Towards a geometric variational discretization of compressible fluids: the rotating shallow water equations

This paper presents a geometric variational discretization of compressible fluid dynamics. The numerical scheme is obtained by discretizing, in a structure preserving way, the Lie group formulation of fluid dynamics on diffeomorphism groups and the associated variational principles. Our framework applies to irregular mesh discretizations in 2D and 3D. It systematically extends work previously made for incompressible fluids to the compressible case. We consider in detail the numerical scheme on 2D irregular simplicial meshes and evaluate the scheme numerically for the rotating shallow water equations. In particular, we investigate whether the scheme conserves stationary solutions, represents well the nonlinear dynamics, and approximates well the frequency relations of the continuous equations, while preserving conservation laws such as mass and total energy.

1. Introduction. This paper develops a geometric variational discretization for compressible fluid dynamics. Geometric integrators form a particular class of numerical schemes which aim to preserve the intrinsic geometric structures of the equations they discretize. As a consequence, such schemes are well-known to correctly reproduce the conservation laws and the global behavior of the underlying dynamical system, [19].
One efficient way to produce geometric integrators is to exploit the variational formulation of the continuous equations and to mimic this formulation at the spatial and/or temporal discrete level. For instance, in classical mechanics, a time discretization of the Lagrangian variational formulation allows for the derivation of numerical schemes, called variational integrators, that are symplectic, exhibit good energy behavior, and inherit a discrete version of Noether's theorem which guarantees the exact preservation of momenta arising from symmetries, see [25].
An extension of this approach to the context of certain partial differential equations may be made through an appeal to their spacetime variational formulation resulting in multisymplectic schemes, [26], [21], see, e.g., [12], [13], [18] for recent developments in variational multisymplectic integrators.
The development of geometric variational integrators for the partial differential equations of incompressible fluid dynamics has been initiated in [27] for the Euler equations of a perfect fluid. This approach exploits the geometric formulation due to [1] which interprets the motion of the ideal fluid as a geodesic curve on the group of volume preserving diffeomorphisms of the fluid domain. As a consequence of this interpretation, the fluid equations arise in the Lagrangian description from the Hamilton variational principle on the group of volume preserving diffeomorphisms, for the Lagrangian given by the kinetic energy of the fluid. By making use of the relabelling symmetry of the fluid, this principle naturally induces a variational principle in the Eulerian description on the Lie algebra of this group, namely, the space of divergence free vector fields. In [27] this variational geometric formulation is implemented on a finite dimensional Lie group discretization of the group of volume preserving diffeomorphisms. A main feature of the discrete level is the occurrence of nonholonomic constraints which require the use of the Lagranged'Alembert principle, a variant of Hamilton's principle applicable to nonholonomic systems. The spatially discretized Euler equations emerge from an application of this principle on the finite dimensional Lie group approximation. This approach was extended in [16] to various equations of incompressible fluid dynamics with advected quantities, such as MHD and liquid crystals. The development of this geometric method for rotating and stratified fluids for atmospheric and oceanic dynamics was given in [14]. Improvements of this variational method in efficiency, generality, and controllability were achieved in [24]. In [6], the geometric variational discretization was extended to anelastic and pseudo-incompressible fluids on 2D irregular simplicial meshes.
In the present paper, we develop this geometric variational discretization towards the treatment of compressible fluid dynamics. This extension is based on a suitable Lie group approximation of the group of (not necessarily volume preserving) diffeomorphisms of the fluid domain, accompanied with an appropriate right invariant nonholonomic constraint obtained by requiring that Lie algebra elements are approximations of continuous vector fields. The spatial discretization of the compressible fluid equations is obtained by an application of the Lagrange-d'Alembert principle on the Lie group of discrete diffeomorphisms, for a semidiscrete Lagrangian which is assumed to have the same relabelling symmetries as the continuous Lagrangian. From these symmetries, one deduces the Eulerian version of this principle and gets the discrete equations by computing the critical curves. This geometric setting is independent of the choice of the mesh discretization of the fluid domain.
As we will see later in the paper, extending the discrete diffeomorphism approach from the incompressible to the compressible case requires nontrivial steps. Simply removing the incompressibility condition in the definition of the group of discrete volume preserving diffeomorphisms defined in [27] is not enough, since this results in a group that is too large as it contains discrete diffeomorphisms that have no physical significance. This difficulty is overcome by the introduction of a nonholonomic constraint that imposes the Lie algebra elements to correspond to a discrete vector field. Note that this extension to the compressible case was not needed for the variational discretization of anelastic and pseudo-incompressible fluids in [6]. Indeed, as shown in [6], in the continuous case these equations can be derived from the Hamilton principle written on a group of diffeomorphisms that preserve a modified volume form. Hence, as opposed to the present case, the variational discretization can still be done via a slight modification of the group of volume preserving discrete diffeomorphisms introduced in [27]. Plan of the paper. We end this introduction by giving a quick overview of the geometric variational discretization of incompressible fluids. In Sect. 2, we describe the geometric setting for the discretization of compressible fluids. In particular, for a given mesh on the fluid domain, we introduce the associated group of discrete diffeomorphisms, its Lie algebra, the nonholonomic constraint, and we identify the appropriate dual spaces and projections that allow us to apply the Lagranged'Alembert principle. In Sect. 3, we derive the discrete compressible fluid equations on 2D irregular simplicial grids. In particular, we identify the expression of the discrete Lie derivative of one-form densities on such meshes. In Sect. 4, we present a structure preserving time discretization (of Crank-Nicolson type) of the scheme and suggest a solving procedure. Moreover, we express the scheme explicitly in terms of velocity and fluid depth. In Sect. 5, we perform numerical test and evaluate if the scheme is capable of conserving stationary solutions and whether it presents well the nonlinear dynamics, frequency relations, and the conservation laws such as mass and total energy. We conclude this introduction by quickly reviewing the approach of [27] based on a discretization of the group of volume preserving diffeomorphisms. Variational discretization of incompressible fluids. Variational discretizations in mechanics always start with a proper understanding of the continuous Lagrangian description of the mechanical system, namely, the identification of its configuration manifold Q and of its Lagrangian, defined on the tangent bundle of Q, to which the Euler-Lagrange equations of motion are associated. In the case of the motion of an ideal fluid on a manifold M , following [1], the configuration space is the infinite dimensional Lie group Diff vol (M ) of volume preserving diffeomorphisms of M . The Lagrangian is defined on its tangent bundle and is given by the L 2 kinetic energy of the Lagrangian velocity of the fluid. The motion of the fluid in the material description thus formally corresponds to L 2 geodesics on Diff vol (M ). An essential property of the configuration manifold Diff vol (M ) is its group structure, which allows one to understand the Eulerian description of fluid dynamics as a "symmetry reduced" version of the material description, associated to the relabelling symmetry of the Lagrangian. One of the main features of the approach undertaken in [27] is that it allows one to preserve this symmetry at the spatially discretized level.
Let us assume that the fluid domain M is discretized as a mesh M of N cells denoted C i , i = 1, ..., N . The mesh is not assumed to be regular. We define the N × N diagonal matrix Ω with diagonal elements Ω ii = Vol(C i ), the volume of cell C i . It is shown in [27] that an appropriate discrete version of the group Diff vol (M ) 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 [27]. Consider the linear action of the group Diff vol (M ) on the space F(M ) of functions on M given by composition, i.e., This linear map has the following two properties: it preserves the L 2 inner product of functions and preserves constant functions. The discrete diffeomorphism group (1) is obtained by imposing that its linear action on discrete functions satisfies these two properties. If one chooses a discrete function to be represented by a vector F ∈ R N , whose value F i on cell i is regarded as the cell average of the continuous function on cell i, then a discrete approximation of the L 2 inner product of functions is given by With this choice, we get the conditions q · 1 = 1 and q T Ωq = Ω in (1). The spatial discretization of the incompressible Euler equations is then obtained by applying a variational principle on the discrete diffeomorphism group D vol (M) for an appropriate spatially discretized right invariant Lagrangian L = L(q,q) and with respect to appropriate nonholonomic constraints. This approach directly follows from a variational discretization of the geometric description of the Euler equations given in [1] that we briefly mentioned above. Being associated to a right invariant Lagrangian, the variational principle can be equivalently rewritten on the Lie algebra of the matrix group D vol (M), given by the space of Ω-antisymmetric, row-null N × N matrices d vol (M) = {A ∈ gl(N ) | A · 1 = 0 and A T Ω + ΩA = 0}.
As shown in [27], a matrix A ∈ d vol (M) represents the discretization of a divergence free vector field u through the identification of a matrix element with a weighted flux, i.e., where D ij is the hyperface common to cells C i and C j and n ij is the normal vector field on D ij pointing from C i to C j . This representation shows that only matrix elements associated to neighboring cells can be non-zero, which is understood as a nonholonomic constraint S ⊂ d vol (M) imposed on the Lie algebra elements and appropriately used in the variational principle. For later use, we recall that in the context of the discrete diffeomorphism group approach, a discrete zero-form (i.e., a discrete function) on M is a vector F ∈ R N . The components of such a vector are regarded as the cell averages of a continuous scalar field f ∈ C 0 (M ), i.e., . The discrete diffeomorphism group approach was developed towards applications to rotating stratified fluids in [14]. In [6] appropriate discrete diffeomorphism groups were defined to develop variational discretization of the equations of anelastic and pseudo-incompressible fluids.
2. Variational Lie group discretization of compressible fluids. In this section, we shall appropriately extend the structure preserving spatial discretization of [27] to the compressible case. As we shall see, the treatment of compressible fluids requires the inclusion of an additional nonholonomic constraint.
Group of discrete diffeomorphisms. Given a mesh M on the fluid domain M , an evident candidate for a discretization of the group of all (i.e., not necessarily volume preserving) diffeomorphisms is obtained by removing the volume preserving condition q T Ωq = Ω in (1), thereby obtaining the matrix Lie group The following lemma characterizes the dual space to d(M) relative to the duality pairing on gl(N ) given by L, A = Tr(L T ΩA).
Lemma 2.1. The dual space to d(M) with respect to the pairing (6) can be identified with the space of N × N matrices with zero diagonal: In particular, given L ∈ gl(N ), we have L, A = 0, for all A ∈ d(M) if and only if Q(L) = 0, for the projector The dual space is therefore identified with the quotient space gl(N )/d(M) • , which is clearly isomorphic to the space (7).
From the preceding Lemma, the coadjoint operator ad where [ , ] is the commutator of matrices, i.e., [A, B] = AB − BA.
As we shall see below, an element in d(M) * does not represent necessarily a discrete momentum, since A does not necessarily represent a discrete velocity.
We shall denote by {M h } h>0 a family of meshes on the fluid domain M , indexed by h = max{h Ci |C i ∈ M h }, where h Ci is the diameter of cell C i . Exactly as in the incompressible case considered in [27], given a family {M h } h>0 of meshes on M , we say that a family Consider a time dependent diffeomorphism ϕ(t) ∈ Diff(M ), fix a function f 0 on M and define the time dependent function f (t) := f 0 • ϕ(t) −1 . The time derivative of f isḟ (t) = −df (t)·u(t), the derivative of f (t) in the direction u(t). Suppose that a family q h (t) of discrete flows approximates ϕ(t) ∈ Diff(M ). From the definition of the discrete diffeomorphism group, the discrete version of Left multiplication by the matrix A h (t) thus corresponds to (minus) the discrete derivative along the discrete vector field A h (t). This motivates the following definition, see Def. 3. of [27]. Given a family {M h } h>0 of meshes on the fluid domain, we say that a family The element A ij of the matrix A(t) =q(t)q(t) −1 describes the infinitesimal exchange of fluid particles between cells C i and C j . We thus assume the same nonholonomic constraint as in the incompressible case, namely that A ij is non-zero only if cells C i and C j share a common boundary. This leads to the constraint where N (i) denotes the set of all indices (including i) of cells sharing a hyperface with cell C i . We shall now use (9) to show explicitly how the elements of the matrix A ∈ S approximate the continuous vector field u as h → 0. To do this, we shall assume some standard conditions on the family of meshes {M h } h>0 , see, e.g., [15]. Recall that the family where ρ Ci is the diameter of the largest ball that can be inscribed in C i . The family {M h } h>0 is quasi-uniform if it is shape-regular and there exists a constant γ independent of h such that We shall also consider the following condition on the family of meshes max Ci,Cj ∈M h , i∈N (j) for some α ≥ 1, where λ is independent of h, x k is the barycenter of cell C k , and x ij is the barycenter of hyperface D ij . In two dimensions, (13) is equivalent to the following condition: Every pair of adjacent triangles C i ∪ C j forms an O(h 1+α ) approximate parallelogram. That is, the lengths of any two opposite edges of C i ∪ C j differ by O(h 1+α ). This assumption is sometimes used in the finite element literature; see, for instance [2], [22], [10].
Lemma 2.2. Consider a family {M h } h>0 of meshes on M . Assume that this family is quasi-uniform and satisfies (13). and where N h is the number of cells in M h and n ij is the normal vector field on D ij pointing from C i to C j .
holds for every f ∈ C 2 (M ).
Proof. By using (14), (15), and Gauss' Theorem, we have, for Ωii Ci f (x) dx and where we drop the index h on F for simplicity.
By the Poincaré-Wirtinger inequality, we have (17) and where | · | W k,∞ (Ci) denotes the W k,∞ semi-norm on C i . We write e 1 as By the Poincaré-Wirtinger inequality, the first term is bounded by By the shape-regularity assumption, we have The second term is bounded by The last factor in (19) is then written as where x i is the barycenter of cell C i and x ij is the barycenter of hyperface D ij . Each of these terms is then estimated via Taylor expansion, giving the bound where the hypothesis (13) is used in the treatment of the third term in (20). When used in (19) this estimation has to be combined with the shape-regular assumption (11) and the quasi-uniform assumption (12) to (19) and hence The combination of all these estimations gives which proves the result.

From this Lemma, it follows that if
is an approximation of the vector field u, i.e., it satisfies (9), then A h − A u h → 0, as h → 0. Formulas (14)-(15) thus give the relation between the Lie algebra element A and the vector field u it approximates. In particular (A h ) ij → 0 for j / ∈ N (i) as h → 0, i.e., A h satisfies the sparsity constraint S in an approximate sense.
Note that the expressions (14)- (15) are consistent with the condition A · 1 = 0 in (5). From these expressions, we also deduce that the matrices A ∈ d(M) have to satisfy, in addition to the constraint A ∈ S imposed earlier, the constraint Ω ii A ij = −Ω jj A ji , for all j = i, i.e., A T Ω + ΩA is a diagonal matrix. We include this in an additional nonholonomic constraint given by This constraint is equivalently described by saying that A decomposes as A = A a + A d , where A a ∈ d vol (M) and A d is diagonal. For A ∈ R, the diagonal part is found from the equality A T Ω + ΩA = 2ΩA d . From the condition A · 1 = 0, we get A d ii = − j A a ij . Semidiscrete variational equations for compressible fluids. The derivation of the spatially discretized equations for compressible fluids is based on the following proposition.
Proposition 1 (Discrete momenta). The dual space to the constraint space R ⊂ d(M) can be identified with the space of discrete one-forms, i.e., The second part of the proposition easily follows from the first part.
The right action of a discrete diffeomorphism q on a discrete function F ∈ Ω 0 d (M) is given by matrix multiplication, i.e., F · q = q −1 F . The space of discrete densities Den d (M) is defined as the dual space to Ω 0 d (M) with respect to the pairing (3). The right action on a density D ∈ Den d (M) is defined by duality as D · q, F 0 := D, F · q −1 0 , for all F ∈ Ω 0 d (M). It is explicitly given by The associated infinitesimal action of a Lie algebra element A ∈ d(M) on a discrete density D ∈ Den d (M) is found by taking the derivative of the action with respect to q at the neutral element. We get Let us assume that a spatially discretized Lagrangian L D0 : T D(M) → R is defined on the tangent bundle T D(M) of the discrete diffeomorphism group. This Lagrangian depends on the initial density D 0 of the fluid in Lagrangian description, hence the notation L D0 . We assume that this Lagrangian preserves the symmetry of the continuous Lagrangian, namely, there exists a Lagrangian = (A, D) : d(M) × R N → R in Eulerian coordinates such that for any given initial density D 0 , we can write The spatially discretized equations follow from the Lagrange-d'Alembert variational principle (see, e.g., [7]) applied to the Lagrangian L D0 and with respect to the nonholonomic constraint S ∩ R ⊂ d(M). This variational principle reads whereqq −1 ∈ S ∩ R and with respect to variations δq such that δqq −1 ∈ S ∩ R and with δq(0) = δq(T ) = 0. In the Eulerian description, using (26), this variational principle can be rewritten as where A ∈ S ∩ R and for variations δA and δD such that We call the principle (28)- (29) an Euler-Poincaré-d'Alembert principle. The expression of the variation δA is obtained from the definition The expression of δD is obtained from the definition D = D 0 · q −1 and from (24). The passage from the Lagrange-d'Alembert principle (27) to the Euler-Poincaré-d'Alembert principle (28)-(29) is rigorously justified by employing the Euler-Poincaré reduction theory, see [20], extended to include the nonholonomic constraint S ∩ R.
where P : g(N ) → Ω 1 d (M) is the projection (23). This is the semidiscrete balance of momenta for compressible fluids. The semidiscrete continuity equation reads for all δA ∈ d(M) and δD ∈ R N . Application of the variational principle (28) yields Isolating the matrix B, integrating by parts, and using for all B ∈ S ∩ R. The result then follows from Proposition 1 and by noting that since B ∈ S, the equations are verified only for neighboring cells, i.e., j ∈ N (i). The semidiscrete continuity equation follows from the definitions Remark 1 (Coadjoint operator and discrete Lie derivative). We have seen that the coadjoint operator for the group D(M) with respect to the pairing (6) has the expression The discrete advection term in (30) is thus related to the coadjoint operator as follows where we recall that L (A) := 1 2 (L−L T ) denotes the skew-symmetric part of a matrix. We shall compute explicitly this advection operator on simplicial grids below and obtain a discretization of the Lie derivative operator on one-form densities.
Discrete Lagrangian for barotropic fluids. The Lagrangian of a compressible rotating barotropic fluid on a manifold M with Riemannian metric g has the general form with u the fluid Eulerian velocity, ρ the mass density, and ε(ρ) the internal energy density. The vector field R is the vector potential of the angular velocity of the Earth. The first and second term of the Lagrangian involve the flat operator associated to the Riemannian metric. The intrinsic expression of the equations for the compressible fluid on Riemannian manifolds is Figure 1. Notation and indexing conventions for the 2D simplicial mesh.
where ∇, grad, and div are, respectively, the covariant derivative, the gradient, and the divergence associated to the Riemannian metric g, Ω is the 2-form defined by 2Ω := dR , and i u Ω = Ω(u, ). The pressure is given by p = ∂ε ∂ρ ρ − ε. On M = R 3 , we recover the classical Coriolis term 2(i u Ω) = 2Ω × u (see, e.g., [5]). If, in addition to ρ, the internal energy depends on other non dynamical fields, like the bottom topography, then (33) has the more general expression A main ingredient in the definition of a discretization of the Lagrangian (32) on d(M) × Den d (M) is a consistent discretization of the Riemannian metric or, equivalently, of the flat operator, on a given mesh M. Since only the skew symmetric part of the discrete momenta counts (i.e., the discrete momenta are in so(N )), we can still use the same flat operators as in [27]. Given such a discrete flat operator , the discrete Lagrangian associated to (32) is 3. Compressible rotating fluids on simplicial grids. In this section we shall use Theorem 2.3 to derive a structure preserving variational discretization of 2 dimensional rotating compressible fluids on irregular simplicial grids. We shall then consider the special case of the rotating shallow water equations.
Simplicial grid. We consider a 2D simplicial mesh on the fluid domain, as described in Fig. 1. We adopt the following notations: fij : = length of a primal edge, triangle edge located between triangle i and triangle j; hij : = length of a dual edge that connect the circumcenters of triangle i and triangle j; Ωii : = area of a primal simplex (triangle) Ti.
The flat operator on a 2D simplicial mesh is defined by the following two conditions, see [27], where e denotes the node common to triangles T i , T j , T k and ζ e denotes the dual cell to e. The vorticity ω of a discrete one-form L ∈ Ω 1 d (M) is defined by ω(L), ζ e := hmn∈∂ζe L mn and the constant K e i is defined as where |ζ e ∩ T k | denotes the area of the intersection of ζ e and T k . Variational discretization of compressible fluids. A main ingredient in the variational discretization (30) is the expression of the discrete advection term in the momentum equation. We shall compute it separately in the next lemma, before giving the semidiscrete compressible fluid equations.
Lemma 3.1 (Discrete Lie derivative on simplicial grids). On a 2D simplicial grid the discrete advection term for A, B ∈ d(M) and D ∈ Den d (M) is given by In this formula, the indices i, j correspond to two neighboring cells T i and T j , the indices i ± and j ± refer to cells as indicated on Fig. 1. The discrete vorticity associated to B at the nodes ± is defined by In the last term of the third line of (38), D·A denotes the infinitesimal action of A ∈ d(M) on D ∈ Den d (M), see (25).
In the last line we introduce the notation L d A (DB ) to indicate that this expression is a discrete version of the Lie derivative, see Remark 2.
Proof. Using the second equation in (36) and 13 We also compute the diagonal matrix elements

GEOMETRIC VARIATIONAL DISCRETIZATION OF COMPRESSIBLE FLUIDS
The final formula (38) is obtained by using the three expressions above, the formula P(X) ij = 1 2 (X ij − X ji − X ii + X jj ), and noting several cancellations and rearrangements.
Remark 2 (Continuous and discrete Lie derivative). Given two vector fields u, v on a manifold and a density ρ, the Lie derivative L of the one-form density ρv is given by where d is the exterior derivative and £ is the Lie derivative of a one-form. Denoting by ω := dv the vorticity of v, we note the strict analogy between this formula and its discrete counterpart (38). The link between the discrete and continuous objects is established by v ∼ B, u ∼ A, ρ ∼ D. The first two terms in (38) correspond to the term ρ i u dv = i ρu ω, the third term in (38) corresponds to ρ d(v · u) and the fourth term corresponds to div(ρ u)v .
Proof. The functional derivatives of the discrete Lagrangian (35) with respect to the pairing , 0 and , 1 are The system (40) results from a series of computations. For the first term in (30), we have d dt

WERNER BAUER AND FRANÇ OIS GAY-BALMAZ
where in the last equality we use the discrete continuity equation (31). For the second term we use Lemma 3.1 with B = A + R . For the last term in (30), we compute Adding these three terms and noting some cancellations, we get (40).
Remark 3. The momentum equations (33) can equivalently be written in the space of one-forms as It is this expression of the compressible equations that appears in a discretized form in the variational discretization (40).
Remark 4. The advection term in (40) consists of two parts associated to either node + or node −: the weighted sums of the absolute vorticity ω + with those elements A ii+ and A jj+ that contribute to ω + and another weighted sum of ω − with A ii− and A jj− . This form follows naturally from the variational principle. In particular, it does not require the reconstruction of tangential velocities out of the prognostic normal ones, in contrast to standard C-grid schemes.
In standard C-grid discretizations of the shallow water equations (see, e.g., [3,31,35]), the advection term is different. Usually, it is a product of a reconstructed tangential velocity value and an averaged absolute vorticity value evaluated at the triangles' edge midpoints. However, it is a nontrivial task to develop reconstructions and suitable average methods that provide well-behaving, conservative discretizations, as the resulting schemes might be inconsistent [36].
The case of the rotating shallow water equations. For the rotating shallow water equations, the density variable is the fluid depth denoted h. The Lagrangian is where B is the bottom topography and where h + B describes the free surface elevation of the fluid. The above developments directly apply to this case by choosing the internal energy ε(h) i.e., the rotating shallow water equations. The formulation (41) becomes The discrete rotating shallow water equations are obtained from the first equation in (40) by replacing the last term with where D i is the discretization of the fluid depth h. In this case (40) becomes the discrete form of the formulation (42) of the rotating shallow water momentum equation.
Remark 5. Note that the first equation in (40) follows from (30) by using the advection equation (31). Without making use of this, we get We end this section by giving in Table 1 a summary that enlightens the correspondence between the continuous and discrete objects.
Lie algebra action on functions Lie algebra action on discrete functions Lie algebra action on densities Lie algebra action on discrete densities Eulerian velocity and density Eulerian discrete velocity and discrete density Discrete compressible Euler equations Form I: These are explicit representations of A introduced in Lemma 2.2 on the simplicial mesh M. Here, div(V ) is a standard FV divergence operator on a triangular mesh, see [3] for instance. Semi-discrete equations. In terms of V ij , the semidiscrete shallow water momentum equation (40) becomes where we defined for valuesR mn related to R mn by R ij = − 1 2Ωii f ijRij , analogously to the relation between V mn and A mn given by (44).
For later use we define the curl-operator and a discrete tangential gradient operator where D ± is the fluid depth at the node ±, see Fig. 1. Moreover, we define Later on, we will apply the f -plane approximation, i.e., for a Coriolis parameter f (defined further below) whose values f | ζ± = f are identical on each node. The application of the variational scheme to the rotating shallow water system on the sphere is given in [9]. Here, we include an explicit representation of the continuity equation (40) by means of V and D to enable comparisons to standard methods. Using the fact that where D ij = Di+Dj 2 . Hence, the continuity equation (40) can be written by means of a standard FV divergence operator div defined on triangles (cf. e.g. [3]). Time discretization. The explicit representations (45) and (48) permit one to apply standard time discretization methods, such as Runga-Kutta or Crank-Nicolson schemes while applying operator splitting methods. Here we proceed differently. Since the spatial discretization has been realized by variational principles in a structure preserving way, a temporal variational discretization can be implemented by following the discrete (in time) Euler-Poincaré-d'Alembert approach, analogously to what has been done in [16] and [14], to which we refer for a detailed treatment. Following [8], the discrete Euler-Poincaré approach is based on the introduction of a local approximant to the exponential map of the Lie group, chosen here as the The temporal scheme consists of the following two steps. First, we compute the advected quantities, here the fluid density D applying the Cayley transform τ . The update equation is then given by D t+1 = τ (∆tA t )D t for the time t and a time step size ∆t. This equation, in particular τ , can be represented as with I the identity matrix (cf. [14] for more details). The elements of the matrix A in terms of V ij for the simplicial mesh M are given in (44). Second, the following update equations for the momentum equation has to be solved: We solve this implicit nonlinear momentum equation by fixed point iteration for all edges ij. To enhance readability, we skip the corresponding subindices in the following. The solving algorithm reads: 1. Start loop over k = 0 with initial guess at t: V * k=0 = V t ; 2. Calculate updated velocity V * k+1 from the explicit equation: then set k + 1 = k; 3. Stop loop over k if ||V * k+1 − V * k || < for a small positive . In case of convergence V * k+1 → V t+1 , this algorithm solves the momentum equation (50). Remark 6. In general, a structure preserving time discretization for the equations of the Euler-Poincaré type (30) is obtained by applying a discrete analogue of the variational principle (28). This is the point of view followed in [16] and [14], to which we refer for the complete treatment. As explained in these papers, by an appropriate choice of the Cayley transform and by dropping cubic terms, it results in a Crank-Nicolson type time update for the momentum equation (30), i.e., (43) for 2D simplicial grids. While in absence of these cubic terms the resulting temporal scheme is no more variational, it was checked that this simplification does not significantly affect the behavior and properties of the solution on the tested configuration. We have used above this time update directly on the momentum equation as reformulated in (40)  The simulations are performed on regular and irregular triangular meshes, as shown in Fig. 2. We denote their resolutions by N = 2(N 1D ) 2 , the number of triangular cells, in which N 1D is the number of subintervals in x-or y-direction. This number is defined with respect to a regular mesh by N 1D := L x /f uni = 2L y / √ 3f uni for the uniform triangle edge length f uni = f ij for all edges ij.
The irregular mesh exhibit a refined region in the center of the domain. The mean edge lengths of the cells in the center are about twice as small as of those in the outer regions. The refinement procedure starts from a regular mesh and is controlled by a monitor function such that the topology of the mesh is conserved and entanglement of the mesh is avoided (r-adaptivity). For more details on how to construct such irregular meshes, we refer the reader to [4] and references therein.
These irregular meshes mimic realistic situations in operational forecasting, in which locally refined meshes are employed to better resolve local small scale features. It is important that these local refinement areas do not impact on the quality of the global large scale flows. We employ these irregular meshes hence to illustrate that the variational RSW scheme is capable to provide excellent results even on meshes with locally refined areas consisting of cells that are possibly strongly non-regular.
Choice of spatial and temporal resolution. We use test cases that are in the geostrophic regime in which the flow is dominated by the geostropic balance. In this context, the Rossby deformation radius L D (64) describes the length scale at which effects caused by rotation are as important as those by gravity. For the test cases studied, L D is at the order of 10 3 km (cf. Section 5.2). Our choice of domain size and spatial resolutions of 2 · 32 2 , 2 · 64 2 , 2 · 128 2 , and 2 · 256 2 throughout all simulations guarantees that L D is well resolved and geostrophic effects equally well represented as gravitational ones.
Despite Crank Nicolson is an implicit time scheme and is, as such theoretically unconditionally stable, in practice the condition number of the implicit system decreases with larger time steps until the iterative solver fails to converge. This imposes an upper bound on the time step also for implicit schemes. To evaluate the ability of the scheme to handle large time steps, we use the gravity Courant number (or CFL number) with gravitational constant g = 7.32 · 10 7 km days −2 and water depth H 0 [km], where c = √ gH 0 is the speed of the fastest traveling wave, ∆x min := min ij (h ij ) is the shortest dual edge length of the mesh, and ∆t is the time step size. In contrast to explicit schemes with necessarily C max ≤ 1, implicit schemes might reach some multiples of this. For the test cases studied below, our implicit time integrator achieves a maximal Courant number of about C max = 3.
Individually for each test case, we choose one fixed time step ∆t (unless indicated otherwise). It is bounded via (51) by C max = 3, the largest water depth H 0 applied, and ∆x min of the irregular mesh with highest resolution used. Besides guaranteeing that the iterative solver converges for all meshes applied and all flow regimes studied, the fixed time step per test case permits us to distinguish between error sources related to the spatial and to the temporal discretizations. Quantities of interest. For all test cases we are particularly interested in studying the time evolution of the relative errors in the following quantities of interest (QOI): mass, total energy, mass-weighted potential vorticity, and potential enstrophy. As these values are conserved quantities in time of the continuous shallow water equations, we study if the corresponding discrete values are conserved too. These quantities can be calculated as follows. The total mass m(t) follows as an integral of the fluid depth h(x, t) over the domain M and is approximated by The total energy E tot (t) is the sum of kinetic E kin (t) and potential energy E pot (t) which are given, respectively, by Defining the absolute vorticity ω a := curl u + f , the potential vorticity q := ωa h , and the relative potential vorticity q rel := curl u h , the conserved quantities of massweighted potential vorticity P V and potential enstrophy P E are given by 5.1. Well-balancedness and frequency representation. By means of two test cases, a lake at rest and a lake at rest with small disturbance in the surface elevation but trivial bottom topography, we investigate whether the variational integrator is capable to preserve stationary solutions of the shallow water equations without generating spurious oscillations and without generating or destroying mass. In particular we study if the scheme is well-balanced with respect to the lake at rest steady state. In addition, we numerically determine the frequency spectrum of occurring surface waves and compare it with the theory. 5.1.1. Lake at rest. This test case serves us to illustrate that our implementation can perfectly handle nontrivial bottom topography. Initializing with constant surface elevation h + B = const. and zero velocity u = 0 (cf. [23]), we expect the scheme to conserve this steady solution in case of flat and nontrivial bottom topography without exciting spurious modes. We expect further that the QOI are preserved too. Initialization. The setup for the lake at rest consists in a fluid in rest with zero initial velocity, V ij = 0 for all edges, and with a surface elevation that coincides with the background depth, here H 0 = 750 m, such that D i + B i = H 0 . The function for bottom topography B(x) describes an underwater island, positioned next to the domain center, that is given by for B = 100 m, 40 L x , and σ y = 3 40 L y . We obtain the discrete function for B i by sampling (55) at the centers of the triangles T i . Here and henceforth, we denote the discrete fields such as B i , D i and (q rel ) i with B(x, y), D(x, y), and q rel (x, y), respectively, when considering them as 2D fields that depend on the x-and y-directions (cf. Fig. 3, for instance).
Here and for the frequency test case below, we use a time step of ∆t = 60 s. According to (51) with ∆x min = 5.183 km of the irregular mesh with 2 · 64 2 cells and H 0 = 750 m, the Courant number is C = 0.99 which is below C max = 3.  Fig. 3, respectively. Clearly, even in case of a nontrivial bottom topography (Fig. 3,  left), the surface elevation is preserved at machine precision for both mesh types. In addition, all QOI are preserved at machine precision too (not shown). This allows us to conclude that the scheme perfectly satisfy the well-balanced property and does not generate spurious modes in case of nontrivial bottom topography.

Frequency spectrum of linearized shallow water equations.
Here we check if the occurring wave frequencies agree with the theory. We assume trivial bottom topography. Linearizing the equations around the undisturbed fluid depth H 0 and zero velocity and inserting plane wave solutions of the form h(x, y, t) = H 0 exp(i(kx + ly − ωt)), there follow the solutions for c = √ gH 0 , wave frequency ω, and wave numbers k, l in x, y-direction, respectively. From this relation we thus have either a stationary solution (ω = 0) or waves with frequencies greater than the Coriolis frequency f , i.e., ω ≥ f . The case ω = f , i.e., k = l = 0, corresponds to inertial oscillations which do not propagate. Because of the double periodic boundary conditions, these waves are not excited here. The case ω > f , corresponds to inertia-gravity (or Poincaré) waves, cf. [38]. Since we have a bounded, double periodic computational domain, the permitted wave numbers for inertia-gravity waves are k = nx2π Lx and l = ny2π Ly for n x , n y = 0, 1, 2, . . . with n x + n y > 0. This gives a minimum wavenumber in each direction. Hence, all frequencies are greater than f but there is no maximal wavenumber. Initialization of surface disturbance. To obtain an initialization that is close to the valuesh(x, t) = H 0 andū(x, t) = 0 around which we linearized the shallow water equations, we superimpose on the lake at rest a small disturbance of magnitude H = 7.5 m in the surface elevation. Hence, we apply the fluid depth

WERNER BAUER AND FRANÇ OIS GAY-BALMAZ
using the periodic extensions The center of the perturbation is positioned at x c1 = 1 2 L x , y c1 = 1 2 L y . To obtain a circular initial (negative) surface evaluation, we use only one value for sigma, i.e. σ x = σ y = 3 40 L y . Note that in terms of implementation, we initialize the fluid depth D i by sampling (57) at each cell center and we set all velocity values V ij to zero. We set all B i = 0 to apply trivial bottom topography.
Using the same time step ∆t = 60 s and the same meshes (i.e. ∆x min = 5.183 km) as for the lake at rest, the Courant number for case (i) with H 0 = 750 m is again C = 0.99 but for case (ii) with H 0 = 1267.5 m it is C = 1.29, both well below C max .
Results of the frequency spectrum study. Recall that we use a gravitational constant of g = 7.32 · 10 7 km days −2 and a double periodic domain with wave numbers k = nx2π Lx and l = ny2π Ly for n x , n y = 0, 1, 2, . . . . According to the dispersion relation in Eqn. (56) (right), we find for two sets of parameters, namely case (i) f = 5.31 days −1 , H 0 = 750 m, and case (ii) f = 6.903 days −1 , H 0 = 1267.5 m, the following frequencies ω(n x , n y ) in units of rad days −1 for some combinations of n x and n y : The parameters for f and H 0 have been chosen such that the flow remains for both case (i) and case (ii) in the quasi-geostrophic regime (i.e. Bu ≈ 1, cf. (64)).
To verify if these theoretical values are well represented by the variational shallow water scheme, we determine the frequencies occurring during the simulations. To this end, we numerically calculate at the center of the domain the Fourier transforms of a time series of the fluid depth D(x, z, t) for the time interval t ∈ [0, 10 days] with a sample frequency of 0.01 days. The resulting spectra for the two choices of parameters are shown in Fig. 4, left for case (i) and right for case (ii).
Besides some small background noise of waves occupying all possible wave numbers, we clearly distinguishes sharp peaks in the spectra exactly at the predicted wave numbers for both parameter sets. For the illustrated combinations of k, l ≤ 2, this perfect match of expected and numerically determined values can easily be seen when comparing the values from the table with those from Fig. 4, while for combinations with larger wave numbers the associated peaks might overlap (not shown). The overlap of the frequencies ω(2, 1) and ω(0, 2) reflects itself in a nearly doubled magnitude of the associated peak. Neither at case (i) nor (ii) we observe unphysical solutions at the frequencies f = 5.31 days −1 or f = 6.903 days −1 , respectively. In addition, we notice the dependency of the spectra on the parameter f when comparing the left and the right spectra. In the latter, the peaks are shifted slightly to higher wave numbers because of the greater Coriolis parameter, in agreement with (56).

Conservation of exact, steady solution of an isolated vortex.
Here we test if the variational shallow water scheme preserves the stationary solution of an isolated vortex. We perform long-term simulations up to 100 days and evaluate alongside the conservation properties of mass, total energy, mass-weighted potential vorticity, and potential enstrophy. For the long-term simulation we apply either regular or irregular computational meshes with only 2 · 64 2 triangular cells because potential instabilities usually occur earlier for coarser mesh resolutions. Comparing the numerical solutions at day 1, for instance, with the initial state allows us further to determine the solutions spatial convergence behavior since, being stationary, every deviation is due to numerical errors. Initialization. A stationary vortex solution of the rotating shallow water equations with trivial bottom topography has the velocity u(x) and fluid depth h(x) given by where V (r) and H(r) verify the gradient wind balance (cf. [34], for instance). This condition complies with the construction method for steady state solutions of the RSW equations suggested by [32]. Here, we apply relation (60) to construct a stationary solutions and consider a test case with trivial bottom topography B(x) = 0. We consider the following radial function to describe the velocity (resp. streamfunction): The function V (r) results from choosing the exponential vorticity profile suggested by [34] for α = 2 combined with the geophysically relevant scaling discussed in [32]. Its integration with respect to r gives the streamfunction Ψ(r). The corresponding radial function for the fluid depth, which follows from (60), reads where u 0 describes the maximal velocity and r 0 is a scaling constant, both to be determined further below. Topography can be easily included, as this steady solution is consistent with the construction method of [32] which allows for such modification. Moreover, our steady solution is an alternative example to that suggested by [32] of a steady isolated vortex. It provides a Cartesian analog for the famous test case 2 of [37] of a steady solution of the shallow water equations in spherical geometry which is frequently applied to measure a scheme's ability to preserve large-scale geostrophic balance (see [31], for instance).
To position the vortex in the center of the domain, we use the definitions x 1 = (x − x c1 ) and y 1 = (y − y c1 ) for x c1 = 1 2 L x , y c1 = 1 2 L y , similarly to (58) but omitting the periodic extension. We take a sufficiently large domain and a corresponding scaling parameter r 0 so that the fluid is at rest at the boundaries. We initialize the fluid depth D i by sampling (62) at each cell center. For the velocity, we have two options to map the analytical initial conditions to the mesh. We either sample (i) the velocity field u of (59) and (61, left) at each triangle edge midpoint before we project it onto the edge's normal direction to obtain V ij . Alternatively, we sample (ii) the streamfunction Ψ of (61, right) at the triangles' vertices and calculate the normal velocities as V ij = k × G ⊥ (Ψ) ij for k = (0, 0, 1), using the tangential gradient operator (46). Both options lead to very similar results, in particular when comparing the numerical solutions visually. However, a comparison in terms of L 2 and L ∞ error norms reveals that initialization (i) leads to slightly smaller error values, in particular on coarse meshes (not shown). Hence, in the following we only present results obtained using the velocities initialization (61, left). Parameter choice and flow regimes. We consider a set of dimensionless parameters to characterise the flow resulting from (61) and (62). The characteristic velocity is described by with characteristic length scale d = 4r 0 and H as maximal deviation of the surface elevation from the background depth H 0 . We consider further the Rossby number Ro, Froude number Fr, and Burger number Bu: with Rossby deformation radius L D = √ gH0 f . For this study, we want to consider fluids in geostrophic regime in which the flow is dominated by the geostrophic balance. This requires Ro 1. The geostrophic regime can further be classified in: (i) semi-geostrophic regime for Bu 1, (ii) quasi-geostrophic regime for Bu ≈ 1, and (iii) incompressible regime for Bu 1 (cf. [11,28], for instance). Because Fr describes the stratification of the fluid -with strong stratification in case of small Fr -the choice of Bu allows us to describe shallow water flows with different degree of compressibility: with (i) and (ii) for compressible and (iii) for almost divergence free flows. As suggested by [17]  We study both the isolated vortex test case and the dual vortex interaction of Sect. 5.3.1 in these three different flow regimes. In particular, we apply the same characteristic values as suggested by [17], which allows us to compare qualitatively as well as quantitatively the different numerical solutions. Hence, we assume for both test cases the same characteristic length of d = 4r 0 = 4σ with σ = 1 2 (σ x + σ y ) for σ x = 3 40 L x and σ y = 3 40 L y (cf. (69)). As maximum velocity we use the characteristic velocity, i.e. we use u 0 = U = 2 gH f d . Given this parameter choice, the isolated vortex is stable, cf. [34].
To assure the convergence of the iterative solver (i.e. C ≤ 3), we use a time step of ∆t = 48 s for the long-term simulations, which yields the Courant number C = 2.90 for ∆x min = 5.183 km of the irregular 2 · 64 2 mesh and H 0 = 10 km, the largest water depth studied. Because we use irregular meshes with resolutions up to 2 · 256 2 cells with ∆x min = 1.313 km, the time step for the convergence study with water depth up to H 0 = 10 km is only ∆t = 12 s, which corresponds to a Courant number of C = 2.86. Results of the long-term simulations. Being a stationary solution of the rotating shallow water equations, we expect the variational integrator to exactly preserve the initial distributions of fluid depth D and relative potential vorticity q rel of the isolated vortex even for long integration times. Here and consistently throughout the manuscript, we illustrate the quantity q rel rather than, e.g., the absolute vorticity ω a or the conserved potential vorticity q. This is because q rel highlights the positive and negative regions of the vorticity distribution and it allows us to compare further below our results with those obtained by [17].
As it can be inferred from Fig. 5 and 6 in which we compare solutions after 100 days of integration on a regular (middle) and an irregular (right) mesh with the initial conditions (left), our variational scheme performs very well because it preserves in fact both fields for both mesh types very well without generation spurious modes. In particular for the regular mesh, the position, extend, and magnitude of the vortex at initial and end states very much agree while we realize for the irregular case a slight oval shape of the initially round vortex in both D and q rel . This deviation is due to numerical errors that are caused by strongly deformed mesh cells, which is particularly apparent on coarse mesh resolution as here for a grid with only 2 · 64 2 cells. However, as it can be inferred from the convergence study below (cf. Fig. 8), this deviation reduces with, at least, 1 st -order with increasing mesh resolution.
Despite the isolated vortex is a steady state solution of the RWS equations, in which any quantity is conserved because of no time dependence, for the numerical scheme such solutions are stationary only up to numerical errors. As such, it is interesting to discuss also here the conservation properties of the QOI. In Fig. 7 we show the time evolutions of the relative errors (determined as ratio of current values at time t over initial value at t = 0) of the total energy E, determined on a mesh with 2 · 64 2 cells (upper row) and with 2 · 32 2 cells (lower row), for fluids in semi-geostrophic (first and second column), in quasi-geostrophic (third and forth column), and in incompressible (fifth and sixth column) regimes. As above, we compare results for regular (first, third, fifth column) with irregular (second, forth, sixth column) meshes. Here, and for all other cases studied, mass m is preserved up to machine precision (not shown).
Comparing the three different flow regimes, we note that the errors in total energy depend on the fluid depth, with smallest values of 10 −10 in the incompressible case. In the other regimes that allow compressibility, total energy is also very Figure 5. Isolated vortex test case: fluid depth D(x, y) at initial time t = 0 (left) and at t = 100 days on a regular (center) and an irregular (right) mesh with 2 · 64 2 triangular cells. Contours between 0.698 km and 0.752 km with interval of 0.003 km. Figure 6. Isolated vortex test case: relative potential vorticity q rel (x, y) at initial time t = 0 (left) and at t = 100 days on a regular (center) and an irregular (right) mesh with 2·64 2 triangular cells. Contours between −1.5 days −1 km −1 and 12.5 days −1 km −1 with interval of 1 days −1 km −1 .
well preserved at the order of 10 −8 . All energy error values relate to the time step ∆t = 48 s and reduce further at 1 st -order convergence with smaller ∆t. But they are more or less independent from the spatial resolution (cf. Fig. 7). The indicated marginal trends of loss in total energy, visible in the energy plots of the irregular mesh cases, diminish with higher spatial resolution for all three flow regimes (compare lower and upper rows of Fig. 7).
The relative errors in P E show a similar dependency on the flow regime as the energy errors (hence not shown). Compared to the values presented in Fig. 7 for ∆t = 48 s, the P E error values are about two orders of magnitude larger while they are more or less independent from the time step size. P V is conserved at machine precision for all cases studied (analogously to the time series presented further below for the nonstationary cases). Results of convergence study. We compute L 2 and L ∞ error norms of the numerical solutions for fluid depth D and relative potential vorticity q rel to study the spatial convergence behavior of solutions of the variational RSW integrator. For the steady state case, any deviation of the numerical solutions from the initial fields is considered as numerical error. Hence, we define the L 2 and L ∞ error measures for a discrete function f i (t) with respect to its initial values f i (0) over all triangles T i by These errors are determined on regular and irregular meshes with resolutions of 2 · 32 2 , 2 · 64 2 , 2 · 128 2 , and 2 · 256 2 triangles. Moreover, we use for all simulations one fixed time step of ∆t = 12 s to consider only the spatial convergence behavior. Fig. 8 illustrates the error values of the numerical solutions for D and q rel after 1 day with respect to the corresponding initial states. As q rel differs from q only by f /h, which is almost constant here, the convergence rates for q are very similar to those of q rel (hence not shown). Also not shown are error values for the absolute vorticity ω a which are, too, very similar to those of q rel . Using numerical solutions for later times provided qualitatively the same results, only the absolute error values would be larger. We compare a fluid in semi-geostrophic (left), quasi-geostrophic (middle), and incompressible (right) regimes. Similarly to the conservation properties of the QOI, the absolute error values are the smallest in the incompressible regime while we realize a slightly increase of the errors for semi-geostrophic and quasi-geostrophic flows. Independently from the regime, all solutions for D and q rel show in both error norms convergence rates between 2 nd -and 1 st -order.
Considering the absolute error values for D, one realizes that for all regimes both L 2 and L ∞ errors on irregular meshes are significantly smaller than on regular meshes with same resolution N , i.e. with the same number of cells. In particular for L 2 [D], the error values on an irregular mesh with N = 2( 1 2 N 1D ) 2 cells are close to those with N = 2(N 1D ) 2 cells on a regular mesh. This agrees well with the fact that on irregular meshes the central region has cells with halved triangle edge length providing effectively a doubled resolution.
In contrast to the improvements for D, the error values for q rel are higher on irregular meshes. Here, the irregular cells might trigger numerical noise in form of additional small scale eddies earlier than on regular meshes leading to the increased error values. Nevertheless, the numerical solutions for q rel converge also here with, at least, 1 st -order to the correct solution.

Nonlinear dynamics.
Let us focus next on flows that are dominated by nonlinear processes. By means of two test cases we study whether our variational integrator is able to correctly represent the general dynamical behavior while conserving the quantities of interest discussed above.
In the first test case we study the evolution of two interacting corotating vortices in different regimes, i.e. we study semi-geostrophic, quasi-geostrophic, and incompressible flows. This allows us to examine how accurate the scheme represents flows that are in advection and/or divergence dominated regimes. In the second test case, the time evolution of a shear flow in quasi-geostropic regime is studied. We examine the general flow pattern, such as position and magnitude of the vortex cores in the first test case or the growth rate of the instability in the second test case, and compare the solutions with literature, in particular with [3,17,29]. As above, we apply also here regular and irregular computational meshes for the simulations.

Vortex pair interaction.
In this test case, the flow evolution of two interacting corotating vortices in the inviscid case is studied. This vortex pair problem is described in [17,29,30,33] for instance. Here, we give a brief description of the test case and its initialization according to [17]. Initialization. We choose the initial conditions in geostrophic equilibrium by prescribing the fluid depth h by an analytic solution while determining the velocity by the constraint of geostrophic balance, i.e. f k × u = −g∇h with k = (0, 0, 1) T , cf. [17] for more details. We use the initial fields where we apply for i = 1, 2 the periodic extensions The centers of the vortices and σ x , σ y are given by using o = 0.1 and H = 75 m. As in Sect. 5.2, we map the analytical function (66) onto the mesh by sampling it at each cell center to obtain D i . Considering h as stream function, we initialize the normal velocity values at the cell faces by using the discrete tangential gradient operator (46), i.e. V ij = − g f G ⊥ (h) ij , rather than initializing the velocity directly via (67). This approach leads to discrete fields for fluid depth and velocity that are in geostrophic balance. The initial fields are shown in Fig. 9. Again, we examine the variational scheme with respect to the three flow regimes: (i) H 0 = 450 m for semi-geostrophic, (ii) H 0 = 750 m for quasi-geostropic, and (iii) H 0 = 10 km for incompressible flows. To allow for comparisons to [17], we apply also here a test case with trivial bottom topography.
For all dual vortex simulations, we apply a time step of ∆t = 12 s and use regular and irregular meshes with 2 · 256 2 cells giving a minimal edge length of ∆x min = 1.313 km. Similarly to the convergence test, we obtain hence for the largest water depth of H 0 = 10 km a Courant number C = 2.86. Results. Let us first discuss the flow evolution of the two interacting corotating vortices. For the chosen initial distance between the vortex cores exhibiting an area of negative vorticity in between them (cf. Fig. 9), the vortex cores are too far apart to allow for a merger (see e.g. [4] for more details). Instead, the two cores are mutually repelled due to nonlinear effects. The corresponding time evolution of the relative potential vorticity field is shown in Fig. 10 for the incompressible case, but the flow evolves very similar for all flow regimes studied (cf. Fig. 12).
Comparing in Fig. 10 simulations performed on either a regular (upper row) or an irregular mesh (lower row), we notice that both relative potential vorticity fields agree very well, in particular the speed of the mutual repulsion of the two vortex cores and the thin filaments between them. Also both fields of fluid depth agree very well (not shown). This very good match between the corresponding fields is also given in case of semi-geostrophic and quasi-geostrophic flows (also not shown).
The snapshots of relative potential vorticity and fluid depth at day 10, presented in Fig. 11 for regular meshes, allows us further to compare our simulations with those performed in [3,17]. For both mesh types and for all flow regimes, i.e. semigeostrophic (right column), quasi-geostrophic (middle column), and incompressible (right column), we obtain solutions that are very close to those determined in [17] with a conventional triangular C-grid discretization of the shallow water equations. In particular the fields agree very well when considering the magnitude of q rel and D fields and the position of the two vortex cores. This good match is also given when comparing our results with those in [3] obtained by a corresponding hexagonal C-grid scheme.
The relative errors of the QOI are shown in Fig. 12. The first and second column correspond to the semi-geostrophic case, the third and forth to the quasi-geostrophic case, and the fifth and sixth to the incompressible case. The relative errors in total energy (upper row) for flows in semi-and quasi-geostrophic regimes are at the order of 10 −7 for both mesh types. In the incompressible case, these errors are about two orders of magnitudes smaller, just as in the isolated vortex test case. They are related to the time step ∆t = 12 s and decrease further at 1 st -order when using smaller time step sizes. The error values are very much independent from the spatial resolution and they exhibit the expected oscillatory behavior of a symplectic time integrator.
Here, and for all simulations performed, the mass-weighted potential vorticity P V is conserved at the order of machine precision, just as mass m, for both regular and irregular meshes and for all temporal and spatial resolutions studied. Considering the relative errors of P E on regular meshes, they show a similar dependency on the flow regime as E, with an accuracy at the order of 10 −4 in the semi-and quasigeostrophic regimes and one order of magnitude smaller in the incompressible case. Our RSW scheme conserves these values well, although the variational discretization method does not treat neither P V nor P E as discrete Casimirs. Hence we do not expect them to be strictly conserved. In fact, here for the vorticity dominated vortex interaction test case, we notice a growth of P E on irregular meshes at the order of 10 −3 for a simulation of 10 days. This growth rate is however rather small and in the shear flow test case studied below, even much smaller. Finally we point out that also the error values for P E are more or less independent from spatial and temporal resolutions (not shown).

5.3.2.
Shear flow in semi-geostrophic regime. In the second test case with dominantly nonlinear effects, we study the evolution of a shear flow in quasi-geostrophic regime. The flow is initialized so that it is in unstable equilibrium. During the flow evolution, perturbations that are superimposed on an initial zonal jet in x-direction (cf. Fig. 13) evolve towards two pairs of counter-rotating vortices (cf. Fig. 14). As the development and growth rate of this instabilities essentially depend on nonlinear effects, we thus evaluate the accurate representation of these effects by the nonlinear terms of the variational RSW scheme.
Initialization. An initialization according to [17] will allow us also for this test case a direct comparisons of the corresponding numerical results. Hence, we initialize the shear flow in quasi-geostrophic regime by the following fluid depth h while enforcing the geostrophic balance. There follows  h(x, y, 0) = H 0 − H y σ y e − y 2 2σ 2 y x = x L x , y = 1 π sin π L y y − L y 2 , y = 1 2π sin 2π L y y − L y 2 , for the parameters λ x = 1 2 , σ y = 1 12 , κ = 0.1. To obtain an inviscid flow in quasigeostrophic (compressible) flow regime with Bu ≈ 1, we choose H 0 = 1076 km and H = 30 m, cf. [17]. We set B(x) = 0. The analytic fields are mapped to the mesh exactly as in Sect. 5.3.1. The corresponding initial fields are shown in Fig. 13.
For the shear flow simulations, we apply a time step of ∆t = 36 s and use regular and irregular meshes with 2 · 256 2 cells and ∆x min = 1.313 km. For the water depth of H 0 = 1.0760 km, the Courant number here is C = 2.82.
Results. Fig. 14 shows the time evolution of the shear flow in unstable equilibrium. The initial perturbations grow within the first three days to an instability with a dominant wave number of two. This instability develops further to two pairs of counter-rotating vortices that are well developed at about day six. The filaments between these vortex cores get thinner for later times and reach scales beyond spatial resolution, which causes the noisy pattern visible at day 10.
A comparison of the snapshots of Fig. 14 for regular (upper row) and irregular meshes (lower row) confirms that (i) the scheme is capable to produce accurate solutions even on very deformed, irregular meshes and that (ii) these solutions agree well with literature, i.e. with a triangular [17] and a hexagonal [3] C-grid discretization of the shallow water equations. In more detail, all schemes' solutions Figure 13. Initial fields of fluid depth (left) and relative potential vorticity (right) in geostrophic balance for the shear flow test case on a regular mesh with N = 2 · 256 2 cells. Contours for q rel between −11 days −1 km −1 and 11 days −1 km −1 with interval of 1 days −1 km −1 , and for D between −0.06 km+H 0 and 0.04 km+H 0 with interval of 0.002 km.
for q rel and D provide very similar general flow pattern -compared at days 3, 6, and 10 -consisting in the number of vortex pairs, their magnitudes and positions, and the structure of the thin filaments. Fig. 16 shows the variational integrator's relative errors in total energy (upper row), mass-weighted potential vorticity P V (middle row), and potential enstrophy P E (lower row) for the shear flow simulations in quasi-geostrophic regime performed on a regular (left column) and an irregular (right column) grid. For both mesh types, the variational integrator conserves the total energy at the order of 10 −7 for a time step size of ∆t = 36 s. As above, the errors decreases at 1 st -order rate for smaller ∆t and are independent from the spatial resolution. Again, mass (not shown) and P V are preserved at the order of machine precision, independently from the chosen mesh type or from the spatial and temporal resolutions. Also the errors in P E do not significantly depend on this choice and remain preserved at the order of 10 −3 . Here in case of the shear flow, we do not observe a growth of P E in case of irregular meshes unlike the vortex interaction test case. 6. Conclusions and outlook. This study provides a first step in the development of geometry preserving numerical integrators for compressible fluids. We derived a variational discretization for compressible fluids by extending the variational discretization framework developed in [27] for incompressible ideal fluids, based on a Lie group approximation of the group of volume preserving diffeomorphisms. This extension was achieved by relaxing the volume preserving condition on the discrete diffeomorphisms and imposing appropriate nonholonomic constraints that naturally follow from the relation between the discrete and continuous velocities. Given a semidiscrete Lagrangian, the semidiscrete equations followed by applying the Lagrange-d'Alembert principle in reduced Eulerian form, i.e., the Euler-Poincaré-d'Alembert equations. The resulting semidiscrete equations are valid on any mesh discretization of the fluid domain, in 3D and 2D. We derived them explicitly for 2D irregular simplicial meshes. In particular, a discrete Lie derivative operator was obtained on such meshes. We then specialized our study by focusing on the case of the rotating shallow water equations.  For this case, we numerically verified that our variational discretization (i) preserves stationary solutions of the shallow water equations, (ii) conserves very accurately quantities of interest such as mass, total energy, mass-weighted potential vorticity, and potential enstrophy, (iii) correctly represents nonlinear dynamics, and (iv) correctly represents (inertia-gravity) waves. In more detail, simulating the time evolution of a lake at rest and of a steady isolated vortex, we showed that the shallow  Figure 16. Shear flow test case: relative errors of total energy E(t) (upper row), of mass-weighted potential vorticity P V (t) (middle row), and potential enstrophy P E(t) (lower row) for a fluid in quasi-geostrophic regime for regular (left) and irregular (right) meshes with 2 · 256 2 cells.
water scheme conserves these stationary solutions to a high degree even for longterm simulations. In addition, the numerical solutions of a stationary isolated vortex converge with, at least 1 st -order, towards the exact solutions for both regular or irregular computational meshes. For all test cases studied, mass-weighted potential vorticity and mass were conserved at machine precision while the total energy shows excellent long terms conservation properties with a maximal error at the order of 10 −7 that decreases further at 1 st -order rate for smaller time step sizes. Moreover, the scheme correctly represents (inertia-gravity) waves as a comparison between numerically determined and theoretically predicted wave spectra confirmed. The study of the nonlinear dynamics with respect to a dual vortex interaction and a shear flow test case showed further that the scheme presents very accurately the dynamics triggered by nonlinear interaction. The quality of the results is similarly good for regular and irregular computational meshes. The correctness of our simulations are underpinned by a comparison to literature. Providing here a variational integrator for the two dimensional shallow water equations, object of current and future work is the extension of the variational discretization framework to derive structure-preserving discretizations for fully three dimensional compressible flows.