2.1. Synopsis of the reduced system

Equations (2.25), (2.17), (2.27) and (2.23) and (2.20), (2.21), (2.22) and (2.19) comprise a novel, multiscale reduced system. For ease of reference, this system is reproduced here, where, for brevity of notation, the numeric subscripts have been omitted and $\overline {\boldsymbol {\nabla }}$ and $\boldsymbol {\nabla }'$ are used in lieu of $\boldsymbol {\nabla }_{\boldsymbol {x}_\bot }$ and $\boldsymbol {\nabla }_{\boldsymbol {\chi }_\bot }$, respectively.

Mean dynamics

(2.28)\begin{gather} \left(\partial_t+\bar{\boldsymbol{u}}_\bot\boldsymbol{\cdot}\overline{\boldsymbol{\nabla}}+\bar{W}\partial_z\right)\bar{\boldsymbol{u}}_\bot ={-}\overline{\boldsymbol{\nabla}}\bar{p}-\partial_z\left(\overline{W'\boldsymbol{u}'_\bot}\right) +\frac{1}{Re_b}\partial_z^{2}\bar{\boldsymbol{u}}_\bot+\bar{\boldsymbol{f}}_\bot, \end{gather}

(2.29)\begin{gather}0={-}\partial_z\bar{p}+\bar{b}, \end{gather}

(2.30)\begin{gather}\left(\partial_t+\bar{\boldsymbol{u}}_\bot\boldsymbol{\cdot}\overline{\boldsymbol{\nabla}}+\bar{W}\partial_z\right)\bar{b} ={-}\bar{W}-\partial_z\left(\overline{W'b'}\right)+\frac{1}{Pr Re_b}\partial_z^{2}\bar{b}, \end{gather}

(2.31)\begin{gather}\overline{\boldsymbol{\nabla}}\boldsymbol{\cdot}\bar{\boldsymbol{u}}_\bot+\partial_z\bar{W}=0. \end{gather}

Fluctuation dynamics

(2.32)\begin{gather} \left(\partial_\tau+\bar{\boldsymbol{u}}_\bot\boldsymbol{\cdot}\boldsymbol{\nabla}'\right)\boldsymbol{u}'_\bot+W'\partial_z\bar{\boldsymbol{u}}_\bot ={-}\boldsymbol{\nabla}' p'+\frac{Fr}{Re_b}\left(\nabla'^{2}+\partial_z^{2}\right)\boldsymbol{u}'_\bot, \end{gather}

(2.33)\begin{gather}\left(\partial_\tau+\bar{\boldsymbol{u}}_\bot\boldsymbol{\cdot}\boldsymbol{\nabla}'\right)W'={-}\partial_z p'+b'+ \frac{Fr}{Re_b}\left(\nabla'^{2}+\partial_z^{2}\right)W', \end{gather}

(2.34)\begin{gather}\left(\partial_\tau+\bar{\boldsymbol{u}}_\bot\boldsymbol{\cdot}\boldsymbol{\nabla}'\right)b'+W'\partial_z\bar{b} ={-}W'+\frac{Fr}{Pr Re_b}\left(\nabla'^{2}+\partial_z^{2}\right)b', \end{gather}

(2.35)\begin{gather}\boldsymbol{\nabla}'\boldsymbol{\cdot}\boldsymbol{u}'_\bot+\partial_z W'=0. \end{gather}

2.2. Generalized quasilinear (QL) structure of the reduced system

Before attempting to integrate the multiscale reduced system numerically, it is instructive to examine the structure of these equations. In the absence of the Reynolds-stress and buoyancy-flux divergence terms, the mean equations (2.28)–(2.31) are readily seen to reduce to the hydrostatic primitive equations. Heuristically, this set of equations governs the comparably slowly evolving, large-scale layer-like motions that are routinely observed in the LAST regime of strongly stratified turbulence. The correlations arising in these equations, which account for the collective effects of the small-scale flows on the strongly anisotropic large-scale motions, generally must be modelled phenomenologically. By exploiting the scale separation associated with the chosen distinguished limit, however, here we are able to avoid the usual closure difficulties via explicit derivation of the leading-order equations governing the evolution of the fluctuation fields (i.e. (2.32)–(2.35)). Taken together, (2.28)–(2.35) comprise a closed reduced system. Direct calculation confirms that this system conserves energy in the absence of dissipation and forcing.

Inspection of (2.32)–(2.34) reveals that the fluctuations are advected horizontally by the mean flow and can interact with vertical gradients of the mean buoyancy and horizontal velocity fields. Equation (2.33) shows that the fluctuation dynamics are non-hydrostatic, as expected. Moreover, the fluctuation dynamics are ‘quasilinear’ with respect to the slowly varying mean fields; that is, fluctuation/fluctuation nonlinearities are absent from the fluctuation equations, but their feedback on the coarse-grained fields is retained. Indeed, a central outcome of the present work is that the increasingly popular QL approximation for anisotropic turbulent flows (Constantinou et al. Reference Constantinou, Lozano-Durán, Nikolaidis, Farrell, Ioannou and Jiménez2014b; Fitzgerald & Farrell Reference Fitzgerald and Farrell2014; Thomas et al. Reference Thomas, Farrell, Ioannou and Gayme2015; Fitzgerald & Farrell Reference Fitzgerald and Farrell2018; Kim, Billant & Gallaire Reference Kim, Billant and Gallaire2020) and the associated statistical formulations (variously referred to as CE2 and SSST: Tobias, Dagon & Marston Reference Tobias, Dagon and Marston2011; Srinivasan & Young Reference Srinivasan and Young2012; Tobias & Marston Reference Tobias and Marston2013; Constantinou, Farrell & Ioannou Reference Constantinou, Farrell and Ioannou2014a; Ait-Chaalal et al. Reference Ait-Chaalal, Schneider, Meyer and Marston2016; Constantinou, Farrell & Ioannou Reference Constantinou, Farrell and Ioannou2016; Farrell et al. Reference Farrell, Ioannou, Jimenez, Constantinou, Lozano-Duran and Nikolaidis2016; Farrell, Gayme & Ioannou Reference Farrell, Gayme and Ioannou2017) can be formally justified for strongly stratified (and perhaps other) shear flows via multiscale asymptotic analysis. Because the mean fields are independent of the fast coordinates, the fluctuation equations are autonomous in $\boldsymbol {\chi }_\bot$ and there is no direct coupling among ‘fast’ Fourier modes; instead, these modes are coupled only through their contribution to the modification of the mean fields. By retaining slow spatiotemporal variability of the mean fields, our formulation in fact extends the usual QL reduction: precisely this sort of multiscale analysis provided inspiration for the so-called ‘generalized quasilinear’ or GQL approximation (Marston, Chini & Tobias Reference Marston, Chini and Tobias2016), which has been demonstrated to improve the accuracy of QL-based predictions significantly for only a modest increase in model complexity (Child et al. Reference Child, Hollerbach, Marston and Tobias2016; Tobias & Marston Reference Tobias and Marston2017). The reduced equations systematically derived here also can be compared to the one-dimensional (1-D) phenomenological model for fluctuating motions in stratified turbulence introduced by Rorai et al. (Reference Rorai, Minini and Pouquet2014) and extended by Feraco et al. (Reference Feraco, Marino, Pumir, Primavera, Mininni, Pouquet and Rosenberg2018). In particular, we find that while the linear coupling and damping of the temperature and vertical velocity are adequately captured by the phenomenological model, neither the instability nor the nonlinear saturation terms in the fluctuation dynamics, which are fundamental to interactions with the large scales, are properly characterized in the 1-D model.

A final point concerns the role of diffusive effects in the multiscale reduced system. Owing to the anisotropic large-scale layering, vertical diffusion of momentum and buoyancy arises at leading-order in the mean equations. Although asymptotically, i.e. as a function of $Fr$ as $Fr\to 0$, $Re_b={O}(1)$, the effective diffusivity ($\propto 1/Re_b$) can be varied numerically to investigate different dynamical regimes captured by the reduced model. Indeed, many DNS studies have suggested that important regime transitions occur at $Re_b\gtrsim 10$ (Bartello & Tobias Reference Bartello and Tobias2013; Maffioli & Davidson Reference Maffioli and Davidson2016; Lucas & Caulfield Reference Lucas and Caulfield2017; Garanaik & Venayagamoorthy Reference Garanaik and Venayagamoorthy2019; Lang & Waite Reference Lang and Waite2019). In contrast to the mean dynamics, the leading-order fluctuation equations are non-dissipative. Nevertheless, in writing (2.32)–(2.34), we have included formally higher-order Laplacian diffusion terms as a simple means of regularizing the fluctuation dynamics, should large $z$-gradients (e.g. possibly associated with critical layers) emerge. Although justifiable as a composite asymptotic approximation, this approach reintroduces the small parameter $Fr$ into the limit system, and a more careful analysis would be required to ascertain whether other formally weak physical processes also may contribute to the dominant balance of terms in ultra-thin regions in which $z$-gradients become large. With this limitation understood, we proceed using the regularized reduced system (2.28)–(2.35). Nonetheless, as shown explicitly in § 4, the 2-D non-dissipative fluctuation equations can be re-expressed as the Taylor–Goldstein (TG) equation (Craik Reference Craik1985) because, formally, the mean fields are frozen during the fast evolution of the fluctuations. Thus, by inspection, it is clear that the reduced system captures all of the linear instabilities admitted by the TG equation, including the classical Kelvin–Helmholtz and Holmboe instabilities.

3.1. Single time-scale formulation

One straightforward way to simulate the reduced (3.1)–(3.4) is to revert to a single time-scale formulation by making the non-asymptotic replacement $\partial _t=(1/Fr)\partial _\tau$ and by reinterpreting the fast ($\boldsymbol {\chi }_\bot$,$\tau$) average as a strict horizontal average (only). This reformulation sacrifices certain advantages accrued via the multiples scales asymptotic analysis. In particular, the resulting system is numerically stiff owing to the reintroduction of $Fr$, and both the fast fluctuations and the slowly evolving mean fields have to be numerically co-evolved on the fast time scale $\tau$. The advantage of this approach is that there is a clear protocol for discretizing this system in both space and time. In fact, the resulting equations are identical to those that would be obtained via an ad hoc QL approximation of the (2-D) Boussinesq equations, in which flow fields are decomposed into a horizontal (‘streamwise’) average plus a fluctuation about that mean (Fitzgerald & Farrell Reference Fitzgerald and Farrell2018, Reference Fitzgerald and Farrell2019). Regarding the spatial discretization, owing to the QL structure, any set of horizontal Fourier modes can be included. Nevertheless, because this set must be specified a priori, standard grids equispaced in physical or Fourier space generally are employed. As described in the following subsection, a chief virtue of the multiple time-scale formulation is that the intrinsic dynamics of the slow–fast QL system self-selects those wavenumbers and associated Fourier modes to be included. Consequently, in that formulation, the fast spatial domain is effectively infinite in extent; there is no a priori quantization of Fourier modes imposed by the seemingly arbitrary specification of $l_x$.

3.2. Multiple time-scale formulation

In our view, an algorithm that explicitly enforces the time-scale separation between the mean and fluctuation fields in (3.1)–(3.4) is preferable. The method for integrating slow–fast QL systems with fast instabilities introduced by Michel & Chini (Reference Michel and Chini2019) leverages the linearity and autonomy of the fluctuation dynamics on the fast time scale, which here implies that both $\psi '$ and $b'$ depend on the fast time $\tau$ only through a term of the form $e^{\sigma \tau }$, where $\sigma$ is the (complex-valued, linear) growth rate. For the Reynolds stress divergence to remain finite in the mean field equations (3.1) and (3.2), either the real part of the growth rate or the amplitude of the fluctuations therefore must vanish. The crucial point, as we demonstrate below, is that for forced strongly stratified flows, the fluctuation amplitude (if non-zero) is then slaved to the mean fields to ensure that the real part of the growth rate $\sigma _r = 0$; that is, to ensure that the Reynolds stress divergence modifies the evolution of $\bar {u}$ and $\bar {b}$ so that the mean fields always evolve (slowly) on a marginal-stability manifold. This constrained evolution is a mathematical manifestation of self-organized criticality, a physical attribute increasingly being associated with stratified turbulent shear flows (Salehipour et al. Reference Salehipour, Peltier and Caulfield2018; Smyth et al. Reference Smyth, Nash and Moum2019), although the argument that stratified flows adjust to a marginally stable ‘pseudo-equilibrium’ dates back at least to Turner (Reference Turner1973). The remainder of this section describes in detail the resulting multiscale algorithm.

First, note that the fluctuation system (3.3)–(3.4) can be treated as a linear homogeneous eigenvalue problem by seeking modal solutions with real-valued wavenumber $k(t)$, (complex-valued) growth rate $\sigma (\bar {\boldsymbol {G}},k)$, where $\bar {\boldsymbol {G}}$ denotes dependencies on the mean fields and their derivatives (i.e. $\bar {u}$, $\partial _z^{2} \bar {u}$ and $\partial _z \bar {b}$), and complex-valued amplitude $A$ of slowly varying magnitude $\vert A(t)\vert$:

(3.5)\begin{gather} \psi'(\chi, z, \tau; t) = A(t) \hat{\varPsi}(z;t) \exp({\sigma \tau + \textrm{i} k(t) \chi}) + \mathrm{c.c.}, \end{gather}

(3.6)\begin{gather}b'(\chi, z, \tau; t) = A(t) \hat{b}(z;t) \exp({\sigma \tau + \textrm{i} k(t) \chi}) + \mathrm{c.c.}, \end{gather}

where $\sigma =\sigma _r+i\sigma _i$ for real $\sigma _r$ and $\sigma _i$, $\mathrm {c.c.}$ denotes complex conjugate, and we have indicated explicitly that the vertical eigenfunctions $\hat {\varPsi }$ and $\hat {b}$ also may vary slowly with time. Recalling that the fluctuation horizontal velocity and buoyancy are ${O}(\epsilon )$ relative to the corresponding mean fields, it is clear from these expressions that the uniformity of the posited asymptotic expansions in (2.6) would be lost if $\sigma _r>0$. Indeed, in that situation, the fluctuations would grow exponentially fast – on the fast time scale – while the mean fields remained locally frozen in time. Accordingly, as an asymptotic uniformity condition, we seek a solvability constraint that ensures $\sigma _r\le 0$. In practice, there may be ‘instants’ during the slow evolution where this condition cannot be satisfied, in which case the fluctuations will attain large (but finite) amplitude and the ${O}(1)$ mean fields will respond rapidly; i.e. on the fast time scale. In these situations, temporal scale separation is transiently lost. Although the algorithm we develop can be modified to incorporate this intermittent bursting behaviour properly (Ferraro Reference Ferraro2019), this extension proves unnecessary for the parameter regime explored in the present investigation and is not described here. In contrast, $\sigma _r<0$ presents no conceptual or pragmatic difficulty, as disturbances are exponentially damped on the fast time scale and, therefore, the amplitude $A$ is self-consistently set to zero. If $\sigma _r=0$, then the averaging required for evaluation of the Reynolds stress divergence terms in the mean-field equations (e.g. $\partial _z (\overline {\partial _z \psi ' \partial _\chi \psi '})$ in (3.1)) is well defined. Specifically,

(3.7)\begin{gather} \partial_z\left(\overline{\partial_z \psi' \partial_\chi\psi'}\right) = {\vert A(t) \vert}^{2} \textrm{i}k \partial_z \left( \hat{\varPsi} \partial_z \hat{\varPsi}^{*} - \hat{\varPsi}^{*} \partial_z \hat{\varPsi} \right) \equiv {\vert A(t) \vert}^{2} \mathrm{RS}_u, \end{gather}

(3.8)\begin{gather}\partial_z\left(\overline{b'\partial_\chi\psi'}\right)= {\vert A(t) \vert}^{2} \textrm{i}k \partial_z \left( \hat{\varPsi} \hat{b}^{*} - \hat{\varPsi}^{*} \hat{b} \right) \equiv {\vert A(t) \vert}^{2} \mathrm{RS}_b, \end{gather}

where $\mathrm {RS}_u$ and $\mathrm {RS}_b$, defined above, have been introduced for brevity of notation and an asterisk denotes complex conjugation. Thus, the QL system (3.1)–(3.4) reduces to the hybrid slow-initial-value/fast-eigenvalue problem:

(3.9)\begin{gather} \partial_t \bar{u}= {\vert A(t) \vert}^{2} \mathrm{RS}_u +\frac{1}{Re_b}\partial_z^{2}\bar{u}+\bar{{f}}, \end{gather}

(3.10)\begin{gather}\partial_t\bar{b}= {\vert A(t) \vert}^{2} \mathrm{RS}_b +\frac{1}{Pr Re_b}\partial_z^{2}\bar{b}, \end{gather}

(3.11)\begin{gather}(\sigma + \textrm{i}k \bar{u} )(\partial_z^{2} - k^{2}) \hat{\varPsi} =\textrm{i}k(\partial_z^{2}\bar{u}) \hat{\varPsi}-\textrm{i}k \hat{b}+\frac{{Fr}}{Re_b}(\partial_z^{2} -k^{2})^{2} \hat{\varPsi}, \end{gather}

(3.12)\begin{gather}(\sigma +\textrm{i}k \bar{u}) \hat{b}=\textrm{i}k(1 + \partial_z\bar{b})\hat{\varPsi}+\frac{{Fr}}{Pr Re_b}(\partial_z^{2} -k^{2})\hat{b}, \end{gather}

where ${\vert A(t) \vert }=0$ if $\sigma _r < 0$.

The central challenge is to determine self-consistently an equation for the a priori unknown amplitude function ${\vert A(t) \vert }$, which, of course, is not constrained via the solution of the linear homogeneous eigensystem (3.11)–(3.12). As demonstrated by Michel & Chini (Reference Michel and Chini2019), to derive this constraint on ${\vert A(t) \vert }$ in the case of marginal stability, i.e. when $\sigma _r=0$, the linear eigensystem itself must be differentiated with respect to the slow time variable and an appropriate solvability condition enforced. This slow-time differentiation enables nonlinearity present in the coupled system (3.9)–(3.12) to be incorporated, as required for determination of the fluctuation amplitude.

To facilitate the required analysis, the linear subsystem (3.11)–(3.12) first is recast in the form $\mathcal {L}X = 0$, with $X = [\hat {\varPsi }(z),\hat {b}(z)]^\textrm {T}$ and

(3.13)\begin{equation} \mathcal{L} = \begin{pmatrix} (\sigma +\textrm{i}k \bar{u})(\partial_z^{2}-k^{2}) -\textrm{i}k \partial_z^{2} \bar{u}- \dfrac{Fr}{Re_b}(\partial_z^{2}-k^{2})^{2} & \textrm{i}k\\ -\textrm{i}k(1+\partial_z \bar{b}) & \sigma+\textrm{i}k\bar{u}- \dfrac{Fr}{Pr Re_b}(\partial_z -k^{2}) \end{pmatrix}, \end{equation}

supplemented with periodic boundary conditions in $z$. With respect to the eigenproduct

(3.14)\begin{equation} (X_1 \vert X_2) \equiv \int_0^{l_z} X_1(z) X_2^{*}(z)\,\mathrm{d}z\quad \forall\,(X_1,X_2), \end{equation}

the adjoint operator of $\mathcal {L}$, defined via $(\mathcal {L}X_1\vert X_2) = (X_1\vert \mathcal {L}^{{\dagger} } X_2)$, is given by

(3.15)\begin{equation} \mathcal{L}^{{{\dagger}}} = \begin{pmatrix} (\sigma^{*} \!-\!\textrm{i}k \bar{u})(\partial_z^{2}\!-\!k^{2}) -2\textrm{i}k \partial_z \bar{u} \partial_z\!-\! \dfrac{Fr}{Re_b}(\partial_z^{2}\!-\!k^{2})^{2} & \textrm{i}k(1+\partial_z \bar{b})\\ -\textrm{i}k & \sigma^{*}-\textrm{i}k\bar{u}- \dfrac{Fr}{Pr Re_b}(\partial_z -k^{2}) \end{pmatrix}. \end{equation}

Note that $\mathcal {L}$ is not self-adjoint ($\mathcal {L}\neq \mathcal {L}^{{\dagger} }$). The linear problem $\mathcal {L}X = 0$ is then differentiated with respect to the slow time $t$, yielding

(3.16)\begin{equation} \mathcal{L} \frac{\mathrm{d} X}{\mathrm{d} t}={-}\frac{\mathrm{d} \mathcal{L}}{\mathrm{d} t} X,\end{equation}

where

(3.17)\begin{equation} \frac{\mathrm{d} \mathcal{L}}{\mathrm{d} t}= \begin{pmatrix} (\mathrm{d}\sigma / \mathrm{d}t + \textrm{i}k \partial_t \bar{u}) (\partial_z^{2}-k^{2})-\textrm{i}k (\partial_z^{2} \partial_t \bar{u}) & 0\\ -\textrm{i}k \partial_z\partial_t\bar{b} & \mathrm{d}\sigma/\mathrm{d}t +\textrm{i}k\partial_t \bar{u}\end{pmatrix} {+ \frac{\mathrm{d}k}{\mathrm{d}t} M} \end{equation}

and $M$ is a matrix whose explicit computation turns out not to be necessary. Because $\mathcal {L}$ is singular, the linear system (3.16) is solvable if and only if the right-hand side is orthogonal to the corresponding null adjoint eigenvector $X^{{\dagger} }$ (the Fredholm alternative condition). This requirement can be readily confirmed by forming the inner product of the left-hand side of (3.16) with respect to $X^{{\dagger} }$, i.e.

(3.18)\begin{equation} \left. \left( \mathcal{L} \frac{\mathrm{d} X}{\mathrm{d} t} \right\vert X^{{{\dagger}}} \right) = \left.\left(\frac{\mathrm{d} X}{\mathrm{d} t} \right\vert \mathcal{L}^{{{\dagger}}} X^{{{\dagger}}} \right) = \left.\left( \frac{\mathrm{d} X}{\mathrm{d} t} \right\vert 0 \right) = 0. \end{equation}

Consequently, using (3.16), we obtain

(3.19)\begin{equation} \left.\left( \frac{\mathrm{d} \mathcal{L}}{\mathrm{d} t} X \right\vert X^{{{\dagger}}} \right)=0.\end{equation}

Crucially, this constraint can be made explicit. Using the mean-field equations (3.9) and (3.10) and noting that $X^{{\dagger} }=[\hat {\varPsi }^{{\dagger} }(z), \hat {b}^{{\dagger} }(z)]^\textrm {T}$, the solvability condition becomes

(3.20)\begin{equation} C_1 \frac{\mathrm{d}\sigma}{\mathrm{d}t} = C_2 {\vert A(t) \vert}^{2}+C_3+C_4 \frac{\mathrm{d}k}{\mathrm{d}t}. \end{equation}

As for the matrix $M$, the computation of the coefficient $C_4$ proves unnecessary. The remaining coefficients are given by the following expressions:

(3.21)\begin{gather} C_1 = \frac{1}{\textrm{i}k}\int_0^{l_z}\left[\hat{\varPsi}^{{{\dagger}} *} (\partial_z^{2} -k^{2}) \hat{\varPsi} + \hat{b}^{{{\dagger}} *} \hat{b} \right]\,\mathrm{d}z, \end{gather}

(3.22)\begin{gather}C_2 = \int_0^{l_z}\left\{ \mathrm{RS}_u \left[\hat{\varPsi}(\partial_z^{2} +k^{2}) \hat{\varPsi}^{{{\dagger}} *} + 2 \partial_z \hat{\varPsi} \partial_z \hat{\varPsi}^{{{\dagger}} *} - \hat{b}\hat{b}^{{{\dagger}} *} \right] - \mathrm{RS}_b \left[\hat{b}^{{{\dagger}} *} \partial_z \hat{\varPsi} + \hat{\varPsi} \partial_z \hat{b}^{{{\dagger}} *} \right] \right\} \,\mathrm{d}z, \end{gather}

(3.23)\begin{gather}C_3 = \int_0^{l_z}\left\{ \left( \bar{f} + \frac{\partial_z^{2} \bar{u}}{Re_b} \right) \left[\hat{\varPsi}(\partial_z^{2} + k^{2}) \hat{\varPsi}^{{{\dagger}} *} + 2 \partial_z \hat{\varPsi} \partial_z \hat{\varPsi}^{{{\dagger}} *} - \hat{b} \hat{b}^{{{\dagger}} *} \right]\right.\nonumber\\ \left. \quad - \frac{\partial_z^{2} \bar{b}}{Pr Re_b} \left( \hat{b}^{{{\dagger}} *} \partial_z \hat{\varPsi} + \hat{\varPsi} \partial_z \hat{b}^{{{\dagger}} *} \right)\right\} \,\mathrm{d}z. \end{gather}

Dividing both sides of (3.20) by $C_1$, the evolution of the linear growth rate $\sigma _r$ therefore is given for fixed $k$ by

(3.24)\begin{equation} { \left( \frac{\partial\sigma_r}{\partial t} \right)_{k} }= \alpha_r - \beta_r {\vert A(t) \vert}^{2},\end{equation}

where both $\alpha _r = \mathrm {Re}\{C_3/C_1\}$ and $\beta _r = \mathrm {Re}\{-C_2/C_1\}$ are computable functionals of the slow fields $\bar {u}$ and $\bar {b}$, the forcing $\bar {f}$, the direct and adjoint eigenfunctions $[\hat {\varPsi }, \hat {b}]^\textrm {T}$ and $[\hat {\varPsi }^{{\dagger} }, \hat {b}^{{\dagger} }]^\textrm {T}$, respectively, the parameters $Re_b$ and $Pr$ and the wavenumber $k$. Note that the total temporal variation of the real part of the growth rate $\sigma (\bar {\boldsymbol {G}},k)$,

(3.25)\begin{equation} \frac{\mathrm{d}\sigma_r}{\mathrm{d}t} = \left( \frac{\partial\sigma_r}{\partial t} \right)_{k} + \frac{\mathrm{d}k}{\mathrm{d}t} \left(\frac{\partial\sigma_r}{\partial k} \right)_{\bar{\boldsymbol{G}}}, \end{equation}

includes an additional contribution owing to the time variation of $k$. The second term on the right-hand side of this expression vanishes, however, provided that $k(t)$ is (locally) the most unstable mode, as required here. As discussed in detail by Michel & Chini (Reference Michel and Chini2019), temporal scale separation then requires (cf. (3.24)) that

(3.26)\begin{equation} {\vert A(t) \vert} = \left\{ \begin{array}{@{}ll} \sqrt{\alpha_r / \beta_r} & \text{if}\ \sigma_r=0, \ \alpha_r >0 \ \mathrm{and}\ \beta_r>0;\\ 0 & \text{otherwise}. \end{array} \right. \end{equation}

In contrast to its magnitude, the phase of $A(t)$ remains unconstrained in this formalism and, in fact, is not needed to evolve the reduced dynamics. Indeed, the reduced system (3.9)–(3.12) is invariant under any such phase change. In the fully nonlinear Boussinesq equations, this quantity would be fixed by the spatial phases of those fluctuation modes of generally small initial amplitude that become linearly unstable on the fast time scale; that is, this phase information ultimately depends on the details of the initial conditions, which are filtered in this framework. A primary virtue of this formalism is precisely that numerical time-integration need only be performed on the slow time scale $t$ without the arbitrary reinitialization of the fluctuation fields at each slow time instant. Specifically, at each slow time step, the following algorithm is executed.

1. (i) The direct eigenvalue problem $\mathcal {L}X=0$ is solved with periodic boundary conditions in $z$, and both the eigenvalue $\sigma$ of largest real part and the corresponding eigenvector $X$ are computed.

2. (ii) Step (i) is repeated over adjacent $k$ (i.e. $k- \Delta k$ and $k+\Delta k$, for suitably small horizontal wavenumber increment $\Delta k$) to reach a wavenumber $k$ corresponding to a local maximum of $\sigma _r(k)$.

3. (iii) If $\sigma _r(k)=0$ (numerically, if $\sigma _r(k)\geqslant 0$), the adjoint eigenvalue problem $\mathcal {L}^{{\dagger} } X^{{\dagger} } = 0$ is solved with periodic boundary conditions in $z$, and a non-zero eigenvector $X^{{\dagger} }$ corresponding to the null eigenvalue is returned. The modal amplitude ${\vert A(t) \vert }$ is then computed through (3.26).

4. (iv) The mean fields $\bar {u}$ and $\bar {b}$ are then time-advanced using a suitable discretization of (3.9)–(3.10).

Finally, note that in the present formulation, $k$ is determined by scanning over a range of wavenumbers at each slow time (step (ii) of the algorithm detailed above) to identify a local maximum of $\sigma _r(k)$. In companion work, we are seeking an extension of this algorithm that obviates the need to solve the eigenvalue problem repeatedly for different wavenumbers $k$ by deriving a set of slow equations for both ${\vert A(t) \vert }$ and $\mathrm {d}k/\mathrm {d}t$; see Ferraro (Reference Ferraro2019). Regardless of the algorithmic details, the wavenumber $k$ in the multiple time-scale formulation can smoothly vary on the slow time scale $t$. Crucially, this variation makes possible the accurate computation of the wavenumber that would be realized in the limit of an infinite horizontal domain, thereby ameliorating a recurring issue in the study of, e.g. the nonlinear dynamics of parallel shear flows and transition to turbulence (Tuckerman, Chantry & Barkley Reference Tuckerman, Chantry and Barkley2020).

4.1. Linear stability analysis

We begin by considering the linear stability of $\bar {u}_L(z) = \cos (3z)$ in the presence of the imposed linear stratification (and zero mean perturbation, so that $\bar {b}_L(z)=0$), as governed by the subsystem (3.3)–(3.4). Making a normal mode ansatz for a disturbance mode with $\chi$-wavenumber $k$, $\psi '(\chi,z,\tau )=\hat {\psi }(z)\mbox {e}^{i k(\chi -c\tau )}+\mathrm {c.c.}$, where the complex phase speed $c\equiv \omega /k$ (for angular frequency $\omega$), yields the TG equation in the non-dissipative limit $Fr/Re_b\to 0$:

(4.1)\begin{equation} \frac{\textrm{d}^{2}\hat{\psi}}{\textrm{d}z^{2}} +\left[\frac{\left(1+\partial_z\bar{b}_L\right)}{\left(\bar{u}_L-c\right)^{2}} -\frac{\partial_z^{2}\bar{u}_L}{\left(\bar{u}_L-c\right)}-k^{2}\right]\hat{\psi}=0. \end{equation}

The local-in-time gradient Richardson number $Ri_g=(1+\partial _z\bar {b})/|\partial _z\bar {u}|^{2}$. For the given basic-state profile, $\bar {u}_L(z) = \cos (3z)$ and $\bar {b}_L(z)=0$, so the gradient Richardson number of the laminar state

(4.2)\begin{equation} Ri_{gL}=\frac{1}{m^{2}\sin^{2}(mz)}, \end{equation}

which attains a minimum value $Ri_{gL_{{min}}}=1/m^{2}$ for $z=\pm {\rm \pi}/(2m),\pm 3{\rm \pi} /(2m) \ldots$. Thus, for the chosen parameter values, $Ri_{gL_{{min}}}=1/9<1/4$, so the Miles–Howard criterion necessary for linear instability is satisfied. Of course, the incorporation of diffusion will modify this inviscid criterion, but the effective diffusivity $Fr/Re_b$ is much less than unity in the subsequent numerical simulations and acts only as a regular perturbation to the non-dissipative problem in smooth regions of the flow.

4.2. Nonlinear evolution

To assess the performance of the reduced system quantitatively, a set of numerical simulations is performed for $Re_b=1$ and decreasing values of $Fr$ down to $10^{-2}$. Three different numerical algorithms are compared, each implemented in the computing environment Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). In each case, a pseudo-spectral method is used with a second-order Runge–Kutta time-stepping scheme. Python codes are provided as supplementary material available at https://doi.org/10.1017/jfm.2021.1060.

First, a set of DNSs of the primitive 2-D Boussinesq equations, that is (2.1)–(2.4) rendered two-dimensional, is performed. With $\alpha =Fr$, we note that the dimensionless total energy density is $(u^{2}+Fr^{2} w^{2}+b^{2})/2$. The domain, of vertical size $l_z=2{\rm \pi} /3$ and of horizontal size $L_x=(2{\rm \pi} /k)Fr$, is discretized using a Fourier–Fourier pseudo-spectral method with 128 grid points in each direction. Note that $k$ must be specified initially and remains fixed for the duration of the simulation.

The second algorithm is the single time-scale formulation of the QL system (STQL) introduced in § 3.1. These simulations also are set in a domain of size $l_z=2{\rm \pi} /3$ and $l_x=2{\rm \pi} /k$ using $128$ grid points (cf. Fourier modes) in each direction. As for the DNSs, $k$ must be fixed a priori. Unlike the DNSs, certain terms in the governing Boussinesq equations are consistently neglected in the STQL simulations, in accord with the multiple scales analysis performed in the limit $Fr \rightarrow 0$ with $Re_b$ fixed.

The third numerical scheme integrates the multiple time-scale formulation of the QL system (MTQL). The equations governing the 1-D evolution of the slowly evolving mean fields $\bar {u}(z,t)$ and $\bar {b}(z,t)$, (3.1) and (3.2), respectively, are discretized using Chebyshev polynomials – rather than Fourier modes, for compatibility with the Dedalus eigenvalue solver – with 128 grid points. As detailed in § 3.2, the filtered dynamics of the fast fluctuations is obtained by solving at each slow time step a 1-D eigenvalue problem, which yields the instantaneous linear growth rate $\sigma _r$ and vertical mode structure ($\hat {\varPsi }(z;t$) and $\hat {b}(z;t)$) of the most unstable mode. Once the wavenumber $k$ is updated and determined to correspond to a local maximum of $\sigma _r(k)$ within $\Delta k=10^{-3}$, the amplitude of the fluctuation mode is set in accord with condition (3.26). Again, we emphasize that in the MTQL algorithm the wavenumber $k$ is not fixed, but evolves so that the marginally stable wavenumber is tracked, with modes corresponding to adjacent values of $k$ being linearly stable. In principle, the amplitude $A(t)$ of the fluctuation mode would be updated continually such that $\sigma _r = 0$ remains fixed; in practice, $\sigma _r$ would remain very close to the first positive value computed, which depends on the specific value of the time step $\Delta t$ but, in any case, will be very small compared with unity (e.g. $\sigma _r = 2\times 10^{-4}$ for $Fr=0.02$, $Re_b=1$ and $\Delta t=10^{-4}$). To achieve a steady state in which $\sigma _r$ is independent of the time step, the modal amplitude $A(t)$ is artificially increased during the transient regime such that

(4.3)\begin{equation} \text{if}\ \sigma_r > 10^{{-}6} \text{ and } A>0, \text{ then } A^{2} =\frac{\alpha_r}{\beta_r} + \frac{\sigma_r}{1000 \beta_r {\Delta t}}. \end{equation}

According to (3.24), $\sigma _r$ will, if $A$ is subject to this change, vary during the ($n+1$)th time step according to $\Delta \sigma _r\equiv \sigma _r^{n+1}-\sigma _r^{n}=-\sigma _r^{n}/1000$, and therefore will relax exponentially to zero provided the stated conditions are satisfied. Note that this modification minimally affects the transient regime and is no longer applied once a steady state is reached, because then $\sigma _r < 10^{-6}$.

The same time step $\Delta t = 10^{-4}$ is used in all of the numerical simulations. It should be emphasized, however, that because the MTQL system is evolved on the slow time scale, much larger values of $\Delta t$ could be used for the MTQL simulation. Although a detailed study of the computational savings afforded by the MTQL algorithm is beyond the scope of this investigation, we note here that the steady state depicted subsequently for $Fr =0.02$ and $Re_b=1$ is also quantitatively reproduced with the MTQL algorithm using a time step $\Delta t = 0.01$ (i.e. 100 times larger) for which the DNS algorithm is numerically unstable.

We first discuss the MTQL dynamics for $Re_b=1$ and $Fr=0.02$. The simulation is forced from a rest state, i.e. with zero initial velocity and zero buoyancy deviation from the imposed background linear stratification. During the transient phase of the dynamics evident in figure 2(a), the real part of the maximum growth rate over all horizontal wavenumbers greater than a minimum viscous threshold (see below) remains negative until $t \simeq 0.175$, when it reaches zero: the amplitude of the fluctuation fields $A$ then jumps to a finite value to prevent further growth of $\sigma _r$. This behaviour, which is clearly depicted in the movie included as supplementary movie 1, accounts for the discontinuity in the slow-time evolution of the total energy (see inset in figure 2a). A nonlinear steady state – one type of ECS arising in forced dissipative nonlinear infinite-dimensional dynamical systems (Kawahara & Kida Reference Kawahara and Kida2001; Lucas & Caulfield Reference Lucas and Caulfield2017; Lucas et al. Reference Lucas, Caulfield and Kerswell2017; Parker, Caulfield & Kerswell Reference Parker, Caulfield and Kerswell2019) – supported by a single horizontal mode with a wavenumber $k = 2.515$ eventually is attained. Note that this single-mode structure is emergent, not imposed.

Figure 2. DNS, STQL and MTQL simulations of stratified Kolmogorov flow forced from a randomly perturbed rest state with $Re_b = 1$ and $Fr=0.02$. (a) Total energy of the flow as a function of time. For the MTQL simulation, $\sigma _r$ also is plotted. (b) Maximum fluctuation growth rate $\sigma _r$ as a function of $k$ shown for three different times during the MTQL simulation. Note the saturation of $\sigma _r$ to zero at wavenumbers $k\simeq 2.5$ and the presence of weakly amplified oscillatory modes (ignored in this study), with non-zero $\sigma _i$, at small $k$.

Figure 2(b) shows the evolution of the growth-rate spectrum during the MTQL simulation. At $t=1.5$, the earliest time depicted, all modes with wavenumbers greater than approximately unity are damped. We have confirmed that the linearly unstable oscillatory modes (with $\sigma _i\ne 0$) arising at small $k$ disappear for vanishingly small stratification; in contrast, stratification is not a necessary condition for the instabilities evident at later times for $k\simeq 2.5$, a distinction recently made by Parker et al. (Reference Parker, Caulfield and Kerswell2019). For simplicity, we choose to avoid constraining the unstable modes arising at small $k$ because $\sigma _r$ remains very small – and is even smaller at larger values of $Re_b$ (see the inset in figure 8) – for these modes when the mode with ${O}(1)$ $k$ is saturated. Nevertheless, we emphasize that, in principle, the MTQL algorithm could be modified to marginalize these modes too. At an intermediate time, $t=1.75$, a marginal mode with wavenumber $k\simeq 2.25$ corresponds to the mode with locally maximum growth rate, but by $t=2.5$, the wavenumber of this mode has evolved to $k\simeq 2.5$. For the given configuration, these modes with ${O}(1)$ $k$ have zero phase speed (because $\sigma _i=0$) and correspond to Kelvin–Helmholtz instabilities.

To perform a quantitative comparison of the three algorithms, both a DNS and STQL simulation also are run for the same parameters ($Re_b=1$ and $Fr=0.02$). For these simulations, the horizontal domain length, $L_x=(2{\rm \pi} /k)Fr$ and $l_x=2{\rm \pi} /k$, respectively, is specified using $k=2.515$, i.e. the emergent steady-state wavenumber obtained using the MTQL algorithm. Were it possible to perform the DNS and STQL simulation in the limit $L_x, l_x \rightarrow \infty$, we would expect an ECS with this horizontal wavenumber to be realized. The evolution of the total energy reported in figure 2 shows that, in all three simulations, there is an overshoot of the total energy density, corresponding to the ‘bursting regime’ documented by Michel & Chini (Reference Michel and Chini2019) and observed in 3-D DNS of the full Boussinesq equations (Rorai et al. Reference Rorai, Minini and Pouquet2014; Feraco et al. Reference Feraco, Marino, Pumir, Primavera, Mininni, Pouquet and Rosenberg2018). The MTQL and STQL simulations converge precisely to the same steady state, as can be confirmed by inspection of figures 2 and 3; of course, this convergence is not entirely surprising (although not guaranteed) given that the same equations are being solved in the steady limit (cf. (3.1)–(3.4)).

Figure 3. Comparison of the steady ECSs computed in stratified Kolmogorov flow for $Re_b = 1$ and $Fr=0.02$ using DNS (a), MTQL (b) and STQL (c) algorithms. The total buoyancy field (including the imposed linear profile) is shown in colour, while arrows are used to depict the isotropically scaled velocity field in each case. To accentuate differences between the QL results and DNS, the total buoyancy contours $z+\bar {b}=\pm 0.5$ are highlighted (dashed curves).

In figure 3, which shows the steady ECSs computed using the DNS (a), MTQL (b) and STQL (c) algorithms, colour indicates the total buoyancy field, i.e. including the imposed background linear stratification, while arrows are used to depict the velocity field. To facilitate the comparison of the 2-D structure, note that the DNS velocity field $(u, \epsilon ^{2} W)$ is plotted, while the corresponding fields $(\bar {u}_0+ \epsilon u_1' , \epsilon W_{-1}')$ are shown for the STQL and MTQL simulations (cf. (2.6) and (2.7)); that is, all velocity components are normalized by $U$. The agreement among the three steady-state ECSs lends confidence to the asymptotically reduced system (3.1)–(3.4) and thereby to the novel MTQL algorithm. The small quantitative discrepancy between the final energies of the ECS obtained from the DNS and via the STQL/MTQL algorithms (figure 2a) is partly attributable to the omission of the mean-field correction $\epsilon \bar {u}_1$ from the latter schemes. If desired, this correction could be self-consistently computed by carrying the analysis to higher-order, although we emphasize that the reduced system (3.1)–(3.4) is both asymptotically consistent (apart from the diffusive regularization) and closed.

The agreement between the DNS and MTQL/STQL simulations is expected to improve as $Fr$ is decreased, because the reduced model is derived in the limit $Fr \rightarrow 0$. To assess the asymptotic convergence of the MTQL algorithm in that limit, a suite of DNSs and MTQL simulations with $Re_b=1$ is performed for a set of decreasing values of $Fr$. For the given parameter regime, the MTQL algorithm always converges to a steady state. The long-time dynamics exhibited by the DNS, however, becomes time-dependent for $Fr< 0.012$. Given that the fluctuation-induced mixing is of central importance (see § 4.3), we choose to compare, in figure 4, the energy of the fluctuation fields computed using DNS and the MTQL algorithm for $Fr \geqslant 0.012$. (In particular, the contribution of the fluctuations to the total energy is ${O}(Fr)$, and hence we use this proxy to examine quantitatively the limit $Fr \rightarrow 0$). Indeed, using Parseval's identity, the dimensionless total energy of the primitive 2-D Boussinesq equations (with $\alpha =Fr$),

(4.4)\begin{equation} E_{{tot}} = \int\,\mathrm{d}\chi\,\mathrm{d}z \left(u^{2} + Fr^{2} w^{2} + b^{2}\right)/2, \end{equation}

can be re-expressed as a sum of the energy contained in the various horizontal modes $k$; that is, $E_{{tot}} = \sum _{k\ge 0}E_k$.

Figure 4. Convergence of the steady ECSs computed at $Re_b=1$ using the asymptotically reduced system (3.1)–(3.4) to the corresponding steady solution of the full 2-D Boussinesq equations in the limit $Fr \rightarrow 0$. The red symbols show the relative error computed using the total fluctuation energy obtained from the DNS (the sum of the energy in each non-zero horizontal Fourier mode), i.e. $r_{k>0}=|E_{k>0}^{{{DNS}}}-E_{k>0}^{{{MTQL}}}|/E_{k>0}^{{{DNS}}}$. The blue symbols show the relative error evaluated using the fluctuation energy only in the Fourier mode with the fundamental wavenumber $k$, again as obtained from the DNS, $r_{k}=|E_{k}^{{{DNS}}}-E_{k}^{{{MTQL}}}|/E_{k}^{{{DNS}}}$. For reference, the solid green line has a slope equal to $Fr^{1}$.

As evident in the figure, the relative error in the fluctuation energy is computed two different ways. The red circles show the relative error when the fluctuation energy is computed using the total DNS fluctuation energy summed across all horizontal Fourier modes with non-zero wavenumber. In addition, the blue squares show the relative error based on the DNS fluctuation energy contained only in the fundamental Fourier mode with wavenumber $k$. Clearly, both metrics confirm that solutions of the asymptotically reduced equations converge to those of the full Boussinesq equations in the distinguished limit considered, namely, as $Fr\to 0$ with $Re_b$ fixed, although only the data for the fundamental mode approximates the expected rate of decay in the relative error.

4.3. Exact coherent states for larger $Re_b$

Steady states can be reached in the MTQL simulations at $Re_b=1$ even for $Fr < 0.012$, a regime in which the DNS exhibit persistent unsteady dynamics. Nevertheless, the long-time MTQL dynamics need not be steady or regular on the slow time scale, and, for sufficiently large $Re_b$, we anticipate that more complex dynamics would be exhibited within the two time-scale framework. For the parameter regime considered here, however, the MTQL algorithm – although ostensibly not specifically designed for this purpose – has the virtue of yielding ECSs that would be linearly unstable within the full Boussinesq dynamics. For at least the last twenty years, such states have been known to be able to capture certain characteristic attributes of the turbulent regime; see e.g. Kawahara & Kida (Reference Kawahara and Kida2001). ECSs have been of particular interest in the context of wall-bounded constant-density shear flows, where a general mechanism supporting such states has been identified (Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995). Numerical computation of ECSs from the governing Navier–Stokes or Boussinesq equations typically requires recurrent flow analyses of expensive DNSs to provide suitable initial conditions for sophisticated Newton-hookstep solvers, although more efficient approaches have been proposed (Page & Kerswell Reference Page and Kerswell2020). In contrast, it is clear from the present investigation and related studies (Hall & Sherwin Reference Hall and Sherwin2010; Beaume et al. Reference Beaume, Chini, Julien and Knobloch2015; Montemuro et al. Reference Montemuro, White, Klewicki and Chini2020) that asymptotically reduced systems, derived using multiple scales analysis, retain the dominant interactions that sustain such states while self-consistently filtering other dynamics, yielding more efficient algorithms for ECS computations and simultaneously exposing the underlying physical mechanisms.

Lucas et al. (Reference Lucas, Caulfield and Kerswell2017) and Lucas & Caulfield (Reference Lucas and Caulfield2017) computed ECSs in 3-D stratified Kolmogorov flow with a horizontally varying forcing (in contrast to the present study). Of particular interest is the mixing efficiency associated with these ECSs, which the authors were able to compute for $Re_b \le 95$ and $Fr \ge 0.23$. With a similar aim, we compute by continuation steady-state ECSs in 2-D stratified Kolmogorov flow using the MTQL algorithm for $Fr=0.01$ and $Re_b\in [1,10]$; i.e. steady states obtained at given $Re_b$ values are used as initial iterates in MTQL simulations performed at incrementally larger $Re_b$. We again emphasize that the horizontal domain size (or fundamental wavenumber) is not imposed a priori but rather is an emergent property of our computations.

The steady ECS fields for the two extreme values of $Re_b$ are depicted in figure 5. Even these relatively simple (2-D) steady states exhibit features that are commonly associated with stratified turbulence; in particular, the solutions indicate the spontaneous emergence of two layers of relatively well-mixed fluid centred on $z = {\rm \pi}/6 \simeq 0.5$ and $z = {\rm \pi}/2 \simeq 1.6$, separated by broader ‘streams’ flowing in opposite directions. This feature is even more clearly evident in the mean buoyancy and velocity profiles ($\bar {b}$ and $\bar {u}$) shown in figure 6 and can be interpreted in the light of the linear stability analysis described in § 4.1. The vertical locations of the mixed layers match those of the minima of the gradient Richardson number of the laminar flow; see (4.2). In the inviscid limit, these locations would correspond to critical layers in which $\bar {u}_L-c$ vanishes, potentially leading to discontinuities in both the fluctuation buoyancy and horizontal velocity fields. In practice, the small but finite viscosity smooths these abrupt variations. In the small $Fr$ limit analysed here, the system self-adjusts to an equilibrium state that is approximately marginally stable with respect to the emergent mean profiles – so the critical layer interpretation remains quantitatively valid provided that $\bar {u}_L$ is replaced with $\bar {u}$, and $z$ (the background stratification) is replaced with $z+\bar {b}$.

Figure 5. Comparison of steady ECSs in stratified Kolmogorov flow for $Fr=0.01$ and $Re_b = 1$ (a) and ${Re_b=10}$ (b) computed using the MTQL algorithm. The total buoyancy field (including the imposed linear profile) is shown in colour, while the arrows are used to depict the isotropically scaled velocity field in each case.

Figure 6. Steady-state mean vertical profiles of buoyancy (a) and horizontal velocity (b) computed using the MTQL algorithm for $Fr=0.01$ and $Re_b=1$ (blue) and $Re_b=10$ (red). The dashed lines correspond to the unstable laminar state, with velocity profile $\bar {u}_L=\cos 3z$ and buoyancy perturbation $\bar {b}=0$.

The staircase-like profile of $z+b$ has been analysed thoroughly in the recent 3-D DNS study of Maffioli (Reference Maffioli2019), who measured the dimensionless third-order moment (i.e. the skewness),

(4.5)\begin{equation} S \equiv \frac{\left\langle \left(\partial_z b\right)^{3} \right\rangle}{\left\langle \left(\partial_z b\right)^{2} \right\rangle^{3/2}}, \end{equation}

where $\langle {\cdot } \rangle$ denotes a spatial average over the entire domain, and it should be recalled that $b(\chi,z) =\bar {b}(z)+b'(\chi,z)$ is the buoyancy deviation from the imposed linear profile ($z$). Interestingly, based on the DNS of Maffioli (Reference Maffioli2019), $S$ does not appear to converge to a finite value in the strongly stratified turbulence regime. Instead, the author finds that $S \sim c Fr^{-0.41}$ (for some constant $c$), the dependence on $Re_b$ not having been investigated. Noting that the layering evident in figure 6 is related directly to $\bar {b}$ rather than to $b$, we therefore also define the skewness of the horizontally averaged buoyancy profile

(4.6)\begin{equation} \bar{S} = \frac{\left\langle \left( \partial_z \bar{b}\right)^{3} \right\rangle}{\left\langle \left(\partial_z \bar{b}\right)^{2} \right\rangle^{3/2}}. \end{equation}

Indeed, given the expansion (2.6), $S$ reduces to $\bar {S}$ at leading order. The averaging over $\chi$ is, of course, straightforward in the present study owing to the multiple scales decomposition but would prove challenging in a DNS because of the spatial variation of the layer with $\chi$. As shown in figure 7, $\bar {S}$ increases monotonically with $Re_b$ with no sign of convergence to a finite value as $Re_b$ is increased. A crude piecewise linear approximation of the buoyancy profiles yields an estimate of $\bar {S}$ as a function of the ratio of the height of the streams $h_{st} \sim L \times Fr$ to that of the emergent mixed layer $h_{ML}$ only (Maffioli Reference Maffioli2019). This estimate suggests that $h_{ML}/L$ does not become independent of $Re_b$ in the strongly stratified limit; in particular, $h_{ML}$ is not simply given by the Ozmidov scale $L_O \sim L \times Fr^{3/2}$.

Figure 7. Mean quantities computed from the steady-state ECSs obtained at $Fr=0.01$ as a function of $Re_b$: (a) skewnesses $\bar {S}$; (b) flux Richardson number $Ri_F$.

Another quantity of interest is the mixing efficiency achieved by these steady ECSs. Characterizing stratified mixing in the natural environment is a longstanding question, notably because stratified turbulence in the atmosphere and oceans occurs in a parameter regime in which both $Fr \ll 1$ and $Re_b \gg 1$. As discussed in detail by Gregg et al. (Reference Gregg, D'Asaro, Riley and Kunze2018), an unfortunate additional source of complexity is that several different measures of mixing have been proposed, which can differ significantly from each other in the turbulent regime; see e.g. Venayagamoorthy & Koseff (Reference Venayagamoorthy and Koseff2016). Here, we choose to compute the flux Richardson number $Ri_F$, defined as the ratio of the volume-integrated buoyancy flux to the power injected by the external force:

(4.7)\begin{equation} Ri_F = \frac{\int \tilde{w}\,\tilde{b} \,\mathrm{d}V}{\int \tilde{\boldsymbol{f}} \boldsymbol{\cdot} \tilde{\boldsymbol{u}} \,\mathrm{d}V},\end{equation}

where tildes refer to dimensional variables. For steady ECSs, this measure coincides with other quantities usually introduced to quantify the mixing, e.g. the mixing efficiency $\mathcal {E}$ based on the dissipation rates of kinetic and available potential energy; see Caulfield & Peltier (Reference Caulfield and Peltier2000), Peltier & Caulfield (Reference Peltier and Caulfield2003), Salehipour & Peltier (Reference Salehipour and Peltier2015) and Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016). The DNSs of Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) and Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019) suggest a constant flux coefficient $\varGamma =Ri_F/(1-Ri_F)\simeq 0.2$ in the strongly stratified turbulent limit, whereas the periodic-orbit ECSs of Lucas & Caulfield (Reference Lucas and Caulfield2017) indicate a progressive decrease toward zero as $Re_b$ is increased. Note that the turbulent flux coefficient experiences strong variations as a function of $Fr$ in the range $[0.05, 0.5]$ (Feraco et al. Reference Feraco, Marino, Pumir, Primavera, Mininni, Pouquet and Rosenberg2018) and that in situ measurements in the open ocean indicate a transition from a constant value to a decreasing one at $Re_b \gtrsim 100$ (Lozovatsky & Fernando Reference Lozovatsky and Fernando2013; Monismith et al. Reference Monismith, Koseff and White2018) although there is still outstanding controversy; see Caulfield (Reference Caulfield2020, Reference Caulfield2021). DNSs in the regime $Fr \ll 1$ and $Re_b \gg 1$ clearly are needed to resolve this issue conclusively but are not feasible with the current computing power (figure 1).

Given the scalings assumed here, the flux Richardson number defined in (4.7) is computed for the steady-state ECSs via

(4.8)\begin{equation} Ri_F ={-} \textrm{i}k |A|^{2} \frac{\int_0^{l_z}(\hat{\varPsi}^{*} \hat{b} - \hat{\varPsi}\hat{b}^{*} ) \,\mathrm{d}z}{\int_0^{l_z}\bar{f}\,\bar{u}_0\, \mathrm{d}z}. \end{equation}

The variation of $Ri_F$ with $Re_b$ is shown in figure 7(b). For the range of $Re_b$ considered, $Ri_F$ varies between 0.13 and 0.18. These values correspond to values of the flux coefficient $\varGamma$ ranging from approximately 0.16 to 0.22, in reasonable accord with the upper bound on $Ri_F$ of 0.15 proposed by Osborn (Reference Osborn1980) and corresponding to $\varGamma \le 0.2$. The overshoot of the mixing measure is reminiscent of that reported by both Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) and Lucas & Caulfield (Reference Lucas and Caulfield2017). Unfortunately, as in other investigations, larger values of $Re_b$ would be required to establish the limiting behaviour of $Ri_F$ for these steady ECSs. One virtue of the systematically reduced formulation derived here is that it should facilitate such investigations. Although we leave this task for a future study, referring to figure 7, it is conceivable that $Ri_F$ asymptotes to a value near 0.17 for large $Re_b$ (again, for the given steady 2-D ECSs); intriguingly, this value corresponds to $\varGamma =0.2$, as seen by Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019).

Finally, we discuss the $Re_b$-dependence of the fluctuation growth-rate spectrum for the steady ECSs. As evident in figure 8, the positive linear growth rate regime at small wavenumber is suppressed as $Re_b$ is increased, partly justifying our choice to disregard these modes in this investigation because our reduced modelling efforts ultimately are aimed at accessing large values of $Re_b$. Interestingly, at $Re_b=10$, there is a second local maxima in the growth rate (ignoring the low-wavenumber mode) around $k\simeq 2.5$ that plausibly could spawn the emergence of a second linearly unstable mode for which marginal stability also would have to be enforced. This eventuality could be treated by introducing a second non-zero modal amplitude $B(T)$ and vertical eigenfunctions $[\hat {\varPsi }_B(z), \hat {b}_B(z)]$ corresponding to this wavenumber, labelled $k_B$, generalizing the procedure introduced in § 3.2: instead of the single equation $(\partial \sigma _{Ar}/\partial t)|_{k_A} =\alpha _{Ar} - \beta _{Ar} |A|^{2}$ for the growth rate $\sigma _{Ar}$ (with associated coefficients $\alpha _{Ar}$, $\beta _{Ar}$) of the sole marginal mode with amplitude $A$, coupled equations for both $\sigma _{Ar}$ and $\sigma _{Br}$ would be derived and solved simultaneously. In this manner, the formalism enables scale-selective adaptivity, in that modes representing newly emergent length scales are introduced only as needed rather than democratically. Again, this procedure can be implemented with no artificial quantization of wavenumbers imposed by the domain size.

Figure 8. Maximum fluctuation-field growth rate as a function of $k$ for the steady ECSs computed at ${Fr=0.01}$ and $Re_b=1$ (blue) and $Re_b=10$ (red). The inset shows the positive growth rates that arise for small $k$; comparing with figure 2(b) confirms that these low-wavenumber modes are viscous instabilities that vanish as $Re_b\to \infty$.

In fact, the potential emergence of a second marginal mode is not the reason we limited this investigation of steady ECSs to $Re_b\le 10$. Rather, for $Re_b>10$ with ${Fr=0.01}$, the steady states become dynamically unstable even within the MTQL algorithm. Ultimately, at a certain time during the evolution, $\sigma _r=0$, $\alpha _r>0$ and $\beta _r<0$. Physically, this situation corresponds to marginally stable fluctuations that would become linearly unstable (because $\alpha _r>0$) and whose feedback on the mean field $\{ \bar {u}, \bar {b} \}$ would render the mean field even more unstable (because $\beta _r<0$). In practice, the fluctuation amplitude would reach such large values that the mean field would be forced to respond on the fast time scale, invalidating the ansatz of time-scale separation (and possibly quasi-linearity). We emphasize, however, that preliminary studies confirm our expectation that this dynamic (a ‘burst’) persists only during a short transient, after which marginal stability is re-established and the MTQL algorithm can be restarted. This issue is described by Michel & Chini (Reference Michel and Chini2019), and some promising preliminary results involving modification of the MTQL algorithm to account for such bursting regimes are reported by Ferraro (Reference Ferraro2019). Notwithstanding this current limitation, we emphasize that the value $Re_b=10$ is attained with the MTQL algorithm for a realistically small Froude number, viz. $Fr = 0.01$. In their 3-D DNS, the smallest Froude number reached by Lucas & Caulfield (Reference Lucas and Caulfield2017) is $0.008$, for which $Re_b=0.25$, and that by Maffioli & Davidson (Reference Maffioli and Davidson2016) is $0.02$, for which $Re_b=17$. Although our study has been restricted to 2-D stratified Kolmogorov flow, the ECS computations have been performed on a laptop computer, demonstrating the potential for asymptotically reduced modelling of strongly stratified flows with as yet inaccessibly large $Re_b$ and small $Fr$. In contrast, the recent DNSs of the 2-D Boussinesq equations carried out by Kumar, Verma & Sukhatme (Reference Kumar, Verma and Sukhatme2017) required high-performance computing resources to achieve $Re_b=130$ (respectively $Re_b= 300$) for their smallest Froude numbers $Fr=0.16$ (respectively $Fr=0.31$), which suggests that only $O(1)$ values of $Re_b$ could be attained for $Fr=0.01$ with similar computational resources.