Abstract

Zonostrophic instability leads to the spontaneous emergence of zonal jets on a β plane from a jetless basic-state flow that is damped by bottom drag and driven by a random body force. Decomposing the barotropic vorticity equation into the zonal mean and eddy equations, and neglecting the eddy–eddy interactions, defines the quasilinear (QL) system. Numerical solution of the QL system shows zonal jets with length scales comparable to jets obtained by solving the nonlinear (NL) system.

Starting with the QL system, one can construct a deterministic equation for the evolution of the two-point single-time correlation function of the vorticity, from which one can obtain the Reynolds stress that drives the zonal mean flow. This deterministic system has an exact nonlinear solution, which is an isotropic and homogenous eddy field with no jets. The authors characterize the linear stability of this jetless solution by calculating the critical stability curve in the parameter space and successfully comparing this analytic result with numerical solutions of the QL system. But the critical drag required for the onset of NL zonostrophic instability is sometimes a factor of 6 smaller than that for QL zonostrophic instability.

Near the critical stability curve, the jet scale predicted by linear stability theory agrees with that obtained via QL numerics. But on reducing the drag, the emerging QL jets agree with the linear stability prediction at only short times. Subsequently jets merge with their neighbors until the flow matures into a state with jets that are significantly broader than the linear prediction but have spacing similar to NL jets.

1. Introduction

Zonal flows are banded, anisotropic, weakly fluctuating alternating jets that form spontaneously and persist indefinitely in an otherwise turbulent plasma or planetary fluid (Diamond et al. 2005; Vasavada and Showman 2005). The subject started with Rhines’ discovery that freely evolving barotropic β-plane turbulence transfers energy into zonal shear modes with zero frequency (Rhines 1975). Also in 1975, experiments by Whitehead showed that forcing, without the exertion of azimuthal torque, in a rapidly rotating basin produces prograde jets; in this context the curved upper surface provides an analog of the β effect. We follow Galperin et al. (2006) in referring to the development and persistence of these anisotropic planetary flows as “zonation.” Williams (1978) showed that zonation occurs in statistically steady forced-dissipative flows on the sphere and proposed this as an explanation of the banded structure of the planetary circulations of Jupiter and Saturn.

Figure 1 shows a typical example of fully developed, forced and dissipative zonation obtained by numerical solution of (3) below. The main features of the statistically steady flow, such as the sharp eastward jets, the broader westward return flows, and the sawtooth relative vorticity, are familiar from many earlier studies of statistically steady, stochastically forced, dissipative β-plane turbulence in doubly periodic geometry (Danilov and Gurarie 2004; Danilov and Gryanik 2004; Maltrud and Vallis 1991; Smith 2004; Vallis and Maltrud 1993) and on the sphere (Williams 1978; Nozawa and Yoden 1997; Huang and Robinson 1998; Scott and Polvani 2007).

Fig. 1.

Nonlinear (NL) zonal jets. (left) A snapshot of the zonally averaged velocity U(y, t) obtained from a solution of (3) in a doubly periodic domain 2πL × 2πL with kfL = 32, where kf is the dominant wavenumber of the forcing ξ. (right) A snapshot of the vorticity ζ, with overlaid zonally averaged vorticity −Uy(y, t) (solid white curve). The parameter values for this run are μ* = 0.018 24 and β* = 1.0. The snapshot is at 2μt = 25 with spinup from rest.

Fig. 1.

Nonlinear (NL) zonal jets. (left) A snapshot of the zonally averaged velocity U(y, t) obtained from a solution of (3) in a doubly periodic domain 2πL × 2πL with kfL = 32, where kf is the dominant wavenumber of the forcing ξ. (right) A snapshot of the vorticity ζ, with overlaid zonally averaged vorticity −Uy(y, t) (solid white curve). The parameter values for this run are μ* = 0.018 24 and β* = 1.0. The snapshot is at 2μt = 25 with spinup from rest.

To establish the notation used in this study, we start with the equations of motion for a forced and dissipative barotropic flow u = (u, υ):

 
formula
 
formula

where f = f0 + βy is the β-plane Coriolis frequency. The flow is energized by a solenoidal (incompressible) force generated by the function a(x, y, t). There is no loss of generality in taking the force to be solenoidal: any compressive component of the external force is balanced by the pressure gradient. Damping is provided by a combination of drag μ and hyperviscosity νn (with n = 4 in numerical simulations, and n = 1 in development of theory).

The incompressible velocity field (2) admits a streamfunction ψ(x, y, t) with (u, υ) = (−ψy, ψx), and relative vorticity ζ = ψxx + ψyy. Eliminating the pressure from (1), one obtains the β-plane vorticity equation

 
formula

The vorticity forcing ξ on the right-hand side of (11) is the curl of the solenoidal force in the momentum equation—that is,

 
formula

We assume that the forcing, a in (1) and ξ in (3), is a rapidly decorrelating, isotropic, spatially homogeneous, random process. Thus energy and enstrophy are injected into a narrow band of wavenumbers centered on a “forced wavenumber” kf (see  appendix A for details of the implementation). This model of exogenous stochastic forcing, first proposed by Lilly (1969), is now a standard protocol used in many barotropic and shallow-water studies of forced-dissipative zonation. The physical interpretation of the forcing and the choice of its spatial structure vary somewhat in literature. Considering ξ to be a representation of baroclinic eddies, Williams (1978) chose the forced wavenumbers to lie in a narrow rectangular band, with the zonal extent of the band equal to the baroclinic deformation radius. Scott and Polvani (2007) and Smith (2004) interpreted the rapidly decorrelating, narrowband, isotropic forcing as a model of small-scale three-dimensional convection. Another possibility is that ξ is a representation of the bubble-cloud forcing used by Whitehead (1975) in the laboratory. Below, in the discussion surrounding (10) we give yet another interpretation of ξ.

We have found no studies that establish any particular forcing protocol as being a reasonable physical representation of three-dimensional small-scale eddies acting on a barotropic flow. However, despite the two strong modeling approximations, namely quasi-geostrophy and the choice of the forcing, some features of Jovian jets, such as the jet width, are approximately captured by simplified models (Smith 2004; Vasavada and Showman 2005). A different model forcing in Showman (2007) uses nonsolenoidal physical space mass forcing to represent moist convection in a shallow water system. Despite the different choice of forcing, Showman’s results on planetary zonal jets are broadly consistent with those obtained by Smith (2004). In light of this fact, and since the barotropic quasigeostrophic (QG) system cannot represent mass forcing, we do not address these issues further.

A common theme in all the studies mentioned above is a separation of scales between the forcing length scale and the width of the emergent jets. Indeed, the spacing of the jets in Fig. 1 is significantly greater than , which is an indication either of the inverse cascade or of a spectrally nonlocal transfer of energy (Huang and Robinson 1998).

A striking feature of β-plane zonation is that the translational symmetry, yy + a, of the equation of motion (3) is spontaneously broken: the locations of the eastward maxima in Fig. 1 are an accident of the initial conditions and of the random number generator used to create ξ. But after the jets form, they remain in the same position, apparently forever. Once these robust quasi-steady jets are in place, their dynamics can be discussed in mechanistic terms using concepts such as potential vorticity (PV) mixing, the resilience of transport barriers at the velocity maxima, radiation stress, and shear straining of turbulent eddies (Rhines and Young 1982; Dritschel and McIntyre 2010). But the primary question addressed here is why the jets form in the first place, given that ξ does not select particular locations. Following earlier investigations of this phenomenon (Farrell and Ioannou 2007; Manfroi and Young 1999), we show that zonation can be understood as symmetry-breaking instability of an isotropic, spatially homogeneous, and jetless β-plane flow.

In section 2 we introduce the eddy-mean decomposition and discuss a statistical method, previously used by Farrell and Ioannou (1993b, 2003, 2007), Marston et al. (2008), and Tobias et al. (2011), which is the basis of our linear stability analysis of zonostrophic instability. This method amounts to forming quadratic averages of the equations of motion and then discarding third-order cumulants. Farrell and Ioannou (2003, 2007) refer to this method as stochastic structural stability theory (SSST), while Marston et al. (2008) call it the second-order cumulant expansion, or CE2. SSST and CE2 are completely equivalent, and only one name is required. We have therefore adopted the more descriptive CE2 terminology of Marston et al. (2008).

In section 3 we present a physical space reformulation of CE2, which has analytic advantages over earlier numerically oriented formulations. Within the context of CE2, section 4 provides a complete analytic description of zonostrophic instability obtained by linearizing around an exact isotropic and homogeneous solution with no jets. As in Farrell and Ioannou (2007), zonation is understood as a linear instability of CE2: part of the linearly unstable eigenmode is a zonal flow. This linear stability problem is characterized by two control parameters, a nondimensional drag parameter μ* and a nondimensional planetary parameter β*, and we determine the CE2 zonostrophic stability boundary in the (β*, μ*)-parameter plane. An important property of CE2 zonostrophic instability is that the most unstable wavenumber, which determines the meridional scale of the exponentially growing jets, is well away from zero. Because the instability unfolds around a nonzero wavenumber, CE2 zonostrophic instability is not properly a negative-viscosity instability. This point is reinforced in section 5 by showing that the CE2 eddy viscosity is identically zero. Section 6 is a comparison between the analytic results and direct numerical simulations of the nonlinear system. Section 7 is the discussion and conclusions. The more technical aspects of the paper are in five appendices.

2. The eddy-mean decomposition and quasilinear dynamics

We use an eddy-mean decomposition

 
formula

where the overbar denotes a zonal average; we also denote the zonal mean velocity as . Applying this average to (3) results in the zonal mean momentum equation

 
formula

and the eddy vorticity equation

 
formula

In (7), the eddy–eddy nonlinearity (EENL) is

 
formula

In addition to the zonal average, the overbar includes a running time average over a short interval so that does not appear on the right-hand side of (6). In presenting equations subsequently used to obtain analytic results, we use n = 1 for the viscosity.

a. Quasilinear dynamics

The main results in this paper are obtained with a quasilinear (QL) system, which is defined by taking

 
formula

in (7). Figure 2 shows a QL solution at the same parameter values as the fully nonlinear (NL) solution as Fig. 1. Because of the coupling between the mean and the eddies, the QL system is nonlinear, and Fig. 2 shows that QL dynamics still results in the spontaneous formation of quasi-steady zonal jets.

Fig. 2.

Quasilinear (QL) zonal jets. (left) A snapshot of the zonally averaged velocity U(y, t) obtained by integrating the QL system (6), (7), and (9). (right) A snapshot of the QL vorticity ζ, with overlaid zonally averaged vorticity −Uy (solid white curve). The parameters for this run are the same as the nonlinear solution in Fig. 1 (i.e., μ* = 0.0182, β* = 1, and kfL = 32). The snapshot is at 2μt = 40 after spinup from rest.

Fig. 2.

Quasilinear (QL) zonal jets. (left) A snapshot of the zonally averaged velocity U(y, t) obtained by integrating the QL system (6), (7), and (9). (right) A snapshot of the QL vorticity ζ, with overlaid zonally averaged vorticity −Uy (solid white curve). The parameters for this run are the same as the nonlinear solution in Fig. 1 (i.e., μ* = 0.0182, β* = 1, and kfL = 32). The snapshot is at 2μt = 40 after spinup from rest.

Comparing the left panels in Figs. 1 and 2, one sees that the QL jets are faster and wider than NL jets, and the jet profiles are different: QL jets are distinctly more east–west symmetric than NL jets. Nonetheless, we show in section 6 that the QL jets in Fig. 2 do have a small east–west QL asymmetry, and at other points in the (β*, μ*)-parameter space, QL jets are strongly east–west asymmetric.

Because the QL jets are faster, the QL system is more zonostrophically unstable than the NL system. In Figs. 1 and 2, quasi-steady jets evolve spontaneously from an initially jetless state, as shown in the Hovmöller diagram of the zonal mean flow in Fig. 3. Comparing Figs. 3a and b shows that the QL system has significantly longer adjustment times than the NL system.

Fig. 3.

(a) Hovmöller diagram of the zonal mean velocity U(y, t) obtained by solution of the full NL system in (3). (b) Hovmöller diagram of the zonal mean velocity U(y, t) obtained by solution of the QL system. (c) A comparison of the zonal mean energy fraction, zmf(t) defined in (78), for QL and NL runs. The time-averaged fractions are 〈zmf〉NL = 0.3 and 〈zmf〉QL = 0.51. This figure shows the time evolution of the runs in Figs. 1 and 2.

Fig. 3.

(a) Hovmöller diagram of the zonal mean velocity U(y, t) obtained by solution of the full NL system in (3). (b) Hovmöller diagram of the zonal mean velocity U(y, t) obtained by solution of the QL system. (c) A comparison of the zonal mean energy fraction, zmf(t) defined in (78), for QL and NL runs. The time-averaged fractions are 〈zmf〉NL = 0.3 and 〈zmf〉QL = 0.51. This figure shows the time evolution of the runs in Figs. 1 and 2.

O’Gorman and Schneider (2007) made the QL approximation (9) in an atmospheric general circulation model and showed by comparison with the full nonlinear version of the model that several important features of the flow are unaffected by complete removal of the eddy–eddy nonlinearity as in (9). Comparing Figs. 1 and 2 we reach a similar conclusion for the more idealized model studied here. This preliminary conclusion is supported by a detailed comparison between NL and QL solutions in section 6.

There are several ways of motivating QL dynamics. The QL system conserves both energy and enstrophy and has the same zonal mean equation and symmetries as the NL system. Thus, arguments based on quadratic integral invariants apply equally to the QL and the NL system (Salmon 1998). Nonetheless, because EENL is discarded, the QL system cannot exhibit a true Batchelor–Kraichnan inverse energy cascade: in the QL model all nonlinear interactions require participation of the zonal mean flow. Because U(y, t) has a larger length scale than the eddy field, all these nonlinear QL interactions are spectrally nonlocal. Fig. 2 shows that the spectrally local Batchelor–Kraichnan inverse cascade is not necessary for zonation.

PV is not materially conserved by the QL system, and consequently nonquadratic functions of PV are not conserved by the nonlinear terms remaining in QL. Thus, Fig. 2 also shows that strict material conservation of PV is not necessary for zonation.

Thus, at the most basic level, the QL system is instructive as an indication of the physically essential processes necessary for zonation.

b. Stochastic closure versus cumulant expansion

A main motivation of the QL system is that using the statistical method pioneered in meteorology by Farrell and Ioannou (1993b, 2003, 2007) one can compute important average quadratic properties of the QL flow, such as Reynolds stress and the eddy enstrophy and energy. However, rather than (9), Farrell and Ioannou (2007) adopt a stochastic closure, which amounts to replacing the eddy–eddy nonlinearity with a combination of random forcing and dissipation:

 
formula

(see also DelSole 2001). The intent of (10) is that the random forcing ξEENL(x, t) and the dissipation μEENLζ′ should be chosen to match the evolution of the NL system. The terms in (10) then augment the exogenous forcing and dissipation on the right-hand side of (7).

However, there is probably no reliable a priori method of determining the right-hand side of (10). Heeding the principle to first do no harm, we prefer the QL alternative (9). This has the advantage that one can then make a specific comparison between QL and NL solutions (e.g., as in Figs. 1 and 2) and assess the role of EENL.

Our point of view, which follows Marston et al. (2008) and Tobias et al. (2011), is to regard the QL system as an approximation to the NL system. In fact, (9) in tandem with the method of Farrell and Ioannou, is precisely the second-order cumulant expansion CE2 of Marston et al. (2008). It is from this perspective that in section 3 we develop a physical-space formulation of CE2, which is suitable for analytic solution.

3. Dynamics of correlations: CE2

In the QL approximation, the eddy vorticity equation can be written as

 
formula

where is the Rayleigh–Kuo operator:

 
formula

In this section we obtain a closed deterministic evolution equation for the two-point correlations function of vorticity ζ′ and streamfunction ψ′. This correlation equation [(23) below] is coupled to the evolution of the zonal mean flow via the Reynolds stresses, and the Reynolds stresses can be obtained by evaluating derivatives of the correlation function at zero separation. Thus one obtains the zonal mean evolution equation in (34) below.

a. Correlation functions: Kinematics

We assume that the external forcing, a(x, t) in (1) and ξ(x, t) in (11), has two-point, two-time correlation functions of the form

 
formula
 
formula

where the overbar above denotes an ensemble average. The dependence of the spatial correlation functions A and Ξ only on the difference

 
formula

indicates that the forcing is zonally homogeneous. We do not assume (yet) that the forcing is isotropic and meridionally homogenous.

Because derivatives commute with the ensemble average, relation (4) implies that A and Ξ are related by

 
formula

where the Laplacian acting on the coordinates of point n is

 
formula

Notice that that we have changed notation: undecorated x in (15) is the zonal difference coordinate. We also use the shorthand , etc. Strictly speaking, we should decorate all the variables in (11) with the subscript n = 1 or 2 to explicitly indicate whether we refer to the eddy vorticity equation at the point x1 = (x1, y1) or at the point x2 = (x2, y2). We forbear from doing so.

We assume “ergodicity” so the overbar is also equivalent to the zonal average of a single realization. We desire the single-time two-point correlation functions

 
formula

and

 
formula

The analog of (16) is

 
formula

Given the streamfunction correlation Ψ(x, y1, y2, t), one can obtain the velocity correlation tensor as

 
formula

Because the choice of denoting one point as x1 and the other as x2 is arbitrary, all correlation function have an important “exchange” symmetry

 
formula

and likewise for A, , Ξ, etc.

b. Correlation functions: Dynamics

It follows from (11) and (14) that the correlation functions evolve as

 
formula

where the Rayleigh–Kuo operator at point n is .

To derive (23), multiply the equation for by and add this to the equation multiplied by . The sum is then ensemble averaged, and after the average all fields depend on x1 and x2 only through the combination x = x1x2. Because of this zonal homogeneity . Thus, for example,

 
formula

A crucial simplification is that the forcing is rapidly decorrelating in time, as expressed by the δ(t1t2) in (14). Considerations summarized in  appendix B (amounting to a simple proof of Ito’s formula) show that

 
formula

The result above is the origin of the first term on the right-hand side of (23).

c. Collective coordinates

As alternatives to y1 and y2, there are advantages in using the “collective coordinates”

 
formula

Eventually we will restrict attention to homogenous and isotropic forcing, and at that point we take Ξ in (14) to be a function only of the two-point separation

 
formula

Collective coordinates are then essential for analytic progress.

In terms of y and , the Laplacians are

 
formula

where is the “separation” Laplacian. Thus, for instance,

 
formula
 
formula

Using the coordinates in (26), the correlation equation (23) becomes

 
formula

where now and and .

d. The zonal mean flow equation

One advantage of collective coordinates is that mean square quantities, such as the enstrophy, are obtained by evaluating correlation functions at zero separation [i.e., by setting (x, y) = 0]. For example, if one possesses then the eddy enstrophy is .

A key result, obtained by evaluating

 
formula

at (x, y) = 0, is that the Reynolds stress is

 
formula

Thus, the mean flow equation (6) can be written as

 
formula

The mean flow equation (34), coupled with the correlation equation (31), is a closed system for the ensemble-averaged properties of QL dynamics.

4. Zonostrophic instability of a spatially homogeneous and isotropic base-state flow

a. The spatially homogeneous basic state

We now suppose that the forcing is statistically homogenous and isotropic, namely that the correlation function of the forcing has the particular form

 
formula

where r is the two-point separation defined in (27). Because Ξ does not depend on , there is a simple exact solution to (31) and (34). This solution is spatially homogeneous and isotropic and has no mean flow, U = 0. With these simplifications the correlation equation (31) collapses to

 
formula

The subscript H emphasizes that is spatially homogeneous (i.e., independent of ). The streamfunction correlation function ΨH(x, y) can be obtained from by solving . It is remarkable that in (36) is independent of β: an isotropic and spatially homogenous forcing drives an isotropic and spatially homogeneous flow, despite the anisotropy of Rossby wave propagation.

We now apply the Fourier integral theorem,

 
formula
 
formula

to (36). We use the notation

 
formula

so that after the Fourier transform , and the streamfunction spectrum is related to the forcing spectrum by

 
formula

We emphasize that in (40) is not singular as h → 0. To see this, we recall (16), which in this homogeneous and isotropic case implies that Ξ = ∇4A, and therefore . In terms of then, the streamfunction spectrum in (40) becomes

 
formula

Since a(x, t) is stationary, the spectra and are finite as h → 0 (provided that μ ≠ 0).

Later we will need two integral constraints on the vorticity forcing correlation function Ξ:

 
formula

and also

 
formula

These follow from Ξ = ∇4A, and the assumption that the correlation function A(r) decays faster than r−1 as r → ∞. The constraints above are satisfied by the standard forcing protocol described in  appendix A, which has zero spectral density around h = 0.

b. The dispersion relation of inviscid and isotropic flow

The linear stability of the spatially homogeneous solution in (41) is determined by imposing small initial disturbances and examining evolution in time. The perturbation variables are added to the base-state variables, (0, ΨH, ) and substituted into (31) and (34). The total “flow,” with mean and imposed small perturbations, is

 
formula

where m is the meridional wavenumber of the disturbances and s is the growth rate, with growing perturbations corresponding to ℜ(s) > 0. Retaining terms linear in the perturbation variables , one has the equations governing the evolution of small perturbations to the homogeneous solution. The details of the subsequent solution are in  appendix C, and a main result of that analysis is the dispersion relation

 
formula

where is the forcing spectrum in (41), and the function Q(χ, n) is defined by the angular integral

 
formula

The dispersion relation (45) applies to the special case of isotropic forcing, A = A(r) and ν = 0; a more general expression of the dispersion relation is in  appendix C.

Dr. George Carnevale has shown that the dispersion relation in (45) and (46) is also obtained from (5.13) in Carnevale and Martin (1982). The field-theoretic approach of Carnevale and Martin (1982) is different from the approximation used to obtain the CE2 system in (31) and (34); for instance, CE2 contains terms such as , which Carnevale and Martin (1982) consider to be fourth order in wave amplitude and therefore negligible. However, after linearization of CE2 around a basic state with U = 0, these terms are neglected. Therefore, the linearized version of CE2 in this section is equivalent to the weak-turbulence limit (5.13) in Carnevale and Martin (1982). This consistency provides confidence in (45) and (46).

c. Ring forcing

In most previous investigations of zonation, the forcing is limited to an annulus of wavenumbers in Fourier space. Typically the annulus of forced modes has a mean radius h = kf and has thickness 2δkkf. This is the “narrow-band forcing” described in  appendix A. We idealize this choice further by considering “ring” forcing corresponding to the limit δk → 0. In other words, we consider a random flow, driven isotropically by injecting energy on the circle h = kf in wavenumber space. This corresponds to

 
formula

where J0 is the Bessel function of order zero. Notice that . With ν = 0—as we assume in (45)—the spatially homogeneous base-state solution in (41) is

 
formula

The parameter ɛ above, with dimensions of watts per kilogram, is the rate of working of the force that sustains the base-state (48) flow against dissipation.

With δ(hkf) in (47), the h integral in (45) is trivial. Before proceeding, however, it is convenient to write the various parameters in the nondimensional form using the length scale and time scale . These scales lead to the control parameters

 
formula

The nondimensional wavenumber and growth rate are

 
formula

The zonostrophic dispersion relation in nondimensional variables is then

 
formula

with the function Q defined in (46). We now lighten the notation by dropping the asterisk on nondimensional variables m and s. We have obtained the growth rate by solving (51) numerically for s = sr + isi. This numerical solution indicates that modes with sr > 0 are found only if 0 < m2 < 1, and these unstable modes have si = 0. We have been unable to obtain a satisfactory nonnumerical proof of these two important properties of zonostrophic instability.

Figure 4 shows some examples of the growth rate s(m) plotted as a function of m for various values of β*, with μ* = 0.15 in all cases. If s > 0 for some values of m (e.g., β* = 0.15 and 1 in Fig. 4), then the homogeneous flow is unstable and zonal jets will grow from very small initial amplitude. Also shown in Fig. 4 are two marginally stable situations β* = 0.0634 and β* = 2.571. These are defined by the condition that the most unstable disturbance has s = 0:

 
formula

A main conclusion resulting from our analysis is that zonostrophic instability is suppressed if β* is either too large or too small (e.g., in Fig. 4 the flow is stable if β* > 2.571 or if β* < 0.0634).

Fig. 4.

The growth rate s as a function of m for μ* = 0.15 and five values of β* indicated on the curves. The variables in this figure are nondimensionalized according to (49) and (50). These modes have si = 0 (i.e., s is real). The curves β* = 2.571 and 0.0634 correspond to the marginally stable situation defined by (52).

Fig. 4.

The growth rate s as a function of m for μ* = 0.15 and five values of β* indicated on the curves. The variables in this figure are nondimensionalized according to (49) and (50). These modes have si = 0 (i.e., s is real). The curves β* = 2.571 and 0.0634 correspond to the marginally stable situation defined by (52).

The marginal stability condition in (52), which is equivalent to the requirements

 
formula

defines a “critical curve” in the (β*, μ*) parameter plane. This curve, is shown in the upper panel of Fig. 5. The solution is linearly unstable in the region below the critical curve. The peak of the critical curve is . This peak defines the “most unstable” point in the (β*, μ*) parameter space (i.e., the largest value of drag μ* at which the homogeneous solution loses stability). The lower panel of Fig. 5 shows the wavenumber mc(β*) of the incipient instability [i.e., the wavenumber determined by simultaneously satisfying the two equations in (53)].

Fig. 5.

(a) The critical curve (solid line); linear zonostrophic instability occurs in the region below the critical curve. (b) The wavenumber on the critical curve (solid line; i.e., the most linearly unstable wavenumber). Asymptotic approximations for the critical curve and the most unstable wavenumber based on β* ≫ 1 (dash-dotted) and β* ≪ 1 (dashed) are shown in both panels.

Fig. 5.

(a) The critical curve (solid line); linear zonostrophic instability occurs in the region below the critical curve. (b) The wavenumber on the critical curve (solid line; i.e., the most linearly unstable wavenumber). Asymptotic approximations for the critical curve and the most unstable wavenumber based on β* ≫ 1 (dash-dotted) and β* ≪ 1 (dashed) are shown in both panels.

d. Approximations to the neutral curve with large and small β*

Also shown in Fig. 5 are analytic approximations in the complementary limits β* ≪ 1 and β* ≫ 1. These results are obtained via asymptotic analysis of the integral Q(χ, n) in (46), and simplification of the dispersion relation [(51); see  appendix E]. If β* ≪ 1 then the critical curve is

 
formula

with wavenumber

 
formula

In the complementary limit β* ≫ 1, the approximation to the critical curve is

 
formula

with wavenumber

 
formula

The lower panel of Fig. 5 shows that linear zonostrophic instability is spectrally nonlocal only in the limit β* → ∞: in that case the most unstable wavenumber is much less than the forced wavenumber kf, implying a scale separation between the scales at which energy is injected and the scale at which jets initially form. In the other limit β* → 0 the linearly unstable wavenumber is close to kf.

e. The small wavenumber structure of the growth rate

The structure of the dispersion relation at small m provides insight into the nature of zonostrophic instability. Looking at Fig. 4, we anticipate that

 
formula

where η2 > 0 might explain the increase in s that results in the instability with s > 0. This would be a “negative-viscosity instability,” which is the interpretation offered by Farrell and Ioannou (2007), Bakas and Ioannou (2011, manuscript submitted to J. Atmos. Sci., hereafter BakIo), and Bakas and Ioannou (2011).

However there is a small surprise: from (E17) we find that the expansion of the dispersion relation (51) around m = 0 is

 
formula

That is, the term η2 in (58), corresponding to viscosity, is identically zero. Instead, the instability is associated with a destabilizing hyperviscous term, namely the Reynolds stresses are related to the zonal mean flow by

 
formula

leading to the small-m growth rate in (59). We analyze this curious situation further in section 5 and show that η2 = 0 follows from the assumed isotropy of the forcing [i.e., η2 = 0 is not a special property of the particular model in (48)]. The conclusion is that zonostrophic instability requires antifrictional momentum fluxes, and in the small-m limit this antifriction is hyperviscous.

In recent work BakIo and Bakas and Ioannou (2011) reach a different conclusion, namely that the antifrictional effect resulting in nonzero Reynolds stress is equivalent to nonzero and positive η2, and that the hyperviscous coefficient η4 is negative and therefore stabilizing. We believe that these differences may result from a different choice of forcing Ξ. Bakas and Ioannou (2011) and BakIo use an anisotropic forcing, while our conclusion above is specifically for isotropic forcing. The importance of isotropy to our conclusion is underscored in the section 5.

5. Isotropy and zero eddy viscosity

In the discussion surrounding (59) we observed that the term in the zonostrophic dispersion relation corresponding to the eddy viscosity is zero. This result emerges from the analysis of a complicated dispersion relation and surely deserves a more fundamental explanation, or at least another explanation. Thus in this section we more directly obtain the eddy viscosity of an isotropically forced QL β-plane shear flow and show that the result is identically zero.

The eddy viscosity is obtained by calculating the Reynolds stresses in a situation where there is good scale separation between a shear flow and eddies. The best possible scale separation is achieved by considering a Couette flow, Un = γyn, and in this case the CE2 correlation equation (31) collapses to

 
formula

For the moment we assume general forcing [i.e., there is no restriction to isotropic forcing (yet)].

The eddy viscosity νe is defined by the relation

 
formula

The goal is to solve (61) and obtain the Reynolds stress by evaluating Ψxy at zero separation [e.g., as in (33)]. The eddy viscosity then follows from definition (62).

We expect that νe defined above is equal to the coefficient η2 in (58). In the m → 0 limit, the modal solution in (44) varies on the length scale m−1, which is much greater than the length scale of the forcing, namely . Thus, on the scale of the forcing, the growing zonal disturbance resembles the Couette flow1 (except at the “shearless” points, where Uy = 0). By calculating the Reynolds stress in this situation one can anticipate the low-wavenumber structure of the dispersion relation. This reasoning is identical to methods in kinetic theory by which the molecular shear viscosity is calculated.

a. A solution of the correlation equation

We can simplify (61) with (i.e., by looking for a solution independent of ):

 
formula

This exact reduction is surprising because the β effect is removed from the problem. Equation (63) can be solved straightforwardly as an ordinary differential equation in x. However, to make contact with a large literature on sheared disturbances, it is instructive to consider the initial-value problem

 
formula

with the initial condition

 
formula

The solution of the steady problem (63) is then obtained as

 
formula

Thus, solving the initial-value problem for F, the vorticity correlation function is written as the time integral of a sheared disturbance:

 
formula

b. The Reynolds stresses

To obtain the correlation function Ψ from we must solve the two-dimensional biharmonic equation , which is accomplished with the Green’s function defined by ∇4G = δ(x)δ(y), or , or

 
formula

With G(x, y) in hand, we have

 
formula

The Reynolds stress now follows by evaluating Ψxy at zero separation, or

 
formula

This is a very convenient and general expression for the Reynolds stresses in terms of the vorticity correlation function .

Substituting (67) into (70) results in a triple integral. To disentangle this, exchange the order so that the t integral is last, and in the inner x and y integrals “unshear” the correlation function with the coordinate change:

 
formula

After these maneuvers the Reynolds stress is

 
formula

where

 
formula

Now we restrict attention to isotropic forcing; that is,

 
formula

Then in polar coordinates, the double integral in (73) factors:

 
formula

The final line follows from constraint (42) and implies that . That is, the eddy viscosity is zero.

We remark that the constraints in (42) and (43) are required so that correlation function Ψ on the left of (69) decays as r → ∞, despite the r → ∞ divergence of the Green’s function G(r) in (68). In the convolution integral on the right-hand side of (69), the large r divergence of G is shielded by zero integrals of the vorticity correlation function , which follows from the integral constraints on Ξ in (42) and (43).

There are two important caveats associated with the conclusion that νe = 0: the stochastic forcing is isotropic and dissipation is provided only by Ekman drag. Relaxing either or both of these assumptions might result in nonzero νe.

c. The kinetic energy density

The energy power integral for the β-plane Couette flow problem considered here is obtained by first rewriting (63) as

 
formula

Canceling a Laplacian, and evaluating the result at zero separation, one obtains2

 
formula

The left-hand side is the transfer of energy between the eddies and the Couette flow, which is zero because Therefore the statistically energy balance is between dissipation due to drag and the rate of working of the random force that drives the eddies. Remarkably, because the Reynolds stresses are zero, the eddy kinetic energy of the statistically steady flow is ɛ/(2μ), independent of both β and γ.

d. Discussion

To a certain extent the result νe = 0 is anticipated in the literature on sheared disturbances. Shepherd (1985) showed that an isotropic initial distribution of Rossby waves maintains a constant energy density, despite shearing by a Couette flow; see also Farrell and Ioannou (1993a) and Holloway (2010). The solution in (67), with the isotropic initial condition in (65), is essentially a time integral of Shepherd’s solution of the sheared-disturbance problem with an isotropic initial condition.

Via direct numerical simulation (but with β = 0), Cummins and Holloway (2010) have recently shown that the eddy–eddy nonlinearity is essential in producing nonzero Reynolds stresses from Couette-sheared eddies. Cummins and Holloway (2010) identify the essential role of EENL as restoration of isotropy at high wavenumbers. Moreover, as a result of nonlinearly restored isotropy, νe is robustly positive and thus cannot serve as an explanation of zonostrophic instability. Whatever the sign of νe, an unfortunate consequence of (9) is that restoration of isotropy at small scales is absent in QL dynamics and not represented in the ensemble-averaged dynamics CE2.

6. Zonation in QL and NL solutions

We now turn to numerical solutions for a comparison of the full nonlinear system, the quasilinear system, and the predictions of CE2. In these calculations the resolution is 512 × 512, and we use the ETDRK4 time-stepping scheme (Cox and Matthews 2005). In addition to the control parameters β* and μ* defined in (49), there is a third control parameter, which is the size of the domain relative to the forced wavenumber kf: in our computations the domain is a doubly periodic square 2πL × 2πL, with kfL = 32. Thus there is scale separation between the forcing and the domain.

We have obtained about 150 QL and NL numerical solutions, with the planetary vorticity gradient in the range

 
formula

and the drag parameter in the range

 
formula

In this section we use these solutions to compare QL and NL solutions and assess the validity of the CE2 linear stability analysis.

a. The onset of zonation in NL and QL solutions

Shown in the bottom panel of Fig. 3 is the evolution of the fraction of kinetic energy in the zonal mean flow

 
formula

where denotes the area integral over the entire domain. The index zmf(t) is a gross measure of the strength of the zonal mean flow. The time average, denoted by 〈zmf〉, is computed by averaging over an interval t1 < t < t1 + 10/μ, where typically 2μt1 > 40. This long spinup ensures that statistical equilibrium has been achieved and is consistent with the equilibration time suggested by Galperin et al. (2006).

The index 〈zmf〉 is used to classify the flow. Figure 6 summarizes a suite of QL and NL calculations in which the drag parameter is varied at fixed β*. The onset of zonation is indicated by the increase in 〈zmf〉. The dotted lines marked correspond to the critical curve in the upper panel of Fig. 5; these analytic predictions compare well with the increase in 〈zmf〉 in the QL numerical solutions. The dotted lines marked in Fig. 6 are eyeball estimates of the onset of NL zonation.

Fig. 6.

The time-averaged zonal mean energy fraction 〈zmf〉 as a function of μ*, with β* fixed as indicated in the bottom-right corner of each panel. QL simulations are indicated by a degree sign and NL solutions by an asterisk.

Fig. 6.

The time-averaged zonal mean energy fraction 〈zmf〉 as a function of μ*, with β* fixed as indicated in the bottom-right corner of each panel. QL simulations are indicated by a degree sign and NL solutions by an asterisk.

The onset of zonostrophic instability requires significantly smaller values of μ* in the NL case than in the QL case: in Fig. 6 the ratio is as large as 5. Thus the QL system is much more unstable than the NL system. Regarding this quantitative difference between NL and QL, we recall that QL (and CE2) is an approximation based on dropping the eddy–eddy nonlinearity. This approximation is most defensible when the mean flow is very strong (i.e., in cases where the zonal mean flow contains almost all of the energy). Therefore, CE2 is not likely to be quantitatively accurate near the linear stability boundary, where the zonal mean flow is weak or nonexistent. The comparison in Fig. 6 is thus a worst-case test of CE2. How, or if, CE2 might be improved to account for the missing eddy–eddy nonlinearity in this weak mean-flow regime is an open issue.

b. Zonostrophically stable NL solutions

Figure 7 shows two NL solutions that are zonostrophically stable; that is, these solutions have

 
formula

and 〈zmf〉 ≈ 0. In the left panel of Fig. 7 the drag is so heavy that the approximate dominant balance in (1) is μζξ and the vorticity field closely resembles a snapshot of the forcing ξ.

Fig. 7.

Snapshots of the vorticity ζ(x, y, t) with overlaid zonally averaged vorticity −Uy(y, t) (solid white curve) with (a) μ* = 0.309 and (b) μ* = 0.0545. Both snapshots are at nondimensional time 2μt = 25, after spinup from rest, and β* = 1.

Fig. 7.

Snapshots of the vorticity ζ(x, y, t) with overlaid zonally averaged vorticity −Uy(y, t) (solid white curve) with (a) μ* = 0.309 and (b) μ* = 0.0545. Both snapshots are at nondimensional time 2μt = 25, after spinup from rest, and β* = 1.

Figure 8 compares energy spectra of statistically steady QL and NL solutions. With strong drag (i.e., μ* = 0.309) only the directly forced wavenumbers are significantly excited. As μ* is reduced there is transfer of energy to small wavenumbers. In the NL case the transfer of energy to wavenumbers smaller than kf is the due to the inverse energy cascade. In the QL case the excitation of small wavenumbers is due only to shearing by the zonal mean flow. Comparing QL and NL solutions at the same value of μ*, one sees from Figs. 8b and 8d that there is significantly more low-wavenumber eddy energy in the NL cases. Yet the zonal mean energy is always stronger in the QL case. There is no clear association between the inverse energy cascade and zonation.

Fig. 8.

(a),(c) The zonal spectrum EZ(ky/kf) for QL and NL solutions with β* = 1. (b),(d) The residual spectrum ER(k/kf), defined as the angularly averaged spectrum after removal of the “zonal modes” with kx = 0. The largest peak in EZ(ky/kf) defines the wavenumber mZ, even if there are no quasi-steady zonal jets [e.g., as in the NL simulation with μ* = 0.0545 in (a)].

Fig. 8.

(a),(c) The zonal spectrum EZ(ky/kf) for QL and NL solutions with β* = 1. (b),(d) The residual spectrum ER(k/kf), defined as the angularly averaged spectrum after removal of the “zonal modes” with kx = 0. The largest peak in EZ(ky/kf) defines the wavenumber mZ, even if there are no quasi-steady zonal jets [e.g., as in the NL simulation with μ* = 0.0545 in (a)].

The NL solution shown in right panel of Fig. 7 with μ* = 0.0545 has an eddy energy spectrum in Fig. 8b exhibiting the beginning of range. However this solution has 〈zmf〉 ≈ 0 and thus serves as example of an isotropic, spatially homogeneous, weakly turbulent, β-plane flow, without jets. To activate zonostrophic instability the drag must be reduced (e.g., to μ* = 0.0182 in Figs. 1 and 8).

c. The jet scale

If zonation occurs, as evinced by significantly nonzero values of 〈zmf〉, then by counting the number of distinct jets one can reliably estimate3 a jet wavenumber mJ. For example, in Fig. 1 there are seven jets and therefore .

However we noticed that there are cases without jets in which the zonal energy spectrum EZ(ky/kf) has a strong peak: an example is the μ* = 0.0545 solution in Fig. 7b: the corresponding zonal energy spectrum in Fig. 8a has a distinct peak even though there are no zonal jets. In cases like this, we report a wavenumber mZ that is the peak of the zonal spectrum EZ(ky/kf). In cases where there are strong jets we invariably find that mZmJ. It is interesting to compare mJ and mZ with a Rhines wavenumber defined as

 
formula

where the root-mean-square velocity is

 
formula

We investigated other choices for the velocity in the Rhines wavenumber; for example, Rhines (1975) advocated the RMS of υ′. We found, however, that VRMS gave the best estimate of the NL jet spacing at small values of μ*. An advantage of VRMS is that the energy power integral4 can be used to express VRMS in terms of external parameters as

 
formula

The relation above applies with an error (due to hyperviscous dissipation) of less than 5% in our simulations. Substituting (82) into (80) one has

 
formula

Figure 9 compares the zonal wavenumber obtained from QL and NL solutions with the Rhines wavenumber on the right-hand side of (83), and with the most unstable wavenumber obtained from the linear stability analysis of section 4. In Fig. 9 we show only the β* = 1 and β* = 0.5 solutions: solutions with other values of β* exhibit a broadly similar dependence of mZ on μ*.

Fig. 9.

A summary of zonal wavenumbers (jet scales) for solutions with (a) β* = 0.5 and (b) β* = 1. The dot-dashed curve is the Rhines wavenumber defined in (83). The solid curve labeled QLS is most unstable wavenumber calculated from the dispersion relation (51). The NL solutions are indicated by asterisks and the QL solutions by degree signs.

Fig. 9.

A summary of zonal wavenumbers (jet scales) for solutions with (a) β* = 0.5 and (b) β* = 1. The dot-dashed curve is the Rhines wavenumber defined in (83). The solid curve labeled QLS is most unstable wavenumber calculated from the dispersion relation (51). The NL solutions are indicated by asterisks and the QL solutions by degree signs.

At large values μ* only the directly forced modes are excited, and consequently mZkf in both the QL and NL cases. At the critical value in Fig. 9, the QL solutions undergo zonostrophic instability, and close to this transition (e.g., at μ* = 0.2 and 0.165 in Fig. 9a) the QL mZ agrees with the analytic result from Fig. 5. In this regime the NL solutions start to develop an inverse cascade (but without exciting zonal jets) and the NL mZ begins to decrease.

There is an interesting transition at in Fig. 9. At this point the QL and NL zonal wavenumbers are equal, and as μ* is reduced the QL and NL wavenumbers are locked together. At in Fig. 9 the NL solutions finally become zonostrophically unstable, resulting in NL jets and significantly nonzero values of NL 〈zmf〉. At the smallest value of μ* in Fig. 9, which corresponds to the runs in Figs. 1 and 2, the QL and NL wavenumbers are almost equal and are estimated roughly by mRh.

In Fig. 9 the analytic result QLS agrees with the observed QL jet scale only when μ* is not too far from the linear stability boundary In the strongly unstable regime, with μ* significantly less than the observed wavenumber is much smaller than the most unstable wavenumber predicted by linear theory. This increase in the QL jet scale is a result of merging jets that initially appear with a spacing, which is well predicted by the linear theory. This phenomenology begins at about and is illustrated in Fig. 10.

Fig. 10.

Hovmöller diagrams for the (a) NL and (b) QL runs with β* = 1.0 and μ* = 0.0545. The NL run corresponds to the vorticity snapshot shown in Fig. 7b and shows zonal “streaks.” In (b) the QL jets initially appear at a wavenumber predicted by linearization of CE2. Then successive mergers result in an increase in jet spacing.

Fig. 10.

Hovmöller diagrams for the (a) NL and (b) QL runs with β* = 1.0 and μ* = 0.0545. The NL run corresponds to the vorticity snapshot shown in Fig. 7b and shows zonal “streaks.” In (b) the QL jets initially appear at a wavenumber predicted by linearization of CE2. Then successive mergers result in an increase in jet spacing.

Figure 10a shows the Hovmöller diagram of the jetless NL solution from Fig. 7b. There is no zonation and U(y, t) shows “streaks” rather than jets. These streaks are not strong relative to the turbulent eddy field (i.e., 〈zmf〉 ≈ 0). The corresponding zonal energy spectrum in Fig. 8a exhibits a strong peak, which is a signature of these transient zonal steaks.

Figure 10b shows the QL case in which jets initially appear with a relatively small meridional spacing predicted by linear theory, followed by a sequence of mergers so that the mature flow has mZ much less than the most linearly unstable wavenumber. The QL jet-merging phenomenology, which is effectively a one-dimensional inverse cascade, is very similar to the “Cahn–Hilliard” solutions obtained by Manfroi and Young (1999) from a model of deterministically forced zonation.

d. The small drag regime

The flows in Figs. 1 and 2 have relatively light damping and both flows have organized jets containing a substantial fraction of the total kinetic energy. Figure 11 shows the time-averaged zonal mean flow 〈U〉 and the corresponding PV gradient β* − 〈Uyy〉. In Fig. 2 the QL jets are almost symmetrical in the zonal direction, in contrast to the NL jets.5 But the QL jets are not perfectly symmetric: the PV gradient in Fig. 11b reveals the QL east–west asymmetry. The NL PV gradient is positive for all y and thus the NL jets are stable according to the Rayleigh–Kuo criterion. The QL PV gradient in Fig. 11b reverses sign on the flanks of the eastward jet, and also at the centers of the westward jets. Nonetheless the QL zonal mean flow shows no indication of barotropic instability [i.e., the deep spikes with β* − 〈Uyy〉 < 0 are permanent features of the QL zonal mean flow even after time averaging].

Fig. 11.

Comparison of zonal mean velocity profiles of the β* = 1 NL and QL runs in Figs. 1 and 2.

Fig. 11.

Comparison of zonal mean velocity profiles of the β* = 1 NL and QL runs in Figs. 1 and 2.

Via integration of their SSST system, Farrell and Ioannou (2007) report equilibrated zonal mean flows with much stronger east–west asymmetry than the QL flow in Fig. 11 (e.g., see their Figs. 8 and 9). There are at least6 two nondimensional parameters β* and μ*, and the jet profile depends on both of these. We will not attempt to characterize this variation systematically. However, to make some contact with the strong-forcing limit considered by Farrell and Ioannou (2007), we consider the QL solution in Figs. 10b and 12a and increase the energy injection rate ɛ by a factor of 1000, while holding β, μ and kf approximately fixed. Then from (49), β* and μ* are each reduced by a factor of 10. The time-averaged zonal mean profile of this strongly forced solution is shown in Fig. 12b and exhibits the parabolic velocity profile seen in the NL run in Fig. 1: there are fast eastward jets with sharp gradients and broad westward jets with smaller PV gradients. Also, the time-averaged QL jet profile in Fig. 12b is more asymmetric than the weaker forced QL jet shown in Fig. 12a, which has a forcing that is a factor of 10 smaller. To quantify the jet asymmetry, we use the ratio

 
formula

where Umax and Umin are the maximum and minimum values attained in the zonal mean velocity profile. By increasing the forcing strength by a factor of 1000, the jet asymmetry increases from α = 1.25 in Fig. 12a to α = 1.56 for the profile in Fig. 12b. This is smaller than the “ideal,” marginally stable (i.e., βUyy = 0 everywhere except at the eastward jet where the PV jumps) profile considered in Danilov and Gurarie (2004), which has α = 2.

Fig. 12.

Comparison of time-averaged zonal mean velocity profiles (thin lines) of QL runs. (a) The solution from Fig. 10b with β* = 1.0 and μ* = 0.054; (b) the strongly forced solution with β* = 0.1 and μ* = 0.005. Also plotted are the corresponding PV gradients (thick curves).

Fig. 12.

Comparison of time-averaged zonal mean velocity profiles (thin lines) of QL runs. (a) The solution from Fig. 10b with β* = 1.0 and μ* = 0.054; (b) the strongly forced solution with β* = 0.1 and μ* = 0.005. Also plotted are the corresponding PV gradients (thick curves).

Thus, although a detailed study of QL jet asymmetry is not a focus of the present work, our QL numerical solutions are generally consistent with the equilibrated SSST jets presented in Farrell and Ioannou (2007).

e. Discussion of the eddy–eddy nonlinearity

An important effect of eddy–eddy nonlinearity is the stirring of PV, producing an exponential-in-time reduction in the length scale of vorticity fluctuations. Eddy-driven stirring is removed from the QL system by (9): shearing by U(y, t) is the only scale-reduction mechanism acting on the QL eddy vorticity. The small-scale structure evident in the QL PV gradient in the right panel of Fig. 11 may reflect the relative inefficiency of shearing by U(y, t) at removing vorticity fluctuations.

Further differences in the jet structure evident in Fig. 11 can be explained by meandering of the NL jets, so that the zonal average reduces the sharpness of the NL PV gradient. The spectral signature of the NL jet meanders is a high energy mode at in the two-dimensional NL spectrum; this same mode is only weakly excited in the QL spectrum. The excitation of almost-zonal modes, with a small but nonzero of value of kx, is a well-known aspect of zonation. These are called a “satellite modes” by Danilov and Gurarie (2004), and they correspond to a domain-scale meander of the NL jets, which is not present in the QL case.

7. Discussion and conclusions

A contribution of this work is the analytic development of the linearized theory of zonostrophic instability within the context of the second-order cumulant expansion (CE2) of Marston et al. (2008) and the stochastic structural stability theory (SSST) of Farrell and Ioannou (2003, 2007). These statistical formulations are equivalent to the correlation dynamics derived in section 3, and that physical-space formulation, in terms of partial differential equations for the correlation functions Ψ and , provides some insight into the mathematical structure of CE2/SSST.

In the top panel of Fig. 5 we display the curve of neutral zonostrophic stability in the (β*, μ*)-parameter plane obtained by solution of linearized CE2 dynamics. We have shown that with isotropic forcing zonostrophic instability is not a negative-viscosity instability: the hallmark of a negative-viscosity instability is that at the stability boundary the most unstable wavenumber is zero. The deterministic model of anisotropically forced β-plane zonation analyzed by Manfroi and Young (1999) provides a bona fide example of the negative-viscosity case. Instead, for the isotropically and stochastically forced model analyzed here, the onset of zonostrophic instability is at the nonzero meridional wavenumber shown in the bottom panel of Fig. 5; only at large β* does this wavenumber approach zero. Moreover, in section 5 we showed that with isotropic forcing the CE2 eddy viscosity νe is identically zero.

Comparison of QL and NL numerical solutions indicates that the CE2 linear stability boundary does not provide an accurate estimate of the onset of zonostrophic instability for NL flows. This quantitative failure of CE2 is not surprising: neglect of the eddy–eddy nonlinearity is most plausible in cases where most of the energy is in the zonal mean flow: close to the stability boundary the zonal mean flow is only incipient. An outstanding open problem is improving CE2 to account for the missing physics in the eddy–eddy nonlinearity. Another important problem is obtaining analytic insight into the solution of the CE2 system in the regime where CE2 is likely to be valid, namely in the strongly unstable regime where the drag μ* is much less than the critical drag μc and the fraction of energy in the zonal mean flow is substantial.

Acknowledgments

This work was supported by the National Science Foundation under OCE1057838. We thank Nikos Bakas, Oliver Bühler, George Carnevale, Brian Farrell, Petros Ioannou, Brad Marston, Rick Salmon, and Steve Tobias for discussion of these results.

APPENDIX A

Implementation of the Random Forcing ξ(x, t)

For numerical purposes we model the δ-correlated forcing ξ(x, y, t) in (3) using a discrete approximation. The goal is to construct a statistically isotropic and narrowband forcing localized close to a radial wavenumber kf. Thus, the forcing is confined to an annulus in wavenumber space, where wavenumbers k in satisfy the inequality

 
formula

We take δk = kf/8, so that is tightly around the average radius kf. We use a fourth-order Runge–Kutta scheme, with time step δt. Implementing the Runge–Kutta scheme requires the value of the forcing not just at points in time separated by the time step δt but also at the midpoints. Some care must be exercised here, though, since the Runge–Kutta scheme assumes a certain degree of smoothness of the solution. To ensure this, we use a forcing that during the nth time step, when (n − 1)δt < t < nδt, has the form

 
formula

where varies linearly from zero to one during the nth time step. The coefficient ξi(n) above is

 
formula

where is the number of wave vectors in , and is the rate of energy injection. The dependence ensures that the forcing is δ correlated in the limit δt → 0. The phase ϕi(n) is a random variable, chosen from a uniform distribution in [0, 2π]; the phase is set independently for each wave vector ki and resets at the start of each time step.

The narrowband forcing is described by the correlation function

 
formula

which has the physical-space form

 
formula

A comparison of Ξb(r) and its idealized form, Ξ(r) ∝ J0(kfr) is shown in Fig. A1.

Fig. A1.

Comparison of the ring forcing δk → 0 (solid line) and the narrow-band forcing with δk = kf/4 (dashed line).

Fig. A1.

Comparison of the ring forcing δk → 0 (solid line) and the narrow-band forcing with δk = kf/4 (dashed line).

APPENDIX B

Rapid Temporal Decorrelation: Derivation of (25)

The crucial assumption leading to (25) is that the temporal decorrelation of the forcing is rapid, as indicated by the δ(t1t2) on the right-hand side of (14). Operationally, this means that we might integrate (11) during the first time step from t = 0 to t = δt as

 
formula

where AOT(x, 0) indicates “all other terms” in (11), evaluated at t = 0. Also in (B1), is a spatial random field with correlation function

 
formula

The forcing “renovates” during each time step; that is, in the nth time step one creates a new independent realization of , but always with the same spatial correlation function Ξ. According to this recipe the magnitude of is independent of δt as δt → 0, and therefore as δt → 0. As demanded by this argument, notice that ξi(n) in (A3) is proportional to .

Before ensemble averaging, we can multiply (B1) evaluated at x1 with (B1) evaluated at x2 to obtain

 
formula

where the subscript n indicates evaluation at xn; for example,

 
formula

Upon ensemble averaging, the terms of order on the right-hand side of (B3) vanish because is independent of ζ′(x, 0). Thus,

 
formula

As δt → 0 the left-hand side is the time derivative of the vorticity correlation function. The terms in (B3), which prohibit a differentiable limit, are nulled in (B5) by the ensemble average. Thus, taking the limit δt → 0 in (B5), we obtain the deterministic differential equation (23) for the evolution of the correlation function .

APPENDIX C

Derivation of the Dispersion Relation (45)

The linearized equations resulting from substitution of (44) into (34) and (31) are

 
formula
 
formula
 
formula

where in (C1)

 
formula

A key intermediate step on the path to (C1) is noting and , so that with U(y, t) in (44) one has

 
formula

Applying the Fourier transform in (38) to (C1) and (C3), one has

 
formula
 
formula
 
formula

In (C6)(C8) we use the notation

 
formula
 
formula
 
formula
 
formula

Eliminating between (C6) and (C7), we obtain the dispersion relation

 
formula

where the functions Λ+(s′, m) and Λ(s′, m) are defined by the integral

 
formula

Changing variables with p → −p and q → −q, and using the exchange symmetry in (22), one finds that

 
formula

so that the right-hand side of (C13) is equal to 2mΛ. Then with the change of variables q′ = qm/2 in the Λ integral, and using (41), one can write the dispersion relation (C13) as

 
formula

where .

If the forcing is isotropic then , and the integral on the right-hand side of (C16) can be simplified using polar coordinates (p, q) = h(cosθ, sinθ):

 
formula

where S is the function

 
formula

One can show that S(χ, n) = −S(−χ, n) = −S(χ, −n), and therefore S(0, n) = S(χ, 0) = 0. These symmetries are important for further work, and they are not manifest from the definition of S in (C18). Thus we seek an alternative form with more obvious properties. The change of variables θθ + π results in

 
formula

The average of (C18) and (C19) is then

 
formula

The function Q(χ, n) is manifestly an even function of n, and θ → −θ shows that Q is also an even function of χ.

Substituting (C20) into (C17) gives the dispersion relation in the form

 
formula

If ν = 0 then s′ = s + 2μ, and we obtain the dispersion relation in (45).

APPENDIX D

The Function Q(χ, n)

In this appendix we summarize some properties of the function Q(χ, n) defined in (C20).

We first note that

 
formula

If 0 ≤ n2 ≤ 1 then

 
formula

where Pr refers to the Cauchy principal value.

The case β* ≪ 1 requires the approximation of Q(χ, n) in the limit χ → ∞. One can expand the integrand in inverse powers of χ and integrate term by term. The first two nonzero terms are

 
formula

The case β* ≫ 1 requires the approximation of Q(χ, n) in the limit χ → 0. A somewhat laborious “range-splitting” calculation shows that

 
formula

Finally, if χ → ∞, with n = 0 then from either (D1) or (D3) we obtain

 
formula

APPENDIX E

The Neutral Curve

The neutral curve in the (β*, μ*)-parameter plane is defined by the conditions in (53). For the dispersion relation in (51), these take the form

 
formula
 
formula

with χ0 = 2μ/ above. (For brevity, in this appendix, we drop the asterisks indicating nondimensional variables.) An examination of the numerical values of χ0 on the neutral curve motivates the possibility that χ0 → ∞ as β → 0 and χ0 → 0 as β → ∞. We now use this numerical observation to derive analytical approximations for the neutral curve in the complementary limiting cases β → 0 and β → ∞. To this end, we use the approximations to Q(χ, m) summarized in  appendix D.

a. Approximation of the marginal curve, β ≫ 1 and χ0 ≪ 1

First in the case of β → ∞, we have from (D4)

 
formula

Clearly, neglecting the O(0) term must be justified post facto once a consistent dominant balance is found. Substituting (E3) in the neutral curve equations, (E1) and (E2) and keeping in mind that that β−1 ≪ 1, we get

 
formula
 
formula

The only consistent balance in the m equation corresponds to m3 ~ 3β−3 and, consequently,

 
formula

Since m ~ O(β−1) the only consistent balance in (E1) is μ ~ 2β−2. A higher-order estimate for μ can be derived by substituting for m to get

 
formula

These approximate expressions are superimposed on the numerically obtained neutral curve in the Fig. 5 with agreement once β > 2.

b. Approximation of the marginal curve, β ≪ 1 and χ0 ≫ 1

A similar analysis is used when β → 0 for which we can write, from (D3),

 
formula

and therefore

 
formula

Substituting into (E1) and (E2) one has

 
formula
 
formula

Above, h.o.t. refers to the higher-order terms that have been neglected from the above equation and can be justified to be small post facto. In (E10), assuming a dominant balance between the term on the left and the first term on the right-hand side, we have

 
formula

Similarly from (E12)

 
formula

which on substituting for μ from (E12), gives us m ~ 1 and

 
formula

Proceeding to the next order, we obtain after some algebra

 
formula

As in the case before, we plot this estimate for m against the numerical estimate in Fig. 5 and the agreement is excellent. In fact, the approximations to mc(β) practically span the entire range of the neutral curve.

c. The small-m expansion

When m → 0 with all other parameters fixed, we observe in Fig. 4 that s → −μ and therefore in (51) the first argument of Q is

 
formula

To compute the small-m expansion of the dispersion relation we write s = −μ + s1 and use (D5) to approximate Q. Thus the dispersion relation (51) becomes

 
formula

This produces the small-m version of the dispersion relation in (59).

REFERENCES

REFERENCES
Bakas
,
N. A.
, and
P.
Ioannou
,
2011
:
Structural stability theory of two-dimensional fluid flow under stochastic forcing
.
J. Fluid Mech.
,
682
,
332
361
.
Carnevale
,
G. F.
, and
P. C.
Martin
,
1982
:
Field theoretical techniques in statistical fluid dynamics: With applications to nonlinear wave dynamics
.
Geophys. Astrophys. Fluid Dyn.
,
20
,
131
164
.
Cox
,
S. M.
, and
P. C.
Matthews
,
2005
:
Exponential time differencing for stiff systems
.
J. Comput. Phys.
,
176
,
430
455
.
Cummins
,
P. F.
, and
G.
Holloway
,
2010
:
Reynolds stress and eddy viscosity in direct numerical simulation of sheared-two-dimensional turbulence
.
J. Fluid Mech.
,
657
,
394
412
.
Danilov
,
S.
, and
V. M.
Gryanik
,
2004
:
Barotropic beta-plane turbulence in a regime with strong zonal jets revisited
.
J. Atmos. Sci.
,
61
,
2283
2295
.
Danilov
,
S.
, and
D.
Gurarie
,
2004
:
Scaling, spectra and zonal jets in beta-plane turbulence
.
Phys. Fluids
,
16
,
2592
2603
.
DelSole
,
T.
,
2001
:
A theory for the forcing and dissipation in stochastic turbulence models
.
J. Atmos. Sci.
,
58
,
3762
3775
.
Diamond
,
P. H.
,
S.-I.
Itoh
,
K.
Itoh
, and
T. S.
Hahm
,
2005
:
Zonal flows in plasmas—A review
.
Plasma Phys. Contr. Fusion
,
47
,
R35
,
doi:10.1088/0741-3335/47/5/R01
.
Dritschel
,
D. G.
, and
M. E.
McIntyre
,
2010
:
Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers
.
J. Atmos. Sci.
,
65
,
855
874
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
1993a
:
Stochastic forcing of perturbation variance in unbounded shear and deformation flows
.
J. Atmos. Sci.
,
50
,
200
211
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
1993b
:
Stochastic forcing of the linearized Navier–Stokes equations
.
Phys. Fluids A
,
5
,
2600
2609
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
2003
:
Structural stability of turbulent jets
.
J. Atmos. Sci.
,
50
,
2101
2118
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
2007
:
Structure and spacing of jets in barotropic turbulence
.
J. Atmos. Sci.
,
64
,
3652
3665
.
Galperin
,
B.
,
S.
Sukoriansky
,
N.
Dikovskaya
,
P.
Read
,
Y.
Yamazaki
, and
R.
Wordsworth
,
2006
:
Anisotropic turbulence and zonal jets in rotating flows with a β-effect
.
Nonlinear Processes Geophys.
,
13
,
83
98
.
Holloway
,
G.
,
2010
:
Eddy stress and shear in 2D flows
.
J. Turbul.
,
11
,
1
14
.
Huang
,
H. P.
, and
W. A.
Robinson
,
1998
:
Two-dimensional turbulence and persistent zonal jets in a global barotropic model
.
Nonlinear Processes Geophys.
,
55
,
611
632
.
Lilly
,
D. K.
,
1969
:
Numerical simulation of two-dimensional turbulence
.
Phys. Fluids
,
12
,
II-240
,
doi:10.1063/1.1692444
.
Maltrud
,
M. E.
, and
G. K.
Vallis
,
1991
:
Energy spectra and coherent structures in forced two-dimensional and beta-plane turbulence
.
J. Fluid Mech.
,
228
,
321
342
.
Manfroi
,
A. J.
, and
W.
Young
,
1999
:
Slow evolution of zonal jets on the beta plane
.
J. Atmos. Sci.
,
56
,
784
800
.
Marston
,
J. B.
,
E.
Conover
, and
T.
Schneider
,
2008
:
Statistics of an unstable barotropic jet from a cumulant expansion
.
J. Atmos. Sci.
,
65
,
1955
1966
.
Nozawa
,
T.
, and
S.
Yoden
,
1997
:
Formation of zonal-band structure in forced two-dimensional turbulence on a rotating sphere
.
Phys. Fluids
,
9
,
2081
2093
.
O’Gorman
,
P. A.
, and
T.
Schneider
,
2007
:
Recovery of atmospheric flow statistics in a general circulation model without nonlinear eddy–eddy interactions
.
Geophys. Res. Lett.
,
34
,
L22801
,
doi:10.1029/2007GL031779
.
Panetta
,
R. L.
,
1993
:
Zonal jets in wide baroclinically unstable regions: persistence and scale selection
.
J. Atmos. Sci.
,
50
,
2073
2106
.
Rhines
,
P. B.
,
1975
:
Waves and turbulence on a beta-plane
.
J. Fluid Mech.
,
69
,
417
443
.
Rhines
,
P. B.
, and
W. R.
Young
,
1982
:
Homogenization of potential vorticity in planetary gyres
.
J. Fluid Mech.
,
122
,
347
367
.
Salmon
,
R.
,
1998
:
Lectures on Geophysical Fluid Dynamics. Oxford University Press, 379 pp
.
Scott
,
R. K.
, and
L. M.
Polvani
,
2007
:
Forced-dissipative shallow-water turbulence on the sphere and the atmospheric circulation of the giant planets
.
J. Atmos. Sci.
,
64
,
3158
3176
.
Shepherd
,
T. G.
,
1985
:
Time development of small disturbances to plane Couette flow
.
J. Atmos. Sci.
,
42
,
1868
1871
.
Showman
,
A. P.
,
2007
:
Numerical simulations of forced shallow-water turbulence: Effects of moist convection on large-scale circulation of Jupiter and Saturn
.
J. Atmos. Sci.
,
64
,
3132
3157
.
Smith
,
K. S.
,
2004
:
A local model for planetary atmospheres forced by small-scale convection
.
J. Atmos. Sci.
,
61
,
1420
1433
.
Tobias
,
S. M.
,
K.
Dagon
, and
J.
Marston
,
2011
:
Astrophysical fluid dynamics via direct statistical simulation
.
Astrophys. J.
,
727
,
127
,
doi:10.1088/0004-637X/727/2/127
.
Vallis
,
G. K.
, and
M. E.
Maltrud
,
1993
:
Generation of mean flow and jets on a beta plane and over topography
.
J. Phys. Oceanogr.
,
23
,
1346
1362
.
Vasavada
,
A. R.
, and
A. P.
Showman
,
2005
:
Jovian atmospheric dynamics: An update after Galileo and Cassini.
Rep. Prog. Phys.
,
68
,
1935
1996
.
Whitehead
,
J. A.
,
1975
:
Mean flow generated by circulation on a beta-plane: An analogy with the moving flame experiment
.
Tellus
,
27
,
358
363
.
Williams
,
G. P.
,
1978
:
Planetary circulations: 1. Barotropic representation of Jovian and terrestrial turbulence
.
J. Atmos. Sci.
,
35
,
1399
1426
.

Footnotes

1

There is also uniform advection by the zonal flow. But that sweeping is eliminated by the difference U1U2 in the correlation equation (31) and is therefore inconsequential to Reynolds stresses.

2

The rate of energy injection is ɛ = ½∇2A|0; see, for example, the model forcing in (47).

3

In certain cases the system may be transitioning between a state with n and n + 1 jets. Following Panetta (1993), we then count jets; no other fractional values are permitted.

4

From (3), the exact energy power integral is 〈ψξ + μ|ψ|2 + (−1)n−1νn|n−1ζ|2〉 = 0, where 〈〉 is both a domain integral and a time average.

5

If β = 0 then the equations of motion are invariant under y → −y and ψ → −ψ. This symmetry, which induces uu, is broken in both QL and NL by nonzero β. This explains the characteristic east–west asymmetry of U(y, t) on the β plane.

6

Farrell and Ioannou (2007) also employ a forcing with a different correlation function than our isotropic choice.