Variational integrators for anelastic and pseudo-incompressible flows

TL;DR

This paper develops structure-preserving numerical schemes for anelastic and pseudo-incompressible flows using a variational approach based on discrete Lie groups, ensuring good energy behavior and applicability to irregular meshes.

Contribution

It introduces a novel variational discretization method for these flow systems, preserving geometric structures and circulation properties.

Findings

Schemes exhibit excellent long-term energy conservation.

Applicable to irregular meshes with preserved divergence conditions.

Preliminary tests demonstrate scheme stability and accuracy.

Abstract

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.

Tables

Table 1. Table 5.1: Parallel between the continuous and discrete forms for the three models. Note that in the Euler-Poincaré form given in the sixth row of the first column, one has to compute the variational derivatives with respect to the three different weighted pairings in order to get the three models. The last row of the first column presents the continuous equations in a form that corresponds to the discrete forms obtained by variational discretization on 2D simplicial meshes. Note that these expressions are not the standard form of the models given in ( 2.4 ), ( 2.7 ), ( 2.8 ).

Continuous diffeomorphisms	Discrete diffeomorphisms
Boussinesq:   ${\mathrm{Diff}}_{\mu }\left(M\right)$	Boussinesq:   $𝖣\left(𝕄\right)$
Anelastic:   ${\mathrm{Diff}}_{\overline{\rho }\mu }\left(M\right)$	Anelastic:   ${𝖣}_{\overline{\rho }}\left(𝕄\right)$
Pseudo-incompressible:   ${\mathrm{Diff}}_{\overline{\rho }\overline{\theta }\mu }\left(M\right)$	Pseudo-incompressible:   ${𝖣}_{\overline{\rho }\overline{\theta }}\left(𝕄\right)$
Lie algebras	Discrete Lie algebras
${𝔛}_{\mu }\left(M\right),{𝔛}_{\overline{\rho }\mu }\left(M\right),{𝔛}_{\overline{\rho }\overline{\theta }\mu }\left(M\right)$	$𝔡\left(𝕄\right),{𝔡}_{\overline{\rho }}\left(𝕄\right),{𝔡}_{\overline{\rho }\overline{\theta }}\left(𝕄\right)$
Euler-Poincaré form	Discrete Euler-Poincaré form
${\partial }_t\frac{\delta \ell }{\delta 𝐮}+{\pounds }_{𝐮}\frac{\delta \ell }{\delta 𝐮}+\frac{\delta \ell }{\delta \theta }𝐝\theta =-𝐝p$,	Equation (4.12)
Common form for the three models	Common discrete form for the three models
	Form independent of the mesh
Expression corresponding to the	Discrete form on 2D simplicial grids
discrete form on 2D simplicial grids
Boussinesq:	Discrete Boussinesq:
${\partial }_t{𝐮}^{♭}+{𝐢}_{𝐮}{\mathrm{𝐝𝐮}}^{♭}=-z𝐝b-𝐝\tilde{p}$	Equation (5.2)
Anelastic:	Discrete Anelastic:
${\partial }_t{𝐮}^{♭}+{𝐢}_{𝐮}{\mathrm{𝐝𝐮}}^{♭}=c_p\overline{\pi }𝐝\theta -𝐝\tilde{p}$	Equation (5.6)
Pseudo-incompressible:	Discrete Pseudo-incompressible:
${\partial }_t\left(\frac{1}{\theta }{𝐮}^{♭}\right)+\frac{1}{\theta }{𝐢}_{𝐮}{\mathrm{𝐝𝐮}}^{♭}=-\frac{1}{{\theta }^2}\left(gz-\frac{1}{2}{\left|𝐮\right|}^2\right)𝐝\theta -𝐝\tilde{p}$	Equation (5.8)

Equations

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+\frac{1}{\rho }\nabla p=-g\mathbf{z},{\partial }_t\rho +\mathrm{div}(\rho \mathbf{u})=0,{\partial }_t\theta +\mathbf{u}\cdot \nabla \theta =0,$$

$$\pi =(p/p_0{)}^{R/c_p},$$

$$\frac{1}{\rho }\nabla p=c_p\theta \nabla \pi .$$

$$\theta (x,y,z,t)=\bar{\theta }(z)+{\theta }^{\prime }(x,y,z,t),\pi (x,y,z,t)=\bar{\pi }(z)+{\pi }^{\prime }(x,y,z,t),$$

$$c_p\bar{\theta }\frac{d\bar{\pi }}{dz}=-g.$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+c_p\theta \nabla {\pi }^{\prime }=g\frac{{\theta }^{\prime }}{\bar{\theta }}\mathbf{z},{\partial }_t\rho +\mathrm{div}(\rho \mathbf{u})=0,{\partial }_t\theta +\mathbf{u}\cdot \nabla \theta =0.$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+c_p{\theta }_0\nabla {\pi }^{\prime }=g\frac{{\theta }^{\prime }}{{\theta }_0}\mathbf{z},\mathrm{div}\mathbf{u}=0,{\partial }_t\theta +\mathbf{u}\cdot \nabla \theta =0,$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+\nabla P_{\mathrm{b}}^{\prime }=b^{\prime }\mathbf{z},\mathrm{div}\mathbf{u}=0,{\partial }_t{b}^{\prime }+\mathbf{u}\cdot \nabla b^{\prime }+N^{2}w=0,$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+\nabla P_{\mathrm{b}}=b\mathbf{z},\mathrm{div}\mathbf{u}=0,{\partial }_{t}b+\mathbf{u}\cdot \nabla b=0,$$

$${\partial }_{t}E+\mathrm{div}((E+P_{\mathrm{b}})\mathbf{u})=0.$$

$$\mathrm{div}(\bar{\rho }\mathbf{u})=0,$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+c_p{\theta }_0\nabla {\pi }^{\prime }=g\frac{{\theta }^{\prime }}{{\theta }_0}\mathbf{z},\mathrm{div}(\bar{\rho }\mathbf{u})=0,{\partial }_t\theta +\mathbf{u}\cdot \nabla \theta =0.$$

$${\partial }_{t}E+\mathrm{div}((E+P_{{\mathrm{a}}_0})\mathbf{u})=0$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+c_p\bar{\theta }\nabla {\pi }^{\prime }=g\frac{{\theta }^{\prime }}{\bar{\theta }}\mathbf{z}.$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+\nabla (c_p\bar{\theta }{\pi }^{\prime })=g\frac{{\theta }^{\prime }}{\bar{\theta }}\mathbf{z},\mathrm{div}(\bar{\rho }\mathbf{u})=0,{\partial }_t\theta +\mathbf{u}\cdot \nabla \theta =0.$$

$${\partial }_{t}E+\mathrm{div}((E+P_{\mathrm{a}})\mathbf{u})=0,$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}+c_p\theta \nabla {\pi }^{\prime }=g\frac{{\theta }^{\prime }}{\bar{\theta }}\mathbf{z},\mathrm{div}(\bar{\rho }\bar{\theta }\mathbf{u})=0,{\partial }_t\theta +\mathbf{u}\cdot \nabla \theta =0.$$

$${\partial }_{t}E+\mathrm{div}((E+P_{\mathrm{pi}})\mathbf{u})=0$$

$$\mathrm{div}(\bar{\sigma }\mathbf{u})=0,$$

$${\pounds }_{\mathbf{u}}\mu ={\mathrm{div}}_{\mu }(\mathbf{u})\mu ,$$

$${\mathfrak{X}}_{\mu }(\mathcal{D})=\{\mathbf{u}\in \mathfrak{X}(\mathcal{D})\mid {\mathrm{div}}_{\mu }(\mathbf{u})=0,\mathbf{u}\Vert \partial \mathcal{D}\}.$$

$${\mathrm{div}}_{\mu }(\bar{\sigma }\mathbf{u})=\bar{\sigma }{\mathrm{div}}_{\bar{\sigma }\mu }(\mathbf{u}).$$

$${\pounds }_{\mathbf{u}}\alpha =\mathbf{d}\left({\mathbf{i}}_{\mathbf{u}}\alpha \right)+{\mathbf{i}}_{\mathbf{u}}\mathbf{d}\alpha ,\mathbf{d}(\alpha \wedge \beta )=\mathbf{d}\alpha \wedge \beta +(-1{)}^k\alpha \wedge \mathbf{d}\beta ,$$

$${\mathbf{i}}_{\mathbf{u}}(\alpha \wedge \beta )={\mathbf{i}}_{\mathbf{u}}\alpha \wedge \beta +(-1{)}^k\alpha \wedge {\mathbf{i}}_{\mathbf{u}}\beta .$$

$${\mathrm{div}}_{\mu }(\bar{\sigma }\mathbf{u})\mu$$

$$=\left({\mathbf{i}}_{\mathbf{u}}\mathbf{d}\bar{\sigma }\right)\mu -{\mathbf{i}}_{\mathbf{u}}\left(\mathbf{d}\bar{\sigma }\wedge \mu \right)+\bar{\sigma }{\mathrm{div}}_{\mu }\mathbf{u}=\left(\mathbf{d}\bar{\sigma }\cdot \mathbf{u}\right)\mu +\bar{\sigma }{\mathrm{div}}_{\mu }\mathbf{u}.$$

$$\bar{\sigma }{\mathrm{div}}_{\bar{\sigma }\mu }(\mathbf{u})\mu ={\pounds }_{\mathbf{u}}(\bar{\sigma }\mu )=\left(\mathbf{d}\bar{\sigma }\cdot \mathbf{u}\right)\mu +\bar{\sigma }{\pounds }_{\mathbf{u}}\mu =\left(\mathbf{d}\bar{\sigma }\cdot \mathbf{u}\right)\mu +\bar{\sigma }{\mathrm{div}}_{\mu }(\mathbf{u}).$$

$$G={\mathrm{Diff}}_{\bar{\rho }\mu }(\mathcal{D})\text{and}G={\mathrm{Diff}}_{\bar{\rho }\bar{\theta }\mu }(\mathcal{D}),$$

$${\mathfrak{X}}_{\bar{\rho }\mu }(\mathcal{D})$$

$${\mathfrak{X}}_{\bar{\rho }\bar{\theta }\mu }(\mathcal{D})$$

$$\delta {\int }_0^T{L}_{{\Theta }_0}(\phi ,\dot{\phi })dt=0,$$

$$\mathbf{u}(t,\phi (t,X))=\dot{\phi }(t,X)\text{and}\theta (t,\phi (t,X))={\Theta }_0(X).$$

Full text

1 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 waves 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, and pseudo-incompressible approximations of the full 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 these approximations of the Euler equations 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 Lagrangians, 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 equations follow from the Euler-Poincaré variational principle, for appropriate Langrangians. In Section 4 we recall the variational discretization framework introduced by [14] and extend it to suit anelastic and pseudo-incompressible equations. The corresponding discretizations on 2D simplicial meshes are presented in Section 5, and preliminary numerical tests are performed in 6. In Section 7 we draw conclusions and provide an outlook.

2 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 pseudo-incompressible 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

[TABLE]

where $\mathbf{u}$ is the three-dimensional velocity vector, $\rho$ is the mass density, $p$ is the pressure, $g$ is the gravitational acceleration, $\mathbf{z}$ is the unit vector directed opposite to the gravitational force. The variable $\theta$ is the potential temperature, defined by $\theta=T/\pi$, in which $T$ is the temperature and $\pi$ is the Exner pressure

[TABLE]

with $R$ the gas constant for dry air, $c_{p}$ the specific heat at constant pressure, and $p_{0}$ a constant reference pressure. Using the equation of state for a perfect gas, $p=\rho RT$, we have the relation

[TABLE]

The equations (2.1) correspond to conservation of momentum, mass, and entropy, respectively.

Let us write

[TABLE]

in which $\bar{\theta}(z)$ and $\bar{\pi}(z)$ characterize a vertically varying reference state in hydrostatic balance, that is,

[TABLE]

In terms of the perturbations $\theta^{\prime}$ and $\pi^{\prime}$, the equations (2.1) can be equivalently written as

[TABLE]

We introduce now three frequently applied approximations to these equations.

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

[TABLE]

in which $\theta_{0}$ is a constant reference potential temperature. These equations can equivalently be written as

[TABLE]

where $P^{\prime}_{\rm b}=c_{p}\theta_{0}\pi^{\prime}$, $N^{2}=\frac{g}{\theta_{0}}\partial_{z}\bar{\theta}$ is the Brunt-Väisälä frequency, and $b^{\prime}=g\frac{\theta^{\prime}}{\theta_{0}}$ is the buoyancy. Making use of the full buoyancy $b=g\frac{\theta}{\theta_{0}}=g\frac{\bar{\theta}+\theta^{\prime}}{\theta_{0}}$, we can write the system as

[TABLE]

where $P_{\rm b}:=P^{\prime}_{\rm b}+\frac{g}{\theta_{0}}\int_{0}^{z}\bar{\theta}(z)dz$.

The total energy is conserved since the energy density $E=\frac{1}{2}|\mathbf{u}|^{2}-bz=\frac{1}{2}|\mathbf{u}|^{2}-g\frac{\theta}{\theta_{0}}z$ verifies the continuity equation

[TABLE]

The requirement for nondivergent flow is easily justified only for liquids, and the errors incurred approximating the true mass conservation relation by $\operatorname{div}\mathbf{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

[TABLE]

where $\bar{\rho}(z)$ is the vertically varying density of the reference state.

In the original anelastic system presented by [13], the reference state is isentropic so that $\bar{\theta}(z)=\theta_{0}$ is constant, which results in the approximation

[TABLE]

The energy density can be written as $E=\bar{\rho}\left(\frac{1}{2}|\mathbf{u}|^{2}-g\frac{\theta}{\theta_{0}}z\right)=\bar{\rho}\left(\frac{1}{2}|\mathbf{u}|^{2}+c_{p}\bar{\pi}\theta\right)$, where $\bar{\pi}(z)=-\frac{g}{c_{p}\theta_{0}}z$ verifies the hydrostatic balance (2.2) for $\bar{\theta}(z)=\theta_{0}$. The total energy is conserved since $E$ verifies the continuity equation

[TABLE]

with $P_{\rm a_{0}}:=\bar{\rho}(c_{p}\theta_{0}\pi^{\prime}+gz)$

In the subsequent work [15], the reference potential temperature $\bar{\theta}$ was allowed to vary in the vertical, leading to the momentum equation

[TABLE]

The resulting system is however not energy conservative. In order to restore energy conservation, [11] considered the approximate momentum equation

[TABLE]

In this case, the energy density $E=\bar{\rho}\left(\frac{1}{2}|\mathbf{u}|^{2}+c_{p}\bar{\pi}\theta\right)$, where $\bar{\pi}(z)$ is such that $c_{p}\frac{\partial\bar{\pi}}{\partial z}=-\frac{g}{\bar{\theta}}$, satisfies the continuity equation

[TABLE]

with $P_{\rm a}:=\bar{\rho}(c_{p}\bar{\theta}\pi^{\prime}+gz)$.

Pseudo-incompressible approximation.

To obtain this approximation developed in [5], one defines the pseudo-density $\rho^{\ast}=\bar{\rho}\bar{\theta}/\theta$ and enforces mass conservation with respect to $\rho^{\ast}$ as $\partial_{t}\rho^{\ast}+\operatorname{div}(\rho^{\ast}\mathbf{u})=0$. When combined with $\partial_{t}\theta+\mathbf{u}\cdot\nabla\theta=0$, it yields $\operatorname{div}(\bar{\rho}\bar{\theta}\mathbf{u})=0$. These last two equations can be used with the momentum equation in (2.3) to yield the energy conservative system

[TABLE]

We note that the balance of momentum is equivalently written as $\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}+c_{p}\theta\nabla\pi=-g\mathbf{z}$, where $\pi=\bar{\pi}+\pi^{\prime}$, with $c_{p}\frac{\partial\bar{\pi}}{\partial z}=-\frac{g}{\bar{\theta}}$. The energy density $E=\rho^{\ast}\left(\frac{1}{2}|\mathbf{u}|^{2}+gz\right)$ verifies the continuity equation

[TABLE]

for $P_{\rm pi}:=c_{p}\rho^{\ast}\theta\pi$.

3 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 $\mathbf{u}(t,\mathbf{x})$:

[TABLE]

for a given strictly positive function $\bar{\sigma}(\mathbf{x})>0$ on the fluid domain $\mathcal{D}$.

We assume that the fluid domain $\mathcal{D}$ is a compact, connected, orientable manifold with smooth boundary $\partial\mathcal{D}$. In our examples, $\mathcal{D}$ is a 2D domain in the vertical plane $\mathbb{R}^{2}\ni\mathbf{x}=(x,z)$ or a 3D domain in $\mathbb{R}^{3}\ni\mathbf{x}=(x,y,z)$.

We fix a volume form $\mu$ on $\mathcal{D}$, i.e., an $n$-form, $n=\operatorname{dim}\mathcal{D}$, with $\mu(\mathbf{x})\neq 0$, for all $\mathbf{x}\in\mathcal{D}$ . If $\mathcal{D}$ is a domain in $\mathbb{R}^{3}\ni(x,y,z)$, one can take $\mu=dx\wedge dy\wedge dz$ to be the standard volume of $\mathbb{R}^{3}$ restricted to $\mathcal{D}$. We shall denote by $\operatorname{div}_{\mu}(\mathbf{u})$ the divergence of $\mathbf{u}$ with respect to the volume form $\mu$. Recall that the divergence is the function $\operatorname{div}_{\mu}(\mathbf{u})$ defined by the equality

[TABLE]

in which $\pounds_{\mathbf{u}}$ denotes the Lie derivative with respect to the vector field $\mathbf{u}$, see, e.g., [1]. When $\mu$ is the standard volume, one evidently recovers the usual divergence operator $\operatorname{div}$ on vector fields.

Diffeomorphism groups.

Let us denote by $\operatorname{Diff}_{\mu}(\mathcal{D})$ the group of all smooth diffeomorphisms $\varphi:\mathcal{D}\rightarrow\mathcal{D}$ that preserve the volume form $\mu$, i.e., $\varphi^{\ast}\mu=\mu$. The group structure on $\operatorname{Diff}_{\mu}(\mathcal{D})$ is given by the composition of diffeomorphisms. The group $\operatorname{Diff}_{\mu}(\mathcal{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 $\operatorname{Diff}_{\mu}(\mathcal{D})$ is given by the space $\mathfrak{X}_{\mu}(\mathcal{D})$ of all divergence free (relative to $\mu$) vector fields on $\mathcal{D}$, parallel to the boundary $\partial\mathcal{D}$:

[TABLE]

Given the strictly positive function $\bar{\sigma}>0$ on $\mathcal{D}$, we consider the new volume form $\bar{\sigma}\mu$ with associated diffeomorphism group and Lie algebra denoted $\operatorname{Diff}_{\bar{\sigma}\mu}(\mathcal{D})$ and $\mathfrak{X}_{\bar{\sigma}\mu}(\mathcal{D})=\{\mathbf{u}\in\mathfrak{X}(\mathcal{D})\mid\operatorname{div}_{\bar{\sigma}\mu}(\mathbf{u})=0,\;\;\mathbf{u}\,\|\,\partial\mathcal{D}\}$, respectively. In the next Lemma, we rewrite the condition $\operatorname{div}_{\bar{\sigma}\mu}(\mathbf{u})=0$ by using exclusively the divergence operator $\operatorname{div}_{\mu}$ associated to the initial volume form $\mu$.

Lemma 3.1

Let $\mathcal{D}$ be a manifold endowed with a volume form $\mu$ and let $\bar{\sigma}>0$ be a strictly positive smooth function on $\mathcal{D}$. Then we have

[TABLE]

Proof:

We will use the following properties of the Lie derivative $\pounds_{\mathbf{u}}$, the exterior differential $\mathbf{d}$, and the inner product $\mathbf{i}_{\mathbf{u}}$ on differential forms (see, e.g., [1]): for a $k$-form $\alpha$, an $n$-form $\beta$, and a vector field $\mathbf{u}$, we have

[TABLE]

On the one hand, we have

[TABLE]

On the other hand, we have

[TABLE]

This proves the result. $\square$

From this Lemma, we deduce that the appropriate Lie groups associated to the anelastic and pseudo-incompressible systems are given by

[TABLE]

respectively. Indeed, from the preceding Lemma, it follows that the Lie algebras of these groups can be written as

[TABLE]

respectively. They correspond to the anelastic and pseudo-incompressible constraints on the fluid velocity. We will continue to uses the subscript $\bar{\sigma}\mu$ when referring to both Lie groups and both Lie algebras.

Euler-Poincaré variational principles.

The diffeomorphism group $\operatorname{Diff}_{\bar{\sigma}\mu}(\mathcal{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 $\varphi(t,\_\,)\in\operatorname{Diff}_{\bar{\sigma}\mu}(\mathcal{D})$: a particle located at a point $X\in\mathcal{D}$ at time $t=0$ travels to ${\bf x}=\varphi(t,X)\in\mathcal{D}$ at time $t$. Exactly as in classical mechanics, the Lagrangian of the system is defined on the tangent bundle $T\operatorname{Diff}_{\bar{\sigma}\mu}(\mathcal{D})$ of the configuration manifold. We shall denote it by $L_{\Theta_{0}}:T\operatorname{Diff}_{\bar{\sigma}\mu}(\mathcal{D})\rightarrow\mathbb{R}$. The index $\Theta_{0}$ indicates that this Lagrangian parametrically depends on the potential temperature $\Theta_{0}(X)$ that is expressed here in the Lagrangian description.

The equations of motion in the Lagrangian description follow from the Hamilton principle

[TABLE]

over a time interval $[0,T]$, for variations $\delta\varphi$ with $\delta\varphi(0)=\delta\varphi(T)=0$.

In the Eulerian description, the variables are the Eulerian velocity $\mathbf{u}(t,\mathbf{x})$ and the potential temperature $\theta(t,\mathbf{x})$. They are related to $\varphi(t,X)$ and $\Theta_{0}(X)$ as

[TABLE]

We assume that the Lagrangian $L_{\Theta_{0}}$ can be rewritten exclusively in terms of these two Eulerian variables, and we denote it by $\ell(\mathbf{u},\theta)$. This assumption means that $L_{\Theta_{0}}$ is right-invariant with respect to the action of the subgroup

[TABLE]

of all diffeomorphisms that keep $\Theta_{0}$ invariant.

By rewriting the Hamilton principle (3.1) in terms of the Eulerian variables $\mathbf{u}$ and $\theta$, we get the Euler-Poincaré variational principle

[TABLE]

where $\mathbf{v}(t,\mathbf{x})$ is an arbitrary vector field on $\mathcal{D}$ parallel to the boundary and with $\operatorname{div}_{\mu}(\bar{\sigma}\mathbf{v})=0$, (i.e., $\mathbf{v}\in\mathfrak{X}_{\bar{\sigma}\mu}(\mathcal{D})$ by Lemma 3.1), and with $\mathbf{v}(0,\mathbf{x})=\mathbf{v}(T,\mathbf{x})=0$. The bracket $[\mathbf{u},\mathbf{v}]$, locally given by $[\mathbf{u},\mathbf{v}]^{i}:=u^{j}\partial_{j}v^{i}-v^{j}\partial_{j}u^{i}$, is the Lie bracket of vector fields.

The expressions for $\delta\mathbf{u}$ and $\delta\theta$ in (3.3) follow by taking the variation of the first and second equalities in (3.2) and defining $\mathbf{v}(t,\mathbf{x})$ as $\mathbf{v}(t,\varphi(t,X))=\delta\varphi(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 $\mathfrak{X}_{\bar{\sigma}\mu}(\mathcal{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 given a vector space $V$, a space in nondegenerate duality with $V$ is a vector space $V^{\prime}$ together with a bilinear form $\left\langle\;\,,\;\right\rangle:V^{\prime}\times V\rightarrow\mathbb{R}$ such that $\left\langle\alpha,v\right\rangle=0$, for all $v\in V$, implies $\alpha=0$ and $\left\langle\alpha,v\right\rangle=0$, for all $\alpha\in V^{\prime}$, implies $v=0$.

Lemma 3.2

The space $\Omega^{1}(\mathcal{D})/\mathbf{d}\Omega^{0}(\mathcal{D})$ of one-forms modulo exact forms is in nondegenerate duality with the space $\mathfrak{X}_{\bar{\sigma}\mu}(\mathcal{D})$, the Lie algebra of $\operatorname{Diff}_{\bar{\sigma}\mu}(\mathcal{D})$. The nondegenerate duality pairing is given by

[TABLE]

where $[\alpha]$ denotes the equivalence class of $\alpha$ modulo exact forms.

Proof:

It is well-known that if $g$ is a Riemannian metric, with $\mu_{g}$ the associated volume form on $\mathcal{D}$, then

[TABLE]

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 $\mu$ and $\bar{\sigma}\mu$ 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 $\mu_{g}$. Let $f$ be the function defined by $\bar{\sigma}\mu=f\mu_{g}$. Since $\mathcal{D}$ is orientable and connected, we have either $f>0$ or $f<0$ on $\mathcal{D}$. We can rewrite the duality pairing (3.4) as

[TABLE]

By successive applications of Lemma 3.1, we have

[TABLE]

where the last equality follows since $\mathbf{u}\in\mathfrak{X}_{\bar{\sigma}\mu}(\mathcal{D})$. This proves that $\mathbf{v}=f\mathbf{u}\in\mathfrak{X}_{\mu_{g}}(\mathcal{D})$. We can thus write the duality pairing $\left\langle\;\,,\;\right\rangle_{\bar{\sigma}}$ in terms of the nondegenerate duality pairing (3.5) as $\left\langle[\alpha],\mathbf{u}\right\rangle_{\bar{\sigma}}=\left\langle[\alpha],f\mathbf{u}\right\rangle$, which proves that it is nondegenerate. $\square$

In a similar way to (3.4), we shall identify the dual to the space of functions $\mathcal{F}(\mathcal{D})$ with itself by using the nondegenerate duality pairing

[TABLE]

Given a Lagrangian $\ell:\mathfrak{X}_{\bar{\sigma}\mu}(\mathcal{D})\times\mathcal{F}(\mathcal{D})\rightarrow\mathbb{R}$, the functional derivatives of $\ell$ are defined with respect to the parings (3.4) and (3.6) and denoted

[TABLE]

Proposition 3.3

The variational principle (3.3) yields the partial differential equation

[TABLE]

where $\pounds_{\mathbf{u}}$ denotes the Lie derivative acting on one-forms, given by $\pounds_{\mathbf{u}}\alpha=\mathbf{d}(\mathbf{i}_{\mathbf{u}}\alpha)+\mathbf{i}_{\mathbf{u}}\mathbf{d}\alpha$. This equation is supplemented with the advection equation

[TABLE]

which follows from the definition of $\theta$ in (3.2).

Proof:

By definition of the functional derivatives, we have

[TABLE]

Using the expression for $\delta\mathbf{u}$ in (3.3), integrating by parts, and using the equalities $\pounds_{\mathbf{u}}\mathbf{v}=[\mathbf{u},\mathbf{v}]$ and $\mathbf{d}(\alpha\cdot\mathbf{v})\cdot\mathbf{u}=(\pounds_{\mathbf{u}}\alpha)\cdot\mathbf{v}+\alpha\cdot(\pounds_{\mathbf{u}}\mathbf{v})$, the first term reads

[TABLE]

We can write $\mathbf{d}\left(\frac{\delta\ell}{\delta\mathbf{u}}\cdot\mathbf{v}\right)\cdot\mathbf{u}=\operatorname{div}_{\bar{\sigma}\mu}\left(\frac{\delta\ell}{\delta\mathbf{u}}\cdot\mathbf{v}\,\mathbf{u}\right)-\frac{\delta\ell}{\delta\mathbf{u}}\cdot\mathbf{v}\operatorname{div}_{\bar{\sigma}\mu}(\mathbf{u})=\operatorname{div}_{\bar{\sigma}\mu}\left(\frac{\delta\ell}{\delta\mathbf{u}}\cdot\mathbf{v}\,\mathbf{u}\right)$, since $\operatorname{div}_{\bar{\sigma}\mu}(\mathbf{u})=0$. Then, by the Gauss Theorem,

[TABLE]

since $\mathbf{u}\,\|\,\partial\mathcal{D}$. Combining these results, we thus get

[TABLE]

for all $\mathbf{v}\in\mathfrak{X}_{\bar{\sigma}\mu}(\mathcal{D})$. By Lemma 3.2, it follows that the one-form $\partial_{t}\frac{\delta\ell}{\delta\mathbf{u}}+\pounds_{\mathbf{u}}\frac{\delta\ell}{\delta\mathbf{u}}+\frac{\delta\ell}{\delta\theta}\mathbf{d}\theta$ is exact, i.e., there exists a function $p$ such that this expression equals $\mathbf{d}p$. $\square$

Remark 3.4

We note that the statement of Proposition 3.3 does not need the introduction of a Riemannian metric $g$ on $\mathcal{D}$. Only a volume form $\mu$ is fixed, together with a strictly positive function $\bar{\sigma}$. 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 $\mu$ or $\bar{\sigma}\mu$ equals $\mu_{g}$). In this case, identifying one-forms and vector fields via the flat operator $\mathbf{u}\in\mathfrak{X}(\mathcal{D})\rightarrow\mathbf{u}^{\flat}=g(\mathbf{u},\_\,)\in\Omega^{1}(\mathcal{D})$, the space $\Omega^{1}(\mathcal{D})/\mathbf{d}\Omega^{0}(\mathcal{D})$ can be identified with the space of vector fields $\mathfrak{X}(\mathcal{D})$ modulo gradient (with respect to $g$) of functions. The nondegenerate duality pairing (3.4) thus reads

[TABLE]

In terms of this duality pairing, the equations (3.7) are equivalently written as

[TABLE]

where $\nabla$ acting on a vector field is the covariant derivative associated to the Riemannian metric $g$, $\nabla$ acting on a function is the gradient relative to $g$, and $\nabla\mathbf{u}^{\mathsf{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 $\mathcal{D}$ is a subset of the vertical plane $\mathbb{R}^{2}\ni(x,z)$ or of the space $\mathbb{R}^{3}\ni(x,y,z)$, and has a smooth boundary $\partial\mathcal{D}$. We fix a volume form $\mu$ on $\mathcal{D}$.

1) Anelastic equations.

For the anelastic equation, we take $\bar{\sigma}(z)=\bar{\rho}(z)$, the reference mass density. The Lagrangian is given by

[TABLE]

where $\bar{\pi}(z)$ is such that $c_{p}\frac{\partial\bar{\pi}}{\partial z}=-\frac{g}{\bar{\theta}}$ and the norm is computed relative to the standard inner product on $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$.

Relative to the pairings (3.8) and (3.6) we get

[TABLE]

so that the Euler-Poincaré equations (3.9) read $\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathsf{T}}\cdot\mathbf{u}-c_{p}\bar{\pi}\nabla\theta=-\nabla 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 $\pi^{\prime}$, we note that $\nabla\mathbf{u}^{\mathsf{T}}\cdot\mathbf{u}=\frac{1}{2}\nabla|\mathbf{u}|^{2}$ and that $-c_{p}\bar{\pi}\nabla\theta$ differs from $-g\frac{\theta^{\prime}}{\bar{\theta}}\mathbf{z}$ by a gradient term, indeed:

[TABLE]

Therefore, with $\pi^{\prime}$ defined in terms of $p$ by the equality $c_{p}\bar{\theta}\pi^{\prime}=p+\frac{1}{2}|\mathbf{u}|^{2}-g\mathbf{z}-c_{p}\bar{\pi}\theta$, the Euler-Poincaré equations yield the anelastic equations (2.7).

2) Pseudo-incompressible equations.

In this case we take $\bar{\sigma}(z):=\bar{\rho}(z)\bar{\theta}(z)$ and the Lagrangian is given by

[TABLE]

As before, the kinetic energy is computed relative to the standard inner product on $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$. Relative to the pairings (3.8) and (3.6) we get

[TABLE]

so that the Euler-Poincaré equations (3.9) read

[TABLE]

After some computations, using the relation $-\frac{g}{\bar{\theta}}=c_{p}\partial_{z}\bar{\pi}$, these equations recover the pseudo-incompressible system (2.8) with $c_{p}(\bar{\pi}+\pi^{\prime})=p+\frac{1}{\theta}(\frac{1}{2}|\mathbf{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.

Theorem 3.5

Consider a domain $\mathcal{D}$ with smooth boundary $\partial\mathcal{D}$ and volume form $\mu$. The anelastic system with reference density $\bar{\rho}$, resp., the pseudo-incompressible system with reference density $\bar{\rho}$ and reference potential temperature $\bar{\theta}$ can be derived from an Euler-Poincaré variational principle for the Lie group

[TABLE]

with Lagrangian (3.10), resp., (3.12).

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

[TABLE]

where $c_{t}=\varphi(t,c_{0})$ is a loop advected by the fluid flow $\varphi(t,\_\,)$ and $\oint_{c_{t}}\alpha$ denotes the circulation of the one-form $\alpha$ along the loop $c_{t}$. Using the equation $\partial_{t}\theta+\mathbf{d}\theta\cdot\mathbf{u}=0$, one also deduces another useful form, namely,

[TABLE]

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 $c_{p}\bar{\pi}\mathbf{d}\theta=c_{p}\mathbf{d}(\bar{\pi}\theta)-c_{p}\mathbf{d}\bar{\pi}\theta=c_{p}\mathbf{d}(\bar{\pi}\theta)+g\frac{\theta}{\bar{\theta}}\mathbf{z}$, expression (3.15) yields the equivalent forms

[TABLE]

of the circulation theorem. For the pseudo-incompressible system, using (3.13) and the equalities $\frac{1}{\theta}\left(\frac{1}{2}|\mathbf{u}|^{2}-gz\right)\mathbf{d}\theta-\theta\mathbf{d}p=\mathbf{d}\left(\frac{1}{2}|\mathbf{u}|^{2}-gz\right)-c_{p}\mathbf{d}\pi\theta$, the expression (3.16) yields

[TABLE]

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.

4 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 $\mathbb{M}$ on the fluid domain $\mathcal{D}$ with cells $C_{i}$, $i=1,...,N$, define a diagonal $N\times N$ matrix $\Omega$ consisting of cell volumes: $\Omega_{i}=\operatorname{Vol}(C_{i})$. The discretization of the group $\operatorname{Diff}_{\mu}(\mathcal{D})$ of volume preserving diffeomorphisms of $\mathcal{D}$ is the matrix group

[TABLE]

where $\operatorname{GL}(N)^{+}$ is the group of invertible $N\times N$ matrices with positive determinant, and $\mathbf{1}$ denotes the column $(1,...,1)^{\mathsf{T}}$ so that the first condition reads $\sum_{j=1}^{N}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 $\operatorname{Diff}_{\mu}(\mathcal{D})$ on the space $\mathcal{F}(\mathcal{D})$ of functions on $\mathcal{D}$, given by

[TABLE]

This linear map has two key proporties:

$(1)$ it preserves the $L^{2}$ inner product of functions;

$(2)$ it preserves constant functions $C$ on $\mathcal{D}$: $C\circ\varphi^{-1}=C$.

In the discrete setting, a function is discretized as a vector $F\in\mathbb{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

[TABLE]

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^{\mathsf{T}}\Omega q=\Omega$ and $q\cdot\mathbf{1}=\mathbf{1}$ are indeed the discrete analogues of the conditions $(1)$ and $(2)$ above, respectively.

The Lie algebra of $\mathsf{D}(\mathbb{M})$, denoted $\mathfrak{d}(\mathbb{M})$, is the space of $\Omega$-antisymmetric, row-null matrices:

[TABLE]

The component $A_{ij}$ of the matrix $A$ is the weighted flux of the vector field $\mathbf{u}$ through the face common to the cells $C_{i}$ and $C_{j}$. This relation induces a nonholonomic constraint on the Lie algebra $\mathfrak{d}(\mathbb{M})$ as only the fluxes through adjacent cells are non-zero:

[TABLE]

in which $N(i)$ is the set of all indices of cells adjacent to cell $C_{i}$.

Once a discrete Lagrangian $\ell_{d}:\mathfrak{d}(\mathbb{M})\rightarrow\mathbb{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 $\mathbb{M}$ on $\mathcal{D}$, a reference density $\bar{\rho}(z)$, and a reference potential temperature $\bar{\theta}(z)$ on $\mathcal{D}$, we define, respectively, the diagonal matrices $\Omega^{\bar{\rho}}$ and $\Omega^{\bar{\rho}\bar{\theta}}$ of $\bar{\rho}$ -weighted and $\bar{\rho}\bar{\theta}$ -weighted volumes as

[TABLE]

The discrete versions of the diffeomorphism groups $\operatorname{Diff}_{\bar{\rho}\mu}(\mathcal{D})$ and $\operatorname{Diff}_{\bar{\rho}\bar{\theta}\mu}(\mathcal{D})$ in (3.14) are therefore

[TABLE]

with Lie algebras

[TABLE]

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 $\mathbf{u}$, by taking into account the weights $\bar{\rho}$ and $\bar{\theta}$. For $\mathbf{u}\in\mathfrak{X}_{\bar{\rho}\mu}(\mathcal{D})$, resp., $\mathbf{u}\in\mathfrak{X}_{\bar{\rho}\bar{\theta}\mu}(\mathcal{D})$, we get

[TABLE]

in which $D_{ij}$ is the boundary common to $C_{i}$ and $C_{j}$, and $\mathbf{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 $\mathfrak{d}_{\bar{\sigma}}(\mathbb{M})$, namely,

[TABLE]

This constraint induces a right-invariant linear constraint on the Lie group $\mathsf{D}_{\bar{\sigma}}(\mathbb{M})$: at $q\in\mathsf{D}_{\bar{\sigma}}(\mathbb{M})$, the constraint is defined as

[TABLE]

In addition to the matrix $A\in\mathfrak{d}_{\bar{\sigma}}(\mathbb{M})$ which discretizes the Eulerian velocity $\mathbf{u}$, we introduce the discrete potential temperature $\Theta\in\mathbb{R}^{N}$ whose component $\Theta_{i}$ is the average of the potential temperature $\theta$ 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_{\Theta_{0},d}:T\mathsf{D}_{\bar{\sigma}}(\mathbb{M})\rightarrow\mathbb{R}$ defined on the tangent bundle of the Lie group $\mathsf{D}_{\bar{\sigma}}(\mathbb{M})$ and being an approximation of the Lagrangian in (3.1). The parameter $\Theta_{0}\in\mathbb{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

[TABLE]

vanishing at $t=0,T$, and with $\dot{q}\in\mathcal{S}_{\bar{\sigma}}(q)$.

In a similar way with the continuous case, the relation between the Lagrangian variables $q(t)\in\mathsf{D}_{\bar{\sigma}}(\mathbb{M})$, $\Theta_{0}\in\mathbb{R}^{N}$, and the Eulerian variables $A(t)\in\mathfrak{d}_{\bar{\sigma}}(\mathbb{M})$, $\Theta(t)\in\mathbb{R}^{N}$, is given by the formulas

[TABLE]

The discrete Lagrangian $L_{\Theta_{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

[TABLE]

The Eulerian version of the Lagrange-d’Alembert principle (4.6) is found to be

[TABLE]

in which $A\in\mathcal{S}_{\bar{\sigma}}$ and where $Y$ is an arbitrary time dependent matrix in $\mathcal{S}_{\bar{\sigma}}$ vanishing at the endpoints. The expressions for the variations $\delta A$ and $\delta\Theta$ in (4.8) are obtained by using the two relations (4.7). In particular, we have $Y=\delta qq^{-1}$, which is therefore an arbitrary time dependent matrix in $\mathcal{S}_{\bar{\sigma}}$ 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 $\mathfrak{d}_{\bar{\sigma}}(\mathbb{M})$ and the space of discrete functions $\Omega_{d}^{0}(\mathbb{M})=\mathbb{R}^{N}$. To this end, we recall from [7] that in the context of discrete diffeomorphism groups, a discrete one-form on $\mathbb{M}$ is identified with a skew-symmetric $N\times N$ matrix. The space of discrete one-forms is denoted by $\Omega^{1}_{d}(\mathbb{M})$. The discrete exterior derivative of a discrete function $F\in\Omega_{d}^{0}(\mathbb{M})$ is the discrete one-form $\mbox{d}F$ given by

[TABLE]

Then, given a strictly positive function $\bar{\sigma}(\mathbf{x})>0$, the discrete version of the $L^{2}$ pairing (3.4) is defined as

[TABLE]

By repeating the arguments of Theorem 2.4 of [7] for the discrete $L^{2}$ pairing (4.9) with weight $\bar{\sigma}$, we get the identification

[TABLE]

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

[TABLE]

A direct application of the principle (4.8) yields the following result.

Proposition 4.1

A curve $(A(t),\Theta(t))\in\mathfrak{d}_{\bar{\sigma}}(\mathbb{M})\times\mathbb{R}^{N}$ is critical for the principle (4.8) if and only if there exists a discrete function $P\in\mathbb{R}^{N}$ such that the following equation holds

[TABLE]

for all $j\in N(i)$, where the functional derivatives $\frac{\delta\ell_{d}}{\delta A}$ and $\frac{\delta\ell_{d}}{\delta\Theta}$ are computed with respect to the pairings (4.9) and (4.11).

Equation (4.12) yields a structure-preserving spatial discretization of the Euler-Poincaré equation (3.7) on the mesh $\mathbb{M}$. For the anelastic, resp., pseudo-incompressible equations, we will choose $\bar{\sigma}=\bar{\rho}$, resp., $\bar{\sigma}=\bar{\rho}\bar{\theta}$ and use in (4.12) suitable approximations

[TABLE]

of the Lagrangians (3.10), resp., (3.12).

Anelastic system.

The discrete Lagrangian associated to (3.10) is

[TABLE]

where the first, resp., the second duality pairing is given in (4.9), resp., (4.11), and $\bar{\Pi}\in\mathbb{R}^{N}$ is a discretization of the reference value $\bar{\pi}(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 $\mathcal{D}$ and is based on a suitable flat operator $A\in\mathcal{S}_{\bar{\rho}}\mapsto A^{\flat}\in\Omega_{d}^{1}(\mathbb{M})$ associated to the mesh $\mathbb{M}$, see [14].

The functional derivatives of $\ell_{d}$ with respect to the pairings $\left\langle\!\left\langle\,,\right\rangle\!\right\rangle_{\bar{\rho}}$ and $\left\langle\,,\right\rangle_{\bar{\rho}}$ are, respectively,

[TABLE]

Using them in (4.12) with $\bar{\sigma}=\bar{\rho}$, we get the structure-preserving spatial discretization of the anelastic system on the mesh $\mathbb{M}$ as

[TABLE]

Pseudo-incompressible system.

The discrete Lagrangian associated to (3.12) is

[TABLE]

in which $Z\in\mathbb{R}^{N}$ is a discretization of the height $z$. Note that we are now using the volumes $\Omega^{\bar{\rho}\bar{\theta}}_{i}$.

The functional derivatives of $\ell_{d}$ with respect to the pairings $\left\langle\!\left\langle\,,\right\rangle\!\right\rangle_{\bar{\rho}\bar{\theta}}$ and $\left\langle\,,\right\rangle_{\bar{\rho}\bar{\theta}}$ are, respectively,

[TABLE]

Using them in (4.12) with $\bar{\sigma}=\bar{\rho}\bar{\theta}$, we get the structure-preserving spatial discretization of the pseudo-incompressible system on the mesh $\mathbb{M}$ as

[TABLE]

5 Variational integrator 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):

[TABLE]

The flat operator on a 2D simplicial mesh is defined by the following two conditions, see [14],

[TABLE]

in which $e$ denotes the node common to triangles $T_{i},T_{j},T_{k}$ and $\zeta^{2}_{e}$ denotes the dual cell to $e$. In (5.1), the vorticity $\omega(K)$ of a discrete one-form $K\in\Omega^{1}_{d}(\mathbb{M})$ is defined by

[TABLE]

where the sum is taken over the dual edges in the boundary $\partial\zeta_{e}^{2}$ counterclockwise around node $e$. The constant $K_{j}^{e}$ is defined as

[TABLE]

where $|\zeta_{e}^{2}|$ and $|\zeta^{2}_{e}\cap T_{j}|$ denote, respectively, the areas of $\zeta^{2}_{e}$ and $\zeta_{e}\cap T_{j}$. Note that the matrix $A^{\flat}$ defined in (5.1) is skew-symmetric, hence $A^{\flat}\in\Omega_{d}^{1}(\mathbb{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 $\operatorname{div}_{\mu}(\mathbf{u})=0$ and that the Boussinesq Lagrangian is given by

[TABLE]

The discrete Lagrangian is therefore chosen as $\ell_{d}:\mathfrak{d}(\mathbb{M})\times\mathbb{R}^{N}\rightarrow\mathbb{R}$,

[TABLE]

where $B\in\mathbb{R}^{N}$ is the discrete buoyancy and $Z$ is the discrete height function, i.e., $Z_{i}$ denotes the $z$-coordinate of the circumcenter of cell $C_{i}$.

Using the Boussinesq Lagrangian and the flat operator (5.1), the discrete Euler-Poincaré equation (4.12) yields

[TABLE]

where $\Omega_{i}A_{ij}=-\Omega_{j}A_{ji}$, for all $i,j$, and $\sum_{j\in N(i)}A_{ij}=0$, for all $i$, and where $\tilde{P}$ is related to $P$ in (4.12) via

[TABLE]

We note that the momentum equation (5.2) corresponds to the discretization of the following form of the Boussinesq equation:

[TABLE]

where, similarly to (5.3), $\tilde{p}=\mathbf{i}_{\mathbf{u}}\mathbf{u}^{\flat}+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 $P_{\rm b}=zb+p+\frac{1}{2}|\mathbf{u}|^{2}=zb+\tilde{p}-\frac{1}{2}|\mathbf{u}|^{2}$.

Anelastic flow.

The continuous and discrete anelastic Lagrangians are given in (3.10) and (4.14). Recall that in this case $\operatorname{div}_{\mu}(\bar{\rho}\mathbf{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

[TABLE]

in which $(\cdot)^{(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

[TABLE]

where $\Omega_{i}^{\bar{\rho}}A_{ij}=-\Omega_{j}^{\bar{\rho}}A_{ji}$, for all $i,j$, and $\sum_{j\in N(i)}A_{ij}=0$, for all $i$, and where $\tilde{P}$ is related to $P$ in (4.12) and (4.15) as before via the formula (5.3). We note that the momentum equation (5.6) corresponds to the discretization of the following form of the anelastic equation:

[TABLE]

where, similarly to (5.3), $\tilde{p}=\mathbf{i}_{\mathbf{u}}\mathbf{u}^{\flat}+p$, with $p$ the pressure function arising in the Euler-Poincaré formulation (3.7). The form (5.7) was shown in §3 to be equivalent to the standard form (2.6).

Pseudo-incompressible flow.

The continuous and discrete pseudo-incompressible Lagrangians are given in (3.12) and (4.16). Recall that in this case $\operatorname{div}_{\mu}(\bar{\rho}\bar{\theta}\mathbf{u})=0$. We take the flat operator (5.1) with the first line modified as in (5.5)

The general discrete pseudo-incompressible equations (4.17) yield

[TABLE]

​​where $\Omega_{i}^{\bar{\rho}\bar{\theta}}A_{ij}=-\Omega_{j}^{\bar{\rho}\bar{\theta}}A_{ji}$, for all $i,j$, and $\sum_{\in N(i)}A_{ij}=0$, for all $i$, and where $\tilde{P}$ is related to $P$ in (4.12) and (4.17) by the fomula

[TABLE]

We note that the momentum equation (5.8) corresponds to the discretization of the following form of the pseudo-incompressible equation:

[TABLE]

where, similarly to (5.9), $\tilde{p}=\frac{1}{\theta}\mathbf{i}_{\mathbf{u}}\mathbf{u}^{\flat}+p$, with $p$ the pressure function arising in the Euler-Poincaré formulation (3.7). The form (5.10) was shown in §3 to be equivalent to the standard form (2.8).

We present in Table 5.1 a parallel between the continuous and discrete variational formulations for the three models.

Time integration.

Since the spatial discretization has been 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. This approach is based on the use of the Cayley transform, a local approximant of the exponential map. For the general discrete Euler-Poincaré system (4.12) and for a given time step $\Delta t$, it results in the following scheme

[TABLE]

where $(A_{ij}^{k-1},\Theta^{k-1}_{i})$ and $(A_{ij}^{k},\Theta_{i}^{k})$ are the values at the consecutive time steps $k-1$ and $k$.

6 Numerical tests

In this section we present preliminary numerical tests for the variational schemes. We will focus on hydrostatic adjustment processes and make for each model a quantitative evaluation of the discrete dispersion relation of the emitted internal gravity waves. The simulations are performed on a regular and an irregular triangular mesh.

Description of the meshes.

The regular mesh consists 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. 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}:=\frac{\sqrt{3}}{{2}}f$. Given a domain size of $L_{x}\times 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\cdot N_{x}\times 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 the irregular mesh, we start from the regular one and randomly move the regularly distributed internal vertices – i.e. vertices that do not belong to boundary cells – from point ${\bf x_{i}}=(x_{i},z_{i})$ to ${\bf x_{i}}+\delta{\bf x_{i}}$ within the bounds $|\delta x_{i}|<c\cdot f_{x}\cdot r$ and $|\delta z_{i}|<c\cdot f_{z}\cdot r$, for a positive constant $c$ and some random number $r\in[-0.5,0.5]$. Although not necessary, we leave the boundary triangles regular as this eases the implementation.

The distortion of the irregular mesh can be quantified using a grid quality measure introduced in [2] that measures the distortion of the dual cells: \Delta h({\bf x}):=\frac{\max_{j}\ h_{ij}}{\min_{j}\ h_{ij}}\, in which $h_{ij}$ is the length of dual edge $j$ of dual cell $\zeta^{2}_{i}$ that contains point $\bf x$. High values of $\Delta h$ indicate strongly deformed cells. For our studies we use $c=0.2$ which leads to a mesh with $\max_{{\bf x}\in\Omega}\Delta h({\bf x})\approx 7$ indicating strongly deformed dual mesh cells.

We use a computation domain of dimension $(x,z)\in\mathcal{D}=[0,L_{x}]\times[0,L_{z}],L_{x}=24\,{\rm m},L_{z}=1\,{\rm m}$, while imposing periodic boundary conditions in $x$-direction and free-slip boundary conditions at the upper and lower boundaries of the domain. Both regular and irregular computational meshes have a resolution of $2\cdot 384\times 20$ triangular cells (cf. Figure 6.1).

Description of the hydrostatic adjustment test case.

The derivations of Boussinesq, anelastic, and pseudo-incompressible models rely on 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.

Applying this test case, we study the schemes’ dynamical behavior, long term energy and mass conservation properties, and their discrete dispersion relations. 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].

Initialization.

Analogously to [4], we initialize the Boussinesq scheme on the basis of a hydrostatic equilibrium, given by $u_{\rm eq}(x,z)=w_{\rm eq}(x,z)=0$ and $b_{\rm eq}(x,z)=-N_{b}^{2}z=:\bar{b}(z)$, on which at $t=0$ a localized positive buoyancy disturbance $\tilde{b}(x,z)$ with compact support is superimposed. Hence, the initial buoyancy field $b(x,z,0)=\bar{b}(z)+\tilde{b}(x,z)$ with Brunt-Väsälä frequency $N_{b}=1/{\rm s}$ is given by the function

[TABLE]

with parameters $r_{0}=0.2\cdot L_{z}$ and $\beta_{b}=0.3\cdot L_{z}$. Note that $[z]=\text{m}$, hence the choice of $N_{b}=1/{\rm s}$ suggests further to set $g=1\,\rm{m/s^{2}}$ and $\theta_{0}=1\,{\rm K}$. Given these analytical functions, the discrete function $B=\{B_{i}|\ \text{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})$ (cf. Figure 6.2).

For the anelastic and pseudo-incompressible schemes, we aim for an initialization that produces results comparable to the Boussinesq scheme and that meets the requirements of constant $N$ and $\sigma_{a}$ or $\sigma_{\rm pi}$ (discussed later in more detail). To this end, the hydrostatic equilibrium is set up by $u_{\rm eq}(x,z)=w_{\rm eq}(x,z)=0$ and a reference state $\bar{\theta}(z)=e^{z+c}$ with constant $c$, on which at $t=0$ a negative potential temperature perturbation $\tilde{\theta}(x,z)$ is superimposed. The initial potential temperature field $\theta(x,z,0)=\bar{\theta}(z)+\tilde{\theta}(x,z)$ is hence given by

[TABLE]

with parameters $r_{0}=0.2\cdot L_{z}$ and $\beta_{a}=0.2\cdot L_{z}$. To obtain a potential temperature field with comparable magnitude (in the order of $\theta_{0}=1\,$K) to the buoyancy field, we set $c=-L_{z}$ giving $0.4\,$K at the bottom and $1\,$K at the top of the domain. The choice of $\beta_{a}$ results in an oscillation comparable in magnitude to the Boussinesq case. For this $\theta$, the Brunt-Väisälä frequency is $N^{2}=\frac{g}{\bar{\theta}}\frac{d\bar{\theta}}{dz}=1/{\rm s}^{2}$, where we set $g=1\,{\rm m/s}^{2}$. The requirement that $\sigma_{a}$ and $\sigma_{\rm pi}$ have to be constant restricts our choice of the stratified density field $\bar{\rho}$ to be either a constant or an exponential function; here we use the profile $\bar{\rho}(z)=e^{-z}$ which mimics a realistic stratification of the atmosphere.

The initialization of the anelastic scheme requires, in addition, to define a discrete Exner pressure $\bar{\Pi}$. The relation between $\bar{\pi}$ and $\bar{\theta}$, see (2.2), allows us to initialize the Exner pressure by the potential temperature field via

[TABLE]

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 $\Theta_{i}(0)=\theta(x_{i},z_{i},0)$, $\bar{\Pi}_{i}=\bar{\pi}(x_{i},z_{i})$, and $\bar{\rho}_{i}=\bar{\rho}(x_{i},z_{i})$ for all triangles $T_{i}$ with cell centers at position $(x_{i},z_{i})$ (cf. Figure 6.2).

We integrate for a time interval of $100\,$s (in correspondence to [4]) and use a fixed time step size of $\Delta t=0.25\,$s for all schemes.

Conserved quantities.

In Section 2, we considered soundproof models that provide energy conserving approximations of the Euler equations. In the following we study if the variational schemes conserve discrete versions of the associated total energies too.

In the same vein, we study if discrete versions of mass are conserved quantities in time also. We note that mass conservation in the Boussinesq case is given by

[TABLE]

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)$. For the anelastic equations, mass conservation is given by

[TABLE]

as $\int_{\mathcal{D}}\bar{\rho}\partial_{t}\theta d\mathbf{x}=-\int_{\mathcal{D}}\mathbf{d}\theta\cdot\bar{\rho}\mathbf{u}d\mathbf{x}=-\int_{\mathcal{D}}\operatorname{div}(\bar{\rho}\mathbf{u}\theta)d\mathbf{x}+\int_{\mathcal{D}}\operatorname{div}(\bar{\rho}\mathbf{u})\theta d\mathbf{x}=0$, which follows by the anelastic constraint $\operatorname{div}(\bar{\rho}\mathbf{u})=0$ and by the choice of boundary conditions, i.e. ${\bf u}\cdot{\bf n}=0$, on $\partial\mathcal{D}$. Following a similar argumentation, mass conservation for the pseudo-incompressible equations is given by

[TABLE]

Results on the dynamics.

Before discussing the quantities of interest, let us first have a look at the general dynamical behavior of the variational schemes. Figure 6.3 shows snapshots at times $t=5\,$s and $t=8\,$s of the buoyancy field $b(x,z,t)$ of the Boussinesq scheme for the central region $[11\,{\rm m},13\,{\rm m}]\times[0,1\,{\rm m}]$ of the regular (left column) and the irregular (right column) mesh. 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 the irregular mesh, in particular visible at the velocity field that is not completely symmetric with respect to the axis $x=10\,{\rm m}$, the results obtained using either the regular or the irregular mesh are very similar.

Analogously, we show in Figure 6.4 snapshots of the potential temperature $\theta(x,z,t)$ of the anelastic scheme. The snapshots for the pseudo-incompressible scheme are very similar and hence not shown. Comparing with Figure 6.3, the wave structure on the velocity and potential temperature fields are rather similar, for both time instances and both mesh types, to those obtained with the Boussinesq scheme, noticing that the magnitude of displacement of $\theta$ from equilibrium is more enhanced in the anelastic and pseudo-incompressible case. Again, the irregular mesh (right column) triggers 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 the regular mesh.

Results on the conservation properties.

Figure 6.5 illustrates 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(t)$ (upper panels) and mass $M(t)$ (lower panels) of the Boussinesq scheme for the regular (left column) and the irregular (right column) mesh. Analogously, Figure 6.6 shows these relative error values for the anelastic scheme and Figure 6.7 for the pseudo-incompressible scheme.

For all three schemes and on both mesh types, the total energy shows an oscillatory behavior while being very well conserved in the mean for long integration times. The magnitudes of these oscillations are at the order of $10^{-6}$, but they depend on the time step size; here we used $\Delta t=0.25\,{\rm s}$. Reducing the time step size by a factor of $10$ decreases simultaneously the magnitude of the relative errors in total energy by the same factor (not shown). Hence, all three variational schemes show the expected $1$st-order convergence rate with time (cf. time scheme derivation in [14]).

In case of the Boussinesq scheme, mass is conserved at the order of $10^{-14}$ for both the regular and the irregular mesh. In the anelastic and pseudo-incompressible cases, mass is conserved at the order of $10^{-13}$ for both mesh types. On the irregular mesh though we observe a slight growth in the anelastic, and a slight decline in the pseudo-incompressible case, but within the order of $10^{-13}$ on a very acceptable level.

Investigation of the frequency representation.

We study the frequency spectra of the occurring internal gravity waves for all three schemes. Consider the Boussinesq system in hydrostatic equilibrium with a reference buoyancy $\bar{b}(z)$ and a pressure $P_{b}$ balance like $\frac{\partial P_{b}}{\partial z}=\bar{b}$. When out of equilibrium, the system tends to a hydrostatic balance by emitting internal gravity waves that obey the dispersion relation

[TABLE]

with wave vector ${\bf k}=(k_{x},k_{y})\in\mathbbm{Z}\setminus 0$, in which $N_{b}^{2}:=\frac{d\bar{b}}{dz}=\frac{g}{\theta_{0}}\frac{d\bar{\theta}}{dz}$, assumed to be a constant, denotes the Brunt-Väsälä frequency for the case of Boussinesq equations.

For the anelastic equations, we assume that the reference states $\bar{\rho}(z)$ and $\bar{\theta}(z)$ are such that

[TABLE]

are constant numbers. Then, the dispersion relation takes the simple form

[TABLE]

Constant values for $N$ and $\sigma_{a}$ are obtained by taking $\bar{\theta}(z)=\alpha e^{\frac{N^{2}}{g}z}$ and $\bar{\rho}(z)=\beta e^{Kz}$, in which case $\sigma_{a}=\frac{K^{2}}{4}$.

Similarly for the pseudo-incompressible equations in hydrostatic equilibrium, we assume that the reference states $\bar{\rho}(z)$ and $\bar{\theta}(z)$ are such that

[TABLE]

are constant numbers. The dispersion relation takes the simple form

[TABLE]

For all three models, one observes that the frequency spectra of the internal gravity waves are anisotropic and bound from above by $N_{b}$, respectively $N$, in the Boussinesq, respectively anelastic or pseudo-incompressible case. To see this, consider the extremes of (6.9), for instance, but the same reasoning works for the other cases too: the lower bound at $\min(\omega)=0$ results from $k_{x}=0$ for any $k_{y}>0$ or $\sigma_{a}$, while the upper bound $\max(\omega)=N$ from $k_{x}\gg k_{y},\sigma_{a}$.

Results.

To study numerically the dispersion relations of our discrete schemes, we determine the Fourier transforms of time series of the buoyancy field $b(x,z,t)$ and of the potential temperature fields $\theta(x,z,t)$ of the anelastic and pseudo-incompressible schemes for the time interval $t\in[0,100\,{\rm s}]$ at various locations of the computation domain $\mathcal{D}$ (similar to those chosen by [4]). The resulting spectra are presented in Figure 6.8 for the Boussinesq, Figure 6.9 for the anelastic, and Figure 6.10 for the pseudo-incompressible schemes; left blocks for the regular, and right blocks for the irregular mesh.

For all selected sample points, these spectra show an anisotropy manifested by the fact that the frequencies lie between zero and $\max(\omega)=N_{b}=N=1/{\rm s}$ with a sharp drop in the spectra at the maximal frequencies $\max(\omega)$. Hence, the spectra are bound from above by $\max(\omega)$ as theoretically expected. Considering the central panel of each block, the spectrum is pronounced for values of $N_{b}$ in agreement with (6.7) or of $N$ in agreement with (6.9) or (6.11): for waves with frequency near $N_{b}$ or $N$, 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 $\mathcal{D}$ has been found by [4]. The simulations on the irregular mesh give very similar frequency spectra. Hence, for all cases the spectra reflect very well the properties of the analytical dispersion relations.

7 Conclusion

In this paper we derived variational integrators for the anelastic and pseudo-incompressible models by exploiting the variational discretization framework introduced in [14] for the discretization of incompressible fluids. 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 (pseudo-incompressible), we could identify the appropriate groups of diffeomorphisms for the two models.

Based on these results, we defined suitable discrete versions of these 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 pairings required to derive the functional derivatives of the discrete Lagrangian that leads to the corresponding discrete equations of motion for the anelastic and pseudo-incompressible models, valid on any mesh discretization of the fluid domain. We then considered in detail the case of irregular 2D simplicial meshes for these two models. For completeness, we also considered the case of the Boussinesq equations on irregular 2D simplicial meshes, thereby extending the results of [4]. For each case, we discussed the form of the equations that appears in discrete form, which is not the standard form in which these equations are usually written, see Table 5.1.

We then tested the obtained variational integrators for both regular and irregular triangular meshes by focusing on hydrostatic adjustment processes. These preliminary tests showed that our variational integrators capture very well the characteristics of the corresponding dispersion relation, in particular the upper and lower bounds of permitted wave numbers. In all cases studied, both mass and energy are conserved to a high degree, following from the structure-preserving nature of our variational integrators.

Acknowledgments

The authors thank D. Cugnet, M. Desbrun, F. Lott, and V. Zeitlin for helpful discussions and valuable feedback. Both authors were partially supported by the ANR project GEOMFLUID, ANR-14-CE23-0002-01. W. Bauer has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 657016.