Emergence of non-zonal coherent structures Nikolaos A. Bakas Department of Physics, University of Ioannina, Ioannina, Greece Petros J. Ioannou∗ Department of Physics, National and Kapodistrian University of Athens, Athens, Greece 5 1 0 Planetary turbulence is observed to self-organize into large-scale structures such 2 as zonal jets and coherent vortices. One of the simplest models that retains the n relevant dynamics of turbulent self-organization is a barotropic flow in a beta-plane a channel with turbulence sustained by random stirring. Non-linear integrations J of this model show that as the energy input rate of the forcing is increased, the 1 homogeneity of the flow is first broken by the emergence of non-zonal, coherent, 2 westward propagating structures and at larger energy input rates by the emergence ] of zonal jets. The emergence of both non-zonal coherent structures and zonal jets h is studied using a statistical theory, Stochastic Structural Stability Theory (S3T). p S3T directly models a second-order approximation to the statistical mean turbulent - o state and allows the identification of statistical turbulent equilibria and study of a their stability. Using S3T, the bifurcation properties of the homogeneous state in . s barotropic beta-plane turbulence are determined. Analytic expressions for the zonal c andnon-zonallarge-scalecoherentflowsthatemergeasaresultofstructuralinstability i s are obtained and the equilibration of the incipient instabilities is studied through y numerical integrations of the S3T dynamical system. The dynamics underlying h p the formation of zonal jets are also investigated. It is shown that zonal jets form [ from the upgradient momentum fluxes that result from the shearing of the eddies by the emerging infinitesimal large-scale flow. Finally, numerical simulations of the 1 v nonlinear equations confirm the characteristics (scale, amplitude and phase speed) of 0 the structures predicted by S3T, even in highly non-linear parameter regimes such as 8 the regime of zonostrophic turbulence. 2 5 0 . 1 0 Atmospheric and oceanic turbulence is commonly ob- et al. 2014) and non-zonal coherent structures (Bakas and 5 served to be organized into spatially and temporally co- Ioannou 2014). This argues that the emergence of coherent 1 herent structures such as zonal jets and coherent vortices. structures in a homogeneous background of turbulence is a : v A simple model that retains the relevant dynamics, is a bifurcation phenomenon. i barotropic flow on a β-plane with turbulence sustained by An advantageous method to study such a phenomenon, X random stirring. Numerical simulations of the stochasti- is to adopt the perspective of statistical state dynamics of r a callyforcedbarotropicvorticityequationonthesurfaceofa the flow, rather than look into the dynamics of sample real- rotatingsphereoronaβ-plane,haveshownthecoexistence izations of direct numerical simulations. This amounts to ofrobustzonaljetsandoflarge-scalewestwardpropagating studythedynamicsandstabilityofthestatisticalequilibria coherent structures that are referred to as satellite modes arising in the turbulent flow, which are fixed points of the (Danilov and Gurarie 2004) or zonons (Sukariansky et al. equationsgoverningtheevolutionoftheflowstatistics. This 2008; Galperin et al. 2010). Emergence of these coherent approach is followed in the Stochastic Structural Stability structuresinbarotropicturbulencehasalsoanotherfeature. Theory (S3T) (Farrell and Ioannou 2003) or Second Order As the energy input of the stochastic forcing is increased or Cumulant Expansion theory (CE2) (Marston et al. 2008). dissipation is decreased, there is a sudden onset of coher- This theory is based on two building blocks. The first is ent zonal flows (Srinivasan and Young 2012; Constantinou to do a Reynolds decomposition of the dynamical variables 1 into the sum of a mean value that represents the coherent this assumption. The second goal is to study in detail the flow and fluctuations that represent the turbulent eddies eddy-mean flow dynamics underlying the S3T instability andthenformthecumulantscontainingtheinformationon focusing on the example of jet formation. And the third the mean values (first cumulant) and on the eddy statistics goal is to compare the characteristics of the structures that (higher order cumulants). The second building block is emergeinS3Tagainstnon-linearsimulations,eveninhighly to truncate the equations governing the evolution of the non-linear regimes that at first glance present a challenging cumulants at second order by either parameterizing the test for the restricted dynamics of S3T. terms involving the third cumulant (Farrell and Ioannou 1993a,b,c;DelSoleandFarrell1996;DelSole2004)orsetting 1. Formulation of Stochastic Structural Stability the third cumulant to zero (Marston et al. 2008; Tobias Theory under a generalized average et al. 2011; Srinivasan and Young 2012). Restriction of Consider a nondivergent barotropic flow on a β-plane the dynamics to the first two cumulants is equivalent to with cartesian coordinates x = (x,y). The velocity field, neglectingtheeddy-eddyinteractionsinthefullynon-linear u=(u,v), is given by (u,v)=(−∂ ψ,∂ ψ), where ψ is the dynamics and retaining only the interaction between the y x streamfunction. Relative vorticity ζ(x,y,t)=∆ψ, evolves eddies with the instantaneous mean flow. While such a according to the non-linear (NL) equation: second order closure might seem crude at first sight, there √ is strong evidence to support it (Bouchet et al. 2013). (∂ +u·∇)ζ+βv =−rζ−ν∆2ζ+ εfe, (1) t Previous studies employing S3T have already addressed the bifurcation from a homogeneous turbulent regime to where ∆=∂2 +∂2 is the horizontal Laplacian, β is the xx yy a jet forming regime in barotropic β-plane turbulence and gradient of planetary vorticity, r is the coefficient of linear identifiedtheemergingjetstructuresbothnumerically(Far- dissipationthattypicallyparameterizesEkmandraginplan- rellandIoannou2007)andanalytically(BakasandIoannou etaryatmospheresandν isthecoefficientofhyper-diffusion 2011;SrinivasanandYoung2012)aslinearlyunstablemodes that dissipates enstrophy flowing into unresolved scales. tothehomogeneousturbulentstateequilibrium. Itwasalso The exogenous forcing term fe, parameterizes processes shownthattheresultingdynamicalsystemfortheevolution such as small scale convection or baroclinic instability, that of the first two cumulants linearized around the homoge- are missing from the barotropic dynamics and is necessary neous equilibrium possesses the mathematical structure of to sustain turbulence. We assume that fe is a temporally the dynamical system of pattern formation (Parker and delta correlated and spatially homogeneous random stir- Krommes 2013). Comparison of the results of the stability ring injecting energy at a rate ε and having a two-point, analysis with direct numerical simulations have shown that two-time correlation function of the form: the structure of zonal flows that emerge in the non-linear simulationscanbepredictedbyS3T(SrinivasanandYoung (cid:104)fe(x ,y ,t )fe(x ,y ,t )(cid:105)=δ(t −t )Ξ(x ,x ,y ,y ), 1 1 1 2 2 2 2 1 1 2 1 2 2012; Constantinou et al. 2014). However, these research (2) efforts, have assumed that the ensemble average employed where the brackets denote an ensemble average over the inS3Tisequivalenttoazonalaverage,asimplificationthat different realizations of the forcing. treats the non-zonal structures as incoherent and cannot S3T describes the statistical dynamics of the first two address their emergence and effect on the jet dynamics. In same time cumulants of (1). The equations governing the addition, the eddy-mean flow dynamics underlying the S3T evolution of the first two cumulants are obtained as follows. instability even in the jet formation case, that involve only We decompose the vorticity field into the averaged field, the interactions of small scale waves with the large-scale Z =T[ζ], defined as a time average over an intermediate coherent structures are not clear. time scale and deviations from the mean or eddies, ζ(cid:48) = So the goals in this article are the following. The first ζ−Z. The intermediate time scale is larger than the time goal is to adopt a more general interpretation of the ensem- scale of the turbulent motions but smaller than the time ble average, in order to address the emergence of coherent scale of the large scale motions. With this decomposition non-zonal structures. We adopt the more general inter- we rewrite (1) as: pretation that the ensemble average is a Reynolds average over the fast turbulent motions (Bernstein 2009; Bernstein (∂t+U·∇)Z+βV =−∇·T[u(cid:48)ζ(cid:48)]−rZ−ν∆2Z, (3) and Farrell 2010). With this definition of the ensemble where U = [U,V] = [−∂ Ψ,∂ Ψ] and u(cid:48) = [u(cid:48),v(cid:48)] = mean, we obtain the statistical dynamics of the interaction y x [−∂ ψ(cid:48),∂ ψ(cid:48)] are the mean and the eddy velocity fields of the coarse-grained ensemble average field, which can y x respectively. The mean vorticity is therefore forced by the be zonal or non-zonal coherent structures represented by divergence of the mean vorticity fluxes. The eddy vorticity their vorticity, with the fine-grained incoherent field represented by the vorticity second cumulant and we revisit the structural stability of the homogeneous equilibrium under 2 ζ(cid:48) evolves according to: where (∂t+U·∇)ζ(cid:48)+(β+∂yZ)v(cid:48)+u(cid:48)∂xZ = A =−U ·∇ −(β+∂ Z)∂ ∆−1+∂ Z∂ ∆−1−r−ν∆2, i i i yi xi i xi yi i i =−rζ(cid:48)−ν∆2ζ(cid:48)+fe+T[u(cid:48)·∇ζ(cid:48)]−u(cid:48)·∇ζ(cid:48), (4) (11) (cid:124) (cid:123)(cid:122) (cid:125) governs the dynamics of linear perturbations about the in- fnl stantaneousmeanflowU. Therighthandsideof(10)isthe where fnl is the term involving the non-linear interactions correlation of the external forcing with vorticity, which for among the turbulent eddies. A closed system for the sta- deltacorrelatedstochasticforcingisindependentofthestate tistical state dynamics is obtained by first neglecting the oftheflowandisequalatalltimestotheprescribedforcing √ eddy-eddy term fnl to obtain the quasi-linear system, covariance: ε(cid:104)feζ(cid:48) +feζ(cid:48)(cid:105) = ε(cid:104)fefe(cid:105) = εΞ. Therefore 1 2 2 1 1 2 The second cumulant evolves then according to: (∂ +U·∇)Z+βV =−∇·T[u(cid:48)ζ(cid:48)]−rZ−ν∆2Z, (5) t (∂ +U·∇)ζ(cid:48)+(β+∂ Z)v(cid:48)+u(cid:48)∂ Z = ∂ C =(A +A )C+εΞ, (12) t y x t 1 2 √ =−rζ(cid:48)−ν∆2ζ(cid:48)+ εfe, (6) and forms with Eq. (9) the closed autonomous system of In order to obtain the statistical dynamics of the quasi- S3T theory that determines the statistical dynamics of the linear system (5)-(6) we adopt the general interpretation flow approximated at second order. that the ensemble average over the forcing realizations TheS3Tsystemhasboundedsolutions(cf.AppendixA) is equal to the time average over the intermediate time and the fixed points ZE and CE, if they exist, define sta- scale (Bernstein 2009; Bernstein and Farrell 2010). Under tistical equilibria of the coherent structures with vorticity, this assumption, the slowly varying mean flow Z is also ZE, in the presence of an eddy field with second order the first cumulant of the vorticity Z = (cid:104)ζ(cid:105), where the cumulant or covariance, CE. The structural stability of brackets denote the ensemble average. The time mean of these statistical equilibria addresses the parameters in the the vorticity flux is equal to the ensemble mean of the physical system which can lead to abrupt reorganization of flux T[u(cid:48)ζ(cid:48)] = (cid:104)u(cid:48)ζ(cid:48)(cid:105). The fluxes can be related to the the turbulent flow. Specifically, when an equilibrium of the second cumulant C(x ,x ,t)≡(cid:104)ζ(cid:48)(x )ζ(cid:48)(x )(cid:105), which is the S3T equations becomes unstable as a physical parameter 1 2 1 2 correlation function of the eddy vorticity between the two changes, the turbulent flow bifurcates to a different attrac- points x = (x ,y ), i = 1,2. We hereafter indicate the tor. In this work, we show that coherent structures emerge i i i dynamic variables that are functions of points x =(x ,y ) as unstable modes of the S3T system and equilibrate at i i i with the subscript i, that is ζ(cid:48) ≡ ζ(cid:48)(x ). By making the finite amplitude. The predictions of S3T regarding the i i identification that the fluxes at point x are equal to the emergence and characteristics of the coherent structures value of the two variable function (cid:104)u(cid:48)ζ(cid:48)(cid:105) evaluated at the are then compared to the non-linear simulations of the 1 2 same point x=x =x , we write the fluxes as: stochastically forced barotropic flow. 1 2 (cid:104)u(cid:48)ζ(cid:48)(cid:105)=(cid:104)u(cid:48)ζ(cid:48)(cid:105) . (7) 1 2 x1=x2 2. S3T instability and emergence of finite ampli- Expressing the velocities in terms of the vorticity [u(cid:48),v(cid:48)]= tude large-scale structure [−∂ ∆−1,∂ ∆−1]ζ(cid:48),where∆−1istheintegraloperatorthat y x The homogeneous equilibrium with no mean flow inverts vorticity into the streamfunction field, we obtain the vorticity fluxes as a function of the second cumulant, Ξ ZE =0, CE = , (13) in the following manner: 2r (cid:104)u(cid:48)ζ(cid:48)(cid:105)=(cid:2)(cid:104)u(cid:48)1ζ2(cid:48)(cid:105)x1=x2,(cid:104)v1(cid:48)ζ2(cid:48)(cid:105)x1=x2(cid:3) is a fixed point of the S3T system (9) and (12) in the =(cid:104)−(cid:10)∂ ∆−1ζ(cid:48)ζ(cid:48)(cid:11) ,(cid:10)∂ ∆−1ζ(cid:48)ζ(cid:48)(cid:11) (cid:105) absence of hyperdiffusion (cf. Appendix B). The linear y1 1 1 2 x1=x2 x1 1 1 2 x1=x2 stability of the homogeneous equilibrium can be addressed =(cid:104)−(cid:0)∂ ∆−1C(cid:1) ,(cid:0)∂ ∆−1C(cid:1) (cid:105). (8) byperforminganeigenanalysisoftheS3Tsystemlinearized y1 1 x1=x2 x1 1 x1=x2 aboutthisequilibrium. Theeigenfunctionsinthiscasehave Consequently, the first cumulant evolves according to: the plane wave form ∂tZ+UZx+V(β+Zy)+rZ+ν∆2Z = δZ =Z einx+imyeσt , δC =C (x˜,y˜)einx+imyeσt, nm nm =∂ (cid:0)∂ ∆−1C(cid:1) −∂ (cid:0)∂ ∆−1C(cid:1) . (9) (14) x y1 1 x1=x2 y x1 1 x1=x2 where x˜ = x −x , x = (x +x )/2, y˜ = y −y , y = 1 2 1 2 1 2 Multiplying (A5) for ∂tζ1(cid:48) by ζ2(cid:48) and (A5) for ∂tζ2(cid:48) by ζ1(cid:48), (y1+y2)/2, n and m are the x and y wavenumbers of the adding the two equations and taking the ensemble average eigenfunction and σ =σ +iσ is the eigenvalue with σ = r i r yields the equation for the second cumulant C: Re(σ), σ = Im(σ) being the growth rate and frequency i √ ∂ C−(A +A )C = ε(cid:104)feζ(cid:48) +feζ(cid:48)(cid:105), (10) of the mode respectively. The eigenvalue σ satisfies the t 1 2 1 2 2 1 3 non-dimensional equation: minimum energy input rate required is ε˜ =67 and occurs c at β˜ =3.5, while the minimum zonostrophy parameter ε˜ (cid:90) ∞ (cid:90) ∞ min dk˜d˜l(1−N˜2/K˜2)Ξˆ(k˜,˜l)× required for the emergence of coherent flows is Rβ =0.82 2πr3L2f −∞ −∞ and occurs for β˜→ 0. For β˜≤ β˜min, the structures that (m˜k˜−n˜˜l)(cid:104)n˜m˜(k˜2 −˜l2)+(m˜2−n˜2)k˜ ˜l (cid:105) first become marginally stable are zonal jets (with n=0). + + + + The critical input rate increases as ε˜ ∼β˜−2 for β˜→0 and × (cid:16) (cid:17) = c iβ˜ k˜K˜2−(k˜+n˜)K˜2 +(σ˜+2)K˜2K˜2 the homogeneous equilibrium is structurally stable for all s s excitation amplitudes when β˜=0. However, the structural =(σ˜+1)N˜2−in˜β˜, (15) stability for β˜=0 is an artifact of the assumed isotropy of the excitation and the assumption of a barotropic flow. In where L is a characteristic length scale, σ˜ = σ/r and f the presence of even the slightest anisotropy (Bakas and (n˜,m˜)=L (n,m) are the non-dimensional eigenvalue and f Ioannou 2011, 2013b), or in the case of a stratified flow wavenumbersrespectively,ε˜=ε/(r3L2)isthenon-dimensional f (Parker and Krommes 2015), zonal jets are S3T unstable energy injection rate of the forcing, β˜=βLf/r is the non- and are expected to emerge even in the absence of β. For dimensional planetary vorticity gradient, β˜> β˜ , the marginally stable structures are non-zonal min and ε˜ grows as ε˜ ∼ β˜1/2 for β˜ → ∞. Since the critical 1 (cid:90) ∞ (cid:90) ∞ c c Ξˆ(k,l)= Ξ(x˜,y˜)e−ikx˜−ily˜dx˜dy˜, (16) forcing for the emergence of zonal jets (also shown in fig- 2π −∞ −∞ ure 1), increases as ε˜ ∼β˜2 for β˜→∞, for large values of c βñon-zonalstructuresfirstemergeandonlyatsignificantly is the Fourier transform of the forcing covariance, K˜2 = higher ε˜zonal jets are expected to appear. Investigation of k˜2+˜l2,K˜2 =(k˜+n˜)2+(˜l+m˜)2,N˜2 =n˜2+m˜2,k˜ =k˜+n˜/2 s + these results with other forcing distributions revealed that and ˜l+ =˜l+m˜/2 (cf. Appendix B). For a forcing with the the results for β˜(cid:29) 1 are independent of the structure of mirror symmetry Ξˆ(k,−l)=Ξˆ(k,l) in wavenumber space the forcing (Bakas et al. 2015). and for n˜ (cid:54)=0, the eigenvalues satisfy the relations: The parameter regime of S3T instability is now related to the results of previous studies and to geophysical flows. σ˜ =σ˜∗ , and σ˜ =σ˜ , (17) (−n˜,m˜) (n˜,m˜) (n˜,−m˜) (n˜,m˜) Previous studies have identified a parameter regime which is distinguished by robust, slowly varying zonal jets as implying that the growth rates depend on |n˜| and |m˜|. As well as propagating, non-dispersive, non-zonal coherent a result, the plane wave δZ =cos(nx+my) and an array structures(Galperinetal.2010). Thisregimethatistermed of localized vortices δZ =cos(nx)cos(my), have the same as zonostrophic, is in a region in parameter space in which growth rate, despite their different structure. For zonally the zonostrophy parameter is large (R ≥ 2.5) and the symmetric perturbations with n˜ =0, only the second rela- β scale k = 0.5(β3/ε)1/5 in which anisotropization of the tionin(17)holdsand(15)reducestotheeigenvaluerelation β turbulent spectrum occurs is sufficiently larger than the derived by Srinivasan and Young (2012) for the emergence forcing scale (k /K ≤ 1/4). This regime is shown in of jets in a barotropic β-plane. β f Weconsiderthecaseofaringforcingthatinjectsenergy figure 2 to be highly supercritical for all β˜. In addition, at rate ε at the total wavenumber K : Bakas and Ioannou (2014) calculated indicative order of f magnitude values of β˜ and ε˜for the Earth’s atmosphere (cid:112) Ξˆ(k,l)=2K δ( k2+l2−K ), (18) andoceanaswellasfortheJovianatmosphere. Fromthese f f values we calculated the relevant zonostrophy parameter and obtain the eigenvalues σ˜ by numerically solving (15). R and indicated the three geophysical flows in figure 2. β For small values of the energy input rate, σ˜ < 0 for all We can see that all three cases are supercritical: the Jovian r wavenumbers and the homogeneous equilibrium is stable. atmosphere is highly supercritical and is well within the At a critical ε˜ the homogeneous flow becomes S3T unsta- zonostrophic regime, while the Earth’s atmosphere and c ble and exponentially growing coherent structures emerge. ocean are slightly supercritical (at least within the context The critical value, ε˜ , is calculated by first determining the of the simplified barotropic model). c energyinputrateε˜(n˜,m˜)thatrenderswavenumbers(n˜,m˜) We now examine the growth rate and dispersion prop- t (cid:0) (cid:1) neutral σ˜ =0 , and then by finding the minimum erties of the unstable modes for ε˜> ε˜ and consider first r(n˜,m˜) c energy input rate over all wavenumbers: ε˜ = min ε˜. the case β˜=1, with ε˜=2ε˜ . The growth rate of the maxi- c (n˜,m˜) t c Thecriticalenergyinputrateε˜c asafunctionofβ˜isshown mally growing eigenvalue, σ˜r, and its associated frequency in figure 1. In addition, the corresponding critical zonos- of the mode, σ˜i, are plotted in figure 3(a) as a function of trophy parameter R = 0.7(ε˜ β˜2)1/20 which was used in |n˜| and |m˜|. We observe that the region in wavenumber β c previous studies to characterize the emergence and struc- space defined roughly by 0<|n˜|<1/2, and 1/2<|m˜|<1 ture of zonal jets in planetary turbulence (Galperin et al. is unstable, with the maximum growth rate occurring for 2010), is shown as a function of β˜ in figure 2. The absolute zonal structures (n˜ = 0) with |m˜| (cid:39) 0.8. The frequency 4 (a) 1 0.4 0.8 0.4 0.3 0.3 0.2 0.6 105 σ˜rjets≥σ˜rnon−zon σ˜rjets≤σ˜rnon−zon ββ˜˜22 ˜m||0.4 0.1 00..12 0.2 104 β˜−2 00 0.1 0.2 0.3 0|n˜.4| 0.5 0.6 0.7 0.8 0 ˜ǫc (b) 103 1 7 6 0.8 0.4 0.2 5 102 β˜1/2 0.6 0.6 4 ˜m| |0.4 3 2 10110−2 10−1 100 101 102 103 0.2 1 β˜ 00 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 |n˜| Figure 1: The critical energy input rate ε˜c for structural Figure 3: Dispersion relation of the unstable modes for instability(thicksolidline)andthecriticalenergyinputrate β˜=1(panela)andβ˜=10(panelb). Thecontoursshowthe forstructuralinstabilityofzonaljets(solidline)asafunction of β˜. The behavior of these critical values for β˜ (cid:28) 1 and growthrateσ˜r andtheshadingshowsthefrequencyσ˜i ofthe β˜(cid:29)1isindicatedwiththedashedasymptotes. Inthelight unstablemodes. ForβÕ(1),stationaryzonaljetsaremore unstable and for β˜ (cid:29) 1, westward propagating non-zonal grayregiononlynon-zonalcoherentstructuresemerge,while structures are more unstable. For both panels, the energy inthedarkgrayregionbothzonaljetsandnon-zonalcoherent structures emerge. The thin dotted vertical line β˜=β˜min inputrateisε˜=2ε˜c. separatestheunstableregion: forβ˜<β˜min zonalstructures grow the most, whereas for β˜>β˜min non-zonal structures growthemost. of the unstable modes is zero for zonal jet perturbations (n˜ =0)andnon-negativeforallotherwavenumbers(n˜ =(cid:54) 0). Using the symmetries (17), this implies that real unstable mean flow perturbations δZ propagate in the retrograde direction if n˜ (cid:54)= 0 and are stationary when n˜ = 0. As ε˜ increases the instability region expands and roughly covers the sector 1/2 < |N˜| < 1, with zonal structures having 4.5 a larger growth rate compared to non-zonal structures, a Jupiter result that holds for any ε˜when β˜<β˜ . 4 min For β˜ > β˜ the non-zonal structures have always min 3.5 zonostrophic regime largergrowthrate. Thisisillustratedinfigure3(b),showing the growth rates and frequencies of the unstable modes 3 for β˜ = 10. For larger β˜ values there is a tendency for Rβ2.5 the frequency of the unstable modes to conform to the corresponding Rossby wave frequency 2 ocean 1.5 βñ˜ atmosphere σ˜ = , (19) 1 R n˜2+m˜2 0.5 atendencythatdoesnotoccurforsmallerβ˜. Acomparison 10−2 10−1 100 101 102 103 β˜ betweenthefrequencyoftheunstablemodesandtheRossby wave frequency is shown in figure 4 in a plot of σ˜ /σ˜ . For Figure 2: The critical zonostrophy parameter R = i R β slightly supercritical ε˜, the ratio is close to one and the 0.7(ε˜cβ˜2)1/20 for structural instability (thick line) and the correspondingcriticalparameterforstructuralinstabilityof unstable modes satisfy the Rossby wave dispersion relation. zonaljets(thinline)asafunctionofβ˜. Theshadedregionde- At higher supercriticalities though, σ˜ departs from the i notesthezonostrophicregimeforwhichboththeinequalities Rossby wave frequency (by as much as 40% for the case of R ≥2.5 and k /K ≤1/4 are satisfied. The stars denote β β f ε˜=50ε˜ shown in figure 4(b)). thepositionoftheEarth’satmosphereandoceanaswellas c theJovianatmosphereintheR ,β˜parameterspace. β 5 and Young 2012). Therefore another mechanism should be (a) 1 responsible for producing the upgradient fluxes in the case 0.8 0.99 of an isotropic forcing. 0.98 ˜m|0.6 In order to investigate the eddy-mean flow dynamics | 0.97 underlyingtheS3Tinstability,wecalculatethevorticityflux 0.4 0.96 divergence that is induced when the statistical equilibrium 0.2 0.95 0 0.1 0.2 0.3 0.4 0.5 (13) is perturbed by an infinitesimal coherent structure δZ. |n˜| (b) For an S3T unstable structure, the induced flux divergence 1.5 tends to enhance the coherent structure δZ producing the 1.4 1 positivefeedbackrequiredforinstability. Sothegoalofthis 1.2 ˜m| sectionistoilluminatetheeddy-meanflowdynamicsleading | 1 0.5 to this positive feedback and to understand qualitatively 0.8 0 0.6 why the homogeneous equilibrium is more stable for small 0 0.2 0.4 0.6 |n˜| 0.8 1 1.2 1.4 and large values of β˜. For zonal mean flows (9), (12) are simplified to: Figure 4: Ratio of the frequency of the unstable modes tσ˜hieovsaemrtehweacvoernreusmpbonerdsinσ˜gRfraetqu(ae)ncε˜y=of2aε˜cRaonsdsb(yb)wεa˜v=ew50itε˜hc ∂tU =−∂y(cid:104)u(cid:48)v(cid:48)(cid:105)−rU =∂y(cid:0)∂y22x1∆−11C(cid:1)x1=x2−rU, (20) whenβ˜=100. Valuesofonedenoteanexactmatchwiththe and Rossbywavefrequency. ∂ C =(A +A )C+Ξ, (21) t 1 2 where 3. Analysisoftheeddy-meanflowdynamicsunder- lying jet formation A =−U ∂ −(β−∂2 U)∆−1∂ −r, (22) i i xi yiyi i xi In this section, we investigate the eddy-mean flow dy- respectively. Asaresultthezonalmeanflowisdrivenbythe namics leading to jet formation. These dynamics should momentum flux divergence of the eddies. The perturbation have the property of directly channeling energy from the in vorticity covariance δC that is induced by the mean turbulent motions to the coherent flow without the pres- flow perturbation δU can be estimated immediately by enceofaturbulentcascade. Previousstudieshaveidentified assuming that the system (20)-(21) is very close to the such mechanisms for the maintenance of zonal jets. Huang stability boundary, so that the growth rate is small. In this and Robinson (1998) showed that shear straining of the case the mean flow evolves slow enough that it remains in turbulent field by the jet produced upgradient momentum equilibrium with the eddy covariance, that is dδC/dt(cid:39)0. fluxes that maintained the jet against dissipation. A sim- BakasandIoannou(2013b)showedthattheensemblemean ple case that clearly illustrates the physical picture for the momentumfluxinducedbyaninfinitesimalsinusoidalmean mechanism of shear straining is to consider the evolution of flow perturbation δU = (cid:15)sin(my), where (cid:15) (cid:28) 1 (i.e the eddies in a planar, inviscid constant shear flow. The eddies eigenfunction of (B4)), is equal in this case to the integral are sheared by the mean flow into thinner elliptical shapes, over time and over all zonal wavenumbers of the responses while their vorticity is conserved. For an elongated eddy to all point excitations in the y direction: this implies that the eddy velocities decrease and the eddy energy is transferred to the mean flow through upgradient δ(cid:104)u(cid:48)v(cid:48)(cid:105)= 1 (cid:90) ∞ (cid:90) ∞ (cid:90) ∞u(cid:48)v(cid:48)(t)dtdξdk, (23) momentum fluxes. This mechanism can operate when the 2π −∞ −∞ 0 time required for the eddies to shear over is much shorter than the dissipation time scale. The reason is that in this where u(cid:48)v(cid:48)(t) is the momentum flux at time t produced by: limit even the eddies with streamfunctions leaning against G(k,y−ξ)=B(k)h(y−ξ)eikx+il0(y−ξ). (24) the shear that initially widen significantly gaining momentum, have the necessary time to shear over, elongate and The Green’s function G has the form of a wavepacket with surrender their momentum to the mean flow. Given that an amplitude B(k) and a carrier wave with wavenumbers for an emerging jet the characteristic shear time scale is (k,l )thatismodulatedintheydirectionbythewavepacket 0 necessarily infinitely longer than the dissipation time scale, envelope h(y). The characteristics of the amplitude, the it needs to be shown that shear straining can produce up- wavenumber and the envelope depend on the forcing char- gradient momentum fluxes in this case as well. In addition, acteristics, but in any case the calculation of the ensemble previousstudieshaveshownthatshearingofisotropiceddies mean momentum fluxes is reduced to calculating the mo- on an infinite domain does not produce any net momentum mentum fluxes over the life cycle of wavepackets that are fluxes (Shepherd 1985; Farrell 1987; Holloway 2010) and initially at different latitudes and then adding their relative should have no effect on the S3T instability (Srinivasan contributions. 6 As the wavepacket propagates in the latitudinal direc- and slightly changes the amplitude of the fluxes as well as tion, its meridional wavenumber and frequency are going slightly speeds up or slows down the wavepacket. The sum to change due to shearing by the mean flow and due to of these two effects will produce the induced momentum the change of the mean vorticity gradient β −U . The fluxes. yy resulting time variable momentum flux u(cid:48)v(cid:48)(t) can be calculated using ray tracing. According to standard ray tracing a. Thelimitofsmallscalewavepacketswithashortpropagation arguments, the wave action is conserved along a ray (in the range absence of dissipation) leading to the momentum flux: In order to clearly illustrate the behavior of the eddy fluxes, we consider the limit of β˜=βL /r (cid:28)1, where L u(cid:48)v(cid:48)(t)=−|B|2A (t)e−2rt|h(y−η(t))|2, (25) f f M is the scale of the wavepackets and in addition we assume where A (t)=kl /(k2+l2)2 is the momentum flux of the that the scale of the mean flow, 1/m, is much larger than M t t the scale of the wavepackets mL (cid:28) 1. In this limit, the carrier wave that determines the amplitude of the fluxes of f wavepackets are dissipated before propagating far from the the wavepacket and l , η(t) are the time dependent merid- t initial position and the effect of the change in the mean ional wavenumber and position of the wavepacket respec- vorticity gradient is higher order. As a result, Bakas and tively(Andrewsetal.1987). Becauseofthesmallamplitude Ioannou(2013b)showthatl andη decreasemonotonically of the mean flow perturbation δU, the wavenumber and 1 1 with time with rates independent of δU and proportional position of the packet vary slowly on a time scale O((cid:15)t) yy to the shear δU (ξ) at the initial position ξ: compared to the dissipation time scale 1/r and the domi- y nant contribution to the time integral in (23) comes from (cid:18) (cid:19) dA small times. We can therefore seek asymptotic solutions of l =−δU (ξ)kt , η =−βδU (ξ) M kt2. (28) 1 y 1 y dl the form t l0 Thatis,theamplitudeofthefluxA andthegroupvelocity l =l +(cid:15)l +··· , η(t)=ξ+c t+(cid:15)η (t)+··· , (26) M t 0 1 0 1 of the packets change only due to the shearing of the phase where c = 2βkl /(k2 +l2)2 is the group velocity in the lines of the carrier wave according to the local shear. 0 0 0 absence of a mean flow and calculate the integral of u(cid:48)v(cid:48)(t) Consider in this limit the first term, u(cid:48)v(cid:48)S, arising from the small amplitude change. Since the wavepacket is dis- over time from the leading order terms. Substituting (26) sipated before it propagates away, we can ignore to first in (25) we obtain: order propagation: u(cid:48)v(cid:48)(t)=−|B|2A (0)e−2rt|h(y−ξ−c t)|2− M 0 (cid:18) (cid:19) dA (cid:124) (cid:123)(cid:122) (cid:125) u(cid:48)v(cid:48) =−(cid:15)|B|2 M l (t)e−2rt|h(y−ξ−c t)|2 −(cid:124) (cid:15)|B|2(cid:18)ddAlMt (cid:19)ul(cid:48)0v(cid:48)lR1(t(cid:123))(cid:122)e−2rt|h(y−ξ−c0t)|(cid:125)2 S (cid:39)|B|2δUy(ξ)kdtlt(cid:18)ddAll0Mt 1(cid:19)l0e−2rt|h(y−ξ)|2,0 (29) u(cid:48)v(cid:48)S so that the packet grows/decays in situ. Since the wave d −(cid:15)|B|2A (0)η (t)e−2rt |h(y−ξ−c t)|2. packet is rapidly dissipated, the integrated momentum flux M 1 dy 0 over its life time will be given to a good approximation by (cid:124) (cid:123)(cid:122) (cid:125) u(cid:48)v(cid:48)β the instantaneous change in the flux1 that is proportional (27) to (dA /dl ) . Figure 5 illustrates the amplitude of the M t l0 momentum flux as a function of the angle θ =arctan(l /k) t t The first term, u(cid:48)v(cid:48) , arises from the momentum flux pro- of the phase lines of the carrier wave of the packet with the R duced by a wavepacket in the absence of a mean flow. y-axis. Itisshownthatthemomentumfluxofawavepacket Because A (0) = kl /(k2 +l2)2 is odd with respect to with |θ |<π/6 (that is with phase lines close to the merid- M 0 0 0 wavenumbers, this term does not contribute to the en- ional direction) excited in regions II or III, will increase semble averaged momentum flux when integrated over all within the dissipation time scale. Compared to an un- wavenumbers and will be hereafter ignored. The second sheared wavepacket, this process leads to upgradient mo- term, u(cid:48)v(cid:48) , arises from the small change in the amplitude mentum flux. The opposite occurs for waves excited in S of the flux A over a dissipation time scale. The third regions I and IV (with |θ |>π/6) that produce downgradi- M 0 term, u(cid:48)v(cid:48) , arises from the small change in the position of ent flux, as their momentum flux decreases. β the packet η compared to a propagating packet in the ab- We now consider the second term, u(cid:48)v(cid:48) arising from β senceofameanflow. Tosummarize,theinfinitesimalmean theeffectofpropagationonthemomentumflux. Thegroup flow refracts the wavepacket due to shearing by the mean 1occurringoverthedissipationtimescale1/r thatisincremental flow and due to the change of the mean vorticity gradient insheartimeunits 7 (a) 0.5 0.2 I II 0.4 0.1 0.3 uv 0 0.2 −0.1 0.1 −0.2 −3 −2 −1 0 1 2 3 M 0 y A x 10−5 (b) −0.1 5 −0.2 0 −0.3 uvβ −5 −0.4 III IV −0.5 −10 −80 −60 −40 −20 0 20 40 60 80 −3 −2 −1 0 1 2 3 θt=arctan(lt/k) y Figure 5: Amplitude of the momentum fluxes, AM(t), of Figure 6: (a) Comparison of the momentum fluxes of an wavepackets as a function of the angle θt = arctan(lt/k) unshearedwavepacketexcitedinregionsII(thicksolidline) between the phase lines of the central wave and the y-axis. and III (solid line) to the momentum fluxes of a sheared Theverticallinesseparatetheregionswith|θt|<π/6(IIand wavepacketshownbythecorrespondingdashedlines,when III)and|θt|>π/6(IandIV). only the change in propagation is taken into account. A snapshotofthefluxesatt=0.2/r isshown. Theplanetary vorticitygradientisβ=0.1,thewavepackethasinitialvortic- (cid:113) velocityisgivenbycg =2βAM inthiscaseandasaresulta ityh(y)=e−y2, k2+l02=1,|θ0|=π/10and|B|=1. (b) wavepacketstartinginregionIII,willpropagatetowardsthe Thedifferenceinmomentumfluxesbetweenashearedandan north(c.f. figure5). Becauseshearingslowsdownthewaves unshearedwavepacketcalculatedovertheirlifecycle. in region III (η ∼−(dA /dl )), the wavepacket will flux 1 M t its momentum from southern latitudes compared to when proportional to the third derivative of δU yielding a hyper- it moved in the absence of the shear flow. This is shown in diffusive momentum flux divergence that tends to reenforce figure6(a)illustratingthedistributionofmomentumfluxof the mean flow and is therefore destabilizing. These destabi- anunshearedanda shearedperturbationwhoseamplitudes lizing fluxes are proportional to β˜2 and as a result, the en- are constant. Figure 6(b) plots this difference, u(cid:48)v(cid:48)β, and ergyinputraterequiredtoformzonaljetsincreasesas1/β˜2 shows that the flux is downgradient in this case. The same inthislimit. Itisworthnotingthatthefirsttermintegrates happens for waves excited in region II, while the waves to zero only for the special case of the isotropic forcing, as excited in regions I and IV produce upgradient flux. even the slightest anisotropy yields a non-zero contribution The net momentum fluxes produced by an ensemble of from u(cid:48)v(cid:48) . For example consider the forcing covariance S wavepackets, will therefore depend on the spectral charac- Ξ(x1,x2,y1,y2)=cos(k(x1−x2))e−(y1−y2)2/δ2 that mim- teristics of the forcing that determine the regions (I-IV), in ics the forcing of the barotropic flow by the most unstable which the forcing has significant power. Bakas and Ioannou baroclinic wave, which has zero meridional wavenumber. In (2013b) show that for the isotropic forcing (18): this case the forcing that is centered at l =0 in wavenum- 0 ber space, injects significant power in a band of waves in δ(cid:104)u(cid:48)v(cid:48)(cid:105)= regions II and III and therefore u(cid:48)v(cid:48) yields upgradient 1 (cid:90) ∞ (cid:90) ∞ 1 (cid:90) ∞ (cid:90) ∞ S = u(cid:48)v(cid:48) dξdk+ u(cid:48)v(cid:48) dξdk fluxes. 2π S 2π β −∞ −∞ −∞ −∞ 3ε˜β˜2r d3δU b. The effects of the change in the mean vorticity gradient and (cid:39)0− . (30) the finite propagation range 32πK4 dy3 f In order to take into account the effect of the change in The first integral is zero, because the gain in momentum the vorticity gradient, we retain higher order terms with occurring for |θ0|<π/6 (waves excited in regions II, III) is respect to mLf (cid:28) 1 in l1 and η1. In this case it can be fully compensated by the loss in momentum for |θ0|>π/6 shown that l1 decreases with time at a rate proportional (waves excited in regions I, IV) since for the isotropic forc- to U (ξ)+U (ξ) (Bakas and Ioannou 2013b). Since the y yyy ing all possible wave orientations are equally excited. The local shear and the local change in the vorticity gradient net momentum fluxes are therefore produced by the u(cid:48)v(cid:48)β have different signs, the wavepacket is ’sheared less’ and as term and are upgradient, because the loss in momentum a result we expect reduced momentum fluxes compared to occurringfor|θ0|<π/6,isovercompensatedbythegainin the limit discussed in section 4.2. Indeed, for the isotropic momentum for |θ |>π/6. The momentum fluxes are also 0 8 forcing: (a) x 10−4 x 10−4 (b) 6 5 3 (cid:32) (cid:33) 3ε˜β˜2r d3δU 1 d5δU δ(cid:104)u(cid:48)v(cid:48)(cid:105)(cid:39)− − . (31) 4 2 32πK4 dy3 4K2 dy5 y 0 E f f 2 1 −5 That is, the change in the mean vorticity gradient has a 0 0 stabilizing effect. 0 x 5 0 50t0 1000 We finally relax the assumption that β˜ (cid:28) 1. In this (c) x 10−3 (d) 6 1040 case, l and η are affected by an integral shear and mean 1 1 5 1030 vorticitygradientovertheregionofpropagation. Forlarger 4 β˜, the wavepacket will encounter regions of both positive y 0 t1020 2 and negative shear and as a result, the momentum fluxes −5 1010 thatarequalitativelyproportionaltotheintegralshearover 0 1000 the propagation region will be reduced. In the limit β˜(cid:29)1, 0 x 5 0 2 x 4 6 the region of propagation is the whole sinusoidal flow with consecutive regions of positive and negative shear and the Figure 7: EquilibrationoftheS3Tinstabilities. (a)Stream- function of the initial perturbation. (b) Energy evolution integralshearalongwiththefluxeswillasymptoticallytend of the initial perturbation shown in panel (a) as obtained to zero. As a result, the energy input rate required for from the integration of the S3T equations (9) and (12) structural instability of zonal jets increases with β˜ in this (dashedline)andfromtheintegrationoftheensemblequasi- limit. linear(EQL)system(4)-(3)withNens=10(solidline)and Nens=100(dash-dottedline)ensemblemembersthatisdis- cussedinsection6. (c)SnapshotofthestreamfunctionΨeq 4. Equilibration of the S3T instabilities of the traveling wave structure and (d) Hovmöller diagram ofΨeq(x,y=π/4,t)forthefiniteequilibratedtravelingwave. We now investigate the equilibration of the instabilities Thethickdashedlineshowsthephasespeedobtainedfrom bystudyingtheS3Tsystem(9),(12)discretizedinadoubly thestabilityequation(15). Theenergyinputrateisε˜=4ε˜c periodic channel of size 2π × 2π. We approximate the andβ˜=100. monochromatic forcing (18), by considering the narrow band forcing form to the harmonic Z =cos(x)cos(5y) that propagates √ (cid:26) westward. This is illustrated in the Hovmöller diagram of K 1, for | k2+l2−K |≤∆K Ξˆ(k,l)= f √ f f , (32) ψ(x,y =π/4,t) shown in 7(d). The sloping dashed line in ∆K 0, for | k2+l2−K |>∆K f f f thediagramcorrespondstothephasespeedofthetraveling where k, l assume integer values, that injects energy at rate wave, which is found to be approximately the phase speed εinanarrowringinwavenumberspacewithradiusK and of the unstable (n,m)=(1,5) eigenmode. f width∆Kf. Weconsiderthesetofparametervaluesβ =10, Considernowtheenergyinputrateε˜=10ε˜c. Whilethe r =0.01, ν =1.19·10−6, K =10 and ∆K =1, for which maximum growth rate still occurs for the (|n|,|m|)=(1,5) f f β˜ = 100. The integration is therefore in the parameter non-zonal structure, zonal jet perturbations are unstable as well. If the S3T dynamics are restricted to account only region of figure 1 in which the non-zonal structures are for the interaction between zonal flows and turbulence by more unstable than the zonal jets. The growth rates of the employing a zonal mean rather than an ensemble mean, coherent structures for integer values of the wavenumbers, an infinitesimal jet perturbation will grow and equilibrate nandmarecalculatedfromthediscreteversionofequation at finite amplitude. To illustrate this we integrate the (15) obtained by substituting the integrals with sums over S3T dynamical system (20)-(21) restricted to zonal flow integer values of the wavenumbers (Bakas and Ioannou coherent structures. The energy of the small zonal jet 2013a). perturbation δZ =0.1cos(4y) is shown in figure 8 to grow We first consider the supercritical energy input rate and saturate at a constant value and the streamfunction ε˜=4ε˜ . Fortheseparametersonlynon-zonalmodesareun- c of the equilibrated jet is shown in the left inset in figure stable, with the perturbation with (n,m)=(1,5) growing 8. However, in the context of the generalized S3T analysis the most. At t=0, we introduce a small random perturba- that takes into account the dynamics of the interaction tion, whose streamfunction is shown in figure 7(a). After between coherent non-zonal structures and jets, we find a few e-folding times, a harmonic structure of the form that these S3T jet equilibria can be saddles: stable to zonal Z = cos(x)cos(5y) dominates the large-scale flow. The jet perturbations but unstable to non-zonal perturbations. energy of this large scale structure shown in figure 7(b), To show this, we consider the evolution of the same jet increases rapidly and eventually saturates. At this point perturbation δZ = 0.1cos(4y) under the generalized S3T the large-scale flow gets attracted to a traveling wave finite dynamics(9),(12)andfindthattheflowfollowsthezonally amplitude equilibrium structure (cf. figure 7(c)) close in 9 x 10−3 10−2 1 6 0.02 10−3 4 0 y 0.8 2 −0.02 E0.6 46 0.01 6 Ez10−4 00 x 5 6 0.02 0 0.01 4 y 0 2 4 y 0.4 00 5 −0.01 y2 −00.01 10−5 02 −0.02 0.2 x 0 0 x 5 0 5 x 0 10−6 0 500 1000 1500 2000 2500 3000 3500 0 200 400 600 800 1000 t t Figure 8: Energy evolution of an initial jet perturbation Figure 9: Zonalenergyevolutionofarandomzonalpertur- δZ = 0.1cos(4y) for the zonally restricted S3T dynamics bationimposedonthenon-zonaltravelingwaveequilibrium (20)-(21)(dashedline)andthegeneralizedS3Tdynamics(9), shownintheleftinset. Thestreamfunctionoftheequilibrated (12) (thin line). The insets show a snapshot of the mean structureisshownintherightinset. Theenergyinputrate flowstreamfunctionatt=700(left)andthestreamfunction isε˜=30ε˜c andβ˜=100. of the equilibrated structure at t = 3500 (right) under the generalizedS3Tdynamics. Theparametersareε˜=10ε˜c and β˜=100. compared to the mixed state in figure 9. Embedded in it are non-zonal vortices with the same meridional scale and with about 14% the energy of the zonal jet. These vortices restricted S3T dynamics and equilibrates to the same finite that are shown in figure 10(b) to have approximately the amplitude zonal jet (cf. figure 8). At this point we insert a compact support structure Ψ = cos(x)cos(4y) propagate small random perturbation to the equilibrated flow. Soon westwardasshownintheHovmo¨llerdiagraminfigure10(c). after, non-zonal undulations grow and the flow transitions to the stable Z =cos(x)cos(5y) traveling wave state that 5. Comparison to ensemble mean quasi-linear and is also shown in figure 8. As a result, the finite equilibrium non-linear simulations zonal jet structure is S3T unstable to coherent non-zonal perturbations and is not expected to appear in non-linear a. Comparison to an ensemble of quasi-linear simulations simulations despite the fact that the zero flow equilibrium Withinthecontextofthesecondordercumulantclosure, is unstable to zonal jet perturbations. We will elaborate the S3T formulation allows the identification of statistical more on this issue in the next section. turbulent equilibria in the infinite ensemble limit, in which Finally,considerthecaseε˜=30ε˜ . Atthisenergyinput c the fluctuations induced by the stochastic forcing are av- rate,thefiniteamplitudenon-zonaltravelingwaveequilibria eraged to zero. However, these S3T equilibria and their become S3T unstable. To show this, we consider the non- stability properties are manifest even in single realizations zonal traveling wave equilibrium obtained by the evolution of the turbulent system. For example, previous studies us- ofthesmallnon-zonalperturbationδZ =0.01cos(x)cos(5y) ing S3T obtained zonal jet equilibria in barotropic, shallow to the homogeneous state that is shown at the left inset water and baroclinic flows in close correspondence with in figure 9 and impose a small random zonal perturbation. observed jets in planetary flows (Farrell and Ioannou 2007, 2 The evolution of the zonal energy E =(1/2)U , where the z 2008,2009b,a). Inaddition,previousstudiesofS3Tdynam- overbar denotes a zonal average, is shown in figure 9. After ics restricted to the interaction between zonal flows and an initial transition period, the zonal perturbations grow turbulence in a β-plane channel showed that when the en- exponentiallyandtheflowtransitionstothejetequilibrium ergy input rate is such that the zero mean flow equilibrium state shown at the right inset in figure 9. Note however, is unstable, zonal jets also appear in the non-linear simu- that the jet equilibrium structure is not zonally symmetric. lations with the structure (scale and amplitude) predicted This is a new type of S3T equilibrium: it is a mix between by S3T (Srinivasan and Young 2012; Constantinou et al. a zonal jet and a non-zonal traveling wave with the same 2014). meridional scale. These mixed equilibria appear to be the A very useful intermediate model that retains the wave- attractors for larger energy input rates as well. This is mean flow dynamics of the S3T system while relaxing the illustrated in figure 10 showing the structure of the mixed infinite ensemble approximation is the quasi-linear system equilibrium at ε=50ε . The equilibrium structure consists c (5)-(6). Under the ergodic assumption, this can be inter- of a large amplitude zonally symmetric jet with larger scale 10

Emergence of non-zonal coherent structures

Similar Emergence of non-zonal coherent structures