In this thesis, we explore the use of stochastic Navier-Stokes equations through the Dynamically Orthogonal (DO) methodology developed at MIT in the Multidisciplinary Simulation, Estimation, and Assimilation Systems Group. Specifically, we examine the effects of the Reynolds number on stochastic fluid flows behind a square cylinder and evaluate computational schemes to do so. We review existing literature, examine our simulation results and validate the numerical solution. The thesis uses a novel open boundary condition formulation for DO stochastic Navier-Stokes equations, which allows the modeling of a wide range of random inlet boundary conditions with a single DO simulation of low stochastic dimensions, reducing computational costs by orders of magnitude. We first test the numerical convergence and validating the numerics. We then study the sensitivity of the results to several parameters, focusing for the dynamics on the sensitivity to the Reynolds number. For the method, we focus on the sensitivity to the: resolution of in the stochastic subspace, resolution in the physical space and number of open boundary conditions DO modes. Finally, we evaluate and study how key dynamical characteristics of the flow such as the recirculation length and the vortex shedding period vary with the Reynolds number.