Werner Bauer

Jörn Behrens

Colin J Cotter

A structure-preserving split finite element discretization of the rotating shallow water equations in split Hamiltonian form

. The splitting introduces two chains of compatible FE spaces such that the differential operators hold strongly. This leads to a separation of the set of equations into topological prognostic and metric-dependent closure equations. Consequently, the structure preservation is related to the topological equations, namely the conservation of the Hamiltonian (i.e. energy) follows from the antisymmetry of the Poisson bracket and the conservation of mass, potential vorticity and potential enstrophy from the bracket's Casimirs. These structure-preserving properties are not affected by the metric closure equations that, in turn, carry metric information and are therefore responsible for the schemes' accuracy, stability, convergence and discrete dispersion properties. We verify the clear separation of structure preservation and numerical properties analytically and show numerical examples of geophysical relevance for a low order single mesh implementation.

Introduction

Compatible FE methods for geophysical fluid dynamics allow for derivations of structure-preserving discretizations that preserve important invariants of the continuous equations of motion. Conservation of such invariants are of particular interest for geophysical flow models in order to avoid, for instance, biases in the statistical behavior of the solutions in case certain invariants are not conserved by the discrete schemes, cf. [START_REF] Dubinkina | Statistical mechanics of Arakawa's discretizations[END_REF][START_REF] Bendall | Statistical properties of an enstrophy conserving discretisation for the stochastic quasi-geostrophic equation[END_REF].

These FE methods apply suitable compatible FE approximations of the variables in use. In the framework of finite element exterior calculus (FEEC) [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF][START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF], the combination of functional analysis and differential geometry permits to find suitable FE spaces while avoiding the problem of instabilities that might occur in conventional mixed FE methods. In the context of geophysical fluid dynamics, several contributions suggest corresponding compatible FE discretizations of various geophysical flow models such as the rotating shallow water (RSW) equations, e.g. [START_REF] Cotter | A finite element exterior calculus framework for the rotating shallow-water equations[END_REF][START_REF] Cotter | Compatible finite elements for numerical weather prediction[END_REF].

Both standard FE and FEEC methods apply integration by parts to address the regularity properties of the FE spaces in use. However, this introduces additional errors and certain operators, such as the co-derivative, to be non-local. [START_REF] Lee | Local coderivatives and approximation of Hodge Laplace problems[END_REF] address this issue by localizing the bilinear form in order to ensure a local volume preservation of the quantities of interest. Alternatively, [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF] introduced a novel approach based on the split form of the equations of GFD [START_REF] Bauer | A new hierarchically-structured n-dimensional covariant form of rotating equations of geophysical fluid dynamics[END_REF], in which pairs of FE spaces are used such that integration by parts is avoided.

In this manuscript, we extend the split FE framework of [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF] to cover the nonlinear RSW equations. We develop a FE discretization method that explicitly uses the split form (topological and metric parts) of the equations of GFD [START_REF] Bauer | A new hierarchically-structured n-dimensional covariant form of rotating equations of geophysical fluid dynamics[END_REF]. The particularity of the approach is that this splitting is preserved during the discretization procedure, resulting in topological prognostic and metric-dependent diagnostic closure equations that both can be expressed in matrix-vector form.

Moreover, the split FE framework will be combined with the Hamiltonian framework which permits one to derive structure-preserving discretizations systematically from a Hamiltonian functional and a Poisson bracket [START_REF] Bauer | Energy-enstrophy conserving compatible finite element schemes for the shallow water equations on rotating domains with boundaries[END_REF]. In fact, the split equations of GFD follow from a split Hamiltonian form in which the topological structure of the equations is encoded in a topological Poisson bracket and the metric information is carried by a Hamiltonian [START_REF] Eldred | Variational and Hamiltonian Formulations of Geophysical Fluids using Split Exterior Calculus[END_REF]. Here, we do not focus on this latter point but rather on the resulting structure-preserving discretizations and their properties.

There are various approaches of structure-preserving discretization methods (see e.g. [START_REF] Bauer | Variational integrators for anelastic and pseudoincompressible flows[END_REF][START_REF] Bauer | Towards a variational discretization of compressible fluids: the rotating shallow water equations[END_REF][START_REF] Desbrun | Variational discretization for rotating stratified fluids[END_REF][START_REF] Gawlik | Geometric, variational discretization of continuum theories[END_REF]). Our method using differential forms shares in particular some basic ideas with mimetic discretizations (see e.g. [START_REF] Da Veiga | Mimetic finite difference methods for Hamiltonian wave equations in 2D[END_REF][START_REF] Bochev | Principles of mimetic discretizations of differential operators, Compatible Spatial Discretizations[END_REF][START_REF] Bossavit | Generalized finite differences" in computational electromagnetics[END_REF][START_REF] Desbrun | Discrete exterior calculus[END_REF][START_REF] Palha | Physicscompatible discretization techniques on single and dual grids, with application to the Poisson equation of volume forms[END_REF]) in which the PDEs are sometimes formulated by differential forms. Although some of these formulations even introduce clear distinction between topological and metric parts [START_REF] Kreeft | Mimetic spectral element method for generalized convection-diffusion problems[END_REF][START_REF] Cotter | A finite element exterior calculus framework for the rotating shallow-water equations[END_REF], none of them associates a proper FE space to each variable, as suggested by our framework.

The paper is organized as follows. In Section 2, a one-dimensional (1D) y-independent slice model of the nonlinear RSW equations is introduced. It is formulated such that an extension to higher dimensional tensor products is straightforward. In Section 3 the Hamiltonian structure with Poisson bracket for the 1D slice model is introduced and a split FE discretization for dual pairs of compatible FE spaces with arbitrary order is introduced. Proofs for the conservation of energy and Casimirs (mass, potential vorticity and enstrophy) are given. As a concrete realization of this general framework, we provide in Section 4 a family of low order split FE schemes and we discuss their properties, in particular how they relate to the topological and metric structures. In Section 5 we perform numerical simulations to test convergence, consistency, and accuracy of the split schemes. In Section 6 we draw conclusions and provide an outlook for ongoing and future work.

A split 1D slice model of the RSW equations

Based on the split form of the equations of GFD [START_REF] Bauer | A new hierarchically-structured n-dimensional covariant form of rotating equations of geophysical fluid dynamics[END_REF], we introduce a 1D y-independent slice model of the RSW equations (e.g. [START_REF] Staniforth | Analysis of a mixed finite-element pair proposed for an atmospheric dynamical core[END_REF]) and discretize it using an extension of the split FE discretization method of [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF]. We do not specify here the FE spaces concretely nor their order of accuracy, only their mutual relations in terms of exact sequences connected by the Hodge star is required. It is only in Section 4 where we choose concrete realizations of FE spaces.

A split 1D slice RSW model

We introduce a 1D slice model of the rotating shallow-water equations on an f -plane in split form, i.e. we apply straight and twisted differential forms and separate the equations into topological and metric terms, according to [START_REF] Bauer | A new hierarchically-structured n-dimensional covariant form of rotating equations of geophysical fluid dynamics[END_REF]. As derived in Appendix A, this model reads ∂u (1) ∂t -˜ q (0) F (0) v + d B (0) = 0, ∂ v (1) ∂t + ˜ q (0) F (0) u = 0, ∂ h (1) ∂t + d F (0) u = 0, (2.1) u (0) = u (1) , v (0) = v (1) , h (1) = h (0) , (2.2) using the definitions F (0) u := h (0) u (0) , F (0) v := h (0) v (0) (mass fluxes), and (2.3) B (0) := gh (0) + 1 2 ( u (0) ) 2 + 1 2 (v (0) ) 2 (Bernoulli function), and (2.4) q (0) h (1) = dv (0) + f dx (potential vorticity (PV)), (2.5) where u (1) (x, t) = u(x, t)dx and v (1) (x, t) = v(x, t) dx are straight, respectively twisted velocity 1forms with coefficient functions u(x, t) and v(x, t), respectively, and h (0) (x, t) is the fluid height with coordinate x in x-direction and t for time. f is the constant Coriolis parameter on the f -plane.

For later discussions we introduce the phase velocity c = √ gH. Note that for this manuscript we employ periodic boundary conditions (BCs) h (0) (0, t) = h (0) (L, t), u (0) (0, t) = u (0) (L, t), v (0) (0, t) = v (0) (L, t), ∀t ∈ [0, T ] where T denotes the total integration time.

The corresponding twisted quantities for velocity and fluid height are defined with respect to the twisted Hodge-star operator : Λ k → Λ (1-k) (resp. Λ k → Λ (1-k) ) that maps from straight (resp. twisted) k-forms to twisted (resp. straight) (1 -k)-forms with k = 0, 1, in 1D. We denote the space of (time-dependent) straight one-forms with (Λ 1 ; 0, T ) and of twisted one-forms such as h (1) with ( Λ 1 ; 0, T ); but usually we simply write Λ 1 or Λ 1 . Similarly, h (0) (x, t) ∈ (Λ 0 ; 0, T ) is a straight 0-form (straight function) and u (0) (x, t) ∈ ( Λ 0 ; 0, T ) a twisted 0-form. The index (k) denotes the degree, and Λ k , Λ k the space of all k-forms. Recall that straight forms do not change their signs when the orientation of the manifold changes in contrast to twisted forms.

Finally, the exterior derivative d is a mapping d :

Λ k → Λ k+1 (or d : Λ k → Λ k+1 ).
Here in 1D, it is simply the total derivative of a smooth function [START_REF] Bauer | A new hierarchically-structured n-dimensional covariant form of rotating equations of geophysical fluid dynamics[END_REF] for full details. Diagram (2.6) illustrates the relations between the operators and spaces.

g (0) ∈ Λ 0 , d g (0) = ∂ x g(x)dx ∈ Λ 1 , cf.
h (0) , v (0) ∈ Λ 0 d ----→ Λ 1 u (1) ˜     ˜ h (1) , v (1) ∈ Λ 1 d ← ----Λ 0 u (0) (2.6)
The split covariant formulation requires the topological equations (2.1) not to mix straight and twisted differential forms (DF); they are related only via the metric equations (2.2), cf. [START_REF] Bauer | A new hierarchically-structured n-dimensional covariant form of rotating equations of geophysical fluid dynamics[END_REF][START_REF] Bossavit | Generalized finite differences" in computational electromagnetics[END_REF]. This requirement is clearly satisfied by both momentum and continuity equations.

Remark 2.1. The covariant split form of the 1D slice model is equivalent to standard formulations as introduced in Appendix B. This easily follows when fixing an orientation of L, representing all 1-forms in coordinates and using the flat operator (cf. [START_REF] Bauer | A new hierarchically-structured n-dimensional covariant form of rotating equations of geophysical fluid dynamics[END_REF]) to map the covariant to the vector-invariant form; here in 1D this is simply to omit the basis dx and consider only the coefficient functions.

Remark 2.2. The combination of a straight and a twisted DF to describe u (1) and v (1) , respectively, for our 1D setup matches well the regularity conditions that are required for higher dimensional tensor product FE methods such that they represent schemes on staggered meshes. For instance in 2D and low order FE, this would result in a velocity that is in H(div), cf. [START_REF] Eldred | Dispersion analysis of compatible Galerkin schemes on quadrilaterals for shallow water models[END_REF].

In general, the split covariant equations, in particularly the use of straight and twisted DF, naturally suggest suitable and compatible FE approximation spaces for the quantities of interest, according to the split FEM approach.

Galerkin discretization of the split slice model

In standard FE methods, the equations of motion are usually formulated first in weak form to perform integration by parts. This is required to guarantee the existence of the derivatives of all functions.

As the split FEM avoids integration by parts, we trivially project the split equations (2.1) on the corresponding FE spaces which gives immediately discrete equations in matrix-vector form.

To this end, we introduce the compatible finite element spaces

H 1 d ----→ L 2   π 0   π 1 Λ 0 h , Λ 0 h d ----→ Λ 1 h , Λ 1 h (2.7) Λ 0 h d ----→ Λ 1 h ˜ 0 h     ˜ 1 h Λ 1 h d ← ----Λ 0 h (2.8)
with commuting, bounded, surjective projections (π 0 , π 1 ). Here for the (1D) tensor product construction, we consider

Λ 0 h , Λ 0 h = CG k and Λ 1 h , Λ 1 h = DG k-1
with polynomial order k. The discrete Hodge star operators ˜ 0 h , ˜ 1 h map between straight and twisted spaces and are allowed to be non-invertible.

Definition 2.3. The split FE discretization of the split equations (2.1)-(2.5) seeks solutions u

(1)

h , v (1)
h , h

h ∈ (Λ 1 h (L), Λ 1 h (L), Λ 1 h (L) (1)
) of the topological equations (in form of trivial projections)

χ (1) h , ∂ ∂t u (1) h -χ (1) h , ˜ q (0) h F (0) v h + χ (1) h , d B (0) h = 0, ∀χ (1)
h ∈ Λ 1 h , (2.9)
χ (1) h , ∂ v (1) h ∂t + χ (1) h , ˜ q (0) h F (0) u h = 0, ∀ χ (1) h ∈ Λ 1 h , (2.10) χ (1) h , ∂ h (1)
h ∂t + χ (1) h , d F (0) u h = 0, ∀ χ (1)
h ∈ Λ 1 h , (2.11) ˜ φ (0) h , q (0) h h (1) h
-˜ φ (0) h , d v (0) h -˜ φ (0) h , f dx = 0, ∀ φ (0) h ∈ Λ 0 h , (2.12)
subject to the metric closure equations (in form of non-trivial Galerkin projections (GP))

u (0) h = u 1-i u (1) h by u (0) h (x), τ i (x) = u (1) h (x), τ i (x) , ∀τ i ∈ Λi h , i = 0, 1, v (0) h = u 1-i v (1) h by v (0) h (x), τ i (x) = v (1) h (x), τ i (x) , ∀τ i ∈ Λi h , i = 0, 1, h (0) h = h 1-j h (1) h by h (0) h (x), τ j (x) = h (1) h (x), τ j (x) , ∀τ j ∈ Λj h , j = 0, 1, (2.13)
with pointwise definitions of fluxes and Bernoulli functions as in (2.3) and (2.4) or with projections as defined further below. Note that •, • := L • ∧ ˜ • is the L 2 the inner product on the domain L. Moreover, (2.13) are nontrivial Galerkin projections between the coefficient functions of zero and one-forms, hence the notation using x. They represent discrete (potentially non-invertible) versions of the Hodge-star operator and are realized exactly as in [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF].

Remark 2.4. We distinguish between discrete and continuous Hodge star operators. The discrete ones in Eqn. (2.13) are realized as Galerkin projections between 0-and 1-forms and we denote them in the following frequently by ˜ h . The continuous Hodge star ˜ is used in the inner product in order to guarantee that k-forms of the same degree are multiplied. The latter translates, for instance, to a standard inner product of coefficient functions, i.e. ˜ u

(1) h , u (0) h = u (1) h (x), u (0) h (x) .

Discretization via the split (Hamiltonian) FE method

To discretize the above 1D slice model in a structure preserving manner while keeping the splitting conserved, we combine here the Hamiltonian framework with the split FE framework of [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF]. Namely, we dedicate to both straight and twisted DF proper FE spaces such that both momentum and continuity equations hold strongly while the metric equations are realized by approximating the Hodge star operator by nontrivial Galerkin projections [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF]. Casting the split equations in Hamiltonian form will allow us to identify the invariants of the equations, such as energy, mass, potential vorticity and enstrophy from the structure of the equations consisting of a Hamiltonian functional and an almost Poisson bracket.

Hamiltonian structure of discrete split nonlinear equations

To find structure preserving discretizations of the split equations (2.1)-(2.5), we apply finite element methods on a formulation that combines the Hamiltonian formulation with the splitting of the equations.

Definition 3.1. Let H : Λ 1 × Λ 1 × Λ 1 → R be the Hamiltonian functional of the 1D split RSW model, defined by H[u (1) h , v (1)
h , h (1)
h ] = 1 2 u (1) h , ˜ h (0) h u (0) h + 1 2 v (1) h , ˜ h (0) h v (0) h + h (1)
h , ˜ gh

(0) h , (3.1)
with metric equations:

u (0) h = u (0) h [u (1)
h ] = ˜ h u (1) h , v (0)
h = v (0) h [ v (1)
h ] = ˜ h v (1) h , h (0)
h = h (0) h [ h (1)
h ] = ˜ h h (1) h , (3.2)
in which the squared brackets [ ] indicate the dependency of a function from another function. We define the (almost) Poisson bracket {, } as

{F, G} := - δF δ h (1) h , d ˜ δG δu (1) h - δF δu (1) h , d ˜ δG δ h (1) h + δF δu (1) h
, ˜ q

(0) h ˜ δG δ v (1) h - δF δ v (1) h , ˜ q (0) h ˜ δG δu (1) h (3.3)
with potential vorticity (PV) q

(1)

h defined implicitly by ˜ φ (0) h , q (0) h h (1) h + d φ (0) h , v (1) h
-˜ φ (0) h , f dx = 0, ∀ φ (0) h ∈ Λ 0 h , (3.4)
which corresponds to the weak form of equation (2.12). Then, the dynamics for any functional F :

Λ 1 × Λ 1 × Λ 1 → R is given by the almost Poisson bracket formulation d dt F[u (1)
h , v (1)
h , h

h ] = {F, H}.

The discrete Hodge star operators ˜ h in (3.2) are realized by the nontrivial Galerkin projections (2.13). In order to proof this proposition, we have to calculate the functional derivative of H. For instance, the functional derivatives of the functional F[u

(1) h ] with respect to u (1)
h is defined weakly for an arbitrary test functions χ (1)

∈ Λ 1 h by δF[u (1)
h ; χ (1) ] = ˜ h χ (1) , ˜ δF δu

(1) h := lim →0 1 F[u (1)
h + χ (1) ] -F[u

h ] (3.5) in which we used the L 2 inner product f, g := L dx f • g for some functions (0-forms) f, g; which is equivalent to

•, • := L • ∧ ˜ •.
The latter provides an inner product also for 1-forms.

Lemma 3.3. The functional derivatives of H with respect to u

(1)

h , v (1)
h , h

h for the test function χ

(1)

h ∈ Λ 1 h , χ (1)
h ∈ Λ 1 h read δH[u (1)
h ; χ

h ] = ˜ h χ (1) h , ˜ δH δu (1) h = ˜ h χ (1) h , h (0) h u (0) h = ˜ h χ (1) h , F (0) u h , δH[ v (1)
h ; χ

h ] = ˜ h χ (1) h , ˜ h δH δ v (1) h = ˜ h χ (1) h , h (1)
h v (0) h = ˜ h χ (1) h , F (0) v h , δH[ h (0)
h ; χ

h ] = ˜ h χ (1) h , ˜ h δH δ h (1) h = ˜ h χ (1)
h , gh h and using the symmetry of the inner product, we find

(0) h + 1 2 u (0) h 2 + 1 2 v (0) h 2 = ˜ h χ (1) h , B (0)
δH[u (1) h ; χ (1) ] = 1 2 ˜ h δu (1)
h , h

(0) h u (0) h + 1 2 ˜ h u (1) h , h (0)
h δ u (0) h (3.6) = 1 2 ˜ h χ (1) , h (0) h u
(0) h + 1 2 ˜ h u (1) h , h (0)
h ˜ h χ (1) (3.7) = 1 2 ˜ h χ (1) , h (0) h u
(0) h + 1 2 ˜ h χ (1) , h (0)
h ˜ h u (1) h (3.8)
= ˜ h χ (1) , h

h u (0) h , (0)
noting that δ u

(0) h [u (1)
h ] = lim →0 1 ( u (0) h [u (1)
h + χ (1) ] -u

(0) h [u (1)
h ]) = ˜ h χ (1) . To map ˜ h u

(1)

h to u (0)
h from line three to line four, we use the metric closure equations (2.13), i.e. for i = 1 this map reads χ

(1)

h (x), u (1)
h (x) = χ (1) h (x), u (0) h (x) ∀χ (1) h (x) ∈ Λ 1
h and for i = 0 it is given by φ

(0) h (x), u (1)
h (x) = φ (0) h (x), u (0) h (x) ∀φ (0) h (x) ∈ Λ 0
h . These discrete versions of a (non-invertible) Hodge star map elements from the space Λ

1 h → Λ 0 h | ker , ˜ h : χ (1) h (x) → φ (0) h (x)| ker(˜ h ) i.e
. from the DG k-1 to the CG k space in which the kernel of this mapping had been removed. This mapping is applied to both the velocity and the test function, hence the functional derivative holds pointwise in φ

(0) h (x)| ker(˜ h ) . The functional derivative δH[ v (1)
h ; χ (1) ] follows from exactly the same line of reasoning. Next, we derive the functional derivative of H with respect to h [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] h for the test functions χ (1) ∈ Λ 1 h . Similarly to above, it can be calculated from the first variation of H and the symmetry of the inner product like

δH[ h (1) h ; χ (1) ] = 1 2 ˜ h u (1)
h , δh

(0) h u (0) h + 1 2 ˜ h v (1)
h , δh

(0) h v (0) h + 1 2 ˜ h δ h (1)
h , gh

(0) h + 1 2 ˜ h h (1)
h , gδh 1) , gh

(0) h = 1 2 ˜ h u (1) h , ˜ h χ (1) u (0) h + 1 2 ˜ h v (1) h , ˜ h χ (1) v (0) h + 1 2 ˜ h χ (
(0) h + 1 2 ˜ h h (1)
h , g˜ h χ (1) = ˜ h χ (1) ,

1 2 ˜ h u (1) h u (0) h + 1 2 ˜ h v (1) h v (0) h + ˜ h χ (1) , 1 2 gh (0) h + 1 2 g˜ h h (1) h = ˜ h χ (1) , 1 2 u (0) h u (0) h + 1 2 v (0) h v (0) h + gh (0) h (3.10) with δh (0) h [ h (1)
h ] = lim →0 1 (h (0)
h [ h (1) + χ (1) ] -h

(0) h [ h (1)
h ]) = ˜ h χ (1) . From the third to the forth line we used again the metric closure equations, now for both velocity and height fields to project ˜ h h (1)

h to h (0) h and ˜ h u (1) h to u (0)
h . Similarly to above, the functional derivative holds pointwise in φ (0) h (x)| ker(˜ h ) for any choice of metric equations. This concludes the proof.

Proof of Proposition 3.2. To derive the continuity equation, consider the functional

F[ h (1) h ; χ (1) h ] = L χ (1) h ∧ ˜ h (1)
h with test function χ [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] h and functional derivatives

δF[ h (1) h ; χ (1) h ] δ h (1) h = χ (1) h , δF[ h (1) h ; χ (1) h ] δu (1) h = 0.
The dynamics of F follows from the Poisson bracket (3.3)

d dt F[ h (1)
h ] = - δF δ h (1) h , d ˜ δH δu (1) h ⇒ χ (1) h , ∂ ∂t h (1) h = -χ (1) h , d F (0) u h ∀ χ (1) h , (3.11)
which gives the continuity equation (2.11). Analogously, to derive the momentum equations we consider

F[u (1)
h ; χ

h ] = χ (1) h , u (1)
= L χ (1) h ∧ ˜ u (1) h
h dx and F[ v (1)
h ; χ

h ] = χ (1) h , v (1)
= L χ (1) h ∧ ˜ v (1) h
h dx with functional derivatives

δF[u (1) h ; χ (1) h ] δ h (1) h = 0, δF[u (1)
h ; χ

h ] δu

(1) h = χ (1) h , δF[ v (1) h ; χ (1) h ] δ v (1) h = χ (1) h , δF[u (1) h ; χ (1) h ] δ v (1) h = δF[ v (1)
h ; χ

h ] δu

(1) h = 0.
Then, the dynamics of the functionals read

d dt F[u (1)
h ] = - δF δu (1) h , d ˜ δH δ h (1) h
+ δF δu (1) h , ˜ q (0) h ˜ δH δ v (1) h (3.12)
⇒ χ

(1)

h , ∂ ∂t u (1) h = -χ (1) h , d B (0) h + χ (1) h , ˜ q (0) h F (0) v , (3.13)
d dt F[ v (1)
h ] = - δF δ v (1) h , ˜ q (0) h ˜ δH δu (1) h (3.14) ⇒ χ (1) h , ∂ ∂t v (1) h = -χ (1) h , ˜ q (0) h F (0) u , (3.15)
yielding the momentum equations (2.9) and (2.10), respectively. Considering also the remaining closure equations, the bracket formulation results in the equations of motion in Definition 2.3.

Properties associated to the topological and metric parts

We discuss the schemes properties in relation to the splitting of the schemes into topological and metric parts.

Properties of the split schemes associated to the topological part

We first discuss the topology dependent properties that are related to structure preservation.

Proposition 3.4. All split schemes preserve total energy independently from the metric equations.

Proof. The Poisson bracket (3.3) does not depend on the discrete Hodge star ˜ h and is therefore antisymmetric for all choices of metric closure equations. Then, the energy functional H is invariant in time:

d dt H[u (1) h , v (1)
h , h

h ] = {H, H} = 0. (3.16) Proposition 3.5. The functionals C = h (1) h , ˜ F ( q (0) h ) with straight function F [ q (0) (1)
h ] depending on the twisted function q (0)

h for F = 1 (mass M), F = q (0)
h 1 (mass-weighted potential vorticity (PV)), and F = ( q (0) h ) 2 potential enstrophy (PE) are Casimirs of the Poisson bracket and conserved for all split schemes.

Proof. We start by showing that C = h

(1) h q (0) h , ˜ q (0) h
is a Casimir of the bracket in (3.3), i.e. PE is conserved. From the definition of potential vorticity in (3.4) and assuming δf = 0, there follows its first variation δ( q

(0) h h (1) h ), ˜ φ (0) h = -δ v (1) h , d φ (0) h ∀ φ (0) h ∈ Λ 0 h . (3.17)
Note that we have chosen the weak form (3.4) with respect to the velocity instead of (2.12) in order to allow to calculate the functional derivatives of the prognostic variables v

h and h

h . Then, for the first variations of C, we find

δC = ˜ 2δ q (0) h , q (0) h h (1) h + ˜ ( q (0) h ) 2 , δ h (1) h (3.18) = ˜ 2 q (0) h , δ( q (0) h h (1) h ) -˜ ( q (0) h ) 2 , δ h (1) h (3.19) = -2 d q (0) h , δ v (1) h -˜ ( q (0) h ) 2 , δ h (1) h (3.20) using ˜ 2 q (0) h , δ( q (0) h h (1) h ) = ˜ 2δ q (0) h , q (0) h h (1) h + ˜ 2( q (0) h ) 2 , δ h (1) h
in the second line and (3.17) in the last line while choosing q

(0) h = φ (0)
h as test function. Consequently, the functional derivatives of C are given by δC δ v

= -2 d q (0) h , χ (1) h , χ (1) h
h ∈ Λ 1 h while δC δu (1) h ∀ χ (1)
and δC δ h

(1) h , χ (1) h
= -˜ ( q (0) h ) 2 , χ (1) h ∀ χ
(1)

h ∈ Λ 1 h . (3.22)
Then, C is a Casimir, i.e. the Poisson bracket vanishes for any functional G:

{C, G} = - δC δ h (1) h , d ˜ δG δu (1) h - δC δu (1) h , d ˜ δG δ h (1) h =0 + δC δu (1) h , ˜ q (0) h ˜ δG δ v (1) h =0 - δC δ v (1) h , ˜ q (0) h ˜ δG δu (1) h = ˜ ( q (0) h ) 2 , d ˜ δG δu (1) h + 2 q (0) h d q (0) h , ˜ ˜ δG δu (1) h = -d( q (0) h ) 2 , ˜ ˜ δG δu (1) h + d( q (0) h ) 2 , ˜ ˜ δG δu (1) h = 0, (3.23)
noting that in 1D we can reformulate the action of the exterior derivative on functions

f ∈ Λ 0 h , g ∈ Λ 0 h like ˜ f, d g = f dx, ∂ x gdx = f, ∂ x g = -∂ x f, g = -d f, g
where we used ˜ 1 = dx. Consequently, {C, G} = 0 for any G, and in particular for H and, therefore, Ċ = {C, H} = 0. Analogously we show that PV, defined by C = h

(1)

h q (0) h , ˜ 1 , is a Casimir of the bracket in (3.
3). This follows immediately from setting in (3.17 h , ˜ 1 . These results hold for all split schemes as the bracket are independent of ˜ h and as the metric equations are not involved in the proof. Remark 3.6. A direct proof of the conservation of Casimir without referring to the Hamiltonian structure and the Poisson brackets applies directly the equations of motions (2.10) for v [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] h and (2.11) for h

(1) h . To this end, we calculate the time derivative of (3.4) noting that ∂ t f = 0:

∂ t ( q (0) h h (1) h ), ˜ φ (0) h = -∂ t v (1) h , d φ (0) h ∀ φ (0) h ∈ Λ 0 h . (3.24) Then, C = h (1) h q (0) h , ˜ q (0) h is conserved, because ∂ t C = ˜ 2∂ t q (0) h , q (0) h h (1) h + ˜ ( q (0) h ) 2 , ∂ t h (1) h (3.25) = ˜ 2 q (0) h , ∂ t ( q (0) h h (1) h ) -˜ ( q (0) h ) 2 , ∂ t h (1) h (3.26) = -2 d q (0) h , ∂ t v (1) h -˜ ( q (0) h ) 2 , ∂ t h (1) h (3.27) = + 2 d q (0) h , ˜ q (0) F (0) u + ˜ ( q (0) h ) 2 , d F (0) u (3.28) = + d( q (0) h ) 2 , ˜ F (0) u + ˜ ( q (0) h ) 2 , d F (0) u (3.29) = + d( q (0) h ) 2 , ˜ F (0) u -d( q (0) h ) 2 , ˜ F (0) u = 0. (3.30)
From line three to line four we used (2.10) for v

h and (2.11) for h

h that hold pointwise, while from line five to line six we applied integration by parts (see proof above).

In sum, structure preservation, in particular conservation of energy, mass, potential vorticity and enstrophy, can be associated to the bracket structure of the Hamiltonian formulations which is independent from the realization of the metric equations, hence it holds for all split schemes. In general, the topological equations incorporate the compatibility conditions of the spaces to each other which are valid for a whole class of compatible FE spaces, but they do not require an actual choice of spaces. The metric equations, in contrast, are actual realizations of these spaces and provide the properties of the schemes associated with metric, as discussed next.

Properties of the split schemes associated to the metric equations

Whereas the topological properties hold for all FE spaces that fulfill the commuting diagram (2.8), the metric properties are associated to a certain choice of these spaces. In this paper, we do not address metric properties that hold for a whole class of spaces but rather study only one family of P0-P1 low-order schemes (cf. Section 4). Therein, we also provide the corresponding metric dependent properties of each family member for these low order P0/P1 realizations. Here, we present only a list of properties that will hold for any split FE scheme and explain their dependency on the metric.

The following properties depend on the metric structure:

• Dispersion relation: the dispersion relation depends on the metric because the physical distance between DoF and their mutual positioning determines the wave dispersion on grids. For instance, the dispersion relations shown in Figure 4.2 all explicitly depend on ∆x;

• Stability: the stability of discretizations can be determined in terms of the inf-sub stability condition, cf. [START_REF] Boffi | Mixed finite element methods and applications[END_REF], that depends on the norm and hence on the metric;

• Convergence and accuracy: the choice of FE spaces determines both convergence and accuracy of the approximated functions (and hence of the entire scheme), all measured in terms of norms.

Alternative model and structure-preserving approximations

As an alternative to Definition (3.1), we introduce a "smoothed" Hamiltonian

H sm [u (1) h , v (1)
h , h

h ] = 1 2 u (1) h , ˜ h (1) h (x) u (0) h + 1 2 v (1) h , ˜ h (1) h (x)v (0) h + h (1) h , ˜ gh (0) h , (1)
in which we use the DG k-1 regular coefficient function h

(1) h (x) instead of h (0)
h . This gives the functional derivatives of H sm with respect to u (1)

h and v (1) h as ˜ h χ (1) h , ˜ δH sm δu (1) h = ˜ h χ (1) h , h (1) h (x) u (0) h = ˜ h χ (1) h , F (0) u h | sm ∀˜ h χ (1) h , ˜ h χ (1) h , ˜ δH sm δ v (1) h = ˜ h χ (1) h , h (1)
h (x)v (0) h = ˜ h χ (1) h , F (0) v h | sm ∀˜ h χ (1)
h .

Note that in contrast to above, here the smoothed fluxes F (0)

u h | sm and F (0) v h | sm are defined only weakly in Λ 0 h | ker or Λ 0 h | ker .
The functional derivative with respect to h

h reads either

˜ h χ (1) h , ˜ δH sm δ h (1) h = ˜ h χ (1)
h , gh

(0) h + 1 2 u (0) h 2 + 1 2 v (0) h 2 = ˜ h χ (1) h , B (0) h
as above or, alternatively when projecting u (0)

h to ˜ h u (1)
h in Eqn. (3.10), it is given by

˜ h χ (1) h , ˜ δH sm δ h (1) h = ˜ h χ (1) h , gh (0) h + 1 2 u
(1)

h (x) 2 + 1 2 v (1) h (x) 2 = ˜ h χ (1) h , B (0)
h | sm . (3.32)
In the latter smoothed Bernoulli function, we apply the DG k-1 regular coefficient functions u

χ (1) h , ∂ ∂t u (1) h -χ (1) h , ˜ q (0) h F (0) v h | sm + χ (1) h , d B (0) h | sm = 0, ∀χ (1)
h ∈ Λ 1 h , (3.33) χ (1) h , ∂ v (1)
h ∂t + χ (1) h , ˜ q (0) h F (0) u h | sm = 0, ∀ χ (1) h ∈ Λ 1 h , (3.34)
χ (1) h , ∂ h (1) h ∂t + χ (1) h , d F (0) u h | sm = 0, ∀ χ (1) h ∈ Λ 1 h , (3.35) φ (0) h , h (1)
h (x) u (0) h -φ (0) h , F (0) u h | sm = 0, ∀ φ (0) h ∈ Λ 0 h , (3.36) φ (0) h , h (1)
h (x)v (0) h -φ (0) h , F (0) v h | sm = 0, ∀φ (0) h ∈ Λ 0 h , (3.37)
with q (0) h defined via (2.12) or (3.17) and subject to the variational metric equations (2.13). Note that the fluxes hold weakly in the function spaces Λ 0 h | ker , Λ 0 h | ker in which the kernel of ˜ h has been removed. Remark 3.7. Note that Propositions 3.4 and 3.5 also hold for the smoothed model, i.e. the energy given by H sm of (3.31) is conserved and likewise the Casimirs. The former follows immediately by the antisymmety of the Poisson bracket which has not been altered when using the smoothed fluxes or Bernoulli function. The latter is a direct consequence of the proof to 3.5 which does not apply concrete realizations of the functional derivatives of the Hamiltonian, hence is valid also for the smoothed model.

Family of P0-P1 low-order split FE schemes

Here, we provide realizations of the split schemes discussed above in form of P0-P1 low-order split FE discretizations. This can be seen as an extension of the split FE schemes of the linear wave equations introduced in [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF] to the nonlinear, rotating case. As such, we proceed similarly as therein, use the same notation, while here, we provide only a brief summary of the discretization steps for completeness.

We use FE spaces V 0 h with a piecewise linear basis {φ l (x)} N l=1 and V 1 h with a piecewise constant basis {χ m (x)} N m=1 to approximate the coordinate functions of the continuous 0-and 1-forms; e.g. we approximate 1-form coordinate functions as u (1) h (x, t) = N m=1 u m (t)χ m (x) and 0-form coordinate functions as u (0) h = N l=1 u l (t)φ l (x); we proceed analogously for the other quantities. Note that the subindex l denotes degrees of freedom (DoFs) at the nodes and m at the cell centers of the mesh. As we are on a one-dimensional domain with periodic boundaries, both have N independent DoFs such that l, m = 1, ..., N , cf. [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF].

The resulting mass and stiffness matrices are defined as follows. M nn and M ee are (N × N ) mass-matrices with metric-dependent coefficients

(M nn ) ll = L φ l (x)φ l (x)dx, (M ee ) mm = L χ m (x)χ m (x)dx. (4.1)
Later, we will use that fact that M ee = Id ee (∆x e ) T , in which () T denotes the transpose of a matrix.

The vector ∆x e = (∆x 1 , . . . ∆x m , . . . ∆x N ) containing metric information of the mesh. M ne is a (N × N ) matrix with metric-dependent coefficients

(M ne ) ml = L χ m (x)φ l (x)dx, (M en ) lm = L φ l (x)χ m (x)dx, (4.2)
with M en = (M ne ) T . Similarly to above, we separate M ne = P ne (∆x e ) T , into the metric-dependent part ∆x e and a metric-free part P ne which is a (N × N ) matrix representing an averaging operator Λ 1 h → Λ 0 h and Λ 1 h → Λ 0 h ; here it corresponds to determining node values from averaging the neighboring cell values. Such separation can be similarly done for M en . D ne is the (N × N ) stiffness matrix with metric-independent coefficients

(D ne ) ml = L χ m (x) dφ l (x) dx dx, (D en ) lm = L dφ l (x) dx χ m (x)dx, (4.3)
with D en = (D ne ) T .

The schemes in matrix-vector formulations

Topological discrete equations. Let us first evaluate the momentum equation (3.13) and (3.15)

χ (1) h , ∂ ∂t u (1) h + χ (1) h , d B (0) h -χ (1) h , ˜ h q (0) h F (0) v h = 0 ∀χ (1) h , (4.4)
χ (1) h , ∂ ∂t v (1) h + χ (1) h , ˜ h q (0) h F (0) u h = 0 ∀χ (1)
h . (4.5)

Using the approximation spaces we get for (4.4) and (4.5)

∂ ∂t u m (t) L χ m (x)χ m (x)dx + B l (t) L dφ l (x) dx χ m (x)dx -q l (t)F v l (t) L φ l (x)χ m (x) dx = 0 ∀m , ∂ ∂t ṽm (t) L χ m (x)χ m (x)dx -q l (t) F u l (t) L φ l (x)χ m (x) dx = 0 ∀m ,
in which we applied the expansions F (0)

u h = F u l (t)φ l (x) and F (0) v h = F v l (t)φ l (x)
for the fluxes, B (0) h = B l (t)φ l (x) for the Bernoulli function, and q (0) h = q l (t)φ l (x) for the potential vorticity that follows from the solution of the vorticity equation (shown further below). In matrix-vector form this reads

∂ ∂t u 1 e + D en B 0 n -M en (q 0 n • F v n ) = 0, ∂ ∂t ṽ1 e + M en (q 0 n • Fu n ) = 0, (4.6)
in which u 1 e = M ee u e is a discrete one-form associated to the vector array u e = {u m (t)|m = 1, ...N }, analogously for ṽ1 e . Note that • denotes the Hadamard (entrywise) product of vectors or matrices. Finally, the evaluation of the integrals in the continuity equation (3.11) yields

χ (1) h , ∂ ∂t h (1) h + χ (1) h , d Fu n = 0 ∀ χ (1) h , (4.7)
∂ ∂t h m (t) L χ m (x)χ m (x) + F u l (t) L dφ l (x) dx χ m (x)dx = 0 ∀m . (4.8)
Using the notion of 1-forms h1 e = M ee he for he = { h m (t)|m = 1, ...N }, we find the matrix-vector equation ∂ ∂t

h1 e + D en Fu n = 0. (4.9)
As each term of both momentum and the continuity equations corresponds to a 1-form, we can consider them as metric-free realization of the corresponding topological equations. Finally, the discrete potential vorticity follows from the (weak) vorticity equation

δ( q (0) h h (1)
h ), ˜ φ

(0) h + δ v (1)
h , d φ (0) h -f dx, ˜ φ (0) = 0 ∀ φ (0) h , (4.10)
q l (t) h m (t) L φ l (x)χ m (x)φ l (x)dx + ṽm (t) L χ m (x) dφ l (x) dx dx -f L 1φ l (x)dx = 0 ∀l .
In matrix-vector form, this reads

M nn (q 0 n • h1 n ) + D ne ṽ1 e -f (∆x e ) T = 0, (4.11)
in which h1 n := A ne h1 e is the height values at cell midpoints averaged to the node values and where A ne is a mean average operator from cell to node values.

Choices for fluxes and Bernoulli function. We will discuss two realizations for the Bernoulli functions and two for the fluxes.

1.) The P0-P1 discretization of the original Hamilton H in Eqn. (3.1) reads in matrix-vector form as

H = 1 2 (u 1 e ) T M en (h 0 n • ũ0 n ) + (ṽ 1 e ) T M en (h 0 n • v 0 n ) E k + 1 2 g( h1 e ) T M en h 0 n Ep . (4.
12)

The associated mass fluxes follow from the projections

φ (0) h , F (0) u h -φ (0) h , h (0)
h u (0) h = 0 ∀φ (0) h , (4.13)
φ (0) h , F (0) v h -φ (0) h , h (0)
h v (0) h = 0 ∀φ (0) h , (4.14)
yielding the discrete mass fluxes in vector array form:

Fu n = h 0 n • ũ0 n and F v n = h 0 n • v 0 n . (4.15)
Moreover, we result in a vector array for the discrete Bernoulli function like

B 0 n = 0.5ũ 0 n • ũ0 n + 0.5v 0 n • v 0 n + gh 0 n . (4.16)
Note that for both fluxes and Bernoulli function, we use only functions in P1.

2.) We have introduced in Section 3.3 an alternative model applying smoothed fluxes and a smoothed Bernoulli function, which provides an alternative structure-preserving approximations of the split 1D slice model of the RSW equations (2.1)-(2.5). In terms of matrix-vector representation of the P0-P1 realization, the matrix-vector form of the smoothed and original model agree except of the following two points:

(i) the discretization of the smoothed Hamiltonian (3.31) reads

H sm = 1 2 (u 1 e ) T h 1 e • M en ũ0 n + (ṽ 1 e ) T h 1 e • M en v 0 n + 1 2 g( h1 e ) T M en h 0 n (4.17)
which leads to the smoothed discrete fluxes

Fu n | sm = A ne h 1 e • ũ0 n and F v n | sm = A ne h 1 e • v 0 n . (4.18)
This provides a set of discrete equations that conserve H sm and the Casimirs of the original system;

(ii) the discretization of the smoothed Bernoulli function (3.32) reads

B 0 n | sm = 0.5A ne (u 1 e • u 1 e ) + 0.5A ne (ṽ 1 e • ṽ1 e ) + gh 0 n . (4.19)
It is compatible with the original Hamiltonian H (3.1) and with the smoothed one (3.31) while for both combinations, the same invariants of the original system of Section 3 are conserved.

As we will see in the next section, using the smoothed fluxes and/or smoothed Bernoulli function instead of the original ones, we can avoid instabilities in the GP1 u -GP0 h and GP0 u -GP1 h schemes (see Remark 4.1) occurring in case when shocks develop.

Metric dependent discrete closure equations. Next, we study the discrete metric closure equations that descend from (3.2), or analogously from (2.13), as matrix-vector representations in form of the above introduced matrices. That is, the evaluation of the nontrivial Galerkin projections for two accuracy levels P0 and P1 provide approximations of the Hodge star operator like

u (0) h = ˜ h u (1) h ⇒ ũ0 n = (M nn ) -1 P ne =: u 1 u 1 e or ũ0 n = (M en ) -1 Id ee =: u 0 u 1 e , v (0)
h = ˜ h v (1) h ⇒ v 0 n = (M nn ) -1 P ne =: u 1 ṽ1 e or v 0 n = (M en ) -1 Id ee =: u 0 ṽ1 e , h (1)
h = ˜ h h (0) h ⇒ h1 e = (P ne ) -1 M nn =: h 1 h 0 n or h1 e = (Id ee ) -1 M en =: h 0 h 0 n ,
where u 1 and h 1 denote the use of P1 as test function in (2.13), and u 0 and h 0 of P0. For more details on the derivation, we refer the reader to [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF]. We use the same Hodge star for both velocity components, but allow this to be different from the Hodge star used for the height fields.

Note the difference between these diagnostic metric equations and the defining diagnostic equation (4.11). The former allow for a free choice of FE spaces that connect the topological equations, the latter has to be in agreement with the topological equations such that potential vorticity and enstrophy are Casimirs of the Poisson equations.

Remark 4.1. The Hodge star operator is discretized via nontrivial Galerkin projections (GP) onto either the piecewise constant (GP0) or piecewise linear (GP1) space. Therefore, we denote the schemes as: (i) GP1 u -GP1 h , (ii) GP0 u -GP0 h , and (iii) GP1 u -GP0 h or GP0 u -GP1 h , respectively, in analogy to the conventional notation for mixed P1-P1 and P1-P0 schemes, cf. [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF] for more details.

Family of split P0-P1 schemes. The introduction of double pairs of compatible FE spaces enriches the choice of potential schemes. For the low order P0-P1 double pairs, we find the following family of split low-order (P0-P1) FE schemes, consisting of one set of topological equations and 4 combinations of metric closure equations: topological momentum eqn.: High accuracy split scheme applies GP1 u -GP1 h ; Medium accuracy split schemes apply GP1 u -GP0 h /GP0 u -GP1 h ; Low accuracy split schemes applies GP0 u -GP0 h , cf. [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF].

∂ ∂t u 1 e + D en B 0 n -M en (q 0 n • F v n ) = 0, ∂ ∂t ṽ1 e + M en (q 0 n • Fu n ) = 0, h 0 n ∈ Λ 0 h ⊂ P 1 D en ----→ Λ 1 h ⊂ P 0 u 1 e , ṽ1 e
GP1 h : M nn h 0 n =P ne h1 e     GP1u: M nn ũ0 n =P ne u 1 e & M nn v 0 n =P ne ṽ1 e GP0 h : M en h 0 n = h1 e     GP0u: M en ũ0 n =u 1 e & M en v 0 n =ṽ 1 e h1 e ∈ Λ 1 h ⊂ P 0 D en ← ----Λ 0 h ⊂ P 1 ũ0 n , v 0 n topological continuity eqn.: ∂ ∂t h1 e + D en Fu n = 0.
Crank Nicolson time discretization. We use the Crank Nicolson (CN) time discretization scheme, which is not fully energy conservation for nonlinear Hamiltonian, such as the one of the RSW equations. Schemes using CN do usually not show a trend in energy for long integration times, but they show oscillations around a mean with values significantly larger than machine precision. Theoretically CN schemes converge at first order with the time step size. For our studies, this behavior is favorable as it allows us to study and distinguish structure preservation related to either the spatial or the temporal discretization. For instance, quantities transported by steady states should be conserved to machine precision also for larger time step sizes, while for non-steady solutions, these quantities are not conserved to machine precision but their conservation properties should converge to machine precision with smaller time step size. We will exploit these properties in the numerical analysis section to thoroughly test the split schemes. An explicit presentation of this time discretization can be found in [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF]. Note that besides this implicit, also explicit time integrators such as Runga-Kutta schemes are usable, but we will not address such time integrators here. Remark 4.2. Note that the split FE method leads to standard matrix-vector formulations in the spatial discretization. As such, they can be implemented straight forwardly to efficient schemes using standard numerical analysis methods.

Properties of the family of schemes

Topological properties. Similarly to Section 3, we distinguish between topological and metric dependent properties for the family of split schemes. We introduce the invariants of the discrete equations in matrix-vector form according to the above representation.

The discrete Hamiltonian (4.12) agrees with the total energy E = E k + E p consisting of kinetic E k and potential energy E p . According to Proposition 3.4, the total energy E is conserved. According to Proposition 3.5, the following Casimirs are conserved in time too:

• mass at elements m e = (∆x) T h1 e or at nodes m n = (∆x) T h 0 n ,

• mass-weighted potential vorticity: PV := (P ne h1 e ) T M nn q0 n ,

• potential enstrophy: PE := (P ne h1 e ) T M nn (q 0 n ) 2 .

In Section 5 we verify the conservation of these invariants numerically.

Remark 4.3. Note that these definitions of invariants descend directly from the matrix-vector representation of the Hamiltonian (3.1) and the corresponding matrix-vector form of the Poisson bracket (3.3), hence the Hamiltonian framework provides the invariants immediately and avoids ambiguity in their definitions.

Remark 4.4. Note the splitting of the schemes' properties associated to topology or metric! The invariants themselves are defined with respect to a choice of metric, i.e. they carry the metric informations. The fact that they are conserved is, however, a purely topological consequence, guaranteed by the Poisson bracket which is defined to be metric independent. This latter property is confirmed by the fact that the resulting prognostic equations are written in terms of (covariant) 1-forms.

Metric properties. Here, we provide two examples of metric dependent properties in order to concretize the abstract discussion in Section 3.2: I.) Consider the dispersion relations of the split schemes for the linearized split RSW model, as derived in Section C in the limit of f = 0, hence for pure gravity waves. This case agrees with the models studied in [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF] in which the derivations of the dispersion relations can be found. In sum, we find for the four choices of metric equations three different dispersion relations. When using the pair of Galerkin projections (GP1 u -GP1 h ) the dispersion relation reads

c d = ω 00 k = ± gH sin(k ∆x 2 ) k ∆x 2 3 cos(k ∆x 2 ) (2 + cos(k∆x)) ;
when using the projections (GP1 u -GP0 h or GP0 u -GP1 h ) they read

c d = ω 01 k = ± gH sin(k ∆x 2 ) k ∆x 2 3 (2 + cos(k∆x)) 1 2
; and finally, when using the projections (GP0 u -GP0 h ), they are This illustrates the metric dependency of the system by the metric closure equations in two ways. First, the dispersion relations themselves changes with different choices of metric equations. And second, each dispersion relation itself is metric dependent, clearly visible by its explicit dependency gH, ω 00 in green for high accuracy GP1 u -GP1 h , ω 01 in blue for medium accuracy GP1 u -GP0 h and GP0 u -GP1 h , and ω 11 in red for low accuracy GP0 u -GP0 h . on the mesh size ∆x. Hence, the metric closure equations impose metric dependent properties on the system of equations.

c d = ω 11 k = ± gH tan(k ∆x 2 ) k ∆x 2 .
II.) Consider the convergence behavior as a second example. As we will show later in the numerical analysis section, the order of the finite element spaces determine the solutions' convergence rates; piecewise constant functions show 1 st , and piecewise linear ones 2 nd order convergence. Convergence is relative to the reduction in the size of the mesh cells ∆x, hence clearly a metric dependent property.

Whereas the concrete realizations of the FE spaces is metric dependent in the sense that, with different FE spaces, properties depending on metric might change, the mutual relations between the compatible FE pairs -determined by the projection operator and the exterior derivative -does not change metric dependent properties of the system.

Numerical analysis

For the numerical analysis of the split schemes we apply three different test cases to study (i) convergence of the schemes, (ii) the long term behavior of the conserved quantities, and (iii) the representation of nonlinear effects.

Suit of test cases. The choice of the corresponding three test cases (TC) is motivated as follows.

For TC 1, we use steady state solutions to study the convergence properties of the split schemes because all deviations of numerical from steady state solutions are caused by numerical errors. However, this is not a demanding problem to study conservation properties because for steady state solutions the quantities of interest (QOI) are usually well preserved given the fact that the solutions do not evolve in time. Therefore, we study the conservation properties of the QOI with TC 2.

In TC 2, we use a flow in (linear) geostrophic balance, i.e. not all terms are fully balanced, in contrast to TC 1. The terms are linearly balanced while nonlinear effects are comparably small. Therefore, this case allows for long integration times while the time dependency of the solutions makes it a more demanding test case for conservation properties than TC 1. We study these quantities for long integration times (about 10 cycles, see explanation in (5.2)) to illustrate the structure-preserving nature of the split schemes.

In TC 3, the initial height distribution is only partly linearly balanced; inertial gravity waves develop and then travel over the domain. Nonlinearity has a larger effect here such that the simulations build shocks during relatively short integration times (at about 0.2 cycles). Therefore, this test case is not suited for long integration times. Instead, we are able to visualize and evaluate the effects of the nonlinear advection term on the development of wave fronts and small scale waves around them.

Flow regime and spatial and temporal resolutions. We perform all simulations in a 1D domain [0, L] with L = 5000 km. We apply periodic boundary conditions and a constant Coriolis parameter of f = 5.3108 days -1 (6.147 • 10 -5 s -1 ) which corresponds to a latitude of 25 • . To account for the test cases' dimensions and enhance readability, we use in the following units of km and days.

We study flows that are dominated by the geostrophic balance, hence in the geostrophic regime in which effects caused by rotation are as important as those by gravity. The corresponding length scale is describes by the Rossby deformation radius L D = √ gH 0 f . Our choice of domain size and spatial resolutions of ∆x = L/N , for N = 32, ..., 2048 elements, for all simulations ensures that L D is well resolved and geostrophic effects equally well represented as gravitational ones. To consider fluids in geostrophic regime in which the flow is dominated by the geostrophic balance, we requires the Rossby number Ro 1. For this study, we distinguish further: (i) quasi-geostrophic regime for Bu ≈ 1 and (ii) incompressible regime for Bu 1 (cf. [START_REF] Cullen | A Mathematical Theory of Large-scale Atmosphere/ocean Flow[END_REF][START_REF] Pedlosky | Geophysical fluid dynamics[END_REF], for instance). We fix H = 75 m which gives a Rossby number of Ro ≈ 0.199. Then, the choice of the background depth H 0 allows us to model flows in (i) quasi-geostrophic regime for H 0 = 750 m and in (ii) incompressible regime for H 0 = 10 km (cf. [START_REF] Bauer | Towards a variational discretization of compressible fluids: the rotating shallow water equations[END_REF] for details).

As mentioned in Section (4.1) we use a Crank Nicolson time integrator, which is an implicit time scheme and is, as such theoretically unconditionally stable. In practice, however, 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. In [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF], we provided a discussion on the upper bound for the different (linear) split schemes. For the nonlinear case, the behavior is similar, but the time step size is a little smaller given the nonlinearities. Here, we will provide only the time step sizes used for the simulations and refer the reader to [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF] for more details on these upper bounds. In general, we choose on fixed time step per test case which will permit us to distinguish between error sources related to the spatial and to the temporal discretizations.

Quantities of interest and error measures. The quantities we are interested in are those presented in Section (4.2), i.e. total energy E(t), mass m e (t) and m c (t), potential vorticity pv(t) and enstrophy pe(t). We expect these quantities to be preserved in time and measure their errors with respect to the relative error measures r f (t) := f (t)-f (0)

f (0)
for any time dependent function f (t). The convergence of the schemes will be evaluated with respect to the following relative L 2 error measures for functions f and velocity components u, v. Let I[f ] = L f (x)dx denote the integral with Gaussian quadrature over the entire domain. Then,

l 2 f = I [(f (T ) -f (0)) 2 ] I [f (0) 2 ] l 2 u = I [(u(T ) -u(0)) 2 ] + I [(v(T ) -v(0)) 2 ] I [u(0) 2 ] + I [v(0) 2 ] (5.1)
in which T denotes the steady state solution at the final time compared to the initial state.

Case 1: Convergence study with steady state solution

We study the convergence in space of the split schemes with respect to a steady state solution of the split 1D slice RSW model. Having such steady state solutions, all errors measured are associated to numerical errors and we study if these errors reduce with smaller grid cell size.

Initialization. We use the steady state solution for the slice 1D RSW equations:

u(x) = 0, v(x) = g4π f L ∆H sin 4πx L , h(x) = H 0 -∆H cos 4πx L , (5.2)
with the choice of parameters as introduced above and project them onto the corresponding FE spaces. This solution is a steady state solution of the 1D RSW equations, cf. Appendix B. We run the simulations for 1 cycle, meaning that the (analytical) wave solutions have traveled once over the entire domain. This simulation time is sufficient to study the spatial influence of the discretization while keeping the errors for the time discretizations low. The initial conditions are indicated in Fig. 5.1 as dashed lines.

For the time step size we use for (i) the quasi geostrophic regime ∆t = 1.6870 • 10 -04 days for GP1 u -GP0 h and GP0 u -GP1 h /GP0 h -GP1 u and ∆t = 1.6870 • 10 -06 days for GP0 u -GP0 h and for (ii) the incompressible regime ∆t = 4.6201•10 -05 days for GP1 u -GP0 h and GP0 u -GP1 h /GP0 h -GP1 u and ∆t = 4.6201 • 10 -07 days for GP0 u -GP0 h .

Results. In Fig. 5.1, we show on the example of the GP0 u -GP0 h scheme the steady state solutions for velocity and height after 1 cycle: left, for the quasi-geostrophic and middle and right, for the incompressible flow regimes. As it can be inferred by the figures, the split scheme preserves the steady state solutions for both flow regimes while the errors decrease with the expected order (shown in Fig. 5.2) for high mesh resolutions. This property is shared by all split schemes. Only for the GP0 u -GP0 h , we see oscillations occurring in the u velocity field at high mesh resolutions (Fig. 5.1, right) for the original Hamiltonian H. However, when using H sm this oscillations can be avoided and convergence can be ensured. In Fig. 5.2 we show the convergence behavior of all schemes for both flow regimes, left for quasi geostrophic and right for incompressible flows. All fields show the expected convergence behavior, namely, P0 fields converge with 1 st order, and P1 fields with 2 nd order towards the analytical solutions. Only for the GP0 u -GP0 h scheme based on the original H, the convergence is not guaranteed for high resolutions (cf. the outlier at N = 512) given the oscillations in the u field (cf. Fig. 5.1, right). The GP0 u -GP0 h scheme based on the smoothed H sm shows, in contrast, the expected convergence behavior. Note that the corresponding relative errors of the different split schemes are very close to each other. Fig. 5.3 shows the dependency of E and pe from the mesh resolution. As we apply a Crank-Nicolson rather then a Poisson time integrators, we do not expect that Energy (or PE) is conserved at machine precision, but we expect a better representation of the invariants for higher resolution. As indicated by the curves, both E and pe converge for high mesh resolutions almost to machine precision in case of incompressible flows (left) and to the order of 10 -12 for the quasi-geostrophic case. Note that these quantities oscillate with the corresponding order of accuracy around their long term means that do not shown significant trends (cf. next test case). Here and henceforth, mass and pv are conserved at machine precision (hence not shown).

Case 2: Conservation behavior for flow in geostrophic balance

We evaluate the structure preserving nature of the split FE schemes. We illustrate that the split 1d slice schemes are capable of conserving the geostrophic balance to a high degree for long integration times while preserving structure, i.e. conservation of E, m, pv and pe. Using the parameter choice introduced above, this test case provides a scenario of geostrophic relevance, with realistic geophysical dimensions. Being linearly balanced, the nonlinear terms trigger small disturbances that propagate over the entire domain. As these nonlinearities are small, we can run the test case for long integration times (here for 10 cycles) to illustrate the good long-term conservation properties of the split schemes.

Initialization. We initialize the surface elevation by

h 0 (x) = H + ∆He -∆w 2π sin π L (x-xc) 2 (5.3)
with x c = L/2 denoting the center of the domain. The velocity fields are initialized by enforcing the (linear) geostrophic balance. To this end, we sample h 0 (x) for each mesh point (and write h 0 n (0) for the corresponding vector) and enforce the geostrophic balance and apply H = 0.75 km to study the (i) quasi-geostrophic, and H = 10 km the (ii) incompressible regimes.

M en v 0 n (0) = g f D en h 0 n (
We use the same time step sizes as for TC1 and run all simulations for 10 cycles except for the GP0 u -GP0 h for only 5 cycles given its the very small time step size.

Results. In Fig. 5.4, we show on the example of the GP1 u -GP1 h scheme the initial conditions (dashed lines) and the solutions after an integration of 10 cycles. The results for the other schemes are very similar, hence not shown. We see that the split schemes conserve the flow in geostrophic balance for long integration times. However, being only linearly balances, gravity waves (visible in the u fields) are emitted from the domain center which travel over the entire domain during the simulation.

For the flow in quasi geostrophic (compressible) regime we realize that the differences in the u velocity field for the initial u(0) = 0 is larger that for the incompressible case. This indicates that structure preservation is more demanding for compressible flows. This can also be seen from Fig. 5.5 in which E is more than two orders of magnitudes better preserved for incompressible than for compressible flows. Similarly, also pe is better preserved in the incompressible case. Note that in none of the cases the time series show a significant positive nor negative trend and that both m and pv are preserved at machine precision. In Fig. (5.6) we illustrate the structure preserving nature of all split schemes for two different mesh resolutions (left: N = 512; right: N = 1024). When comparing the relative errors in the QOI of the different schemes for the same mesh resolutions with each other, we deduce that the errors are more or less independent from the choice of metric equations (compare rows) and that with increasing resolution (compare left and right columns), the errors in the QOI converge to machine precision, similarly to Fig. (5.3). We choose to present here only the more demanding compressible case, while the corresponding QOI in the incompressible case are not shown as they are much better preserved (all better then about 10 -13 for N ≥ 512).

Given its very small CFL number (cf. [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF]), we use for the GP0 u -GP0 h a mesh resolutions of N = 128 (left) and N = 256 (right) with time step ∆t = 1.6870 • 10 -05 days for a simulation time of 5 cycles (shown in the lower row of Fig. (5.6)). Only because of the use of courser resolutions, the conservation of QOI is not as high as for the other cases. Namely, when using N = 1024 elements, E and pe are preserved at the order of about 10 -11 and 10 -10 , respectively, and when using N = 512 elements, they are preserved at the order of about 10 -10 and 10 -09 , respectively, (determined via short term runs). Hence, they behave analogously to the other split schemes. Note that for all cases shown (and studied), the time series do not show a significant trend verifying the structure preserving nature of these low order split FE schemes.

Case 3: Nonlinear advection

By keeping the surface elevation only partly in geostrophic balance while the surface elevation points here upwards, surface waves with half the magnitude of the initial surface elevation are emitted and propagate over the domain. In case of quasi geostrophic, hence compressible flows, these waves develop shock-like structures because of nonlinear effects. We evaluate the schemes' representation of the nonlinear processes by studying the wave fronts and the small scale waves that develop around these fronts in accordance with the discrete dispersion relation.

Initialization. The initialization is equivalent to TC 2 with the exception that the velocity field v is only partly in geostrophic balance, i.e. we use

M en v 0 n (0) = 1 2 • g f D en h 0 n (0) and ũ0 n (0) = 0 (5.5)
for the initial surface elevation (5.3). We use the same spatial and time resolutions as above but integrate only for up to 0.2 cycles, as the nonlinear effects lead to shocks that cannot be resolved properly for much later times.

Results. In Fig. 5.7-5.9 we illustrate that the nonlinearity leads to the building of shocks in case of compressible flows. For the given parameter choice, the wave fronts pile up around the fronts after Namely, for the GP1 u -GP1 h scheme in Fig. 5.7 the highest wave numbers are not propagated (cf. ω 00 in Fig. 4.2), hence the oscillations occur after the fronts; for GP1 u -GP0 h /GP0 u -GP1 h in Fig. 5.9 and 5.10, the highest wave numbers are a little faster then the analytically expected values (cf. curves for ω 10 and ω 01 ), therefore, the oscillations are a little faster then the fronts; finally for GP0 u -GP0 h in Fig. 5.8, the highest wave number tends to infinitely fast wave speed (cf. ω 00 ), hence the oscillations are even much further ahead of the fronts compared to the GP1 u -GP0 h /GP0 u -GP1 h cases. These results agree well with the expected behavior imposed on the schemes by their discrete dispersion relations (and agree well with the discussion in [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF]). Note that, as the development of shocks is only visible in the height and the u velocity field, we only show these plots here and skip the illustration of the v fields.

The results in Fig. 5.9 are obtained by using the H sm from Eqn. (3.31) and in Fig. 5.10 by using the smoothed Bernoulli function B 0 n | sm from Eqn. (4.19). In more detail, whereas the GP1 u -GP1 h and the GP0 u -GP0 h schemes permit to use the original version, the GP0 u -GP1 h is not stable for the original Bernoulli function (4.16) and the GP1 u -GP0 h is not stable for the original flux terms (4.15) here for this test case. Because these are exactly the cases in which the dispersion relations do not support spurious modes, we expect the source of these instabilities to be in the kernel projections in the metric equations between potentially not fully compatible chain complexes of FE spaces. However, using the suggested smoothed fluxes and Bernoulli functions allows us to stabilize the schemes, as illustrated. A more profound study on the stability of the schemes will be subject of future work.

For all split schemes, the relative errors of E and pe here for this demanding, strongly nonlinear test case are at the order of about 10 -04 , or better, on a mesh with N = 512. We use the same time step sizes as above. For higher spatial resolutions, the conservation properties get better, analogously equations and the corresponding split FE discretization method, a pure separation of topological and metric equations can be translated from the continuous to the discrete level. Moreover, this splitting accounts also for the schemes' properties, namely the conservation of invariants descend from the topological structure of the Poisson bracket while the metric properties relate to a concrete choice of FE spaces that, in turn, determine the approximation of the Hamiltonian and its metric-dependent functional derivatives.

(GP 0 u , GP 1 h ) u (0) n (GP 0 u , GP 1
Currently, we further extend the split Hamiltonian idea to other sets of equations of GFD, similarly to [START_REF] Eldred | Variational and Hamiltonian Formulations of Geophysical Fluids using Split Exterior Calculus[END_REF]. Moreover, we will study higher dimensional tensor product schemes as direct extensions of the abstract schemes introduced in Section 3 and we will also address general FE spaces on 2D and 3D meshes. Alongside, alternative approximations of the metric equations will be studied to introduce various structure-preserving primal-primal or primal-dual split schemes of different order of accuracy. Here, the split FE approach provides an ideal framework as the structure preservation is not affected by such modifications in the metric equations.

2d = 2d h, u

2d = 2d u

(1) 2d

(A.1) with twisted potential vorticity function q defined via q h

(2)

2d = d u (1)
2d + f dx ∧ dy with Coriolis parameter f and with B (0) := gh + 1 2 (u

2d ) 2 the Bernoulli function. Note that bottom topography b can be easily included in B (0) via g(h-b), but for simplicity, we do not consider it here. Then, the momentum equations in a 2D coordinate representation read ∂ ∂t u(x, y, t)dx -q(x, y, t)h(x, y, t)v(x, y, t) dx + ∂ ∂x B (0) (x, y, t)dx = 0, (A.2) ∂ ∂t v(x, y, t)dy + q(x, y, t)h(x, y, t)u(x, y, t) dy + ∂ ∂y B (0) (x, y, t)dy = 0, (A.3)

where we used the relations -dx = 2d dy and dy = 2d dx and where x and y denote the spatial, and t the temporal coordinates. To derive the 1D slice model, assume that all coefficient functions depend on x and t only. Then, To find a covariant representation of these coordinate versions, we note that for the Hodge star in 1D 1d = we have dx = 1 and we introduce the 1-forms u (1) = u(x, t)dx, v (1) = v(x, t) dx, and h (1) := h(x, t) dx. This yields ∂ ∂t u (1) -q (0) F (0) v h + d B (0) = 0, ∂ ∂t v (1) + q (0) F (0) u h = 0, ∂ ∂t h (1) + d F (0) u h = 0, (A.8)

∂
for F (0)

v h := h (0) h v (0)
h and F (0)

u h := h (0)
h u (0) . As usual for the splitting, the metric equations close the set of equations. Here they read u (0) = u (1) , v (1) = v (0) , h (1) = h (0) . (A.9)

B Standard and Hamiltonian form of the 1D slice RSW model

For the reader unfamiliar with the covariant or split form of the RSW equations, we introduce the 1D slice model in (I) standard coordinate representation, cf. [START_REF] Staniforth | Analysis of a mixed finite-element pair proposed for an atmospheric dynamical core[END_REF], and (II) in standard Hamiltonian form. I.) For velocity components u(x, t), v(x, t), and height h(x, t), the 1D slice model in coordinate representation reads with the definition of q as given above. The latter formulations can equivalently be written in an almost Poisson bracket, as done for the split schemes in the main text.

Remark B.2. For a fixed orientation and by substituting the metric into the topological equations, it is straightforward to show that the split equations, derived in Section A agree with the standard equations (B.4).

C Linearization of the split RSW equations

We derive the linearization of the split RSW equations from the Hamiltonian, linearized around background height H and u

(1)

h = v (1)
h = 0, i.e.

H[u

h , v

h , h

h ] = u

h , H d u

(0) h ∀ χ (1)
h .

(C.1)

Proposition 3 . 2 .

 32 The equations in Definition 2.3 imply the Poisson bracket formulation in Definition 3.1.

h

  as test function giving δC = 0. Hence (3.23) vanishes and PV is conserved. The same argument holds for mass with C = h (1)

  h (x). These functional derivatives follow from Lemma 3.3 which also holds for the smoothed fluxes F (0) u h and F (0) v h and the smoothed Bernoulli function B (0) h | sm . Similarly to the proof of Proposition 3.2, the smoothed version of the equations of motion (in weak form) follows from these functional derivatives and the Poisson bracket 3.3, i.e.

Figure 4 . 1 :

 41 Figure 4.1: Family of split low-order (P0-P1) schemes consisting of one set of topological equations with original or smoothed fluxes and Bernoulli function as well as four realizations of metric equations:High accuracy split scheme applies GP1 u -GP1 h ; Medium accuracy split schemes apply GP1 u -GP0 h /GP0 u -GP1 h ; Low accuracy split schemes applies GP0 u -GP0 h , cf.[START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF].

Figure 4 .

 4 Figure 4.2 illustrates the three different dispersion relations for the four split schemes.This illustrates the metric dependency of the system by the metric closure equations in two ways. First, the dispersion relations themselves changes with different choices of metric equations. And second, each dispersion relation itself is metric dependent, clearly visible by its explicit dependency

Figure 4 . 2 :

 42 Figure 4.2: Dispersion relations: analytic (black) for c =√ gH, ω 00 in green for high accuracy GP1 u -GP1 h , ω 01 in blue for medium accuracy GP1 u -GP0 h and GP0 u -GP1 h , and ω 11 in red for low accuracy GP0 u -GP0 h .

Figure 5 . 1 :

 51 Figure 5.1: Steady velocity and height fields after 1 cycle for the GP0 u -GP0 h scheme; left in the quasi geostrophic regime, middle and right in the incompressible regime for a mesh with 512 elements. Initial fields are shown in dashed lines. In the middle panel, H sm is used which avoids oscillations occurring for the original H at highly resolved meshes (right).

Figure 5 . 2 :

 52 Figure 5.2: Convergence for the steady state test case in quasi geostrophic regime (left) and incompressible regime (right) for a simulation time of 1 cycle.

Figure 5 . 3 :

 53 Figure 5.3: Relative errors in energy (E) and potential enstrophy (pe) errors for the steady state test case in quasi geostrophic regime (left) and incompressible regime (right) for a simulation time of 1 cycle.

4 )

 4 For the linearized equations in (C.1), these are steady state solutions with∂ ∂t h1 e = ∂ ∂t u 1 e = ∂∂t ṽ1 e = 0. We use the same parameters as above, i.e. L = 5000 km, g = 7.32 • 10 7 km days -2 , f = 5.3108 days-1

Figure 5 . 4 :

 54 Figure 5.4: Solutions after 10 cycles of the GP1 u -GP1 h scheme for (i) quasi geostrophic (left) and (ii) incompressible regime (right) on a mesh with N e = 512 elements.

11 Figure 5 . 5 :

 1155 Figure 5.5: Time series of the quantities of interest for 10 cycles of the GP1 u -GP1 h scheme for (i) quasi geostrophic (left) and (ii) incompressible regime (right) on a mesh with N e = 512 elements.

9 Figure 5 . 6 :

 956 Figure 5.6: Time series of the quantities of interest for 10 cycles (resp. 5 cycles) of the GP1 u-GP1 h , GP1 u -GP0 h /GP0 u -GP1 h (resp. GP0 u -GP0 h )schemes for quasi geostrophic flows on mesh with N e = 512 (left) and N e = 1028 elements (right) (resp. N e = 128 (left) and N e = 256 (right)).

Figure 5 . 10 :

 510 Figure 5.10: Fields with oscillations at the wave fronts for the GP0 u -GP1 h scheme with B 0 n | sm in the quasi geostrophic regime on a mesh with N e = 512 elements and after a simulation time of t = 0.225T .

7 Acknowledgements

 7 WB has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 657016. JB acknowledges support by the DFG excellence cluster CliSAP (EXC177).A Derivation of the split 1D y-independent slice RSW modelThe split covariant RSW equations in 2D read

1 . 3 )

 13 The derivatives of the steady state solution (independent of time t) of Eqns. (5Putting them into (B.1) shows that (5.2) is in fact a steady state solutions of the 1D slice model. II.) Alternatively to the advection form (B.1), a curl-form of the RSW equations exists (cf.[START_REF] Cotter | Compatible finite elements for numerical weather prediction[END_REF]): h(x, t), velocity u(x, t), v(x, t), Coriolis parameter f , and gravity g. Based on this curl-form, a standard Hamiltonian form can be introduces as follows. For a Hamiltonian motion for the 1D slice model follow from the structure

  h , with metric equations(3.2). This give the functional derivatives: h ), zero else.The dynamics are given by the almost Poisson bracket (3.3) for the linearized q

  Figure 5.9: Fields with oscillations at the wave fronts for the GP1 u -GP0 h scheme with H sm F in the quasi geostrophic regime on a mesh with N e = 512 elements and after a simulation time of t = 0.225T .

					1.4
		1
			h(x, 0)		1.3								0	u(x, 0)
			h (1) e (GP 1 u , GP 0 h )		1.2								u (1) e (GP 1 u , GP 0 h )
		0.95	h (0) n (GP 1 u , GP 0 h )	km	1.1								u (0) n (GP 1 u , GP 0 h )
					1
					0.9
		0.9			0.8								-500	1500
	km	0.85			0.7	0	1000	2000	km	3000	4000	5000	km/day	500 1000
													-1000	km/day	0
													-500
		0.8
													-1000
													-1500
													0	1000	2000	3000	4000	5000
		0.75											km
													-1500
		800	1000	1200			1400			1600		1800	800	1000	1200	1400	1600	1800
					km							km
					1.4
		1
			h(x, 0)		1.3								0	u(x, 0)
			h (1) e (GP 0 u , GP 1 h )		1.2								u (1) e
		0.95	h (0) n (GP 0 u , GP 1 h )	km	1.1
					1
					0.9
		0.9			0.8								-500
	km	0.85			0.7	0	1000	2000	km	3000	4000	5000	km/day
													-1000
		0.8
		0.75
													-1500
		800	1000	1200			1400			1600		1800	800	1000	1200	1400	1600	1800
					km							km

INRIA Rennes, France. werner.bauer.email@gmail.com

Imperial College London, United Kingdom.

CEN/Department of Mathematics, Universität Hamburg, Germany.

to Fig. 5.3. Again, mass and potential vorticity are preserved at machine precision for all schemes. h (1) e (GP 1 u , GP 1 h ) h (0) n (GP 1 u , GP Fields with oscillations at the wave fronts for the GP0 u -GP0 h scheme in the quasi geostrophic regime on a mesh with N e = 512 elements and after a simulation time of t = 0.225T .

Conclusions

We introduced a split one-dimensional y-independent slice model of the nonlinear rotating shallow water equations and a corresponding split Hamiltonian form. Using the idea of the split FE method of [START_REF] Bauer | A structure-preserving split finite element discretization of the split wave equations[END_REF], we first introduced a general split FE discretization of arbitrary order for which we showed conservation of energy, mass, potential vorticity and enstrophy. As concrete realization, we derived a family of low order split P0-P1 schemes and we studied their properties related to topology and metric for test cases of geophysical relevance. Analytical calculations and simulations showed that the structure preservation is carried by the topological equations independently from the metric equations. The latter, in turn, determine the schemes accuracy, convergence, stability, and discrete dispersion properties. In this vein, we could control occurring instabilities for a strongly nonlinear test case by using an alternative, smoothed Hamiltonian. This modification did not negatively impacting on the conservation properties.

Besides these locally conservative discretizations that can easily be extended to higher dimensional tensor product schemes, the major outcome of the manuscript is the following: by using the split


Farid Medjedoub
email: farid.medjedoub-renexter@renault.com

Gilles Dour

Sabine Le Roux
email: leroux@enstimac.fr

Pascal Lamesle
email: pascal.lamesle@enstimac.fr

Mehdi Salem
email: mehdi.salem@enstimac.fr

Patrick Hairy
email: hairy@ctif.com

Farhad Rézaï-Aria

Experimental conditions and environment effects on thermal fatigue damage accumulation and life of die-casting steel X38CrMoV5 (AISI H11)

Keywords: Thermal Fatigue, TF, crack and initiation, chromium steels, fatigue-oxidation interaction fatigue

Thermal Fatigue (TF) is a life-limiting factor in die-casting dies. The effect of the maximum temperature of thermal cycle, T max , and the heating (or heat flux density) of thermal fatigue (X38CrMoV5 AISI H11) at 47 HRC is investigated. Two microscopic and macroscopic heat-checking cells are identified. The microscopic heat checking is limited to the double-layer oxide scales formed on the external surface of TF specimens. An engineering Paris-type law using K max , can describe TF surface crack propagation of all tests examined. TF life based on σ max is used to rationalise TF micro-and macroscopic cracking.

Introduction

Thermal Fatigue (TF) causes gradual deterioration and eventual cracking of materials by alternate heating and cooling during which free thermal expansion/contraction is partially or completely constrained [START_REF] Rézaï-Aria | Thermal fatigue of hot work tool steels[END_REF]. Hot-forming tools work under transient thermal cycling and so are prone to TF crack initiation and propagation (as in rolling [START_REF] Da Silva | Thermal cracking of multicomponent white cast iron[END_REF], forging, die casting, etc.). TF is a life-limiting factor in die-casting dies, in particular high-pressure die casting. TF damage may be activated and/or coupled with environment effects such as oxidation and corrosion (see Figure 1). Surface heat-checking that forms under multiaxial loading is the typical damage network observed in die-casting and hot-working tools. The temperature distribution in a die changes from one region to another [START_REF] Medjedoub | Endommagement par faïençage des moules de fonderie sous pression en fatigue thermique: origines, mécanismes et approaches[END_REF]. The crack initiation and propagation is controlled by the actual stress and thermal conditions, as well as by environmental effects. It is therefore important to know about the industrial thermal cycle characteristics. This contribution deals with the investigations of heat-checking under various thermal fatigue conditions. The heat-checking mechanisms and oxidation effects are discussed. Thermal fatigue crack propagation and life are reported. 2 Experiments and materials

Materials

A double-tempered martensitic steel X38CrMoV5 (AISI H11) with 47HRC (mostly) was investigated. Forged bars were purchased from Aubert and Duval France and double heat-treated using Balzers (Balitherme structure). TF specimens were then machined from the heat-treated bars.

Thermal fatigue reference cycle and test conditions

Recently a special well-designed heat flux gauge (sensor) using a Sapphir fibre and several well-placed thermocouples with small diameters was developed to measure the heat flux density transferred during casting from the surface to the die (Dour et al., 2005a;2006) (see Figure 2). The temperature at the surface of the die is calculated via normalisation and a numerical inverse method by proceeding with a careful calibration of the pyrometer connected to the optical Sapphir fibre for determination of the emissivity of the molten Al-alloy as well as the entire data acquisition chain response, the position of thermocouples in the gauge and also a careful placement of the sensor in the die considering the critical position in a die. As an example, the variation of the calculated temperature via the sensor and estimated thermoelastic stress versus time is reported (see Figure 2).

Details of TF and an in-house lab view data acquisition system using induction heating (25 kW, 100 to 400 kHz) are given in [START_REF] Jean | An investigation on heat checking of hot work tool steels[END_REF] and [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. Tubular specimens (10 mm wall thickness) with no external air forcing and an internal continuous water-cooling are employed (see Figure 3). The external surface is polished. Based on several heat flux measurements on HPDC dies, a typical 'reference thermal fatigue cycle' is used. A dummy specimen was first carefully instrumented by at least five K-type thermocouples spot-welded on the external surface (0.1 mm wire's diameter), and eight thermocouples brazed in drilled holds along the longitudinal direction, z-axis (0.2 mm wires diameter). The measured temperature-time cycles, T-t, were then used to optimise numerical thermal analysis down under ABAQUS. The minimum temperature at the external surface (T min ) is set at 100°C. Two series of TF experiments are performed with: 1) four heating periods (1.2, 2, 3.7, and 6.5 s) to achieve T max = 650°C, 2) five T max (500°C, 550°C, 600°C, 650°C and 685°C) with the same heating period (about 1 s). To achieve such conditions, a constant energy is put on the external surface of specimen. The resulting heat flux density is then estimated by an inverse method well described in Dour et al. (2005a) (see Dour et al., 2006 in this journal). A constant prescribed nominal 'induction power' is therefore input on the external surface of TF specimens for heating. The temperature is continuously recorded via a thermocouple-K spot-welded on the middle part of specimens. This part of specimens is used for monitoring the experiments. The temperature-time cycle measured is used for eventual minor readjustments of induction power while testing. The opposite part of this region is used for optical and SEM observations. Tests were regularly interrupted to examine systematically the formation and evolution of both micro-and macroscopic heat checking patterns (or cell forming) on prescribed regions. The morphology characteristics of cells (density, mean diameter, etc.) were calculated by image analysis of SEM-constructed cartographies.

Results and discussions

Thermal and mechanical analysis

Thermoelastoplastic strains and stresses were calculated by ABAQUS using axisymmetric DCAX4 and CAX4 elements for thermal and mechanical numerical analyses, respectively [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. As boundary conditions in thermal analysis, several heat flux densities Q and a constant heat transfer coefficient were imposed to external (subdivided in several strips) and internal surface nodes, respectively. The heat flux densities were estimated by 1D analytical calculations [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]Dour et al., 2005b). The measured (see Figure 4a) and calculated temperature-time cycles were compared. Then Q was adjusted [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF] for each region to minimise the difference between the estimated and the measured temperatures. A difference less than 7% was achieved for all TF experiments. For FEM thermomechanical numerical analysis, tensile curves measured at prescribed temperatures ranging from RT to 700°C were used according to ABAQUS software requirements for isotropic hardening constitutive calculations. As an example, a mechanical total strain-temperature loop calculated for a selected element located in the centre of the TF specimen is reported in Figure 4(b). As can be observed, this element is put under an extended compressive strain in first heating. Then the element is put under tension during cooling. The maximum stress σ max is achieved at T min . The σ min during the first heating is achieved at an intermediate temperature between T min and T max . The first cycle produces maximum plastic strain amplitude defined as the difference between strain at point A and point D in the first hysteresis loop. This is a typical stress-strain behaviour in the thermal fatigue of metallic materials. During the second heating and cooling cycle, the plastic strain becomes very much less. In addition, contrary to the first cycle, σ min is achieved at T max . This is one of the reasons why several fully stop/restarting operations in hot forming (for example, for changing or repairing of tool) do more damage to a tool experiencing transient temperature conditions compared to a non-stop working operating condition. The thermo-elasto-viscoplatic analyses show that σ min and σ max can evolve [START_REF] Ahmer | Cyclic behaviour simulation of X38CrMoV5 (AISI H11) tempered martensitic hot work tool steel[END_REF]. However, as the plastic strain of concern is small, an elasto-plastic calculation can be relevant to the damage analysis and life modelling [START_REF] Daffos | Endommagement des outillages de forgeage traités par nitruration: Étude et Modélisation[END_REF][START_REF] Oudin | Thermomechanical fatigue behaviour and life assessment of hot work tool steels[END_REF].

Damage mechanisms: oxidation, micro-and macroscopic TF cracking

While the oxide morphology depends upon T max , the oxide scale always consists of tow layers, one rich in Cr beneath the steel and one poor in Cr in direct contact with air that has quasi-identical thickness (see Figure 5). The TF oxidation kinetic seems to obey a classical parabolic law (doxide = K.√t) [START_REF] Medjedoub | Modeling of Cyclic Oxidation at Variable Temperature, Application to X38CrMoV5 Tool Steel[END_REF][START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. The same kinetics is also observed in isothermal LCF [START_REF] Daffos | Endommagement des outillages de forgeage traités par nitruration: Étude et Modélisation[END_REF]) and non-isothermal (thermomechanical) fatigue [START_REF] Oudin | Thermomechanical fatigue behaviour and life assessment of hot work tool steels[END_REF]. Both oxide-scales are cracked by 'heat-checking'. Such damage mechanism is named 'microscopic heat-checking' considering the depthsof short cracks that are comparable to the thickness of both oxide scales (see Figure 5). Macro cracks propagate inward base steel from the short crack tips forming the micro heat-checking cells. By the coalescing of isolated cracks and recracking of already existent cracks, a network of interconnected macro cracks, named 'macroscopic-cell' or 'macroscopic heat-checking' form (see Figure 5). Figure 6 shows that the heat-checking (density) is T max and/or ∆T (= T max -100°C) independent while it is heating-rate (or heat flux density)-dependent when ∆T is constant (T max = 650°C). et al., 2006;[START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. The cartography of selected regions is constructed by assembling nine SEM images taken from a surface (area 0.8*1 mm*mm). The procedure has several steps: the median filtering to reduce the noise in the grey scale image, the contrast enhancement, thresholding routine to detect darken micro cracks and binarisation of the image to finally obtain a network of closed skeletonised cells that is convenient for relevant quantification and image analysis [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF].

Figure 7 presents the variation of heat-checking density ρ as a function of the number of TF cycles for all experiments. The heat-checking network forms very early and revolves drastically according to a sigmoïdal trend to achieve a quasi-asymptotic limit that remains unchanged with further TF cycles [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. This asymptotic limit is called the stabilised, or saturated heat checking cell density (ρstab). The number of cycles required to achieve ρstab is named Nstab. The variation of ρstab as a function of φ presents a quasi-linear trend (see Figure 8). By increasing the heat flux density, ρstab increases. It should be emphasised that ρstab might presumably achieve a certain limit for much higher Φ. However, owing to limitations of our TF rig facility, such heat-checking morphology characteristics were not assessed. The number of cycles for the micro heat-checking initiation (Na) is an arithmetic value between the cycle when the heat-checking is first detected (N2) and the last previous TF test interruption (N1): Na = (N1 + N2)/2. The precision of Na depends principally on intervals of test interruptions. In major experiments, tests were interrupted at each 500 TF cycle prior to cells formation thus giving a relatively high accuracy to Na.

Thermal fatigue crack propagation

Several destructive experiments are required to establish an in-depth crack propagation curves [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. We have measured the mean and maximum crack depths. When cracks are short (less than few millimetres), the mean and/or maximum crack length could reasonably be determined with acceptable statistical estimations. As only few cracks propagate inward, sometimes shorter crack depths could be obtained with a higher number of TF cycles, which is irrelevant. Thermal stresses decrease from the external surface and therefore long cracks tend to stop, making it difficult to achieve the reliable in-depth crack versus the number of cycle curves.

The maximum and the mean surface crack lengths were determined on a circumferential strip of 20 mm width in the centre of specimens (see Figure 9). No macroscopic cracking could be identified with T max lower then 600°C. With increasing T max , the crack propagation rate increases. Under a prescribed TF condition (T max and T min ), the crack propagation curve is heating period-dependent. In fact, by reducing the heating period, the heat flux density input to the external surface increases, leading to enhance the thermal stresses and consequently the crack propagation rate. Figure 10 presents the TF propagation rate as a function of calculating the maximum Stress Intensity Factor (SIF), K max = σ θθmax √π.a mean . This is an 'engineering' definition of SIF, since it does not consider multiaxial loading and also a mean presents an 'equivalent' crack length. This Paris-type relationship seems to better rationalise the TF crack propagation than mechanical strain intensity factor defined as: ∆K max = ∆ε m .√π.a mean [START_REF] Rézaï-Aria | [END_REF][START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF][START_REF] Bergström | High temperature fatigue of tool steels[END_REF]. In general, TF resistance rates of tool materials are compared based on a maximum or a mean crack depth achieved under a prescribed number of TF cycles. The maximum temperature is not the only relevant criterion for the TF life prediction, since depending upon thermal gradient and thermal stains and stresses, different life cycles can be achieved. We have defined several fatigue life cycles corresponding to several prescribed surface crack lengths (for example, 0.1 mm, 0.3 mm, etc.) [START_REF] Medjedoub | Détermination des paramètres influant sur le phénomène d'endommagement par fatigue thermique des moules en fonderie sous pression d'aluminium[END_REF]. Then each TF life was plotted as a function of T max or mechanical parameters such ∆ε m , ∆σ, σ θθmax , etc. Figure 11 presents TF life as a function of the σ θθmax . Such mechanical approaches are more relevant to the comparison between different material and also between TF and isothermal LCF [START_REF] Oudin | Thermomechanical fatigue behaviour and life assessment of hot work tool steels[END_REF] or TMF [START_REF] Le Roux | [END_REF][START_REF] Engler-Pinto | Thermo-mechanical fatigue investigations of SRR99[END_REF] life curves.

Conclusions

The thermal fatigue behaviour of X38CrMoV5 AISI H11 (47HRC) is studied. The effects of the maximum temperature of thermal cycle (from 550°C to 685°C) and the heating period (from about 1.2 to 6.5 s) are assessed.

The micro-and macroscopic heat-checking cells are formed. The microscopic heat-checking mechanism consists of surface oxide-scales cracking. The depth of such cracking is mainly limited to the thickness of oxide scales. The formation of the micro-and macroscopic heat checking is due to the multiaxial mechanical stress and/or strain loadings that are generated by transient thermal fatigue solicitations.

The evolution of the microscopic heat-checking density versus the number of thermal fatigue cycles shows a sigmoïdal trend and achieves a quasi-asymptotic limit. This stabilised density is very much heating period (or heat flux density)-dependent, while it is quasi-independent on maximum temperature of the thermal cycle. The macroscopic heat-checking network is formed from the microscopic heat-checking.

The thermal fatigue strains and stresses are calculated via numerical simulations using ABAQUS. A Paris-type law relating the mean thermal fatigue crack propagation rate to the maximum stress intensity factor, K max , could rationalise the thermal fatigue crack propagation resistance of the steel examined under different thermal fatigue test conditions.

It is found that σ max is a relevant criterion for rationalising both the microscopic heat-checking initiation life and the surface macroscopic cracking life.

Figure 1

 1 Figure 1 A cross-sectioned view of typical cracking of high-pressure Al die-casting die surfaces hat are damaged by the surface heat-checking

Figure 2

 2 Figure 2 Heat flux density gauge placed in HPDCs (a) and calculated temperature and thermoelastic stress versus time (b)

Figure 3

 3 Figure 3 Thermal fatigue tubular specimen

Figure 4

 4 Figure 4 Measured surface temperature time cycle: (a) at various heating period (or heat flux density) T max = 650°C and (b) a typical calculated stress-mechanical strain hysteresis loop by numerical simulations

Figure 5 Figure 6

 56 Figure5The role of the oxidation on the formation of the micro-and macroscopic heat-checking thermal fatigue cracking

Figure 7

 7 Figure 7 Evolution of the microscopic heat-checking cell's density as a function of the number of TF cycles

Figure 8

 8 Figure 8 Variations of stabilised microscopic heat-checking density, ρstab as a function of the heat flux density Φ for all TF conditions examined

Figure 9

 9 Figure 9 Evolution of the surface mean crack length as a function cycle for different thermal fatigue test conditions

Figure 10

 10 Figure 10 Thermal fatigue crack propagation rate as a function of the maximum stress intensity factor

Figure 11

 11 Figure 11 Thermal fatigue life (N cycles) as a function of σ θθmax

Acknowledgements

The authors from Ecole des Mines d'Albi-Carmaux would like to acknowledge the French Casting Centre (CTIF) for partially supporting these activities. The authors also acknowledge Ing. G. Bocquet and Dr. A. Hamasaiid from the Ecole des Mines d'Albi-Carmaux for their valuable contributions to the heat transfer measurements and analyses on actual dies.