2.1. QL approximation

The QL approximation to (1.1), has as its starting point the Reynolds decomposition

(2.1)\begin{equation} \varphi = \bar{ \varphi } + \varphi', \end{equation}

where $\bar { \varphi }$ denotes the zonal mean of a quantity $\varphi$ and $\varphi '$ its non-zonal perturbation. Decomposing the vorticity equation (1.1), and neglecting the non-zonal terms given by $\boldsymbol {\nabla }\boldsymbol {\cdot } (\boldsymbol {u}' \zeta ')- \boldsymbol {\nabla } \boldsymbol {\cdot } (\bar {\boldsymbol {u}' \zeta '})$ in the perturbation equation, results in the QL equations

(2.2)$$\begin{gather} \bar{\zeta}_t ={-}(\overline{ v'\zeta' })_y -\mu \bar{ \zeta } + ({-}1)^{n+1}\nu_{2n}\partial_y^{2n}\bar{ \zeta } + \sqrt{\varepsilon} \bar{ \eta }, \end{gather}$$

(2.3)$$\begin{gather}\zeta'_t ={-} \bar{u} \zeta'_x - v'\bar{ \zeta }_y - \beta v' -\mu\zeta' +({-}1)^{n+1}\nu_{2n}\nabla^{2n}\zeta' + \sqrt{\varepsilon} \eta'. \end{gather}$$

By construction, the QL equations allow for nonlinear interactions only between pairs of wave modes and a single zonal mode, and not between triads of wave modes. This excludes the possibility of a spectrally local inverse cascade of energy, as in the classical Kraichnan–Batchelor 2D turbulence theory (Batchelor Reference Batchelor1967; Kraichnan Reference Kraichnan1967).

2.2. CE2 equations

Following, e.g., Srinivasan & Young (Reference Srinivasan and Young2012) and Marston et al. (Reference Marston, Conover and Schneider2008), the CE2 equations are obtained by introducing an averaging operator $\langle {\cdot } \rangle$, which in different interpretations (see, e.g., Ait-Chaalal et al. Reference Ait-Chaalal, Schneider, Meyer and Marston2016; Marston, Qi & Tobias Reference Marston, Qi and Tobias2019) is an ensemble average over the noise process $\eta$, or a zonal average. In both cases, $\langle {\cdot } \rangle$ satisfies the Reynolds decomposition property $\langle \varphi _1 \varphi _2 \rangle = \langle \varphi _1 \rangle \langle \varphi _2 \rangle - \langle \varphi _1' \varphi _2' \rangle$, where $\varphi '_i= \varphi _i -\langle \varphi _i \rangle$. The CE2 equations are most conveniently formulated for the averaged zonal velocity $U=\langle u \rangle$ and the second cumulant of vorticity

(2.4)\begin{equation} \mathcal{Z}(\boldsymbol{x}_1,\boldsymbol{x}_2,t)=\langle \zeta'(\boldsymbol{x}_1,t) \zeta'(\boldsymbol{x}_2,t) \rangle. \end{equation}

Applying the averaging operator to (1.1), neglecting the non-zonal terms as was done for (2.2)–(2.3), results in the CE2 equations

(2.5)\begin{gather} U_t ={-}\langle u'v' \rangle_y - \mu U - ({-}1)^n\nu_{2n}\partial_y^{2n}U, \end{gather}

(2.6) \begin{gather} \partial_t\mathcal{Z} + \left(U_1-U_2\right) \partial_x \mathcal{Z} ={-}( (\beta - U''_1) \nabla_1^{{-}2} - (\beta - U''_2) \nabla_2^{{-}2}) \partial_x \mathcal{Z}\nonumber\\ - 2\mu \mathcal{Z} - ({-}1)^n \nu_{2n}(\nabla_1^{2n} + \nabla_2^{2n}) \mathcal{Z} + \varepsilon \varPi. \end{gather}

Here $\varPi =\langle \eta '(\boldsymbol {x}_1,t) \eta '(\boldsymbol {x}_2,t) \rangle$ is the covariance of the stochastic process $\eta '$, which is assumed homogeneous so that $\varPi \equiv \varPi ( \boldsymbol {x}_1-\boldsymbol {x}_2,t)$. It then follows that, in (2.6), $\mathcal {Z}\equiv \mathcal {Z}(x,y_1,y_2)$ where $x=x_1-x_2$ is the relative zonal distance between the two arguments $\boldsymbol {x}_1$ and $\boldsymbol {x}_2$. Numerical subscripts in (2.6) refer to the argument variable, i.e. $\nabla _i$ acts on the $\boldsymbol {x}_i$ variable, $\nabla ^{-2}_i$ denotes the inverse of the Laplacian operator again applied to the $\boldsymbol {x}_i$ variable, and $U_i=U(y_i)$.

To close the equations it is necessary to evaluate the Reynolds stress term $\langle u'v' \rangle$ appearing in (2.5) in terms of $\mathcal {Z}$. From the relationship between relative vorticity $\zeta '$ and velocity $\boldsymbol {u}'$ the required result is

(2.7) \begin{equation} \langle u'v' \rangle(y) ={-} [\partial_{y_1}\partial_{x}\nabla^{{-}2}_1\nabla^{{-}2}_2 \mathcal{Z}(x, y_1,y_2)]_{x \to 0,~y_1\to y_2=y}. \end{equation}

The CE2 equations (2.5)–(2.7) are the simplest set in a hierarchy of cumulant models of the statistics of (1.1). A sequence of potentially more accurate, if typically more computationally expensive, cumulant models result from adopting different closures for the wave–wave interaction terms which are neglected in CE2 (see, e.g., Farrell & Ioannou (Reference Farrell and Ioannou2003), Marston et al. (Reference Marston, Qi and Tobias2019) and references therein). The CE2 equations are mathematically equivalent to the SSST equations originally derived and analysed in Farrell & Ioannou (Reference Farrell and Ioannou2003, Reference Farrell and Ioannou2007). In each case, the closure models prioritise non-local zonal energy exchanges, pitching them in opposition to turbulent cascade theories Rhines (Reference Rhines1975) and PV mixing theories (Dritschel & McIntyre Reference Dritschel and McIntyre2008; Dritschel & Scott Reference Dritschel and Scott2011).

2.3. The SY14 result for the momentum flux in a linear shear flow

Next, a re-working of the SY14 solution to (2.6), for the special case in which $U(y)=\gamma y$ is a steady constant shear flow, is presented. The aim is to obtain a relatively simple and easy to evaluate formula for the steady momentum flux $\langle u'v' \rangle$, obtained from the solution to (2.6) as $t \to \infty$, which will be used to understand the outcome of our numerical calculations in the following.

Our focus is restricted to a time-independent homogeneous forcing covariance $\varPi (\boldsymbol {x})$, which takes the argument $\boldsymbol {x}=(x,y)^T=\boldsymbol {x}_1-\boldsymbol {x}_2$, and its Fourier transform $\hat {\varPi }(\boldsymbol {k})$ defined by

(2.8a,b)\begin{equation} \hat{\varPi}(\boldsymbol{k}) = \frac{1}{(2{\rm \pi})^2} \int_{\mathbb{R}^2} \varPi({\boldsymbol{x}}) {\rm e}^{-\mathrm {i} \boldsymbol{k}\, \boldsymbol{\cdot}\, \boldsymbol{x} } \, \mathrm{d}\kern0.7pt \boldsymbol{x}, \quad \varPi(\boldsymbol{x}) = \int_{\mathbb{R}^2} \hat{\varPi}({\boldsymbol{k}}) {\rm e}^{\mathrm {i} \boldsymbol{k}\, \boldsymbol{\cdot}\, \boldsymbol{x} } \, \mathrm{d}\boldsymbol{k}. \end{equation}

From the definition of $\varPi$, it follows that $\hat {\varPi }({\boldsymbol {k}}) = \langle \hat {\eta }({\boldsymbol {k}})\hat {\eta }({-{\boldsymbol {k}}}) \rangle$, where $\hat {\eta }$ is the Fourier transform of the stochastic forcing $\eta$ in (1.1). In addition, because $\eta$ is real, $\hat {\eta }({-{\boldsymbol {k}}})= \hat {\eta }^*(\boldsymbol {k})$, and it follows that $\hat {\varPi }=\langle |\hat {\eta }(\boldsymbol {k})|^2\rangle$ is real and has the symmetry $\hat {\varPi }({\boldsymbol {k}})=\hat {\varPi }(-{\boldsymbol {k}})$. Further, $\eta$ must be normalised so that $\varepsilon$ is the energy injection rate in (1.1), which in Fourier space corresponds to the constraint

(2.9)\begin{equation} \int_{\mathbb{R}^2} \frac{\hat{\varPi}({\boldsymbol{k})}}{2|\boldsymbol{k}|^2}\, \mathrm{d}\boldsymbol{k}=1. \end{equation}

The SY14 analysis exploits the linearity of (2.6) to separate the contribution to $\langle u'v' \rangle$ from forcing at each wavevector $\boldsymbol {k}$. It turns out, because there is no intrinsic length scale associated with a constant shear flow, that the amplitude of the wavevector at which energy is injected does not enter the SY14 analysis. Rather, for forcing at each $\boldsymbol {k}$, the SY14 result depends only on the phase angle $\phi =\tan ^{-1}{(k_y/k_x)}$. Consequently it is possible to express the general result in terms of an energy input density $\rho ^\varepsilon (\phi )$, which is defined to be

(2.10)\begin{equation} \rho^\varepsilon(\phi)= \int_0^\infty \frac{ \hat{\varPi}(k \boldsymbol{e}_k(\phi) ) }{k^2} \, k \,{\rm d}k, \quad \mbox{where}\ \boldsymbol{e}_k(\phi)=(\cos{\phi},\sin{\phi})^{\rm T}, \ -\frac{\rm \pi}{2}<\phi \leqslant \frac{\rm \pi}{2}. \end{equation}

The symmetry $\hat {\varPi }(\boldsymbol {k})=\hat {\varPi }(-\boldsymbol {k})$ means that $\rho ^\varepsilon (\phi )= \rho ^\varepsilon (\phi +{\rm \pi} )$, which explains the restriction to $-{{\rm \pi} }/{2}<\phi \leqslant {{\rm \pi} }/{2}$. The energy input constraint means that

(2.11)\begin{equation} \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \rho^\varepsilon(\phi) \, \mathrm{d}\phi=1, \end{equation}

which means that $\rho ^\varepsilon (\phi )$ can be thought of as a density function describing the distribution of the energy input with respect to the phase angle (cf. a probability density). Note that the case of isotropic forcing, i.e. forcing on an annulus in wavenumber space in the idealised limit of an infinite domain, corresponds to $\rho ^\varepsilon (\phi )=1/ {\rm \pi}$ (const.).

The key result of this section, which is obtained in Appendix A by adapting the general approach of SY14, is that the momentum flux can be expressed as

(2.12)\begin{equation} \langle u'v' \rangle= \frac{\varepsilon}{\gamma} \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \rho^\varepsilon (\phi) K( \phi, \alpha) \, \mathrm{d}\phi, \end{equation}

where $\alpha =2\mu /\gamma$ is a parameter governed by the ratio of the Rayleigh friction rate $\mu$ to the shear rate $\gamma$. The kernel function $K(\phi,\alpha )$ appearing in (2.12) quantifies the contribution of a wave forced at phase angle $\phi$ to the momentum flux. It is given by

(2.13) \begin{align} K(\phi,\alpha)= 1 + \frac{1}{\alpha} |z|^2 \,{\rm Im} \left\{ e^{z}{E}_1(z) \right\}, \enspace \mbox{where}\ z(\phi,\alpha)=\alpha(-\tan{\phi}+{\rm i}), \left(-\frac{\rm \pi}{2}<\phi \leqslant \frac{\rm \pi}{2}\right). \end{align}

Here ${E}_1(z)$ is the complex exponential integral with the branch cut on the negative real axis

(2.14a,b)\begin{equation} {E}_1(z) := \int_z^{\infty} \frac{e^{{-}t}}{t}\mathrm{d}t,\quad |\mathrm{Arg}(z)|<{\rm \pi}. \end{equation}

Note that for fixed $\alpha$ the branch cut of ${\rm E}_1(z)$ is not crossed as $\phi$ varies, and that the restriction $\phi <|{\rm \pi} /2|$ ensures that $\tan {\phi }$ is continuous. Note also that $K(\phi,\alpha )$ is an even function of $\alpha$.

The result (2.12) is an alternative, more accessible presentation of the results in equations (B6)–(B9) of SY14. Our expression has several advantages over those of SY14. First $\mathrm {E}_1(z)$ is a tabulated function which is implemented in most mathematical software packages (e.g. as ExpIntegralE[z,1] in Mathematica and expint(z) in Matlab), meaning that (2.12) is relatively straightforward to evaluate. In addition, $\mathrm {E}_1(z)$ has well-known series expansions for both $|z|\ll 1$ and $|z|\gg 1$, which can be exploited to understand more about the behaviour of $\langle u'v' \rangle$ for any forcing $\varPi$, allowing results obtained by SY14 and Bakas & Ioannou (Reference Bakas and Ioannou2014) for rather specific forcings to be generalised to all $\rho ^\varepsilon (\phi )$.

Finally, the form (2.12) is particularly amenable to the case of forcing applied at discrete wavenumbers, e.g., on the doubly periodic domain typical for numerical simulations. In that case

(2.15)\begin{equation} \hat{\varPi}(\boldsymbol{k})=\sum_{\boldsymbol{k}_j} \rho^\varepsilon_j |\boldsymbol{k}_j|^2 \delta ( \boldsymbol{k}-\boldsymbol{k}_j), \end{equation}

where the sum is over all of the discretised wavenumber vectors $\boldsymbol {k}_j$ generated by the domain, and $\rho ^\varepsilon _j$ is the fraction of energy injected at wavenumber $\boldsymbol {k}_j$. By construction $\sum _{\boldsymbol {k}_j} \rho ^\varepsilon _j=2$. The result (2.12) then simplifies to

(2.16)\begin{equation} \langle u'v' \rangle= \frac{\varepsilon}{\gamma} \sum_{\boldsymbol{k}_j, k_x>0} \rho^\varepsilon_j K( \phi_j, \alpha), \end{equation}

which is a convenient form for accurate comparison with the numerical results in the following. Forcing on the zonal mean (i.e. with $k_x=0$, or at $\phi = \pm {\rm \pi}/2$) self-evidently does not contribute to the wave-driven momentum flux $\langle u'v' \rangle$, which explains why the sum is over wavevectors with $k_x>0$.

Given a specific energy input density $\rho ^\varepsilon (\phi )$, a useful alternative form of (2.12) is

(2.17)\begin{equation} \langle u'v' \rangle = \frac{\varepsilon}{2\mu} \, G\left( \frac{2 \mu}{\gamma} \right) \quad \mbox{where} \ G(\alpha) := \alpha \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \rho^\varepsilon (\phi) K( \phi, \alpha) \, \mathrm{d}\phi. \end{equation}

One reason (2.17) is useful is that it is often helpful to scale the momentum flux $\langle u'v' \rangle$ in terms of the expected equilibrium energy $\mathcal {E}=\varepsilon /2\mu$ in the system. The function $G(\alpha )$ is a non-dimensional odd function of $\alpha$ which is unique to each energy input density $\rho ^\varepsilon (\phi )$.

2.4. General properties of the momentum flux formula

In order to understand the general momentum flux formulae (2.12) and (2.17), and to extend and place in context previous results by SY14, WB19 and others, it is useful to understand the behaviour of the kernel function $K(\phi,\alpha )$. Figure 1 shows $K(\phi,\alpha )$ for several values of $\alpha$, with panel (a) showing the full range $(-{{\rm \pi} }/{2}<\phi \leqslant {{\rm \pi} }/{2})$ and panel (b) showing a scaled region (around $\phi = {{\rm \pi} }/{2})$ to illustrate some of the self-similar behaviours as $\alpha \to 0$. Some key properties of $K(\phi,\alpha )$, established in Appendix B and presented here for the case $\alpha >0$, are:

1. (i) $K(\phi,\alpha ) < 1$ for all $\phi \in (-{{\rm \pi} }/{2},{{\rm \pi} }/{2})$ and all $\alpha \in \mathbb {R}\setminus \{0\}$;

2. (ii) the limiting values at the boundaries of the interval are $K(\pm {{\rm \pi} }/{2},\alpha )=0$;

3. (iii) $K(\phi,\alpha )$ has a single minimum $K_-(\alpha )$; in the limit $\alpha \to 0$, $K_-(\alpha ) \sim -4{\rm \pi} e^{-2}/\alpha$, and the location of the minimum is asymptotically close to $\phi ={\rm \pi} /2-\alpha /2$;

4. (iv) in the limit $\alpha \to 0$,

   (2.18)\begin{equation} K(\phi,\alpha) = 1-\alpha \left(\phi+\frac{\rm \pi}{2}\right) \sec^2\phi+O(\alpha^2 \log{\alpha}); \end{equation}
   note that this formula applies to fixed $\phi$, and therefore cannot be used to approximate $K(\phi,\alpha )$ when $|\phi \pm {\rm \pi}/2| \sim O(\alpha )$, meaning that it cannot describe the minima of $K(\phi,\alpha )$ shown in figure 1;

5. (v) as $\alpha \to \infty$,

   (2.19)\begin{equation} K(\phi,\alpha) ={-}\alpha^{{-}1} \sin{2\phi}+2\alpha^{{-}2}\cos{\phi} \cos{3\phi}+O(\alpha^{{-}3}); \end{equation}

6. (vi) we have

   (2.20)\begin{equation} \int_{-{\rm \pi}/2}^{{\rm \pi}/2} K(\phi,\alpha) \, {\rm d}\phi = 0. \end{equation}

Figure 1. (a) The function $K(\phi,\alpha )$ against $\phi \in [-{\rm \pi} /2,{\rm \pi} /2]$ for various values of $\alpha$. When $\alpha \gg 1$, $K$ is sinusoidal (blue). As $\alpha \to 0$, $K\to 1$ for all $\phi$ except an $O(\alpha )$ region at $\phi ={\rm \pi} /2$, where an $O(\alpha ^{-1})$ minimum is obtained. (b) The same plot as the (a) with rescaled axes illustrates the self-similar behaviour of $K(\phi,\alpha )$ about its minimum. Dotted black plots the $\alpha \to 0$ solution $\alpha K(\phi,\alpha )=-{\rm \pi} (\alpha /\theta )^2\exp (-\alpha /\theta )$ found in Appendix B.

Much of the structure of $K(\phi,\alpha )$ can be understood by considering the momentum fluxes generated by plane waves with different $\phi$, and how those plane waves evolve in a linear shear flow, under the so-called Orr mechanism (Orr Reference Orr1907). When friction is large ($\alpha \to \infty$), plane waves are dissipated before they evolve under the action of the shear flow, and $K(\phi,\alpha )$ is an odd function of $\phi$ determined by the momentum flux of waves at the local phase angle $\phi$ at which they are generated. Hence, $K(\phi,\alpha )$ is negative for positive $\phi$ and positive for negative $\phi$ (see explanation on pp. 516–517 of Vallis Reference Vallis2006). When friction is low ($\alpha \to 0$), by contrast, the waves generated at angle $\phi$ will be long-lived and as time evolves their phase angle will decrease monotonically, as the wave is advected by the shear flow. Waves generated at an initial angle $\phi$ close to ${{\rm \pi} }/{2}$, i.e. with an extreme tilt against the shear, undergo considerable transient growth due to the Orr mechanism as their phase angle approaches zero. This transient growth occurs on a timescale $1/\gamma \theta$, where $\theta = {{\rm \pi} }/{2}-\phi$. As they evolve, these waves are also damped on a timescale $\mu ^{-1}$. The waves which generate the most extreme negative momentum fluxes are therefore those with $\theta \approx \mu /\gamma$, a value which allows for strong transient growth prior to the wave being dissipated, but with sufficient dissipation occurring that the wave is attenuated significantly by the time its phase angle becomes negative, so that the positive momentum fluxes during this stage of its life cycle are insufficient to cancel the earlier negative momentum fluxes.

Properties (i)–(vi) are useful as they make it relatively straightforward to prove and extend previous results regarding $\langle u'v' \rangle$ in a simple setting. The most important of these are as follows.

1. (1) Bounds on the momentum flux: As $K(\phi,\alpha )$ is bounded on $[-{{\rm \pi} }/{2},{{\rm \pi} }/{2}]$ exact bounds on $\langle u'v' \rangle$ follow naturally, which apply to any choice of the energy input density $\rho ^\varepsilon (\phi )$. For positive shear ($\gamma >0$) these are

   (2.21)\begin{equation} \frac{\varepsilon}{\gamma} K_-(\alpha) \leqslant \langle u'v' \rangle \leqslant \frac{\varepsilon}{\gamma} K_+(\alpha),\quad \mbox{where}\ K_+(\alpha)=\sup_{\phi \in \left[-{\rm \pi}/{2},{\rm \pi}/{2}\right]} \{ K(\phi,\alpha) \}, \end{equation}
   and $K_-$ is the corresponding infimum. Note that property (i) ensures that $K_+ < 1$. Each bound can be attained exactly by setting $\rho ^\varepsilon (\phi )=\delta (\phi -\phi _\pm (\alpha ))$ where $\phi =\phi _\pm (\alpha )$ denotes the locations of the supremum and infimum, respectively. Interestingly, the upper bound improves only very slightly on the bound $\langle u'v' \rangle \leqslant \varepsilon /\gamma (1+\alpha )$ found by SY14 from the energy power integral. The lower bound, which satisfies $K_-(\alpha ) \sim 4{\rm \pi} e^{-2}\alpha ^{-1}$ when $\alpha \to 0$, is only useful when $\alpha$ is order unity. Equation (2.13) allows $K_{\pm }$ to be evaluated numerically to high accuracy, and the results are shown in figure 2 along with the SY14 bound.

2. (2) The low friction limit (WB19 result): A key question raised in the introduction is under what circumstances is the low-friction approximation (1.3), $\langle u'v' \rangle \approx \varepsilon / \gamma$ valid? It turns out that a simple necessary condition for (1.3) to be valid in the limit $\alpha \to 0$, is that there is a wave angle cut-off in the forcing density $\rho ^\varepsilon (\phi )$. That is, it is necessary that there exists a cut-off angle $\phi _{\alpha }$ such that $\rho ^\varepsilon (\phi )=0$ for all $\phi$ satisfying $\phi _{\alpha } <|\phi |<{\rm \pi} /2$. For example, such a cut-off occurs naturally for simulations in a finite domain when the forcing is on an annulus in wavenumber space, because in this case $k_x$ is bounded below and $k_y$ is bounded above. In this case (1.3) can be refined by inserting the expansion (2.18) from property (iv) into (2.12), to give (for $\alpha >0$)

   (2.22)\begin{align} \langle u'v' \rangle & = \frac{\varepsilon}{\gamma} \left( 1 - \alpha \int_{-\phi_{\alpha}}^{\phi_{\alpha}} \left(\phi+\frac{\rm \pi}{2} \right) \sec^2{\phi} \, \rho^\varepsilon(\phi) \, \mathrm{d}\phi \right)+ O(\alpha^2 \log{\alpha})\nonumber\\ & \to \frac{\varepsilon}{\gamma} \quad \mbox{as}\ \alpha\to 0. \end{align}
   The phase angle cut-off is evidently necessary in (2.22) in order for the leading-order correction term to avoid the singularities in the integrand at $\phi =\pm {\rm \pi}/2$ and consequently remain bounded.
   

   Figure 2. Plot illustrating the range of $K(\phi,\alpha )$. Plotted are the $\phi$-supremum and $\phi$-infimum, $K_+(\alpha )$ (red) and $K_-(\alpha )$ (blue), respectively. The shaded grey region therefore is the range of $K(\phi,\alpha )$. The dashed black line corresponds to the upper bound for Reynolds stress established in SY14.

   In practice, while formally valid for $\alpha \ll 1$, (2.22) will be an accurate approximation to (2.12) only if $|{{\rm \pi} }/{2}-\phi _{\alpha } | \gg \alpha$, in order that the leading correction term in (2.22) remains $O(\alpha )$. Although this condition will always be formally satisfied if $\phi _{\alpha }$ is fixed and the limit $\alpha \to 0$ is taken, long domains which allow for forcing at wavenumbers with $|k_x| \ll |k_y|$ will have $|{{\rm \pi} }/{2}-\phi _{\alpha } | \ll 1$, and will therefore require $\alpha$ to be very small before (2.22) becomes accurate. A further practical point is that, when applied to slowly varying jets, (2.22) cannot hold everywhere because no matter how small the friction, $\alpha$ will not remain small in boundary layer regions around the jet extrema, where $|U_y|/\mu \lesssim 1$.

3. (3) The high-friction/low-shear limit: The behaviour in the limit $\alpha \to \infty$, corresponding to high friction or weak shear, is always important in flows with jets close to the jet extrema. This limit is also relevant for studying the emergence of $\beta$-plane jets from an infinitesimal shear (see, e.g., Bakas & Ioannou Reference Bakas and Ioannou2013; Bakas, Constantinou & Ioannou Reference Bakas, Constantinou and Ioannou2015). The behaviour of $\langle u'v' \rangle$ in the limit $\alpha \to \infty$ under the local theory is found by inserting the expansion (2.19) from property (v) into (2.12). A cancellation of leading terms follows and the result is

   (2.23)\begin{align} \langle u'v' \rangle={-} \frac{\varepsilon}{2\mu} \left( \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \sin{2\phi} \rho^\varepsilon(\phi) \, \mathrm{d}\phi -\frac{2}{\alpha} \int_{-{\rm \pi}/2}^{{\rm \pi}/2} \cos{\phi} \cos{3\phi} \rho^\varepsilon(\phi) \, \mathrm{d}\phi \right) +O(\alpha^{{-}2}). \end{align}
   Reassuringly, (2.23) shows explicitly that $\langle u'v' \rangle$ is bounded as $\alpha \to \infty$ for all possible $\rho ^\varepsilon (\phi )$, consistent with the findings of SY14. The leading term can be recognised as the momentum flux arising from the high-friction solution of (2.6), in which all terms involving the shear are neglected to give $\mathcal {Z}=\varepsilon \varPi / 2\mu$. If $\rho ^\varepsilon (\phi )$ is taken to be an even function of $\phi$, as is typical in simulations, the first term vanishes and the second term is dominant. In this case $\langle u'v' \rangle$ is linear in the shear $\gamma \equiv U_y$ near to the jet extrema.

4. (4) The isotropic forcing paradox: (2.20) of property (vi) ensures that if the forcing in (2.12) is isotropic, i.e. $\rho ^\varepsilon (\phi )=1/{\rm \pi}$, then

   (2.24)\begin{equation} \langle u'v' \rangle= 0, \quad\mbox{for all}\ \alpha. \end{equation}
   This remarkable result, that an isotropic forcing leads to zero momentum flux, was discovered by Farrell (Reference Farrell1987) and was strongly emphasised by SY14. At first glance, considering the limit $\alpha \to 0$, (2.24) appears to be in contradiction to the WB19 result (1.3). However, isotropic forcing with $\rho ^\varepsilon (\phi )=1/{\rm \pi}$ does not satisfy the wave angle cut-off property which is required for (1.3) to hold. The subtlety here is the remarkable structure of $K(\phi,\alpha )$ which allows properties (i) and (vi) to hold simultaneously as $\alpha \to 0$. Essentially, an increasingly thin boundary layer close to $\phi ={\rm \pi} /2$, illustrated in figure 1(b), allows the integral property (vi) ($\int _{-{\rm \pi} /2}^{{\rm \pi} /2} K(\phi,\alpha ) \, {\rm d}\phi = 0$) to hold even as $K(\phi,\alpha )\to 1$ for all $\phi \in (-{{\rm \pi} }/{2}, {{\rm \pi} }/{2})$. Physically, what is occurring as $\alpha \to 0$ is that only those waves with the most extreme tilt against the shear ($k_x \ll k_y$) can contribute a negative momentum flux, but the momentum flux for these waves becomes increasingly large in magnitude because they become very long-lived when friction is low. In any finite domain, with forcing confined to an annulus in wavenumber space, waves with such extreme tilts will not exist and (1.3) holds as $\alpha \to 0$. This reasoning provides an explanation as to how numerical experiments, which use a discretised ‘isotropic’ forcing in a finite domain (as in, e.g., Scott & Dritschel Reference Scott and Dritschel2012; Srinivasan & Young Reference Srinivasan and Young2012; Bakas & Ioannou Reference Bakas and Ioannou2014; Constantinou et al. Reference Constantinou, Farrell and Ioannou2014), can result in jet flows with decidedly non-zero $\langle u'v' \rangle$ in the context of the local theory only. An alternative perspective on the same issue is given in Frishman (Reference Frishman2017). There, the focus is on the length of time required for the negative momentum fluxes of the ‘extreme tilt’ waves to become established, which becomes increasingly long as $\alpha \to 0$. As a result, switching the order of the limits $\alpha \to 0$ and $t \to \infty$ in the solution of (2.6) results in isotropic forcing giving (1.3) instead of (2.24). Here our focus is on steady equilibria and consequently the limit $t \to \infty$ is always taken first. Another consideration is that, in any realistic scenario where the background shear is not exactly linear, next order corrections to the linear theory will also produce non-zero momentum flux. For example, Srinivasan & Young (Reference Srinivasan and Young2012) find that the correction term in the momentum flux is proportional to $\beta ^2 U'''$, and a physical explanation for this result is given in Bakas & Ioannou (Reference Bakas and Ioannou2013, Reference Bakas and Ioannou2019).

Given a specific energy input density $\rho ^\varepsilon (\phi )$, (2.17) shows that momentum flux is determined by the function $G(\alpha )$. Figure 3 shows $G(\alpha )$ for some simple examples of different wave forcings (WF1, WF2 and WF3 hereafter, see the caption), together with the corresponding limiting forms obtained from (1.3) and (2.22) (valid as $\alpha \to 0$) and (2.23) (valid as $\alpha \to \infty$). The asymptotic expressions are seen to do a good job of approximating $G(\alpha )$ for $\alpha <10^{-1}$ and $\alpha >10^1$, with (2.22) representing a significant improvement on (1.3). However, both the small $\alpha$ and large $\alpha$ expressions are inaccurate in the range $10^{-1} \lesssim \alpha \lesssim 10^1$. In figure 3 there is a notable contrast between $G(\alpha )$ for forcings for which $\rho ^\varepsilon (\phi )$ is an even function of $\phi$ (figure 3c–f) compared with otherwise (figure 3a,b). In the latter case, $G(\alpha )$ tends towards a constant value as $\alpha \to \infty$, consistent with the fact that the forcing in this case generates a non-zero momentum flux even in the absence of a shear flow (see the leading term in 2.23).

Figure 3. Examples of the SY14 momentum flux function $G(\alpha )$ defined in (2.17) for the wave forcings WF1, WF2 and WF3. Panels (b,d, f) show $G(\alpha )$ (purple), the limiting forms (2.22) (red curves) and (2.23) (orange curves) and the WB19 result (1.3, $G(\alpha )=\alpha$) (dashed blue). Panels (a,c,e) show the patterns of the wave forcings WF1, WF2 and WF3 in wavenumber space: WF1, $\phi _1={\rm \pi} /4$; WF2, $\{ \phi _1,\phi _2 \} = \{ -{\rm \pi} /4, {\rm \pi}/4\}$ and $\rho ^\varepsilon _j=1/2$; WF3, $\phi _j=\tan ^{-1}{(\,j/8)}$ for $j=-8,\ldots,0,\ldots,8$ and $\rho ^\varepsilon _j \propto \cos ^2\phi _j$.

2.5. The SY14 momentum flux closure applied to a jet flow

Up to this point the focus has been on constant shear flows with $U(y)=\gamma y$. How relevant is the SY14 expression for the momentum flux when $U(y)$ is instead a smooth jet flow? Intuitively, when the parameter $F$, which determines the ratio of the forcing scale to the jet scale, satisfies $F \ll 1$ one would expect that the eddies ‘see’ only a linear shear locally and SY14 will be accurate. A key question is how small does $F$ have to be in practice for the SY14 formula to hold.

The reason this question is important is that if a regime exists in which SY14 holds, the question of whether the equilibrium jet profile $U(y)$ is ever determined by a purely zonal closure can be addressed. The SY14 zonal closure equation is obtained by inserting (2.17) into (2.5) to give

(2.25)\begin{equation} \mu U ={-}\frac{\varepsilon}{2 \mu} \partial_y G\left(\frac{2\mu}{U_y}\right) - \nu_{2n}(-\partial_{yy})^{n}U, \end{equation}

where the function $G$ is defined in (2.17). In principal (2.25), which becomes a second-order ordinary differential equation for $U$ when hyperdiffusion is neglected, can be solved numerically to obtain an equilibrium jet profile $U(y)$, for any $G$. A special case is $G(\alpha )=\alpha$, which corresponds to the WB19 situation in which $\langle u'v' \rangle = \varepsilon /U_y$ (see (1.3)). For this case the exact solution, neglecting hyperdiffusion, can be given in terms of the inverse of the error function (Frishman Reference Frishman2017)

(2.26)\begin{equation} U(y)=2\sqrt{\mathcal{E}}\textrm{erf}^{{-}1}(y/\sqrt{A^2{\rm \pi}\mathcal{E}}) \quad \mbox{for} \ -A\sqrt{{\rm \pi}\mathcal{E}}< y< A\sqrt{{\rm \pi}\mathcal{E}}, \end{equation}

where $\mathcal {E}=\varepsilon / 2\mu$. Here the existence of a latitude $y=0$ where $U=0$ is assumed, and $A=U_y(0)>0$ is the local (undetermined) shear there. The profile (2.26), which becomes singular at $y= \pm A\sqrt {{\rm \pi} \mathcal {E}}$ has nevertheless been argued by WB19 to be a reasonable fit to the ‘between jet’ zonal wind profiles seen in simulations. Obviously, however, the singularities in (2.26) mean that it cannot be accurate everywhere. One of our main findings is, however, that replacing $G(\alpha )=\alpha$ with the full SY14 expression in (2.25) does not generally result in more accurate jet profiles. Much of what follows is dedicated to explaining why it does not.

A good starting point for investigating the SY14 closure is to perform a scattering experiment. That is, the CE2 equation (2.6) is solved numerically for a fixed steady jet-like wind profile $U(y)$, and the resulting momentum flux $\langle u'v' \rangle$ is compared with the theoretical prediction ((2.17) with $\gamma$ replaced with $U_y$). Note that the CE2 solution is exactly equivalent to that obtained by solving (2.3) with $\bar {u}=U(y)$ held fixed and obtaining $\langle u'v' \rangle$ by statistical averaging, however CE2 is clearly more efficient because statistical error is eliminated. In addition, when $U(y)$ is held fixed, (2.6) becomes a Lyapunov equation (see (A2)), which is much cheaper to solve numerically compared with the full time-dependent CE2 equations. High-resolution solutions of (2.6) can therefore be obtained, allowing for a numerical investigation of the small $F$ limit.

In a $2{\rm \pi}$-periodic domain, the scattering experiments are performed on a fixed flow $U(y)=2\sin y$, with quantisation and zonostrophy parameters $Q=1.22$ and $Z=1.94$, respectively. The value of $Q$ is consistent with a single jet in the domain, and the value of $Z$ allows the boundary layers near jet extrema to be resolved, which becomes prohibitively expensive at higher $Z$. By varying the forcing wavenumber $k_f$ across experiments, we then investigate the effect of varying $F$ in the range 0.04–0.84. Figure 4 shows the resulting momentum flux $\langle u'v' \rangle$ for the wave forcing patterns WF2 and WF3, illustrated in figure 4(a,d), and described further in the caption to figure 3. Figures 4(b,e) show the main comparison and demonstrate that, away from the jet extrema, in both cases there is good agreement between the CE2 scattering experiments at all three values of $F$ and the SY14 asymptotic result (2.17), and to a lesser extent the simpler expression (1.3).

Figure 4. Results from scattering experiments with $U(y)=2 \sin y$ and $Q=1.22$, $Z=1.94$ and three different values of $F$ (corresponding to forcing wavenumber $k_f=(8,32,128)$ respectively, or $F=(0.677,0.169,0.042)$ in (a–c), and $F=(0.844,0.211,0.053)$ in (d–f)). Panels (a–c) and (d–f) show the results for the wave forcings WF2 and WF3, respectively (see figure 3 caption). Panels (b) and (e) compare calculated and predicted $\langle u'v' \rangle$ across the full domain, and panels (c) and ( f) show a close-up of the situation near the east and west jet cores. The theoretical results ((1.3), WB19, dashed purple), ((2.17), SY14, red dotted line) and ((2.23), $\alpha \to \infty$, dashed blue) are also plotted.

Figures 4(c, f) show a blow-up of the situation close to the jet extrema, where $\langle u'v' \rangle$ varies rapidly in thin boundary layers, in both the CE2 calculations and in the theory (2.17). The theoretical expression (2.17) depends only on $U_y$ and not $\beta$, and therefore is symmetric at both east and west jets. The CE2 solutions include the effect of $\beta$ and therefore the convergence to (2.17) as $F \to 0$ is seen to be rather different for each jet. In figure 4(a–c), in which the forcing is concentrated on just a few waves, even the sign of the momentum flux is opposite to the (2.17) prediction for the larger values of $F$, showing that convergence is slow. When the forcing is distributed over a wider range of wavenumbers, as shown in figure 4(d–f), convergence is more uniform. In both sets of calculations convergence is significantly slower at the east jet, where the PV gradient $\beta - U_{yy}$ is large and positive, compared with the west jet where the magnitude of $\beta - U_{yy}$ is smaller.

In summary, the scattering experiments show clearly that, for sufficiently small but nevertheless physically reasonable (and numerically accessible) values of $F$, the SY14 expression (2.17) can do a good job of predicting the momentum fluxes everywhere, for a steady smooth and, crucially, stable zonal flow $U(y)$. Why, then, does the closure (2.25) not describe the equilibrium jets in our $\beta$-plane turbulence calculations? This question is addressed next.

3.1. Experimental design: separation of momentum flux contributions

A standard approach in $\beta$-plane simulations of jets is to use a (near-)isotropic wave forcing $\eta$, in which the forcing is applied only at wavevectors $\boldsymbol {k}$ which lie within an annulus $(k_-<|\boldsymbol {k}|< k_+)$. However, for fixed $Q$ and $Z$, the equilibrium jet structure in nonlinear simulations is known to be largely insensitive to the details of the forcing mechanism (e.g. Scott & Dritschel Reference Scott and Dritschel2012). Here, the freedom to choose the forcing structure is exploited by instead using the forcing WF3, which is concentrated on a single zonal wavenumber $k_x=k_f$ (see the previous discussion). In all but one simulation, WF3 is augmented (slightly) by also including an extremely weak forcing applied to waves with $k_x \ne k_f$ to allow for the excitation of instabilities. The reason for this is twofold.

1. (i) The first is to allow the momentum fluxes associated with waves directly forced by $\eta$ and those generated by secondary instabilities to be clearly distinguished. Exploiting the properties of the Fourier transform in $x$, it is helpful to first decompose the momentum fluxes into contributions from each zonal wavenumber $k_x$,

   (3.1)\begin{equation} \langle u'v' \rangle_{k_x}(y) = 2\,\mathrm{Re}\left( \tilde{u}_{k_x}(y) \, \tilde{v}_{k_x}^* (y) \right), \end{equation}
   where $\tilde {u}_{k_x}$ and $\tilde {v}_{k_x}$ denote the coefficient of the $k_x$ term in the $x$-Fourier transform of $u'$ and $v'$, respectively. This allows for the decomposition of the momentum flux $\langle u'v' \rangle =\langle u'v' \rangle _{D}+\langle u'v' \rangle _{S}$ into directly forced and secondary components defined to be
   (3.2)$$\begin{gather} \langle u'v' \rangle_{D} = \langle u'v' \rangle_{k_f} \end{gather}$$
   (3.3)$$\begin{gather}\langle u'v' \rangle_{S} = \sum_{k_x \ne k_f} \langle u'v' \rangle_{k_x}. \end{gather}$$
   The above decomposition pre-supposes that no secondary instabilities will occur at $k_x=k_f$, which is the case in our simulations, because $F \ll 1$ and there is a clear scale separation between the jet-scale secondary instabilities and the forcing. Forcing on wavenumbers in an annulus does not allow for a clean decomposition, because the zonal wavenumbers associated with the secondary instabilities will also be forced directly. In the following, it will also prove useful to further decompose $\langle u'v' \rangle _{S}$ to help distinguish between instabilities on the east and west jets.

2. (ii) The second motivation for our forcing choice is to introduce a long-wave cut-off for the excitation of zonal wave modes. It is now widely accepted that baroclinic instability is the primary mechanism driving the extratropical jets on the giant planets (e.g. Liu & Schneider Reference Liu and Schneider2010; Young, Read & Wang Reference Young, Read and Wang2019; Read et al. Reference Read, Kennedy, Lewis, Scolan, Tabataba-Vakili, Wang, Wright and Young2020) and, if $\eta$ is to represent this process, it should have a long-wave cut-off in the zonal direction (i.e. there should be no stochastic excitation of waves with $|k_x|< k_{min}$ for some $k_{min}>0$, cf. the Eady or Phillips model for baroclinic instability; e.g. Vallis Reference Vallis2006). As discussed previously, the existence of a long-wave cut-off is quantitatively important because very long zonal waves can make a disproportionately large (and, if not treated carefully, an unphysical) contribution to the momentum flux in the local theory (see also Frishman Reference Frishman2017).

With the described forcing, a reference simulation (REF) is performed at parameter settings $(Z,Q,F)= (5.05,2.91,1.01)$. (This value of $F$ corresponds to forcing at $k_f=16$.) For the purposes of comparison and model validation, simulation REF is repeated in the full nonlinear equations (1.1), the QL equations (2.2)–(2.3) and the CE2 equations (2.5)–(2.6). The resolution in the nonlinear simulation is $256^2$ Fourier modes, whereas for QL and CE2 it is $16 \times 256$ because wavenumbers $|k_x| > k_f$ are not necessary. Hyper-diffusivity in all cases is $\nu _4=2.5\times 10^{-8}$, which is sufficient to remove enstrophy at small scales in the nonlinear simulation, but remains insignificant in the energy balance. The results are reported next.

3.2. Momentum flux decomposition in an equilibrated jet flow

Figure 5 shows the spin-up and equilibration of jets in the QL, CE2 and nonlinear (NL) simulations. The left panels are Hovmöller plots of the zonal mean wind $U(y,t)$ and the right panels show a time-mean taken over the final $0.5\ {\rm \mu}^{-1}$ time units in each simulation. The QL and CE2 simulations result in very similar jet profiles, as expected because the latter describes the statistics of the former. The jets in the NL simulation are less symmetric, with a more rounded westward jet, reflecting the fact that at $Z\approx 5$ the NL simulation only approaches the high zonostrophy ($Z \gg 1$) regime in which QL is expected to be a good approximation to (1.1). Nevertheless figure 5 confirms the usefulness of QL and CE2 as simplified models of the full nonlinear behaviour.

Figure 5. Left panels: Hovmöller plots of zonal mean flow $U(y,t)$. (a) QL simulation with wave forcing WF3 with the ‘seed forcing’ omitted. (b) QL simulation with wave forcing WF3 including the seed forcing. (c) CE2 simulation with WF3 including the seed forcing. (d) Fully nonlinear simulation with WF3. Right panels: time average of $U(y,t)$ over the last $0.5\ {\rm \mu}^{-1}$ time period. In all simulations $Z=5.05$, $Q=2.91$ and $F=1.01$.

Figure 5(a) shows a QL simulation which gives a simple illustration of the importance of secondary instabilities for the equilibrated jet profiles. In this simulation the ‘seed forcing’ in WF3 is switched off, meaning that waves with zonal wavenumber $k_x \ne k_f$ cannot be excited. Despite the fact that, energetically speaking, the forcing is effectively unchanged, because the fraction of energy input into the seed forcing into waves with $k_x \ne k_f$ is less than $10^{-4}$ of the total, the outcome of the simulation is radically different. Thin jets are formed with widths far less than the Rhines scale seen in the other simulations. Further QL simulations (not shown) demonstrate that the equilibrated jets in the QL simulation (second row) are independent of the amplitude used for the seed forcing, providing further evidence that the role of the seed forcing is to excite instability.

To understand the mechanisms at play it is informative to look first at the momentum flux decomposition for the CE2 and QL simulations, shown in figure 6(a–h), because these are somewhat cleaner than their NL counterpart and give near-identical results. Figures 6(a,e) show the time-mean zonal wind $U(y)$ for reference, and figures 6(b–d, f–h) show the contributions to the momentum flux convergence, i.e. the wave-induced force on the zonal flow, for different zonal wavenumbers. Results for QL are obtained using a ‘jet-following’ averaging procedure to compensate for any gradual evolution in jet position which would otherwise smear out the results. The jet-following procedure resets in the origin in the $y$-direction prior to averaging, by phase-shifting all quantities in Fourier space to the phase of the jets, as determined by the second Fourier coefficient of $U(y,t)$. Using this method, averages are calculated over a long period ($50\ {\rm \mu}^{-1}$) of equilibrium, to obtain good statistical convergence. CE2 adopts an alternative averaging procedure, because equilibrium solutions for CE2 undergo small oscillations (Marston et al. Reference Marston, Qi and Tobias2019) about a steady jet configuration. In the solution presented the oscillations are relatively small, but a temporal average is taken anyway once a steady state is reached (from $t=49\ {\rm \mu}^{-1}$ to $50\ {\rm \mu}^{-1}$).

Figure 6. Long-time equilibrium mean quantities from the (a–d) CE2, (e–h) QL and (i–l) NL simulations reported in figure 5. Panels (a,e,i) show the average mean wind profile $U(y)$. The other panels plot $-\mu ^{-1}\partial \langle u'v'\rangle _{k_x}$ for significant modes $k_x$. The CE2 and QL results are similar, clearly identifying distinct instabilities at the westward (shaded red) and eastward (shaded blue) jets which counteract the jet-sharpening contribution from the forcing wave $k_x=k_f$. Analysis of the NL simulation reveals wave ranges performing similar roles as the QL and CE2 simulations: e.g. $k_x=1,2$ damps westward jet growth, $k_x=3\unicode{x2013}7$ all have a similar structure counteracting the eastward jet sharpening and the waves $k_x>13$ perform the jet sharpening role of the forcing wave $k_f$.

The directly forced contribution to the momentum flux convergence $-\partial _y \langle u'v' \rangle _{D}$ (for $k_x=k_f=16$) is shown in figure 6(d,h). The secondary instability contribution $-\partial _y \langle u'v' \rangle _{S}$ is found to be entirely dominated by zonal wavenumbers $k_x=3, 4$ and $6$ and these contributions are plotted in figure 6(b,c, f,g). A striking feature of figure 6 is that at the jet cores there is a cancellation between two large terms: the direct wave induced force which acts to accelerate both the westward and eastward jets, and the contributions from the secondary instabilities which are decelerating. At the westward jets, in particular, the direct and secondary momentum flux convergences are about an order of magnitude larger than the frictional force in the momentum balance in (2.5). The main contributions at the eastward and westward jets have somewhat different scales ($k_x=6$ and $k_x=3,4$, respectively), which we show in the following is characteristic of the scale of barotropic instability on each jet. The secondary instabilities are non-local, however their influence is not felt across the full domain, and there are wide regions in the jet flanks where $\langle u'v' \rangle _S\ll \langle u'v' \rangle _D$. Figure 7 shows that, within these wide regions, SY14 formula (2.12) is an excellent prediction for the Reynolds stress.

Figure 7. Further analysis of the CE2 simulation in figure 6. The half domain $[-{\rm \pi} /4,3{\rm \pi} /4]$ is reported to view a single eastward and westward jet. (a) Mean zonal wind profile (blue) and the PV gradient (orange). (b) Momentum flux quantities scaled by $\mathcal {E}$. These are the CE2 results for $\langle u'v' \rangle$ (dotted blue), $\langle u'v' \rangle _D$ (black) and $\langle u'v' \rangle _S$ (green) and the SY14 $\langle u'v' \rangle$ (red dotted line). Sufficiently distant from the points where $\beta -U_{yy}\approx 0$ it is observed that $\langle u'v' \rangle _D\approx 0$ and the SY14 solution agrees with the CE2.

The momentum flux decomposition for the NL simulation is shown in figure 6( j–l). There are striking similarities between NL and QL/CE2 in the patterns of the momentum flux convergences at each jet, but also significant differences in both the magnitudes and latitudinal scales, as well as the breakdown by zonal wavenumber. The differences can be accounted for by the fact that $Z$ is finite in the NL simulation, and agreement between NL and QL is expected to improve at higher $Z$. As $Z$ increases, away from jet extrema the zonal eddy-eddy interactions are suppressed (e.g. Tobias & Marston Reference Tobias and Marston2013.) Indeed, in some small region around jet extrema, it is not clear why the QL approximation should hold; however, such a region is expected to be vanishingly narrow with increasing $Z$. The simulations indicate that finite $Z$ affects the NL flow statistics in two important ways.

1. (i) First, there is significant wave–wave interaction in the nonlinear simulation, which leads to wave energy being scattered in wavenumber space. This affects the direct wave-induced force $-\partial _y \langle u'v' \rangle _{D}$ by spreading its contribution over a range of wavenumbers centred on $k_x=k_f$, and similarly the spectrum of $-\partial _y \langle u'v' \rangle _{S}$ is broadened to a larger range of $k_x$.

2. (ii) Second, the average variance of the mean wind profile about its temporal mean is roughly five times greater in the NL case compared with QL (both cases are calculated using the ‘jet following’ procedure). The increased variance is exhibited as (a) less steady relative jet positions, which smears the statistics to give weaker and broader patterns for the momentum flux convergence, and (b) increased variability in the magnitude of $U$ at each jet, meaning that there are periods when the Raleigh–Kuo condition for stability is broken for a broader range of $k_x$, which results in a broader zonal spectrum of emergent waves.

These effects account for the difference between the CE2/QL and NL momentum flux decompositions in figure 6, because the latter can be viewed as a ‘smeared out’ version of the former, but is otherwise qualitatively similar. From figure 6 it is evident that the QL and CE2 statistics are almost identical. An exception is the QL $k_x = 3, 4$ contribution near the eastward jet which is absent in CE2: a small reminder that CE2 is not an exact representation of the QL statistics as its averaging procedure cannot fully model the mean flow fluctuations which occur in realisations of the QL flow. Curiously, the NL $k_x = 1, 2$ modes provide the opposite (sharpening) Reynolds stress at the eastward jets compared with $k_x = 3\unicode{x2013}7$ which is larger. Hence, the overall contribution from $\langle u'v' \rangle _S$, as in CE2/QL, is to counteract the sharpening from $\langle u'v' \rangle _D$.

In summary, momentum fluxes from the directly forced waves act to accelerate both the eastward and westward jets, and opposing these are momentum fluxes from relatively long ($k_x=1$–$7$) jet-scale waves which appear to derive from secondary instabilities. To investigate the origin of these long waves more thoroughly, a linear instability analysis of the time-mean flow is conducted next.

3.3. Linear stability analysis of the equilibrium jet flow

To investigate the origin of the emergent long waves that drive the secondary momentum fluxes $\langle u'v' \rangle _S$, a linear stability analysis of the time-mean flow in the CE2 simulation is presented next. The CE2 simulation is chosen for analysis because the time-variability of $U(y,t)$ is significantly lower than for the stochastic simulations, meaning that its time average is an excellent statistical representation of the actual state at any fixed time. One of the pitfalls of a single numerical linear stability calculation of a time-mean state of a system near marginal stability is that the results can return large number of wave modes with near-zero growth rates, making the results hard to interpret. The goal of our analysis, therefore, is to investigate whether there are states nearby to the time-mean state (in a sense to be described) in which there are specific wave modes with significant growth rates. The idea is that, as $U(y,t)$ evolves it will spend a proportion of its time in these more unstable states, exciting the strongly unstable wave modes, which will persist because they are relatively weakly damped when the system is outside the unstable regime.

The linear stability problem is formulated by seeking a solution to (2.3) of the form

(3.4)\begin{equation} \psi'(x,y,t)=\mathrm{Re}\left\{ \varPsi(y) \exp\left( \mathrm{i}k_x(x-ct) \right) \right\}, \end{equation}

resulting in the generalised eigenvalue problem

(3.5)\begin{equation} \mathcal{L} \varPsi = c \mathcal{M} \varPsi, \end{equation}

for linear operators

(3.6) \begin{equation} \left.\begin{gathered} \mathcal{L} =(\mathrm{i}k_xU+\mu-({-}1)^{n+1}\nu_{2n}(\partial_{yy}-k_x^2)^{n})(\partial_{yy}-k_x^2) + \mathrm{i}k_x(\tilde{\beta}-U_{yy}),\\ \mathcal{M} =\mathrm{i}k_x(\partial_{yy} - k_x^2). \end{gathered}\right\} \end{equation}

Here $\tilde {\beta }=\beta +\delta \beta$, where $\delta \beta$ is a perturbation to the value of $\beta$ used in the CE2 simulations, which has been introduced as a device to investigate the stability of ‘nearby’ states to the CE2 time-mean flow. Our assumption here is that the stability properties of nearby states generated by varying $\delta \beta$ are representative of those of the nearby states generated by fluctuations in $U(y,t)$. We believe this is reasonable as growth rates are largely determined by the width and magnitude of the PV gradient reversal region in which $\beta -U_{yy}$ changes sign, as suggested by the Rayleigh–Kuo necessary criterion for instability, which states that a sign change must be present for instability in the inviscid system. Note that the reference time mean profile, with $\delta \beta =0$, just satisfies the Rayleigh–Kuo criterion, with small regions of opposite sign PV gradient located at the westward jet extrema and in flanks of the eastward jet.

The generalised eigenvalue problem (3.5) is discretised on a grid of $1024$ points by replacing $\mathcal {L}$ and $\mathcal {M}$ with the matrices obtained when standard centred-difference approximations replace $y$-derivatives. The calculated growth rate $k_x c_i$ of the fastest growing mode, obtained by taking the imaginary part of the computed eigenvalues $c$ of (3.5), is plotted in figure 8 as a function of $(k_x, \delta \beta )$. As $\delta \beta$ is reduced and the system becomes more unstable, growing waves are seen to emerge at around $k_x\approx 3.7$ and $k_x\approx 6.1$, which correspond to instabilities at the westward and eastward jets, respectively. Notably, there is a significant asymmetry between strong positive growth rates at negative $\delta \beta$ and weak decay rates at positive $\delta \beta$, supporting the idea that the system need spend only a relatively small fraction of time in the unstable regime to support the emergence of these waves.

Figure 8. Linear stability analysis for the CE2 profile in figure 6. The contour plot shows the maximum eigenmode growth rate $kc_i$ as $\delta \beta$ and $k_x$ are varied. Here, $\delta \beta =\tilde {\beta }-\beta$ measures the deviation of the Coriolis parameter $\tilde {\beta }$ used for the stability analysis compared to the actual $\beta$ value used in the CE2 simulation. Two distinct unstable regions are visible at $k_x = 3, 4$ and $k_x = 6$, corresponding to instabilities at the westward and eastward jets as labelled. The unstable wavenumbers agree with those identified in figure 6. The right-hand side line plots show the normalised momentum flux divergence (solid blue) of the most unstable mode at the points indicated by the green and blue dots. For reference the PV gradient $\tilde {\beta }-U_{yy}$ is also shown (dotted red).

The latitudinal structure of the momentum flux convergence associated with each unstable waves can be calculated from the corresponding eigenvector $\varPsi$. Figures 8(b,c) show the calculations for a typical unstable wave on the westward jet (green dot) and eastward jet (blue dot), respectively. The pattern of the latitudinal structure in each case is seen to be close to those calculated in the equilibrium QL/CE2 and NL simulations, shown in figure 7. This correspondence in the momentum flux structures at each jet, together with the close matches in the respective emergent zonal wavenumbers ($k_x=3,4$ for the westward jet and $k_x=6$ for the eastward jet) provides conclusive evidence that each jet is independently marginally stable to barotropic instability, and that these barotropic instabilities are the source of the secondary instabilities which contribute to $\langle u'v' \rangle _S$ in the equilibrium simulations reported previously.

4.1. A stable monotonic zonal flow

To generate a stable monotonic zonal flow, the Rayleigh friction term in the original governing equation (1.1) is modified so that the flow is relaxed towards a prescribed profile $U_0(y)=\textrm {erf}(y)$, modelling, e.g., a wind stress term in an oceanic flow. The boundary conditions are periodic in the zonal $x$-direction, as before, but for this problem periodicity in $y$ is replaced with solid wall boundaries at $y=\pm {\rm \pi}$ at which the usual (Phillips) boundary conditions are applied. To minimise wall effects, the stochastic excitation is confined to the shear zone, so the governing equation becomes

(4.2)\begin{equation} \zeta_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \zeta + \beta v ={-} \mu (\zeta-\zeta_0) -\nu_{4}\nabla^{4}\zeta + \sqrt{\varepsilon f(y)}\eta, \end{equation}

with $\zeta _0 = -2\exp (-y^2)/\sqrt {{\rm \pi} }$. The local energy injection rate is now $\varepsilon f(y)$ with $f(y)=\exp ( - y^2/ 2 \sigma _n^{2} ) / \sqrt {2 {\rm \pi}} \sigma _n$ and $\sigma _n={\rm \pi} /8$. As the flow is monotonic it is not necessary to use the smoothed closure (SY14$\sigma$ described previously), therefore results from (4.2) can be compared with those from the unsmoothed SY14 closure equation which is (cf. (2.25))

(4.3)\begin{equation} \partial_t U ={-} \mu (U-U_0) - \frac{\varepsilon}{2\mu}\partial_y \left(f(y) G\left(\frac{2\mu}{U'(y)}\right)\right) + \nu_{4}\partial_y^{4}U, \end{equation}

with $G(\alpha )$ given by (2.17).

Figure 9 shows a comparison between the equilibrated states of the CE2 equations obtained from (4.2) with the equilibrated solution of the SY14 closure equation (4.3). The parameters $(\mu, (2{\rm \pi} )^2\varepsilon /\mu, \beta )=(0.002, 4, 2)$ have here been chosen in order that the stochastic forcing is rather weak, in the sense that the equilibrated flow remains close to the relaxation profile $U_0(y)$. It is clear (see figure 9(b)) that SY14 does an excellent job of describing the deviation from $U_0(y)$ induced by the stochastic forcing in the CE2 model (the small differences here can be attributed to finite $F$ effects). Evidently, a class of channel flows exist in which zonal local closure theories are entirely successful, replicating previous success with vortex condensate flows (Laurie et al. Reference Laurie, Boffetta, Falkovich, Kolokolov and Lebedev2014).

Figure 9. (a) Plots of $U(y)$ from CE2 (blue) and LCT (purple) for a channel flow linearly relaxed to $U_0=\textrm {erf}(y)$ (red). (b) Similar to (a), but plots the difference $U-U_0$. (c) The forcing profile $f(y)$. The forcing profile is also indicated on the two left-most panels with a blue gradient indicating the forcing magnitude.

Further experiments (not shown) reveal that, even in this simple setting, this regime in which the SY14 local zonal closure theory remains accurate is rather restricted. Increasing the forcing (higher $\varepsilon$) or broadening the forcing region (higher $\sigma _n$) tends to lead to the formation of local extrema with $U_y=0$, which lead to secondary instabilities by a mechanism to be described in the next subsection. For the reasons described previously, SY14 is sensitive to finite $F$ effects at local extrema, and consequently local extrema deserve special attention, in the following flow scenario.

4.2. A stable zonal flow with alternating jets

To generate a stable alternating jet flow, (1.1) is modified to include a ‘radiative relaxation’ term, as has been used by, e.g., Scott & Polvani (Reference Scott and Polvani2008) to model the effect of large scale radiative damping on the jets of the giant planets. The governing equation, solved on the doubly periodic domain, becomes

(4.4)\begin{equation} \zeta_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \zeta + \beta v = r\left( \bar{\psi} - \varPsi_{rad} \right) - \mu \zeta - \nu_{4}\nabla^{4}\zeta + \sqrt{\varepsilon}\eta. \end{equation}

Here the streamfunction $\psi =\nabla ^{-2} \zeta$, and $\varPsi _{rad}=\cos (y)$ so that the zonal flow is effectively relaxed towards a sinusoidal ‘deep jet’ profile (cf. Thomson & McIntyre Reference Thomson and McIntyre2016), given here by

(4.5)\begin{equation} U_0(y)=\frac{1}{1+(\mu/r)}\sin(y). \end{equation}

To maintain the stability of the equilibrium flow, it was found to be necessary to use a relatively strong relaxation and damping. This set-up thus allows an exploration of strongly relaxed jets without excessive damping of eddy activity at the forcing scale.

The parameter settings for our main calculation, using the CE2 truncation of (4.4), solved numerically on a grid of 256 points, are $(\beta,\mu,r,\varepsilon )=(22, 0.05, 0.05, 0.05/(2{\rm \pi} )^2)$. The large value of $\beta$ was found to be necessary to suppress secondary instabilities and, in fact, because $\beta \gg |U_{yy}|$ throughout the domain, has the effect of making the behaviour at the east and west jets almost symmetric. Useful insight into the expected behaviour at the jet cusps is provided by the scattering experiments of § 3, where typical momentum fluxes for a fixed sinusoidal jet were plotted for different values of the forcing scale parameter $F$. Strong sensitivity to $F$ is therefore expected, and because it is not computationally feasible to explore the $F \to 0$ limit in these interactive CE2 calculations, we take $k_f=16$ and aim to compare our results with the SY14$\sigma$ closure (i.e. (4.3) modified for this system, with the smoother (4.1) applied), which can capture the qualitative effects of finite $F$.

Figure 10 shows the comparison between the CE2 calculation and solutions obtained from the SY14$\sigma$ closure equation (solved numerically on a 2048 grid). The SY14$\sigma$ results show sensitivity to the smoothing parameter $\sigma$, reflecting the sensitivity to $F$ in the CE2 momentum flux patterns seen in figure 4. For the present calculation, an optimal value of $\sigma =\sigma _*=0.237/k_f$ results in good agreement between CE2 and SY14$\sigma$ throughout the domain. Admittedly, this good agreement requires a flow-specific parameter fit for $\sigma$, and more research is required to determine how best finite $F$ effects can be captured in a modified or extended SY14.

Figure 10. CE2 and LCT quantities from the radiatively damped experiment. Panel (a) compares the mean wind profile with $U_0$ for CE2 and SY14$\sigma$ (with the optimal $\sigma ^*=0.237/k_f$). Panel (b) investigates results in more detail by plotting the profile deviation from the radiatively relaxed profile, $U-U_0$. SY14$\sigma$ results are given for $\sigma = [0.125/k_f, \sigma ^*, 0.5/k_f]$. In all cases $k_f=16$.

The eddy-induced changes to $U$ seen in figure 10 serve to illustrate why the regime in which $\langle u'v' \rangle _S=0$, and the purely zonal closure based on SY14$\sigma$ is applicable, is restricted to a narrow region of parameter space in which the stochastic forcing is weak and the relaxation of the jet is strong. Because $\langle u'v' \rangle _D$ peaks strongly near the jet cores, even rather modest eddy forcing, which has little effect on the flow elsewhere, has a strong effect at the jet tips, acting to sharpen them. The jet curvature $U_{yy}$ at the westward jet tip and eastward jet flanks increases rapidly with forcing strength, leading to $\beta -U_{yy}$ changing sign and, thus, the potential for secondary instability. The effect is more pronounced as $F$ is reduced, or $\sigma$ is reduced in SY14$\sigma$, as seen in figure 10. Through this sharpening mechanism, stochastic forcing at a jet cusp has a much stronger tendency to lead to secondary instability than elsewhere, explaining why the stable alternating jet regime discovered here occupies such a narrow region of parameter space, especially so at low $F$. Much more important are the canonical equilibrated jets, which have active secondary instabilities, to be discussed next.

4.3. A typical equilibrated jet flow in $\beta$-plane turbulence

Finally, we are ready to address the question of whether the SY14$\sigma$ closure can be adapted to parameterise the stochastic term in the canonical NL $\beta$-plane turbulence simulation described in § 3.2 (see figure 5). Recall that SY14$\sigma$ can capture only the directly forced momentum flux $\langle u'v' \rangle _D$, and not the secondary flux $\langle u'v' \rangle _S$ due to barotropic instability, which we know from § 3 is equally important in the equilibrated momentum balance. In the absence of a (known) means to parameterise $\langle u'v' \rangle _S$, our approach is simply to allow unstable waves to develop spontaneously in the parameterised flow, by solving the deterministic nonlinear equation

(4.6)\begin{equation} \zeta_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \zeta + \beta v ={-}\mu \zeta + ({-}1)^{n+1} \nu_{2n}\nabla^{2n}\zeta - \frac{\varepsilon}{2 \mu} \left[ \partial_{yy} G\left(\frac{2\mu}{\bar{\zeta}}\right)\right]_\sigma, \end{equation}

where the square brackets $[ {\cdot } ]_\sigma$ denote that the smoother (4.1) is applied. Equation (4.6) is simply (1.1) with the stochastic term replaced by the deterministic SY14$\sigma$ parameterisation.

Figure 11 compares PV and vorticity snapshots in the NL simulation of § 3.2 with those from a corresponding integration of the parameterised equation (4.6). The parameter settings and numerical configurations are identical, and both simulations are integrated for the same length of time until an equilibrated state is reached. The value $\sigma =0.237/k_f$, found in the alternating jet experiment above, is used for the smoothing parameter in SY14$\sigma$.

Figure 11. (a,c) Snapshots of PV $\zeta +\beta y$ and relative vorticity $\zeta$ from the NL simulation reported in figures 5 and 6. (b,d) Snapshots of the same quantities from an SY14$\sigma$ simulation with the same parameters. In the SY14$\sigma$ model the stochastic forcing is replaced with the deterministic forcing term described in the main text. The zonal jet structure and scaling are comparable in each case, as are the non-zonal instabilities which emerge in the flow.

There are several striking similarities between the two simulations apparent in figure 11.

1. (i) Excellent correspondence between the structure of the equilibrated jets in the two simulations.

2. (ii) Relatively short waves $(k_x=4\unicode{x2013} 6$) propagating on the PV barriers at the core of the eastward jets. These correspond to the eastward jet instability identified in the linear analysis of figure 8.

3. (iii) Longer waves $(k_x=1\unicode{x2013} 3$) propagating on the PV barriers located at the flanks of the westward jets. Note that, although these waves are only clearly visible in snapshots of the parameterised run, their presence in the NL simulation is evident from the momentum flux decomposition of figure 6. They correspond to the waves generated by the westward jet instability identified in figure 8.

From the compiled evidence, it is clear that SY14$\sigma$ is successful because it captures the direct momentum flux $\langle u' v' \rangle _D$ sufficiently accurately to permit the flow to reach a marginally stable equilibrium in which emergent waves generate the correct secondary flux $\langle u' v' \rangle _S$. Further analysis of the momentum budget (not shown) supports this interpretation. In these simulations, the main drawback of SY14$\sigma$ is that it does not give a computational saving compared with the NL simulation, in fact the computational cost is almost identical as both runs are on the same grid over the same time period. However, this need not be true in general, particularly as $F$ is reduced. This is because, as inspection of figure 11 confirms, the length scale of waves in the parameterised simulation is set by the jet spacing, not by the forcing scale. Therefore, a scenario with $F \ll 1$ in which the parameterised run is many orders of magnitude cheaper than NL would not be difficult to set up.