## Abstract

Turbulent flows are often observed to be organized into large-spatial-scale jets such as the familiar zonal jets in the upper levels of the Jovian atmosphere. These relatively steady large-scale jets are not forced coherently but are maintained by the much smaller spatial- and temporal-scale turbulence with which they coexist. The turbulence maintaining the jets may arise from exogenous sources such as small-scale convection or from endogenous sources such as eddy generation associated with baroclinic development processes within the jet itself. Recently a comprehensive theory for the interaction of jets with turbulence has been developed called stochastic structural stability theory (SSST). In this work SSST is used to study the formation of multiple jets in barotropic turbulence in order to understand the physical mechanism producing and maintaining these jets and, specifically, to predict the jet amplitude, structure, and spacing. These jets are shown to be maintained by the continuous spectrum of shear waves and to be organized into stable attracting states in the mutually adjusted mean flow and turbulence fields. The jet structure, amplitude, and spacing and the turbulence level required for emergence of jets can be inferred from these equilibria. For weak but supercritical turbulence levels the jet scale is determined by the most unstable mode of the SSST system and the amplitude of the jets at equilibrium is determined by the balance between eddy forcing and mean flow dissipation. At stronger turbulence levels the jet amplitude saturates with jet spacing and amplitude satisfying the Rayleigh–Kuo stability condition that implies the Rhines scale. Equilibrium jets obtained with the SSST system are in remarkable agreement with equilibrium jets obtained in simulations of fully developed *β*-plane turbulence.

## 1. Introduction

Large-scale coherent jets for which no obvious jet-scale forcing can be ascribed are often observed in turbulent flows; examples include the banded winds of the gaseous planets (Ingersoll 1990; Vasavada and Showman 2005) and the earth’s midlatitude jets. This phenomenon of spontaneous jet organization in turbulence has been extensively investigated in observational and in theoretical studies (Williams 1979, 2003; Panetta 1993; Nozawa and Yoden 1997; Huang and Robinson 1998; Lee 2005; Manfroi and Young 1999; Vallis and Maltrud 1993; Cho and Polvani 1996; Read et al. 2004; Galpirin et al. 2004). A comprehensive theory for the formation and maintenance of jets in turbulence must account not only for the emergence of the jets but also for their structure and persistence. The mechanism by which jets are maintained is upgradient eddy momentum flux and, specifically, flux by small-scale eddies widely separated in scale from the jet scale. This mechanism is distinct from mechanisms based on a turbulent cascade, which takes place locally in spectral space. Furthermore, at equilibrium the flux is carried by the continuous spectrum of shear waves as distinct from the discrete jet modes (Huang and Robinson 1998; Panetta 1993; Vallis and Maltrud 1993). This upgradient momentum transfer mechanism has been found to maintain jets both in barotropic forced dissipative models (Huang and Robinson 1998) and in baroclinic two-layer free turbulence models (Panetta 1993; Williams 2003).

This mechanism of interaction between jets and turbulence is quasi-linear in the sense that the continuous-spectrum momentum fluxes can be obtained from a linear model if the jet structure is known (Huang and Robinson 1998). The forcing of the eddies that produce the fluxes can be traced either to explicitly exogenous short-time-scale processes, such as convection (as in the case of the Jovian jets), or to endogenous turbulence generated internally (as in the example of the earth’s midlatitude jets). The wave forcing can be modeled stochastically because these processes have short time scales compared to the jet persistence time scale. It follows that the central component in a theory for jets in turbulence is a method to obtain the structure of the turbulence and the associated fluxes given the jet. This problem has been solved using stochastic turbulence modeling (Farrell and Ioannou 1993a, c, 1996; DelSole and Farrell 1996; Newman et al. 1997; Zhang and Held 1999; DelSole 2004), which provides an analytic method to obtain the distribution of momentum flux arising from the turbulent wave field associated with a given jet structure. Once these fluxes are known the equilibrium states of balance between the large-scale forcing or relaxation to a mean flow and the turbulent momentum flux divergence can be identified, from which a theory for jets in turbulence proceeds.

There are three length scales in barotropic dynamics on a *β* plane: the scale imposed by the boundaries, the Rhines scale *L _{β}* =

*U*/

*β*, and, if a shallow water approximation has been made, the Rossby radius

*L*=

_{R}*gH*/

*f*, in which

*f*is the Coriolis parameter,

*β*is its meridional gradient,

*g*is the acceleration of gravity,

*H*is the depth of the fluid, and

*U*is a characteristic velocity. Emergence of any of these length scales in the solution to turbulence in the barotropic equations cannot be used to argue for any particular physical mechanism for maintaining the jets because the intrinsic scales of the problem are common to all mechanisms and can at best be used to argue consistency with the dynamics. To establish that a physical mechanism is at work in the jet structure problem requires diagnosing that mechanism in the evolution dynamics and at equilibrium. For instance, the nonlocal nature of the momentum transfer from the small eddy scale directly to the jet scale argues against a turbulent cascade arrested at the Rhines scale being the mechanism, despite the often-remarked coincidence of the Rhines radius in observations of jets and in the theory of barotropic

*β*-plane turbulence. The role of

*β*in the jet structure problem has also been advanced as required for breaking the symmetry of the

*f*plane and establishing a preferred direction for the jets, but this role can be assumed by the boundaries and examples of equatorial jets maintained by turbulence, such as the quasi-biennial oscillation (QBO), and eddy-driven jets in nonrotating systems (Plumb and McEwan 1978) suggest that the planetary vorticity gradient may not be essential for the maintenance of turbulent jets.

In this paper we study the structure and dynamics of jets in turbulence by joining the equation governing the stochastically forced perturbation field with the equation for the zonal mean to form a coupled wave–mean flow evolution system. This nonlinear coupled equation system is the basic tool in stochastic structural stability theory (SSST) analysis (Farrell and Ioannou 2003, hereafter FI03).

## 2. Dynamics of the zonally averaged velocity in turbulent flows

### a. Formulation

A theory for jet dynamics in turbulent flow was developed in FI03. We now review the salient ideas of this theory called SSST.

Consider a turbulent rotating barotropic fluid of uniform density confined to the *x–y* plane and let *U*(*y, t*) be the latitudinally dependent (*y*) and time-dependent (*t*) mean velocity averaged in the zonal (*x*) direction. Assume isotropy of the planar geometry is broken either by the flow being on a *β* plane or by the flow being confined to a zonal channel. Let *ψ*(*x, y, t*) be the streamfunction perturbation to this mean flow. The mean flow obeys the Reynolds-averaged equation:

in which overbars denote zonal averaging and the mean flow acceleration arises from an imbalance between the northward zonally averaged eddy vorticity flux , which is equal to the convergence of the northward momentum flux

and linear relaxation to a state of no flow. In the preceding equations the zonal and meridional perturbation velocity are nondivergent with streamfunction *ψ:*

The perturbation vorticity is then

The perturbation streamfunction is written as a Fourier sum of zonal harmonics:

The total eddy vorticity flux in Eq. (1) is the sum of the eddy vorticity flux due to each Fourier component:

where * denotes complex conjugation.

The continuous operators are discretized and the dynamical operators are approximated by finite dimensional matrices in which the state *ψ** _{k}* is represented by a column vector with entries of the complex value of the streamfunction at collocation points.

Each Fourier component of the perturbation streamfunction evolves according to the stochastically forced linear equation:

where 𝗨 is the diagonal matrix formed from the mean velocity *U*(*y*, *t*) at the collocation points and 𝗔* _{k}*(𝗨) is the operator that governs linear dynamics:

with 𝗜 the identity matrix and

the matrix representation of the Laplacian (𝗗^{2} is the matrix representation of *d*^{2}/*dy*^{2}). In writing this perturbation equation we have included in the linear operator parameterization of the nonlinear terms in the complete perturbation equation as stochastic forcing and added dissipation (Farrell and Ioannou 1993a, 1996; DelSole and Farrell 1996; Newman et al. 1997; DelSole 2004). Care must be taken that the structure of this forcing matrix, 𝗙* _{k}*, does not bias the response. We select forcing matrices that are homogeneous in the latitudinal direction with Gaussian autocorrelation function about

*y*proportional to exp[−(

_{i}*y*−

*y*)

_{i}^{2}/

*δ*

^{2}]. Also, the forcing matrices will be the same for all zonal wavenumbers.

The forcing is delta-correlated white noise with zero mean and unit variance:

where 𝗜 is the identity matrix and ^{†} denotes Hermitian transposition. The angle brackets denote an ensemble average, that is, an average over realizations of the forcing.

With these assumptions the time evolution of the ensemble-average covariance matrix of the perturbation field 𝗖_{k} = 〈*ψ*_{k}*ψ*^{†}_{k}〉 can be obtained as follows:

Taking advantage of the explicit expression for the state

we obtain that

under the assumption that the initial state of the system is uncorrelated with the noise, so that 〈*ψ*_{k}(0)** ξ**(

*t*)

^{†}〉 = 0 and the property of the delta functions ∫

^{t}

_{0}〈

**(**

*ξ**s*)

**(**

*ξ**t*)

^{†}〉

*ds*= 𝗜/2. We therefore obtain that the covariance obeys the deterministic equation:

with 𝗤_{k} = 𝗙_{k}𝗙^{†}_{k} the spatial covariance matrix of the stochastic forcing (FI03). For our purposes the matrix 𝗤* _{k}* includes all the relevant characteristics of the forcing; the meridional distribution of the forcing is given by the diagonal elements of this matrix and the autocorrelation function by the rows.

The *j*th component of the ensemble average of the zonally averaged eddy vorticity flux from (6) is

(no sum over *j*).

When the ensemble average of the zonally averaged eddy vorticity flux equals the zonally averaged eddy vorticity flux, that is, when

then the evolution equations for the covariances, (14), the diagnostic equation for the total eddy vorticity flux, (15), and the Reynolds-averaged equation for the mean zonal flow, (1), define a closed nonlinear system for the evolution of the mean flow under the influence of its consistent field of turbulent eddies. This system is globally stable and its stationary solutions are equilibria in which the mean flow is maintained with constant structure by eddy forcing.

This ideal limit in which the ensemble average is equal to the zonal average is of particular interest. In this limit, although the effect of the ensemble-average turbulent fluxes is retained in the solution, the fluctuations associated with the turbulent eddy dynamics are suppressed by the averaging, and the zonally averaged flow dynamics become deterministic. In this limit the mean zonal flow equilibria emerge with great clarity. This ideal limit is approached when the autocorrelation scale of the perturbation field, *l*, is much smaller than the zonal extent, *L*, of the channel. In that case taking the zonal average is equivalent to averaging *N* = *O*(*L*/*l*) statistically independent realizations of the forcing, which, as *N* → ∞, converges to the ensemble average (FI03). Examples of physical systems in which this limit applies include the Jovian upper atmosphere and the solar convection zone, both forced by relatively small-scale convection. When this limit is not approached sufficiently closely the system exhibits stochastic fluctuations about the ideal equilibria. The ideal equilibrium may nevertheless be detected underlying the noisy observations using statistical methods (Koo et al. 2002). One physical system that exhibits such behavior is the polar jet in the earth’s atmosphere. In that case the eddy dynamics are dominated by global wavenumbers 4–10. As a result the number of independent systems that are averaged over a latitude circle is at most order 10, which is too small for the ideal dynamics to be realized closely. Accordingly, while the underlying order is revealed by the ideal dynamics, the zonal jet exhibits stochastic fluctuations about the ideal jet (FI03).

The assumption that the stochastic forcing is independent of the mean and eddy fields limits our analysis. While this is probably a good assumption for the Jovian atmosphere because the forcing is thought to arise from internally generated convection, this is a crude assumption for the earth’s polar jet where the forcing distribution parameterized by 𝗙𝗙^{†} is influenced by the mean flow and eddy amplitude itself. Obtaining the forcing consistent with the turbulence supported by a jet is equivalent to obtaining a closure of the turbulent system. Progress on this problem has been made by DelSole (1999), and his or a similar closure could be used in the present study. However, while it is an attractive avenue for future study, such a closure is not necessary for understanding the basic dynamics underlying the emergence of zonal jets in turbulence. In the interests of simplicity and clarity, for the present we retain a spatially uniform forcing and structure.

### b. Scaling of the equations and boundary conditions

We nondimensionalize the equations, choosing the earth day for the time scale, *T* = 1 day, and length scale *L _{y}* = 10

^{6}m. All variables will henceforth be considered nondimensional. The calculations are performed in a channel 40 units wide in the meridional directions. We impose periodic boundary conditions in the meridional direction. Under spatially uniform forcing the periodic boundary conditions in the meridional direction enforce translational symmetry in the problem. The size of the channel has been selected wide enough to accommodate a number of jets.

The strength of the stochastic forcing, 𝗤, is measured by the forcing density Θ = trace(𝗤)*δ _{y}*/

*L*(where

_{y}*δ*=

_{y}*L*/

_{y}*n*,

*L*is the channel width, and

_{y}*n*is the number of discretization points in the meridional direction). We will generally give the strength of the forcing in dimensional units (mW kg

^{−1}). To obtain the nondimensional value of trace(𝗤)/

*n*divide the values given by

*L*

^{2}

_{y}/

*T*

^{3}. Unless otherwise stated the number of discretization points is

*n*= 128.

## 3. Structural stability of the state of no zonal flow

Consider SSST applied to the case of an eddy perturbation field with a single zonal wavenumber, *k*. The coupled nonlinear system is

With **U** we denote the diagonal elements of the velocity matrix 𝗨. Further consider imposing meridionally homogeneous forcing and periodic boundary conditions. With these assumptions the absolute vorticity flux vanishes, and it follows that the state of no zonal flow, **U** = 0, is a fixed point of the nonlinear Eqs. (17a)–(17b). The question arises whether this equilibrium is stable. Its stability can be determined by considering perturbations *δ***U** and *δ*𝗖 to the equilibrium, (**U*** _{E}* = 0, 𝗖

*), that satisfies the Lyapunov equation 𝗔(𝗨*

_{E}_{E})𝗖

_{E}+ 𝗖

_{E}𝗔(𝗨

_{E})

^{†}+ 𝗤 = 0. The perturbation SSST equations are

in which the summation convention has been used, and where *E* refers to the equilibrium state consisting of the mean velocity 𝗨* _{E}* and eddy covariance 𝗖

*with associated equilibrium absolute vorticity flux*

_{E}_{E}. As discussed in FI03, the nonanalytic projection operator that selects the imaginary part of

*δ*𝗖 requires that the evolution of the real and imaginary part be separated so that the stability equations written symbolically are

with 𝗟 a real linear operator and with the superscripts *R* and *I* denoting the real and imaginary parts (respectively) of the linear operator obtained by linearizing about the equilibrium mean flow **U*** _{E}*. Eigenanalysis of 𝗟 determines the stability of the coupled equilibria.

^{1}

Stability of an equilibrium state (**U*** _{E},* 𝗖

*) depends on the strength of the stochastic forcing, 𝗤, measured by the forcing density Θ = trace(𝗤)*

_{E}*δ*/

_{y}*L*. Consider an equilibrium flow,

_{y}**U**

*, which is linearly stable, in the sense that the operator 𝗔(𝗨*

_{E}*) is stable. If 𝗤 = 0 then 𝗖*

_{E}*= 0 (at equilibrium there are no eddies), so that the stability Eqs. (18a)–(18b) inherit the stability of 𝗔(𝗨*

_{E}*), implying the stability of 𝗟. When forcing is introduced 𝗖*

_{E}*becomes nonzero and mean velocity perturbations can induce changes in the covariance*

_{E}*δ*𝗖, as can be seen from (18a), with the result that the stability of 𝗟 departs from that of 𝗔(𝗨

*), and 𝗟 can become unstable even if 𝗔(𝗨*

_{E}*) is stable. In contrast to previous work on wave–mean flow interaction dynamics in which the eddy variance is not forced and a mode freely develops as a marginal instability (Pedlosky 1972), the presence here of a realistic turbulent background variance allows the whole spectrum of the stable 𝗔(𝗨*

_{E}*) to contribute to the absolute vorticity fluxes and to participate in the destabilization of the zonal flow as revealed through the perturbation stability operator 𝗟. This is an example of a fundamentally new physical mechanism for destabilizing fluid flows, which is essentially emergent from the interaction between turbulence and the mean flow.*

_{E}As an example, in Fig. 1 we present typical results of the maximum growth rate of perturbations from equilibrium **U*** _{E}* = 0 as a function of the strength of the stochastic forcing for the case of eddy forcing at

*k*= 1.885 corresponding to the global zonal wavenumber 12 in a channel of length 40. It can be easily demonstrated that, at Θ = 0, the system (18a)–(18b) has eigenvalues consisting of the pairs

*λ*+

_{i}*λ**

_{j}, where

*λ*are the eigenvalues of the operator 𝗔(𝗨

_{i}*) and in addition,*

_{E}*n*eigenvalues equal to −

*r*, the rate of relaxation of the mean flow to zero. Consequently, in the absence of stochastic forcing the state of zero mean flow is stable if 𝗔(𝗨

_{e}*) is stable. As Θ increases the state of zero flow is progressively destabilized, becoming unstable in the specific example when the critical Θ*

_{E}*= 0.65 mW kg*

_{c}^{−1}is reached. At first this instability grows exponentially, but eventually it reaches a stable equilibrium. As discussed in FI03, the SSST equations are globally stable and their dynamics are remarkably regular despite their high dimensionality. For Θ < Θ

*the zero mean flow is globally stable and any initial mean flow will relax to the state of zero flow.*

_{c}Each eigenfunction of 𝗟 has two components: a perturbation covariance and a perturbation mean flow. The eigenfunctions are harmonic in the meridional direction because of the periodic boundary conditions and the homogeneity of the 𝗟 operator in *y*. The most unstable eigenfunction of 𝗟 is shown in Fig. 2 for forcing amplitude Θ = 1.5 mW kg^{−1} and the other parameters as in Fig. 1. This eigenfunction corresponds to a perturbation with global meridional wavenumber 8 for a channel that extends 40 nondimensional units in the meridional direction. It can be seen from Fig. 2 that the vorticity flux associated with the perturbation covariance part of the eigenfunction (dashed line) is proportional to the velocity perturbation itself, without any phase lag, indicating that this eigenfunction grows or decays without meridional translation and that the corresponding eigenvalue is purely real. Since the forcing Θ is large enough compared to the zonal flow perturbation dissipation, *r _{e}*, the acceleration induced by the absolute vorticity flux dominates the relaxation toward zero flow and the global meridional wavenumber-8 zonal flow perturbation grows exponentially.

The growth rate for this specific example as a function of meridional wavenumber is shown in Fig. 3. The maximum growth rate is attained at the global meridional wavenumber 8 for these parameter values. The wavenumber of maximum growth rate depends on the perturbation dissipation parameter *r*: as the dissipation parameter is reduced the maximum shifts to higher global meridional wavenumber. We have also shown in Fig. 3 the estimated growth rate of meridional modes based on the adiabatic approximation. In the adiabatic approximation the eddy covariance is assumed to be in asymptotic equilibrium with perturbation *δ***U** at all times [lhs of (18a) set to zero]. Estimates of growth rate made using this adiabatic approximation agree identically with the time-dependent eigenfunctions of (18a)–(18b) only at the zeroes of *λ* (this observation is the basis for an efficient method for determining the parameter values for neutral stability of the zero state), but they provide a good estimate of growth rate otherwise.

For eddy forcings confined to single wavenumbers smaller than a critical *k _{c}*, all the eigenfunctions of 𝗟 produce downgradient perturbation vorticity fluxes, and consequently the zero zonal flow is stable for single zonal wavenumber forcing of any amplitude at this and lower values of

*k*. This critical wavenumber is

*k*≈ 1.25 for

_{c}*r*= 0.25 and is shifted to smaller wavenumbers as the dissipation rate

*r*is decreased.

In summary, with *r* fixed the zero zonal mean state is unstable to eddy components with zonal wavenumbers *k* > *k _{c}* and for forcing exceeding a critical strength, but it is stable for sufficiently small zonal wavenumbers at all forcing strengths. The zero state can also be shown to be globally attracting for those wavenumbers for which it is locally stable. Consequently, if a nonzero mean flow is to be maintained, sufficient eddy forcing must be supplied and the eddy field must include appropriately high wavenumber components.

## 4. Equilibration of the SSST eigenstates

Consider equilibration of the unstable zero mean state with eddy field restricted to zonal wavenumber *k* = 1.885 treated in the previous section and with dissipation parameters *r* = 0.25, *r _{e}* = 0.1. Consider first the case with Θ = 1.5 mW kg

^{−1}, a rather small but destabilizing forcing (recall that the critical forcing density to produce instability in this example is Θ

*= 0.65 mW kg*

_{c}^{−1}; cf. Fig. 1). At

*t*= 0 introduce a mean velocity perturbation with a rich meridional spectrum (left panel in Fig. 4); by

*t*= 40 the most unstable global meridional wavenumber-8 mean flow perturbation has emerged and this structure grows exponentially and eventually equilibrates at

*t*≈ 150 (see Fig. 5). This finite-amplitude equilibrium has the same global meridional wavenumber as the most unstable eigenfunction of the perturbation operator 𝗟. The final equilibrated mean flow,

**U**

*, is associated with an asymptotically stable 𝗔(𝗨*

_{E}*) although it satisfies the Rayleigh–Kuo necessary condition for instability as the equilibrated flow has reversals of the sign of*

_{E}*β*−

*U*″

_{E}. It turns out that the flow is unstable at zonal wavenumbers not included in the calculation. This often occurs for quantized zonal wavenumbers: that instability is avoided by adjustment of the mean flow so that unstable wavenumbers are not included in the set of retained zonal wavenumbers. The stability of this equilibrated-state

**U**

*can be investigated by eigenanalysis of the operator 𝗟. The equilibrated state is linearly stable if the neutral eigenvalues that correspond to meridional translations of the equilibrated jets are excepted.*

_{E}Other mean flow equilibria can be constructed by following other unstable meridional mean flow perturbations from the zero mean state to equilibrium. For the chosen strength of forcing these other meridional mean flow perturbations also equilibrate to finite-amplitude states with the same meridional structure as the initial perturbation, and these finite-amplitude equilibria are locally linearly stable. However, their domain of attraction is limited and sufficiently large-amplitude zonal flow perturbations provoke reorganization to the most stable global meridional wavenumber-8 mean zonal state.

## 5. Maintenance of zonal jets on a *β* plane

Consider an eddy field composed of global zonal wavenumbers 1–14 in a reentrant zonal channel of zonal length 40. The corresponding zonal wavenumbers, *k*, take now the discrete values 2*πj*/40, where *j* is the index of the global zonal wavenumber zonal 1 ≤ *j* ≤ 14. Each mode is forced identically with forcing density Θ, and the dissipation parameters are *r* = 0.25 and *r _{e}* = 0.1. As discussed in the previous section, high zonal wavenumber forcing typically produces absolute vorticity fluxes in phase with mean flow perturbations, which destabilize the zero flow equilibrium mean flow, while low zonal wavenumber forcing produces out-of-phase absolute vorticity fluxes that tend to relax perturbations to the zero flow equilibrium. In the presence of the assumed rich spectrum of zonal waves, we have determined [using the adiabatic approximation to the stability Eqs. (18a)–(18b) that can, as previously mentioned, be used to find nonoscillatory neutrally stable states] that the critical Θ

*in the presence of these waves is Θ*

_{c}*= 0.084 mW kg*

_{c}^{−1}(total forcing 14Θ

*), and the most unstable zonal mean flow perturbation eigenfunction has meridional wavenumber 8. The critical forcing required for maintaining a nonzero zonal flow depends crucially on the eddy perturbation dissipation parameter*

_{c}*r*. As

*r*increases, the zonal wavenumber for which the potential vorticity flux becomes upgradient shifts to higher wavenumbers and the critical forcing increases. For example, a 4-fold increase in

*r*leads to a 20-fold increase in the critical forcing, which becomes Θ

*= 2 mW kg*

_{c}^{−1}. As a result, we conclude that the emergence of zonal jets is inhibited as eddy dissipation increases. For the dissipation parameters used in this study the critical forcing for which meridional wavenumber

*n*is destabilized is given in Table 1.

The equilibrated meridional wavenumber-8 jet flow that results from the slightly supercritical forcing amplitude Θ = 0.105 mW kg^{−1} is shown in Fig. 6. The equilibrated flow is stable with the sign of the mean absolute vorticity *β* − *U* ″_{E} everywhere positive (center panel), and the eddy momentum flux everywhere upgradient and proportional to the local shear, − = *DdU _{E}*/

*dy*, where

*D*is the eddy exchange coefficient that has the negative value −1/15. The equilibrated velocity profile is nearly of the sinusoidal form of the most unstable eigenfunction of (18a)–(18b). It is instructive to recover this eigenfunction alternatively by using the above (negative) diffusive approximation of the absolute vorticity flux in (1) to obtain the approximate balance at equilibrium:

implying that

with *U _{m}* determined by higher-order balance in (18a)–(18b). Because Θ = 0.105 mW kg

^{−1}is only a slightly supercritical forcing the amplitude of the equilibrated flow is small and it has nearly the same structure as the eigenfunctions of the stability operator 𝗟.

When the forcing is slightly supercritical the equilibrated mean flow is weakly constrained, and at equilibrium it can have any of the meridional wavenumbers that are unstable for the zero-state flow, although the most unstable meridional wavenumber is favored. We have now obtained a complete understanding of the formation and structure of jets in the limit of weak forcing. However, as the forcing strength increases the mean velocity also increases and the process of equilibration of zonal flow perturbations to the zero mean state becomes more complex. An example of the evolution of the maximum zonal velocity as a function of time for the strongly forced eddy field with Θ = 41.6 mW kg^{−1} is shown in Fig. 7. At first the most unstable eigenfunction of the 𝗟 operator, which happens to also have global zonal wavenumber 8, emerges and grows exponentially to high amplitude. However, this mean flow (first panel in Fig. 8) is not a fixed point of the equilibrium of the coupled system (17a)–(17b): the associated 𝗔 operator is unstable and the absolute vorticity fluxes do not balance the relaxation of the mean flow to zero. A process of readjustment occurs in which the zonal flow moves to a smaller meridional wavenumber, repeatedly adjusting the maximum velocity and the number of jets until equilibrium is reached, which for this case occurs at global meridional wavenumber 2. The operator 𝗔 associated with this equilibrium is stable to infinitesimal eddy perturbations.

The meridional structure of the mean flow that results when the eddy field is forced with various supercritical Θ is shown in Fig. 9. The equilibrated flows are organized in the four panels of Fig. 9 according to the meridional wavenumber of the equilibrated flow. As the forcing increases the maximum wavenumber of the stable equilibrium decreases. For each forcing amplitude there may be multiple stable equilibria with different meridional wavenumbers. When the forcing is strong the mean flows equilibrate near the stability boundary and the equilibrium with the largest meridional wavenumber, *n*_{0}, has the largest basin of attraction, although additional stable equilibria exist with meridional wavenumbers *n* < *n*_{0}. It is remarkable that the finite-amplitude equilibria can each be traced back to an instability of the zero-state mean flow: if the eigenfunction mean flow perturbation with meridional wavenumber *n* is unstable, a stable equilibrated flow with this meridional wavenumber may exist if the meridional wavenumber is less or equal to *n*_{0}. But instability of the zero state to perturbations with meridional wavenumber *n* does not necessarily imply the existence of a finite-amplitude equilibrium with this meridional wavenumber. For example, consider the case of *n* = 1 shown in Fig. 9a. The zero-state mean flow is unstable to *n* = 1 perturbations for Θ > 3.15 mW kg^{−1}. However, no *n* = 1 finite-amplitude equilibria were found for forcing strengths that ranged from the slightly supercritical to Θ = 10.2 mW kg^{−1}.

The structures shown in the bottom right-hand panel of Fig. 8 and in Fig. 9 exhibit the east–west flow asymmetry characterized (for *β* positive) by sharp westerly flow maxima and comparatively broad easterly flow maxima characteristic of jets in strongly forced *β*-plane simulations (cf. Panetta 1993). The equilibrated flows are structurally stable^{2} and also eddy-perturbation stable: the operator 𝗔(𝗨* _{E}*) is stable for all zonal wavenumbers included in the calculation. The maximal growth rate of 𝗔(𝗨

*) as a function of zonal wavenumber,*

_{E}*k*, for some of the equilibrated flows of Fig. 9 are shown in Fig. 10. While these mean flows are stable for zonal wavenumbers 1–14 and with the chosen value of the coefficient of Rayleigh dissipation (the lowest six wavenumbers used in the calculations are indicated by black dots), the equilibrated flows may be unstable to longer waves that are not included in the eddy field. This may be particularly so when the eddy field is strongly forced, as the quantization of global zonal wavenumbers that results from the imposition of periodic boundary conditions excludes the unstable set. The mean absolute vorticity gradient for the equilibrated flow that results with Θ = 5.2 mW kg

^{−1}is shown in the top panel of Fig. 11. There are small regions of reversal in the sign of the mean absolute vorticity gradient that would allow for an instability if there were no dissipation or if other zonal wavenumbers were included in the eddy field structure. The maintained equilibrated mean flow has an approximately sawtooth shape in which regions of almost constant shear are joined together and the momentum flux is everywhere upgradient leading to an eddy exchange coefficient,

*D*(

*y*), such that − =

*D*(

*y*)

*dU*/

*dy*, which is everywhere negative and of approximate value −0.4 for the case shown in Fig. 11 (bottom panel). The meridional structure of the eddy field shown in Fig. 12 is in good qualitative agreement with the structures found in the two-layer simulations of Panetta (1993).

It follows that, in the presence of *β*, the Rayleigh–Kuo stability condition provides the useful if not tight upper bound on the magnitude of the mean flow *U _{max}* and the scale of the jets,

*L*:

When the velocity maximum approaches this upper bound we have the Rhines scaling of the equilibrium jets.

The extent to which the Rhines scaling is observed is shown in Fig. 13, in which the scale of the jet is taken to be *L* = *L _{y}*/(2

*πn*), where

*n*is the global meridional wavenumber of the mean flow. It is evident from this graph that the scaling is observed for the cases of strong forcing. At each meridional wavenumber,

*n*, the equilibria, indicated by squares, progress from right to left as the forcing, Θ, is increased. The equilibrium for the highest forcing is indicated by the square farthest to the left. This equilibrium occurs close to but generally slightly above the stability boundary, and therefore just above the dashed line that indicates the Rhines scaling.

For weak forcing the resulting weak mean flow equilibrates with the mean flow dissipation far from the stability boundary, and for these cases (rightmost squares for each *n* in Fig. 13) the structure and the resulting mean flow are determined not by the Rayleigh–Kuo stability condition but by dissipative equilibrium of the most unstable eigenfunction of the operator 𝗟 about the zero-state mean flow, and therefore these weakly forced equilibria do not follow Rhines scaling. This regime has been observed in 2D *β*-turbulence simulations under weak forcing (cf. Fig. 10 in Vallis and Maltrud 1993).

### a. The limit β → 0

It is interesting to examine jet stability and equilibration in the limit *β* → 0. In the absence of *β* and with periodic boundary conditions the channel becomes isotropic so that no zonal direction can be distinguished. However, we may assume that a zonal direction is imposed by the orientation of the meridional boundaries or by the presence of a small amount of *β*. Using the same parameters and forcing, except that the calculations are performed with *β* = 0, we find that equilibrated zonal jets are still maintained in the presence of adequate forcing. These are shown for a variety of forcing amplitudes in Fig. 14. Note that for the case Θ = 5.2 mW kg^{−1} two equilibrium jets were found. The equilibrated mean flows are both structurally stable and also perturbation stable. They have a sawtooth velocity profile, and as expected do not have east–west asymmetry. The meridional structure of the shear, the momentum flux, and the curvature are shown in Fig. 15 for the equilibrated flow that results with Θ = 10.4 mW kg^{−1}. The eddy momentum flux at equilibrium is everywhere upgradient (has the same sign as the shear) and is nearly constant over the sections of the sawtooth that have constant shear, as would also be the case if band-limited shear waves were forced in an unbounded constant shear flow (Farrell and Ioannou 1993b). We can define a negative equivalent eddy exchange coefficient, *D*, connecting the shear with the momentum flux through the relation − ≈ *DdU*/*dy*. From Fig. 15, *D* is approximately constant with a value of −1. Note that, in the regions of nearly constant shear, the curvature varies linearly with meridional distance indicating that the velocity profile to first order in *r _{e}* is

consistent with a local balance between perturbation momentum flux divergence, −*Dd*^{2}*U _{e}*/

*dy*

^{2}, and mean flow dissipation

*r*≈

_{e}U_{E}*αr*.

_{e}y## 6. Conclusions

Steady coherent large-scale zonal jets are commonly observed in both rotating and nonrotating turbulent fluids. These jets emerge spontaneously from the turbulence, and once established are often apparently stable. Although laminar flow instabilities and turbulent cascades have been considered in the search for a detailed theoretical explanation of jet emergence from turbulence, observations suggest that a transfer of momentum directly to the jet scale by small-scale shear waves that are nonlocal in wavenumber is responsible. Accepting that the primary physical mechanism maintaining the jets is this systematic eddy momentum flux convergence by small-jet relative-scale eddies, a detailed dynamical explanation for the required systematic organization of the eddy fluxes has been lacking. In this paper we provided this theory, building on results from linear stochastic turbulence modeling and specifically SSST. One advantage of SSST is that it allows an analytic solution for the eddy covariance in statistical equilibrium with the jet so that feedback between the jet and the turbulence can be diagnosed in detail. Central to this model is coupling the evolution equation for the eddy forcing obtained from the statistical equilibrium stochastic turbulence model with the evolution equation for the zonal jet. This nonlinear coupled set of equations exhibits robust and relatively simple behavior. Using this model we find that the state of zero zonal flow is a stationary state of the coupled system, but that this stationary state is unstable to zonal perturbations if the turbulence is sufficiently strong, and these zonal perturbations evolve into jets that grow and adjust in structure until reaching finite-amplitude equilibrium.

With weak forcing the equilibrated jets do not follow the Rhines scaling; however, if the flow is sufficiently strongly forced, equilibrium occurs near a modal stability boundary for the jet. This modal stability boundary coincides with the intrinsic length scale in the problem that is the Rhines scale, but it is unrelated to arrested turbulent cascades. Instead, the stability boundary is related to the Rhines scale through the Rayleigh–Kuo necessary condition for instability because jet growth is due to transports produced by the smoothly distributed continuous spectrum of shear waves and cannot be in equilibrium with the local critical-layer transport associated with weakly unstable modes. For this reason strongly forced unstable jets grow and reorganize until reaching equilibrium near the stability boundary and therefore near the Rhines scale. Perturbation of the zero state reveals the coupled system to be unstable for a range of zonal and meridional wavenumbers but with a maximum growth rate at a preferred scale, which is a function of the strength of the stochastic forcing and the eddy dissipation. The finite-amplitude equilibria reached starting from these instabilities are found generally to be individually stable but with a most stable and therefore a preferred meridional wavenumber that may be expected to eventually prevail as the statistically preferred structure. This preferred equilibrium jet structure for geophysical typical values of *β* exhibits pronounced asymmetry, with the westerly jet being sharper than the easterly jet because the positive curvature of the westerly jet is limited by the Rayleigh stability requirement. Moreover, the structure of the jets at finite amplitude can be understood from the upgradient transport of momentum proportional to shear obtained in the analytic solution for band-limited excitation of free shear flow (Farrell and Ioannou 1993b). This ubiquitous instability of barotropic turbulence to jet formation can be viewed as arising because the equilibrium stochastic spectrum of wave modes in a turbulent flow is always available to be organized by a perturbation zonal jet to produce a flux proportional to the jet amplitude. This local upgradient flux due to the continuous spectrum in shear flow results in an exponential jet growth rate because it produces a flux proportional to jet amplitude. This is a new mechanism for destabilizing turbulent flows that is essentially turbulent in nature. As a further consequence, at finite amplitude the jet sides have nearly constant shear with only that small variation in velocity that is necessary, consistent with the proportionality of flux and shear, to balance the imposed relaxation of the jet to zero. For nearly constant shear flow regions on the jet flank this requires a deviation from constant shear proportional to *y*^{3} to produce the requisite linearly proportional flux divergence.

The wavenumber of a weakly forced jet is that of the most unstable eigenmode of the linearized SSST system, and the amplitude is obtained from the balance between upgradient eddy momentum flux proportional to mean shear flux and dissipation.

This work provides a detailed physical theory for jet structure in barotropic turbulence. Baroclinic turbulence is similar to barotropic turbulence in that momentum flux convergence into the upper-level jet plays a similar role in the maintenance of both baroclinic and barotropic jets, but in baroclinic turbulence there is the possibility for heat flux and related mean meridional flows to also play a role in redistributing momentum. In addition, the forcing in baroclinic jets is primarily intrinsically related to baroclinic eddy growth processes in the jet rather than extrinsically imposed by processes external to the jet dynamics. The dynamics of baroclinic jet formation in quasigeostrophic turbulence will be examined using the methods of structural stability in a separate paper.

## Acknowledgments

This work was supported by NSF ATM-0432463. Petros Ioannou was supported in part by an ELKE grant from the University of Athens. We acknowledge the comments of the two anonymous reviewers.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

## Footnotes

C*orresponding author address:* Brian Farrell, Department of Earth and Planetary Sciences, Harvard University, Geological Museum 452, 24 Oxford Street, Cambridge, MA 02138. Email: farrell@deas.harvard.edu

^{1}

The size of the matrix 𝗟 is *O*(n^{2} × *n*^{2}) where *n* is the number of discretization points in the meridional direction, and for reasonable discretizations eigenanalysis requires exploiting the sparseness of 𝗟. The calculations to be presented were performed with 128 points.

^{2}

As previously discussed there is a zero eigenvalue of the stability operator corresponding to meridional translation of the mean flow that is not pertinent to flow stability. This translational mode can be suppressed by imposing a small amount of spatially varying Rayleigh friction that breaks the meridional homogeneity of the flow.