Datasets

All direct numerical simulations (DNS) are performed under the following conditions.

The density field $\rho$ satisfies \begin{equation} \rho\left(\mathbf{x},t\right)=\rho_{0}+\bar{\rho}\left(z\right)+\rho’\left(\mathbf{x},t\right),\ \ \bar{\rho}\left(z\right)=-Az, \end{equation} where $z$ is taken to point upward in the opposite direction of gravity, $\rho_0$ denotes a reference density, $A$ is a constant describing the linear background vertical density gradient, and $\rho’$ denotes (turbulent) fluctuations away from the background state.

The Navier-Stokes equations are then solved under the Boussinesq approximation; specifically, variations in density away from the background reference value $\rho_0$ are assumed to be small and such density variations are thus only accounted for when computing the buoyancy force. The resulting equations of motion that are solved numerically are:

Equation (2) is the incompressibility condition arising from mass conservation under the Boussinesq approximation, with $\mathbf{u}=\left(u,v,w\right)$ describing the fluid velocity. Equation (3) is the momentum equation, where the pressure field has been decomposed as $p\left(\mathbf{x},t\right)=\bar{p}\left(z\right)+p’\left(\mathbf{x},t\right)$ such that $\partial\bar{p}/\partial z=-\rho_{0}g$ (i.e. $\bar{p}\left(z\right)$ is the hydrostatic pressure arising from the ambient density $\rho_0$). The gravitational acceleration is $g$, the viscosity $\mu$, and $\mathbf{F}$ denotes a body force that is used in some simulations to sustain turbulence in a steady-state (i.e. statistically-stationary turbulence). Equation (4) is an advection-diffusion equation governing density perturbations $\rho’$ away from the linear background stratification, with $\kappa$ denoting the molecular diffusivity.

A fractional-step pseudospectral method is used to solve equations (2-4) in dimensional form in a three-dimensional box. The variables $\mathbf{u}=\left(u,v,w\right)$ and $\rho’$ all satisfy periodic boundary conditions. Values of the following constants are provided with each dataset in a text file “params.txt”:

Note that parameter values need only be chosen such that the dimensionless groups of relevance in characterizing the turbulent flow (i.e. $Re$, $Fr$ etc, described below) emerge to be the desired values. Therefore, it is often simpler to use scaled values of the parameters (such as $g=1$ instead of $g=9.81 \ m/s^2$, for example) when running the simulations.

The equations of motion (2-4) may be non-dimensionalized using

\begin{equation} \hat{x}=x/L,\ \hat{u}=u/U,\ \hat{t}=t\left(U/L\right),\ \hat{p}=p/\left(\rho_{0}U^{2}\right),\ \hat{\rho}=\rho/\left(\left[d\bar{\rho}/dz\right]L\right), \end{equation} where hats indicate non-dimensional quantities. The equations of motion in non-dimensional form are then

where \begin{equation} Re=\frac{UL}{\nu},\ \ Fr=\frac{U}{LN},\ \ Pr=\frac{\nu}{\kappa}, \ \ N=\sqrt{-\frac{g}{\rho_{0}}\frac{\partial\bar{\rho}}{\partial z}}. \end{equation}

The three key dimensionless groups that emerge are the Reynolds ($Re$), Froude ($Fr$), and Prandtl ($Pr$) numbers. The parameter $N$ appearing in the Froude number is the buoyancy frequency $N=2\pi f$ charactizing the vertical oscillation frequency $f$ of a parcel of fluid that is perturbed in a local background density gradient $\partial\bar{\rho} / \partial z$. We note that in characterizing our datasets, we typically use a scaled version of the Froude number $Fr=\frac{2\pi U}{NL}$, as it is $f$ (with dimensions of $1/s$) that characterizes the oscillation frequency of the fluid parcel.

The above dimensionless groups may be combined to form other informative groups, and can be physically interpreted in terms of ratios of length or time scales characterizing the flow, as summarized in the table below.

The figure below illustrates the range of relevant length scales that must be resolved in a DNS of stratified turbulence.

Simulation outputs are stored as binary files of the form “{var}.{#}”, where {var}={u,v,w,r} corresponds to the three velocity components $(u,v,w)$ and the perturbed density field $\rho’$ away from the background density gradient. The {#} corrsponds to the timestamp of the simulation, which may be converted into a physical time based on the reported timestep.

The three-dimensional grid for a given variable and timestep is stored as a one-dimensional string of single-precision numbers in order of $x$ (fastest changing index), $y$, $z$ (slowest changing index). To easily detect errors in reading and reshaping the raw data into a three-dimensional grid, two zeros are padded onto the end of the $x$-dimension. The datafiles thus have size $\left(N_{x}+2\right)\times N_{y}\times N_{z}$. Python codes are provided to read and save a given binary file as a three-dimensional array (e.g. as a NumPy array). For larger datasets, it is advisable to use a memory map which allows one to easily extract slices of the fields without loading the entire three-dimensional array into memory.