Variational integrators for anelastic and pseudo-incompressible flows

The anelastic and pseudo-incompressible equations are two well-known soundproof approximations of compressible flows useful for both theoretical and numerical analysis in meteorology, atmospheric science, and ocean studies. In this paper, we derive and test structure-preserving numerical schemes for these two systems. The derivations are based on a discrete version of the Euler-Poincar\'e variational method. This approach relies on a finite dimensional approximation of the (Lie) group of diffeomorphisms that preserve weighted-volume forms. These weights describe the background stratification of the fluid and correspond to the weighed velocity fields for anelastic and pseudo-incompressible approximations. In particular, we identify to these discrete Lie group configurations the associated Lie algebras such that elements of the latter correspond to weighted velocity fields that satisfy the divergence-free conditions for both systems. Defining discrete Lagrangians in terms of these Lie algebras, the discrete equations follow by means of variational principles. Descending from variational principles, the schemes exhibit further a discrete version of Kelvin circulation theorem, are applicable to irregular meshes, and show excellent long term energy behavior. We illustrate the properties of the schemes by performing preliminary test cases.


Introduction
Numerical simulations of atmosphere and ocean on the global scale are of high importance in the field of Geophysical Fluid Dynamics (GFD). The dynamics of these systems are frequently modeled by the full Euler equations using explicit time integration schemes (see, e.g., [6]). These simulations are however computationally very expensive. Besides highly resolved meshes to capture important small scale features, the fast traveling sound waves have to be resolved too, by very small time step sizes, in order to guarantee stable simulations [6]. As these sound wave are assumed to be negligible in atmospheric flows, soundproof models, in which these fast waves are filtered out, are a viable option that permits to increase the time step sizes and hence to speed up calculations significantly.
Frequently applied soundproof models are the Boussinesq, anelastic, pseudo-incompressible approximations of the Euler equations [5,11,13]. There exist elaborated discretizations of these equations in literature. However, these discretizations often do not take into account the underlying geometrical structure of the equations. This may result in a lack of conserving mass, momentum, energy, or to the fact that the Helmholtz decomposition of vector fields or the Kelvin-Noether circulation theorem are not satisfied. Structure-preserving schemes descending from Euler-Poincaré variational methods [14], [7], [4] conserve these quantities, as they arise from a Lagrangian formulation, in which these conserved quantities are given by invariants of the Lagrangian under symmetries, [8].
With this paper, we contribute to develop variational integrators in the area of GFD by including the anelastic and pseudo-incompressible schemes into the variational discretization framework developed by [14]. To use this framework, we first have to describe the approximation of the Euler equations for a perfect fluid in terms of the Euler-Poincaré variational method [8]. The central idea is to use volume forms that are weighted by the corresponding background stratifications such that they match the divergence-free conditions of the correspondingly weighed velocity fields associated to either the anelastic or the pseudo-incompressible approximations. This will allow us to identify for these approximations the appropriate Lie group configuration with corresponding Lie algebras. Using the latter to define appropriate Lagrangian, the equations of motion follow by Hamilton's variational principle of stationary action.
The definition of appropriate discrete diffeomorphism groups will be based on the idea to use weighted meshes that provide discrete counterparts of the weighted volume forms. The corresponding discrete Lie algebras will incorporate the required divergence-free conditions on the weighed velocity fields. Defining appropriate weighted pairings required to derive the functional derivatives of the discrete Lagrangians, the flat operator introduced in [14] is directly applicable and we can thus avoid to discuss this otherwise delicate issue. Mimicking the continuous theory, the discretizations of anelastic and pseudo-incompressible equations follow by variations of appropriate discrete Lagrangians.
We structure the paper as follows. In Section 2 we recall the standard formulations of Boussinesq, anelastic, and pseudo-incompressible approximations of the Euler equations for perfect fluids. In Section 3 we show that these formulations follow for appropriate Langrangians from the Euler-Poincaré variational principle. In Section 4 we recall the variational discretization framework introduced by [14] and extend it to suit anelastic and pseudo-incompressible equations. In Section 5 we derive corresponding discretizations on 2D simplicial meshes and provide explicit formulations in terms of velocity and buoyancy or potential temperature fields. In Section 6 we perform numerical tests to demonstrate the structure-preserving properties of the variational integrators. In Section 7 we draw conclusions and provide an outlook.

Anelastic and pseudo-incompressible systems
In this section we review the three approximations of the Euler equations of a perfect gas that will be the subject of this paper, namely, the Boussinesq, the anelastic, and the pseudoincompressible approximations (see, e.g., [6] for more details).
The Euler equations for the inviscid isentropic motion of a perfect gas can be expressed in the form ∂ t u + u · ∇u + 1 ρ ∇p = −gz, ∂ t ρ + div(ρu) = 0, ∂ t θ + u · ∇θ = 0, (2.1) where u is the three-dimensional velocity vector, ρ is the mass density, p is the pressure, g is the gravitational acceleration, z is the unit vector directed opposite to the gravitational force. The variable θ is the potential temperature, defined by θ = T /π, in which T is the temperature and π is the Exner pressure π = (p/p 0 ) R/cp , with R is the gas constant for dry air, c p is the specific heat at constant pressure, and p 0 is a constant reference pressure. Using the equation of state for a perfect gas, p = ρRT , we have the relation 1 ρ ∇p = c p θ∇π.
Boussinesq approximation. This approximation is obtained by assuming a nondivergent flow and by neglecting the variations in potential temperature except in the leading-order contribution to the buoyancy. We thus get, from (2.3), the system ∂ t u + u · ∇u + c p θ 0 ∇π = g θ θ 0 z, div u = 0, ∂ t θ + u · ∇θ = 0, in which θ 0 is a constant reference potential temperature. These equations can equivalently be written as where P b = c p θ 0 π , N 2 = g θ 0 ∂ zθ is the Brunt-Väisälä frequency, and b = g θ θ 0 is the buoyancy. Making use of the full buoyancy b = g θ θ 0 = gθ +θ θ 0 , we can write the system as The total energy is conserved since the energy density E = 1 2 |u| 2 − bz = 1 2 |u| 2 − g θ θ 0 z verifies the continuity equation ∂ t E + div((E + P b )u) = 0. (2.5) The requirement for nondivergent flow is easily justified only for liquids, and the errors incurred approximating the true mass conservation relation by div u = 0 can be quite large in stratified compressible flows. In this case, the anelastic and pseudo-incompressible models have to be considered, which better approximate the true mass continuity equation.
Anelastic approximation. The anelastic system approximates the continuity equation as div(ρu) = 0, whereρ(z) is the vertically varying density of the reference state.
The resulting system is however not energy conservative. In order to restore energy conservation, [11] considered the approximate momentum equation In this case, the energy density E =ρ 1 2 |u| 2 + c pπ θ , whereπ(z) is such that c p ∂π ∂z = − ḡ θ , satisfies the continuity equation with P a :=ρ(c pθ π + gz).

Variational formulation
We shall now formulate the anelastic and pseudo-incompressible equations in Euler-Poincaré variational form. Euler-Poincaré variational principles are Eulerian versions of the classical Hamilton principle of critical action. We refer to [8] for the general Euler-Poincaré theory based on Lagrangian reduction and for several applications in fluid dynamics. An Euler-Poincaré formulation for anelastic systems was given in [3]. We shall develop below a slightly different Euler-Poincaré approach, well-suited for the variational discretization, by putting the emphasis on the underlying Lie group of diffeomorphisms associated to these systems. As we have recalled above, the anelastic and pseudo-incompressible equations are based on a constraint of the following type on the fluid velocity u(t, x): div(σu) = 0, for a given strictly positive functionσ(x) > 0 on the fluid domain D.
We assume that the fluid domain D is a compact, connected, orientable manifold with smooth boundary ∂D. In our examples, D is a 2D domain in the vertical plane R 2 x = (x, z) or a 3D domain in R 3 x = (x, y, z).
We fix a volume form µ on D, i.e., an n-form, n = dim D, with µ(x) = 0, for all x ∈ D . If D is a domain in R 3 (x, y, z), one can take µ = dx ∧ dy ∧ dz to be the standard volume of R 3 restricted to D. We shall denote by div µ (u) the divergence of u with respect to the volume form µ. Recall that the divergence is the function div µ (u) defined by the equality in which £ u denotes the Lie derivative with respect to the vector field u, see, e.g., [1]. When µ is the standard volume, one evidently recovers the usual divergence operator div on vector fields.
Diffeomorphism groups. Let us denote by Diff µ (D) the group of all smooth diffeomorphisms ϕ : D → D that preserve the volume form µ, i.e., ϕ * µ = µ. The group structure on Diff µ (D) is given by the composition of diffeomorphisms. The group Diff µ (D) can be endowed with the structure of a Fréchet Lie group, although in this paper we shall only use the Lie group structure at a formal level. The Lie algebra of the group Diff µ (D) is given by the space X µ (D) of all divergence free (relative to µ) vector fields on D, parallel to the boundary ∂D: Given the strictly positive functionσ > 0 on D, we consider the new volume formσµ with associated diffeomorphism group and Lie algebra denoted Diffσ µ (D) and Xσ µ (D) = {u ∈ X(D) | divσ µ (u) = 0, u ∂D}, respectively. In the next Lemma, we rewrite the condition divσ µ (u) = 0 by using exclusively the divergence operator div µ associated to the initial volume form µ.
Lemma 3.1 Let D be a manifold endowed with a volume form µ and letσ > 0 be a strictly positive smooth function on D. Then we have div µ (σu) =σ divσ µ (u).
Proof. We will use the following properties of the Lie derivative £ u , the exterior differential d, and the inner product i u on differential forms (see, e.g., [1]): for a k-form α, an n-form β, and a vector field u, we have On the one hand, we have On the other hand, we havē This proves the result.
From this Lemma, we deduce that the appropriate Lie groups associated to the anelastic and pseudo-incompressible systems are given by G = Diffρ µ (D) and G = Diffρθ µ (D), respectively. Indeed, from the preceding Lemma, it follows that the Lie algebras of these groups can be written as respectively. They correspond to the anelastic and pseudo-incompressible constraints on the fluid velocity. We will continue to uses the subscriptσµ when referring to both Lie groups and both Lie algebras.
Euler-Poincaré variational principles. The diffeomorphism group Diffσ µ (D) plays the role of the configuration manifold for these fluid models. The motion of the fluid is completely characterized by a time dependent curve ϕ(t, ) ∈ Diffσ µ (D): a particle located at a point X ∈ D at time t = 0 travels to x = ϕ(t, X) ∈ D at time t. Exactly as in classical mechanics, the Lagrangian of the system is defined on the tangent bundle T Diffσ µ (D) of the configuration manifold. We shall denote it by L Θ 0 : T Diffσ µ (D) → R. The index Θ 0 indicates that this Lagrangian parametrically depends on the potential temperature Θ 0 (X) that is expressed here in the Lagrangian description.
The equations of motion in the Lagrangian description follow from the Hamilton principle over a time interval [0, T ], for variations δϕ with δϕ(0) = δϕ(T ) = 0.
In the Eulerian description, the variables are the Eulerian velocity u(t, x) and the potential temperature θ(t, x). They are related to ϕ(t, X) and Θ 0 (X) as u(t, ϕ(t, X)) =φ(t, X) and θ(t, ϕ(t, X)) = Θ 0 (X). (3.2) We assume that the Lagrangian L Θ 0 can be rewritten exclusively in terms of these two Eulerian variables, and we denote it by (u, θ). This assumption means that L Θ 0 is right-invariant with respect to the action of the subgroup of all diffeomorphisms that keep Θ 0 invariant. By rewriting the Hamilton principle (3.1) in terms of the Eulerian variables u and θ, we get the Euler-Poincaré variational principle where v(t, x) is an arbitrary vector field on D parallel to the boundary and with div(σv) = 0, (i.e., v ∈ Xσ µ (D) by Lemma 3.1), and with v(0, x) = v(T, x) = 0. The bracket [u, v], locally given by [u, v] The expressions for δu and δθ in (3.3) follow by taking the variation of the first and second equalities in (3.2) and defining v(t, x) as v(t, ϕ(t, X)) = δϕ(t, X). A direct and efficient way to obtain these expressions, or the variational principle (3.3), is to apply the general theory of Euler-Poincaré reduction on Lie groups, see [8].
In order to compute the associated equations, one needs to fix an appropriate space in nondegenerate duality with Xσ µ (D). This is recalled in the next Lemma, which follows from the Hodge decomposition and shall play a crucial role in the discrete setting later. Recall that if V is a vector space, a space in nondegenerate duality with V is a vector space V to together with a bilinear form , : V × V → R such that α, v = 0, for all v ∈ V , implies α = 0 and α, v = 0, for all α ∈ V , implies v = 0.
Lemma 3.2 The space Ω 1 (D)/dΩ 0 (D) of one-forms modulo exact forms is in nondegenerate duality with the space Xσ µ (D), the Lie algebra of Diffσ µ (D). The nondegenerate duality pairing is given by
Proof. It is well-known that if g is a Riemannian metric, with µ g the associated volume form on D, then , : is a nondegenerate duality pairing, see e.g., [12, §14.1]. This result follows from the Hodge decomposition of 1-forms, which needs the introduction of a Riemannian metric g. In our case, the volume forms µ andσµ are not necessarily associated to a Riemannian metric. We shall thus introduce a Riemannian metric g uniquely for the purpose of this proof, with associated Riemannian volume form µ g . Let f be the function defined byσµ = f µ g . Since D is orientable and connected, we have either f > 0 or f < 0 on D. We can rewrite the duality pairing (3.4) as (3.5) By successive applications of Lemma 3.1, we have where the last equality follows since u ∈ Xσ µ (D). This proves that v = f u ∈ X µg (D). We can thus write the duality pairing , σ in terms of the nondegenerate duality pairing (3.5) as [α], u σ = [α], f u , which proves that it is nondegenerate.
In a similar way to (3.4), we shall identify the dual to the space of functions F(D) with itself by using the nondegenerate duality pairing ∈ Ω 1 (D)/dΩ 0 (D), for δ δu ∈ Ω 1 (D), and δ δθ ∈ F(D).

Proposition 3.3
The variational principle (3.3) yields the partial differential equation where £ u denotes the Lie derivative acting on one-forms, given by £ u α = d(i u α) + i u dα. This equation is supplemented with the advection equation which follows from the definition of θ in (3.2).
Proof. By definition of the functional derivatives, we have Using the expression for δu in (3.3), integrating by parts, and using the equalities We can write d δ Combining these results, we thus get for all v ∈ Xσ µ (D). By Lemma 3.2, it follows that the one-form ∂ t δ δu + £ u δ δu + δ δθ dθ is exact, i.e., there exists a function p such that this expression equals dp. Remark 3.4 We note that the statement of Proposition 3.3 does not need the introduction of a Riemannian metric g on D. Only a volume form µ is fixed, together with a strictly positive functionσ. It can be however advantageous to formulate the equations (3.7) in terms of a Riemannian metric g (note that we do not suppose that µ orσµ equals µ g ). In this case, identifying one-forms and vector fields via the flat operator u ∈ X(D) → u = g(u, ) ∈ Ω 1 (D), the space Ω 1 (D)/dΩ 0 (D) can be identified with the space of vector fields X(D) modulo gradient (with respect to g) of functions. The nondegenerate duality pairing (3.4) thus reads In terms of this duality pairing, the equations (3.7) are equivalently written as where ∇ acting on a vector field is the covariant derivative associated to the Riemannian metric g, ∇ acting on a function is the gradient relative to g, and ∇u T denotes the transpose with respect to g.
We shall now apply this setting to the anelastic and the pseudo-incompressible equations. The fluid domain D is a subset of the vertical plane R 2 (x, z) or of the space R 3 (x, y, z), and has a smooth boundary ∂D. We fix a volume form µ on D.
1) Anelastic equations. For the anelastic equation, we takeσ(z) =ρ(z), the reference mass density. The Lagrangian is given by (3.10) whereπ(z) is such that c p ∂π ∂z = − ḡ θ and the norm is computed relative to the standard inner product on R 2 or R 3 .
Relative to the pairings (3.8) and (3.6) we get δ δu = u and δ δθ = −c pπ , (3.11) so that the Euler-Poincaré equations (3.9) read ∂ t u + u · ∇u + ∇u T · u − c pπ ∇θ = −∇p, in terms of the pressure p. To permit a comparison of these anelastic equations with those given in the standard form of (2.7) in terms of Exner pressure π , we note that ∇u T · u = 1 2 ∇|u| 2 and that −c pπ ∇θ differs from −g θ θ z by a gradient term, indeed: Therefore, with π defined in terms of p by the equality c pθ π = p + 1 2 |u| 2 − gz − c pπ θ, the Euler-Poincaré equations yield the anelastic equations (2.7).
2) Pseudo-incompressible equations. In this case we takeσ(z) :=ρ(z)θ(z) and the Lagrangian is given by As before, the kinetic energy is computed relative to the standard inner product on R 2 or R 3 . Relative to the pairings (3.8) and (3.6) we get so that the Euler-Poincaré equations (3.9) read After some computations, using the relation − ḡ θ = c p ∂ zπ , these equations recover the pseudoincompressible system (2.8) with c p (π + π ) = p + 1 θ ( 1 2 |u| 2 − gz). Based on these results, we can formulate the following statement that will allow us to derive the variational discretization of these two models by the discrete diffeomorphism group approach.
Kelvin-Noether circulation theorems. The Euler-Poincaré formulation is well adapted for a systematic derivation of the circulation theorems, see [8]. From (3.7) one indeed deduces the following general form of the circulation theorem where c t = ϕ(t, c 0 ) is a loop advected by the fluid flow ϕ(t, ) and ct α denotes the circulation of the one-form α along the loop c t . Using the equation ∂ t θ + dθ · u = 0, one also deduces another useful form, namely, We shall not present the derivation of (3.15) and (3.16) since they follow from similar arguments with those explained in details in [8]. For the anelastic system, using (3.11) and the equalities of the circulation theorem. For the pseudo-incompressible system, using (3.13) and the equali- As shown further below, these conservation laws of the continuous equations, here presented with explicit formulas, are also preserved by the discrete variational discretizations that we will derive in the following section.

Variational discretizations
In this section we first quickly review from [14] the discrete diffeomorphism group approach in the incompressible case. Then, based on the results of Theorem 3.5, we show that an appropriate adaptation of this approach allows us to derive a variational discretization of the anelastic and pseudo-incompressible systems valid on a large class of mesh discretizations of the fluid domain.
Review of the discrete diffeomorphism group approach in the incompressible case.
The spatial discretization of the equations is accomplished by considering the finite dimensional approximation of the group of volume preserving diffeomorphisms developed in [14], which we roughly recall below. Given a mesh M on the fluid domain D with cells C i , i = 1, ..., N , define a diagonal N × N matrix Ω consisting of cell volumes: Ω i = Vol(C i ). The discretization of the group Diff µ (D) of volume preserving diffeomorphisms of D is the matrix group where GL(N ) + is the group of invertible N × N matrices with positive determinant, and 1 denotes the column (1, ..., 1) T so that the first condition reads N j=1 q ij = 1 for all i = 1, ..., N . The main idea behind this definition is the following (see [14] for the detailed treatment). Consider the linear action of Diff µ (D) on the space F(D) of functions on D, given by This linear map has two key proporties: (1) it preserves the L 2 inner product of functions; (2) it preserves constant functions C on D: In the discrete setting, a function is discretized as a vector F ∈ R N whose value F i on cell C i is regarded as the cell average of the function. Accordingly, the discrete L 2 inner product of two discrete functions is defined by The discrete diffeomorphism group (4.1) is such that its action on discrete functions by matrix multiplication, is an approximation of the linear map (4.2). The conditions q T Ωq = Ω and q · 1 = 1 are indeed the discrete analogues of the conditions (1) and (2) above, respectively. The Lie algebra of D(M), denoted d(M), is the space of Ω-antisymmetric, row-null matrices: The component A ij of the matrix A is the weighted flux of the vector field u through the face common to the cells C i and C j . This relation induces a nonholonomic constraint on the Lie algebra d(M) as only the fluxes through adjacent cells are non-zero: in which N (i) is the set of all indices of cells adjacent to cell C i . Once a discrete Lagrangian d : d(M) → R has been selected, the derivation of the spatial variational discretization then proceeds by applying an Euler-Poincaré variational principle to this Lagrangian that takes into account the nonholonomic constraint. This approach has been developed in [14] for the incompressible homogenous ideal fluid and extended to several models of incompressible fluids with advection equations in [7] and to rotating and/or stratified Boussinesq flows in [4].
Discrete diffeomorphism groups for the two models. The results obtained in Theorem 3.5 make it possible to extend this approach to treat the anelastic and pseudo-incompressible systems. The main idea consists in defining weighted versions of the volume of the cells in order to permit the use of the results recalled above in the incompressible case.
Given a mesh M on D, a reference densityρ(z), and a reference potential temperatureθ(z) on D, we define, respectively, the diagonal matrices Ωρ and Ωρθ ofρ -weighted andρθ -weighted volumes as The discrete versions of the diffeomorphism groups Diffρ µ (D) and Diffρθ µ (D) in (3.14) are therefore Dρ(M) : = q ∈ GL(N ) + | q · 1 = 1 and q T Ωρq = Ωρ , with Lie algebras In order to treat the anelastic and pseudo-incompressible cases, we also need to appropriately modify the relation between the components A ij and the velocity vector fields u, by taking into account the weightsρ andθ. For u ∈ Xρ µ (D), resp., u ∈ Xρθ µ (D), we get in which D ij is the boundary common to C i and C j , and n ij is the normal vector field on D ij pointing from C i to C j . The same nonholonomic constraint as before needs to be imposed, but this time on dσ(M), namely, This constraint induces a right-invariant linear constraint on the Lie group Dσ(M): at q ∈ Dσ(M), the constraint is defined as In addition to the matrix A ∈ dσ(M) which discretizes the Eulerian velocity u, we introduce the discrete potential temperature Θ ∈ R N whose component Θ i is the average of the potential temperature θ on cell C i .
Variational discretization for the two models. The variational discretization is carried out by mimicking the Euler-Poincaré approach of Theorem 3.5. Consider a discrete Lagrangian L Θ 0 ,d : T Dσ(M) → R defined on the tangent bundle of the Lie group Dσ(M) and being an approximation of the Lagrangian in (3.1). The parameter Θ 0 ∈ R N is the discrete potential temperature in the Lagrangian description.
Hamilton's principle has to be appropriately modified to take into account the nonholonomic constraint, namely, we apply the Lagrange-d'Alembert principle stating that the action functional is critical with respect to variations subject to the constraint. In our case it reads vanishing at t = 0, T , and withq ∈ Sσ(q).
In a similar way with the continuous case, the relation between the Lagrangian variables q(t) ∈ Dσ(M), Θ 0 ∈ R N , and the Eulerian variables A(t) ∈ dσ(M), Θ(t) ∈ R N , is given by the formulas The discrete Lagrangian L Θ 0 ,d is assumed to have the same right-invariance as its continuous counterpart in (3.1), hence it can be exclusively written in terms of the Eulerian variables in (4.7). We thus get the reduced Lagrangian The Eulerian version of the Lagrange-d'Alembert principle (4.6) is found to be in which A ∈ Sσ and where Y is an arbitrary time dependent matrix in Sσ vanishing at the endpoints. The expressions for the variations δA and δΘ in (4.8) are obtained by using the two relations (4.7). In particular, we have Y = δqq −1 , which is therefore an arbitrary time dependent matrix in Sσ vanishing at t = 0, T . The principle (4.8) is a nonholonomic version of the Euler-Poincaré variational principle.
In order to derive the equations associated to the principle (4.8), we need to introduce appropriate dual spaces to the Lie algebra dσ(M) and the space of discrete functions Ω 0 d (M) = R N . To this end, we recall from [7] that in the context of discrete diffeomorphism groups, a discrete one-form on M is identified with a skew-symmetric N × N matrix. The space of discrete one-forms is denoted by Then, given a strictly positive functionσ(x) > 0, the discrete version of the L 2 pairing (3.4) is defined as By repeating the arguments of Theorem 2.4 of [7] for the discrete L 2 pairing (4.9) with weight σ, we get the identification which is the discrete analogue of the identification in Lemma 3.2.
Concerning functions, the discrete analogue of the pairing (3.6) is given by A direct application of the principle (4.8) yields the following result.
Anelastic system. The discrete Lagrangian associated to (3.10) is where the first, resp., the second duality pairing is given in (4.9), resp., (4.11), andΠ ∈ R N is a discretization of the reference valueπ(z) of the Exner pressure. The first term in (4.14) is the discretization of the kinetic energy associated to a given Riemannian metric on D and is based on a suitable flat operator A ∈ Sρ → A ∈ Ω 1 d (M) associated to the mesh M, see [14]. The functional derivatives of d with respect to the pairings , ρ and , ρ are, respectively, Using them in (4.12) withσ =ρ, we get the structure-preserving spatial discretization of the anelastic system on the mesh M as Pseudo-incompressible system. The discrete Lagrangian associated to (3.12) is in which Z ∈ R N is a discretization of the height z. Note that we are now using the volumes Ωρθ i . The functional derivatives of d with respect to the pairings , ρθ and , ρθ are, respectively, Using them in (4.12) withσ =ρθ, we get the structure-preserving spatial discretization of the pseudo-incompressible system on the mesh M as (4.17) Figure 5.1: Notation and indexing conventions for the 2D simplicial mesh.

Numerical scheme on irregular simplicial meshes
In this section, we shall use the general results of §4, valid for any kind of reasonable (i.e. non-degenerated) meshes, to deduce the variational discretization on 2D simplicial meshes. On such meshes, we adopt the following notations (cf. Figure 5.1): f ij : = length of the primal edge, located between triangle i and triangle j; h ij : = length of the dual edge, between triangle i and triangle j; The flat operator on a 2D simplicial mesh is defined by the following two conditions, see [14], in which e denotes the node common to triangles T i , T j , T k and ζ 2 e denotes the dual cell to e. In (5.1), the vorticity ω(K) of a discrete one-form K ∈ Ω 1 d (M) is defined by where the sum is taken over the dual edges in the boundary ∂ζ 2 e counterclockwise around node e. The constant K e i is defined as where |ζ 2 e | and |ζ 2 e ∩ T k | denote, respectively, the areas of ζ 2 e and ζ e ∩ T k . Note that the matrix A defined in (5.1) is skew-symmetric, hence A ∈ Ω 1 d (M).

Boussinesq flow
Variational discretization of the Boussinesq fluid on regular Cartesian grids has been carried out in [4]. Here we shall derive from (4.12) the variational scheme on irregular 2D simplicial grids. Recall that in this case div µ (u) = 0 and that the Boussinesq Lagrangian is given by The discrete Lagrangian is therefore chosen as d : where B ∈ R N is the discrete buoyancy and Z is the discrete height function, i.e., Z i denotes the z-coordinate of the center of cell C i . Using the Boussinesq Lagrangian and the flat operator (5.1), the discrete Euler-Poincaré equation (4.12) yields where Ω i A ij = −Ω j A ji , for all i, j, and j∈N (i) A ij = 0, for all i, and whereP is related to P in (4.12) viaP We note that the momentum equation (5.2) corresponds to the discretization of the following form of the Boussinesq equation: where, similarly to (5.3),p = i u u +p, with p the pressure function arising in the Euler-Poincaré formulation (3.7). The form (5.4) is easily seen to be equivalent to the standard form (2.4) with Recall that the Lie algebra element A ∈ d(M) and the velocity u are related as We introduce the discrete velocity V ij such that We note that since A · 1 = 0, we get j∈N (i) f ij V ij = 0, which corresponds to the divergence free condition. We used everywhere the index b, referring to Boussinesq, to distinguish from the expressions occurring further below for the other two models.

Anelastic flow
The anelastic Lagrangian is given in (3.10) and the discrete Lagrangian d : dρ(M) × R N → R in (4.14). Recall that in this case div µ (ρu) = 0. The flat operator (5.1) has to be slightly modified in order to produce a skew-symmetric matrix, namely, we modify the first line in (5.1) to in which (·) (A) denotes the skew-symmetric part. For the Boussinesq model, this definition recovers (5.1), since the matrix M is in this case already skew-symmetric. One checks that this definition still satisfies the properties of a flat operator in [14]. The general discrete anelastic equations (4.15) yield where Ωρ i A ij = −Ωρ j A ji , for all i, j, and j∈N (i) A ij = 0, for all i, and whereP is related to P in (4.12) and (4.15) as before via the formula (5.3). We note that the momentum equation (5.7) corresponds to the discretization of the following form of the anelastic equation: where, similarly to (5.3),p = i u u +p, with p the pressure function arising in the Euler-Poincaré formulation (3.7). The form (5.8) was shown in §3 to be equivalent to the standard form (2.6).
Recall that the Lie algebra element A ∈ dρ(M) and the velocity u are related as We thus introduce the velocity V ij such that A ij = − 1 2Ωρ i f ijρij V ij , in whichρ ij :=ρ i +ρ j 2 , with ρ i := Ωρ i /Ω i the mean value ofρ(z) over cell C i . We get the relation In terms of V , the semidiscrete anelastic momentum equation (5.7) be- where we defined

Pseudo-incompressible flow
The pseudo-incompressible Lagrangian is given in (3.12) and the discrete Lagrangian d : dρ(M) × R N → R is given in (4.16). Recall that in this case div µ (ρθu) = 0. We take the flat operator (5.1) with the first line modified as in (5.6) The general discrete pseudo-incompressible equations (4.17) yield 10) where Ωρθ i A ij = −Ωρθ j A ji , for all i, j, and ∈N (i) A ij = 0, for all i, and whereP is related to P in (4.12) and (4.17) by the fomulã We note that the momentum equation (5.10) corresponds to the discretization of the following form of the pseudo-incompressible equation: where, similarly to (5.11),p = 1 θ i u u + p, with p the pressure function arising in the Euler-Poincaré formulation (3.7). The form (5.12) was shown in §3 to be equivalent to the standard form (2.8).
Recall that the Lie algebra element A ∈ dρθ(M) and the velocity u are related as We thus introduce the velocity V ij such that , with ρθ i := Ωρθ i /Ω i the mean value ofρ(z)θ(z) over cell C i . We get the relation In terms of V , the semidiscrete pseudo-incompressible momentum equation (5.10) becomes where we defined

Time integration and Poisson solvers
Since the spatial discretization has be realized in a structure-preserving way, a corresponding temporal variational discretization follows by applying the general discrete (in time) Euler-Poincaré approach, as it has be done in [7], [4] to which we refer for a detailed treatment. Hence, the temporal scheme consists of the following two steps. Given the time t and a time step size ∆t, we first compute Θ t+1 = τ (hA t )Θ t (with Θ = B in case of Boussinesq), in which τ is the Cayley transform. To this end we represent τ by the equation: (5.14) with I the identity matrix (cf. [4] for more details). In the second step, the following update equations are solved. They are given in explicit form for the Boussinesq/anelastic schemes by and for the pseudo-incompressible equation by This time update scheme, derived from discrete Euler-Poincaré equations, is analogous to a Crank-Nicolson time update scheme. To solve these implicit equations, we proceed as follows. Noting that Θ t+1 has been calculated first and is thus known, we rearrange the equations into known and unknown parts as and similarly for the pseudo-incompressible equations: We can summaries the edge values into vector arrays, which we will indicate by omitting the notation ij. Then, we solve these vector equations for Boussinesq, anelasic, and pseudoincompressible schemes by a fixed point iteration, given by the following algorithm: 1. Start loop over k = 0 with initial guess at t: V * k=0 = V t and P * k=0 = 0 2. Calculate tentative velocity V * k : • Boussinesq/anelastic: 3. Calculate pressure P * k+1 by solving the Poisson equation: 4. Update velocity and set k + 1 = k: • Boussinesq/anelastic: Substituting the tentative velocity V * k into the velocity update equation (step 4), we realize that in case of convergence, i.e. V * k+1 → V t+1 and P * k+1 → P t+1 , this algorithm indeed solves equations (5.17) and (5.18).
Being derived both spatially and temporally via the Euler-Poincaré approach, the schemes (5.15) and (5.16) admit discrete version of the Kelvin circulation theorem. We refer to [7] and [4] for detailed discussions of the discrete Kelvin circulation theorems and its application to several examples. Since the application to the present models is similar we shall not repeat it here.

Numerical tests
We study numerically the behavior of the variational discretization schemes for the Boussinesq, anelastic, and pseudo-incompressible equations for regular and irregular computational meshes. On various test cases we investigate the models' dynamical behavior and their conservation of mass and energy. Descending from variational principles, these schemes are supposed to conserve the latter quantities to a high degree. Moreover as these models support internal gravity waves, we are interested in a quantitative evaluation of their discrete dispersion relations.
For all ensuing test cases we perform the simulations on regular and irregular triangular meshes. The regular meshes consist of equilateral triangles of constant edge length f = |f ij |, where f ij , j = 1, 2, 3, denote the edges of triangle T i (cf. notation of Section 5). The distance between neighboring vertices in x-direction is given by f x := f while the height of the triangles in z-direction is given by f z := √ 3 2 f . Given a domain size of L x × L z , in which L x and L z denote the domain's length in x-and z-directions, respectively, the mesh resolution, denoted by 2 · N x × N z for N x := L x /f x and N y := L z /f z , corresponds to the number of triangular cells.
To construct irregular meshes, we start from a regular one and randomly move the regularly distributed internal vertices -i.e. vertices that do not belong to boundary cells -from point x i = (x i , z i ) to x i + δx i within the bounds |δx i | < c · f x · r and |δz i | < c · f z · r, for a positive constant c and some random number r ∈ [−0.5, 0.5]. Leaving the boundary triangles regular is not necessary in terms of the discretization schemes, however it eases the way of how we implemented the boundary conditions (a more general treatment is possible, but we leave it to future work).  shown in [2], we apply for the advection dominated test cases in Section 6.2 the irregular mesh (i), where we restrict the irregularity of the grid to max x∈Ω ∆h(x) ≤ 1.5 by setting c = 0.1. In contrast to such rather smooth grid, we apply for the hydrostatic adjustment test cases in Section 6.1 the irregular mesh (ii) for c = 0.2 with max x∈Ω ∆h(x) ≤ 7 (cf. Figure 6.1) in order to illustrate that our variational schemes are also stable on strongly deformed meshes while conserving mass and total energy to a high degree.

Internal gravity waves by hydrostatic adjustment
An essential ingredient of the derivation of Boussinesq, anelastic and pseudo-incompressible models is the assumption of a vertically varying reference state that is in hydrostatic balance, i.e. the gravitational and pressure terms compensate each other (cf. Section 2). When out of equilibrium, the system tends to a balanced state by the so-called hydrostatic adjustment process [10] by emitting internal gravity waves.
Here, we investigate the above derived variational integrators using the hydrostatic adjustment process. We study their dynamical behavior, long term energy and mass conservation properties, and their discrete dispersion relations. To this end, we initialize the Boussinesq scheme as in [4], and adapt the therein suggested test case to suit also for the anelastic and pseudo-incompressible schemes. This will allow us to compare quantitatively the simulation results of our schemes with each other and with those of [4].
We use for the simulations in this section a computation domain of dimension (x, z) ∈ D = [0, L x ] × [0, L z ], L x = 24 m, L z = 1 m, while imposing periodic boundary conditions in x-direction and free-slip boundary conditions at the upper and lower boundaries of the domain. We integrate for a time interval of 100 s (in correspondence to [4]) with a time step size of ∆t = 0.25 s. Both regular and irregular computational meshes have a resolution of 2 · 384 × 20 triangular cells (cf. Figure 6.1).

Hydrostatic adjustment of the Boussinesq model
Consider the Boussinesq system in hydrostatic equilibrium between a stable reference buoyancȳ b(z) and the pressure P b of Section 2, i.e., When out of equilibrium, the system tends to a hydrostatic balance by emitting internal gravity waves that obey the dispersion relation dz , assumed to be a constant, denotes the Brunt-Väsälä frequency for the case of Boussinesq equations. One observes that the frequencies of the internal gravity waves are bounded from above by N b .
We further note that mass conservation is given by Being implicitly related to the density, we refer to this quantity, and the upcoming similar ones for anelastic and pseudo-incompressible equations, generally as mass M (t). Hence, we study whether the discrete integral of b over the domain D, denoted as mass M (t), remains constant in time and whether the discrete total energy, analogously to (2.5), is conserved.
Initialization. Analogously to [4], we initialize the system on the basis of a hydrostatic equilibrium, given by u eq (x, z) = w eq (x, z) = 0 and b eq (x, z) = −N 2 b z =:b(z), on which at t = 0 a localized positive buoyancy disturbanceb(x, z) with compact support is superimposed. Hence, the initial buoyancy field b(x, z, 0) =b(z) +b(x, z) with Brunt-Väsälä frequency N b = 1/s is given by the function with parameters: Note that [z] = m, hence the choice of N b = 1/s suggests further to set g = 1 m/s 2 and θ 0 = 1 K. Given these analytical functions, the discrete function B = {B i | for all triangles T i } is obtained by setting B i (0) = b(x i , z i , 0) for all triangles T i with cell centers at position (x i , z i ).
Results. To study numerically the dispersion relation of our discrete model, we determine the Fourier transforms of a time series of the buoyancy field b(x, z, t) for the time interval  [4]). The resulting spectra are presented in Figure 6.3. For all selected sample points, these spectra show an anisotropy manifested by the fact that they are bound from above by N b = 1/s. For the central point, the spectrum is pronounced for values of N b , in well agreement with equations (6.2): for waves with frequency near N b , the group velocity tends to zero leaving the corresponding waves trapped in the center of the domain. A very similar distribution of frequency spectra within the domain D has been found by [4]. The simulations on irregular meshes give very similar frequency spectra. For these early times, before waves that are reflected by the boundaries reach the center, one clearly observes the internal gravity waves, caused by the buoyancy perturbation, that propagate from the center along the channel in x-direction. Besides of small irregularities of the solutions on irregular meshes, in particular visible at the velocity field that is not completely symmetric with respect to the axis x = 10 m, the results obtained either on regular or irregular meshes are very similar. Figure 6.5 illustrates for regular (left column) and irregular (right column) meshes the time evolution of the relative errors (determined as ratio of current values at t over initial value at t = 0) of total energy E tot (t) (upper panels) and mass M (t) of equations (6.3) (lower panels). The middle panels show the energy transfer between kinetic E k (t) and potential E p (t) energy such that the sum of both energy sources, E tot (t), is conserved. One realizes that the initial surplus of potential energy, caused by the superimposed disturbance on the equilibrium state, triggers an oscillatory transfer between both energy sources.
As one can further deduce from these plots, the total energy shows an oscillatory behavior while being very well conserved in the mean for long integration times. Being symplectic, our variational integrator hence fulfills the expected behavior. The magnitude of oscillation is in the order of 10 −6 , but depends on the time step size; here we used ∆t = 0.25 s. Reducing the time step size by a factor of 10 decreases simultaneously the magnitude of the relative errors in  Figure 6.5: Boussinesq scheme: relative errors of total energy E tot (t) and mass M (t), and energy transfer between kinetic E k (t) and potential E p (t) energy, for regular (left block) and irregular (right block) meshes.
total, kinetic, and potential energy by the same factor (not shown). Hence, our scheme shows the expected 1st-order convergence rate with time (cf. time scheme derivation in [14]). The mass M (t) is conserved with an order of 10 −15 . The quantities of interest for the simulations on irregular meshes are equally well conserved as on regular meshes.

Hydrostatic adjustment of the anelastic model
We consider the anelastic equations in hydrostatic equilibrium. By assuming that the reference statesρ(z) andθ(z) are such that Constant values for N and σ a are obtained by takingθ(z) = αe N 2 g z andρ(z) = βe Kz , in which case σ a = K 2 4 . Analogously to equations (6.2), also this dispersion relation shows an anisotropic spectra that is bound from above by N . To see this, consider the extremes of (6.6): the lower bound at min(ω) = 0 results from k x = 0 while the upper bound max(ω) = N results for k x k y , σ. For the anelastic equations, mass conservation is given by as Dρ ∂ t θdx = − D dθ ·ρudx = − D div(ρuθ)dx + D div(ρu)θdx = 0, which follows by the anelastic constraint div(ρu) = 0 and by the choice of boundary conditions, i.e. u · n = 0, on ∂D. We study if our anelastic discretization scheme conserves the latter quantity and if the total discrete energy, a discrete version of E =ρ 1 2 |u| 2 + c pπ θ is well conserved, too.
Initialization. We aim for an initialization that produces results comparable to the Boussinesq scheme and that meets the requirements of constant N and σ a in equations (6.5). To this end, the hydrostatic equilibrium is set up by u eq (x, z) = w eq (x, z) = 0 and a reference statē θ(z) = e z+c for a constant c, on which at t = 0 a negative potential temperature perturbatioñ θ(x, z) is superimposed. The initial potential temperature field θ(x, z, 0) =θ(z) +θ(x, z) is hence given by with parameters r 0 = 0.2 · L z and β a = 0.2 · L z . To obtain a potential temperature field with comparable magnitude (in the order of θ 0 = 1 K) to the buoyancy field -but now with increasing value with height z as it would be in the atmosphere -we set c = −L z = −1 (corresponding to α = 1/e from above) giving 0.4 K at the bottom and 1 K at the top of the domain. Note that β a is slightly smaller than β b of equations (6.4); it is chosen such that at the domain center x c , z c , there is θ(x c , z c , 0) ≈ −b(x c , z c , 0) which results in similar magnitudes of oscillation. With this choice of θ, the Brunt-Väisälä frequency is N 2 = ḡ θ dθ dz = 1/s 2 , where we set g = 1 m/s 2 . This setup will allow us to compare the results obtained by the anelastic with those of the Boussinesq scheme from above and with those in [4].
The requirement that σ a is constant, restricts our choice of the stratified density fieldρ to be either a constant or an exponential function. We study our scheme using two profiles: (i) ρ(z) = e −z , which is very similar to a realistic stratification, and (ii)ρ(z) = e −8z providing a very strongly stratified density field. To remain concise, we only present the frequency spectra for both profiles, while we restrict our representation elsewhere to the realistic profile and only comment on results obtained when using the strong density stratification. The relation betweenπ andθ, see (2.2), allows us to initialize the Exner pressure by the potential temperature field viaπ (z) = g c p e −(z+c) . (6.9) for any values of specific heat at constant pressure c p (here we set c p = 1), as it will cancel out in the anelastic equations. Given these functions, the discrete ones are obtained by setting Θ i (0) = θ(x i , z i , 0),Π i = π(x i , z i ) andρ i =ρ(x i , z i ) for all triangles T i with cell centers at position (x i , z i ).
Results. Similarly to the Boussinesq scheme, we determine the frequency spectra by taking Fourier transforms of a time series of the potential temperature field θ(x, z, t) for the time interval t = [0, 100 s] at locations in D similar to those in Section 6.1.1. Figure 6.6 shows the various frequency spectra for initialization (i)ρ(z) = e −z ; left block for the regular, and right block for the irregular meshes. These spectra consist of frequencies lying between zero and max(ω) = N = 1/s, which agrees well with the frequency pattern in Figure 6.3 and reflects very well the properties of the analytical dispersion relation (6.6).
Using strong stratification (density initialization (ii)ρ(z) = e −8z ), the spectra shown in Figure 6.7 for both mesh types are similar to Figure 6.6, in particular their lower and upper boundaries. The distribution of occupied wave numbers though slightly differs from initialization (i): the upper panels show in general higher magnitude, the central panels show increased magnitude of small wave numbers, while the peaks at N = 1/s are not as distinct as for the weaker density stratification. In case of strong stratification, the magnitudes of oscillation are strongest at the upper boundary, i.e. in the region with lowest density. This causes small scale waves (noise) that are distributed over the domain and that occupy the small waves numbers in the spectra.
In Figure 6.8, we show snapshots of the potential temperature θ(x, z, t) at time t = 5 s and t = 8 s for the region [11 m, 13 m] × [0, 1 m] using the realistic density profile (i). Comparing with Figure 6.4, the velocity and potential temperature fields for both time instances and both mesh types are rather similar to those obtained with the Boussinesq scheme, noticing that the magnitude of displacement of θ from equilibrium is more enhanced in the anelastic case. Again, the irregular meshes (right column) trigger solutions that are slightly non axis-symmetric with respect to x = 10 m, but agree in general very well with the internal gravity wave propagations obtained on regular meshes. Figure 6.9 shows the relative errors of total energy (upper panels) and mass of equations (6.7) (lower panels) for regular (left column) and irregular (right column) meshes. Again, the central panels show the energy transfer between kinetic and potential energy. In all cases, these quantities are as well conserved as for the Boussinesq scheme, in particular their mean is very well conserved for long integration times. The conservation of mass shows a slight growth (in particular for the irregular mesh), but within the order of 10 −13 on a very acceptable level. For these plots we used a time step size of ∆t = 0.25 s. Similar to above, the anelastic scheme shows the expected 1st-order convergence rate with time (cf. Section 6.1.1 for more details). Remark 6.1 By settingρ = ρ 0 andθ = θ 0 as constants, the analytical Boussinesq and anelastic models agree. By initializing the discrete schemes in equations (5.5) and (5.9) accordingly, i.e. by setting θ = −b,ρ = 1,θ = θ 0 = 1, and c p = 1, we could exactly recover the numerical results of the Boussinesq by the anelastic scheme (not shown).

Hydrostatic adjustment of the pseudo-incompressible model
We consider the pseudo-incompressible equations in hydrostatic equilibrium. By assuming that the reference statesρ(z) andθ(z) are such that Therefore, the dispersion relation (6.11) is anisotropic and bound from above by max(ω) = N just as for the Boussinesq and anelastic models, which results from the same reasoning as done for equation (6.6).
For the pseudo-incompressible equations, we note that the mass conservation is given by d dt Dρ (z)θ(z)θ(x, z, t)dx = 0, (6.12) which follows by a similar argumentation as presented for equation (6.7).
Initialization. The pseudo-incompressible scheme is initialized exactly as the anelastic scheme with the only difference that here we do not need to define the Exner pressureπ(z) (cf. equation (5.13)). In particular, by choosing exponential functions for the potential temperature and the density fields, that isθ(z) = e z and (i)ρ(z) = e −z or (ii)ρ(z) = e −8z , we guarantee that the conditions put on N and σ pi to be constant, required to be able to calculate the dispersion relation (6.11), are fulfilled. As above, we study the dynamics, the conservation of the quantities of interest, and the dispersion relations for the two different density stratifications (i) and (ii).
Results. Analogously to Section 6.1.2, we determine the frequency spectra by taking Fourier transforms of a time series of θ(x, z, t) for the time interval t = [0, 100 s] at various locations in D. Figure 6.10 shows the resulting frequency spectra for initialization (i)ρ(z) = e −z ; left block for the regular, and right block for the irregular mesh. These spectra show for both regular and irregular meshes a very similar pattern as in Figure 6.6, in particular the maximal frequencies at each panel drop strongly at max(ω) = 1/s as theoretically expected. When using initialization (ii), the spectra for both regular and irregular meshes are similar to those presented in Figure 6.7 for the anelastic case; one observes also here an increase of small wave numbers occupying the frequency spectra, in particular on the upper panels (not shown). In Figure 6 both time instances and both mesh types very similar dynamical behavior as for Boussinesq and anelastic cases. Figure 6.12 shows the relative errors of total energy and mass of equation (6.12) and the energy transfer between kinetic and potential energy. On both mesh types, the errors are in the same order as for the Boussinesq and anelastic schemes. Also here, the time step size is ∆t = 0.25 s and the scheme shows the expected 1st-order convergence rate with time. The time evolution of mass shows a slight decline, but only in the order of 10 −14 , hence it is very well conserved.

Rising and falling bubble test cases
On advection dominated test cases, namely (i) a rising hot bubble in a cold environment and (ii) a cold drop of air in a warm environment, we illustrate the dynamics of the anelastic and pseudoincompressible schemes. We point out the similarities and differences between the underlying approximations and we test whether our variational integrators share these properties.

Rising bubble test cases
Using the test case of a rising hot bubble of air in a cold environment, we will verify in the following that in case of small potential temperature variations, the anelastic and pseudoincompressible schemes behave very similarly. In order to enable quantitative comparisons also with literature, we set up the test case similarly to [9]. The potential temperature distribution consisting of a constant background stateθ = 300 K and some initial disturbance θ (x, z, 0) is given by θ(x, z, 0) = 300 K + δθ cos 2 π 2 r for r ≤ 1, r = 5 with values for δθ = 2 K, L = 10 km, and z off = 0.03 L. This results in an Exner pressurē π(z) = − g cpθ z withπ(z = 0) = 0 according to (2.2). The hydrostatic background density stratification shall be described by (cf. [13])  [9]). Using this analytical functions we initialize the discrete ones as described above.
Using this initialization, we expect to obtain results very similar to those of Klein [9]. As the integration time with 1000 s is such that the bubble does not interact with the boundaries and as the initialization is axis-symmetric to x c = 10 km (cf. Figure 6.13), the difference in the choice of boundary conditions for our setup and that in [9] does not influence the results.
Results. Figure 6.14 shows the potential temperature θ(x, z, t) for the anelastic (left column) and pseudo-incompressible schemes (right column) for regular (upper row) and irregular (lower row) meshes after an integration of 1000 s. The initial data are symmetric with respect to vertical axis x c = 10 km. As the initial maximum of potential temperature is located on this symmetry line and as the boundary conditions are mirror-symmetric too, the maximum value should stay on the symmetric line and it's value should be conserved during the integration.
Choosing the initialization of Klein [9], we are able to directly compare the results. First we note that the results for the anelastic and pseudo-incompressible schemes under small variations in potential temperature (i.e. θ max /θ min = 302 K/300 K) are almost indistinguishable, in agreement with [9]. This is true for both regular and irregular meshes. Moreover, both schemes present well the advection speed of the bubbles reaching a height of 8 km after an integration of 1000 s, in well agreement with literature.
As to the symmetry of the solutions, anelastic and pseudo-incompressible schemes on regular meshes show nicely the required symmetry behavior while the maxima stay indeed on the symmetry line, but increases to θ max = 2.5 K. This is in contrast to [9], in which the solution develops two local extrema located symmetrically some finite distance away from the symmetry line while the maximal value of θ reaches, in their best case simulations, θ max = 1.73K. These latter results had been obtained by using a scheme that applies flux limiters that avoid the development of too steep gradients of θ. As our schemes do not apply such flux limiters, or any kind of smoothing methods, such steep gradients develop, in particular for latter times when the bubble is strongly stretched, leading to the increased maximal potential temperature values.
Using irregular meshes for the variational anelasic and pseudo-incompressible schemes, the general dynamical behavior of the rising bubble and, in particular, the symmetry of the solutions are well represented. Only the shapes of the stretched bubbles show some wiggles that result from the distorted triangles on the irregular meshes. Here, the maximal values of potential temperature is about θ max = 2.8 K.
From Figure 6.15 we deduce that for regular and irregular meshes both anelasic and pseudo-incompressible schemes conserve very well the total energy in the order of 10 −6 and that, for each scheme, there is almost no difference of the results obtained on either regular or irregular meshes. For all cases, the mass as defined in equations (6.7) and (6.12) is conserved in the order of 10 −14 (not shown).

Cold bubble test case
Finally we study the different behavior of the variational anelastic and pseudo-incompressible schemes in case of large potential temperature variations. Again, we set up the test case similarly to [9] which allow us to compare the results.
Initialization. The computational domain is given by ( , L x = 20 m, L z = 10 m, with periodic boundary condition in x-direction and free-slip boundary conditions in z-direction. Also here, we calculate on regular and irregular meshes with a resolution of 2 · 321 × 186 triangles. We use a time step size of ∆t = 0.0125 s. The potential temperature distribution consisting of a constant background stateθ = 300 K and some initial disturbance θ (x, z, 0) is given by θ(x, z, 0) = 300 K − δθ 2 (1+cos(πr 2 )) for r ≤ 1, r = 5 with values for δθ = 30 K and 210 K, L = 10 m, and z off = −0.02 L. The Exner pressurē π(z) required for the anelastic, and the background density stratificationρ(z) required for both schemes are the same as in the warm bubble test case. Using this analytical functions we initialize the discrete ones as described above.
Analogously to the warm bubble case, by using this initialization we expect to obtain very similar results to those presented in Klein [9], as also for this test case the bubble does not interact with the boundaries and as the initialization is axis-symmetric to x c = 10 m (cf. Figure 6.13).
Results. Figure 6.16 shows the potential temperature fields of the falling cold air bubble with δθ = 30 K (upper row) and δθ = 210 K (lower row) for the anelastic (left column) and pseudoincompressible schemes (right column) on regular meshes. We realize that for small potential density variations, the anelastic and pseudo-incompressible schemes give very similar results (upper row), as already verified above for the warm bubble test case. For larger density variations though, the corresponding potential temperature fields obviously differ (cf. lower row). The results obtained on irregular meshes are almost indistinguishable from those presented in Figure 6.16 (hence not shown).
The differences between the anelastic and pseudo-incompressible models are caused by the different pressure gradient terms, i.e by the, in time, linear term ∇(θπ) for the anelastic and by the nonlinear term θ∇π for the pseudo-incompressible models, as explained in the following. The shape of the bubble in the anelastic case (left column of Figure 6.16) is mainly cause by the horizontal variation of buoyancy force, with maximum value at the symmetry axis x c = 10 m. There, the air is the heaviest and will fall the fastest, with decreasing values to either the right or left hand side. This gradient in fall velocity will cause negative vorticity on the right and positive vorticity on the left of x c , stretching the initial round bubble into a thin sheet. This, so  called, first baroclinic effect, which merely changes even in case of high potential temperature variations, causes the stretching of the bubble when it falls down and is the reason for the stretching of the bubble in the anelastic case. The first baroclinic effect is prominent in both the anelastic and the pseudo-incompressible models. In the latter model, however, a second baroclinic effect counteracts the first one, in particular for large potential temperature variations (see lower right panel of Figure 6.16). Here, the nonlinear pressure gradient term causes a vorticity production proportional to ∇θ × ∇π. This effect is not captured by the anelastic (or Boussinesq) model as ∇ × ∇(θπ) = 0. With a maximum pressure value at the bottom of the droplet which decreases with height, and with an outward pointing potential temperature gradient, there develop hence in the pseudoincompressible case positive vorticity on the right and negative vorticity on the left of the symmetry axis, counteracting the stretching. Therefore, in the upper right panel with a relative small potential temperature variation, the stretching is dominate, but in case of a larger variation (lower right panel) the second baroclinic effect counteracts the stretching.
In general, our variational schemes represent well the first, and in case of the pseudo-incompressible scheme also the second baroclinic effect for both mesh types. Our results agree well with those of [9] for low and for large potential temperature variations up to θ min /θ max = 90 K/300 K. As we do not apply flux limiters or any smoothing in our schemes (see also discussion in Section 6.2.1), we did not result in smooth solutions for the cold bubble test case with even higher potential temperature variations (as studied for instance by Klein [9]) when integrating the pseudo-incompressible scheme for more than 1.5 s. Here, too steep gradients in θ, which are strongly under-resolved on the front of the bubble, destroy the bubble's round shape. Moreover, this effect causes in both schemes, similarly to the warm bubble case, a decrease of the minimum θ value at the bubble's front of about 5 K for a variation of θ min /θ max = 270 K/300 K and of about 20 K for θ min /θ max = 90 K/300 K. Figure 6.17 shows, for regular meshes, the relative errors of total energy for the anelastic (left column) and pseudo-incompressible schemes for potential temperature profiles with δθ = 30 K (upper row) and δθ = 210 K (lower row). The values for total energy are well conserved while also here we realize for the pseudo-incompressible scheme a decrease in accuracy when the potential temperature gradients are too steep. For all cases, the values for mass, defined in equations (6.7) and (6.12), are very well conserved in the order of 10 −14 (not shown). The relative errors obtain on irregular meshes are in the same order (also not shown).

Conclusion
In this paper, we derived variational integrators for the anelastic and pseudo-incompressible models by applying the variational discretization framework introduced in [14]. In order to enable the use of this framework, we first described the anelastic and pseudo-incompressible approximations of the Euler equations of a perfect gas in terms of the Euler-Poincaré variational method. Applying the idea of weighted volume forms, i.e. weighted in terms of the background stratifications of density (anelastic) or of density times potential temperature (pseudoincompressible), we could identify the appropriate Lie group configuration space for the two models. The corresponding equations of motion follow by Hamilton's variational principle of stationary action.
The inclusion of the anelastic and pseudo-incompressible models into the theory of Euler-Poincaré variational method allows us to apply the variational discretization approach of Pavlov et al. [14]. In fact, we extended the latter by suitable discrete diffeomorphism groups, that incorporate the idea of weighted meshes as discrete counterparts of the weighted volume forms, in order to match the divergence-free conditions of the corresponding weighted velocity fields. Alongside, we defined appropriate weighted pairing required to derive the functional derivatives of the discrete Lagrangian that leads to the corresponding discrete equations of motion. Following this pathway, we could directly apply the flat operator of [14], which is a delicate issue otherwise. In this vein, we derived structure-preserving discretizations on irregular 2D simplicial meshes for the Boussinesq, anelastic, and pseudo-incompressible models.
We verified the structure-preserving properties of these variational integrators for both regular and irregular triangular meshes with respect to two test cases. In the hydrostatic adjustment process, we considered an atmosphere in which the vertically varying reference states are in hydrostatic balance. When out of this equilibrium, the system tends to a balanced state by emitting internal gravity waves subject to a model-dependent dispersion relation. Numerical simulations showed that each of our variational integrators capture very well the characteristics of the corresponding dispersion relation, in particular the upper and lower bounds of permitted wave numbers. To test also the nonlinear behavior, we performed advection dominated tests using a warm rising, resp., cold falling bubble in a cold, resp., warm environment. Our variational integrators perform very well when compared with reference literature, in particular they capture nicely the shape and advection speed of the bubbles. In all these cases studied, both mass and energy are conserved to a high degree confirming the structure-preserving nature of our variational integrators.
Subject of our current work is the study of variational integrators for nonlinear rotating shallow-water equations. In this vein, it is our goal to extend the variational discretization framework in order to derive structure-preserving discretization for fully two and three dimensional compressible flows.