## Abstract

Zonal jets and nonzonal large-scale flows are often present in forced–dissipative barotropic turbulence on a beta plane. The dynamics underlying the formation of both zonal and nonzonal coherent structures is investigated in this work within the statistical framework of stochastic structural stability theory (S3T). Previous S3T studies have shown that the homogeneous turbulent state undergoes a bifurcation at a critical parameter and becomes inhomogeneous with the emergence of zonal and/or large-scale nonzonal flows and that these statistical predictions of S3T are reflected in direct numerical simulations. In this paper, the dynamics underlying the S3T statistical instability of the homogeneous state as a function of parameters is studied. It is shown that, for weak planetary vorticity gradient *β*, both zonal jets and nonzonal large-scale structures form from upgradient momentum fluxes due to shearing of the eddies by the emerging infinitesimal flow. For large *β*, the dynamics of the S3T instability differs for zonal and nonzonal flows but in both the destabilizing vorticity fluxes decrease with increasing *β*. Shearing of the eddies by the mean flow continues to be the mechanism for the emergence of zonal jets while nonzonal large-scale flows emerge from resonant and near-resonant triad interactions between the large-scale flow and the stochastically forced eddies. The relation between the formation of large-scale structure through modulational instability and the S3T instability of the homogeneous state is also investigated and it is shown that the modulational instability results are subsumed by the S3T results.

## 1. Introduction

Atmospheric turbulence is commonly observed to be organized into slowly varying large-scale structures such as zonal jets and coherent vortices. Prominent examples are the banded jets and the Great Red Spot in the Jovian atmosphere (Ingersoll 1990; Vasavada and Showman 2005). Laboratory experiments as well as direct numerical simulations of turbulent flows have shown that these coherent structures appear and persist for a very long time despite the presence of eddy mixing (Vallis and Maltrud 1993; Weeks et al. 1997; Read et al. 2004; Espa et al. 2010; Di Nitto et al. 2013).

A model that exhibits many aspects of turbulent self-organization into coherent structures yet is simple enough to extensively investigate is a barotropic flow on the surface of a rotating planet or on a beta plane with turbulence sustained by random stirring. Numerical simulations of this model have shown that robust zonal jets coexist with large-scale westward-propagating coherent waves (Sukoriansky et al. 2008; Galperin et al. 2010). These waves were found to either obey a Rossby wave dispersion or form nondispersive packets that are referred to as satellite modes (Danilov and Gurarie 2004) or zonons (Sukoriansky et al. 2008). In addition, the formation of these coherent structures was shown to be a bifurcation phenomenon. As the energy input of the stochastic forcing is increased, the flow bifurcates from a turbulent, spatially homogeneous state to a state in which zonal jets and/or nonzonal coherent structures emerge and are maintained by turbulence (Bakas and Ioannou 2013a; Constantinou et al. 2014). In this work, we address the eddy–mean flow dynamics underlying the emergence of both zonal and nonzonal structures.

Since organization of turbulence into coherent structures involves complex nonlinear interactions among a large number of degrees of freedom, which erratically contribute to the maintenance of the large-scale structure, an attractive approach is to study the statistical state dynamics (SSD) of the turbulent flow rather than single realizations of the turbulent field. Recently, the SSD of barotropic and baroclinic atmospheres has been studied by truncating the infinite hierarchy of cumulant equations to second order. Stochastic structural stability theory (S3T) is such a second-order Gaussian approximation of the full SSD, in which the third cumulant is parameterized as the sum of a known correlation function and a dissipation term (Farrell and Ioannou 2003). This is equivalent to a parameterization of the eddy–eddy nonlinearity as an exogenous random forcing with the required dissipation to remove the energy injected by the forcing. Such a representation is strongly supported by the results of previous studies. Linear inverse modeling studies showed that this parameterization is the best linear representation of the eddy–eddy nonlinear interactions in planetary turbulence (DelSole and Farrell 1996; DelSole 1996; DelSole and Hou 1999; DelSole 2004), while earlier studies have shown that the turbulent transport properties (heat and momentum fluxes) of the midlatitude transient climatology are accurately obtained as the stochastic response of the large-scale flow to stochastic forcing (Farrell and Ioannou 1994, 1995; Whitaker and Sardeshmukh 1998; Zhang and Held 1999). In addition, Bouchet et al. (2013) have shown that, in the limit of weak forcing and dissipation, the formal asymptotic expansion of the probability density function of the Euler equations around a mean flow that is assumed to only have a singular spectrum of modes comprises the second-order S3T closure with an additional stochastic term forcing the mean flow. Therefore, S3T formally describes the statistical equilibrium mean flow and the eddy statistics in this case, as the additional stochastic term only produces fluctuations around this statistical equilibrium. Similar to the S3T closure of the full SSD is the second-order cumulant expansion (CE2) closure in which the third-order cumulant is neglected without parameterization (Marston et al. 2008; Marston 2010, 2012; Tobias and Marston 2013). It has been shown that the predictions of S3T (or CE2) simulations are reflected in corresponding nonlinear simulations (O’Gorman and Schneider 2007; Srinivasan and Young 2012; Tobias and Marston 2013; Constantinou et al. 2014).

The second-order closure results in a nonlinear, autonomous dynamical system that governs the evolution of the mean flow and its consistent second-order perturbation statistics. Its fixed points define statistical equilibria, whose instability brings about structural reconfiguration of the mean flow and of the turbulent statistics. Previous studies employing S3T addressed the bifurcation from a homogeneous turbulent regime to a jet-forming regime in barotropic beta-plane turbulence and identified the emerging jet structures as linearly unstable modes of the homogeneous turbulent state equilibrium (Farrell and Ioannou 2003, 2007; Bakas and Ioannou 2011; Srinivasan and Young 2012; Parker and Krommes 2013, 2014). The S3T stability analysis of the homogeneous equilibrium was further advanced with the introduction of the continuum formulation by Srinivasan and Young (2012), who derived a compact analytic expression for the growth rate and frequency of the unstable structures. Interestingly, Carnevale and Martin (1982) using field theoretic techniques have arrived at the same stability equation for the statistical description of fluctuations about a homogeneous state.

Comparisons of the jet structure predicted by S3T with direct numerical simulations have shown that the structure of zonal flows that emerge in the nonlinear simulations can be predicted by S3T (Srinivasan and Young 2012; Tobias and Marston 2013; Constantinou et al. 2014). However, Srinivasan and Young (2012) found quantitative differences between the predictions of S3T as seen in the bifurcation diagram for the emergence of jets and the corresponding diagram obtained from the nonlinear simulations, calling into question the validity of the S3T (or CE2) approximations when the mean flow is very weak. Constantinou et al. (2014) demonstrated that this discrepancy was due to the prior emergence of nonzonal coherent structures that modified the background equilibrium spectrum and showed that S3T predictions were accurate when this modification in the background spectrum was accounted for.

The nonzonal structures were treated in these studies as incoherent because of the assumption that the ensemble average over the forcing realizations is equivalent to a zonal average and therefore their emergence and effect on the jet dynamics could not be directly addressed. By making the alternative interpretation of the ensemble mean as a Reynolds average over the fast turbulent motions that was introduced in earlier studies of atmospheric blocking (Bernstein 2009; Bernstein and Farrell 2010), Bakas and Ioannou (2013a, 2014) addressed the emergence of the nonzonal coherent structures in barotropic beta-plane turbulence in terms of the parameters and , where *β* is the gradient of the planetary vorticity, is the length scale of the forcing, *ε* is the energy input rate of the forcing, and is the dissipation time scale. Characteristic values of these parameters for Earth’s midlatitude atmosphere and oceans and the Jovian atmosphere are given in Table 1. It was found that for isotropic forcing the homogeneous statistical equilibrium becomes unstable when the energy input rate exceeds a critical value that depends on as shown in the stability regime diagram in Fig. 1. In marginally unstable flows with zonal jets first emerge, while for westward-propagating nonzonal structures first emerge and equilibrate to finite-amplitude traveling waves. At larger energy input rates, the finite-amplitude nonzonal traveling states are unstable and the flow equilibrates to mixed zonal jet–traveling wave states that consist of strong zonal jets with weaker traveling nonzonal structures embedded in them. These predictions of the S3T stability analysis were verified by direct numerical simulations of turbulent barotropic flow (Bakas and Ioannou 2014).

The S3T dynamics that underlie the formation of large-scale structure cannot depend on turbulent anisotropic inverse cascade processes because local-in-wavenumber-space eddy–eddy interactions are absent in S3T. In S3T, large-scale structure emerges from a cooperative instability arising from the nonlocal in wavenumber space interaction between the large-scale mean flow and the forced, small-scale turbulent eddies. The eddy–mean flow dynamics of this cooperative instability has been investigated by Bakas and Ioannou (2013b) for the case of zonal jet emergence in the limit of . It was shown that shear straining of the small-scale eddies by the local shear of an infinitesimal sinusoidal zonal jet, as described by Orr dynamics in a beta plane, produces upgradient fluxes that intensify that zonal jet. In this work we extend the study of this cooperative eddy–mean flow instability to address not only zonal jet formation but also formation of nonzonal coherent structures, and also we will address the formation of these coherent structures for a wide range of values of . We will show that for the eddy–mean flow dynamics of nonzonal structures is also dominated by shearing of the eddies, whereas for resonant and near-resonant interactions play an important role in the dynamics. The importance of near-resonant interactions in the formation of large-scale flows has been previously discussed by Lee and Smith (2007). However, this effect is rigorously quantified in this work. Finally, we discuss the connection between the modulational instability of plane Rossby waves (Lorenz 1972; Gill 1974) and the S3T stability of the homogeneous turbulent equilibrium.

This paper is organized as follows: In section 2 we derive the S3T system for a barotropic flow and the resulting eigenvalue problem addressing the stability of the homogeneous statistical equilibrium. In section 3 we transform the eigenvalue problem to a rotated frame of reference so that the formation of zonal jets and nonzonal structures can be studied under a uniform framework. In section 4 we identify the eddy–mean flow dynamics underlying the S3T instability for isotropic stochastic forcing and in section 5 we study the effect of the forcing anisotropy to the S3T instability. The results are summarized in sections 6 and 7.

## 2. Formulation of S3T dynamics and emergence of nonzonal coherent structures

Consider a barotropic flow on an infinite beta plane with *x* and *y* Cartesian coordinates along the zonal and the meridional direction respectively and with planetary vorticity gradient, . The nondivergent velocity field with components is expressed in terms of the streamfunction, *ψ*, as , where is the unit vector normal to the plane of the flow. The vorticity of the fluid , with , evolves as

where and *J* is the two-dimensional Jacobian, . The flow is dissipated with linear damping at a rate *r*, which typically models Ekman drag in planetary atmospheres. Turbulence is maintained by the external stochastic forcing *ξ*, which models exogenous processes, such as turbulent convection or energy injected by baroclinic instability. We assume that is a temporally delta-correlated and spatially homogeneous random stirring that injects energy at a rate *ε*. We nondimensionalize (1) using the dissipation time-scale and the typical length scale of the stochastic excitation . In these units , , , , , and , where the asterisks denote nondimensional variables. We hereinafter drop the asterisks for simplicity.

To construct the S3T dynamical system in the continuous formulation of Srinivasan and Young (2012), we proceed as follows:

- The averaged fields, denoted with uppercase letters, are calculated by taking a time average, denoted with , over an intermediate time scale, larger than the time scale of the turbulent motions but smaller than the time scale of the large-scale motions. Deviations from the mean (eddies) are denoted with dashes and lowercase letters. For example, the vorticity field is split as , where . The equations for the mean and the eddies that derive from (1) are where is the linear perturbation operator about the instantaneous mean flow and . Neglecting the nonlinear term in (2b) we obtain the quasi-linear system:
- The quasi-linear system (4) under the ergodic assumption that the time average over the intermediate time scale is equal to an ensemble average produces the S3T system: where
*C*is the ensemble-mean eddy-vorticity spatial covariance between points and ;*Q*is the spatial covariance of the delta-correlated and spatially homogeneous forcing, defined by and is the ensemble-mean vorticity forcing of the large scales by the eddy field, given by The subscript*a*(*b*) in the operators indicates that the coefficients of the operator are evaluated on*a*(*b*) and that the operator acts only on the variable (). The subscript indicates that any function of and is evaluated at the same point, .

The S3T system (5) is a closure of the statistical dynamics of (2) at second order. Being autonomous, it may posses statistical equilibria , the stability of which is addressed by considering small perturbations and performing an eigenanalysis of the linearized S3T equations about these equilibria.

For any spatially homogeneous forcing *Q*, there is always the homogeneous S3T equilibrium

with no mean flow and a homogeneous eddy field—that is, with a translationally invariant covariance . The stability of the homogeneous equilibrium (9) is determined from eigenalysis of the linearized S3T equations about this equilibrium:

with [obtained by setting in (3)] and . It can be shown from (10) that the homogeneous equilibrium is S3T stable for and becomes unstable when *ε* exceeds the critical value that depends on *β* and on the structure of *Q*. In this work we focus our analysis close to the instability threshold and identify the physical processes underlying the S3T instability. We follow Srinivasan and Young (2014) and consider a ring stochastic forcing of waves of total wavenumber , with power spectrum:^{1}

where

with and . The parameter *μ* modulates the anisotropy of the spectrum of the forcing and takes values in order that the spectrum is everywhere positive and therefore physically realizable. Example realizations of the spatial structure of stochastic excitations at different *μ* are shown in Fig. 2. When , the forcing is isotropic and could model the forcing of the Jovian atmosphere at cloud level from turbulent convection. When , the stochastic excitation favors small- Fourier components as the baroclinic forcing of the upper-level jet in the midlatitude atmosphere. When , the forcing favors the almost-zonal Fourier components around .

## 3. Emergence of nonzonal structures as zonal flows in a rotated frame

The eigenfunctions of the S3T stability equations shown in (10) are specified by their two components: the mean-flow component and the covariance component . Because the stability equations given in (10) are linearized about the homogeneous equilibrium depicted in (9), the eigenfunction structure simplifies significantly and assumes the form

where is the wavevector that characterizes the eigenfunction and is the homogeneous component of the covariance eigenfunction. The mean-flow component of each eigenfunction has the form of a zonal flow when and the form of a nonzonal flow when . However, nonzonal mean-flow perturbations can be rendered zonal through a rotation of the frame of reference.

For an eigenfunction with wavenumber , clockwise rotation of the axes by an angle transforms the components of to

with and the components of the planetary vorticity gradient to . Correspondingly, in the rotated frame the eigenfunction has only the mean-flow component , which is of the form of a zonal jet in the direction (i.e., the wavevector has zero component; cf. Fig. 3), and . The eigenvalue problem (10) about the homogeneous equilibrium transforms to

with and . The equilibrium covariance in the rotated frame is defined as , where are the components of in the rotated frame. The perturbation vorticity flux is given in terms of as

In writing the S3T eigenvalue problem in the rotated frame and by transforming a nonzonal perturbation into a zonal jet perturbation, there is a twofold gain. The first is that we can use the methods that were previously developed by Bakas and Ioannou (2013b) in the context of the emergence of zonal jets in order to understand the mechanisms responsible for the emergence of nonzonal structures. The second is that we can directly address the eddy–mean flow dynamics that give rise to zonal jets with constant topographic vorticity gradient but in a direction other than the meridional (Boland et al. 2012).

From here on we will study the S3T instability of the homogeneous equilibrium (9) in the rotated frame. While in the rotated frame all eigenfunctions have the form of a zonal jet, we distinguish the perturbations as zonal when and nonzonal when (i.e., as they manifest in the unrotated frame). A direct implication of (15a) is that that the vorticity flux induced by eigenfunction , must be proportional to and it therefore can be written as

with *f* determining the amplitude and relative phase of the vorticity flux feedback induced by the mean-flow eigenfunction with eigenvalue *σ*. When the real part of *f* is positive, the induced vorticity flux is upgradient, and when it is negative, the flux is downgradient. It is shown in appendix A that

with , , and . The power spectrum of the stochastic forcing given in (11) in the rotated frame takes the form

With this notation, the eigenvalues *σ* of (15) satisfy the equation

For , (20) reduces to the eigenvalue relation of Srinivasan and Young (2012) that governs the stability of the homogeneous equilibrium (9) to zonal jet perturbations in the unrotated frame. For mirror symmetric forcing [i.e., with covariance satisfying or ], like (11), and for we find numerically^{2} that the eigenvalue corresponding to unstable zonal jet eigenfunctions have ; that is, the unstable zonal jets grow in situ. For the above expression produces the eigenvalue relation obtained by Bakas and Ioannou (2014) for the growth rate of nonzonal perturbations with wavenumbers in the unrotated frame. The growing eigenfunctions in this case are numerically found to be propagating () and at marginal stability for their frequency becomes the Rossby wave frequency

for a plane wave with wavevector .

From (20), we obtain that a necessary condition for S3T instability is that the real part of the vorticity flux feedback factor must be positive. To illuminate the eddy–mean flow dynamics underlying the S3T instability, we study the behavior of for energy input rates close to . Near the stability boundary and under the assumption that at marginal stability^{3 }, the feedback on the mean flow for the delta function forcing (19) can be written as

where (cf. appendix A) is the contribution to from Fourier components of the forcing with wavevectors and (see Fig. 3). When , the induced vorticity fluxes are upgradient and the critical energy input rate is . The integrand can be alternatively interpreted as the contribution of the stochastically forced waves or eddies to the vorticity fluxes. These forced waves have a total wavenumber and are characterized only by the angle *ϑ* between their phase lines and the axis. We can isolate the dependence of the feedback factor on *β* by writing with

and, as shown in appendix A, functions , , and independent of *β*.

In the following sections we will determine the contribution of the various waves to the vorticity flux feedback and identify the angle *ϑ* that produces the most significant contribution to this feedback. We will also calculate as a function of the mean-flow wavenumber *n* for . We will limit our discussion to the emergence of mean flows with (i.e., with scale larger than the scale of the forcing). In section 4 the analysis is mostly focused to isotropic forcing () while the effect of anisotropy is discussed in section 5.

## 4. Eddy–mean flow dynamics leading to formation of zonal and nonzonal structures for isotropic forcing

### a. Induced vorticity fluxes for

We expand the integrand of (22) in powers of *β*:

with . The leading-order term is the contribution of each wave with wavevector to the vorticity flux feedback in the absence of *β* and is shown in Fig. 4a. For , the dynamics are rotationally symmetric and for isotropic forcing is independent of *φ*. Therefore all zonal and nonzonal eigenfunctions with the same wavenumber *n* grow at the same rate. Upgradient fluxes () to a mean flow with wavenumber *n* are induced by waves with phase lines inclined at angles satisfying (cf. appendix B). This implies that all waves with necessarily produce upgradient vorticity fluxes to any mean flow with wavenumber , while waves with produce upgradient fluxes for any mean flow with large-enough wavenumber (cf. Fig. 4a). The eddy–mean flow dynamics was investigated in the limit of by Bakas and Ioannou (2013b). It was shown that the vorticity fluxes can be calculated by time averaging the fluxes over the life cycle of an ensemble of localized stochastically forced wavepackets initially located at different latitudes. For , the wavepackets evolve in the region of their excitation under the influence of the infinitesimal local shear of and are rapidly dissipated before they shear over. As a result, their effect on the mean flow is dictated by the instantaneous (with respect to the shear time scale) change in their momentum fluxes. Any pair of wavepackets having a central wavevector with phase lines forming angles with the *y* axis surrender instantaneously momentum to the mean flow and reinforce it, whereas pairs with gain instantaneously momentum from the mean flow and oppose jet formation. Therefore, anisotropic forcing that injects significant power into Fourier components with (such as the forcing from baroclinic instability that primarily excites Fourier components with ) produces robustly upgradient fluxes that asymptotically behave antidiffusively. That is, for a sinusoidal mean-flow perturbation we have with *K* positive and proportional to the anisotropy factor *μ* [cf. appendix B and Bakas and Ioannou (2013b)].

For isotropic forcing, the net vorticity flux produced by shearing of the perturbations vanishes (i.e., ), given that the upgradient fluxes produced by waves with exactly balance the downgradient fluxes produced by the waves with . However, a net vorticity flux feedback is produced and asymptotically behaves as a negative fourth-order hyperdiffusion with coefficient for [cf. (25) and Bakas and Ioannou (2013b)]. In appendix B it is shown that the feedback factor for isotropic forcing in the limit with is

which is accurate even up to , as shown in Fig. 5. To understand the contribution of *β* to the vorticity flux feedback, we plot for a zonal (Fig. 4b) and a nonzonal perturbation (Fig. 4c) as a function of the mean-flow wavenumber *n* and wave angle *ϑ*. We choose to scale by because in (25) increases as . Consider first the case of a zonal jet. It can be seen that at every point, has the opposite sign to , implying that *β* tempers both the upgradient (for roughly ) and the downgradient (for ) fluxes of . However, in the sector the values of are much larger than in the sector and the net fluxes integrated over all angles are upgradient, as in (25) for the isotropic case.

The asymptotic analysis of Bakas and Ioannou (2013b), which is formally valid for , offers understanding of the dynamics that lead to the inequality and to the positive net contribution of (i.e., to ). Any pair of wavepackets with wavevectors at angles instantaneously gain momentum from the mean flow as described above (i.e., for ), but their group velocity is also increased (decreased) while propagating northward (southward). This occurs because shearing changes their meridional wavenumber and consequently their group velocity. The instantaneous change in the momentum fluxes resulting from this speed up (slow down) of the wavepackets is positive in the region of excitation leading to upgradient fluxes (). The opposite happens for pairs with [cf. Fig. 3 in Bakas and Ioannou (2013b)]; however, the downgradient fluxes produced are smaller than the upgradient fluxes, leading to a net positive contribution when integrated over all angles. Figure 4b shows that this result is valid for larger mean-flow wavenumbers as well.

Consider now the case of a nonzonal perturbation (Fig. 4c). We observe that the angles for which the waves have significant positive or negative contributions to the vorticity flux feedback are roughly the same as in the case of zonal jets. In addition, the vorticity flux feedback factor decreases with the angle *φ* of the nonzonal perturbations [cf. (25)]. As a result, zonal jet perturbations always produce larger vorticity fluxes compared to nonzonal perturbations and are therefore the most unstable in the limit . Additionally, these results show that, for , the mechanism for structural instability of the nonzonal structures is the same as the mechanism for zonal jet formation, which is shearing of the eddies by the infinitesimal mean flow.

### b. Induced vorticity fluxes for

Consider first the emergence of nonzonal structures in the limit . The contribution of each Fourier component of the forcing to the vorticity flux feedback *f _{r}* for the case of nonzonal structures at is shown in Fig. 6a. In contrast to the cases with [or , discussed in section 4c], there is only a small band of Fourier components that contribute significantly to the vorticity flux feedback, as indicated with the narrow tongues in Fig. 6a. The reason for this selectivity in the response is that for the components that produce appreciable fluxes, as seen from (23), are concentrated on the curves that satisfy (shown in Fig. 6b) or equivalently for the that satisfy the resonant condition [cf. (A9)]. This is the resonant condition satisfied when a Rossby wave with wavevector and frequency forms a resonant triad with the nonzonal structure with wavevector and frequency . We concentrate our analysis to these “resonant contributions” because they dominate the vorticity flux feedback of nonzonal perturbations for .

Resonant triads do not occur for all mean-flow perturbations . For in region D of Fig. 7a, has no roots and therefore there are no Fourier components with that form a resonant triad with and is determined by the sum over the nonresonant contributions as illustrated in Fig. 7b. In region B of Fig. 7a, there are only two resonant angles *ϑ*. The resonant and nonresonant contribution for a typical case in region B is shown in Fig. 7c. Note that it is the resonant contributions that determine the vorticity flux feedback. However, they produce a negative vorticity flux feedback (a downgradient tendency), which is stabilizing, a result that holds for all in region B. In regions A and C, there exist four resonant angles that dominate the vorticity flux. In C, all resonant contributions are stabilizing and therefore C is also a stable region. In region A, which at most extends to (cf. appendix B), two of the four resonances give positive contributions to (cf. Figs. 7d,e). Therefore, only for in region A does a destabilizing vorticity flux feedback occur. The largest destabilizing feedback occurs when the positively contributing resonances are near coalescence (i.e., as in Fig. 7d), which occurs for close to the curve separating regions A and B. The reason is that when the resonances are apart, as in Figs. 7c and 7e, the significant contributions come from near-resonant waves with angles within a band of around the resonant angles and the integrated resonant contributions to the vorticity flux are . However, when the resonances are near coalescence, as for the case shown in Fig. 7d, the band of near-resonant waves contributing significantly increases as the integrand assumes a double-humped shape and, as shown in Appendix B, the destabilizing vorticity flux feedback becomes . Note that as , the width over which we have significant contributions diminishes and therefore unless an infinite amount of energy is injected exactly at the resonant angles (as is assumed in modulational instability studies).

It can be shown (cf. appendix B) that the resonant contribution for asymptotically approaches

where the subscript *j* indicates the value of the functions at the *j*th of the *M* roots of and . The values , , and are all , whereas is always positive and the only quantity that has dependence on *β*. It is only for just above the separating boundaries of regions A and B and regions B and D in Fig. 7a yielding and is elsewhere yielding , as is also qualitatively described above. The sign of the *j*th resonant contribution to the total vorticity flux feedback depends only on the sign of . For just above the boundary separating regions B and D, and attains its minimum value, which corresponds to the largest stabilizing tendency. This is illustrated in Fig. 8, showing as a function of *n*. For just above the boundary separating regions A and B, coalescence of the two positive contributing resonances occurs and attains its maximum value, which corresponds to the largest destabilizing tendency. For small mean-flow wavenumbers *n* (corresponding to region D), the vorticity flux feedback is negative and owing to the absence of resonant contributions.

An interesting exception to the results discussed above occurs for the important case of zonal jet perturbations (). In that case, in (26) as the roots of and coincide and the resonant contribution (26) is exactly zero. As shown in Fig. 9, positive vorticity flux feedback is obtained from a broad band of the nonresonant Fourier components with , corresponding to waves with lines of constant phase nearly aligned with the *y* axis (remember that, for smaller *β*, the region that produces destabilizing fluxes extends up to ). For large *β*, is always destabilizing for all zonal jet perturbations with , as shown by (B18) and Fig. 8, and the largest destabilizing vorticity flux feedback, , is obtained for jets with the largest allowed scale. The reason for the weak fluxes and the preference for the emergence of jets of the largest scale in this limit is understood by noting that the stochastically forced eddies for propagate with group velocities. Therefore, in contrast to the limit of in which they evolve according to their local shear, the forced waves will respond to the integrated shear of the sinusoidal perturbation over their large propagation extent, which will be very weak.

To summarize, although zonal jets and most nonzonal perturbations induce fluxes that decay as for large *β*, resonant and near-resonant interactions arrest the decay rate of certain nonzonal perturbations by a factor of , leading to fluxes that decay as . This makes the nonzonal perturbations to be the most S3T unstable perturbations for . Also in contrast to when is positive for all *n* and *φ* (cf. Fig. 5), the vorticity flux feedback is negative for in regions B and D of Fig. 7a. As a result, the mean flows that produce negative fluxes and are by necessity S3T stable are interestingly in the interior of the dumbbell shape shown in Fig. 10. The largest destabilizing fluxes occur in the narrow region adjacent to the outer boundaries of the dumbbell shape, which demarcates the boundary separating regions A and B. Because of the selectivity of the resonances these results do not depend on the forcing anisotropy as we will see in the next section.

### c. Induced vorticity fluxes for

We have seen that in the singular case of isotropic forcing the only process available for the emergence of mean flows is the fourth-order antidiffusive vorticity flux feedback induced by the variation of the group velocity of the forced eddies due to the mean flow shear. For , the waves interact with the local shear producing fluxes proportional to . As *β* increases, this growth is reduced since the waves interact with an effective integral shear within their propagation extent, which is weak and, eventually, as we have seen in the previous section, for , the fluxes decay as . Therefore, the fluxes attain their maximum at an intermediate value of *β*. This occurs for , as can be seen in Fig. 11a where the maximum over all is shown. It will be demonstrated in the next section that this intermediate *β* maximizes the S3T instability for all forcing spectra.

While the eddy–mean flow interaction of both zonal and nonzonal perturbations is dominated by the same dynamics when , for the eddy–nonzonal flow interaction is dominated by resonances that do not occur for zonal flow perturbations. The resonant interactions lead to the possibility of arrested decay of the vorticity flux at the rates of and instead of the decay in the absence of resonances. The vorticity flux attains its maximum at an intermediate value for nonzonal mean flows as well, which is nonetheless large enough for the resonant contributions to reinforce the contribution from the shearing mechanism. Figure 12 shows the contribution to the vorticity flux feedback induced by the various wave components that are excited for two values of *β* ( and ) in the case of zonal jets () and nonzonal perturbations (). As *β* increases, the resonant contributions start playing an important role for nonzonal perturbations as there is enhanced contribution to the vorticity flux feedback in the vicinity of the curves, indicated by the white dashed lines. These resonant contributions enhance the vorticity fluxes relative to the fluxes obtained for zonal jets and render the nonzonal structures more unstable compared to zonal jets when (Bakas and Ioannou 2014).

## 5. Effect of anisotropic forcing on S3T instability

In this section we investigate the effect of the anisotropy of the excitation on the S3T instability. The maximum vorticity flux feedback for three cases of anisotropy ( and ) and for isotropic forcing () is shown in Fig. 11a. For , the main contribution to for zonal jet perturbations comes from forced waves with nearly meridional constant phase lines (angles near ; cf. Fig. 9). Therefore, attains larger (smaller) values for a stochastic forcing that injects more (less) power in waves with angles near —that is, for positive (negative) anisotropicity factor *μ* (cf. Fig. 2). The maximum value of over all wavenumbers *n* depends in this case linearly on *μ* (cf. appendix B),

For nonzonal perturbations, the main contribution comes from forced waves satisfying the resonant condition and depends only on the sum of the resonant contributions. The sign of that determines whether the resonant contribution is positive or negative [cf. (26)] depends only on the sign of and not on *μ* [cf. (A7c)]. The anisotropicity affects only the magnitude of . For any it is found that the resonances giving positive contribution occur at angles for which . A stochastic excitation, which injects more power near () will give larger positive resonant contributions and therefore increases with *μ*. However, the effect on the maximum vorticity feedback is weak, as the spectral selectivity of the resonances renders the characteristics of the most unstable nonzonal structure independent of the spectrum of the forcing. That is, that correspond to the maximum asymptotically approaches to , (marked with star in Fig. 7a) as , a result that is very weakly dependent on *μ* (cf. Figs. 11b,c).

For , the characteristics of the S3T instability are dependent on the anisotropy of the stochastic forcing. The vorticity flux feedback is at leading order proportional to *μ*:

This shows that there can be upgradient vorticity fluxes leading to S3T instability for as long as . For , the maximum is achieved by zonal jets (), while for any nonzonal perturbation with can grow, with the maximum achieved for when the nonzonal perturbations assume the form of jets in the *y* direction (meridional jets) (cf. Fig. 11c).

It is worth noting that Srinivasan and Young (2014) also find that that the eddy momentum fluxes are proportional to *μ* when a constant shear flow is stochastically forced with power spectrum (11). This result is intriguing as the two studies address two different physical regimes. This work treats the limit appropriate for emerging structures in which the shear time scale is far larger than the dissipation time scale with the fluxes determined by the instantaneous response of the eddies on the shear. Srinivasan and Young (2014) study the opposite limit in which the mean-flow shear is finite and the shear time scale is much shorter than the dissipation time scale with the fluxes determined by the integrated influence of the shear on the eddies over their whole life cycle, which may include complex effects such as reflection and absorption at critical levels.

In summary,

the S3T instability of the homogeneous state is a monotonically increasing function of

*μ*for all*β*,the forced waves that contribute most to the instability are structures with small

*γ*(i.e., waves with phase lines nearly aligned with the*y*axis, as in Fig. 2a), andthe anisotropy of the excitation affects prominently the S3T stability of the homogeneous state only for .

## 6. Discussion

In this work we addressed the dynamics underlying the onset of the S3T instability leading to the formation of large-scale structure but not the nonlinear development and equilibration of the instability. The emergent structure may be susceptible to either hydrodynamic or structural secondary instabilities as it reaches finite amplitude [cf. Farrell and Ioannou 2003, 2007; Parker and Krommes (2014) for zonal jets and Bakas and Ioannou (2014) for nonzonal flows]. For example, the most unstable jet structure for marginally unstable parameters is at the scale of the forcing [for small *β*, maximum instability occurs at a scale slightly larger than for isotropic forcing and close to for anisotropic forcing^{4} while for larger *β*, maximum instability occurs at about (cf. Fig. 11b)]. Most geophysical flows are far from marginal stability and the jet scale predicted at marginal stability does not characterize the scale of the actual jets. S3T predicts that the emergent jets through a series of mergers usually equilibrate at a much larger scale (Farrell and Ioannou 2007; Parker and Krommes 2014). These predictions of S3T have been shown to be accurately reflected in sample nonlinear simulations (Srinivasan and Young 2012; Constantinou et al. 2014).

Although in this work we examined the statistical dynamical instability of a homogeneous state of turbulence in the presence of forcing and dissipation, the results bear a relation to the deterministic barotropic hydrodynamic instability of nonzonal flows on a *β* plane in the absence of forcing and dissipation. Parker and Krommes (2015) have recently shown that in the inviscid limit the modulational instability of a Rossby wave (Lorenz 1972; Gill 1974; Connaughton et al. 2010) and the S3T instability of a homogeneous turbulent state with equilibrium vorticity power spectrum corresponding to the Rossby wave: obey the same stability condition. This equivalence is formal because physically the two problems are very different. In the problem of Lorenz (1972), the stability of a basic state in the form of a coherent Rossby plane wave is studied, while S3T addresses the statistical stability of an incoherent homogeneous state with the power spectrum of the Rossby wave. In that sense, as noted also by Parker and Krommes (2015), S3T stability analysis embeds the modulational instability results into a more general physical framework. In appendix C we extend the result of Parker and Krommes (2015) and show the formal equivalence between the modulational instability of any solution of the barotropic equation, which may be time dependent in general but has stationary power spectrum, with the S3T instability of the homogeneous state with the same power spectrum. Such a nonlinear solution of the inviscid barotropic vorticity equations is for example a superposition of any number of Rossby waves:

all with the same total wavenumber, , that forms a nondispersive structure moving westward (cf. appendix C). If we assume a zonal jet perturbation superimposed on this nonlinear solution, then the results in this work show further that the dynamics underlying the instability of this structure can be interpreted in the limit of as shearing of the finite-amplitude solution by the weak shear of the jet perturbation.

## 7. Conclusions

The mechanism for formation of coherent structures in a barotropic beta plane under a spatially homogeneous and temporally delta-correlated stochastic forcing was examined in the framework of stochastic structural stability theory (S3T). In this framework, a second-order closure for the dynamics of the flow statistics is obtained by ignoring or parameterizing the eddy–eddy nonlinearity. The resulting deterministic system for the joint evolution of the coherent flow and of the second-order turbulent eddy covariance admits statistical equilibria.

For a spatially homogeneous forcing covariance, a homogeneous state with no mean coherent structures is such an equilibrium solution of the S3T dynamical system. When a critical energy input rate of the forcing is exceeded, this homogeneous equilibrium is unstable and propagating nonzonal coherent structures and/or stationary zonal jets emerge in agreement with direct numerical simulations. To identify the processes that lead to the formation of coherent structures, the vorticity fluxes induced by a plane wave-mean flow, which is the eigenfunction of the linearized S3T system around the homogeneous equilibrium, were calculated close to the bifurcation point and closed form asymptotic expressions for these fluxes were obtained. Upgradient fluxes in this limit are consistent with S3T instability and coherent structure formation.

The induced fluxes were calculated in a rotated frame of reference, in which the plane wave-mean flow corresponds to a zonal jet evolving in a beta plane with a nonmeridional planetary vorticity gradient. This was done because in this rotated frame of reference the intuition gained by previous studies for the eddy–mean flow dynamics underlying zonal jet formation can be utilized to clarify the dynamics underlying nonzonal wave formation, or formation of zonal jets when the effect of topography is equivalent to a nonmeridional planetary vorticity gradient.

In the limit of a weak planetary vorticity gradient (), the eddy–mean flow dynamics are similar for both zonal jets and nonzonal structures. The stochastically forced eddies that propagate with slow group velocities in this limit are rapidly dissipated as they are sheared over by the infinitesimal mean flow. Their effect on the mean flow is therefore determined at leading order by the instantaneous, with respect to the shear time scale, change in their momentum fluxes and to second order by the instantaneous change in their group velocity. The waves with constant phase lines that form angles with the meridional direction instantaneously surrender momentum to the mean flow and lead to upgradient fluxes that reinforce the mean flow for an anisotropic forcing. For an isotropic forcing, this leading-order effect produces no net fluxes when integrated over all forced waves and the instability is controlled by the second-order effect that the instantaneous change of the waves’s group velocity has on the momentum fluxes. In this case, the group velocity of waves that form angles with the meridional direction is instantaneously increased (decreased) for waves propagating northward (southward) because of refraction. The difference in momentum fluxes resulting from this change in group velocity is positive in the region of their excitation leading to upgradient fluxes. As a result, the anisotropy of the forcing has a significant effect on the induced fluxes and the S3T instability in this limit. In any case, the effect of the eddies on the mean flow due to shearing is larger for zonal jets compared to nonzonal perturbations and consequently zonal jets are more unstable in this limit.

In the limit of strong planetary vorticity gradient , the eddy–mean flow dynamics producing upgradient vorticity fluxes are different for zonal and nonzonal perturbations, but in both cases the fluxes decrease with *β*. Zonal jets continue to induce upgradient vorticity fluxes through wave shearing, which decrease as . The reason is that in this limit the waves that can propagate in the meridional direction are influenced by the integrated shear over the sinusoidal flow, which is very small. However, the nonzonal mean-flow perturbations can sustain fluxes that decrease only as . The reason for these larger fluxes is that resonant and near-resonant interactions dominate the dynamics in this limit [cf. section 3.26 in Pedlosky (1992)]. Resonance occurs between the emerging structure, which, close to the stability boundary, satisfies the Rossby wave dispersion, and the stochastically forced waves satisfying the Rossby wave frequency resonant condition. The resonant interactions that occur for nonzonal structures may produce upgradient or downgradient net fluxes and it was found that upgradient fluxes cannot be induced by nonzonal flows with wavenumbers in a region of wavenumber space in the shape of a dumbbell. Maximum upgradient fluxes occur for both zonal and nonzonal flows for . In this regime, shearing of the forced waves by the infinitesimal nonzonal flows is reinforced by fluxes from the resonant interactions, enhancing the vorticity fluxes and rendering the nonzonal structures more unstable compared to zonal jets when . In contrast to the limit , these results were found to be insensitive to the anisotropy of the forcing.

Finally, the relation of the S3T instability and modulational instability of finite-amplitude Rossby waves was discussed. Parker and Krommes (2015) showed that the growth rates obtained when three Rossby waves interact with the primary finite-amplitude Rossby wave match exactly in the inviscid limit that the growth rates obtained by the S3T stability analysis for the homogeneous equilibrium with the vorticity covariance produced by the primary Rossby wave. It was shown in this work that this agreement can be found for more general cases (e.g., when the covariance is produced by any linear combination of Rossby waves with the same total wavenumber). Such an agreement occurs because retaining only the interaction between four waves in modulational instability is equivalent to neglecting the eddy–eddy nonlinearity in S3T. The equivalence of the dynamics underlying modulational and S3T instability in this case shows that S3T stability analysis generalizes modulational instability analysis to stochastically forced dissipative flows. However, contrary to modulational instability, the underlying S3T dynamics can capture both the emergence of large-scale structure and its equilibration. In addition, the dynamics underlying modulational instability can be interpreted under the alternative eddy–mean flow view adopted in this work.

## Acknowledgments

Nikolaos Bakas is supported by the AXA Research Fund and Navid Constantinou acknowledges the support of the Alexander S. Onassis Public Benefit Foundation. Nikolaos Bakas would also like to thank Izumi Saito for fruitful discussions on the results of this work.

### APPENDIX A

#### Eddy-Vorticity Flux Response to a Mean-Flow Perturbation

In this appendix we study the eigenvalue problem (15), which determines the S3T stability of jet perturbations to the homogeneous turbulent equilibrium (9) in the rotated frame of reference. The eigenfunction corresponding to eigenvalue *σ* has the spatial structure

The power spectrum of the homogeneous part of the covariance eigenfunction, , is determined from (15b) to be

with , , , , , and the Fourier transform of the forcing covariance (19). The vorticity flux induced by this eigenfunction is

which is proportional to . By using the symmetry , which implies that , and by changing the integration variable in (A3) to , we obtain the following expression for the feedback factor *f*:

with and .

Introducing (A3) into (15a) we obtain the stability equation (20) that determines *σ*, which can be shown to be exactly the stability equation obtained by Bakas and Ioannou (2014). The stability equation can be written in terms of the real and imaginary part of *σ* as

The real part of the vorticity flux feedback contributes to the growth rate of the mean flow and the imaginary part determines the departure of the phase speed of the mean flow from the Rossby wave frequency . For the first term in (A5b) is while for marginally unstable eigenfunctions is at most of . As will be shown, the critical *ε* increases as or as *β* and therefore the frequency of the marginally unstable waves is approximately equal to the Rossby phase frequency.

We focus on the real part of the feedback gain near marginal stability (). Setting in (A5a) for the marginally unstable structures and for the ring forcing in the rotated frame, , takes the form

with

Positive values of indicate that the vorticity flux induced by the stochastic forcing at marginal stability on the mean flow with wavenumber *n* and nonzonality parameter *φ* is upgradient, and the marginal energy input rate is .

Note that , which is defined in (23) as the rhs of (A6), is unchanged when the angle *φ* is shifted by () or when there is a simultaneous shift of and . As a result, it suffices to only consider cases with .

As Parker and Krommes (2015) noted, the stability equation (20) can be written in coordinate independent form as

where is the Rossby frequency of a wave with wavenumber [defined in (21)]. As a result, in coordinate-free form

and the roots of on the plane satisfy the resonant condition

### APPENDIX B

#### Asymptotic Expressions for the Induced Vorticity Flux Feedback

In this appendix we calculate in closed form asymptotic expressions for the vorticity flux feedback induced by a mean-flow perturbation in the form of a zonal jet in the rotated frame of reference with wavenumber *n* for and .

##### a. Case

When and for *n* satisfying , we expand