1. Introduction
Owing to the associated enhanced rates of irreversible scalar mixing, stratified turbulence is a critical process in the Earth's atmosphere and oceans, impacting both weather and climate, and in the interiors of stars and gaseous planets, affecting their long-term evolution. Assuming that the buoyancy of the fluid, whether liquid or gaseous, is controlled by a single scalar field, which could be temperature or the concentration of a single solute, the dimensionless Boussinesq equations governing the fluid motions (Spiegel & Veronis Reference Spiegel and Veronis1960) are
$$\begin{gather} \frac{\partial{\boldsymbol{u}}}{\partial{t}} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}={-} \boldsymbol{\nabla} p + \frac{b}{Fr^2} \boldsymbol{e}_z + \frac{1}{Re} \nabla^{2} \boldsymbol{u} + \boldsymbol{F}_h, \end{gather}$$
$$\begin{gather}\frac{\partial{b}}{\partial{t}} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} b + w = \frac{1}{Pe} \nabla^2 b, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0. \end{gather}$$
Here, 
$\boldsymbol {u}=(u,v,w)$ is the velocity field expressed in units of 
$U^*$ (where 
$U^*$ is a characteristic horizontal velocity of the large-scale flow), 
$t$ is the time variable in units of 
$L^*/U^*$ (where 
$L^*$ is a characteristic large horizontal scale of the flow), 
$p$ is the pressure fluctuation away from hydrostatic equilibrium in units of 
$\rho _m^* U^{*2}$ (where 
$\rho _m^*$ is the mean density of the fluid), and 
$b$ is the deviation of the buoyancy field away from a linearly stratified background, expressed in units of 
$L^* N^{*2}$ (where 
$N^*$ is the buoyancy frequency of the stable stratification). The 
$w$ term in the buoyancy equation thus represents the vertical advection of the background stratification. Note that here and throughout this paper, starred quantities are dimensional while non-starred quantities are non-dimensional. The flow is assumed to be driven by a non-dimensional divergence-free horizontal force 
${\boldsymbol F}_h$, which only varies on large spatial scales and long time scales. The unit vector 
$\boldsymbol {e}_z$ points in the direction opposite to gravity. We note that the validity of the Boussinesq approximation in the context of gaseous atmospheric and astrophysical flows (Spiegel & Veronis Reference Spiegel and Veronis1960) must be verified a posteriori, by checking that the characteristic vertical scale of the flow remains much smaller than a pressure, temperature or density scale height. It is assumed here that 
$U^*$ is always much smaller than the sound speed.
The usual dimensionless governing parameters of the flow emerge; namely, the Reynolds number 
$Re = U^*L^*/\nu ^*$, the Péclet number 
$Pe = U^*L^*/\kappa ^*$ and the Froude number 
$Fr = U^*/N^*L^*$, where 
$\nu ^*$ is the kinematic viscosity, and 
$\kappa ^*$ is the buoyancy diffusivity, while the Prandtl number, 
$Pr =\nu ^*/\kappa ^*= Pe/Re$ is a property of the fluid. Typically, 
$Pr \sim O(1)$ in air and water, but is very small in astrophysical fluids (of order 
$10^{-2}$ in degenerate plasmas and liquid metals, and much smaller in non-degenerate stellar plasmas, see Lignières Reference Lignières2020).
As the stratification increases (
$Fr \rightarrow 0$), vertical motions are increasingly suppressed and restricted to small characteristic vertical scales 
$l_z = O(\alpha )$, where the emergent aspect ratio 
$\alpha = l_z^* / L^*$ is an increasing function of 
$Fr$ but could also depend on 
$Re$ and 
$Pe$. In the limit 
$(\alpha,Fr) \rightarrow 0$, asymptotic analysis has successfully been used to derive reduced equations for stratified turbulence and to gain insight into its properties. Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), following Billant & Chomaz (Reference Billant and Chomaz2001), proposed an asymptotic reduction in which the vertical coordinate 
$z$ is rescaled as 
$\zeta = z / \alpha$ (with 
$\alpha \ll 1$), while the horizontal coordinates 
${\boldsymbol x}_h=(x,y)$ remain of order unity. Accordingly, all dependent variables 
$q$ are expressed as 
$q({\boldsymbol x}_h,\zeta,t;\alpha )$. As in Klein (Reference Klein2010), we refer to this type of model, which is designed to capture the essence of a scale-specific process, as a single-scale asymptotic (SSA) model. This vertical rescaling, when used in (1.1), reveals the importance of the emergent buoyancy Reynolds and Péclet numbers, defined here as
\begin{equation} Re_b = \alpha^2 Re \quad \text{and} \quad Pe_b = \alpha^2 Pe, \end{equation}
respectively. Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) showed that balancing the mass continuity equation in the limit 
$\alpha \rightarrow 0$ requires 
$w = O(\alpha )$. When 
$Re_b$ and 
$Pe_b$ are at least 
$O(1)$, dominant balance in the buoyancy equation implies 
$b = O(\alpha )$. Finally, the vertical component of the momentum equation reduces to hydrostatic equilibrium when 
$\alpha \rightarrow 0$, yielding 
$\alpha = Fr$, as first argued in the inviscid and non-diffusive case by Billant & Chomaz (Reference Billant and Chomaz2001).
More recently, Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024) noted that in the limit of 
$Pr \ll 1$, it is possible to have a regime in which 
$Pe_b \ll 1 \le Re_b$. They demonstrated that, in this case, the SSA model and corresponding asymptotic expansion reveal instead that 
$w = O(\alpha )$ and 
$b = O( \alpha Pe_b)$, with 
$\alpha = (Fr^2/Pe)^{1/4}$ (see also Lignières Reference Lignières2020; Skoutnev Reference Skoutnev2023).
Crucial to the SSA theory is the notion that every component of the flow is strongly anisotropic, with large horizontal scales and a small vertical scale. In this model, therefore, the vertical fluid motions are primarily driven by the divergence of the horizontal flow, as illustrated schematically in figure 1. Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), however, noted that the SSA theory ignores the possibility that isotropic motions with small horizontal scales may also exist and, in fact, are commonly seen in numerical simulations of stratified turbulence at sufficiently large 
$Re_b$ (cf. Maffioli & Davidson Reference Maffioli and Davidson2016; Cope, Garaud & Caulfield Reference Cope, Garaud and Caulfield2020; Garaud Reference Garaud2020). They proposed a new asymptotic reduction that explicitly incorporates two horizontal scales and two time scales, such that all dependent variables are expressed as
\begin{equation} q({\boldsymbol x}_f,{\boldsymbol x}_s,\zeta,t_f,t_s;\alpha),\quad \text{where } {\boldsymbol x}_s = {\boldsymbol x}, {\boldsymbol x}_f = {\boldsymbol x}/\alpha,\ t_s = t,\ t_f = t/\alpha , \end{equation}
where the subscripts 
$s$ and 
$f$ are used to denote slow and fast scales, respectively. Again following Klein (Reference Klein2010), we refer to the resulting reduced equations as a multiscale asymptotic (MSA) model. Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) showed that these definitions imply that the large-scale motions remain strongly anisotropic with an aspect ratio 
$\alpha$, as in the SSA model, but can coexist with isotropic small-scale motions that evolve on the fast time scale 
$t_f$, and vary on the small vertical coordinate 
$\zeta$ and the ‘fast’ horizontal coordinate 
${\boldsymbol x}_f$. These small-scale motions are self-consistently driven by an instability of the local vertical shear emergent from the larger-scale horizontal flow (see figure 1), and are gradually stabilized as the stratification increases at fixed 
$Re$. They essentially disappear beyond a certain threshold, at which point the MSA model naturally recovers the SSA model and its predicted scalings. For the sake of clarity, however, we refer in what follows to the SSA model and its scalings whenever small horizontal scales are dynamically negligible, and to the MSA model and its scalings whenever they are dynamically important, even though the MSA model does in fact naturally cover both cases.
Illustrations and summary of the SSA and MSA model predictions for 
$w$ and 
$l_z$ in both the non-diffusive and diffusive regimes. Horizontal eddies are shown in red, and vertical eddies are shown in blue.
In the asymptotic limit where 
$Re_b \ge O(1)$ and 
$Pe_b \ge O(1)$, Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) found that 
$w = O(\alpha ^{1/2})$, 
$b = O(\alpha )$ and 
$\alpha = Fr$ when 
$\alpha \rightarrow 0$. Their scaling prediction for 
$w$ thus deviates substantially from that of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), but recovers that of Riley & Lindborg (Reference Riley and Lindborg2012) albeit using different arguments (for details see Shah et al. Reference Shah, Chini, Caulfield and Garaud2024). That scaling has been tentatively validated by Maffioli & Davidson (Reference Maffioli and Davidson2016) in run-down direct numerical simulations (DNS) of stratified turbulence.
Extending the MSA theory to the low-
$Pr$ case, Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024) recovered the results of Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) when 
$Pe_b \geq {O}(1)$. They also found that 
$\alpha = Fr$ and 
$w = O(\alpha ^{1/2})$ both continue to hold in an ‘intermediate’ regime where 
$O(\alpha ) \le Pe_b \ll 1$. However, when 
$Pe_b \ll \alpha$, a new fully diffusive regime emerges in which 
$w = O(\alpha ^{1/2})$, 
$b = O(Pe_b \alpha ^{1/2})$ and 
$\alpha = (Fr^2/Pe)^{1/3}$. As for the 
$Pr={O}(1)$ scenario, the scaling predictions of the MSA theory differ substantially from those emerging from the low-
$Pe_b$ limit of the SSA theory but recover them when small scales are absent. The various theories and their predicted scalings for 
$Re_b \ge O(1)$, with 
$Pe_b \ge O(\alpha )$ or 
$Pe_b \ll \alpha$, respectively, are summarized in figure 1.
Therefore, an interesting question is whether evidence for these scaling laws can be found in DNS data. Recently, two series of DNS were presented by Cope et al. (Reference Cope, Garaud and Caulfield2020) and Garaud (Reference Garaud2020), respectively, which solved (1.1) with 
$\boldsymbol {F}_h = \sin (y) \boldsymbol {e}_x$ (where 
$\boldsymbol {e}_x$ is a unit vector in the streamwise, i.e. 
$x$, direction; see § 2 for further details). In their 
$Pe < 1$ simulations (where by construction 
$Pe_b \ll \alpha$), Cope et al. (Reference Cope, Garaud and Caulfield2020) found that the vertical length scale of the turbulent motions scales as 
$\alpha = (Fr^2/Pe)^{1/3}$, validating the predictions of Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024) in that limit. This scaling, however, was not as clearly evident in the high-
$Pe$ but low-
$Pe_b$ data of Garaud (Reference Garaud2020). One potential explanation is that 
$Re_b$ is relatively low in these simulations (which have 
$Pr = 0.1$, so 
$Re_b = 10Pe_b$), implying viscous effects are not necessarily negligible. In the limit of high 
$Pe_b$, Garaud (Reference Garaud2020) was unable to find evidence for the 
$\alpha = Fr$, 
$w = O(Fr)$ scaling of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) and was unaware at the time of the scaling 
$w = O(Fr^{1/2})$ obtained by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), proposing instead on empirical grounds that 
$w \propto \alpha = Fr^{2/3}$ provides the best fit to the data. The apparent discrepancy between Garaud's data and previous models therefore prompts us to analyse some new DNS results and to revisit the available data from Cope et al. (Reference Cope, Garaud and Caulfield2020) and Garaud (Reference Garaud2020) in the light of the MSA models of stratified turbulence recently derived by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) and Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024).
2. Comparison of theory with DNS
The SSA and MSA theories differ primarily in their predictions for the characteristic vertical length scale of the flow (or equivalently, 
$\alpha$) and for the characteristic vertical velocity. We are therefore interested in comparing these predictions to the data. In practice, however, the characteristic vertical length scale is a relatively difficult quantity to extract from the DNS, as there is no unique and universally accepted definition. Consequently, we focus solely on comparing the theoretical predictions for the characteristic vertical velocity of the flow to the root-mean-square (r.m.s.) of the 
$w$ field because that quantity is both well-defined and easy to compute.
In what follows, we extend and re-analyse the datasets presented in Cope et al. (Reference Cope, Garaud and Caulfield2020) and Garaud (Reference Garaud2020). Both studies performed DNS of the set of non-dimensional equations (1.1) with 
$\boldsymbol {F}_h = \sin (y) \boldsymbol {e}_x$ in a triply periodic domain of size 
$L_x = 4{\rm \pi}, L_y = 2{\rm \pi}, L_z = 2{\rm \pi}$ using the PADDI code (Traxler et al. Reference Traxler, Stellmach, Garaud, Radko and Brummell2011). With this choice, 
$L^*$ is simply the inverse horizontal wavenumber 
$k^*_f$ of the forcing, and the dimensional velocity scale 
$U^*$ is 
$(F_0^*/\rho ^*_m k^*_f)^{1/2}$, where 
$F^*_0$ is the dimensional amplitude of the forcing. In both studies, the streamwise length of the domain 
$L_x = 4{\rm \pi}$ was chosen to be close to the period of the fastest-growing mode of the horizontal shear instability of the Kolmogorov flow driven by that body force (see e.g. Cope et al. Reference Cope, Garaud and Caulfield2020). This ensures the natural generation of large-scale flows in the horizontal direction, whose nonlinear evolution then generates flows on 
$O(1)$ scales in both the 
$x$ and 
$y$ directions. As illustrated in the Appendix, we have verified that the results in the non-viscous, non-diffusive regime are essentially independent of 
$L_x$ and 
$L_y$ as long as the domain is large enough to allow that primary mode of horizontal shear instability to grow.
We focus on the highest-Reynolds-number simulations of Cope et al. (Reference Cope, Garaud and Caulfield2020) and Garaud (Reference Garaud2020), which were performed for 
$Re = 600$. Note that here and in these papers 
$Re$ is defined using the inverse wavenumber of the forcing, and thus is a factor of 
$2{\rm \pi}$ smaller than that of simulations which use the box size as the unit of length instead. Similarly, 
$Fr$ is a factor of 2
${\rm \pi}$ larger here than if we had used the domain size instead. Cope et al. (Reference Cope, Garaud and Caulfield2020) presented a range of DNS for 
$Pe \le 0.1$ and low 
$Fr$ (using the parameter 
$B = Fr^{-2}$ to characterize the stratification). They also ran a few simulations in the asymptotically low-
$Pe$ regime (Lignières Reference Lignières1999), called the LPN regime hereafter, where the buoyancy equation is replaced by 
$w = Pe^{-1} \nabla ^2b$. Garaud (Reference Garaud2020) presented DNS with 
$Re = 600$, 
$Pe = 60$ (i.e. 
$Pr = 0.1$) and low 
$Fr$. For each distinct value of 
$(Re,Pe)$, a first simulation with 
$Fr = 0.33$ was initialized from 
$b = 0$ and 
${\boldsymbol u} = \sin (y) {\boldsymbol e}_x$, plus small amplitude white noise (Cope et al. Reference Cope, Garaud and Caulfield2020; Garaud Reference Garaud2020). Subsequent simulations at higher or lower values of 
$Fr$ were restarted from the end-point of that first run, to bypass the long transient required for the primary horizontal shear instability to develop. Each simulation was integrated until a statistically stationary state was reached, lasting at least 100 time units (see the Appendix for sample time series and a justification for this choice). We have confirmed that the results are independent of the initialization provided the simulations are integrated for at least this period of time. In some cases, we had to further extend the original DNS from Cope et al. (Reference Cope, Garaud and Caulfield2020) or Garaud (Reference Garaud2020) to have a sufficiently long stationary time series. The quantity 
$w_{rms} = \langle w^2\rangle _t^{1/2}$, where 
$\langle {\cdot } \rangle _t$ denotes a volume and time average, was then measured in that statistically stationary state using the extended data.
To complement this dataset, we have run additional simulations at 
$Re = 1000$, 
$Pe = 100$. These DNS have twice the spatial resolution of those of Cope et al. (Reference Cope, Garaud and Caulfield2020) and Garaud (Reference Garaud2020) and, thus, have only been integrated for up to 50 time units in the statistically stationary regime. In addition, the full fields are too large to be saved regularly, so we have saved two-dimensional slices through the data in the 
$(x,y)$, 
$(y,z)$ and 
$(x,z)$ planes. These simulations are only used for visualizations and to assess the influence of viscosity by comparison with the 
$Re = 600$ results. Note that for these new 
$Re = 1000$ runs, and for all previously published ones in Cope et al. (Reference Cope, Garaud and Caulfield2020) and Garaud (Reference Garaud2020), we have ensured that the product of the maximum resolved wavenumber and the Kolmogorov scale is always greater than one, ensuring that the flow field (and the buoyancy field, since 
$Pr \le 0.1$ in all cases) is resolved down to the dissipation scale (see Cope et al. Reference Cope, Garaud and Caulfield2020; Garaud Reference Garaud2020, for details).
We compare the 
$w_{rms}$ data and the various theories in figure 2(a,b). Panel (a) shows 
$w_{rms}$ as a function of 
$Fr^{-1}$, measured for the 
$Pe = 60$, 
$Re=600$ runs (green symbols), and for the 
$Pe = 100$, 
$Re = 1000$ runs (orange symbols). Note that 
$Pr = 0.1$ in both cases. We refer to these simulations as ‘non-diffusive’ because 
$Pe$ is large. The fact that the measured values of 
$w_{rms}$ are identical for the two sets of simulations at different 
$Re$ demonstrates that viscous effects are negligible, at least for 
$Fr^{-1} \le 20$. Panel (b) shows the results of the suite of experiments at 
$Pe=0.1$, 
$Re = 600$ (purple symbols), for which 
$Pr = 0.1/600$. We refer to these simulations as ‘diffusive’ runs, because 
$Pe$ is small. Note that some of these runs were actually integrated using the LPN regime equations instead (square symbols). In that case, the relevant input parameters are 
$Re$ and 
$\chi = Pe/Fr^2$ (
$=BPe$ in the notation of Cope et al. Reference Cope, Garaud and Caulfield2020). To obtain the corresponding value of 
$Fr$ for a given 
$\chi$, simply note that 
$Fr = \sqrt {Pe / \chi }$ for a given 
$Pe$. Illustrated as well in the same plots are the various theoretical predictions for the vertical velocity: the red line in each panel corresponds to the SSA scalings, while the blue line corresponds to the MSA scalings.
(a,b) Comparison between the model predictions and the data for non-diffusive simulations with 
$Pr = 0.1$ and two different values of 
$Re$ (a) and diffusive simulations with 
$Re = 600$, 
$Pe = 0.1$ (b). Symbols show 
$w_{rms}$ extracted from the DNS and error bars show the standard deviation of its temporal variability. Squares in (b) denote LPN simulations (see main text for details). The blue and red lines in (a,b) show the MSA and SSA scaling predictions, respectively. (c,d) Regime diagrams for stratified turbulence at 
$Pr = 0.1$ (c) and 
$Pr = 0.1/600 \simeq 0.00017$ (d), adapted from Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024). Grey regions support isotropic motions. White regions are viscously controlled (
$Re_b \le 1$). Green regions support non-diffusive anisotropic stratified turbulence (
$Pe_b \ge \mathit {O}(1)$), and purple regions support diffusive anisotropic stratified turbulence (
$Pe_b \ll \alpha$). The yellow regions are in the ‘intermediate’ regime of Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024) (
$\mathit {O}(\alpha ) \le Pe_b \ll 1$). Horizontal arrows show the transects through the regimes corresponding to the panels above.
Figure 2(c,d) shows for comparison the expected regime diagrams for the corresponding values of 
$Pr$ in each case, based on the asymptotic theory of Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024). The coloured horizontal arrows show the transect taken through parameter space for each series of DNS shown in figure 2(a,b). The background colours show the expected regime: isotropic motions with 
$\alpha \simeq 1$ (grey), non-diffusive anisotropic turbulence (green), intermediate regime (yellow), diffusive anisotropic turbulence (violet) and viscous regime (white). The same background colours in figure 2(a,b) show the expected regime transitions as a function of 
$Fr^{-1}$ at the value of 
$Pe$ corresponding to the transect taken. We note that while Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024) distinguished the non-diffusive and intermediate regimes, these have the same predicted scalings for 
$\alpha$ and 
$w$; in any case, the intermediate regime does not span a large region of parameter space at 
$Pr = 0.1$ and would be difficult to identify even if the scaling laws differed.
Examination of figure 2(a,b) confirms that none of the theories apply when the stratification is weak so the flow is isotropic on all scales (
$\alpha \simeq 1$, grey region), or in the viscous regime (white region), where 
$Re_b \le 1$. This is, of course, as expected. However, we also see that neither the SSA nor the MSA model predictions fit the data in the entire region where they are supposedly valid (i.e. the green/yellow regions in the non-diffusive case, and the purple region in the diffusive case). Instead, we find that the MSA predictions appear to be better at weaker stratifications (higher 
$Fr$) while the SSA predictions appear to be better at higher stratification (lower 
$Fr$), when 
$Re$ is fixed.
3. Turbulent patches vs quiescent flow
To gain insight into the applicability of the predicted scalings, we examine the actual flow field more closely. Figure 3(a–f) show snapshots of 
$u$ and 
$w$ in three different high-resolution DNS at 
$Pe = 100$ and 
$Re = 1000$, with 
$Fr$ decreasing from approximately 
$0.18$ to approximately 
$0.058$. It is clear that while the 
$Fr \simeq 0.18$ case is fully turbulent, the more strongly stratified 
$Fr \simeq 0.058$ case is not as the turbulence is localized to small ‘patches’ (see e.g. the regions of high 
$|w|$). Similar findings were reported by Cope et al. (Reference Cope, Garaud and Caulfield2020) in their low-
$Pe$ simulations in the regime that they named ‘stratified intermittent’ (see their figure 6) and are also evident in Garaud (Reference Garaud2020); see the volume renderings of 
$u$ and 
$w$ at 
$Fr = 0.05$ (
$B =400$) in her figure 1 for instance. Figure 3(g–i) shows the kinetic energy spectra of the horizontal flows (red lines) and of the vertical flow (blue lines) for the same simulations. Different lines correspond to different instants in time, in order to illustrate the intrinsic variability of the spectra.
(a–f) DNS snapshots of 
$u$ and 
$w$ in the 
$y=0$ plane for 
$Re = 1000$, 
$Pe = 100$ and various Froude numbers, with stratification increasing from top to bottom. (g–i) Kinetic energy spectra of the horizontal flows (red lines) and of the vertical flows (blue lines) as a function of the horizontal wavenumber 
$k_h$, for the same three Froude numbers, with stratification increasing from left to right. Each line corresponds to a particular instant in time.
These results clearly illustrate the coexistence of large and small horizontal scales, with 
$u$ increasingly dominated by large scales with subdominant small scales as stratification increases, and 
$w$ increasingly dominated by small scales with subdominant large scales (cf. Riley & Lindborg Reference Riley and Lindborg2012). We see from the spectra in particular that the large scales are highly anisotropic. The horizontal flows have a kinetic energy spectrum proportional to 
$k_h^{-3}$ (where 
$k_h$ is the horizontal wavenumber), consistent with oceanic observations (e.g. Klymak & Moum Reference Klymak and Moum2007; Falder, White & Caulfield Reference Falder, White and Caulfield2016) and the classical empirical ‘Garrett-Munk’ spectrum for internal waves (Garrett & Munk Reference Garrett and Munk1975). It is also consistent with observations of turbulence in the stratosphere on 
$\mathit {O}$(10 km) scales (Lilly & Lester Reference Lilly and Lester1974). Meanwhile, the kinetic energy spectrum of vertical motions is a weakly increasing function of 
$k_h$ at large scales, which peaks at a value of 
$k_h$ that appears to scale as 
$Fr^{-1}$ (with the limited data available). This tentatively shows that 
$l_z^* \simeq Fr L^*$ is indeed the injection scale for the vertical motions in these non-diffusive strongly stratified simulations, consistent with the interpretation that they arise from shear instabilities of the layerwise horizontal flow on that scale.
At smaller scales, we see that the flow becomes much more isotropic. In the 
$Fr \simeq 0.18$ and 
$Fr =0.1$ cases, for which 
$Re_b = Fr^2 Re \simeq 33$ and 
$10$, respectively, the isotropic small-scale flow seems to have a standard 
$k_h^{-5/3}$ spectrum until viscous effects come into play. In the 
$Fr \simeq 0.058$ case, by contrast, the spectrum of the small-scale flow is steeper than 
$k_h^{-5/3}$, indicating that the turbulence is somewhat suppressed. This observation is not surprising, because at these parameter values 
$Re_b \simeq 3.3$, and so viscous effects are likely to suppress inertial range dynamics that are unaffected by stratification, as well as suppressing the local vertical shear instability except in regions where the shear is exceptionally strong, leading to the patchiness of the flow observed in the snapshots.
The snapshots in the more strongly stratified case further reveal that the small horizontal scales are only dominant within the turbulent patches and essentially disappear outside of these patches. As such, the distinct MSA model scalings are only expected to apply within the turbulent patches. In the more orderly layer-like flow outside of those patches, the SSA scalings – which coincide with the MSA model predictions in regions where small scales are not excited – should hold.
To verify this interpretation quantitatively, we sought to identify a reliable diagnostic for the turbulent patches, i.e. regions where the flow exhibits small horizontal scales. It is common to use the enstrophy 
$|{\boldsymbol {\omega }}|^2$ as a diagnostic for turbulence, where 
$\boldsymbol {\omega } = \boldsymbol {\nabla } \times \boldsymbol {u}$ is the flow vorticity. Indeed, the turbulent cascade to small scales implies that enstrophy must be large within the patches. However, enstrophy turns out to be an inappropriate diagnostic for our purpose because it can also be large in the layer-like regions of strong vertical shear outside of the turbulent patches, such as those described by the SSA model. This fact is illustrated in figure 4(a), which shows the enstrophy field in a particular snapshot of a strongly stratified simulation, and can be understood as follows. According to Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) and Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024), 
${\boldsymbol u} = \bar {\boldsymbol u} + {\boldsymbol u}'$ where 
$\bar {\boldsymbol u}$ can be thought of as the large-scale anisotropic component of the flow, which varies on the 
$O(1)$ horizontal scales and 
$O(\alpha )$ vertical scale, as in the SSA model. Meanwhile 
${\boldsymbol u}'$ can be thought of as the small-scale isotropic and turbulent component of the flow, which varies on 
$O(\alpha )$ scales in all directions, as in the MSA model. Furthermore, these authors show that 
$\bar u \sim \bar v \sim O(1)$, while 
$\bar w \sim O(\alpha )$, and 
$u'\sim v' \sim w' \sim O(\alpha ^{1/2})$. Accordingly, we find that the horizontal vorticity components are dominated by the contribution from 
$\bar {\boldsymbol u}$, namely 
$\omega _x \sim \bar \omega _x \sim \omega _y \sim \bar \omega _y \sim O(Fr^{-1})$, while the vertical vorticity component is dominated by the contributions from 
${\boldsymbol u}'$, with 
$\omega _z \sim \omega '_z \sim O(Fr^{-1/2})$. For comparison, note that the vertical vorticity associated with the large-scale forcing, and with the mean horizontal flow 
$(\bar u, \bar v)$, is 
$O(1)$ and therefore negligible in comparison with 
$\omega _z'$. As a result, we argue that 
$\omega _z^2$ is a more reliable diagnostic of the small-scale turbulence than the enstrophy. This assertion is confirmed in figure 4(b), which shows 
$\omega _z^2$ for the same snapshot depicted in figure 4(a). We see that the regions of high 
$\omega _z^2$ only highlight the turbulent patches of the flow.
Snapshots of enstrophy (a) and vertical vorticity squared (b) from a simulation at 
$Re = 600$, 
$Pe = 60$ and 
$Fr = 0.05$. The latter is a better diagnostic of the turbulent patches.
In what follows, we therefore define the following quantities, using weighted averages with the weight function 
$\omega _z^2$ to emphasize the turbulent patches and the function 
$\omega _z^{-2}$ to emphasize the more quiescent regions, respectively:
\begin{equation} w^{turb}_{rms} = \frac{\langle w^2 \omega_z^2\rangle_t^{1/2}}{\langle \omega_z^2\rangle_t^{1/2}} \quad \text{and}\quad w^{no turb}_{rms} = \frac{\langle w^2 \omega_z^{{-}2}\rangle_t^{1/2}}{\langle \omega_z^{{-}2}\rangle_t^{1/2} }. \end{equation}
The first can be viewed as the r.m.s. of 
$w$ taken over the turbulent patches, where the distinct MSA scalings should apply. The second can be viewed as the r.m.s. of 
$w$ taken everywhere other than the turbulent patches, where the SSA scalings should apply. Note that the computation of 
$w^{turb}_{rms}$ and 
$w^{noturb}_{rms}$ requires integrals of 
$w^2$, 
$\omega _z^2$, and their product or ratio over the entire volume, which was not one of the simulation diagnostics originally saved. As such, we are unable to extract these quantities from the 
$Re = 1000$ simulations. However, we can compute them from the full-data snapshots regularly saved in the 
$Re = 600$ simulations in both high- and low-
$Pe$ datasets (of which there are usually between 50 and 100 depending on the simulation). The variance is naturally larger than for the 
$w_{rms}$ data because of the smaller amount of data available. We also note that similar results can be obtained by choosing 
$|\omega _z|$ and 
$|\omega _z|^{-1}$ as the weight functions for the turbulent and quiescent regions, respectively. However, 
$|\omega _z|$ is only a factor of 
$Fr^{-1/2}$ larger in the turbulent regions than in the quiescent ones. Since 
$Fr$ in our simulations is not extremely small, we use 
$\omega _z^2 = O(Fr^{-1})$ instead to more clearly identify the turbulent patches.
We present the results in figure 5, with the non-diffusive 
$Re = 600, Pe =60$ simulations in panel (a) and the diffusive 
$Re = 600, Pe = 0.1$ simulations in panel (b). The background colours are the same as in figure 2. The 
$w_{rms}$ data from figure 2 is again shown in green and purple symbols. We plot the 
$w^{turb}_{rms}$ data using blue symbols in both cases, and the MSA scalings for turbulent regions using a blue line for comparison. Similarly, we plot the 
$w^{noturb}_{rms}$ data using red symbols and the SSA scalings using a red line. We see, quite clearly, that each theory fits the data in its respective region of validity – the distinct MSA scalings being valid in the turbulent patches, and the SSA scalings only being valid outside of the turbulent patches. This shows that the transition observed in figure 2, from simulations that appear to satisfy the turbulent MSA scalings better at low stratification to simulations that appear to fit the SSA scalings better at high stratification, primarily is a consequence of the decrease in the volumetric fraction of the domain occupied by the turbulent patches when 
$Fr^{-1}$ increases.
Comparison between models and data for the characteristic vertical velocity at 
$Re = 600$, 
$Pe = 60$ (a) and at 
$Re = 600$, 
$Pe = 0.1$ (b). Green and purple symbols show 
$w_{rms}$ in (a) and (b), respectively. In (a,b), blue symbols show 
$w^{turb}_{rms}$ and should be compared with the turbulent MSA scalings (blue lines), while red symbols show 
$w^{noturb}_{rms}$ and should be compared with the corresponding SSA scalings (red lines).
As the stratification continues to increase, the buoyancy Reynolds number 
$Re_b = \alpha ^2 Re$ eventually decreases below a critical value 
$Re_{b,{crit}} = O(1)$, where viscous effects become dominant. Assuming 
$Re_{b,{crit}} =1$, we show this transition in figure 2 as the line separating the coloured area from the white region, for the MSA model. In the non-diffusive case (panel (a)), 
$\alpha = Fr$, so 
$Re_{b} =1$ is equivalent to 
$Pe = Pr Fr^{-2}$, which is the edge of the yellow region. We see that the data are consistent with this prediction: beyond the viscous transition, 
$w_{rms}$ rapidly drops to very low values consistent with a viscously dominated flow. For the diffusive case (panel (b)), 
$\alpha = (Fr^2/Pe)^{1/3}$ in the MSA model so 
$Re_b = 1$ is equivalent to 
$Pe = Pr^3 Fr^{-4}$, which corresponds to the edge of the purple region. We see that the effects of viscosity appear to become important at slightly weaker stratification than predicted assuming 
$Re_{b,{crit}} =1$, but plausibly attribute this discrepancy to missing 
$O(1)$ constants in the estimates for 
$\alpha$ and/or 
$Re_{b,{crit}}$.
4. Conclusion
In this paper, we have presented a detailed comparison of DNS data with various theoretical predictions for the characteristic vertical velocity of fluid motions in forced stratified turbulence. In particular, we have studied both moderate- and low-Prandtl-number regimes, resulting in a wide range of Péclet numbers at fixed Reynolds number. When buoyancy diffusion is negligible, our results notably provide compelling evidence for the 
$w \propto Fr^{1/2}$ scaling law for stratified turbulence first proposed by Riley & Lindborg (Reference Riley and Lindborg2012) using heuristic arguments and rigorously derived by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) using multiscale asymptotic analysis. In the latter investigation, this scaling law is intrinsically tied to the existence of small-scale isotropic flow motions driven by the emergent vertical shear between larger-scale primarily horizontal eddies, as corroborated by the results from the DNS presented here. The vertical shear instability is gradually stabilized as the buoyancy Reynolds number decreases towards a critical value of order unity, and the small-scale isotropic component of the turbulence becomes confined to localized patches rather than being domain filling. Outside of these turbulent patches, small horizontal scales disappear, and we find that 
$w \propto Fr$ instead, consistent with the model of Billant & Chomaz (Reference Billant and Chomaz2001) and Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007). As 
$Re_b$ gradually approaches unity from above, therefore, the r.m.s. vertical velocity of the flow computed from an average over the whole domain differs from either of these scaling laws, and additionally depends on the volume filling factor of the small-scale turbulence, whose dependence on stratification and 
$Re_b$ is the subject of ongoing work. The patchiness of the small-scale isotropic flow also explains the incorrect conclusion reached by Garaud (Reference Garaud2020) regarding the possible existence of another regime of stratified turbulence where 
$w \propto Fr^{2/3}$. In hindsight, we understand her empirically inferred intermediate scaling as a consequence of the decrease in the volume filled by turbulent patches with increasing stratification at fixed 
$Re$.
At low 
$Pr$, Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024) revised the predictions of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) and Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) to account for the effects of buoyancy diffusion. They showed that the presence of small isotropic motions implies that 
$w \propto (Fr^2/Pe)^{1/6}$, consistent with an early model and DNS data by Cope et al. (Reference Cope, Garaud and Caulfield2020), while in their absence 
$w \propto (Fr^2/Pe)^{1/4}$, consistent with predictions from Lignières (Reference Lignières2020) and Skoutnev (Reference Skoutnev2023). Revisiting the very-low-
$Pr$ DNS of Cope et al. (Reference Cope, Garaud and Caulfield2020) in this new light, we have confirmed both scaling laws within and outside of the turbulent patches, respectively. Finally, the models of Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) and Shah et al. (Reference Shah, Chini, Caulfield and Garaud2024) also predict where in parameter space viscous effects become important. We have confirmed these predictions, too, with our DNS data.
An important caveat of our conclusions is that the simulations were limited to a specific type of horizontal forcing. We believe that the results ought to apply more generally as long as the forcing drives primarily horizontal flows on long time scales and large length scales, but this conjecture will need to be verified in future work.
In summary, this investigation demonstrates that the combination of rigorous multiscale analysis (Chini et al. Reference Chini, Michel, Julien, Rocha and Caulfield2022; Shah et al. Reference Shah, Chini, Caulfield and Garaud2024) with idealized DNS (Cope et al. Reference Cope, Garaud and Caulfield2020; Garaud Reference Garaud2020, and new simulations presented here) can be a powerful tool to identify and validate scaling laws for stratified turbulence across different regions of parameter space. In future work, we will incorporate the effects of rotation and magnetic fields, which must be taken into account for a more realistic description of stratified turbulence in geophysical and astrophysical settings.
Funding
This work uses the Expanse supercomputer at the San Diego Supercomputing Center. P. Garaud thanks the SDSC support team for their help. The authors gratefully acknowledge the Geophysical Fluid Dynamics Summer School (NSF 1829864), particularly the 2018, 2022 and 2023 programs. K.S. acknowledges funding from the J.S. McDonnell Foundation. G.P.C. acknowledges funding from the U.S. Department of Energy through award DE-SC0024572.
Declaration of interests
The authors report no conflict of interest.
Appendix. Numerical considerations
The results presented in this paper require simulations that have achieved a statistically stationary state lasting at least 100 time units. This duration was chosen to ensure that the time series from which flow statistics are computed are sufficiently uncorrelated in time. Indeed, assuming a mean streamwise velocity of 1, the fluid has time to flow approximately 8 times through the domain in 100 time units when the domain length is 
$4{\rm \pi}$. In practice, the mean streamwise velocity ranges from approximately 2 to 4, depending on the input parameters, and the integration time interval is often larger than 100 time units, so the true number of ‘laps’ is generally much higher. Time series of the instantaneous r.m.s. streamwise and vertical velocity for a few selected simulations in both non-diffusive and diffusive regimes are shown in figure 6.
Sample time series of the instantaneous r.m.s. horizontal velocity 
$u_{rms}$ (red lines) and vertical velocity 
$w_{rms}$ (blue lines), as a function of time for various simulations. Note that the starting points of the simulations have been offset to an arbitrary position for ease of visualization. The grey area shows the time-averaging interval used in figures 2 and 5.
We have also verified that the scaling laws obtained are independent of the selected domain size by performing a few simulations in an 
$8{\rm \pi} \times 4{\rm \pi} \times 2{\rm \pi}$ domain with the same body force. With that choice, the domain is sufficiently wide to accommodate two wavelengths of the applied sinusoidal force in the spanwise direction, and two wavelengths of the fastest-growing mode of horizontal shear instability in the streamwise direction. Because of the heavily increased computational cost, we have only run cases for 
$Re = 600, Pe = 60$, and three values of 
$Fr$, and these simulations have been run for a shorter duration. Yet, as demonstrated in figure 7, the same scaling laws are found in the turbulent patches (
$w^{turb}_{rms} \propto Fr^{1/2}$) and outside of the turbulent patches (
$w^{noturb}_{rms} \propto Fr$), respectively. We note that the prefactor is slightly smaller, but this is not surprising given that the meandering flow structure is allowed to be a little different in the larger domain.
(a) Filled symbols show 
$w^{turb}_{rms}$ and 
$w^{noturb}_{rms}$ extracted from simulations in larger computational domains (of size 
$8{\rm \pi} \times 4{\rm \pi} \times 2{\rm \pi}$) and can be compared to those extracted from simulations in regular-sized domains, shown as open symbols. (b,c) Snapshot of the 
$(x,y)$ plane at some arbitrary value of 
$z$ in the statistically stationary state of a simulation at 
$Re = 600$, 
$Pe = 60$, 
$Fr = 0.1$, in a domain of size 
$8{\rm \pi} \times 4{\rm \pi} \times 2{\rm \pi}$.
Open access
