Parameterisation of small-scale random forcing in β-plane turbulence

Abstract

Zonal jet formation in $\beta$-plane turbulence is investigated with the focus on whether an accurate closure can be developed for the eddy momentum fluxes due to small-scale random forcing. The approach of Srinivasan and Young (J. Atmos. Sci., vol. 71, 2014, pp. 2169–2185) is developed to give a relatively simple expression for the local Reynolds stress, due to a white-in-time random forcing with characteristic length scale much less than the jet spacing. In typical jet flows, however, it is demonstrated that the Srinivasan–Young flux is not the full story because momentum fluxes due to jet-scale waves, present as a result of distinct barotropic instabilities of the eastward and westward jets, respectively, also play a key role in the momentum balance. Numerical simulations that explicitly include the random forcing are then compared with those in which the Srinivasan–Young closure is applied. For typical jet flows, good agreement of the equilibrium zonal flow is found provided that the closure simulation is not truncated to be purely zonal, i.e. jet-scale secondary barotropic instabilities are allowed to develop. Flows in which the geometry or external forcing acts to suppress the development of secondary instabilities are also simulated, and for these flows the Srinivasan–Young closure is shown to be successful as a purely zonal closure. It is argued that vortex condensates in isotropic forced-dissipative 2D turbulence are an example of this latter situation.

1. Introduction

From the atmospheres of giant planets to terrestrial oceans, self-organisation of turbulence into jets is a ubiquitous feature of geophysical fluid dynamics (see, e.g., Galperin & Read (2019) and references therein). The two-dimensional barotropic vorticity equation on a $\beta$-plane, subject to stochastic forcing, has long been considered a key ‘toy’ model for understanding the physics of such jet formation, in particular the long-lived zonal jets of Jupiter and Saturn. The governing equation

(1.1)\begin{equation} \zeta_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \zeta + \beta v ={-} \mu \zeta + ({-}1)^{n+1} \nu_{2n}\nabla^{2n}\zeta + \sqrt{\varepsilon}\eta \end{equation}

describes the evolution of the absolute vorticity field $\zeta +\beta y$, where $\zeta$ is relative vorticity and $\beta y$ is the variable part of the planetary vorticity, under the action of advection by the velocity field $\boldsymbol {u}\equiv (u,v)^{\rm T}=\boldsymbol {\nabla }^\perp (\nabla ^{-2} \zeta )$, a linear drag with rate $\mu$, a regularising hyperdiffusion of order $n$ and coefficient $\nu _{2n}$, and a stochastic forcing $\eta$ which is white-in-time and spatially correlated over a characteristic length scale $L_f$. The forcing $\eta$ is here normalised so that $\varepsilon$ is the rate of its energy input into the system. Numerous studies based on (1.1) and its near relatives have uncovered much about its behaviour. In a doubly periodic domain of dimension $2{\rm \pi} L_d \times 2{\rm \pi} L_d$ the most physically significant non-dimensional parameters have been found to be

(1.2a–c)\begin{equation} Z=\frac{L_{Rh}}{L_\varepsilon},\quad Q= \frac{L_d}{L_{Rh}},\quad F =\frac{L_f}{L_{Rh}}. \end{equation}

Each parameter is a ratio of the ‘jet-spacing’ or Rhines scale $L_{Rh}=(U/\beta )^{1/2}$ (Rhines 1975; Williams 1978), where $U=(\varepsilon /\mu )^{1/2}$ is the turbulent velocity scale, to other length scales in the system, namely the transition length scale $L_\varepsilon = (\varepsilon / \beta ^3)^{1/5}$ of Maltrud & Vallis (1991), the domain scale $L_d$ and the forcing scale $L_f$ respectively.

Phenomenologically, in the equilibrated state, the zonostrophy parameter $Z=\beta ^{1/10}\varepsilon ^{1/20}\mu ^{-1/4}$ (Sukoriansky, Dikovskaya & Galperin 2007; Scott & Dritschel 2012) controls the ratio of energy of the zonal modes (i.e. the jet structures) to the non-zonal modes (i.e. waves and vortices). At low $Z (\ll 1)$ the zonal energy is negligible and the flow is dominated by vortices as in 2D Navier–Stokes turbulence, but as $Z$ increases (${\gtrsim }1$) jet structures form, and almost all energy is contained in the jets at larger $Z$ (${\gtrsim }5$). The quantisation parameter $Q$ simply determines the number of jets which appear in the domain, provided that $Z$ is sufficiently large for jets to form and $F$ has a more subtle effect on the jet dynamics (see, e.g., Scott & Dritschel (2012), where a range of different forcings are considered) with an apparent stronger tendency to ‘potential vorticity (PV) staircase’ (Dritschel & McIntyre 2008) formation at larger $F$. Here, our focus is on the regime of most relevance to the giant planets, $Z\gg 1$, $Q \gg 1$ and $F \ll 1$.

In the high zonostrophy limit $Z \gg 1$, it is well-established (e.g. Bouchet, Nardini & Tangarife 2013; Tobias & Marston 2013) that nonlinear exchanges of energy in (1.1) are dominated by triad interactions between pairs of wave modes and a mean flow mode (defined here as Fourier modes with zonal wavenumber $k_x \ne 0$ and $k_x=0$, respectively). As the high zonostrophy limit is approached (typically once $Z\gtrsim 5$), triad interactions between triplets of wave modes can be neglected, and the dynamics of (1.1) is well captured by a quasi-linear (QL) approximation (see (2.3)). An advantage of the QL system is that equations describing its statistics, known variously as the CE2 equations (Marston, Conover & Schneider 2008; Srinivasan & Young 2012; Ait-Chaalal et al. 2016) or the SSST equations (Farrell & Ioannou 2003, 2007), can be derived and analysed. An important result is the exact solution for the second-order moments in the CE2 system, in the special case where the mean flow is a linear shear flow (Srinivasan & Young 2014, SY14 hereafter). In the small $F$ limit, in which the excited waves are short compared to the scale of the jets and will therefore ‘see’ only a linear shear flow locally, solutions following the SY14 approach appear to hold the promise of a local closure for the ensemble average wave momentum flux $\langle u' v' \rangle$. Such a closure, a long-time goal of researchers, would allow for a single equation describing the evolution of the ensemble mean zonal wind $U(y,t)$.

Recently, using a different but essentially equivalent approach to SY14, Woillez & Bouchet (2019, WB19 hereafter) have shown that in the joint limits $Z\gg 1$, $F \ll 1$, under some mild restrictions on the forcing $\eta$ to be discussed in the following, the momentum flux closure becomes independent of the details of $\eta$ and is given by

(1.3)\begin{equation} \langle u' v' \rangle= \varepsilon / U_y. \end{equation}

A similar closure (modified slightly for polar geometry) has been postulated for the related problem of vortex condensate formation in forced 2D Navier–Stokes turbulence (Laurie et al. 2014), and can be used to predict the radial structure of the vortex condensate, thus demonstrating that a local closure theory can be successful in the right setting. However, in the $\beta$-plane jets problem, as WB19 recognised, it is clear that the closure (1.3) is inadequate. First, (1.3) fails because the equilibrated flow must include jet minima and maxima where $U_y=0$, and where the predicted momentum fluxes of (1.3) are therefore singular. Second, (1.3) does not depend on $\beta$, and therefore it cannot possibly lead to equilibrated jets with spacing on the Rhines scale $L_{Rh}$. These observations, however, do not necessarily preclude a local theory. An extended local theory can be hypothesised, which could allow for higher-order terms in $F$ to introduce dependency on $\beta$ and the flow curvature $U_{yy}$, and which would require a particular focus on boundary layer regions where $|U_y|=O(F)$. In principle, such a theory might resolve both of the previously described issues and lead to a closure that is valid everywhere in the flow.

The present work aims, using a series of carefully designed CE2 and fully nonlinear numerical simulations, to demonstrate the following.

1. (i) That in the right parameter regime ($F\ll 1$, $Z\gg 1$) the SY14 approach does lead to accurate predictions for the Reynolds stresses $\langle u' v' \rangle$ due to the waves induced by the small-scale forcing, including at jet maxima and minima where (1.3) breaks down.

2. (ii) That in a typical equilibrated state, both the eastward and westward jets fluctuate around a state which is marginally unstable to barotropic instability (see also Constantinou, Farrell & Ioannou 2014). Jet-scale barotropic waves due to these instabilities are present on both jets, and these have a distinct characteristic structure on each jet. Momentum fluxes due to these waves, which are here disentangled from momentum fluxes due to the small-scale forcing, have a key role in the momentum balance of the jets.

3. (iii) Because the barotropic waves emerge on the jet scale, i.e. they are global, a purely local closure for the jet profile of the type suggested previously is doomed to failure.

The plan of the work is as follows. In § 2, the results of SY14 are reviewed and reconciled with those of WB19, clarifying the limit in which (1.3) is valid. In § 3 CE2 and QL simulations designed to demonstrate the main points are described and the results are analysed. A careful study of the linear stability properties of both the eastward and westward jets is performed to clarify the role of the waves generated by barotropic instability in the momentum balance of the equilibrated flow. In § 4, the model equations are parameterised by replacing the stochastic term by a momentum flux convergence based on SY14. A purely zonal closure is found to work in some special cases, including a monotonic zonal flow profile and a flow which is strongly relaxed to an alternating jet profile. In the general case, however, the parameterisation is found to work only if non-zonal terms are retained in the equations and the barotropic secondary instabilities are allowed to develop. Finally in § 5 conclusions are presented.

2. The QL approximation, the CE2 equations and the SY14 local theory

In this section, analytical and numerical approaches to understanding the behaviour of (1.1) are reviewed and extended.

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 1967; Kraichnan 1967).

2.2. CE2 equations

Following, e.g., Srinivasan & Young (2012) and Marston et al. (2008), the CE2 equations are obtained by introducing an averaging operator $\langle {\cdot } \rangle$, which in different interpretations (see, e.g., Ait-Chaalal et al. 2016; Marston, Qi & Tobias 2019) 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 (2003), Marston et al. (2019) and references therein). The CE2 equations are mathematically equivalent to the SSST equations originally derived and analysed in Farrell & Ioannou (2003, 2007). In each case, the closure models prioritise non-local zonal energy exchanges, pitching them in opposition to turbulent cascade theories Rhines (1975) and PV mixing theories (Dritschel & McIntyre 2008; Dritschel & Scott 2011).

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 (2014) 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 1907). 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 2006). 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 2013; Bakas, Constantinou & Ioannou 2015). 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 (1987) 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 2012; Srinivasan & Young 2012; Bakas & Ioannou 2014; Constantinou et al. 2014), 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 (2017). 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 (2012) 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 (2013, 2019).

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 2017)

(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. Momentum balance in equilibrated jets

In this section the aim is to explain why purely zonal closures such as (2.25) cannot by themselves describe the equilibrium jet profiles in $\beta$-plane turbulence. Evidence from QL simulations, CE2 calculations and fully nonlinear simulations will be used to show that momentum fluxes due to waves arising from jet-scale instabilities also have an essential role in the momentum balance. It has been previously shown (e.g. Farrell & Ioannou 2017) that equilibrated jets in $\beta$-plane turbulence are marginally stable to barotropic instability. In their study upgradient momentum fluxes from the small scale forcing act to sharpen the jets, whereas emergent secondary instabilities generate downgradient fluxes at the westward jet which counter the sharpening. The secondary instability is associated with the Rayleigh–Kuo stability criterion being violated at the westward jet cores. Our aim here is to demonstrate that a similar secondary instability mechanism, due to a distinct unstable mode with a different wavenumber, also leads to downgradient momentum fluxes at the eastward jet. A specific choice of wave forcing is used in the simulations to facilitate the analysis, which is described 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 2012). 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 2010; Young, Read & Wang 2019; Read et al. 2020) 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 2006). 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 2017).

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. 2019) 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 2013.) 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. Deterministic parameterisation of the stochastic forcing

In the previous section it was established that the SY14 closure for the momentum flux (2.17) can accurately predict the direct contribution $\langle u'v' \rangle _D$ due to stochastic forcing, at least in the limit $F \to 0$. In the equilibrated jet flows, however, an equally important contribution $\langle u'v' \rangle _S$ is present due to barotropic instabilities at both the eastward and westward jets. Next, our aim is to address the question of whether, when and how exactly the SY14 closure can practically be used to model stochastically forced zonal flows, by exploring three settings of increasing complexity.

1. (i) A stable monotonic zonal flow: For this class of flows the SY14 closure and, because there are no turning points, the simple approximation (1.3) can be used to predict $\langle u'v' \rangle _D$. Provided the stochastic forcing is not too strong, the equilibrium flow will remain in the stable regime and there will be no secondary instabilities, i.e. $\langle u'v' \rangle _S=0$. This scenario is therefore closely analogous to the ‘vortex condensate’ situation (e.g. Laurie et al. 2014) discussed in the introduction, in which zonal closures have been shown to be successful.

2. (ii) A stable zonal flow with alternating jets: In this scenario a strong ‘radiative’ damping is used to maintain the flow in a stable state consisting of an eastward and westward jet. This is a tougher test of SY14 than the monotonic flow above, because the flow must adjust to rapid changes in $\langle u'v' \rangle _D$ near the jet cusps. However, it remains simpler than the $\beta$-plane jet flows of § 3, because by design the flow is stable and $\langle u'v' \rangle _S=0$. This regime is key to understanding the mechanism by which the directly forced momentum fluxes $\langle u'v' \rangle _D$ lead to the onset of secondary instability.

3. (iii) A typical equilibrated jet flow in $\beta$-plane turbulence: In this flow we know from the previous section that it is essential to include some representation of barotropic instability so that $\langle u'v' \rangle _S$ is captured. Our approach in the following is simply to replace the stochastic forcing with a deterministic term based on SY14 to model $\langle u'v' \rangle _D$, and to permit barotropically unstable waves to emerge spontaneously in the flow, by including non-zonal perturbations in the initial conditions. Therefore, our parameterisation in this case is not purely zonal, but can still be useful when $F \ll 1$, i.e. there is a large-scale separation between the forcing scale and the jet scale, because the latter determines the scale of the barotropic waves.

A practical consideration when implementing the SY14 closure concerns the large gradients in the calculated $\langle u'v' \rangle _D$ which invariably emerge in thin boundary layers near jet cusps. Numerically these boundary layers can always be resolved with a sufficiently high-resolution grid. However, it is often computationally more practical to apply a smoother to the momentum flux $\langle u'v' \rangle _D$ calculated from (2.17). Here we use a kernel smoother which acts on a function $f(y)$ according to

(4.1)\begin{equation} f_{\sigma}(y)=\int_{-\infty}^\infty f(y') \, K_\sigma(y-y') \, {\rm d} y', \end{equation}

where $K_\sigma (y)$ denotes a smoothing kernel with characteristic length scale. In practice, a Gaussian with variance $\sigma ^2$ is used. It turns out that, for a significant range of $\sigma$, this modified closure (denoted SY14$\sigma$ hereafter) actually improves the comparison with our numerical results, for the simple reason that the smoothing in (4.1) can replicate the effect of finite $F$ (see figure 4).

Results from the three different flow scenarios are reported next.

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. 2014).

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 (2008) 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 2016), 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.

5. Discussion and conclusions

The dynamics of jets in geophysical fluids have been described as a ‘wave–turbulence jigsaw’ (e.g. McIntyre 2008). Here, we have made progress in understanding how the pieces of the jigsaw fit together for the canonical model flow of stochastically forced $\beta$-plane turbulence. The key findings are that the momentum fluxes in the equilibrated jet states can be decomposed into a direct contribution due to the forcing $\langle u' v' \rangle _D$, and a secondary flux $\langle u' v' \rangle _S$ due to emergent jet-scale instabilities at both the eastward and westward jets. Note that our finding of jet-scale barotropic instability at the eastward jet does not support the local theory for the momentum balance there proposed recently by WB19. Moreover, $\langle u' v' \rangle _D$ can be accurately modelled throughout the flow (at least for $F \ll 1$) using the analytical approach of SY14, which has been developed and simplified in § 2. Using these insights a parameterised simulation in which the stochastic forcing is replaced by a deterministic term based on SY14, and which accurately captures the dynamical equilibrium of a stochastically forced flow, was presented in § 4.3.

The success of our parameterised simulation depends upon allowing the jet-scale waves responsible for $\langle u' v' \rangle _S$ to emerge spontaneously, due to the barotropic instability of the equilibrated flow. For a complete theory of $\beta$-plane turbulence these waves also need to be parameterised. However, the development of a closure for unstable flows in forced-dissipative equilibrium remains a (classic) unsolved problem (see, e.g., Marston et al. 2008), and is a fundamentally more difficult problem than parameterising the stochastic forcing, because there is no scale separation between the jets and the waves to exploit. Some flows do exist, nevertheless, in which there is no secondary instability (i.e. $\langle u' v' \rangle _S=0$) and these are most easily found when the flow profile is monotonic (see § 4.1), otherwise strong relaxation towards a stable flow is required (as in § 4.2). For such flows a completely zonal closure theory based on SY14 (or SY14$\sigma$) can be successful. An important example is the vortex condensate flows (Laurie et al. 2014; Kolokolov & Lebedev 2016; Frishman 2017; Frishman & Herbert 2018) found in forced-dissipative isotropic 2D turbulence, for which the simplest local closure theory, based on (1.3), gives excellent predictions for the azimuthal velocity profile of the vortex condensate. A key qualitative finding in this work is just how destabilising even weak turbulent forcing can be near flow extrema, particularly as $F \to 0$, as the momentum flux pattern near jet extrema has a strong tendency to sharpen jets, causing a large increase in flow curvature, and subsequent instability due to the PV gradient reversals which are generated.

Even in the absence of a closure for $\langle u' v' \rangle _S$, the ability to parameterise $\langle u' v' \rangle _D$ has the potential to be an extremely useful for modellers wishing to explore jet flows that are subject to very small-scale forcing (i.e. $F \ll 1$), and may even be of practical use in the development of general circulation models for the giant planets. The main reason is that the parameterised model need only resolve the emergent waves which occur due to the instability, for which $10^1\unicode{x2013} 10^2$ model grid points per jet spacing should be adequate. Resolving the actual energy injection scale $k_f^{-1}$ may easily require many orders of magnitude more resolution. However, to end on a note of caution, the SY14 expressions for $\langle u' v' \rangle _D$ presented here are strictly valid only for (near-)steady jets. Extending the results to the full time-dependent situation will be the subject of future work.