Energy stability of the Charney-DeVore quasi-geostrophic equation for atmospheric blocking
aa r X i v : . [ m a t h . A P ] F e b Energy stability of the Charney–DeVore quasi-geostrophic equation foratmospheric blocking
Zhi-Min Chen a , Xiangming Xiong b, ∗ a School of Mathematics and Statistics, Shenzhen University, Shenzhen 518060, China b The Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ 08544, USA
Abstract
Charney and DeVore [J. Atmos. Sci. 36 (1979), 1205-1216] found multiple equilibrium states as aconsequence of bottom topography in their pioneering work on the quasi-geostrophic barotropicflow over topography in a β -plane channel. In the present paper, we prove that the basicflow is asymptotically stable in a parameter region, including the flat topography situation,which excludes the existence of multiple equilibrium states therein. Moreover, we show that anadditional condition on the average zonal force or the average zonal velocity is indispensable tothe well-posedness of the Charney–DeVore quasi-geostrophic equation. Coexistence of at leastthree equilibrium states is confirmed by a pseudo-arclength continuation method for differenttopographic amplitudes. The stabilities of the equilibrium states are examined by high-resolutiondirect numerical simulations. Keywords: quasi-geostrophic equation, barotropic flow, beta-plane channel model, topographiceffect, nonlinear stability, pseudo-arclength continuation method
1. Introduction
A surface quasi-geostrophic flow is the first order approximation of a geostrophic flow withrespect to the quasi-geostrophic approximation under a small Rossby number [1]. From themathematical formulation of [2], the energy dissipation controlled by the fractional Laplacian κ ( − ∆) α is applied to a quasi-geostrophic flow [3, 4]. Mathematical theory of the dissipativequasi-geostrophic equation has been extensively studied (see, for example, [3, 5–10]). In thepresent study, we are interested in a quasi-geostrophic flow, which is from the understanding ofblocks in atmosphere and involves the weak energy dissipation κ ( − ∆) α with α = 0.Large-scale flows in the atmosphere are strongly influenced by the bottom topography. Whenthe topography is not flat, the small-scale and large-scale components of atmospheric flow in-teract via not only advection but also topographic stress. Charney and DeVore [11] proposed abarotropic channel model to study the influence of topography on the atmospheric flow. Theyfound multiple equilibrium states in the quasi-geostrophic β -plane flow over topography withhighly truncated spectral expansions. In some region of parameter space, there were one un-stable equilibrium state and two stable equilibrium states. The two stable equilibrium states ∗ Corresponding author.
Email addresses: [email protected] (Zhi-Min Chen), [email protected] (Xiangming Xiong)
Accepted by SCIENTIA SINICA Mathematica (in Chinese) omprised a “high-index” flow with a relatively larger zonal flux and a “low-index” flow showingatmospheric blocking. Atmospheric blocking refers to the unusual persistence of a large-scaleweather pattern in an area, and usually occurs when the normal midlatitude eastward flow isinterrupted by topography or strong meridional flow [12, 13]. Atmospheric blocking may bringabout extreme weather events such as floods, droughts, and persistent abnormal temperatures[14].In the work of Charney and DeVore [11], the governing equation of the vertical vorticity ζ was ∂ζ∂t + ∂ψ∂x ∂ ( ζ + h + βy ) ∂y − ∂ψ∂y ∂ ( ζ + h + βy ) ∂x = − kζ − ∂F∂y , (1)where ψ was the stream function satisfying ζ = ∂ ψ/∂x + ∂ ψ/∂y , and h was the bottomtopography. β = 2Ω(cos φ ) /a , where Ω was the angular speed of the earth’s rotation, a was theradius of the earth, and φ was the average latitude of the zonal channel. k was a coefficient offriction exerted by the Ekman layer near the ground. − ∂F/∂y was the vorticity source arisenfrom the zonal body force F .Charney and DeVore [11] assumed the flow was periodic in the zonal direction with a non-dimensional period of 2 π . The meridional boundaries y = 0 and y = π were assumed to benon-penetrative. Therefore, the boundary conditions of the stream function ψ were ψ ( x + 2 π, y, t ) = ψ ( x, y, t ) , ∂ψ∂x (cid:12)(cid:12)(cid:12)(cid:12) y =0 = ∂ψ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) y = π = 0 . (2)However, the governing equation (1) was not well-posed under the boundary conditions (2).Although the average zonal force did not appear in the governing equation (1), it evidentlyinfluenced the atmospheric flow. Therefore, another equation of the average zonal velocity isneeded to take into account this influence. This problem was noticed by Davey [15]. He proposeda quasi-linear theory for steady flow over topography in a periodic channel, which included theequation of the average zonal velocity but neglected the advection due to the mean shear andthe interactions between eddies. He confirmed that there were multiple equilibrium states whenthe amplitude of topography was sufficiently large.The multiple equilibrium states of the Charney–DeVore quasi-geostrophic equation were anal-ysed and calculated by many researchers. Hart [16] found stable equilibrium states when thetopography was anisotropic and mainly varied in the zonal direction. In the fully-nonlinearcalculation performed by Legras and Ghil [17], however, the parameter region where multipleequilibrium states existed was found to be smaller than that in a quasi-linear model. Tung andRosenthal found that the multiple equilibrium states only existed in a parameter region that wasfar from realistic if fully-nonlinear effects were considered and the truncation of modes was notsevere [18]. However, multiple equilibrium states were still found with wind speeds only twiceas large as the observed ones by introducing mid-latitudinal zonal jets in the model [19]. Otherresearchers focused on the more complicated case with the baroclinic model [20] or two layer flow[21]. Multiple equilibrium states from a bifurcation viewpoint of a non-parallel vortex flow wereinvestigated in [22]. One may also refer to [23] for the Hopf bifurcation of a quasi-geostrophicatmospheric flow involving strong energy dissipation.2hen the potential vorticity is proportional to the stream function, the nonlinear Jacobianterm in the equation of potential vorticity will vanish. The inviscid stabilities of the equilibriumstates with potential vorticity proportional to the stream function have already been analysedin a β -plane channel [24, 25]. In the absence of forcing and dissipation, the nonlinear stability ofthe equilibrium states was studied with the conservation of a combinition of the kinetic energyand the potential enstrophy in a double-periodic domain [26, 27]. Zou and Fyfe considered theenergy stability of the equilibrium states of forced and dissipated barotropic flow [28]. Energystability theory of general parallel flows was recently investigated by the authors [29].In [30], we proved the global stability of the equilibrium state in a parameter region, andperformed numerical simulations with 12 Fourier modes in the zonal direction and 65 colloca-tion points in the meridional direction. In this paper, we will prove the global stability of theequilibrium state in a larger parameter region. After that, we will calculate the equilibriumstates with a pseudo-arclength continuation method, and will examine their stabilities with di-rect numerical simulations with higher resolution. The derivation of the governing equations ofthe quasi-geostrophic barotropic flow is presented in Section 2. The proof of the global stabilityis given in Section 3. After briefly introducing the numerical methods in Section 4, we presentthe numerical results of multiple equilibrium states in Section 5. Discussions on the equilibriumstates found by Charney and DeVore are made in Section 6, and conclusions are given in Section7.
2. Mathematical model
The large-scale flow of atmosphere outside the Ekman boundary layer is nearly inviscid, andtherefore the governing equations and boundary conditions are ∇ ∗ · U ∗ = 0 , (3) ∂ U ∗ ∂t ∗ + ( U ∗ · ∇ ∗ ) U ∗ + f ∗ × U ∗ = − ρ ∇ ∗ P ∗ + 1 ρ F ∗ , (4)( U ∗ , P ∗ )( x ∗ + 2 πL e x , t ∗ ) = ( U ∗ , P ∗ )( x ∗ , t ∗ ) , (5) V ∗ | y ∗ =0 = V ∗ | y ∗ = πL = 0 , (6) W ∗ | z ∗ = h ∗ = U ∗ | z ∗ = h ∗ ∂h ∗ ∂x ∗ + V ∗ | z ∗ = h ∗ ∂h ∗ ∂y ∗ , W ∗ | z ∗ = H = 0 , (7)where x ∗ = ( x ∗ , y ∗ , z ∗ ), and x ∗ and y ∗ are the zonal (eastward) and meridional (northward)coordinates in the tangent plane of the earth, respectively. The coordinate in the directionperpendicular to the tangent plane (opposite to the direction of the gravitational acceleration) is z ∗ . The unit vectors in these directions are e x , e y , and e z . Assume the latitude φ to be confinedto a small interval around a middle latitude φ in the Northern Hemisphere, and then we have x ∗ ≈ aθ cos φ and y ∗ ≈ πL/ a ( φ − φ ), where θ , a , and πL are the longitude, the radius of theearth, and the distance between zonal walls, respectively. U ∗ = ( U ∗ , V ∗ , W ∗ ) is the velocity ofthe atmosphere. The density of the atmosphere ρ is assumed to be a constant. H is the heightof the upper free surface, and h ∗ ( x ∗ , y ∗ ) is the lower topography, which satisfies Z Z h ∗ d x ∗ d y ∗ = 0 . β -plane approximation, f ∗ ≈ φ ) e z ≈ ( f ∗ + β ∗ y ∗ ) e z , where Ω is theangular speed of the earth’s rotation, f ∗ = 2Ω sin φ − Ω πL (cos φ ) /a , and β ∗ = 2Ω(cos φ ) /a .Only the local normal component of the earth’s rotation is preserved here [1]. The body forceis F ∗ = F ∗ e x − ρg e z , where g is the magnitude of the gravitational acceleration, and F ∗ is theexternal force acting in the zonal direction. In the above definitions, asterisks denote dimensionalvariables. The constants H , L , a , Ω, ρ , and g are also dimensional quantities.Suppose | h ∗ | ≪ H ≪ πL ≪ a , and therefore the vertical scale of the geophysical flowis far less than its horizontal scale. Neglecting W ∗ , ∂U ∗ /∂z ∗ , and ∂V ∗ /∂z ∗ in the momentumequations, we have ∂U ∗ ∂x ∗ + ∂V ∗ ∂y ∗ + ∂W ∗ ∂z ∗ = 0 , (8) ∂U ∗ ∂t ∗ + U ∗ ∂U ∗ ∂x ∗ + V ∗ ∂U ∗ ∂y ∗ − ( f ∗ + β ∗ y ∗ ) V ∗ = − ρ ∂P ∗ ∂x ∗ + 1 ρ F ∗ , (9) ∂V ∗ ∂t ∗ + U ∗ ∂V ∗ ∂x ∗ + V ∗ ∂V ∗ ∂y ∗ + ( f ∗ + β ∗ y ∗ ) U ∗ = − ρ ∂P ∗ ∂y ∗ , (10)0 = − ρ ∂P ∗ ∂z ∗ − g. (11)The vertical momentum equation (11) implies ∂∂z ∗ ( ∂P ∗ ∂x ∗ ) = ∂∂z ∗ ( ∂P ∗ ∂y ∗ ) = 0 . Assume F ∗ = F ∗ ( x ∗ , y ∗ ), and then all terms in the horizontal momentum equations (9)–(10) areindependent of z ∗ .Integrating the continuity equation (8) from z ∗ = h ∗ to z ∗ = H , we have W ∗ | z ∗ = H − W ∗ | z ∗ = h ∗ = − Z Hh ∗ ( ∂U ∗ ∂x ∗ + ∂V ∗ ∂y ∗ ) d z ∗ = − ( H − h ∗ )( ∂U ∗ ∂x ∗ + ∂V ∗ ∂y ∗ ) . By substituting the boundary condition (7) into the above equation, we obtain ∂∂x ∗ (( H − h ∗ ) U ∗ ) + ∂∂y ∗ (( H − h ∗ ) V ∗ ) = 0 . Let ˜ U ∗ = (1 − h ∗ /H ) U ∗ and ˜ V ∗ = (1 − h ∗ /H ) V ∗ , and then the equations (8)–(10) areapproximated by ∂ ˜ U ∗ ∂x ∗ + ∂ ˜ V ∗ ∂y ∗ = 0 , (12)(1 + h ∗ H ) (cid:16) ∂ ˜ U ∗ ∂t ∗ + ( ˜ U ∗ ∂∂x ∗ + V ∗ ∂∂y ∗ )( ˜ U ∗ + h ∗ H ˜ U ∗ ) − ( f ∗ + β ∗ y ∗ ) ˜ V ∗ (cid:17) = − ρ ∂P ∗ ∂x ∗ + 1 ρ F ∗ , (13)(1 + h ∗ H ) (cid:16) ∂ ˜ V ∗ ∂t ∗ + ( ˜ U ∗ ∂∂x ∗ + V ∗ ∂∂y ∗ )( ˜ V ∗ + h ∗ H ˜ V ∗ ) + ( f ∗ + β ∗ y ∗ ) ˜ U ∗ (cid:17) = − ρ ∂P ∗ ∂y ∗ , (14)4here (1 − h ∗ /H ) − ≈ h ∗ /H is used.The characteristic length, velocity, and time of large-scale atmospheric motions are L ∼ km, U ∗ ∼
10 m / s, and L/U ∗ ∼ s ≈
28 hours, respectively [31, 32]. The average depth ofthe atmosphere and the characteristic elevation of topography are supposed to be H ∼
10 km[32] and h ∗ ∼ φ ≈ ◦ N. Using a ≈ g ≈
10 m / s ,and Ω ≈ π/ − , we have f ∗ ∼ − s − and β ∗ ∼ − m − s − . The Rossby number is Ro = U ∗ / ( f ∗ L ) ∼ .
1. Then we have the following estimates: f ∗ ˜ U ∗ , f ∗ ˜ V ∗ ∼ f ∗ U ∗ ∼ − m / s ,h ∗ H f ∗ ˜ U ∗ , h ∗ H f ∗ ˜ V ∗ ∼ h ∗ H f ∗ U ∗ ∼ − m / s ,β ∗ y ∗ ˜ U ∗ , β ∗ y ∗ ˜ V ∗ ∼ β ∗ LU ∗ ∼ − m / s ,∂ ˜ U ∗ ∂t ∗ , ∂ ˜ V ∗ ∂t ∗ ∼ U ∗ L ∼ − m / s , ˜ U ∗ ∂ ˜ U ∗ ∂x ∗ , ˜ V ∗ ∂ ˜ U ∗ ∂y ∗ , ˜ U ∗ ∂ ˜ V ∗ ∂x ∗ , ˜ V ∗ ∂ ˜ V ∗ ∂y ∗ ∼ U ∗ L ∼ − m / s . It follows that the terms on the right hand side of the equations (13)–(14) are on the order of10 − m / s .Neglecting terms that have magnitudes less than 10 − m / s in the equations (12)–(14), wehave ∂ ˜ U ∗ ∂x ∗ + ∂ ˜ V ∗ ∂y ∗ = 0 , (15) ∂ ˜ U ∗ ∂t ∗ + ˜ U ∗ ∂ ˜ U ∗ ∂x ∗ + V ∗ ∂ ˜ U ∗ ∂y ∗ − ( f ∗ + f ∗ h ∗ H + β ∗ y ∗ ) ˜ V ∗ = − ρ ∂P ∗ ∂x ∗ + 1 ρ F ∗ , (16) ∂ ˜ V ∗ ∂t ∗ + ˜ U ∗ ∂ ˜ V ∗ ∂x ∗ + V ∗ ∂ ˜ V ∗ ∂y ∗ + ( f ∗ + f ∗ h ∗ H + β ∗ y ∗ ) ˜ U ∗ = − ρ ∂P ∗ ∂y ∗ . (17)In the above derivation, we focus on the large-scale flow of atmosphere outside the Ekmanboundary layer. However, the ground exerts friction to the flow of atmosphere through theEkman boundary layer. Denote the zonal velocity, the meridional velocity, and the pressure inthe Ekman boundary layer as ˜ U ∗ E ( z ∗ ), ˜ V ∗ E ( z ∗ ), and P ∗ E ( x ∗ , y ∗ , z ∗ ), respectively. The governingequations and boundary conditions are − f ∗ ˜ V ∗ E = − ρ ∂P ∗ E ∂x ∗ + 1 ρ F ∗ + ν E ∂ ˜ U ∗ E ∂z ∗ , (18) f ∗ ˜ U ∗ E = − ρ ∂P ∗ E ∂y ∗ + ν E ∂ ˜ V ∗ E ∂z ∗ , (19)0 = − ρ ∂P ∗ E ∂z ∗ − g, (20)˜ U ∗ E = ˜ V ∗ E = 0 , z ∗ = 0 , (21)˜ U ∗ E → ˜ U ∗ , ˜ V ∗ E → ˜ V ∗ , P ∗ E → P ∗ , δ ∗ E ≪ z ∗ ≪ H, (22)5here topography is neglected for simplicity. δ ∗ E is the characteristic depth of the boundarylayer, which is to be determined. ν E is the bulk eddy viscosity in the boundary layer. From theequation (20), we have ∂∂z ∗ ( ∂P ∗ E ∂x ∗ ) = ∂∂z ∗ ( ∂P ∗ E ∂y ∗ ) = 0 . As a result, − ρ ∂P ∗ E ∂x ∗ + 1 ρ F ∗ = − ρ ∂P ∗ E ∂x ∗ | z ∗ ≫ δ ∗ E + 1 ρ F ∗ = − ρ ∂P ∗ ∂x ∗ + 1 ρ F ∗ ≈ − f ∗ ˜ V ∗ , (23) − ρ ∂P ∗ E ∂y ∗ = − ρ ∂P ∗ E ∂y ∗ | z ∗ ≫ δ ∗ E = − ρ ∂P ∗ ∂y ∗ ≈ f ∗ ˜ U ∗ , (24)where the geostrophic balance approximation is used. The geostrophic balance approximationresults from the balance of terms that have magnitude of 10 − m / s in the equations (12)–(14).Substituting (23)–(24) to the equations (18)–(19), we have − f ∗ ( ˜ V ∗ E − ˜ V ∗ ) ≈ ν E ∂ ˜ U ∗ E ∂z ∗ , (25) f ∗ ( ˜ U ∗ E − ˜ U ∗ ) ≈ ν E ∂ ˜ V ∗ E ∂z ∗ . (26)The solution is˜ U ∗ E ≈ ˜ U ∗ − ˜ U ∗ exp( − z ∗ δ ∗ E ) cos( z ∗ δ ∗ E ) − ˜ V ∗ exp( − z ∗ δ ∗ E ) sin( z ∗ δ ∗ E ) , (27)˜ V ∗ E ≈ ˜ V ∗ − ˜ V ∗ exp( − z ∗ δ ∗ E ) cos( z ∗ δ ∗ E ) + ˜ U ∗ exp( − z ∗ δ ∗ E ) sin( z ∗ δ ∗ E ) , (28)where δ ∗ E = (2 ν E /f ∗ ) / [31]. Therefore, the shear stresses exerted by the ground are τ ∗ x = − ρν E ∂ ˜ U ∗ E ∂z ∗ | z ∗ =0 ≈ − ρν E ( ˜ U ∗ − ˜ V ∗ ) δ ∗ E , (29) τ ∗ y = − ρν E ∂ ˜ V ∗ E ∂z ∗ | z ∗ =0 ≈ − ρν E ( ˜ U ∗ + ˜ V ∗ ) δ ∗ E . (30)The original equations (15)–(17) are modified as ∂ ˜ U ∗ ∂x ∗ + ∂ ˜ V ∗ ∂y ∗ = 0 , (31) ∂ ˜ U ∗ ∂t ∗ + ˜ U ∗ ∂ ˜ U ∗ ∂x ∗ + V ∗ ∂ ˜ U ∗ ∂y ∗ − ( f ∗ + f ∗ h ∗ H + β ∗ y ∗ ) ˜ V ∗ = − ρ ∂P ∗ ∂x ∗ − ν E ( ˜ U ∗ − ˜ V ∗ ) δ ∗ E H + 1 ρ F ∗ , (32) ∂ ˜ V ∗ ∂t ∗ + ˜ U ∗ ∂ ˜ V ∗ ∂x ∗ + V ∗ ∂ ˜ V ∗ ∂y ∗ + ( f ∗ + f ∗ h ∗ H + β ∗ y ∗ ) ˜ U ∗ = − ρ ∂P ∗ ∂y ∗ − ν E ( ˜ U ∗ + ˜ V ∗ ) δ ∗ E H , (33)6o that the volume integrals in the computational domain of the body forces τ ∗ x /H and τ ∗ y /H areequal to the surface integrals of the shear stresses on the ground.Introducing ˜ P ∗ = P ∗ − ρ ( f ∗ + ν E / ( δ ∗ E H )) ˜ ψ ∗ , where ˜ ψ ∗ is the stream function satisfying˜ U ∗ = − ∂ ˜ ψ ∗ /∂y ∗ and ˜ V ∗ = ∂ ˜ ψ ∗ /∂x ∗ , we have ∂ ˜ U ∗ ∂x ∗ + ∂ ˜ V ∗ ∂y ∗ = 0 , (34) ∂ ˜ U ∗ ∂t ∗ + ˜ U ∗ ∂ ˜ U ∗ ∂x ∗ + V ∗ ∂ ˜ U ∗ ∂y ∗ − ( f ∗ h ∗ H + β ∗ y ∗ ) ˜ V ∗ = − ρ ∂ ˜ P ∗ ∂x ∗ − ν E δ ∗ E H ˜ U ∗ + 1 ρ F ∗ , (35) ∂ ˜ V ∗ ∂t ∗ + ˜ U ∗ ∂ ˜ V ∗ ∂x ∗ + V ∗ ∂ ˜ V ∗ ∂y ∗ + ( f ∗ h ∗ H + β ∗ y ∗ ) ˜ U ∗ = − ρ ∂ ˜ P ∗ ∂y ∗ − ν E δ ∗ E H ˜ V ∗ . (36)Define the non-dimensional coordinates, time, velocity, stream function, pressure, topography,and zonal force as in [11]:( x, y ) = ( x ∗ , y ∗ ) L , t = f ∗ t ∗ , ( U, V ) = ( ˜ U ∗ , ˜ V ∗ ) f ∗ L , ψ = ˜ ψ ∗ f ∗ L , P = ˜ P ∗ ρf ∗ L ,h = h ∗ H , F = F ∗ ρf ∗ L .
Define non-dimensional parameters β and k as β = Lf ∗ β ∗ = 2 L a tan φ − πL , k = δ ∗ E H .
The zonal force F is further decomposed as the sum of the average zonal force F ave and thefluctuating zonal force F ′ , which are defined as F ave = 12 π Z π Z π F d x d y and F ′ = F − F ave .Then we have the non-dimensional governing equations and boundary conditions of the prim-itive variables: ∂U∂x + ∂V∂y = 0 , (37) ∂U∂t + U ∂U∂x + V ∂U∂y − ( h + βy ) V = − ∂P∂x − kU + F ave + F ′ , (38) ∂V∂t + U ∂V∂x + V ∂V∂y + ( h + βy ) U = − ∂P∂y − kV, (39)( U, V, P )( x + 2 π, y, t ) = ( U, V, P )( x, y, t ) , V | y =0 = V | y = π = 0 . (40)Introducing the vertical component of the non-dimensional vorticity ζ = ∇ ψ , where ∇ = ∂ /∂x + ∂ /∂y , we have the governing equation of the vertical vorticity [11]: ∂ζ∂t + ∂ψ∂x ∂ ( ζ + h + βy ) ∂y − ∂ψ∂y ∂ ( ζ + h + βy ) ∂x = − kζ − ∂F ′ ∂y , (41)7nd the boundary conditions of the stream function [11]: ψ ( x + 2 π, y, t ) = ψ ( x, y, t ) , ∂ψ∂x (cid:12)(cid:12)(cid:12)(cid:12) y =0 = ∂ψ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) y = π = 0 . (42)The equations of the primitive variables (37)–(40), however, cannot be recovered from theequations of the vertical vorticity (41)–(42) without the information of the average zonal force F ave . Therefore, the equations of the vertical vorticity need an additional condition on the averagezonal force to be well-posed.This additional condition can also be given in the form with respect to the average velocity U ave = 12 π Z π Z π U d x d y = 1 π Z π U d y. (43)Observing that Z π V dx = Z π ∂ψ∂x dx = 0due to the periodic boundary condition (42), we notice from the equation (38) that the averagezonal velocity and the average zonal force are coupled throughd U ave d t = ( hV ) ave − kU ave + F ave , (44)where( hV ) ave = 12 π Z Z hV d x d y is the average topographic drag [15]. If the average zonal velocity U ave is prescribed to be aconstant, the average zonal force will be F ave = kU ave − ( hV ) ave , which may vary with time. Therefore, a specified average zonal velocity gives the average zonalforce implicitly, making the equations of the vertical vorticity well-posed.The additional condition can even be given as aF ave + bU ave + c = 0 , where a , b , and c are constants, as in the study of the plane Poiseuille flow [33]. This additionalcondition includes the constant F ave condition and the constant U ave condition as special cases,and also permits a class of more general conditions F ave = − ba U ave − ca a, b = 0. In Section 6, we will show that Charney and Devore solved the equations of thevertical vorticity (41)–(42) under this type of additional condition, although they did not realiseit. The additional condition F ave = − ba U ave − ca seems to lack a clear physical meaning, but it is useful in the calculation of multiple equilibriumstates that are unstable in both the constant F ave calculation and the constant U ave calculation.In addition, the ill-posedness of the equations of the vertical vorticity can also be revealedwith the uniqueness of the solution of Poisson’s equation. The equations of the vertical vorticity ζ depend on the stream function ψ , so we need to solve the Poisson’s equation ∇ ψ = ζ in each time step, which has a unique solution ψ under the boundary conditions ψ | x =0 = ψ | x =2 π and ψ | y =0 = ψ | y = π = 0 . Therefore, the general solution of the Poisson’s equation under the boundary condition (42) is ψ = ψ + C + C y , where C and C are constants to be determined. The constant C can bechosen arbitrarily, but C y contributes to the term ∂ψ/∂y in the equation (41) and should bespecified. Specifying C is equivalent to specifying the average zonal velocity U ave because C = 1 π ( ψ | y = π − ψ | y =0 ) = − U ave . Therefore, an additional condition on the average zonal velocity guarantees a unique solutionof Poisson’s equation of the stream function, and makes the equations of the vertical vorticitywell-posed.In this paper, we concentrate on a simple case with F ′ = kC (cid:18) sin y − π (cid:19) , (45)where C is a constant as in [11]. Cases with constant average zonal forces and cases with constantaverage zonal velocities are both considered in the following.
3. Nonlinear stability
In this section, we will prove that there is an asymptotically stable solution to the equations(37)–(40) under the condition of constant average zonal force or constant average zonal velocityfor F ′ given by (45) and h = Cη cos y , where C and η > − F ave is a constant, and F ′ and h are functions of only themeridional coordinate y , the basic flow of the equations (37)–(40) is a parallel flow: U = 1 k ( F ave + F ′ ) , V = 0 , P = − k Z y ( h ( y ′ ) + βy ′ )( F ave + F ′ ( y ′ )) d y ′ + P , (46)9here P is any constant.From the equations (37)–(40), the governing equations and boundary conditions of the dis-turbance ( u, v, p ) = ( U, V, P ) − ( U , , P ) are ∂u∂x + ∂v∂y = 0 , (47) ∂u∂t + U ∂u∂x + v d U d y + u ∂u∂x + v ∂u∂y − ( h + βy ) v = − ∂p∂x − ku, (48) ∂v∂t + U ∂v∂x + u ∂v∂x + v ∂v∂y + ( h + βy ) u = − ∂p∂y − kv, (49)( u, v, p )( x + 2 π, y, t ) = ( u, v, p )( x, y, t ) , v | y =0 = v | y = π = 0 , (50)where ( u, v ) are the velocity components of the disturbance, and p is the pressure perturbation.Introducing the disturbance vertical vorticity ζ ′ = ∂v/∂x − ∂u/∂y , we have ∂ζ ′ ∂t + U ∂ζ ′ ∂x + u ∂ζ ′ ∂x + v ∂ζ ′ ∂y + ( d h d y − d U d y + β ) v = − kζ ′ . Then we havedd t Z π Z π u + v x d y = − Z π Z π uv d U d y d x d y − k Z π Z π ( u + v ) d x d y, (51)dd t Z π Z π ζ ′ x d y = − Z π Z π uv ( d h d y − d U d y ) d x d y − k Z π Z π ζ ′ d x d y, (52)where integration by parts and the boundary condition (50) are used.If we assume F ′ = kC (sin y − /π ) and F ave = k ( C ′ + 2 C/π ), where C and C ′ are constants,then we have U = C sin y + C ′ according to (46). When h = Cη cos y , where η is a constant,(51)–(52) lead todd t Z π Z π ( ζ ′ + ( η − u + v )) d x d y = − k Z π Z π ( ζ ′ + ( η − u + v )) d x d y. (53)which is an extension of the remarkable identity in the case of flat topography [34].When η >
1, (53) implies that the disturbance kinetic energy 0 . RR ( u + v ) d x d y willexponentially decay. However, we are more interested in the stability of the basic flow when thetopography is flat ( η = 0). It is natural to ask when is the integral R π R π ( ζ ′ + ( η − u + v )) d x d y positive definite. Actually, we have the following lemma: Lemma 3.1. If u, v ∈ H ((0 , π ) × (0 , π )) are real functions with period π in the zonal direction,and satisfy ∂u∂x + ∂v∂y = 0 , (54) v | y =0 = v | y = π = 0 , (55) Z π u ( x,
0) d x = Z π u ( x, π ) d x, Z π Z π u d x d y = 0 , (56)10 hen Z π Z π ζ ′ d x d y ≥ Z π Z π ( u + v ) d x d y. (57)To prove Lemma 3.1, we need the following two versions of Wirtinger’s inequality [35]: Lemma 3.2. If f ∈ H (0 , L y ) , f (0) = f ( L y ) , and R L y f d y = 0 , then Z L y ( d f d y ) d y ≥ ( 2 πL y ) Z L y f d y. Lemma 3.3. If f ∈ H (0 , L y ) and f (0) = f ( L y ) = 0 , then Z L y ( d f d y ) d y ≥ ( πL y ) Z L y f d y. Proof of Lemma 3.1.
The Fourier expansions of u , v , and ζ ′ are u ( x, y ) = + ∞ X m = −∞ ˆ u m ( y )e i mx , v ( x, y ) = + ∞ X m = −∞ ˆ v m ( y )e i mx , ζ ′ ( x, y ) = + ∞ X m = −∞ (i m ˆ v m − dˆ u m d y )e i mx , where ˆ u − m = ¯ˆ u m and ˆ v − m = ¯ˆ v m , and the overlines denote complex conjugations. ˆ u , ˆ v ∈ H (0 , π )are real functions. Substituting the Fourier expansions into (54)–(56), we havei m ˆ u m + dˆ v m d y = 0 , ( m ≥ , (58)ˆ v m (0) = ˆ v m ( π ) = 0 , ( m ≥ , (59)ˆ u (0) = ˆ u ( π ) , Z π ˆ u d y = 0 . (60)Then we have ˆ v = 0, ˆ u m ∈ H (0 , π ), and ˆ v m ∈ H (0 , π ) for m ≥
1, and12 π Z π Z π ( u + v ) d x d y = Z π ˆ u d y + X m ≥ Z π (cid:16) m (cid:12)(cid:12)(cid:12) dˆ v m d y (cid:12)(cid:12)(cid:12) + 2 | ˆ v m | (cid:17) d y, (61)12 π Z π Z π ζ ′ d x d y = Z π ( dˆ u d y ) d y + X m ≥ Z π (cid:16) m (cid:12)(cid:12)(cid:12) d ˆ v m d y (cid:12)(cid:12)(cid:12) + 4 (cid:12)(cid:12)(cid:12) dˆ v m d y (cid:12)(cid:12)(cid:12) + 2 m | ˆ v m | (cid:17) d y. (62)Using Lemma 3.2 for ˆ u , we have Z π ( dˆ u d y ) d y ≥ Z π ˆ u d y. (63)11sing Lemma 3.3 for ˆ v m , we have Z π (cid:16) m (cid:12)(cid:12)(cid:12) d ˆ v m d y (cid:12)(cid:12)(cid:12) + 4 (cid:12)(cid:12)(cid:12) dˆ v m d y (cid:12)(cid:12)(cid:12) + 2 m | ˆ v m | (cid:17) d y ≥ Z π (cid:16) (cid:12)(cid:12)(cid:12) dˆ v m d y (cid:12)(cid:12)(cid:12) + ( 2 m + 2 m ) | ˆ v m | (cid:17) d y ≥ Z π (cid:16) m (cid:12)(cid:12)(cid:12) dˆ v m d y (cid:12)(cid:12)(cid:12) + 2 | ˆ v m | (cid:17) d y (64)for m ≥
1. Substituting (63) and (64) into (61) and (62), we arrive at the inequality (57), whichcompletes the proof.
Remark 3.4.
The conclusion of Lemma 3.1 is optimal. When ( u, v ) = (sin x cos y, − cos x sin y ), Z π Z π ζ ′ d x d y = 2 Z π Z π ( u + v ) d x d y. Under the conditions (54)–(56), the integral R π R π ( ζ ′ + ( η − u + v )) d x d y is positivedefinite for η > −
1. A solution to the equations of the disturbance (47)–(50), however, does nothave to satisfy the condition (56) for all t ≥
0. According to the equation (48), we havedd t Z π ( u ( x, − u ( x, π )) d x = − k Z π ( u ( x, − u ( x, π )) d x, (65)dd t Z π Z π u d x d y = − k Z π Z π u d x d y. (66)Although the integrals R π ( u ( x, − u ( x, π )) d x and R π R π u d x d y will not vanish if they do notvanish at t = 0, they will both decay exponentially with time as exp( − kt ). Their influence onthe conclusion of Lemma 3.1 is characterised by the following lemma: Lemma 3.5. If u, v ∈ H ((0 , π ) × (0 , π )) are real functions with period π in the zonal direction,and satisfy (54) and (55), then Z π Z π ζ ′ d x d y ≥ Z π Z π ( u + v ) d x d y − π (cid:16) Z π Z π u d x d y (cid:17) − π − π (cid:16) Z π ( u ( x, − u ( x, π )) d x (cid:17) . (67) Proof.
Performing the same Fourier expansions as in the proof of Lemma 3.1, we still have (61),(62), and the inequality (64). Using Lemma 3.2 for f ( y ) = ˆ u ( y ) − π Z π ˆ u ( y ′ ) d y ′ + ( yπ −
12 )(ˆ u (0) − ˆ u ( π )) , we have Z π ( dˆ u d y ) d y − π (ˆ u (0) − ˆ u ( π )) ≥ Z π ˆ u d y + 8(ˆ u (0) − ˆ u ( π )) Z π ( yπ −
12 )ˆ u d y − π (cid:16) Z π ˆ u d y (cid:17) + π u (0) − ˆ u ( π )) ≥ Z π ˆ u d y − π (cid:16) Z π ˆ u d y (cid:17) − π u (0) − ˆ u ( π )) , Z π ( dˆ u d y ) d y ≥ Z π ˆ u d y − π (cid:16) Z π Z π u d x d y (cid:17) − π − π (cid:16) Z π ( u ( x, − u ( x, π )) d x (cid:17) . (68)Substituting (68) and (64) into (61) and (62), we obtain the inequality (67), which completesthe proof. Theorem 3.6. If F ′ = kC (sin y − /π ) , F ave = k ( C ′ + 2 C/π ) , and h = Cη cos y , where C , C ′ ,and η are constants, and η > − , the basic flow U = ( C sin y + C ′ , is stable in the sense that k U ( t ) − U k H ≤ (cid:16)
31 + η ′ C + 2(2 − η ′ ) π (1 + η ′ ) C + ( π − − η ′ )6 π (1 + η ′ ) C (cid:17) / exp( − kt ) , (69) where η ′ = min { η, } ,C = (cid:16) Z π Z π ( ζ ′ + ( η − u + v )) d x d y (cid:17) t =0 ,C = (cid:16) Z π Z π u d x d y (cid:17) t =0 ,C = (cid:16) Z π ( u ( x, − u ( x, π )) d x (cid:17) t =0 . Proof.
When − < η <
2, we have the following estimate of the H norm of u = ( u, v ): k u k H = Z π Z π ( ζ ′ + u + v ) d x d y ≤
31 + η Z π Z π ( ζ ′ + ( η − u + v )) d x d y + 2(2 − η ) π (1 + η ) (cid:16) Z π Z π u d x d y (cid:17) + ( π − − η )6 π (1 + η ) (cid:16) Z π ( u ( x, − u ( x, π )) d x (cid:17) . When η ≥
2, it is obvious that k u k H ≤ Z π Z π ( ζ ′ + ( η − u + v )) d x d y. Furthermore, noticing (53), (65), and (66), we finish the proof of the theorem.In [30], we proved a weaker version of Theorem 3.6, which stated that k ζ ′ k L would decayexponentially with time provided that η ≥
0. But the proof therein is not strict. The governingequations in [30] need an additional condition on the average zonal force or the average zonalvelocity to be well-posed, as in the work of Charney and DeVore [11].Now we consider the case with a constant average zonal velocity instead of a constant averagezonal force. When the average zonal velocity U ave is a constant, and h is a function of only themeridional coordinate y , the average zonal force F ave = kU ave − ( hV ) ave = kU ave is also a constantaccording to the equation (44), where ( hV ) ave = ( ∂ ( hψ ) /∂x ) ave = 0 is used. Therefore we have13 heorem 3.7. If F ′ = kC (sin y − /π ) , U ave = C ′ + 2 C/π , and h = Cη cos y , where C , C ′ , and η are constants, and η > − , the basic flow U = ( C sin y + C ′ , is stable in the sense that theinequality (69) holds.
4. Numerical methods
The periodic boundary conditions for U , V , and P in the zonal direction permit Fourierexpansions of these primitive variables. Denoting the Fourier expansion in the zonal direction ofa real function f ( x, y, t ) as P m ˆ f m ( y, t )e i mx or P m ( F f ) m e i mx , where ˆ f ∈ R and ˆ f − m = ¯ˆ f m , with¯ˆ f m representing the complex conjugation of ˆ f m . The equations of primitive variables (37)–(40)with constant average zonal force F ave lead toi m ˆ U m + ∂ ˆ V m ∂y = 0 , (70) ∂ ˆ U m ∂t + F ( U ∂U∂x + V ∂U∂y ) m − F ( hV ) m − βy ˆ V m = − i m ˆ P m − k ˆ U m + δ m, F ave + ˆ F ′ m , (71) ∂ ˆ V m ∂t + F ( U ∂V∂x + V ∂V∂y ) m + F ( hU ) m + βy ˆ U m = − ∂ ˆ P m ∂y − k ˆ V m , (72)ˆ V m | y =0 = ˆ V m | y = π = 0 , (73)where δ m, = 1 for m = 0, and δ m, = 0 for 1 ≤ m ≤ M − M is the number of independentFourier modes used in the discretisation. Eliminating ˆ P m from the above equations, we have ∂ ˆ U ∂t + F ( U ∂U∂x + V ∂U∂y ) − F ( hV ) = − k ˆ U + F ave + ˆ F ′ , (74)ˆ V = 0 , (75)and ˆ U m = i m ∂ ˆ V m ∂y , (76) ∂∂t ( ∂ ∂y − m ) ˆ V m + i m F ( U ∇ V − V ∇ U ) m + i m F ( U ∂h∂x + V ∂h∂y ) m + i mβ ˆ V m = − k ( ∂ ∂y − m ) ˆ V m − i m ∂ ˆ F ′ m ∂y , (77)ˆ V m | y =0 = ˆ V m | y = π = 0 , (78)for 1 ≤ m ≤ M − U ave is prescribed as a constant, the average zonal force will be F ave = kU ave − ( hV ) ave according to the equation (44). Instead of the equation (74), we have ∂ ˆ U ∂t + F ( U ∂U∂x + V ∂U∂y ) − F ( hV ) = − k ˆ U + kU ave − ( hV ) ave + ˆ F ′ . (79)14f the additional condition is aF ave + bU ave + c = 0, where a , b , and c are constants, and a, b = 0, we have the following equation instead of the equation (74): ∂ ˆ U ∂t + F ( U ∂U∂x + V ∂U∂y ) − F ( hV ) = − k ˆ U − ba U ave − ca + ˆ F ′ . (80)The Chebyshev-Tau method with N Chebyshev polynomials is used for the discretisation inthe meridional direction. In the temporal discretisation, linear terms except terms containing h are discretised implicitly with the second-order Crank-Nicolson scheme, and other terms arediscretised explicitly with the four-stage third-order Runge-Kutta scheme [36]. Direct numerical simulations are difficult to obtain linearly unstable equilibrium states, be-cause any small-amplitude disturbances will lead the solutions away from the unstable equilib-rium states in the phase space. Therefore, it is more adequate to solve the nonlinear steadyequations directly, instead of calculating equilibrium states as final states after a sufficientlylong time in direct numerical simulations. The nonlinear steady equations are solved with apseudo-arclength continuation method.Removing the time derivative terms in (74)–(78), expanding the nonlinear terms, and elimi-nating ˆ U m , we have − X m ≥ ( i m ¯ˆ V m D ˆ V m − i m ˆ V m D ¯ˆ V m ) + X m ≥ (¯ˆ h m ˆ V m + ˆ h m ¯ˆ V m ) − k ˆ U + F ave + ˆ F ′ = 0 , (81)and 0 = X m + m = mm =0 − mm (D ˆ V m )(D − m ) ˆ V m + mm ˆ V m (D − m )D ˆ V m − i mm m (D ˆ V m )ˆ h m + i m ˆ V m Dˆ h m ! + (i m ˆ U + k )(D − m ) ˆ V m − i m (D ˆ U − β ) ˆ V m − m ˆ U ˆ h m + i m D ˆ F ′ m (82)for 1 ≤ m ≤ M −
1, where D ≡ d / d y .Using the Chebyshev-Tau method, we expand ˆ U and ˆ V m with N Chebyshev polynomials,and denote the column vectors composed of the coefficients as ˆ U and ˆ V m , respectively. Thenthe equations (81) and (82) with the boundary condition (78) constitute a nonlinear equation G ( f , F ave ) = , where f = ( ˆ U T0 Re ˆ V T1 Im ˆ V T1 · · · Re ˆ V T M − Im ˆ V T M − ) T . The superscript “T” denotes the transpose of a vector or a matrix. “Re” and “Im” represent thereal part and the imaginary part of a complex vector, respectively.We use the following predictor-corrector continuation method [37–39] to find a branch ofsolutions of G ( f , F ave ) = . 15 .2.1. Predictor (the first solution) The initial guess of the first solution is chosen as an approximation solution of the equations(81), (82), and (78) at a sufficiently large F ave . For example, if h = h cos 2 x sin y , then anapproximation solution is F ave = 0 . U = ( F ave + ˆ F ′ ) /k , ˆ V = 0 . h sin y , and ˆ V m = 0 for1 ≤ m ≤ M − m = 2. Denote this approximation solution as ( f , F ave , ) or ( f , F ). Thesubscript “ave” will be omitted in the following for conciseness. The iterative scheme for the calculation of the first solution is ∂ G ∂ f ( f n − , F )∆ f n = − G ( f n − , F ) , f n = f n − + ∆ f n , for n ≥
1. If the sequence { f n } n ≥ converges to f , then the first solution is defined as ( f , F ).Otherwise, a more accurate approximation solution is needed in Section 4.2.1. For example, wecan choose the same form of the approximation solution, but at a larger F ave . ν -th solution, ν ≥ ) After the solution ( f ν − , F ν − ) is obtained, we search the next solution ( f ν , F ν ) which satisfies || ( f ν − f ν − , F ν − F ν − ) || = ∆ s , where || ( f , F ) || = f T W f + F / (2 k ). W is a weight matrixso that f T W f is the discretisation of the average kinetic energy E ave = 12 π Z π Z π U + V x d y. (83)An approximation of the ν -th solution can be obtained by a tangent predictor [38]: ∂ G ∂ f ( f ν − , F ν − )∆ ˜ f ν = − ∂ G ∂F ( f ν − , F ν − ) , (∆ f ν , ∆ F ν ) = ± (∆ ˜ f ν , || (∆ ˜ f ν , || ∆ s, ( f ν , F ν ) = ( f ν − , F ν − ) + (∆ f ν , ∆ F ν ) . The plus or minus sign is determined by demanding the angle between (∆ f ν , ∆ F ν ) and ( f ν − − f ν − , F ν − − F ν − ) is less than 90 ◦ when ν ≥
3, where the angle is calculated with the innerproductcorresponding to the norm || · || . When ν = 2, we choose the (∆ f , ∆ F ) with a negative ∆ F to calculate the solutions ( f, F ) with F < F . The other half of the branch of solutions with F > F can also be obtained by choosing a positive ∆ F in a new calculation. ν -th solution, ν ≥ ) The iterative scheme for the calculation of the ν -th solution ( f ν , F ν ) ( ν ≥
2) is (cid:18) ∂ G ∂ f ( f nν , F nν ) ∂ G ∂F ( f nν , F nν )2( f nν − f ν − ) T W ( F nν − F ν − ) /k (cid:19) (cid:18) ∆ f n +1 ν ∆ F n +1 ν (cid:19) = (cid:18) − G ( f nν , F nν ) − A ν ( f nν , F nν ) (cid:19) , ( f n +1 ν , F n +1 ν ) = ( f nν , F nν ) + (∆ f n +1 ν , ∆ F n +1 ν ) , n ≥
1, where A ν ( f , F ) = || ( f − f ν − , F − F ν − ) || − ∆ s . This scheme is inspired by the Taylor expansions of G ( f n +1 ν , F n +1 ν ) = and A ν ( f n +1 ν , F n +1 ν ) = 0at ( f nν , F nν ).If the sequence { ( f nν , F nν ) } n ≥ converges, the ν -th solution ( f ν , F ν ) is defined as the limit ofthe sequence. Otherwise, we change the step from ∆ s = ∆ s to ∆ s = 0 . s , and recalculatethe ν -th solution.The branch of solution ends where new solutions of G ( f , F ave ) = cannot be obtained nomatter how small the step is, which means either there are no more equilibrium states or theequilibrium states need a better resolution (larger M and N ) to be calculated.
5. Results
In Section 3, we have proved the asymptotic stability of the basic flow for a zonal-invarianttopography h = Cη cos y in a region of parameter space. In the following, we perform numericalsimulations for a more realistic topography h = h cos 2 x sin y , where h is a constant. Thistopography was also considered by Charney and Devore [11]. When h = h cos 2 x sin y and F ′ = √ kψ A0 (sin y − /π ), the equations (37)–(40) can be solved in the following three spaces: S = { ( U, V ) ∈ H × H | ∂U∂x + ∂V∂y = 0 , V ( x,
0) = V ( x, π ) = 0 } , S even = { ( U, V ) ∈ S | U ( x + π, y ) = U ( x, y ) , V ( x + π, y ) = V ( x, y ) } , S sym = { ( U, V ) ∈ S even | U ( x + π , π − y ) = U ( x, y ) , V ( x + π , π − y ) = − V ( x, y ) } , where H represents the space of real functions in H ((0 , π ) × (0 , π )) with period 2 π in thezonal direction.The equilibrium states of the equations (37)–(40) are calculated with the pseudo-arclengthcontinuation method in the space S sym for k = 0 . β = 0 .
25, and ψ A0 = 0 .
2. We use 128 evenor odd Chebyshev polynomials in the meridional direction, and use 16 even Fourier modes in thezonal direction. It is equivalent to setting N = 256 and M = 32 in Section 4.2 because of thesymmetry of the space S sym . The maximum step in the continuation method is 0.01. h = 0 . h = 0 .
2, there are multiple equilibrium states with different U ave corresponding tothe same F ave (Figure 1). For example, when F ave = 0 . U ave = 0 . . . ψ = 0, − . πU ave , − . πU ave , − . πU ave , and − πU ave from bottom to top. In this paper, we alwayschoose ψ | y =0 = 0, and then we have ψ | y = π = − πU ave . The thin black lines are the contours ofthe topography h = h cos 2 x sin y for reference, with dashed lines representing negative values.There is a small-scale circulation in the region circled by the closed streamline ψ = 0 inthe dashed rectangle of Figure 2(a), which can be seen clearly in Figure 2(b). The stream17 ave F a v e Continuation method(Charney & DeVore 1979) F ave =0.000747 U ave +0.00167 Fig. 1.
The average zonal force and the average zonal velocity of the equilibrium states for h = 0 .
2. The solidline is composed of the equilibrium states calculated with the continuation method. The three equilibrium statescalculated by Charney and DeVore [11] are shown by squares (Table 1), and satisfy F ave = 0 . U ave + 0 . k = 0 .
01 and ψ A0 = 0 . x / π y / π x / π y / π x / π y / π x / π y / π (c)(a) (b)(d) Fig. 2.
The contours of the stream function ψ (thick coloured lines) and the topography h (thin black lines withdashed lines representing negative values) of the three equilibrium states for h = 0 . F ave = 0 . U ave = 0 .
170 (a), U ave = 0 . U ave = 0 . ψ = 0, − . πU ave , − . πU ave , − . πU ave , and − πU ave from bottom to top. (b) is a close-up view of the dashed rectangle in (a) with velocity vectors plotted ona coarse mesh in the region ψ > − . y = 0. ψ max at the stagnation point ( x, y ) = (0 . π, . π ) in Figure2(b). Therefore, the circumferential flux of the small-scale circulation in the region ψ > φ c = ψ max −
0, where a positive flux represents a clockwise circulation.Actually, four regions contain small-scale circulations in Figure 2(a). Two of them attachingto the boundary y = 0 have the same circumferential flux φ c , and other two regions attachingto the boundary y = π have the same circumferential flux − φ c , because of the symmetry of thevelocity field in the space S sym . The regions of small-scale circulations will shrink when U ave is increasing. When U ave ≥ .
23, small-scale circulations no longer exist. The contours of thestream function when U ave = 0 .
22 and 0 .
23 are plotted in Figure 3. In Figure 3(b), the contours ψ = 0 and ψ = − πU ave coincide with the boundary y = 0 and y = π , respectively. x / π y / π (a) (b) x / π y / π Fig. 3.
The contours of the stream function ψ (thick coloured lines) and the topography h (thin black lines withdashed lines representing negative values) for ( F ave , U ave ) = (0 . , . F ave , U ave ) = (0 . , . h = 0 .
2. The thick coloured contours in each subfigure correspond to ψ = 0, − . πU ave , − . πU ave , − . πU ave , − . πU ave , − . πU ave , and − πU ave from bottom to top. In Charney and DeVore’s numerical simulation, the low-index and the high-index equilibriumstates are stable, whereas the medium-index equilibrium state is unstable [11]. We examine thestabilities of some equilibrium states with direct numerical simulations. The velocity fields ofthe equilibrium states calculated with the pseudo-arclength continuation method are used asthe initial velocity fields in direct numerical simulations. We use 192 even or odd Chebyshevpolynomials and 192 even Fourier modes in the meridional direction and the zonal direction,respectively. It is equivalent to setting N = M = 384 in Section 4.1 because of the symmetry ofthe space S sym . The time step is ∆ t = 0 . F ave = 0 . F ave , U ave ) = (0 . , . S sym , and will evolve into the high-index equilibrium state ( F ave , U ave ) = (0 . , . F = √ kψ A0 sin y and required the “net pressure gradient” to be zero, which was equivalentto setting constant average zonal force to be F ave = 2 √ kψ A0 /π and F ′ = √ kψ A0 (sin y − /π ). Using the same parameters k = 0 . β = 0 . h = 0 .
2, and ψ A0 = 0 .
2, we obtainthree equilibrium states for F ave = 2 √ kψ A0 /π ≈ . U ave = 0 .
127 and 0 . U ave = 0 . ψ = − . − . a v e
400 800 120000.030.060.090.120.150.180.210.24 U ave (0) = 0.170 U ave (0) = 0.0915 U ave (0) = 0.0336 Fig. 4.
The evolution of U ave from the three equilibrium states with ( F ave , U ave ) = (0 . , . F ave , U ave ) = (0 . , . F ave , U ave ) = (0 . , . F ave is fixed to be 0.002 and h = 0 . − .
25 and − .
35 for the equilibrium state U ave = 0 .
127 in Figure 5(a), and plot the contours of ψ = 0 . − . − .
06 and − . U ave = 0 . x / π y / π (a) (b) x / π y / π Fig. 5.
The contours of the stream function ψ (thick coloured lines) and the topography h (thin black lineswith dashed lines representing negative values) of the “high-index” equilibrium state U ave = 0 .
127 (a) and the“low-index” equilibrium state U ave = 0 . F ave = 0 . h = 0 .
2, using the same contour levels asFig. 6 in [30]. ≤ h ≤ . h = 0), the equilibrium state U = ( F ave + F ′ ) /k , V = 0 isglobally stable for each F ave according to Theorem 3.6. The average zonal velocity is U ave = F ave /k from the equation (44). Actually, the same equilibrium state U = U ave + F ′ /k , V = 0 is alsoglobally stable for each U ave according to Theorem 3.7. So there is only one stable equilibriumstate for each F ave or for each U ave when h = 0. When h = 0 .
01, there is only one equilibriumstate for each F ave as in the case of h = 0 (Figure 6). When h is increased to 0.02, there arethree equilibrium states for each F ave when 0 . ≤ F ave ≤ . h = 0 .
04, the stabilities of the equilibrium states change when F ave is nearits extreme values. When h = 0 . h = 0 .
2, however, the stabilities of the equilibriumstates may be influenced by limit cycles which are not considered in the present continuationmethod. For example, when h = 0 .
1, the equilibrium state with ( F ave , U ave ) = (0 . , . U ave oscillating between 0.0407 and 0.0448 (Figure 7(a)). The projection of the phase spacetrajectory on the average kinetic energy–average zonal velocity ( E ave - U ave ) plane in Figure 7(b)shows a stable limit cycle. The period of the limit cycle is T ≈
56 (Figure 7(c)). The contoursof the stream functions when U ave achieves its minimum and maximum are shown in Figure 7(d)and (e), respectively. The small-scale circulations almost do not move in the zonal direction asa travelling wave, but they intrude further into the zonal flow in the meridional direction when U ave is minimum than when U ave is maximum. U ave F a v e h = 0.2 h = 0.1 h = 0.04 h = 0.02 h = 0.01 h = 0DNS (stable)DNS (unstable) Fig. 6.
The average zonal force and the average zonal velocity of the equilibrium states for 0 ≤ h ≤ .
6. Discussion
Charney and DeVore [11] did not realise that the equations of the vertical vorticity (41)–(42)are ill-posed, and approximated the stream function as ψ ( x, y, t ) = √ ψ A ( t ) cos y + 2 ψ K ( t ) cos 2 x sin y + 2 ψ L ( t ) sin 2 x sin y, (84)21 a v e E ave U a v e t < 10001000 < t < 15003500 < t < 4000 U ave E ave x / π y / π x / π y / π (d) (e)(b) (c)(a) Fig. 7.
The evolution of the equilibrium state with ( F ave , U ave ) = (0 . , . F ave is fixed to be 0.0022 and h = 0 .
1. (a) Evolution of U ave . (b) Phase space trajectory(projection on the E ave - U ave plane). (c) Evolution of U ave and E ave (a close-up view). The contours of the streamfunction ψ (thick coloured lines) and the topography h (thin black lines with dashed lines representing negativevalues) when t = 4021 (d) and t = 4049 (e). The thick coloured contours in each subfigure correspond to ψ = 0, − . πU ave ( t ), − . πU ave ( t ), − . πU ave ( t ), and − πU ave ( t ) from bottom to top. h = h cos 2 x sin y and F ′ = √ kψ A0 (sin y − /π ), substituting (84) into (41), and using the Galerkin method, they obtainedd ψ A d t = 8 √ π h ψ L − k ( ψ A − ψ A0 ) , (85)d ψ K d t = − √ π ψ A ψ L + 25 βψ L − kψ K , (86)d ψ L d t = 64 √ π ψ A ψ K − √ π h ψ A − βψ K − kψ L . (87)Table 1 lists the three equilibrium states of the equations of the vertical vorticity calculatedby Charney and DeVore [11] for k = 0 . β = 0 . h = 0 .
2, and ψ A0 = 0 .
2. The average zonalvelocity is U ave = 1 π ( ψ | y =0 − ψ | y = π ) = 2 √ π ψ A , (88)and the average zonal force is F ave = d U ave d t + kU ave − ( hV ) ave = 2 √ π ( d ψ A d t + kψ A ) − h ψ L = ( 323 π − h ψ L + 2 √ π kψ A0 , (89)where (44) and (85) are used. The three equilibrium states calculated by Charney and DeVore[11] have neither the same average zonal velocity nor the same average zonal force (Table 1).According to the equation (85), these equilibrium states satisfy 8 √ h ψ L / (3 π ) − k ( ψ A − ψ A0 ) = 0,which is equivalent to F ave = (1 − π
32 ) kU ave + 3 π √ kψ A0 , (90)where (88) and (89) are used. Table 1
The equilibrium states of the equations (85)–(87) for k = 0 . β = 0 . h = 0 .
2, and ψ A0 = 0 . U ave and F ave arecalculated according to (88) and (89), respectively. ψ A ψ K ψ L U ave F ave . Conclusion In this paper, we show that the equations of the vertical vorticity are not well-posed in thestudy of quasi-geostrophic barotropic flows over topography. There lacks an additional conditionon the average zonal force F ave , which can be given explicitly as a constant or as a function ofthe average zonal velocity U ave , or given implicitly by prescribing U ave to be a constant.We prove that there is an asymptotically stable equilibrium state under the condition ofconstant F ave or constant U ave when F ′ = kC (sin y − /π ) and h = Cη cos y , where C and η > − h = 0, the equilibrium state U = ( F ave + F ′ ) /k , V = 0is asymptotically stable in a constant F ave problem, and the equilibrium state U = U ave + F ′ /k , V = 0 is asymptotically stable in a constant U ave problem, which excludes the existence ofmultiple equilibrium states in the case of flat topography.When the topography is h = h cos 2 x sin y and the fluctuating zonal force is F ′ = √ kψ A0 (sin y − /π ) , we calculate the equilibrium states with a pseudo-arclength continuation method for k = 0 . β = 0 . ψ A0 = 0 .
2, and 0 ≤ h ≤ . F ave appear only for h ≥ .
02. Their stabilities are examined with direct numerical simulations.When h = 0 .
04, along the F ave – U ave curve, the stability of the equilibrium state changes nearthe extreme points of F ave . However, this is not the case for h = 0 . h = 0 .
2, where theremay be Hopf bifurcation. The emergence and evolution of limit cycles deserve further studies.
Acknowledgment . The research was partially supported by NSFC of China (11571240, 11602148).
References [1] J. Pedlosky, Geophysical Fluid Dynamics, 2 ed., Springer, New York, 1987.[2] P. Constantin, A. Majda, E. Tabak, Formation of strong fronts in the 2d quasi-geostrophicthermal active scalar, Nonlinearity 7 (1994) 1495–1533.[3] J. Wu, Quasi-geostrophic-type equations with initial data in morrey spaces, Nonlinearity10 (1997) 1409–1420.[4] P. Constantin, J. Wu, Behavior of solutions of 2d quasi-geostrophic equations, SIAM J.Math. Anal. 30 (1999) 937–948.[5] Q. Chen, C. Miao, Z. Zhang, A new bernstein’s inequality and the 2d dissipative quasi-geostrophic equation, Commun. Math. Phys. 271 (2007) 821–838.[6] A. Kiselev, F. Nazarov, A. Volberg, Global well-posedness for the critical 2d dissipativequasi-geostrophic equation, Inventiones Math. 167 (2007) 445–453.[7] P. Constantin, G. Iyer, J. Wu, Global regularity for a modified critical dissipative quasi-geostrophic equation, Indiana Univ. Math. J. 57 (2008) 2861–2897.248] J. Wu, Global solutions of the 2d dissipative quasi-geostrophic equation in besov spaces,SIAM J. Math. Anal. 36 (2005) 1014–1030.[9] Z. M. Chen, W. G. Price, Stability and instability analyses of the dissipative quasi-geostrophic equation, Nonlinearity 21 (2008) 765–782.[10] Z. M. Chen, Bifurcating steady-state solutions of the dissipative quasi-geostrophic equationin lagrangian formulation, Nonlinearity 29 (2016) 3132–3147.[11] J. G. Charney, J. G. DeVore, Multiple flow equilibria in the atmosphere and blocking, J.Atmos. Sci. 36 (1979) 1205–1216.[12] M. Ghil, S. Childress, Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dy-namo Theory, and Climate Dynamics, Springer, New York, 1987.[13] B. Cushman-Roisin, J.-M. Beckers, Introduction to Geophysical Fluid Dynamics: Physicaland Numerical Aspects, 2 ed., Academic Press, Amsterdam, 2011.[14] F. V. Dolzhansky, Fundamentals of Geophysical Hydrodynamics, 2 ed., Springer, Berlin,2013.[15] M. K. Davey, A quasi-linear theory for rotating flow over topography. part 1. steady ββ