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.
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).
To establish the notation used in this study, we start with the equations of motion for a forced and dissipative barotropic flow u = (u, υ):
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
The vorticity forcing ξ on the right-hand side of (11) is the curl of the solenoidal force in the momentum equation—that is,
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, y → y + 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
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
and the eddy vorticity equation
In (7), the eddy–eddy nonlinearity (EENL) is
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
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.
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.
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:
(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
where is the Rayleigh–Kuo operator:
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
where the overbar above denotes an ensemble average. The dependence of the spatial correlation functions A and Ξ only on the difference
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
where the Laplacian acting on the coordinates of point n is
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
The analog of (16) is
Given the streamfunction correlation Ψ(x, y1, y2, t), one can obtain the velocity correlation tensor as
Because the choice of denoting one point as x1 and the other as x2 is arbitrary, all correlation function have an important “exchange” symmetry
and likewise for A, , Ξ, etc.
b. Correlation functions: Dynamics
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 = x1 − x2. Because of this zonal homogeneity . Thus, for example,
A crucial simplification is that the forcing is rapidly decorrelating in time, as expressed by the δ(t1 − t2) in (14). Considerations summarized in appendix B (amounting to a simple proof of Ito’s formula) show that
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”
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
Collective coordinates are then essential for analytic progress.
In terms of y and , the Laplacians are
where is the “separation” Laplacian. Thus, for instance,
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
at (x, y) = 0, is that the Reynolds stress is
Thus, the mean flow equation (6) can be written as
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
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
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,
to (36). We use the notation
so that after the Fourier transform , and the streamfunction spectrum is related to the forcing spectrum by
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
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 Ξ:
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
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
where is the forcing spectrum in (41), and the function Q(χ, n) is defined by the angular integral
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δk ≪ kf. 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
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 δ(h − kf) 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
The nondimensional wavenumber and growth rate are
The zonostrophic dispersion relation in nondimensional variables is then
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:
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).
The marginal stability condition in (52), which is equivalent to the requirements
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)].
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
In the complementary limit β* ≫ 1, the approximation to the critical curve is
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
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).
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
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
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
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 ):
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
with the initial condition
The solution of the steady problem (63) is then obtained as
Thus, solving the initial-value problem for F, the vorticity correlation function is written as the time integral of a sheared disturbance:
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
With G(x, y) in hand, we have
The Reynolds stress now follows by evaluating Ψxy at zero separation, or
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:
After these maneuvers the Reynolds stress is
Now we restrict attention to isotropic forcing; that is,
Then in polar coordinates, the double integral in (73) factors:
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
Canceling a Laplacian, and evaluating the result at zero separation, one obtains2
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 γ.
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
and the drag parameter in the range
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
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.
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
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.
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 mZ ≈ mJ. It is interesting to compare mJ and mZ with a Rhines wavenumber defined as
where the root-mean-square velocity is
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
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 μ*.
At large values μ* only the directly forced modes are excited, and consequently mZ ≈ kf 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.
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].
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
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.
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.
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.
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
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
where varies linearly from zero to one during the nth time step. The coefficient ξi(n) above is
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
which has the physical-space form
A comparison of Ξb(r) and its idealized form, Ξ(r) ∝ J0(kfr) is shown in Fig. A1.
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 δ(t1 − t2) 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
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 .
where the subscript n indicates evaluation at xn; for example,
Upon ensemble averaging, the terms of order on the right-hand side of (B3) vanish because is independent of ζ′(x, 0). Thus,
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 .
Derivation of the Dispersion Relation (45)
where in (C1)
where the functions Λ+(s′, m) and Λ−(s′, m) are defined by the integral
Changing variables with p → −p and q → −q, and using the exchange symmetry in (22), one finds that
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θ):
where S is the function
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
The function Q(χ, n) is manifestly an even function of n, and θ → −θ shows that Q is also an even function of χ.
If ν = 0 then s′ = s + 2μ, and we obtain the dispersion relation in (45).
The Function Q(χ, n)
In this appendix we summarize some properties of the function Q(χ, n) defined in (C20).
We first note that
If 0 ≤ n2 ≤ 1 then
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
The case β* ≫ 1 requires the approximation of Q(χ, n) in the limit χ → 0. A somewhat laborious “range-splitting” calculation shows that
The Neutral Curve
with χ0 = 2μ/mβ 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)
Clearly, neglecting the O(mχ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
The only consistent balance in the m equation corresponds to m3 ~ 3β−3 and, consequently,
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
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),
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
Similarly from (E12)
which on substituting for μ from (E12), gives us m ~ 1 and
Proceeding to the next order, we obtain after some algebra
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
This produces the small-m version of the dispersion relation in (59).
There is also uniform advection by the zonal flow. But that sweeping is eliminated by the difference U1 − U2 in the correlation equation (31) and is therefore inconsequential to Reynolds stresses.
The rate of energy injection is ɛ = ½∇2A|0; see, for example, the model forcing in (47).
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.
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.
If β = 0 then the equations of motion are invariant under y → −y and ψ → −ψ. This symmetry, which induces u → u, is broken in both QL and NL by nonzero β. This explains the characteristic east–west asymmetry of U(y, t) on the β plane.
Farrell and Ioannou (2007) also employ a forcing with a different correlation function than our isotropic choice.