## Abstract

Stochastic structural stability theory (S3T) provides analytical methods for understanding the emergence and equilibration of jets from the turbulence in planetary atmospheres based on the dynamics of the statistical mean state of the turbulence closed at second order. Predictions for formation and equilibration of turbulent jets made using S3T are critically compared with results of simulations made using the associated quasi-linear and nonlinear models. S3T predicts the observed bifurcation behavior associated with the emergence of jets, their equilibration, and their breakdown as a function of parameters. Quantitative differences in bifurcation parameter values between predictions of S3T and results of nonlinear simulations are traced to modification of the eddy spectrum which results from two processes: nonlinear eddy–eddy interactions and formation of discrete nonzonal structures. Remarkably, these nonzonal structures, which substantially modify the turbulence spectrum, are found to arise from S3T instability. Formation as linear instabilities and equilibration at finite amplitude of multiple equilibria for identical parameter values in the form of jets with distinct meridional wavenumbers is verified, as is the existence at equilibrium of finite-amplitude nonzonal structures in the form of nonlinearly modified Rossby waves. When zonal jets and nonlinearly modified Rossby waves coexist at finite amplitude, the jet structure is generally found to dominate even if it is linearly less unstable. The physical reality of the manifold of S3T jets and nonzonal structures is underscored by the existence in nonlinear simulations of jet structure at subcritical S3T parameter values that are identified with stable S3T jet modes excited by turbulent fluctuations.

## 1. Introduction

Spatially and temporally coherent jets are a common feature of turbulent flows in planetary atmospheres with the banded winds of the giant planets constituting a familiar example (Vasavada and Showman 2005). Fjørtoft (1953) noted that the conservation of both energy and enstrophy in dissipationless barotropic flow implies that transfer of energy among spatial spectral components results in energy accumulating at the largest scales. This argument provides a conceptual basis for understanding the observed tendency for formation of large-scale structure from small-scale turbulence in planetary atmospheres. However, the observed large-scale structure is dominated by zonal jets with specific form and, moreover, the scale of these jets is distinct from the largest scale in the flow. Rhines (1975) argued that the observed spatial scale of jets in beta-plane turbulence results from arrest of upscale energy transport at the length scale , where *β* is the meridional gradient of planetary vorticity and *u* is the root-mean-square velocity in the turbulent fluid. In Rhines’s interpretation this is the scale at which the turbulent energy cascade is intercepted by the formation of propagating Rossby waves. Balk et al. (1991) extended Rhines’s argument by showing that in addition to energy and enstrophy, dissipationless barotropic turbulence conserves a third quadratic invariant called zonostrophy, which constrains the large-scale structures in dissipationless beta-plane turbulence to be predominantly zonal (cf. Balk and Yoshikawa 2009). While these results establish a conceptual basis for expecting large-scale zonal structures to form in beta-plane turbulence, the physical mechanism of jet formation, the structure of the jets, and their dependence on parameters remain to be determined.

One mechanism for formation of jets is vorticity mixing resulting from Rossby wave breaking, which leads to homogenization of vorticity in localized regions and formation of vorticity staircases. The risers of these staircases correspond to thin prograde jets located at the latitudes of steep vorticity gradients separating parabolic retrograde jets corresponding to the well-mixed steps of the staircase (Baldwin et al. 2007; Dritschel and McIntyre 2008). While vorticity staircases have been obtained in numerical simulations (Scott and Dritschel 2012), in many cases mixing is insufficient to produce a staircase structure. Moreover, jets are shown here to form from a bifurcation at infinitesimal perturbation amplitude and in the absence of wave breaking.

Arguments based on equilibrium statistical mechanics have also been advanced to explain emergence of jets (e.g., by Miller 1990 and Robert and Sommeria 1991). This theory is based on the principle that dissipationless turbulence tends to produce configurations that maximize entropy while conserving both energy and enstrophy. These maximum entropy configurations in beta-plane barotropic turbulence assume the form, when observed at large scale, of zonal jets or vortices (cf. Bouchet and Venaille 2012). However, the relevance of these results to planetary flows that are strongly forced and dissipated and therefore out of equilibrium remains to be shown.

An important constraint on theories of jet maintenance is that the primary mechanism by which planetary turbulent jets are maintained is eddy momentum flux systematically directed up the mean velocity gradient. This upgradient momentum flux is produced by a broad spectrum of eddies, implying that the large-scale jets are maintained by spectrally nonlocal interaction between the eddy field and the large-scale zonal jets. This has been verified in observational studies on Jovian atmosphere (Ingersoll et al. 2004; Salyk et al. 2006) and in numerical simulations (Nozawa and Yoden 1997; Huang and Robinson 1998). Wordsworth et al. (2008) studied jet formation in rotating tanks and found strong evidence confirming that jets are maintained by spectrally nonlocal energy transfer.

Laminar instability of a meridional Rossby wave or of a zonally varying meridional flow can generate zonal flows (Lorenz 1972; Gill 1974; Manfroi and Young 1999; Berloff et al. 2009; Connaughton et al. 2010). Equations for the dynamics of these jets in the weakly nonlinear limit were obtained by Manfroi and Young (1999). This instability, referred to as modulational instability, produces spectrally nonlocal transfer to the zonal flow from the forced meridional waves but presumes a constant source of these finite amplitude meridional waves. In baroclinic flows, baroclinic instability has been advanced as the source of these coherent waves (Berloff et al. 2009).

Stochastic structural stability theory (SSST, contracted to S3T) addresses turbulent jet dynamics as a two-way interaction between the mean flow and its consistent field of turbulent eddies (Farrell and Ioannou 2003). Both S3T and modulational instability involve nonlocal interactions in wavenumber space but these theories differ in that in S3T the mean flow is supported by its interaction with a broad turbulence spectrum rather than with specific waves. In fact, S3T is a nonequilibrium statistical theory that provides a closure comprising a dynamics for the evolution of the mean flow together with its consistent field of eddies. In S3T, the dynamics of the turbulence statistics required by this closure are obtained from a stochastic turbulence model (STM), which provides accurate eddy statistics for the atmosphere at large scale (Farrell and Ioannou 1993, 1994, 1995; Zhang and Held 1999). Marston et al. (2008) have shown that the S3T system is obtained by truncating the infinite hierarchy of cumulant expansions to second order and they refer to the S3T system as the second-order cumulant expansion (CE2). In S3T, jets initially arise as a linear instability of the interaction between an infinitesimal jet perturbation and the associated eddy field and finite-amplitude jets result from nonlinear equilibria continuing from these instabilities. Analysis of this jet formation instability determines the bifurcation structure of the jet formation process as a function of parameters. In addition to jet formation bifurcations, S3T predicts jet breakdown bifurcations, as well as the structure of the emergent jets, the structure of the finite amplitude equilibrium jets to which they continue, and the structure of the turbulence accompanying the jets. Moreover, S3T is a dynamics so it predicts the time-dependent trajectory of the statistical-mean turbulent state as it evolves and, remarkably, the mean turbulent state is often predicted by S3T to be time dependent in the sense that the statistical mean state of the turbulence evolves in a manner predicted by the theory (Farrell and Ioannou 2009b). The formation of zonal jets in planetary turbulence was studied as a bifurcation problem in S3T by Farrell and Ioannou (2003, 2007, 2008, 2009a,c), Bakas and Ioannou (2011), and Srinivasan and Young (2012). A continuous formulation of S3T developed by Srinivasan and Young (2012) has facilitated analysis of the physical processes that give rise to the S3T instability and construction of analytic expressions for the growth rates of the S3T instability in homogeneous beta-plane turbulence (Srinivasan and Young 2012; Bakas and Ioannou 2013b).

Relating S3T to jet dynamics in fully nonlinear turbulence is facilitated by studying the quasi-linear model, which is intermediate between the nonlinear model and S3T. The quasi-linear (QL) approximation to the full nonlinear dynamics (NL) results when eddy–eddy interactions are not explicitly included in the dynamics but are either neglected entirely or replaced with a simple stochastic parameterization, so that no turbulent cascade occurs in the equations for the eddies, while interaction between the eddies and the zonal-mean flow is retained fully in the zonal-mean equation. S3T is essentially QL with the additional assumption of an infinite ensemble of eddies replacing the single realization evolved under QL. Although the dynamics of S3T and QL are essentially the same, by making the approximation of an infinite ensemble of eddies, the S3T equations provide autonomous and fluctuation-free dynamics of the statistical-mean turbulent state, which transforms QL from a simulation of turbulence into a predictive theory of turbulence.

A fundamental attribute of QL–S3T is that the nonlinear eddy–eddy cascade of NL is suppressed in these systems. It follows that agreement in predictions of jet formation and equilibration between NL and QL–S3T provides compelling evidence that cascades are not required for jet formation and theoretical support for observations showing that the turbulent transfers of momentum maintaining finite amplitude jets are nonlocal in spectral space.

Previous studies demonstrated that unstable jets maintained by body forcing can be equilibrated using QL dynamics (Schoeberl and Lindzen 1984; DelSole and Farrell 1996; O’Gorman and Schneider 2007; Marston et al. 2008). In contrast to these studies, in this work, we investigate the spontaneous emergence and equilibration of jets from homogeneous turbulence in the absence of any coherent external forcing at the jet scale. S3T predicts that infinitesimal perturbations with zonal jet form organize homogeneous turbulence to produce systematic upgradient fluxes giving rise to exponential jet growth and eventually to the establishment of finite-amplitude equilibrium jets. Specifically, the S3T equations predict initial formation of jets by the most unstable eigenmode of the linearized S3T dynamics. In agreement with S3T, Srinivasan and Young (2012) found that their NL simulations exhibit jet emergence from a homogeneous turbulent state with subsequent establishment of finite-amplitude jets, while noting quantitative differences between bifurcation parameter values predicted by S3T and the parameter values for which jets were observed to emerge in NL. Tobias and Marston (2013) also investigated the correspondence of CE2 simulations of jet formation with corresponding NL simulations and found that CE2 reproduces the jet structure, although they noted some differences in the second cumulant, and suggested a remedy by inclusion of higher cumulants.

In this paper, we use NL and its QL counterpart together with S3T to examine further the dynamics of emergence and equilibration of jets from turbulence. Qualitative agreement in bifurcation behavior among these systems, which is obtained for all the spatial turbulence forcing distributions studied, confirms that the S3T instability mechanism is responsible for the formation and equilibration of jets. Quantitative agreement is obtained for bifurcation parameters between NL and QL–S3T when account is taken of the modification of the turbulent spectrum that occurs in NL but not in QL–S3T. Remarkably, a primary component of this spectral modification can itself be traced to S3T instability, but of nonzonal rather than of zonal form. We investigate the formation and equilibration of these nonzonal S3T instabilities and the effect these structures have on the equilibrium spectrum of beta-plane turbulence. We also investigate circumstances under which nonzonal structures are modified and suppressed by the formation of zonal jets.

A dynamic of potential importance to climate is the possibility of multiple equilibria of the statistical-mean turbulent state being supported with the same system parameters (Farrell and Ioannou 2003, 2007). Multiple jet equilibria were recently related to pattern formation in the context of S3T by Parker and Krommes (2014). We verify existence of multiple equilibria, predicted by S3T, in our NL simulations. Finally, we show that weak jets result from stochastic excitation by the turbulence of stable S3T modes, which demonstrates the physical reality of the stable S3T modes. Turbulent fluctuation induced excitation of these stable S3T jet modes and the weak but zonally extended jets that equilibrate from weak instabilities at slight supercriticality may explain the enigmatic latent jets of Berloff et al. (2011).

## 2. Formulation of nonlinear and quasi-linear barotropic beta-plane dynamics

Consider a beta plane with *x* and *y* Cartesian coordinates along the zonal and the meridional direction, respectively. The nondivergent zonal and meridional velocity fields are expressed in terms of a streamfunction *ψ* as *u* = −∂_{y}*ψ* and *υ* = ∂_{x}*ψ*. The absolute vorticity is *ζ* + 2Ω + *βy*, where *ζ* = Δ*ψ* and . The NL dynamics of this system is governed by the barotropic vorticity equation:

The flow is dissipated with linear damping at rate *r* and hyperviscosity with coefficient *ν*_{4}. Periodic boundary conditions are imposed in *x* and *y* with periodicity 2*πL*. Distances have been nondimensionalized by *L* = 5000 km and time by *T* = *L*/*U*, where *U* = 40 m s^{−1}, so that the time unit is *T* = 1.5 days and *β* = 10 corresponds to a midlatitude value. Turbulence is maintained by stochastic forcing with spatial and temporal structure *F* and amplitude *ε*.

Flow fields are decomposed into zonal-mean components, denoted with an overbar, and deviations from the zonal mean (eddies), which are indicated lowercase with primes. The zonal velocity is *U*(*y*, *t*) + *u*′(*x*, *y*, *t*), with *U* the zonal-mean velocity, the meridional velocity is *υ*′(*x*, *y*, *t*), and the eddy vorticity is *ζ*′(*x*, *y*, *t*). From Eq. (1), equations for the evolution of the zonal-mean flow *U* and the associated eddy field *ζ*′ obtained by zonal averaging are

where is the nonlinear term representing the eddy–eddy interactions. Equations (2) define the nonlinear system, NL. Note that the stochastic forcing appears only in Eq. (2b).

The QL approximation of NL is obtained by setting *F*_{e} = 0, which implies neglect of the eddy–eddy interactions in Eq. (2b), while retaining the Reynolds stress forcing in the zonal-mean flow equation:

## 3. S3T formulation of barotropic beta-plane dynamics

The S3T system governs evolution of the ensemble-mean state of the QL system [Eqs. (3)]. Derivation and properties of S3T can be found in Farrell and Ioannou (2003, 2007). Expressed in matrix form, the S3T system is

In these equations, the eddy fields and the forcing have been expanded in zonal harmonics—that is, —with *k* = 1, … *N*_{k} the zonal wavenumbers. The zonal-mean flow **U** and the *k* Fourier components of the eddy vorticity form column vectors with elements of their corresponding values at the *N*_{y} discretization points *y*_{j} for *j* = 1, … , *N*_{y}. The second-order statistics of the eddy field are specified by the *N*_{k} covariances , with the angle brackets denoting an ensemble average over realizations and the superscript dagger denoting the Hermitian transpose. The operator vecd() returns the column vector of the diagonal elements of matrix and returns the imaginary part. The matrix is

where **Δ**_{k} = ^{2} − *k*^{2}, ^{2} is the discretized ∂_{yy} operator, is the identity matrix, is the inverse of **Δ**_{k}, is the diagonal matrix with vecd() = , and _{yy} is the diagonal matrix with vecd(_{yy}) = ^{2}. The forcing amplitude is controlled by the parameter *ε* and the spatial covariance of the forcing enters Eq. (4b) as , with = *F*_{kp}(*y*_{j}), where the Fourier component of the stochastic forcing is assumed to have the form (cf. appendix B).

## 4. Stochastic forcing structure

Because the S3T instability mechanism that results in jet bifurcation from a homogeneous turbulent state differs for isotropic and nonisotropic turbulence, we consider examples of both isotropic and nonisotropic turbulence forcing. The jet forming instability in the case of homogeneous, nonisotropic forcing arises from the upgradient fluxes induced by shearing of the turbulence by the infinitesimal perturbation jet, while the upgradient fluxes for the case of homogeneous isotropic forcing arise from the refraction of the eddies caused by the variation in the potential vorticity gradient induced by the infinitesimal perturbation jet (Bakas and Ioannou 2013b).

Three stochastic forcing structures will be used in our investigation of the correspondence among S3T, QL, and NL dynamics. The first independently excites a set of zonal wavenumbers. This forcing was first used by Williams (1978) to parameterize excitation of barotropic dynamics by baroclinic instabilities. This forcing was also used by DelSole (2001) in his study of upper-level tropospheric jet dynamics and in the study of jet formation using S3T dynamics by Farrell and Ioannou (2003, 2007) and Bakas and Ioannou (2011). This stochastic forcing is spatially homogeneous but not isotropic and will be denoted as nonisotropic forcing (NIF). The second forcing is an isotropic narrow ring forcing (IRFn) concentrated near a single total wavenumber, *K _{f}*. This forcing structure has been used extensively in studies of beta-plane turbulence (cf. Vallis and Maltrud 1993) and was also used in the recent study of Srinivasan and Young (2012). It was introduced by Lilly (1969) in order to isolate the inverse cascade from the forcing in a study of two-dimensional turbulence.

^{1}The third forcing that we use is an isotropic ring forcing, in which the forcing is distributed over a wide annular region in wavenumber space around the central total wavenumber; it is denoted IRFw.

Specification of these stochastic forcing structures are given in appendix B. Plots of the corresponding power spectra together with instantaneous realizations both in vorticity and streamfunction for the three types of forcing structures are shown in Fig. 1. IRFn is peculiar in that it primarily excites vortices of scale 1/*K*_{f} that are evident in both the vorticity and streamfunction fields, while IRFw produces a streamfunction field dominated by large-scale structure similar to the fields excited by the other broadband forcings.

## 5. Stability in S3T of the homogeneous equilibrium state

The S3T system with homogeneous stochastic forcing and *ν*_{4} = 0 admits the homogeneous equilibrium solution:

which has no jets and eddy covariance proportional to the covariance of the forcing. The linear stability of perturbations to (**U**^{E}, ^{E}) is determined from the linearized equations:

with and . The temporal eigenvalues *σ* of these linear equations, which determine the stability of the equilibrium (**U**^{E}, ^{E}), satisfy Eqs. (C2) in appendix C. We obtain the stability of the homogeneous equilibrium under NIF and IRFn with *r* = 0.01 as a function of the parameter *ε*, which corresponds to the nondimensional rate of energy injection into the system (cf. appendix B). The growth rates, *σ*_{r} = Re(*σ*), as a function of the meridional wavenumber *n* of the zonal-mean flow perturbation and of *ε*, obtained using the method discussed in appendix C, are shown in Fig. 2. For both types of forcing, instability occurs for *ε* > *ε*_{c} and over a band of *n*. In all cases, the *σ* with greatest real part has zero imaginary part implying nontranslating jets have the largest growth rate. In the next section, predictions of S3T stability analysis for the bifurcation structure associated with jet formation will be compared with the corresponding QL and NL simulations. While QL and NL simulations reveal an apparent bifurcation, they cannot provide theoretical predictions of this bifurcation. We wish to examine the circumstances under which the underlying bifurcation structure predicted theoretically by the S3T stability analysis is reflected in the QL and NL simulations.

## 6. Bifurcations predicted by S3T and their reflection in QL and NL simulations

We examine the counterpart in NL and QL simulations of the S3T structural instability by comparing the evolution of the domain averaged energy of the zonal flow, . The amplitude of the zonal flow is measured, as in Srinivasan and Young (2012), with the zonal-mean-flow (zmf) index defined as , where *E*_{m} is the time-averaged energy of the zonal-mean flow, *E*_{m}(*t*), and *E*_{p} is the time average of the domain-averaged kinetic energy of the eddies, .

Zmf indices are shown as a function of *ε* in Fig. 3a for NIF forcing and in Fig. 3b for IRFn forcing, both for the case with *r* = 0.01 presented in Fig. 2. The fundamental qualitative prediction of S3T that jets form as a bifurcation in the strength of the turbulence forcing is verified in these plots. Agreement in the critical value *ε*_{c} for jet emergence is also obtained between S3T and QL while this parameter value is substantially larger in NL. For example, jets emerge in the NL simulations at under NIF forcing and at under IRFn. Similar behavior was noted by Srinivasan and Young (2012). The reason for this difference will be explained in section 7.

S3T dynamics not only predicts the emergence of zonal jets as a bifurcation in turbulence forcing, but also predicts the structure of the finite amplitude jets that result from equilibration of the initial jet formation instability. These finite-amplitude jets correspond to fixed points of the S3T dynamics. An example for IRFn strongly forced with *ε* = 100*ε*_{c} and with damping *r* = 0.01 is shown in Fig. 4. This example demonstrates the essential similarity among the jets in NL, QL, and S3T simulations.

Under strong turbulence forcing, the initial S3T jet formation instability typically reaches final equilibrium as a finite-amplitude jet at a wavenumber smaller than that of the initial instability. An example is the case of IRFn at *ε* = 100*ε*_{c} shown in Fig. 4. In this example, the jets emerge in S3T initially with *n* = 10, which is in agreement with the prediction of the S3T instability of the homogeneous equilibrium, but eventually equilibrate at *n* = 3 following a series of jet mergers, as seen in the Hovmöller diagram. Similar dynamics are evident in the NL and QL simulations. This behavior can be rationalized by noting that if the wavenumber of the jet remains fixed, then as jet amplitude continues to increase under strong turbulence forcing, violation of the Rayleigh–Kuo stability criterion would necessarily occur. By transitioning to a lower wavenumber, the flow is able to forestall this occurrence of inflectional instability. However, detailed analysis of the S3T stability of the finite-amplitude equilibria near the point of jet merger reveals that these mergers coincide with the inception of a structural instability associated with eddy–mean flow interaction, which precedes the occurrence of hydrodynamic instability of the jet (Farrell and Ioannou 2003, 2007).^{2}

## 7. Influence of the turbulence spectrum on the S3T jet formation instability

Both QL and S3T dynamics exclude interactions among eddies and include only the nonlocal interactions between jets, with *k* = 0, and eddies, with *k* ≠ 0. Therefore, there is no enstrophy or energy cascade in wavenumber space in either QL or S3T dynamics and the homogeneous S3T equilibrium state [Eq. (6)] has spectrum , which is determined by the spectrum of the forcing ( is the spectral power of the forcing covariance; cf. appendix B). However, this is not true in NL, which includes eddy–eddy interactions producing enstrophy/energy cascades. For example, in NL an isotropic ring forcing is spread as time progresses, becoming concentrated at lower wavenumbers and forming the characteristic dumbbell shape seen in beta-plane turbulence simulations (cf. Vallis and Maltrud 1993), and consequently the homogeneous turbulent state is no longer characterized by the spectrum of the forcing. We can take account of this modification of the spectrum by performing S3T stability on the homogeneous state under the equivalent forcing covariance,

which maintains the observed NL spectrum in the S3T dynamics. The NL modified eddy vorticity spectrum, , is obtained from an ensemble of NL simulations. Plots of under IRFn are shown in Figs. 6a–c for various *ε* and *r* values. The departure of the NL spectra from the spectra of the QL and S3T equilibria is evident and this departure depends on the amplitude of *ε* and *r*.

We now demonstrate that while the fundamental qualitative prediction of S3T that jets form as a bifurcation in turbulence forcing and in the absence of turbulent cascades is verified in both QL and NL, a necessary condition for obtaining quantitative agreement between NL and both S3T and QL dynamics is that the equilibrium spectrum used in the S3T and QL dynamics be close to the equilibrium spectrum obtained in NL so that the stability analysis is performed on similar states. In the case with IRFn and *r* = 0.01, formation of persistent finite-amplitude zonal jets occurs in the NL simulations at (cf. Fig. 3b). In agreement, S3T stability analysis on the NL modified equilibrium IRFn spectrum (denoted S3Tb and shown in Fig. 6c) predicts instability for *ε* ≥ 2.8*ε*_{c} (cf. Fig. 6e). Moreover, S3T stability analysis with the S3Tb spectrum predicts jet formation at *n* = 6 and in agreement with this prediction, jets emerge in NL with *n* = 6. Hovmöller diagrams demonstrating similar jet evolution in NL under IRFn and in S3T under S3Tb forcing are shown in Fig. 7. We also note that agreement between NL and S3T in predictions of jet amplitude at large supercriticality is also obtained by using the S3Tb spectrum (cf. Fig. 3).^{3}

This influence of the eddy spectrum on jet dynamics is revealed in the case of IRFn at *ε* = 2*ε*_{c}, which is shown in Fig. 8. Although, at this energy input rate, S3T under IRFn is structurally unstable, no jets emerge in NL. We have shown that agreement in bifurcation structure is obtained between NL and S3T when S3T analysis is performed with the S3Tb spectrum. We now examine the development of the NL spectrum toward S3Tb and demonstrate the close control exerted by this evolving spectrum on S3T stability. The evolving spectrum, shown in Figs. 9a–f, is obtained using an ensemble of NL simulations, each starting from a state of rest and evolving under a different forcing realization. A sequence of S3T analyses performed on this evolving ensemble spectrum is shown in Fig. 9g. The weak NL ensemble spectrum at *t* = 1 does not support instability, but by *t* = 20, the ensemble spectrum, having assumed the isotropic ring structure of the forcing, becomes S3T unstable. This structural instability results in the formation of an incipient *n* = 6 jet structure, which is evident by *t* = 50 in the NL simulation shown in Fig. 8. As the spectrum further evolves, the S3T growth rates decrease and no jet structure is unstable for *t* > 120, and decay rates continue to increase until *t* = 250 (cf. Fig. 9g). This example demonstrates the tight control on S3T stability exerted by the spectrum. Furthermore, it shows the close association between S3T instability and the emergence of jet structure in NL.

## 8. Influence of nonzonal structures predicted by S3T on the turbulence spectrum and on jet dynamics

Despite S3T supercriticality, no persistent jets emerge in NL simulations with IRFn in the interval *ε*_{c} < *ε* < 2.8*ε*_{c} (cf. Fig. 3b). Comparisons of NL, QL, and S3T simulations with IRFn forcing at *ε* = 2*ε*_{c} are shown in Fig. 8. Instead of zonal jets, in the NL simulation, prominent nonzonal structures are seen to propagate westward at the Rossby wave phase speed. These nonzonal structures are also evident in the concentration of power in the enstrophy spectrum at (|*k*|, |ℓ|) = (1, 7) (cf. Figs. 10a and 10b). At this forcing amplitude, these structures are essentially linear Rossby waves that, if stochastically forced, would be coherent only over the dissipation time scale 1/*r*. Coherence on the dissipation time scale is observed in the subdominant part of the spectrum as seen in the case of the (3, 6) structure in Fig. 11c. However, the dominant (1, 7) structure remains coherent over time periods far exceeding the dissipation time scale (cf. Hovmöller diagram in Fig. 11b). This case represents a regime in which the flow is dominated by a single nonzonal structure. Both the concentration of power in and the coherence of this structure will be addressed below.

When the forcing is increased to *ε* = 10*ε*_{c}, a (0, 6) jet structure emerges, suppresses the nonzonal (1, 7) structure, and becomes the dominant structure. A prominent phase-coherent nonzonal (1, 5) structure propagating with the Rossby wave speed is also present, as shown in Fig. 12a. A similar regime of coexisting jets and nonzonal structures is also evident at higher supercriticalities. An example is the case of the equilibrium state at *ε* = 100*ε*_{c} (cf. Fig. 4) in which the energy of the flow is shared between the (0, 3) jet and the (1, 3) structure, as shown in Fig. 12b. At this forcing level, the (1, 3) structure is not phase coherent, but its phase speed is still given by the Rossby wave speed. At even higher forcing, similar nonzonal structures, referred to as zonons, have been reported to coexist with zonal jets while propagating phase incoherently at speeds that differ substantially from the Rossby wave speed (Sukoriansky et al. 2008). These cases provide examples of the regime in which jets and nonzonal structures coexist.

To study the dynamics of nonzonal structures within the framework of S3T, a different interpretation of the ensemble mean in the S3T formulation is required: instead of interpreting the ensemble means as zonal means, interpret them rather as Reynolds averages over an intermediate time scale (Bernstein 2009; Bernstein and Farrell 2010; Bakas and Ioannou 2013a, 2014). Analysis of S3T stability of the homogeneous equilibrium state using this broader interpretation (cf. appendix D) reveals that when the energy input rate reaches the value *ε*_{c}, which is the S3T stability threshold for the emergence of zonal jets, the state may already be unstable to nonzonal structures. This can be seen in the stability analysis shown in Fig. 13, which reveals that the maximum growth rate occurs at wavenumbers corresponding to nonzonal structures. In agreement with this stability analysis, the spectrum of the NL simulation shows concentration of power in these most S3T unstable wavenumbers (cf. Fig. 10).

The dominance and persistence of the structures seen in these NL simulations can be understood from this stability analysis and its extension into the nonlinear regime. Because the stochastic forcing is white in time, the energy injection rate is fixed and state independent and, assuming linear damping at rate *r* dominates the dissipation, the total flow energy assumes the fixed and state-independent mean value *E*_{m} + *E*_{p} = *ε*/(2*r*). At finite amplitude, the set of S3T unstable structures equilibrate to allocate among themselves most of this energy, which results in the dominance of a small subset of these structures. However, we find that in this competition, a specific zonal jet structure has primacy so that even if this structure is not the most linearly unstable, it emerges as the dominant structure.

An attractive means for exploring the dynamics of the interaction between jets and nonzonal structures is changing the jet damping rate in Eq. (2a) from *r* to *r*_{m} and allowing it to assume values different from *r* in Eq. (2b). In this way, we can control the relative stability of jets and nonzonal structures, as well as the finite equilibrium amplitude reached by the jet. This asymmetric damping may be regarded as a model for approximating jet dynamics in a baroclinic flow in which the upper-level jet is lightly damped, while the active baroclinic turbulence generating scales are strongly Ekman damped. This asymmetry in the damping between upper and lower levels contributes to making jets in baroclinic turbulence generally stronger than jets in barotropic turbulence (Farrell and Ioannou 2007, 2008). By appropriate choice of *r* and *r*_{m}, a regime can be obtained in which the zonal jet instability appears first as *ε* increases. Because once jets are unstable they dominate nonzonal structures, in this regime zonal jets are the dominant coherent structure and S3T analysis based on the zonal interpretation of the ensemble mean produces very good agreement with NL. For example, a comparison of bifurcation structures among S3T, QL, and NL under NIF and IRFn using the asymmetric damping *r* = 0.1 and *r*_{m} = 0.01 demonstrates that jets emerge at the same critical value in S3T, QL, and NL (cf. Figs. 14a and 14b). This agreement, which has been obtained by asymmetric damping, induced suppression of the nonzonal instability up to *ε*_{c}, implies that in the simulations with symmetric damping the difference between the S3T prediction for *ε* required for the first emergence of jets and the value of *ε* obtained for first jet appearance in NL (cf. Fig. 3) can be attributed to modification of the background spectrum by the prior emergence of the nonzonal structures predicted by S3T. Once unstable, zonal structures immediately dominate nonzonal structures, which explains why S3T dynamics based on the zonal-mean interpretation of the ensemble mean produces accurate results for parameter values for which zonal jets are the first instability to occur.

A comparison of the development of jets in S3T, QL, and NL with this asymmetric damping and NIF forcing, shown in Fig. 15, demonstrates the accuracy of the S3T predictions. S3T analysis predicts that in this case with NIF forcing maximum instability occurs at *n* = 6. When these maximally growing eigenfunctions are introduced in the S3T system, the jets grow exponentially at first at the predicted rate and then equilibrate. Corresponding simulations with the QL and NL dynamics reveal nearly identical jet growth followed by finite-amplitude equilibration (shown in Fig. 15). Similar results are obtained with IRFn. This demonstrates that the S3T dynamics comprises both the jet instability mechanism and the mechanism of finite-amplitude equilibration.

Although no theoretical prediction of this bifurcation behavior can be made directly from NL or QL, they both reveal the bifurcation structure obtained from the S3T analysis. By suppressing the peripheral complexity of nonzonal structure formation by nonzonal S3T instabilities, these simulations allow construction of a simple model example that provides compelling evidence for identifying jet formation and equilibration in NL with the S3T theoretical framework. Moreover, agreement among the NL, QL, and S3T bifurcation diagrams shown in Figs. 14a and 14b provides convincing evidence that turbulent cascades, which are absent in S3T or QL, are not required for jet formation.

While under NIF agreement between NL and S3T, equilibrium jet amplitudes extends to all values of *ε*, under IRFn, the NL and S3T equilibrium amplitudes diverge at larger values of *ε* (cf. Figs. 14a and 14b). This difference between NL and QL-S3T at large *ε* cannot be attributed to nonlinear modification of the spectrum, which is accounted for by use of the S3Tb spectrum (cf. S3Tb response in Fig. 14b). Rather, this difference is primarily due to nonlinear eddy–eddy interactions retained in NL that disrupt the upgradient momentum transfer. This disruption is accentuated by the peculiar efficiency with which IRFn gives rise to vortices, as can be seen in Figs. 1d–f. The more physical distributed forcing structures do not share this property (cf. Fig. 1). We verify that IRFn is responsible for depressing NL equilibrium jet strength at high supercriticality by broadening the forcing distribution to assume the form IRFw (cf. appendix B, as well as Fig. 1 for IRFn–IRFw comparison). Using IRFw while retaining other parameters as in Fig. 14b, we obtain agreement between S3T, QL, and NL simulations, as is shown in Fig. 14c.

## 9. Identification of intermittent jets with stable S3T zonal eigenfunctions

For subcritical forcing, S3T predicts a stable homogeneous statistical equilibrium and a set of eigenfunctions that govern the decay of perturbations to this equilibrium. We wish to show that these eigenfunctions are excited in NL by fluctuations in the turbulence and that this excitation gives rise in NL simulations to the formation of intermittent jets with the form of these eigenfunctions.

As an example, consider the simulation with asymmetric damping and IRFn subcritical forcing shown in Fig. 16. For these parameters, the least damped eigenfunctions are zonal jets, and confirmation that the intermittent jets in NL, shown in Fig. 16a, are consistent with turbulence fluctuations exciting the S3T damped modes is given in Fig. 16c, where the intermittent jets resulting from stochastic forcing of the S3T modes themselves are shown. This diagram was obtained by plotting , with *α*_{n} independent red noise processes, associated with the damping rates |*σ*(*n*)| of the first *N* = 15 least damped S3T modes. These *α*_{n} are obtained from the Langevin equation, *dα*_{n}/*dt* = *σ*(*n*)*α*_{n} + *ξ*(*t*), where *ξ*(*t*) is a *δ*-correlated complex valued random variable.

The fluctuation-free S3T simulations reveal a persistent jet structure only coincident with the inception of the S3T instability, which occurs only for supercritical forcing. However, in QL and NL simulations, fluctuations excite the damped manifold of modes predicted by the S3T analysis to exist at subcritical forcing amplitudes. This observation confirms the reality of the manifold of S3T stable modes.

In NL and QL simulations, these stable modes predicted by S3T are increasingly excited as the critical bifurcation point in parameter space is approached, because their damping rate vanishes at the bifurcation. The associated increase in zonal-mean-flow energy on approach to the bifurcation point obscures the exact location of the bifurcation point in NL and QL simulations compared to the fluctuation-free S3T simulations for which the bifurcation is exactly coincident with the inception of the S3T instability (i.e., Fig. 14).

## 10. Verification in NL of the multiple jet equilibria predicted by S3T

As is commonly found in nonlinear systems, the finite-amplitude equilibria predicted by S3T are not necessarily unique and multiple equilibria can occur for the same parameters. S3T provides a theoretical framework for studying these multiple equilibria, their stability, and bifurcation structure. An example of two such S3T equilibria is shown in Fig. 17 together with their associated NL simulations. As the parameters change, these equilibria may cease to exist or become S3T unstable. Similar multiple equilibria have been found in S3T studies of barotropic beta-plane turbulence (Farrell and Ioannou 2003, 2007; Parker and Krommes 2014) and in S3T studies of baroclinic turbulence (Farrell and Ioannou 2008, 2009c), and the hypothesis has been advanced that the existence of such multiple jet equilibria may underlie the abrupt transitions found in the record of Earth’s climate (Farrell and Ioannou 2003; Wunsch 2003).

## 11. Conclusions

In this work, predictions of S3T for jet formation and equilibration in barotropic beta-plane turbulence were critically compared with results obtained using QL and NL simulations. The qualitative bifurcation structure predicted by S3T for emergence of zonal jets from a homogeneous turbulent state was confirmed by both the QL and NL simulations. Moreover, the finite-amplitude equilibrium jets in NL and QL simulations were found to be as predicted by the fixed-point solutions of S3T. Differences in jet formation bifurcation parameter values between NL and QL–S3T were reconciled by taking account of the fact that the spectrum of turbulence is substantially modified in NL. Remarkably, the modification of the spectrum in NL could be traced in large part to the emergence of nonzonal structures through S3T instability. When account is taken of the modification of the turbulent spectrum resulting substantially from these nonzonal structures, S3T also provides quantitative agreement with the threshold values for the emergence of jets in NL. The influence of the background eddy spectrum on the S3T dynamics was found to be immediate, in the sense that in spinup simulations, jets emerge in accordance with the instability calculated on the temporally developing spectrum. The fact that jets are prominent in observations is consistent with the robust result that when a jet structure emerges, it has primacy over the nonzonal structures, so that even if the jet eigenfunction is not the most linearly S3T unstable eigenfunction, the jet still emerges at finite amplitude as the dominant structure.

These results confirm that jet emergence and equilibration in barotropic beta-plane turbulence results from the cooperative quasi-linear mean flow–eddy instability that is predicted by S3T. These results also establish that turbulent cascades are not required for the formation of zonal jets in beta-plane turbulence. Moreover, the physical reality of the manifold of stable modes arising from cooperative interaction between incoherent turbulence and coherent jets, which is predicted by S3T, was verified in this work by relating observations of intermittent jets in NL and QL to stochastic excitation by the turbulence of this manifold of stable S3T modes.

S3T provides an autonomous, deterministic nonlinear closure of turbulence dynamics at second order that provides an attractive vehicle for further investigation of the dynamics of turbulent flows.

## Acknowledgments

The authors acknowledge discussions with N. Bakas, F. Bouchet, K. Srinivasan, and W. Young. Navid Constantinou acknowledges the support of the Alexander S. Onassis Public Benefit Foundation. Brian Farrell was supported by NSF AGS-1246929 and ATM-0736022. Brian Farrell and Petros Ioannou acknowledge the hospitality during June 2012 of the Aspen Center for Physics (supported by NSF under Grant 1066293) where part of this paper was written. Petros Ioannou acknowledges the generous support of the John S. Latsis Foundation under “Research Projects 2011.”

### APPENDIX A

#### Numerical Details and Parameters

Both the nonlinear (NL) simulations of Eq. (1) and the quasi-linear (QL) simulations of Eqs. (3) were carried out with a pseudospectral Fourier code. The maximum resolved wavenumbers were *k*_{max} = *N*_{x}/2 and ℓ_{max} = *N*_{y}/2 and the maximum resolved total wavenumber was . For the time integration, a fourth-order Runge–Kutta method (RK4) was used together with a Godunov step for integrating the stochastic forcing. In all calculations, hyperviscosity was added for numerical stability with coefficient , where *δt* is the time step. In all calculations, *δt* = 2.5 × 10^{−3} and *N*_{x} = *N*_{y} = 256, which imply *ν*_{4} = 1.86 × 10^{−7}.

### APPENDIX B

#### The Stochastic Forcing Structure and Its Associated Power Spectrum

The stochastic forcing at point (*x*_{i}, *y*_{j}) is

in which *ξ*_{kp} values are temporally *δ* correlated and independent and satisfy 〈*ξ*_{kp}(*t*)〉 = 0, . The stochastic forcing is correlated in *y* by the columns of the matrix . For the nonisotropic forcing (NIF), this meridional structure, in a periodic domain with period 2*π* in *y*, is specified by

We choose and force zonal components *k* = 1, … , 14 . Because the stochastic forcing is *δ* correlated in time, the energy input rate, given by

does not depend on the state of the system and can be independently specified. The normalization constant *c*_{k} in Eq. (B2) is chosen so that each *k* is excited equally and one unit of energy is injected in total. It follows that the total energy input rate in the NL, QL, or S3T simulations is given by *ε*.

The isotropic ring forcing is specified by

with ℓ(*p*) = (*p* − 1) − *N*_{y}/2 and . Meridional wavenumber ℓ extends from −*N*_{y}/2 to *N*_{y}/2 − 1 because only these wavenumbers are resolved when *N*_{y} points in the meridional direction are retained (for *N*_{y} even). For IRFn, *w*_{kp} = 1 for |*K* − *K*_{f}| ≤ *δk*_{f} and *w*_{kp} = 0 otherwise, with *K*_{f} = 14 and *δk*_{f} = 1. For IRFw, with *K*_{f} = 14 and . The normalization constant *c* is chosen for both cases so that the total energy input rate is unity. IRFn and IRFw are both spatially homogeneous and nearly isotropic in a finite doubly periodic domain. They approach exact isotropy as the domain size increases.

The spatial covariance of the forcing, *Q*(*x*_{a}, *x*_{b}, *y*_{a}, *y*_{b}) = 〈*F*(*x*_{a}, *y*_{a}, *t*) *F*(*x*_{b}, *y*_{b}, *t*)〉, being homogeneous in both *x* and *y*, depends only on *x*_{a} − *x*_{b} and *y*_{a} − *y*_{b} and has Fourier expansion

The values are the spatial power spectrum of the stochastic forcing. Fourier coefficients of the forcing covariance for only positive values of zonal wavenumbers *k* [cf. Eq. (B5)], can be related to Fourier expansions in both positive and negative zonal wavenumbers,

through and for *k* > 0. In the derivation of these relations, the symmetry of the forcing covariance under exchange of the two points is used.

### APPENDIX C

#### Determining the S3T Stability of the Homogeneous State

Equation (7) determines the S3T stability of the equilibrium state. Because of the presence of the imaginary part in Eq. (7a), in order to proceed with eigenanalysis of this system, we need to treat the real and imaginary parts of the covariances as independent variables. Writing the covariances as *δ*_{k} = *δ*_{k,}_{R} + *iδ*_{k,}_{I} and , and the operators as and , we obtain the real coefficient system:

in which the subscript *k* in all the variables in Eqs. (C1b) and (C1c) has been omitted. In Eqs. (C1), the coefficient of linear damping of the mean flow *r*_{m} may differ from the coefficient of linear damping of the nonzonal perturbations *r* (cf. section 8). The asymptotic stability of Eqs. (C1) is determined by assuming solutions of the form , with and by determining the eigenvalues *σ* and the eigenfunctions of the system:

In most cases, direct eigenanalysis of this system is computationally prohibitive because it involves eigenanalysis of matrices of dimension if *N*_{y} grid points are used to approximate the functions and *N*_{k} zonal wavenumbers are forced. In this section, we describe an efficient iterative method that can produce solutions to this stability problem for large *N*_{y}. The method is a generalization of the adiabatic approximation used in earlier studies (Farrell and Ioannou 2003, 2007; Bakas and Ioannou 2011).

When Eqs. (C2) have eigenvalues with Re(*σ*) > 0, the equilibrium is S3T unstable. When *σ* is complex the eigenfunctions , , , and will be complex. Realizable, Hermitian solutions can then be formed by superposing the complex conjugate eigenfunction. Note that the covariances are required to be Hermitian but need not be positive definite [for a discussion of eigenvalue problems involving covariances, cf. Farrell and Ioannou (2002)].

Because of the periodic boundary conditions, the mean-flow eigenfunctions are, in general, a superposition of harmonics. However, here we are treating the stability of the homogeneous equilibrium and the eigenfunctions can be shown to be single harmonics, , and Eq. (C2a) becomes

The number of unstable jets, if the equilibrium is unstable, is *n*. Equation (C3) can be regarded as an equation for *σ* given that satisfies the coupled Eqs. (C2b) and (C2c) and is therefore a function of *σ* and *n*. Having transformed Eq. (C3) into an equation for *σ* for a given *n*, the eigenvalues can be determined by iteration.

It is advantageous to solve Eqs. (C2b) and (C2c) for the variables , , which satisfy the decoupled Sylvester equations:

(Note that if the eigenvectors and are complex, despite the notation, then and .)

### APPENDIX D

#### Determining the S3T Stability of the Homogeneous State to Nonzonal Perturbations

The S3T stability of a homogeneous background state to nonzonal perturbations of the form *e*^{σt}*e*^{i(mx+ny)}, in which *m* is the zonal wavenumber of the perturbations and *n* is the meridional wavenumber of perturbations, is determined by the eigenvalues *σ* that solve for each (*m*, *n*) the equation

In the above equation, *N*^{2} = *m*^{2} + *n*^{2}, *K*^{2} = *k*^{2} + ℓ^{2}, , *k*_{+} = *k* + *m*/2, ℓ_{+} = ℓ + *n*/2, *δ*_{ij} is Kronecker’s delta, and values are the Fourier coefficients of the forcing covariance as defined in Eq. (B6). This form of the equation is appropriate for a square domain of length 2*π* and the summations are over the integer wavenumbers *k* and ℓ. The derivation of the above equation can be found in Bakas and Ioannou (2014). For a specified forcing with spectrum , the growth rates are obtained using Newton’s method.

At high supercriticality (i.e., as *ε* → **∞**), the maximal growth rate *σ*_{r} of the (*m*, *n*) large-scale structure asymptotically approaches

with

and the frequency of this eigenstructure assumes asymptotically the Rossby wave frequency, *σ*_{i} = Im(*σ*) = *mβ*/(*m*^{2} + *n*^{2}). This asymptotic expression for the growth rate and phase speed of the large-scale structure is useful for tracing the maximal growth rates as a function of supercriticality using Newton’s iterations.

The asymptotic growth rates depend only on the forcing distribution. The growth rates for the NIF and IRFn used in this paper are shown in Fig. D1. It can be shown that the asymptotic growth rate vanishes for exactly isotropic forcing. Asymptotically, the growth rates do not depend on the damping rate of the mean flow *r*_{m}. NIF favors at least initially jet formation, while IRFn favors formation of nonzonal structures.

## REFERENCES

*β*-plane turbulence.

*Phys. Rev. Lett.,*

**110,**224501, doi:10.1103/PhysRevLett.110.224501.

*J. Fluid Mech.,*

**740,**312–341, doi:10.1017/jfm.2013.663.

*Phys. Plasmas,*

**16,**112903, doi:10.1063/1.3258666.

*Jupiter: the Planet, Satellites, and Magnetosphere,*F. Bagenal, T. E. Dowling, and W. B. McKinnon, Eds., Cambridge University Press, 105–128.

*Geophys. Res. Lett.,*

**34,**L22801, doi:10.1029/2007GL031779.

*Phys. Rev. Lett.,*

**101,**178501, doi:10.1103/PhysRevLett.101.178501.

*Phys. Rev. Lett.,*

**110,**104502, doi:10.1103/PhysRevLett.110.104502.

*Phys. Fluids,*

**20,**126602, doi:10.1063/1.2990042.

## Footnotes

^{1}

In the limit that only wavenumber *K _{f}* is excited by IRFn, this forcing would produce no cascade because any sum of Rossby waves, each of which has the same total wavenumber, is a nonlinear solution. This observation establishes that the jet formation mechanism for IRFn is necessarily that of S3T because nonlinearity other than the QL interaction with the jet, which is retained in S3T, vanishes. However, in the presence of perturbations with other wavenumbers, a cascade ensues leading eventually to an equilibrium spectrum under IRFn forcing. The dynamics and consequences for the jet formation process of the establishment of this spectrum are studied in section 7 (cf. Fig. 9).

^{2}

Jet mergers occur in the Ginzburg–Landau equations that govern the dynamics of the S3T instability of the homogeneous equilibrium state for parameter values for which the system is close to marginal stability (Parker and Krommes 2014). However, these mergers in the Ginzburg–Landau equations are associated with equilibration of the Eckhaus instability rather than equilibration of the inflectional instability associated with violation of the Rayleigh–Kuo criterion, as is the case for mergers of finite-amplitude jets (cf. Fig. 5). Characteristic of this difference is that in the case of the Ginzburg–Landau equations, both the prograde and retrograde jets merge, while in the case of the finite-amplitude jets, only the prograde jets merge. The same phenomenology as in the Ginzburg–Landau equations occurs in the case of the Cahn–Hilliard equations that govern the dynamics of marginally stable jets in the modulational instability study of Manfroi and Young (1999).

^{3}

The spectral peaks near the axis do not directly influence the stability of the NL modified spectrum, which is determined by the distorted and broadened ring spectrum. However, while the spectral peaks do not influence the stability directly, they do influence it indirectly by distorting the incoherent spectrum.