Variational balance models for the three-dimensional Euler-Boussinesq equations with full Coriolis force
aa r X i v : . [ m a t h . A P ] D ec Variational balance models for the three-dimensionalEuler-Boussinesq equations with full Coriolis force
G¨ozde ¨Ozden a,b,c and Marcel Oliver a a School of Engineering and Science, Jacobs University, 28759 Bremen, Germany; b MARUM – Center for Marine Environmental Sciences, Universit¨at Bremen, 28334 Bremen,Germany; c Meteorological Institute, Universit¨at Hamburg, 20146 Hamburg, Germany
ARTICLE HISTORY
Compiled January 1, 2021
ABSTRACT
We derive a semi-geostrophic variational balance model for the three-dimensionalEuler–Boussinesq equations on the non-traditional f -plane under the rigid lid ap-proximation. The model is obtained by a small Rossby number expansion in theHamilton principle, with no other approximations made. In particular, we assumeneither hydrostaticity, nor do we neglect the horizontal components of the Coriolisparameter, i.e., we do not make the so-called “traditional approximation”. The re-sulting balance models has the same structure as the “ L balance model” for theprimitive equations: a kinematic balance relation, the prognostic equation for thethree-dimensional tracer field, and an additional prognostic equation for a scalarfield over the two-dimensional horizontal domain which is linked to the undeter-mined constant of integration in the thermal wind relation. The balance relation iselliptic under the assumption of stable stratification and sufficiently small fluctua-tions in all prognostic fields. KEYWORDS
Boussinesq equation, full Coriolis force, variational asymptotics, balance models
1. Introduction
We describe the derivation of a semi-geostrophic variational balance model for thethree-dimensional Euler–Boussinesq equations on the non-traditional f -plane underthe rigid lid approximation. In non-dimensional variables, the Euler–Boussinesq systemtakes the form ε ( ∂ t u + u · ∇ u ) + Ω × u = − ∇ p − ρ k , (1a) ∇ · u = 0 , (1b) ∂ t ρ + u · ∇ ρ = 0 . (1c)where ε is the Rossby number, assumed small, u is the three-dimensional velocity field, Ω = φ sin φ ≡ cs (2)he full Coriolis vector at constant latitude φ which, without loss of generality, isassumed to lie in the y - z plane, p is the pressure, ρ the density, and k the unitvector in the vertical. For simplicity, we assume periodic boundary conditions in thehorizontal and rigid lid boundary conditions in the vertical on a layer of fluid withconstant depth, D = T × [ −
1, 0] , (3)so that k · u = 0 at z = −
1, 0 . (4)The semi-geostrophic scaling used in (1) makes no smallness assumption for ver-tical length scales. More precisely, the aspect ratio, the ratio of typical horizontaland typical vertical scales is assumed to be of order one. Thus, the model is fullynon-hydrostatic. To be consistent, this necessitates that all components of the Cori-olis vector are considered—the so-called “traditional approximation” where only thevertical component of the Coriolis vector is retained cannot be separated from mak-ing the hydrostatic approximation (e.g. White, 2002). In addition, the Rossby num-ber is taken to scale like the square of the Froude number. For a detailed exposi-tory treatment of scaling assumptions and model hierarchies, we refer the reader toFranzke, Oliver, Rademacher, and Badin (2019).We emphasize that our setting is more general than that of the quasi-geostrophiclimit which is based on the hydrostatic and traditional approximations, a constant ratioof Rossby to Froude numbers, and the assumption of strong density stratification.The quasi-geostrophic equations as the first-order balance model in this limit arewidely used for proof-of-concept studies and can be found in standard textbooks (e.g.Pedlosky, 1987; Vallis, 2006). As we shall see, the assumption of strong stratificationwill re-emerge in our setting as a sufficient condition for solvability of the balancerelation; however, it is not assumed a priori . Our main point is the derivation of abalance relation for non-hydrostatic flow.Our derivation is based on variational asymptotics. Approximations to Hamilton’svariational principle for the equations of rotating fluid flow were pioneered by Salmon(1983, 1985) in the context of the shallow water equations and later, in Salmon (1996),for the primitive equations of the ocean. The basic idea is to consider the “extendedHamilton principle” or “phase space variational principle” and use the leading orderbalance relation (geostrophic balance for the shallow water equations or the thermalwind relation for the primitive equations) to constrain the momentum variables. Thestationary points of the constrained action with respect to variations of the positionvariables give a balance relation that includes the next-order correction to the leadingorder constraint. In principle, this method can be iterated to obtain higher-order bal-ance relations. For example, Allen, Holm, and Newberger (2002) derive second ordermodels and show that the so-called L and L models are numerically well behaved.A second idea, also proposed in Salmon (1985), is the use of a near-identity co-ordinate transformation to bring the resulting variational principle or the equationsof motion into a more convenient form. This transformation may be applied pertur-batively, so that the resulting models coincide only up to terms of a certain order.From an analytic perspective, higher-order terms may matter and it turns out that atransformation to a coordinate system in which the Hamilton equations of motion takecanonical form, the original motivation behind this approach, leads to ill-posedness of2he full system of prognostic equations augmented by the balance relation.An even more general framework is obtained by reversing the steps “constrain” and“transform”. Oliver (2006) noted that is possible to assume an entirely general near-identity change of coordinates and expand the transformed Lagrangian in powers ofthe small parameter. Whenever the transformation is chosen such that the Lagrangianis degenerate to the desired order, the variational principle implies a phase space con-straint. This approach gives rise to a greater variety of candidate models that arenot accessible via Salmon’s approach. In particular, it allows us to retain some im-portant mathematical features such as regularity of potential vorticity inversion. Inthe context of the shallow water equations on the f -plane, it turns out that the L model first proposed by Salmon is already optimal among a larger family of mod-els (Dritschel, Gottwald, & Oliver, 2017). In other configurations, such as when theCoriolis parameter is spatially varying (Oliver & Vasylkevych, 2013), or in the fullythree dimensional situation considered here, the additional flexibility that comes fromputting the change of coordinates first is necessary.In this work, we show that semi-geostrophic variational asymptotics can be donedirectly for the three-dimensional Euler–Boussinesq system with full Coriolis forcewithout any preparatory approximations. Our motivation derives, first, from the roleof the Boussinesq equation as the common parent model of nearly all of geophysicalfluid dynamics. Second, we would like to see how “non-traditional” Coriolis effects,associated with the vertical component of the Coriolis force and horizontal Coriolisforces coming from the vertical velocity, enter the balance dynamics. Non-traditionaleffects are significant in a variety of circumstances (Gerkema & Shrira, 2005a, 2005b;Gerkema, Zimmerman, Maas, & van Haren, 2008) and are also studied in the contextof other reduced models such as layers of shallow water (Stewart & Dellar, 2010, 2012)and models on the β -plane (Dellar, 2011).Our approach is similar, in principle, to the primitive equation case(Oliver & Vasylkevych, 2016), and we find that the resulting balance models have thesame structure: a kinematic balance relation which is elliptic under suitable assump-tions, the prognostic equation for the three-dimensional tracer field, and an additionalprognostic equation for a scalar field over the two-dimensional horizontal domain whichis linked to the undetermined constant of integration in the thermal wind relation. Thisfield, in the present setting, takes the form of a skewed relative vorticity averaged alongthe axis of rotation.On the technical level, the computations are considerably more difficult than forthe primitive equation. Some expressions, such as the thermal wind relation, take anatural form in an oblique coordinate system where the axis of rotation takes the roleof the “vertical”. Other expressions, most notably the top and bottom boundary con-ditions and the gravitational force, are most easily described in the usual Cartesiancoordinates. This incompatibility requires a detailed study of averaging and decompo-sition of vector fields in oblique coordinates, and of the translations between obliqueand Cartesian coordinates.The remainder of the paper is organized as follows. Section 2 recalls the variationalderivation for the Euler–Boussinesq equations in the language of the Euler–Poincar´evariational principle. Section 3 explains the general setting for variational asymptoticsin this framework. In Section 4, we develop a variety of technical tools, starting fromthe mapping to oblique coordinates, averaging along the axis of rotation, and thedecomposition of vector field in oblique coordinates. Section 5 discusses the leadingorder balance, the thermal wind relation, on the mid-latitude f -plane with full Coriolisforce. In Section 6, we derive the first order balance model Lagrangian. To ensure3hat the Lagrangian is maximally degenerate, we carefully decompose the expressionfor kinetic energy into the parts that can be removed by a suitable choice of thetransformation vector field and a residual component which cannot be removed. InSection 7, we derive the balance model Euler–Poincar´e equations from the balancemodel Lagrangian. In Section 8, we find that it can be decomposed into a singleevolution equation for a scalar field in the two horizontal variables, and a kinematicrelationship for all other components. This kinematic relationship is analyzed in moredetail in Section 9, where we show that it is elliptic under the assumption of stablestratification and sufficiently small fluctuations in all prognostic fields. Section 10discusses the reconstruction of the full velocity field from the balance relation and theprognostic variables. We finally give a brief discussion and conclusions.
2. Variational principle for the Boussinesq equations
In this section, we recall the derivation of the equations of motion via the Hamiltonprinciple, here in the abstract setting of Euler–Poincar´e theory. We also set up thederivation of reduced equations of motion via degenerate variational asymptotics.We write Diff µ ( D ) to denote the group of H s -class volume-preserving diffeomor-phisms on D that leave its boundary invariant. The associated “Lie algebra” of vectorfield is the space V div = { u ∈ H s ( D , R ) : ∇ · u = 0, k · u = 0 at z = −
1, 0 } . (5)Here, H s is the usual Sobolev space of order s consisting of functions with squareintegrable weak derivatives up to order s . For s > /
2, Diff µ ( D ) is a smooth infinite-dimensional manifold and V div is its tangent space at the identity (Ebin & Marsden,1970; Palais, 1968).Let η = η ( · , t ) ∈ Diff µ ( D ) denote the time-dependent flow generated by a time-dependent vector field u = u ( · , t ) ∈ V div , i.e.,˙ η ( a , t ) = u ( η ( a , t ), t ) (6)or ˙ η = u ◦ η for short; here and in the following we use the symbol “ ◦ ” to denotethe composition of maps. In the Lagrangian description of fluid flow, x ( t ) = η ( a , t ) isthe trajectory of a fluid parcel initially located at a ∈ D . We retain the letter x forEulerian positions and a for Lagrangian labels throughout.To formulate the variational principle, we note that the continuity equation (1c) isequivalent to ρ ◦ η = ρ , (7)where ρ denotes the given initial density field. The Lagrangian for the non-dimensionalEuler–Boussinesq system is given by (Franzke et al., 2019) L ( η , ˙ η ; ρ ) = Z D R ◦ η · ˙ η + ε | ˙ η | − ρ k · η d a = Z D R · u + ε | u | − ρ z d x ≡ ℓ ( u , ρ ) . (8)4he vector R is a vector potential for the Coriolis vector, i.e., ∇ × R = Ω . Here, onthe f -plane, we can choose R = 12 J x = 12 cz − sysx − cx , (9)where J is the skew-symmetric matrix J = − s cs − c . (10)As we see in the second line of (8), the Lagrangian can be written as a functional of u and ρ alone. In this form, it is referred to as the reduced Lagrangian ℓ . Computingvariations of the full Lagrangian L with respect to the flow map η is equivalent tocomputing variations of the reduced Lagrangian ℓ with respect to u and ρ subject toso-called Lin constraints. This is summarized in the following theorem. Theorem 2.1 (Holm, Marsden, & Ratiu, 1998) . Consider a curve η in Diff µ ( D ) withLagrangian velocity ˙ η and Eulerian velocity u ∈ V div . Then the following are equiva-lent. (i) η satisfies the Hamilton variational principle δ Z t t L ( η , ˙ η ; ρ ) dt = 0 (11) with respect to variations of the flow map δ η = v ◦ η , where v is a curve in V div vanishing at the temporal end points. (ii) u and ρ satisfy the reduced Hamilton principle δ Z t t ℓ ( u , ρ ) dt = 0 (12) with respect to variations δ u and δρ that are subject to the Lin constraints δ u =˙ v + u · ∇ v − v · ∇ u and δρ + v · ∇ ρ = 0 , where v is a curve in V div vanishingat the temporal end points. (iii) m and ρ satisfy the Euler–Poincar´e equation Z D ( ∂ t + L u ) m · v + δℓδρ L v ρ d x = 0 (13) for every v ∈ V div , where L u is the Lie derivative in the direction of u and m = δℓδ u (14) is the momentum one-form.
5n the language of vector fields on a region of R , (13) reads Z D (cid:16) ∂ t m + ( ∇ × m ) × u + ∇ ( m · u ) + δℓδρ ∇ ρ (cid:17) · v d x = 0 (15)for every v ∈ V div . This implies that the term in parentheses must be zero up to agradient of a scalar potential φ .For the Euler–Boussinesq Lagrangian (8), m = δℓδ u = R + ε u and δℓδρ = − z (16)Setting φ = − p − z ρ − ε | u | and using the vector identity ( ∇ × u ) × u = u · ∇ u − ∇ | u | , we see that (15) implies the Euler–Boussinesq momentum equation (1).As the Lagrangian L is invariant under time translation, the corresponding Hamil-tonian H = Z D δℓδ u · u d x − ℓ ( u , ρ ) = Z D ε | u | + z ρ d x (17)is a constant of the motion. A second conservation law arises via the invariance ofthe Lagrangian with respect to “particle relabeling”, i.e., composition of the flowmap with an arbitrary time-independent map in Diff µ ( D ): It can be shown (see, e.g.,Franzke et al., 2019, for a direct verification) that the potential vorticity q = ( ∇ × m ) · ∇ ρ = ( Ω + ε ∇ × u ) · ∇ ρ (18)is materially conserved, i.e., satisfies the advection equation ∂ t q + u · ∇ q = 0 . (19)
3. Variational asymptotics
Our construction of variational balance models is based on the following construction.Suppose that the flow of the balance model is related to the flow of the full model viaa near-identity change of variables. To be explicit, let˙ η ε = u ε ◦ η ε (20)represent Lagrangian resp. Eulerian velocity of the Boussinesq system, and let˙ η = u ◦ η (21)denote the corresponding balance model velocities. We suppose that the change ofcoordinates is the flow of a vector field where ε takes the role of an artificial timeparameter, namely η ′ ε = v ε ◦ η ε , η ε (cid:12)(cid:12) ε =0 = η . (22)6ere and in the following, we use the prime symbol to denote a derivative with respectto ε . Cross-differentiation of (20) and (22) gives u ′ ε = ˙ v ε + u ε · ∇ v ε − v ε · ∇ u ε . (23)Likewise, differentiation in ε of the Lagrangian density gives ρ ′ ε + v ε · ∇ ρ ε = 0 . (24)These two relations can be seen as the Lin constraints for the ε -flow, cf. the Linconstraints for the t -flow stated in Theorem 2.1.At this point, the choice of v ε is completely arbitrary and no assumptions have beenmade. However, we will always restrict ourselves to incompressible transformationswhich leave the domain invariant, so that v ε ∈ V div .We now treat ε as a perturbation parameter, expanding (23) and (24) in ε . Likewise,we expand Lagrangian as a formal power series in ε and use the Lin constraints toeliminate all ε -derivatives of u and ρ . Model reduction is then achieved in the followingsteps.(i) Truncate the expansion of the Lagrangian at some fixed order. Here, we will onlylook at the truncation to O ( ε ), the first nontrivial order.(ii) Choose the expansion coefficients of the transformation vector fields such thatthe resulting Lagrangian is maximally degenerate. When computing the modelreduction to first order, the zero-order transformation vector field v = v ε | ε =0 isthe only choice to be made.(iii) Apply the Euler–Poincar´e variational principle to derive the equations of motion.Since the Lagrangian is degenerate, some degrees of freedom will be kinematic,thus imply a phase-space constraint that can be understood as a so-called Diracconstraint (e.g. Salmon, 1988).The first two steps are done in reverse order relative to the original method ofSalmon (1985) who constrained the system first, using the readily available leadingorder balance relation, and transformed to different coordinates second. In fact, in hisapproach, the second step is optional and it turns out that, for f -plane shallow water,the so-called L model, which skips the transformation, is superior to all other modelsin a more general family (Dritschel et al., 2017). In our viewpoint, the transformationis always necessary. The advantage is that the constraint arises as a consequence ofthe formalism and does not need to be guessed or derived a priori . In simple cases,it is possible to choose the transform such that it vanishes to the order considered;in this case, we reproduce Salmon’s L model. In more complicated cases, in partic-ular in our setting here, it is not possible to cancel all terms in the transformationvector field at the order considered. This means that a direct application of Salmon’smethod would result in a model with additional spurious prognostic variables. Wefinally remark that at least in finite dimensions, the procedure is rigorous and theresulting model is correct to the expected order of approximation (Gottwald & Oliver,2014). For partial differential equations, this question is open; it is clear, though, thatadditional conditions are necessary. 7 . Oblique vertical coordinate and decomposition of vector fields The axis of rotation plays a crucial role because it defines the direction of the charac-teristics of the thermal wind relation. It is a distinguished direction for many aspects ofour derivation. In particular, we need to decompose vector fields into the componentsaligned with and perpendicular to this axis. As it is generally not aligned with thegeometric vertical, we need to collect a number of identities and elementary results tofacilitate working in this oblique coordinate setting.We begin by defining a projector Q = ΩΩ T = c cs cs s (25)in the direction of the axis of rotation; the projector onto the orthogonal complementis then given by P = I − Q = s − cs − cs c . (26)We note that Ω spans the kernel of J and is perpendicular to the range of J . Thus, Q is the orthogonal projector onto Ker J and P is the orthogonal projector onto Range J .In the following, we need to integrate along the axis of rotation. To do so, thefollowing parameterization is convenient: We write x = xy + ζ Ω ≡ A ξ ≡ χ ( ξ ) (27)where ξ ≡ ( x , y , ζ ) T , x and y are the horizontal coordinates of the characteristic lineat the surface, ζ is the arclength parameter along the characteristic line, with ζ = 0being at the surface and ζ = − s − being at the bottom, and A = c s . (28)We note that det A = s . Further,for any scalar function f and with the mapping χ defined in (27), ∇ ( f ◦ χ ) = ( A T ∇ f ) ◦ χ , (29a)( ∇ f ) ◦ χ = A − T ∇ ( f ◦ χ ) , (29b)and, in components, ∂ i ( f ◦ χ ) = ( ( ∂ i f ) ◦ χ for i = x , y ( Ω · ∇ f ) ◦ χ for i = ζ . (29c)8hat, for an arbitrary vector field v , ∇ · ( v ◦ χ ) = ( ∇ · A v ) ◦ χ . (30)In particular, setting v = A − P u for some vector field u , and noting that A − P = SP with S = diag(1, s − , s − ) , (31)(30) turns into ( ∇ · P u ) ◦ χ = ∇ · ( SP u ◦ χ ) . (32)For any function φ , we define ¯ φ as its mean along the axis of rotation, i.e.,¯ φ ◦ χ = s Z − s − φ ◦ χ dζ . (33)Further, we write ˆ φ = φ − ¯ φ to denote the deviation from the mean. The definition forvector fields is analogous. Since, for arbitrary functions φ and ψ , ¯ ψ ◦ χ is independentof ζ and det A = s , we see that Z D ˆ φ ¯ ψ d x = s Z T Z − s − ˆ φ ◦ χ dζ ¯ ψ ◦ χ dx = 0 . (34)In other words, mean and fluctuating components in the sense of (33) are L -orthogonal, which we will use without explicit further mention. The following lemmastates that even though averaging is done along an oblique line, some of the expectedrelations of an averaging operator hold true. Lemma 4.1.
Let φ and ψ be arbitrary functions and v an arbitrary vector field on D . Then (i) ψ ¯ φ = ¯ ψ ¯ φ , (ii) v · ∇ ¯ φ = ¯ v · ∇ ¯ φ . Proof.
For (i), we note that ¯ φ ◦ χ is independent of ζ , so that ψ ¯ φ ◦ χ = s Z − s − ψ ◦ χ ¯ φ ◦ χ dζ = s Z − s − ψ ◦ χ dζ ¯ φ ◦ χ = ¯ ψ ◦ χ ¯ φ ◦ χ .For (ii), we note that Ω · ∇ ¯ φ = 0, so that the z -derivative can be replace by anequivalent y -derivative. Since horizontal derivatives commute with taking the average,part (i) applies and yields the claim.Commutation of vertical derivatives of arbitrary functions with averaging is moresubtle, as the next lemma shows. Here and in the following, we shall write χ (0) and χ ( − s − ) to indicate that the expression is evaluated at the top ( ζ = 0) or at thebottom boundary ( ζ = s − ), as a function of the remaining horizontal variables.9 emma 4.2. Let φ be a function with φ ◦ χ (0) = φ ◦ χ ( − s − ) and v an arbitraryvector field on D . Then (i) ∂ z ¯ φ = ∂ z φ , (ii) ¯ v · ∇ φ = ¯ v · ∇ ¯ φ . Proof.
On the one hand, Ω · ∇ ¯ φ = 0, so that s ∂ z ¯ φ = − c ∂ y ¯ φ . On the other hand, ∂ z φ ◦ χ = s Z − s − ( ∂ z φ ) ◦ χ dζ = Z − s − ( ∂ ζ − c ∂ y )( φ ◦ χ ) dζ = − cs ∂ y ¯ φ ◦ χ + φ ◦ χ (0) − φ ◦ χ ( − s − ) . (35)Under the condition stated, the boundary terms cancel. This implies (i). For (ii), weuse Lemma 4.1 to move ¯ v out of the average. Horizontal derivatives can be moved outof the average without restrictions, for the z -derivative, part (i) applies.We proceed to prove a number of lemmas which describe the splitting of divergence-free vector fields with zero-flux boundary conditions into, on the one hand, mean andmean-free components and, on the other hand, components along and perpendicularto the axis of rotation. Lemma 4.3.
Let v ∈ V div . Then ∇ · S h P ¯ v = 0 . Proof.
For a divergence-free vector field, ∇ · P v = − ∇ · Q v = − Ω · ∇ ( Ω · v ), so that(32) turns into ∇ · ( S h P v ◦ χ ) = − ∂ ζ ( Ω · v ◦ χ ) − ∂ ζ ( S P v ◦ χ ) = − s − ∂ ζ ( v ◦ χ ) (36)where the last equality can be verified by direct computation in coordinates. Integrat-ing in ζ and noting that the right hand side is zero due to the boundary conditions,we obtain the statement of the lemma.The following is a somewhat weaker converse of Lemma 4.3. Lemma 4.4.
Let v be a vector field with ∇ · S h P ¯ v = 0 . Then ∇ · ¯ v = 0 . Proof.
The assumption implies 0 = ∇ · ¯ v − cs ∂ y ¯ v . (37)Since ¯ v ◦ χ is independent of ζ , this implies ∇ · ( A − ¯ v ◦ χ ) = 0. By (30), this impliesthat ∇ · ¯ v = 0. Corollary 4.5. If v ∈ V div , then ∇ · ˆ v = ∇ · ¯ v = 0 . Proof.
Lemma 4.3 followed by Lemma 4.4 yields ∇ · ¯ v = 0. Then ˆ u = u − ¯ u is alsodivergence-free. 10 emma 4.6. Let u ∈ V div . Then Q ˆ u is uniquely determined by P ˆ u and given by theformula Q ˆ u = Ω ˆ g with ˆ g ◦ χ = − s Z − s − ζ ∇ · ( SP ˆ u ◦ χ ) dζ − Z ζ − s − ∇ · ( SP ˆ u ◦ χ ′ ) dζ ′ (38) where we write χ ′ to abbreviate χ ( x , y , ζ ′ ) . Remark 1.
Setting U ◦ χ = Z ζ − s − ˆ u ◦ χ ′ dζ ′ , (39)we can write (38) in the more compact formˆ g = − ∇ · P ˆ U . (40) Proof.
With u = P u + Q u , the divergence condition reads0 = ∇ · u = ∇ · P u + ∇ · Q u . (41)Recalling that Q = ΩΩ T and setting g = Ω · u , we can write Q u = Ω g . Then, Ω · ∇ ˆ g = Ω · ∇ g = − ∇ · P u (42)so that, due to (29c), ∂ ζ (ˆ g ◦ χ ) = ( Ω · ∇ ˆ g ) ◦ χ = − ( ∇ · P u ) ◦ χ = − ∇ · ( SP u ◦ χ ) = − ∇ · ( SP ˆ u ◦ χ ) . (43)We note that the third equality in (43) is due to (32) and the last equality is due to u = ˆ u + ¯ u , where ¯ u ◦ χ is independent of ζ , so that the entire contribution from ¯ u vanishes due to Lemma 4.3. Equation (43) determines g uniquely up to a constant ofintegration on each of the characteristic lines; we write g = ˆ g + ¯ g and choose ¯ g as thisconstant of integration. The condition that ˆ g is mean-free along each characteristicline implies 0 = s Z − s − ˆ g ◦ χ dζ = ˆ g ◦ χ ( − s − ) − s Z − s − ζ ∂ ζ (ˆ g ◦ χ ) dζ . (44)Then, substituting (43) and (44) intoˆ g ◦ χ ( ζ ) = ˆ g ◦ χ ( − s − ) + Z ζ − s − ∂ ζ (ˆ g ◦ χ ′ ) dζ ′ , (45)we obtain (38), which determines Q ˆ u uniquely. Lemma 4.7.
Let u ∈ V div ; we write ˆ u to denote its mean-free component as before.Further, let ˆ w be any mean free vector field. Then Z D ˆ w · Q ˆ u d x = Z D ˆ u · P A [ ˆ w ] d x (46)11 ith A [ ˆ w ] defined by A [ ˆ w ] ◦ χ = − S ∇ Z ζ − s − Ω · ˆ w ◦ χ ′ dζ ′ . (47) Proof.
By Lemma 4.6, Z D ˆ w · Q ˆ u d x = Z D ˆ w · Ω ˆ g d x = s Z T Z − s − Ω · ˆ w ◦ χ ˆ g ◦ χ d ξ , (48)where ˆ g ◦ χ is given by (38). As it is integrated against a mean-free vector field, thefirst term on the right of (38) does not contribute to the integral (48), so that Z D ˆ w · Q ˆ u d x = − s Z T Z − s − Ω · ˆ w ◦ χ Z ζ − s − ∇ · ( SP ˆ u ◦ χ ′ ) dζ ′ d ξ = s Z T Z − s − ∇ · ( SP ˆ u ◦ χ ) Z ζ − s − Ω · ˆ w ◦ χ ′ dζ ′ d ξ = − s Z T Z − s − ˆ u ◦ χ · PS ∇ Z ζ − s − Ω · ˆ w ◦ χ ′ dζ ′ d ξ = Z D ˆ u · P A [ ˆ w ] d x (49)where A [ ˆ w ] is given by (47). We remark that the second inequality is based on integra-tion by parts in ζ , the third equality is due to the divergence theorem. In both cases,the boundary terms vanish due to the mean-free condition on ˆ w . We have further usedthe symmetry of the matrices S and P . Lemma 4.8.
Let u = ˆ u + ¯ u ∈ V div and let ¯ φ be the vertical mean of an arbitraryscalar field. Then Z D ¯ φ ¯ u d x = − Z D ˆ u · P ∇ ¯ φ z d x . (50) Proof.
By Lemma 4.6, Q ˆ u = Ω ˆ g withˆ g ◦ χ ( − s − ) = − s Z − s − ζ ∇ · ( SP ˆ u ◦ χ ) dζ . (51)At the bottom boundary, k · ¯ u = − k · ˆ u = − k · P ˆ u − k · Q ˆ u , so that Z D ¯ φ ¯ u d x = s Z T Z − s − ¯ φ ◦ χ (cid:2) − k · P ˆ u ◦ χ ( − s − ) − s ˆ g ◦ χ ( − s − ) (cid:3) d ξ = Z T ¯ φ ◦ χ (cid:20) − k · P ˆ u ◦ χ ( − s − ) dx + s Z − s − ζ ∇ · ( SP ˆ u ◦ χ ) dζ (cid:21) dx = − Z T ¯ φ ◦ χ k · P ˆ u ◦ χ ( − s − ) dx + s Z T Z − s − ¯ φ ◦ χ (cid:0) ∇ · ( ζ SP ˆ u ◦ χ ) − ∇ ζ · SP ˆ u ◦ χ (cid:1) d ξ . (52)12ince ∇ ζ = k , the second term in the last integral vanishes as the vertical integrationis over a mean free quantity. The first term in the last integral, we integrate by parts.The boundary term from the upper boundary is zero. The boundary term from thelower boundary exactly cancels the integral on the second last line, so that Z D ¯ φ ¯ u d x = − s Z T Z − s − ˆ u ◦ χ · PS ∇ ( ¯ φ ◦ χ ) ζ d ξ = − s Z T Z − s − ˆ u ◦ χ · PS ( A T ∇ ¯ φ z ) ◦ χ d ξ , (53)where the last equality is due to (29a) and z = sζ . Noting that PSA T = P and changingback to Cartesian coordinates, we obtain (50).The next lemma provides a converse statement to Lemma 4.6. Lemma 4.9.
Suppose P ˆ u is a given vector field which is mean-free and contained inthe range of P . Define ˆ g as in Lemma 4.6, i.e., ˆ g ◦ χ = − s Z − s − ζ ∇ · ( SP ˆ u ◦ χ ) dζ − Z ζ − s − ∇ · ( SP ˆ u ◦ χ ′ ) dζ ′ . (54) Then ˆ u = P ˆ u + Q ˆ u with Q ˆ u = Ω ˆ g is mean-free and divergence-free. Proof.
The fact that ˆ g , hence ˆ u , is mean-free is a direct consequence of the choice ofconstant of integration in the proof of Lemma 4.6.To prove that ˆ u is divergence-free, we take the ζ -derivative of (54), ∂ ζ (ˆ g ◦ χ ) = − ∇ · ( SP ˆ u ◦ χ ) . (55)This implies that 0 = ∂ ζ (ˆ g ◦ χ ) + ∇ · ( SP ˆ u ◦ χ )= ( Ω · ∇ ˆ g ) ◦ χ + ( ∇ · P ˆ u ) ◦ χ = ( ∇ · Q ˆ u ) ◦ χ + ( ∇ · P ˆ u ) ◦ χ = ( ∇ · ˆ u ) ◦ χ . (56)Since χ is invertible, we find that ∇ · ˆ u = 0. Lemma 4.10.
Under the conditions of Lemma 4.9, there exists ¯ u , also divergencefree, such that u = ˆ u + ¯ u satisfies the zero-flux boundary condition k · u = 0 at z = 0, − . Proof.
By Lemma 4.9, ˆ u is divergence free. We now choose ¯ u such that the vectorfield u is tangent to the top and bottom boundaries. At z = −
1, 0, we require k · ¯ u = − k · ˆ u . (57)Observe thatˆ u ◦ χ (0) = k · P ˆ u ◦ χ (0) + s ˆ g ◦ χ (0)13 k · P ˆ u ◦ χ (0) + s ˆ g ◦ χ ( − s − ) − s Z − s − ∇ · ( SP ˆ u ◦ χ ) dζ = k · P ˆ u ◦ χ (0) + s ˆ g ◦ χ ( − s − ) − k · P ˆ u ◦ χ (0) + k · P ˆ u ◦ χ ( − s − )= k · P ˆ u ◦ χ ( − s − ) + s ˆ g ◦ χ ( − s − ) . (58)This implies ˆ u takes the same value at the bottom and top boundaries along any linein the direction of the axis of rotation. Consequently, we can use (57) to define k · ¯ u at the bottom, i.e., we set¯ u ◦ χ = − k · P ˆ u ◦ χ ( − s − ) − s ˆ g ◦ χ ( − s − ) . (59)Then the boundary condition (57) is satisfied at z = 0 as well.We now choose ¯ u such that u is divergence-free. Indeed, due to Lemma 4.4, it sufficesto ensure that ∇ · S h P ¯ u = 0, cf. (37) for an explicit expression. Setting ¯ u = ∇ φ , wesee that this implies ∆ φ = cs ∂ y ¯ u , (60)which can be solved as a Poisson equation on T . Lemma 4.11.
Let u , v ∈ V div be such that their domain-mean is zero. Then Z D u · J v d x = Z D ˆ u · J ˆ v d x . (61) Proof.
The mean of u over the domain D is zero if and only if the horizontal mean of ¯ u is zero, and likewise for v . Moreover, by Lemma 4.3, ∇· ( S h P ¯ u ) = 0 and ∇· ( S h P ¯ v ) = 0.Thus, there exist scalar fields ψ and θ such that S h P ¯ u = ∇ ⊥ ψ and S h P ¯ v = ∇ ⊥ θ . Wewrite S − h = s − cs (62)to denote the pseudo-inverse of S h on Range P . It satisfies S h S − h = I , (63a)where I denotes the 2 × S − h S h P = P . (63b)Then, P ¯ u = S − h ∇ ⊥ ψ and P ¯ v = S − h ∇ ⊥ θ . Further, observing that S − T h J S − h = s J , (64)where J denotes the canonical 2 × PJ = J = JP ,14e compute Z D ¯ u · J ¯ v d x = Z D P ¯ u · JP ¯ v d x = Z D S − h ∇ ⊥ ψ · JS − h ∇ ⊥ θ d x = Z D ∇ ⊥ ψ · S − T h JS − h ∇ ⊥ θ d x = − s Z D ∇ ⊥ ψ · ∇ θ d x (65)By orthogonality of gradients and curls, the last integral is zero, which implies (61).The balance relation we are going to derive in Section 7 below is most convenientlyformulated in terms of U , the anti-derivative of ˆ u along the axis of rotation. Lemma 4.12.
Let u = ˆ u + ¯ u ∈ V div and define U via (39) . Then u = − s ∇ · S h P U (66a) and ¯ u = − s ∇ · U . (66b) Proof.
By incompressibility, Ω · ∇ u = c ∂ y u − s ∇ · u = − s ∇ · S h P u = − s ∇ · S h P ˆ u . (67)where the last identity is due to Lemma 4.3. Integration along the axis of rotation anduse of the zero-flux boundary condition at z = − u ◦ χ = − s Z ζ − s − ( ∇ · S h P ˆ u ) ◦ χ ′ dζ ′ . (68)By the definition of U via (39), this implies (66a). Further, noting that ˆ u = Ω · ∇ U ,we obtain ¯ u = − s ∇ · S h P U − Ω · ∇ U . (69)By direct computation, the right hand side of (69) simplifies to − s ∇ · U .
5. Thermal wind relation
The formal leading order balance in the Euler–Boussinesq momentum equation (1a)is Ω × u = − ∇ p − ρ k . (70)15o remove the pressure, we take the curl of this equation, Ω · ∇ u g = (cid:18) −∇ ⊥ ρ (cid:19) . (71)This expression is known as the thermal wind relation . The characteristics of this firstorder equation are the lines in the direction of the axis of rotation. With the notationset up in Section 4, it is easy to integrate (71) along its characteristic lines. Splitting u g = ˆ u g + ¯ u g , we first note that (71) determines only the mean-free component ˆ u g .Indeed, ∂ ζ (ˆ u g ◦ χ ) = ( Ω · ∇ ˆ u g ) ◦ χ = (cid:18) −∇ ⊥ ρ ◦ χ (cid:19) (72)The mean-free component ˆ u g must further satisfy0 = s Z − s − ˆ u g ◦ χ dζ = ˆ u g ◦ χ ( − s − ) − s Z − s − ζ ∂ ζ (ˆ u g ◦ χ ) dζ . (73)Then, ˆ u g ◦ χ = ˆ u g ◦ χ ( − s − ) + Z ζ − s − ∂ ζ (ˆ u g ◦ χ ′ ) dζ ′ = s Z − s − ζ ∂ ζ (ˆ u g ◦ χ ) dζ + Z ζ − s − ∂ ζ (ˆ u g ◦ χ ′ ) dζ ′ . (74)Inserting the thermal wind relation (72), we find that the vertical component of thethermal wind is independent of ζ , so k · ˆ u g = 0, and the horizontal components of thethermal wind are given byˆ u g ◦ χ = − s Z − s − ζ ∇ ⊥ ρ ◦ χ ′ dζ ′ − Z ζ − s − ∇ ⊥ ρ ◦ χ ′ dζ ′ (75)or ˆ u g = ∇ ⊥ θ (76a)where θ ◦ χ = − s Z − s − ζ ρ ◦ χ dζ − Z ζ − s − ρ ◦ χ ′ dζ ′ . (76b)Equation (76) shows that ˆ u g is horizontally divergence-free. Since k · ˆ u g = 0, we alsohave ∇ · u g = 0. 16 . Derivation of the first-order balance model We now follow the procedure outlined in Section 3 to first order in ε . Using the Linconstraints for the ε -flow, (23) and (24), we can write u ε = u + ε ( ˙ v + u · ∇ v − v · ∇ u ) + O ( ε ) (77)and ρ ε = ρ − ε ∇ · ( ρ v ) + O ( ε ) . (78)The transformation vector field v will be specified in the following. By construction,we assume that all flows are volume-preserving and leave the domain invariant, so that u , v ∈ V div . For technical reasons, we also assume that the domain-mean of u is zero.This is not a restriction since the assumption is only removing a steady solid-bodytranslation from the system, which is a constant of the motion of the full Euler–Boussinesq system so that it remains zero if it vanishes initially. The transformationvector field v shall also be chosen as not to generate a solid body translation in balancemodel coordinates.Inserting these relations into the Euler–Boussinesq Lagrangian (8), expanded inpowers of ε and truncated to O ( ε ), we obtain after a short computation L = Z D R ◦ η ε · ˙ η ε + ε | ˙ η ε | d a − Z D ρ ε z d x = Z D u · J x − ρ z d x + ε Z D u · J v + | u | − ρ v · k d x = L + ε L . (79)We note that we have used the divergence theorem to rewrite the potential energy con-tribution to L , where the boundary integral vanishes due to the boundary condition k · v = 0.We proceed to rewrite L . The transformation vector field v must be chosen in aform that makes the Lagrangian at order ε affine. To determine the transformationvector field, we decompose the kinetic energy part of L as Z D | u | d x = Z D | ˆ u | + 2 ˆ u · ¯ u + | ¯ u | d x . (80)We note that the cross term in (80) vanishes due to (34). The square terms are decom-posed into terms that can be written as an L -pairing with P ˆ u , and a final remainderterm which cannot. Starting with the contribution from | ˆ u | , we write Z D | ˆ u | d x = Z D ˆ u · P ˆ u + ˆ u · Q ˆ u d x = Z D ˆ u · P (ˆ u + A [ˆ u ]) d x , (81)using Lemma 4.7 in the final equality. The contribution from | ¯ u | is split differently.Setting S = diag(1, s − , 0), noting that s ( I − S P )¯ u = Ω ¯ u , and noting that ( S P ) T − P is skew so that S P ¯ u · ¯ u − ¯ u · S P ¯ u = 0, we write Z D | ¯ u | d x = Z D ( I + S P )¯ u · ( I − S P )¯ u + S P ¯ u · S P ¯ u d x = − s − Z D ˆ u · P ∇ ( Ω · ( I + S P )¯ u ) z d x + Z D | S P ¯ u | d x . (82)The second equality is based on Lemma 4.8 with ¯ φ = s − Ω · ( I + S P )¯ u . Note thatLemma 4.8 could be used in different ways, so long as we retain a pairing of somefunction ¯ φ with ¯ u . Our chosen splitting is distinguished, due to Lemma 4.3, by thefact that the remainder integral in (82) is proportional to the kinetic energy of a divergence-free two-dimensional vector field. In the following, this remainder term willbe the only term that yields an evolution equation in the variational principle. As S P ¯ u is divergence-free, this evolution equation can always be written in terms of ascalar stream function, i.e., is determined by the evolution of a single scalar field. Anyother splitting would yield a two-component evolution equation except for one specialtilt of the axis of rotation, which is not desirable.Collecting terms, we find Z D | u | d x = Z D ˆ u · P ˆ V d x + Z D | S P ¯ u | d x (83)with V = ˆ u + A [ˆ u ] − s − ∇ (cid:0) Ω · ( I + S P )¯ u (cid:1) z . (84)We note that V is not necessarily mean-free but, due to the pairing with ˆ u , only itsmean-free component ˆ V contributes to the Lagrangian.Since J has a one-dimensional kernel, it is impossible to remove all quadratic termsfrom the L Lagrangian, but choosing the transformation vector field v as a solutionto the equation J v = − P ˆ V , (85)up to terms that only depend on ρ , we can make L affine in ˆ u . Indeed, noting that J is invertible on Range P with pseudo-inverse J T , and further that PJ T = J T = J T P ,we seek v in the form v = P ˆ v + Q ˆ v + ¯ v (86)with P ˆ v = − J T ˆ V + λ J T ˆ V g . (87)In this expression, λ is a free parameter and V g = ˆ u g + A [ˆ u g ] . (88)18e note that V g takes the same form as V , with u replaced by u g . Since u g is notconstrained by the thermal wind relation, we do not include a term that matches thethird term of (84) into the ansatz for V g .The terms Q ˆ v and ¯ v in (86) are chosen in the following such that v becomesdivergence free and tangent to the top and bottom boundaries. Using the constructionfrom Lemma 4.6, we write Q ˆ v ≡ Ω ˆ g , whereˆ g ◦ χ = ˆ g ◦ χ ( − s − ) − Z ζ − s − ∇ · ( SP ˆ v ◦ χ ′ ) dζ ′ . (89)By Lemma 4.9, ˆ v is divergence-free. By Lemma 4.10, ¯ v can be chosen such that thezero flux boundary condition¯ v ◦ χ = − k · P ˆ v ◦ χ ( − s − ) − s ˆ g ◦ χ ( − s − )= − k · P ˆ v ◦ χ (0) − s ˆ g ◦ χ (0) . (90)are satisfied and, moreover, ¯ v is divergence free. Lemma 4.11 shows that the choiceof the horizontal components ¯ v will not enter the computation of the equations ofmotion. Only the contribution from ¯ v will appear, but it is directly a function of ˆ v by (90), so that the final expression will not contain any references to ¯ v .Combining the contributions from rotation and kinetic energy, we find Z D u · J v + | u | d x = Z D ˆ u · JJ T (cid:0) − ˆ V + λ ˆ V g (cid:1) + ˆ u · P ˆ V + | S P ¯ u | d x = Z D λ ˆ u · P ˆ V g + | S P ¯ u | d x = Z D λ ˆ u · ˆ u g + | S P ¯ u | d x (91)where, in the last equality, we have used Lemma 4.7 and the fact that the geostrophicvelocity is exclusively horizontal.Inserting the boundary condition for the transformation vector field (90) and therepresentation of Q ¯ v via (89), we rewrite the potential energy contribution to the L -Lagrangian as follows: − Z D ρ v · k d x = − s Z T Z − s − ρ ◦ χ (cid:0) k · ¯ v ◦ χ + k · P ˆ v ◦ χ + k · Q ˆ v ◦ χ (cid:1) d ξ = s Z T Z − s − ρ ◦ χ (cid:20) k · P ˆ v ◦ χ ( − s − ) + s ˆ g ◦ χ ( − s − ) − k · P ˆ v ◦ χ − s ˆ g ◦ χ ( − s − ) + s Z ζ − s − ∇ · ( SP ˆ v ◦ χ ′ ) dζ ′ (cid:21) d ξ = s Z T Z − s − ρ ◦ χ (cid:20) k · P ˆ v ◦ χ ( − s − ) − k · P ˆ v ◦ χ + Z ζ − s − ∂ ζ ( k · P ˆ v ◦ χ ′ ) dζ ′ + s Z ζ − s − ∇ · ( S h P ˆ v ◦ χ ′ ) dζ ′ (cid:21) d ξ = − s Z T Z − s − ∇ ρ ◦ χ · Z ζ − s − S h P ˆ v ◦ χ ′ dζ ′ d ξ (92)19urther, inserting (87), using the identity S h J T = − s − J P h , (93)and noting that horizontal gradients and composition with χ commute, we obtain − Z D ρ v · k d x = − s Z T Z − s − ∇ ρ ◦ χ · Z ζ − s − S h J T (cid:0) − ˆ V + λ ˆ V g (cid:1) ◦ χ ′ dζ ′ d ξ = s Z T Z − s − ∇ ρ ◦ χ · Z ζ − s − J P h (cid:0) − ˆ V + λ ˆ V g (cid:1) ◦ χ ′ dζ ′ d ξ = s Z T Z − s − ∇ ⊥ ρ ◦ χ · Z ζ − s − P h (cid:0) ˆ V − λ ˆ V g (cid:1) ◦ χ ′ dζ ′ d ξ . (94)Inserting the thermal wind relation (72), integrating by parts with respect to ζ , andchanging variables, we continue the computation: − Z D ρ v · k d x = − s Z T Z − s − ∂ ζ (ˆ u g ◦ χ ) · Z ζ − s − P h (cid:0) ˆ V − λ ˆ V g (cid:1) ◦ χ ′ dζ ′ d ξ = Z D ˆ u g · P h (cid:0) ˆ V − λ ˆ V g (cid:1) d x = Z D ˆ u g · (cid:0) ˆ u − λ ˆ u g (cid:1) d x (95)where, in the last step, we have made use of Lemma 4.7, Lemma 4.8, and the fact that k · u g = 0.Altogether, the L -Lagrangian reads L = Z D λ ˆ u · ˆ u g + | S P ¯ u | +ˆ u g · ( ˆ u − λ ˆ u g ) d x = Z D ν ˆ u · ˆ u g + ¯ u · M ¯ u − λ | ˆ u g | d x , (96)where ν = λ + and M = PS P = − c/s − c/s c /s . (97)In the following, we choose λ such that ν >
7. Derivation of the balance model equations of motion
Taking the variation of (96), we obtain δL = Z D δ u · ( ν ˆ u g + M ¯ u ) + ν u · δ ˆ u g − λ δ ˆ u g · ˆ u g d x Z D δ u · p + δ ˆ u g · ˆ b d x , (98)where p = M ¯ u + ν ˆ u g and ˆ b = ν ˆ u − λ ˆ u g . (99)To rewrite the second term on the right of (98), we insert the expression for thegeostrophic velocity (75), change variables, and recall that horizontal derivatives andcomposition with χ commute, so that Z D δ ˆ u g · ˆ b d x = − s Z D Z ζ − s − ∇ ⊥ δρ ◦ χ ′ dζ ′ · ˆ b ◦ χ d ξ = s Z D ∇ ⊥ δρ ◦ χ · Z ζ − s − ˆ b ◦ χ ′ dζ ′ d ξ = − s Z D δρ ◦ χ ∇ ⊥ · Z ζ − s − ˆ b ◦ χ ′ dζ ′ d ξ = − Z D δρ ∇ ⊥ · B d x , (100)where B is the anti-derivative of ˆ b along the axis of rotation, i.e., B = ν U − λ U g (101)where U are the horizontal components of U defined in (39) and U g is defined likewise,with ˆ u g in place of ˆ u .Inserting (100) into (98), we have δL = Z D δ u · p − δρ ∇ ⊥ · B d x . (102)We recall the general Euler–Poincar´e equation (15), which can be written ∂ t m + ( ∇ × m ) × u + δℓδρ ∇ ρ = ∇ ˜ φ (103)with ˜ φ an arbitrary scalar field. Here, m = δℓδ u = R + ε p and δℓδρ = − z − ε ∇ ⊥ · B , (104)so that, with ˜ φ = − φ − zρ Ω × u + ρ k + ε ( ∂ t p + ( ∇ × p ) × u − ∇ ρ ∇ ⊥ · B ) = − ∇ φ . (105)By direct computation, using the fact that c∂ z = − s∂ y when applied to averagedquantities, we find that ∇ × M ¯ u = Ω ω with ω = s − ∇ ⊥ · M h ¯ u . Hence, ξ ≡ ∇ × p = Ω ω + ν γ (106)21ith γ = ∇ × ˆ u g = − ∇ × curl θ = (cid:18) −∇ ∂ z θ ∆ θ (cid:19) , (107)where we recall that ˆ u g = ∇ ⊥ θ , see (76), write curl f = ∇ × (0, 0, f ) to identify thecurl of a scalar field with the curl of a vector field oriented in the vertical, and use ∆ todenote the horizontal Laplacian. With this notation in place, taking the curl of (105)and noting that ∇ × ( ξ × u ) = u · ∇ ξ − ξ · ∇ u as both u and ξ are divergence-free,we find − Ω · ∇ u + ∇ × ( ρ k ) + ε (cid:0) ∂ t ξ + u · ∇ ξ − ξ · ∇ u + ∇ ρ × ∇ ∇ ⊥ · B (cid:1) = 0 . (108)The corresponding balance model Hamiltonian, via (17), is given by H = Z D ρ z + ε (cid:0) ¯ u · M ¯ u + λ | ˆ u g | (cid:1) d x . (109)The potential vorticity for the balance model reads q = ( ∇ × m ) · ∇ ρ = ( Ω + ε Ω ω + εν γ ) · ∇ ρ = ( s + εs ω + εν ∆ θ ) ∂ z ρ + (Ω + ε Ω ω − εν ∂ z ∇ θ ) · ∇ ρ = − ( s + εs ω + εν ∆ θ ) ∂ z ( Ω · ∇ θ ) − (Ω + ε Ω ω − εν ∂ z ∇ θ ) · ∇ ( Ω · ∇ θ ) (110)where, in the last equality, we have used that ρ = − Ω · ∇ θ . Alternatively, we can write q = (cid:12)(cid:12)(cid:12)(cid:12) − ( s + εs ω + εν ∆ θ ) ∇ ( Ω · ∇ θ )(Ω + ε Ω ω − εν ∂ z ∇ θ ) ∂ z ( Ω · ∇ θ ) (cid:12)(cid:12)(cid:12)(cid:12) . (111)We remark that the relation between q and θ can be seen as a second order nonlineardifferential operator of the form q = F ( ω , D θ ) , (112)which is nonlinearly elliptic, Gilbarg and Trudinger (2001), so long as[ F ij ] = ∂F∂ ( D θ ) (113)is positive (or negative) definite. Direct computation shows that[ F ij ] = εν∂ z ρ − εc ξ − ε ( ξ + ν∂ x ρ ) − εc ξ εν∂ z ρ − c ( c + εξ ) − (cid:0) cs + ε ( cξ + sξ + ν∂ y ρ ) (cid:1) − ε ( ξ + ν∂ x ρ ) − (cid:0) cs + ε ( cξ + sξ + ν∂ y ρ ) (cid:1) − s − εs ξ ,(114)so that − F is positive definite provided its principal minors are positive, i.e., − εν∂ z ρ > εν ∂ z ρ + c (cid:2) c + ε (cid:0) ξ + cξ ν∂ z ρ (cid:1)(cid:3) > F < ∂ z ρ <
0, thedeviations from a steady mean state are sufficiently small, and ε is sufficiently small.We note that these conditions are similar to the conditions under which the linearbalance relation is elliptic, which are discussed in more detail in Section 9 below.
8. Separation of balance relation into dynamic and kinematic components
In the following, we shall show that the balance relation is mostly a kinematic relation-ship between density ρ and balanced velocity field. However, there is one horizontalscalar field, the curl of M ¯ u , which evolves dynamically. To see this, we first look atthe vertical component equation, which reads ∂ t ξ + u · ∇ ξ = ξ · ∇ u + ∇ ρ · ∇ ⊥ ∇ ⊥ · B + ε − Ω · ∇ ˆ u . (116)We remark that the seemingly unbalanced O ( ε − )-term in this equation is actuallyonly an O (1)-contribution because u is an O ( ε ) perturbation of u g whose verticalcomponent is zero.Taking the average of (116) along the axis of rotation and using Lemma 4.1 andLemma 4.2 to commute averaging with directional derivatives where possible, we findthat the evolution equation for ξ reduces to a prognostic equation for ω , ∂ t ω + ¯ u · ∇ ω = ∇ ρ · ∇ ⊥ ∇ ⊥ · B − ν ( ∂ z ∇ θ · ∇ u − ∆ θ ∂ z u + u · ∇ ∆ θ ) . (117)All other contributions to (108) are entirely kinematic. To see this, we multiply (117)by Ω and subtract this expression from (108) to obtain0 = − Ω · ∇ u + ∇ × ( ρ k ) + ε (cid:2) ν ∂ t γ + Ω ˆ u · ∇ ω + ν u · ∇ γ − ξ · ∇ u + ∇ ρ × ∇ ∇ ⊥ · B + Ω ∇ ρ · ∇ ⊥ ∇ ⊥ · B + ν Ω ( γ · ∇ u − u · ∇ ∆ θ ) (cid:3) . (118)To remove the remaining time derivative from (118), we proceed as follows. By thedefinition of θ , equation (76b), Ω · ∇ θ = − ρ . Noting that ˆ u = Ω · ∇ U , and using thecontinuity equation in the form (1c), we obtain Ω · ∇ ˙ θ = u · ∇ ρ = Ω · ∇ U · ∇ ρ + ¯ u · ∇ ρ = Ω · ∇ ( U · ∇ ρ ) − U · ∇ ( Ω · ∇ ρ ) + ¯ u · ∇ ρ . (119)We recall that ˙ θ is mean-free by construction. Thus, integrating (119) in the directionof Ω , we find that ˙ θ = U · ∇ ρ − U · ∇ ρ − ˆΓ + ˆΓ (120)23here Γ ◦ χ = Z ζ − s − ( U · ∇ ( Ω · ∇ ρ )) ◦ χ ′ dζ ′ (121a)and Γ ◦ χ = Z ζ − s − (¯ u · ∇ ρ ) ◦ χ ′ dζ ′ . (121b)We can thus rewrite (118) exclusively in terms of U by inserting ˆ u = Ω · ∇ U and(66a), eliminating ˙ γ via (107) and (120), and replacing B by its definition in (101).We obtain0 = − ( Ω · ∇ ) U + ∇ × ( ρ k ) + ε (cid:2) − ν ∇ × curl( U · ∇ ρ ) + ν ∇ × curl ( U · ∇ ρ )+ ν ∇ × curl(ˆΓ − ˆΓ ) + Ω ( Ω · ∇ U ) · ∇ ω − ν ( Ω · ∇ U ) · ∇∇ × curl θ − ν ¯ u · ∇∇ × curl θ − ξ · ∇ ( Ω · ∇ U ) − ξ · ∇ ¯ u + ν ∇ ρ × ∇ ∇ ⊥ · U − λ ∇ ρ × ∇ ∇ ⊥ · U g + ν Ω ∇ ρ · ∇ ⊥ ∇ ⊥ · U − λ Ω ∇ ρ · ∇ ⊥ ∇ ⊥ · U g + νs Ω ( ∇ × curl θ ) · ∇ ∇ · S h P U − ν Ω ( Ω · ∇ U ) · ∇ ∆ θ (cid:3) . (122)We note that there is no contribution from ¯ u to the last term in (122) as θ is mean-free.To expose that the highest order terms appearing in this expression define a secondorder differential operator, we need the following vector identity, which can be seen asa variation of the well-known curl-curl identity. Lemma 8.1.
For a general vector field U , ( ∇ × curl U ) · ∇ ρ − ∇ ρ × ∇ ∇ ⊥ · U = (cid:18) ∂ z ρ ∇ ∇ · U −∇ ρ · ∇ ∇ · U (cid:19) − ∂ z ρ ∆ U + ∇ ρ · ∇ ∂ z U . (123) Proof.
Starting from the left hand side of (123), using the Einstein summation con-vention and the product formula for the Levi–Civita symbols, we compute ε ikl ε lm δ jn ∂ km U n ∂ j ρ − ε ijk ε mn ∂ km U n ∂ j ρ = (cid:0) δ j ( δ im δ kn − δ in δ km ) + δ jm ( δ in δ k − δ i δ kn ) (cid:1) ∂ km U n ∂ j ρ = ∂ ik U k ∂ ρ − ∂ kk U i ∂ ρ + ∂ j U i ∂ j ρ − δ i ∂ jk U k ∂ j ρ . (124)Canceling k ∂ z ∇ · U ∂ z ρ between the first and the fourth term, and ∂ zz U ∂ z ρ betweenthe second and the third term, we obtain the expression on the right of (123).Since u ∈ V div , we can use Lemma 4.12 to remove ∇ · U from the right hand sideof (123), so that the balance relation (122) takes the form L U ≡ A U + S U + K U = g , (125a)where A is the linear second order differential operator acting on scalar functions via Aφ = − ( Ω · ∇ ) φ + εν ∂ z ρ ∆ φ , (125b)24nd acting diagonally on vectors, S contains second order derivatives that are not in A , i.e, S U = ε (cid:2) − ν ∇ ρ · ∇ ∂ z U − ξ · ∇ ( Ω · ∇ U ) + ν ∇ × curl ( U · ∇ ρ )+ ν Ω ∇ ρ · ∇ ⊥ ∇ ⊥ · U + νs Ω ( ∇ × curl θ ) · ∇ ∇ · S h P U (cid:3) , (125c) K is the first-order linear operator K U = ε (cid:2) − ν ∇ × curl( U · ∇ ρ ) + ν ( ∇ × curl U ) · ∇ ρ + ν ∇ × curl ˆΓ + Ω ( Ω · ∇ U ) · ∇ ω − ν ( Ω · ∇ U ) · ∇∇ × curl θ − ν Ω ( Ω · ∇ U ) · ∇ ∆ θ (cid:3) ,(125d)and g = − ∇ × ( ρ k ) + ε (cid:20) ν s − (cid:18) − ∂ z ρ ∇ ¯ u ∇ ρ · ∇ ¯ u (cid:19) + ν ∇ × curl ˆΓ + ν ¯ u · ∇∇ × curl θ + ξ · ∇ ¯ u + 2 λ ∇ ρ × ∇ ∇ ⊥ · U g + 2 λ Ω ∇ ρ · ∇ ⊥ ∇ ⊥ · U g (cid:21) . (125e)We note that S , via the third term in (125c), and K also contain lower order contribu-tions. However, this is not important for the argument that follows and sorting termsstrictly by order leads to even more complicated expressions.
9. Ellipticity of the balance relation for stably stratified flow
The balance relation written in the form (125a) is complicated, but we shall show in thefollowing that U can be recovered from ρ and ¯ u as the solution of an elliptic equationin the case when the flow is strongly and uniformly stratified. More specifically, weshall assume that ∂ z ρ = − α + ∂ z ˜ ρ (126)where ˜ ρ is small in a suitable norm, to be specified further below. It seems possible toweaken this assumption and allow order one variations in the stratification profile solong as stratification is uniformly stable across the layer, but the simpler assumption(126) makes the structure of the problem more transparent.Under this assumption, we can replace the operator A by A = − ( Ω · ∇ ) − ενα ∆ (127a)and S U = S U + εν ∂ z ˜ ρ ∆ U = ε (cid:2) − αν ∇ × curl ¯ U − ω ( Ω · ∇ ) U + ν ∂ z ˜ ρ ∆ U − ν ∇ ˜ ρ · ∇ ∂ z U + ν ∇ × curl ( U · ∇ ˜ ρ )+ ν Ω ∇ ˜ ρ · ∇ ⊥ ∇ ⊥ · U + νs Ω ( ∇ × curl ˜ θ ) · ∇ ∇ · S h P U (cid:3) (127b)25here ˜ θ is defined by (76b) with ˜ ρ in place of ρ and we note that ∇ × curl ˜ θ = ∇ × curl θ .In the following, we will drop the subscript 1 from our notation with the understandingthat we refer to the modified definitions (127).We first note that the operator L defined in (125a) has an infinite-dimensionalkernel which we must quotient out. Indeed, all reference to U in (122), the precursorexpression for (125a), come in either via ˆ u = Ω · ∇ U which is invariant to arbitrarychanges in ¯ U or via ∇ ⊥ · U which is invariant under changes in ∇ · U as well as changesin U so long as these changes are divergence-free due to the relation between ¯ u and ∇ · U established in Lemma 4.12. Thus, we seek solutions to L U = g in the space X s = { U ∈ H s ( D , R ) : U ( − H ) = U (0) = 0 } (128a)modulo functions in Y s = { W ∈ H s ( D , R ) : ∇ ⊥ · W = 0, ∇ · W = 0, and ˆ W = 0 } . (128b)A vector field W ∈ Y s is mean-free and divergence-free, so that ∇ · W = cs ∂ y W . (129)Moreover, its horizontal components are curl-free, so they can be written as the hori-zontal gradient of a potential, i.e., W = ∇ ¯ φ , so that ¯ φ can be obtained by solving thetwo-dimensional elliptic problem with periodic boundary conditions,∆ ¯ φ = cs ∇ ∂ y W . (130)Hence, W = cs ∇ ∆ − ∂ y W . (131)Since (131) is an operator of order zero, the operator P Y U = (cid:18) c/s ∇ ∆ − ∂ y ¯ U ¯ U (cid:19) (132)is bounded on any H s . Further, endowing H s ( D , R ) with the inner product h U , V i H s = h U , V i L + h U , ( − ∆) s V i L + h ( Ω · ∇ ) s U , ( Ω · ∇ ) s V i L (133)which induces a norm which is equivalent to the canonical H s -norm, we see that P Y is the H s -orthogonal projector onto Y s . Noting that the mean and the mean-free partof a vector field in X s are also orthogonal with respect to the inner product (133), wesee that we can decompose any vector field U ∈ X s into three orthogonal components, U = ˆ U + W + (cid:18) ∇ ⊥ ¯ ψ (cid:19) (134)where W ∈ Y s . Thus, solving L U = g in the quotient space X s /Y s amounts to findinga solution to the problem that has W = P Y U = 0.26e recall that L is elliptic if the bilinear form a ( U , V ) = h A U + S U , V i L (135)is coercive in H . (Note that the notion of ellipticity in terms of positivity of the co-efficient matrix of the principal symbol of the operator is not available here because S involves averaging so that S is not in the classical form of a second order differ-ential operator.) A quick inspection reveals that the contribution from A is clearlycoercive. When solving in the quotient X s /Y s , the first term in (127b) is zero. Allother contributions from S are small relative to the contributions from A provided ω is small and bounds on second derivatives of ˜ ρ are small. We conclude that, underthese assumptions, L is elliptic.Similarly, an inspection of K in (125d) reveals that it is a first-order operator (thesecond-order contributions cancel between the first and second term) with coefficientsthat again depend on either ω or ˜ ρ ; if these are small in suitable norms, K is anoperator with small norm from H s +1 into H s . In that case, we can solve the balancerelation in X s /Y s by iteration. Setting Q Y = 1 − P Y , applying Q Y to the balancerelation (125a), and noting that P Y and A commute, we can write A U + Q Y S U + Q Y K U = Q Y g , (136)so that U = A − Q Y g − A − Q Y ( S + K ) U . (137)In conclusion, when the mean flow is small and the density field is a small deviationfrom a uniform stable stratification, the right hand side is contractive and the iterationwill converge.
10. Reconstruction of the full velocity field
The balance relation (136) so far only provides the mean-free component ˆ u of thevelocity. To find the mean velocity ¯ u , we proceed as follows. First, note that theboundary condition k · u = 0 at z = − H , 0 already provides that ¯ u = − ˆ u at z = − H (or z = 0, but one boundary suffices to fully specify ¯ u ). Alternatively, Lemma 4.12provides the expression ¯ u = − s ∇ · U .To find the horizontal mean velocity ¯ u , note that s ω = ∇ ⊥ · M h ¯ u satisfies theevolution equation (117). Moreover, ∇ · M h ¯ u = 0 by Lemma 4.4 since M h = S h P .Noting that ∇ ⊥ · M h ¯ u = ∇ ⊥ · ¯ u − c/s ∂ x ¯ u , ∇ · M h ¯ u = ∇ · ¯ u − c/s ∂ y ¯ u , and writing¯ u = ∇ ⊥ ψ + ∇ φ in terms of its stream function ψ and velocity potential φ , we see thatwe can reconstruct the horizontal mean velocity by solving the two Poisson equations∆ ψ = s ω + cs ∂ x ¯ u = s ω − c ∂ x ∇ · U , (138a)∆ φ = cs ∂ y ¯ u = − c ∂ y ∇ · U . (138b)The full velocity field u = ˆ u + ¯ u can then be used to propagate ρ and the potentialvorticity. 27
1. Discussion and Conclusion
In this paper, we have achieved a variational model reduction for the full three-dimensional Euler–Boussinesq equation with a full Coriolis force. We have studieda simple setting, namely the f -plane approximation, a layer of fluid of constant depth,and periodic boundary conditions in the horizontal. The picture which emerges isstructurally identical to that of variational balance models for the primitive equationsas derived by Salmon (1996) and generalized in Oliver and Vasylkevych (2016):(i) The prognostic variables of the first-order balance model are the tracer field(density or potential temperature) and the vertically averaged relative vortic-ity, here tilted and averaged along the axis of rotation, which depends on thehorizontal spatial coordinates only.(ii) The velocity field is related to the prognostic fields via a kinematic balancerelation. The balance relation is elliptic if rotation is sufficiently fast and theprognostic fields are small perturbations of a stably stratified equilibrium state.(iii) The balance model conserves energy and has a materially conserved potentialvorticity.Thus, we have verified that the derivation of balance models of semigeostrophic typeextends all the way to one of the most general models of geophysical flows. In particu-lar, the assumption of hydrostaticity and of the “traditional approximation” changesdetails, but does not change the structural features of the semigeostrophic limit.At the same time, we note that a full Coriolis force causes difficulties that are notseen in the simpler cases. Since the axis of rotation is not aligned with the directionof gravitational force, there are two distinguished “vertical” directions. This is anobstacle to using a fully intrinsic geometric formulation of the derivation in the spiritof Arnold and Khesin (1999) or Gilbert and Vanneste (2018), forcing us to resort todetailed coordinate calculations. The resulting equations, therefore, appear to lack therelative simplicity of balance models in more idealized settings; many of the new termssimplify or disappear when the Coriolis force is acting exactly in the horizontal plane.We also note that our derivation requires a nontrivial change of coordinates comein already at O ( ε ), where the transformation vector field depends on the prognosticpart of the mean velocity; cf. the last term in (84). We believe that it is not possible toremove this contribution to the transformation at O ( ε ) since there is no leading-orderconstraint through the thermal wind relation on this component. For the analogouscomputation for the primitive equations, it suffices to set λ = in (87), which formallycancels all terms at O ( ε ). For the Euler–Boussinesq system, it is the last term in (84)that cannot be canceled. We note that this additional contribution to the transforma-tion vector field appears even in the case when the axis of rotation is aligned with thegeometric vertical and only disappears when the hydrostatic approximation is made.Thus, the more straightforward derivation by Salmon (1996), who inserts the thermalwind relation directly into the extended Lagrangian to constrain the variational prin-ciple does not work in this case; our more general transformational setting describedin Section 3 cannot be avoided.When the axis of rotation is aligned with the horizontal, i.e., for the equatorial f -plane, all our final expression remain non-singular, which is surprising given thatsome of the intermediate expressions do contain diverging terms. What fails, however,is ellipticity of the balance relation. Stratification is providing regularization in thehorizontal plane, rotation is providing regularization in the vertical. When the axis ofrotation is tilted into the horizontal, all control in the geometric vertical is lost. We28o not expect that the model maintains any sense of well-posedness in this case.Whether the balance model is well-posed as an initial value problem on the mid-latitude f -plane is not obvious as the right hand side of the balance relation containshigh-order derivatives of the prognostic variables. All we can assert is that the bal-ance relation itself allows, under the conditions stated, stable reconstruction of theslaved components of the velocity field from sufficiently smooth prognostic variables.In practical terms, the question of solving the balance model as a prognostic modelis probably not relevant, as the complexity-accuracy tradeoff is not in its favor. Thebalance relation itself, however, might be useful as a diagnostic relation. Moreover,the structure asserted here gives definite proof that the concept of balance does notrequire conceptual modification if the assumption of hydrostaticity or exact verticalalignment of the axis of rotation are relaxed. Acknowledgment
This paper is a contribution to project M2 (Systematic Multi-Scale Analysis and Mod-eling for Geophysical Flow) of the Collaborative Research Center TRR 181 “EnergyTransfer in Atmosphere and Ocean” funded by the Deutsche Forschungsgemeinschaft(DFG, German Research Foundation) under project number 274762653. Additionalfunding was received through the Ideen- und Risikofund 2020 at Universit¨at Ham-burg.
References
Allen, J. S., Holm, D. D., & Newberger, P. A. (2002). Toward an extended-geostrophic Euler–Poincar´e model for mesoscale oceanographic flow. In J. Norbury & I. Roulstone (Eds.),
Large-scale atmosphere–ocean dynamics (Vol. 1, pp. 101–125). Cambridge University Press.Arnold, V. I., & Khesin, B. A. (1999).
Topological methods in hydrodynamics . Springer.Dellar, P. J. (2011). Variations on a beta-plane: derivation of non-traditional beta-planeequations from Hamilton’s principle on a sphere.
J. Fluid Mech. , , 174–195.Dritschel, D. G., Gottwald, G. A., & Oliver, M. (2017). Comparison of variational balancemodels for the rotating shallow water equations. J. Fluid Mech. , , 689–716.Ebin, D. G., & Marsden, J. (1970). Groups of diffeomorphisms and the motion of an incom-pressible fluid. Ann. of Math. (2) , , 102–163.Franzke, C. L. E., Oliver, M., Rademacher, J. D. M., & Badin, G. (2019). Multi-scale methodsfor geophysical flows. In C. Eden & A. Iske (Eds.), Energy Transfers in Atmosphere andOcean (pp. 1–51). Springer-Verlag.Gerkema, T., & Shrira, V. I. (2005a). Near-inertial waves in the ocean: beyond the ‘traditionalapproximation’.
J. Fluid Mech. , , 195–219.Gerkema, T., & Shrira, V. I. (2005b). Near-inertial waves on the “nontraditional” β -plane. J.Geophys. Res. Oceans , (C1).Gerkema, T., Zimmerman, J. T. F., Maas, L. R. M., & van Haren, H. (2008). Geophysical andastrophysical fluid dynamics beyond the traditional approximation. Rev. Geophys. , (2),RG2004.Gilbarg, D., & Trudinger, N. S. (2001). Elliptic partial differential equations of second order .Springer.Gilbert, A. D., & Vanneste, J. (2018). Geometric generalised Lagrangian-mean theories.
J.Fluid Mech. , , 95–134.Gottwald, G. A., & Oliver, M. (2014). Slow dynamics via degenerate variational asymptotics. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. , (2170), 20140460. olm, D. D., Marsden, J. E., & Ratiu, T. S. (1998). The Euler–Poincar´e equations andsemidirect products with applications to continuum theories. Adv. Math. , (1), 1–81.Oliver, M. (2006). Variational asymptotics for rotating shallow water near geostrophy: atransformational approach. J. Fluid Mech. , , 197–234.Oliver, M., & Vasylkevych, S. (2013). Generalized LSG models with spatially varying Coriolisparameter. Geophys. Astrophys. Fluid Dyn. , , 259–276.Oliver, M., & Vasylkevych, S. (2016). Generalized large-scale semigeostrophic approximationsfor the f -plane primitive equations. J. Phys. A: Math. Theor. , , 184001.Palais, R. S. (1968). Foundations of Global Non-Linear Analysis . W. A. Benjamin, Inc., NewYork-Amsterdam.Pedlosky, J. (1987).
Geophysical fluid dynamics (Second ed.). Springer.Salmon, R. (1983). Practical use of Hamilton’s principle.
J. Fluid Mech. , , 431–444.Salmon, R. (1985). New equations for nearly geostrophic flow. J. Fluid Mech. , , 461–477.Salmon, R. (1988). Semigeostrophic theory as a Dirac-bracket projection. J. Fluid Mech. , , 345–358.Salmon, R. (1996). Large-scale semigeostrophic equations for use in ocean circulation models. J. Fluid Mech. , , 85–105.Stewart, A. L., & Dellar, P. J. (2010). Multilayer shallow water equations with completeCoriolis force. Part 1. Derivation on a non-traditional beta-plane. J. Fluid Mech. , ,387.Stewart, A. L., & Dellar, P. J. (2012). Multilayer shallow water equations with completeCoriolis force. Part 2. Linear plane waves. J. Fluid Mech. , , 16–50.Vallis, G. K. (2006). Atmospheric and oceanic fluid dynamics: fundamentals and large-scalecirculation . Cambridge University Press.White, A. A. (2002). A view of the equations of meteorological dynamics and various ap-proximations. In J. Norbury & I. Roulstone (Eds.),
Large-scale atmosphere–ocean dynamics (Vol. 1, pp. 1–100). Cambridge University Press.(Vol. 1, pp. 1–100). Cambridge University Press.