Analysis of dissipation operators that damp spurious modes while maintaining discrete approximate geostrophic equilibriums for the B-grid staggered scheme on triangular meshes

We consider the lowest order - so called B-grid - discretization on triangular cells which is common in the ocean or atmosphere simulation science, with pressure field at the vertices of the triangles and the velocity field at the cell centers. It is known that this discretization presents spurious modes due to a too large number of velocities unknowns and that these modes should be damped; since preservation of approximate geostrophic equilibrium (GE) may be of primary importance in the simulations, we propose and analyze damping operators that, when applied to the linearized constant coefficient configuration, exactly preserve discrete versions of the GEs. We also analyze their effects in more general configurations and their possible coupling with more standard damping operators. These strategies cover triangular variants of the Low Froude (LF) scheme [1], the Apparent Topography (AT) scheme [2, 3], as well as a new modification of it. The LF scheme has no numerical diffusion terms in the pressure equation and damps the normal velocity jumps between adjacent cells in the momentum equation. In the AT scheme, numerical diffusion is added in the pressure equation in such a way that it does not impact the discrete GE. However, since the AT scheme alone cannot be proved to be damping through energy estimates and does not preserve the space which is orthogonal to the discrete GE, we also study what we call the Modified Apparent Topography (MAT) that satisfies these important properties. We also suggest semi-implicit time discretizations that preserve these properties. Extensions to the full non-linear equations are also provided. Numerical experiments illustrate the performance of the new schemes.


Introduction
Large-scale geophysical flows are those in which the earth's rotation has a significant impact on the dynamics [4]; in such configurations, oceanic and atmospheric circulations are often perturbations of the so-called geostrophic equilibrium (GE) which results from the balance between the pressure gradient and the Coriolis force [4,5].When such circulations are investigated by numerical approximations of the corresponding physical model, it is therefore essential that the numerical schemes do not completely damp out this approximate equilibrium.Although the GE is not a perfect steady-state in real dynamics because of forcing and non-linearities, there is some interest in studying whether numerical schemes preserve it in the linear case, for the reasons exposed below.
We begin with the investigation of the dimensionless shallow water system of equations on the rotating frame, which is given by In System (1) unknowns h and ū respectively denote the water or layer depth and the velocity of the fluid and function b(x) denotes the topography and is a given function.Dimensionless numbers St, F r and Ro respectively stand for the Strouhal, the Froude and the Rossby numbers.As mentioned above, we are interested in large scale flows, so we will focus on the case in which the Froude and Rossby numbers are small.In particularly, we shall assume that with M is a small parameter (typical values lead to M ∼ 10 −2 ).For a Strouhal number of order O( 1 M ), i.e. for small time scale, the solution of system (1) with flat topography is h = h 0 + M r + O(M 2 ) and ū = u + O(M ), where h 0 is a space-time constant and where, up to multiplicative constants, (r, u) satisfies the linear wave equation with Coriolis source term where u = (u, v) T , and u ⊥ = (−v, u) T .The parameters a ⋆ and ω are of order one, respectively related to the wave velocity and the rotating velocity.In case of constant coefficients, System (2) admits as stationary state the GE which is given by with incompressibility of the velocity field As a consequence, a good reason for analyzing schemes that can preserve a discrete version of the GE in the linearized constant coefficient case, is that this is the minimum requirement that a scheme should have in order to get accurate simulations in more general situations around the GE.We shall then mention how the schemes we propose may be extended to more general cases.
There are few works related to numerical schemes that preserve a discrete version of geostrophic equilibriums.Since Godunov type finite volume collocated schemes (all unknowns are located at the centers of the mesh cells and fluxes are computed using a Riemann solver) have a long and successful history for the numerical approximation of systems of conservation laws such as (1) with ω = 0, it is quite natural to first try to extend them to the treatment of the Coriolis source term (ω ̸ = 0).Some success in this direction on Cartesian grids was previously obtained in [6], but we show in a discussion that is postponed to the Appendix that the direct extension of these schemes to triangular meshes seems to lead to a dead end (we recall that triangular grids may be preferred to Cartesian grids for example because they are much more flexible for accurate discretizations of complex detailed geometries and for local mesh refinement).Instead, we investigate in the present contribution the lowest-order (k = 0) staggered arrangement belonging to the well-studied family of schemes with completely discontinuous P disc k velocities and continuous Lagrange P k+1 pressures.As is described in [7], this configuration is very natural for the preservation of a discrete version of (3), since discrete pressure gradients belong to the velocity space.An additional useful property in the lowest order case is that discretely divergence free velocity fields can be characterized as the curls of P 1 Lagrange functions defined by their values at the vertices of the triangular cells.Finally, our choice of the P 0 − P 1 configuration is supported by [8] which stresses that, among a large study of finite element pairs, this is an appealing choice (see the Conclusion of [8] for details).
This, as well as the fact that such a staggered strategy (cell-centered velocities and node based pressures) is proposed and analyzed in the finite volume coastal ocean model (FVCOM) [9]) and actually used for large-scale ocean modeling [10,11], leads us to the main purpose of this work, which is to propose a family of staggered explicit schemes with damping terms on triangular cells in which pressures and velocities that verify discrete versions of (3)-( 4) are not affected by the numerical diffusion of the scheme.This numerical diffusion is necessary for at least two reasons: it maintains the scheme stability when an explicit (or semi-implicit) discretization in time is used, and it helps damp spurious numerical modes that are associated with this staggered configuration, as analyzed for example in [12] (together with many other finite element pairs) and in [13].In the latter, discrete velocity Laplacians are used not only to damp the spurious modes, but also to damp variance on small scales.Full velocity Laplacians have the drawback of destroying GEs in circumstances in which the flow takes place around this approximate equilibrium, but may be needed in other circumstances, for example for turbulent forced-dissipative flows.In the present work, the schemes we propose preserve discrete GEs while damping the spurious modes; in order to perform this, they do not use full velocity Laplacians, but only the normal direction of the velocity jumps at the cell interfaces.We indicate in Section 5 of the article how this could be coupled to diffusion of the tangential velocity jumps to recover standard full Laplacians when they are needed.A by-product of the analysis is that we are able to characterize the spurious modes on general triangular grids: they are velocities which have vanishing discrete curl and divergence, while not being zero.Additionally, the family of schemes we propose also verifies points 2 to 6 mentioned in [14]: mass conservation, curl-free pressure gradient, energy conserving pressure terms, energy conserving Coriolis term and steady geostrophic modes.An additional feature which, to our knowledge, has not been discussed previously in the literature at the discrete level, is the rigorous definition of a splitting between the vector subspace of discrete GEs and an orthogonal subspace containing their perturbations.In the linearized case, some of the schemes we present not only preserve discrete GEs, but also the orthogonal subspace: whenever the initial condition is orthogonal to the space of GEs, then so is the solution at any time; this ensures that perturbations around GEs do not exchange energy with the space of GEs.Finally, we also show the stability of two of the proposed schemes.
This work is organized as follows.We propose some semi-discrete (continuous in time) staggered schemes in Section 2: • The Low Froude staggered scheme (LF), • The Apparent Topography staggered scheme (AT), • The Modified Apparent Topography scheme (MAT).
Particularly, we show that the discrete differential operators that appear in the scheme definition satisfy mimetic properties that lead to the desirable properties mentioned above.The mathematical properties mentioned in Sections 2.4 to 2.6 mainly concern the schemes in the linearized, constant coefficient case; this is interesting per se from a mathematical point of view and gives hints on what may happen in more general configurations.Stress is then added on the damping properties of the schemes we propose.We then take into account in Section 3 the time discretization to introduce appropriate fully discrete schemes that still possess the same properties as the semi-discrete schemes.Some numerical test cases are shown in Section 4 to confirm the analysis led in the theoretical part and we discuss an extension to the case of non constant Coriolis parameter.Then, further extensions to non-flat topography and non-linear configurations are discussed and tested in Section 5, in which it is also shown that combining the diffusive operators we propose with some additional velocity diffusion allows to recover usual full velocity Laplacian operators.
Concluding remarks complete the study in Section 6, and in Appendix A, we briefly present the problem encountered with the collocated Godunov scheme on triangular meshes when trying to check whether the scheme steadystates verify discrete versions of the GE.

Definition of the discrete operators and the semi-discrete staggered scheme
We first define the discrete version of the gradient and divergence operators.As mentioned above, the discrete gradient of a scalar field will be defined over the cells, from values of the scalar field defined at the vertices.The discrete divergence operator will then be defined such that a discrete integration by parts holds; in this way we shall be able to prove stability of the scheme by energy estimates.We shall denote by T i a generic triangular cell used for the discretization of the computational domain which is supposed to be rectangular, periodic with no internal hole, for the sake of simplicity.Let A ij be the common edge of the neighboring cells T i and T j , and n ij be the unit normal vector to A ij pointing from T i to T j .Moreover, we also denote the area of the cell T i by |T i | and the length of the edge A ij by |A ij |.
Let (r k ) k∈ [1,Nr] be a discrete scalar field defined by its values at the vertices of the mesh, where N r is the number of vertices.Let us denote by r h the globally continuous function that is piecewise P 1 on each cell and defined by the values (r k ) k∈ [1,Nr] at the vertices of the mesh.We can define the discrete gradient on the (primal) cells by using the following formula where for any edge A ij , we have denoted by N ij and S ij its extremities (see Figure 1(b)).It is easily checked that ( 5) is the gradient of the P 1 Lagrange function r h , so that this definition is quite standard.Let (u i ) i∈[1,N ] be a discrete vector field defined by its values on the triangular cells.We shall define its divergence on the so-called barycentric dual mesh constructed as follows.Each vertex is associated with a dual cell obtained by joining the barycenters of the cells which share the vertex to the midpoints of the edges (see Figure 1(a)).Let us denote the area of the dual cell D k by |D k |.Then, we can define the discrete divergence (∇ D h • u h ) on the dual cell by the following (also quite standard) formula where for a triangle T i and a vertex k, l ik is the length of the edge of T i that is opposite to vertex k and n ik is the unit normal vector pointing outside T i on this edge.This edge will be named A ik below.Moreover, we can define the discrete curl of the vector field by On the other hand, we can define the discrete scalar product between q 1 h = (r 1 h , u 1 h ) and q 2 h = (r 2 h , u 2 h ) by With these discrete operators, we propose the following family of semi-discrete staggered schemes which can be applied to the linear wave equation with Coriolis source term: where (γ r ) i , (ν r ) i = κra⋆hi

2
, with h i being a length related to triangle T i , e.g. its circumradius, and ν u = κua⋆ 2 represent the parameters that control the diffusion terms and (∇ CR u h • n ) i is given by The classical values for κ r and κ u are either 0 (no damping) or 1 (damping).The term defined in (10) is the standard diffusive upwinding term used in the collocated Godunov scheme in the velocity equation.Note that this diffusion term only acts on the jumps of the normal velocities, contrarily to a full velocity Laplacian diffusion term.This difference is essential in the preservation of discretely divergence free velocity fields that have continuous normal components across cell interfaces, and are thus unaffected by this diffusion term.The notation ∇ CR is here to recall that the righthand side in (10) is the gradient of the Crouzeix-Raviart non-conforming P 1 function with value [(u j −u i )•n ij ] at the midpoint of the edge A ij .We also note that the Low Froude (LF) staggered scheme corresponds to ν r = γ r = 0, the classical Apparent Topography (AT) scheme corresponds to ν r > 0, γ r = 0 and the Modified Apparent Topography (MAT) to ν r = γ r > 0.
Remark 1.In Scheme (9), it is stressed that the additional terms all vanish when applied to a geostrophic equilibrium.This is obvious for the terms that contain ∇ T h r h + ω a⋆ u ⊥ h because the definition of the GE is exactly that this term vanishes.Moreover, the discrete GE implies that u h is the curl of a P 1 continuous Lagrange function on the triangles; these curls have continuous normal components on the edges of the mesh, which implies that the jump term u h • n also vanishes for such a discrete GE.More details are provided in Section 2.4.

Properties of the discrete operators
Proposition 1.With the discrete divergence, curl, gradient and scalar product defined respectively by (6), (7), (5) and (8), we have the following properties for the semi-discrete staggered scheme (9) (These properties are not new findings of this work; they are repeated here without proof for the sake of completeness): i. Discrete integration by parts (energy conservation for the pressure gradient force) ii.Energy conservation for the Coriolis force iii.No vorticity production for the pressure gradient force 2.3.Evolution of the discrete energy Proposition 2. With ν r = γ r and the discrete energy defined with the following expression which means that the Low Froude and Modified Apparent Topography semi-discrete staggered scheme dissipate the discrete energy.
Proof.We take the scalar product of the staggered scheme (9) with q h = (r h , u h ) in order to obtain Classically (see [15]), the term • n ij | 2 ; then, using (11) and (12), and denoting by || • || P the norm associated to the scalar product on the primal mesh, we get which leads to the conclusion.
Remark 2. The fact that γ r = ν r is essential in the calculation above.For this reason, the usual Apparent Topography scheme (ν r > 0, γ r = 0) cannot be proved to dissipate energy with this kind of energy estimation (which does not mean that it does not).

Analysis of the discretized steady-states and their orthogonal subspace
We now define a set of discretized steady-states with staggered variables on triangular meshes by the following expression which is a consistent discretization of the GE (3).Defining further a set A h by we have the following discrete orthogonal decomposition: Proof.For each qh = (r h , ûh ) ∈ E △ ω̸ =0 and arbitrary qh = (r h , ũh ) ∈ R Nr × R 2N , we use the discrete integration by parts formula (11) to obtain Hence, if qh ∈ A h , we obviously have ⟨q h , qh ⟩ = 0 which leads to A h ⊂ E △,⊥ ω̸ =0 .On the other hand, since rh can be arbitrary in R Nr when qh ∈ E △ ω̸ =0 , then the equality rh , rh Remark 3. The discrete orthogonal decomposition allows us to define the following discrete orthogonal projection and we can construct qh by what follows: Let q h = (r h , u h ) be given in R Nr+2N .For all (p h , vh ) ∈ E △ ω̸ =0 , using orthogonality, we have We then use the definition of the discrete steady-states and the discrete integration by parts formula to get As a result, since ph can be arbitrary in R Nr , it is possible to find rh by solving the following linear system Then, by definition of the discrete steady-states, the part of the velocity field in Finally, the orthogonal component is simply given by qh = q h − qh .Moreover, the linear system (18) defines a unique solution since −∇ D h • and ∇ T h are adjoint operators as shown by the discrete integration by parts formula (11).

Well-balanced and orthogonality preserving properties
Definition 1.A semi-discrete scheme is said to be well-balanced with respect to system (2) if Definition 2. A semi-discrete scheme is said to be orthogonality preserving with respect to system (2) if Proposition 4. We have: (i) The semi-discrete staggered type scheme (9) is a well-balanced scheme in the sense that it can capture the discrete steady states (15).
Proof.With the discrete steady state (15), the velocity field can be written as Since we have no vorticity production of the gradient term (see ( 13)), we get On the other hand, since ûh is the curl of a discrete piecewise P 1 , globally continuous function, then its normal jumps through edges vanish: Therefore, the definition of the discrete kernel (15), the divergence free property (20) and the vanishing of the jumps of the velocity field through the edges A ij (21) imply the well-balanced property of the semi-discrete staggered scheme (9).This proves Point (i).
In consideration of the orthogonality preserving property, by taking the discrete scalar product of the semi-discrete staggered scheme with any stationary state qh ∈ E △ ω̸ =0 , we obtain when By using the discrete integration by parts formula and ( 20) , we have Moreover, by simple calculations, we get, thanks to (21): Using a final discrete integration by parts formula and the fact that qh ∈ E △ ω̸ =0 , we get On the other hand, using integration by parts, ∀q h ∈ E △ ω̸ =0 we get Therefore, the semi-discrete staggered scheme is orthogonality preserving.
Remark 4. Proposition 4 tells us that all considered schemes preserve GE because dissipation terms vanish for such steady states.However, the Low Froude and Modified Apparent Topography schemes have an additional advantage over the usual Apparent Topography scheme in that they operate a strict splitting between the GE and their orthogonal subspace.More precisely, by virtue of the orthogonal decomposition of vector fields, it is always possible to decompose the numerical solution into two parts: with qh (t) in the GE space E △ ω̸ =0 and qh in the orthogonal space.Since the Low Froude and Modified Apparent Topography schemes are orthogonality preserving, their orthogonal parts qh (t) are equal to the time evolution of qh (0) for all time t.As a result, we can ensure that qh (t) = P h q 0 h , ∀t ≥ 0 for both schemes.This means that the kernel part of their numerical solutions is always equal to the initial projection P h q 0 h , which is also the final state of the numerical solution when the orthogonal part is damped to zero by the numerical diffusion of the schemes.On the other hand, the classical Apparent Topography scheme is not orthogonality preserving because of its numerical diffusion on the pressure equation.As a consequence, the orthogonal part not only damps out, but also partly moves into the kernel.Therefore, the kernel part of this scheme may be changed at each time step until there is no energy left in the space orthogonal to the kernel.

Low Froude number accuracy
When both a ⋆ and ω scale like 1 M , with M ≪ 1, then the solution of ( 2) is typically a GE (an element in the kernel) plus an orthogonal perturbation of order M .We shall verify under which (sufficient) condition a scheme like (9) mimics this.More precisely, we set Definition 3. Let q ν,h (t) be the solution of the semi-discrete scheme (9) with initial condition q 0 h .This scheme is said to be accurate at low Froude number if Proposition 5.With ν r = γ r , the semi-discrete scheme ( 9) is accurate at low Froude number; this is in particular the case for the LF and the MAT schemes.
Proof.By linearity, the solution of semi-discrete staggered scheme q ν,h (t) can be written as where q a ν,h (t) and q b ν,h (t) are solutions of ( 9) with initial conditions respectively given by q a ν,h (0) = P h (q 0 h ) and q b ν,h (0) = q 0 h − P h (q 0 h ).
Then, we have Moreover, when ν r = γ r , the dissipation of the semi-discrete staggered scheme proved in Lemma 2 leads to the conclusion that ∥q b ν,h (t)∥ ≤ ∥q b ν,h (0)∥.For this reason, the accuracy of the scheme is linked to the behavior of q a ν,h (t).Since the semi-discrete scheme ( 9) is a well-balanced scheme, we obviously have q a ν,h (t) = P h (q 0 h ).Therefore, we obtain ∀t ≥ 0, ∥q ν,h (t) Remark 5. Since it is difficult to prove dissipation of energy for the classical Apparent Topography scheme (see Remark 2), we do not have enough evidence to conclude that this well-balanced scheme is accurate at low Froude number.

Analysis and damping of spurious modes
One of the drawbacks of the triangular B-grid scheme is that the standard scheme without damping (i.e.ν r = γ r = ν u = 0 in (9)) supports spurious numerical modes: the analysis of the numerical dispersion relation on uniform equilateral meshes performed in [13] shows that (non trivially constant) divergence free velocity fields (with vanishing pressure) can oscillate with pulsation ω.This means that there are non trivial velocity fields verifying We analyse in what follows these fields on general triangular grids: taking the ∇ D h • operator on both sides of (23) and using (22), we have that With relationship (7), we conclude that for ω ̸ = 0, it also holds that If u h were a continuous velocity field in a periodic rectangular domain with no holes as is considered here, a vanishing curl and a vanishing divergence would be sufficient to infer that u h is a constant vector field.However, at the discrete level, ( 22) and ( 24) do not imply that u h is a constant vector field: indeed, with N the number of triangles and N r the number of vertices, a cell-centered velocity field has 2N degrees of freedom, while ( 22) and ( 24) sum up to only 2N r = N constraints (this equality holds on a periodic triangular mesh of a periodic rectangular domain with no holes).Therefore there are many such spurious modes and a general construction is given in the associated numerical test in Section 4.3.We shall build on ( 14) to show that these modes are damped to zero by the family of schemes we propose as soon as ν r = γ r > 0.More precisely, with an initial condition r h (0) = 0 and u h (0) such that ( 24) is verified, then by definition (16) and property (17), it holds that (r h (0), u h (0)) ∈ E △,⊥ ω̸ =0 .By virtue of Lemma 4 (ii), we infer that (r h (t), u h (t)) ∈ E △,⊥ ω̸ =0 for all t ≥ 0. This implies by ( 16)-( 17) that Now, we use a general decomposition of discrete velocity fields on periodic triangular meshes with no holes (see [16,Lemma 3.5] with α ≡ 1): there exists two time-dependent discrete fields ψ h (t) and ϕ h (t), with ψ h defined by its values at the vertices of the mesh and ϕ h defined by its values at the edges of the mesh, and a uniform vector field Moreover, in this decomposition, each of the three components is orthogonal to the other two with respect to the scalar product ⟨ , ⟩ P .Now, computing the second term in the right-hand side of ( 14), we obtain, thanks to (26) and the orthogonality of ∇ T h with ∇ CR × and uniform fields where ν min is the minimum value of ν r over the mesh.Developing the first term in the right-hand side of (27) we obtain In order to evaluate the second term in the right-hand side of (28), we first use (25), then a discrete integration by part for the curl operators (that can be proved just like (11) was proved), then (26) and the fact that this decomposition is orthogonal, to obtain Plugging ( 29) into ( 28) and then into (27), and using ( 14) we obtain the following inequalities which implies and thus the exponential damping of the spurious modes with a damping rate (at least) proportional to ν min and thus to h min , where h min is the minimum of the h i over the mesh.

Analysis of fully discrete staggered schemes
We present in this section a time discretization for scheme (9).Of course, one may first think about the splitting scheme where we deal with the homogeneous system in the first step and take the source term into account in the second step.However, this naive splitting scheme is unable to capture the steady states.Therefore, in this section, we focus on the analysis of the so-called one step scheme.

The fully discrete one step scheme
Since a totally explicit discretization of the Coriolis term (u ) would be unstable [1], we need to discretize it implicitly enough.In order to achieve this, let us introduce two new parameters θ 1 and θ 2 corresponding to the time discretization of the Coriolis source term.Let us also introduce τ 1 and τ 2 corresponding to the time discretization of the divergence term in the pressure equation.For convenience, let us denote Then, we introduce the θ-τ scheme which is given by Remark 6.The θ-τ scheme (30) is not totally explicit but we can still calculate the velocity and pressure updates without having to solve any linear system, although the velocity field u n+1 appears in the pressure equation: we first compute the updated velocity field from the momentum equation, and then use it to compute u τ in the pressure equation.
3.2.Well-balanced and orthogonality preserving scheme Proposition 6.It holds that (i) The fully discrete one step scheme (30) is a well-balanced scheme.
(ii) When ν r = γ r , the fully discrete one step scheme (30) is an orthogonality preserving scheme if which means that the velocity field used in the Coriolis source term in the velocity equation and that used to compute the divergence term in the pressure equation should be the same.
Proof.Let us consider at time t n a numerical solution qn h = (r n h , ûn h ) which is in the discrete kernel, which reads (19).We also recall that this implies that ûn h also verifies (20) and (21).We shall show that the numerical solution at time t n+1 is still equal to qn h .First, developing u θ and taking (19) into account to express the pressure gradient, as well as the fact that the jump term vanishes due to (21), the velocity equation of (30) leads to Next, this implies that u τ = ûn and ∇ D h • u τ = 0 by (20).Therefore, the pressure equation of the fully discrete one step scheme (30) reduces to r n+1 k = rn k .This proves (i).Next, let us consider a numerical solution q n h belonging to E △,⊥ ω̸ =0 , so that ⟨q n h , qh ⟩ = 0 for any steady state qh ∈ E △ ω̸ =0 .Taking the scalar product of ( 30) with any qh ∈ E △ ω̸ =0 , and following the proof of Proposition 4, the only terms that remain are We realize that when τ 1 = θ 1 and τ 2 = θ 2 , we get u τ h = u θ h , from which it follows that Therefore, we conclude that q n+1 h ∈ E △,⊥ ω̸ =0 and so the scheme is orthogonality preserving.

Numerical test cases
In this section, we perform several test cases for all proposed schemes: LF, AT and MAT.Moreover, for the sake of comparison, we also run the scheme with the Laplacian diffusion in the pressure equation and with the normal velocity jump diffusion for the velocity equations (PL-VJ scheme): in ( 9), the diffusion in the pressure equation is replaced by −∇ D h • ν r ∇ T h r h and γ r is set to 0. The classical scheme with the full Laplacian diffusion for the velocity equations and no diffusion in the pressure equation (VL scheme) is also tested in this section: in (9), we set ν r = γ r = 0 and then replace the diffusion term (10) in the velocity equation by The crucial difference with (10) is that (10) only involves the normal jump of the velocity field, and this normal jump vanishes on discretely divergence free velocities, as recalled in Remark 1 and Proposition 4.

Well-balanced test case
In this test case, we investigate the behavior of the staggered schemes with a GE as an initial condition.Particularly, we consider the stationary vortex in the square periodic domain T 2 = [−0.5,0.5] × [−0.5, 0.5] with initial pressure r 0 given by r(x, y, t and we construct the discrete initial pressure by interpolating this pressure field at the mesh vertices.Then, we construct the initial velocity field u 0 by using the definition of the discrete kernel (15) so that we can obtain a discrete stationary state.Figure 2 indicates that the classical VL staggered scheme is not well-balanced since it quickly produces some spurious waves in the orthogonal subspace (Fig. 2b); it also strongly damps the kernel part (Fig. 2a).Similarly, The PL-VJ scheme also destroys the geostrophic equilibrium.However, compared to the VL scheme, it produces waves with a smaller amplitude in the orthogonal part and the damping rate of the kernel part is much slower.On the contrary, the LF, AT and MAT staggered schemes are well-balanced as expected.This can be checked because the orthogonal parts of their solutions remain equal to zero (Figure 2b) during the computation and their kernel parts remain constant (Figure 2a).Additionally, Figure 3 shows that the pressure contours of the VL and, to a much smaller extent, PL-VJ schemes are different from those of the other schemes.This is another evidence that they are unable to capture steady-states, contrarily to the LF, AT and MAT schemes.

Orthogonality preserving test case
In this test case, we consider periodic boundary conditions and an initial vector field given by u(x, y, t = 0) = 1 2 exp − 4x 0.4 .5]and we construct the discrete initial velocity by interpolating this velocity field at the cell centers.Then the initial pressure r(x, y, t = 0) is constructed at the mesh vertices by using the definition of the discrete orthogonal subspace ( 16)- (17).In all cases, we choose θ 1 = θ 2 = 1 2 for the time discretization of the Coriolis force in (30) and the value of (τ 1 , τ 2 ) may vary from test to test. Figure 4a shows that the classical VL, PL-VJ and AT schemes are not orthogonality preserving since the kernel components of these schemes are not equal to zero.Moreover, since the kernel part is updated at each time step, it is not a constant in time.Figure 4b shows that the damping rates of the orthogonal parts of the VL, PL-VJ, AT and MAT schemes are larger than that of the Low Froude scheme.Although the PL-VJ and the AT schemes are not orthogonality preserving, these strategies on triangular grid create waves in the kernel that have much smaller amplitudes than those on Cartesian grids (see [6]).Next, Figure 4a was obtained with τ 1 = τ 2 = 1 for the fully discrete LF and MAT schemes, and this generates waves with very small amplitudes in the kernel part; this underlines that in order to obtain schemes that are really orthogonality preserving, we need to use a velocity field in the Coriolis source term that is equal to that used for the computation of the divergence term in the pressure equation.Indeed, Figures 5a and 5b clearly show that the LF-θ-τ and the MAT-θ-τ schemes with 2 are perfectly orthogonality preserving since these strategies do not create waves in the kernel subspace during the computational process, as proved by Lemma 6.

Spurious mode test case
In this test case, we consider an initial condition given by a spurious mode of the triangular B-grid when no damping is used (i.e.(9) with ν r = γ r = ν u = 0), and we shall observe that this mode is now correctly damped to 0. As described in Section 2.7, this initial condition satisfies In order to construct such a spurious mode, we recall the discrete Hodge decomposition on periodic triangular grids that is given by ( 26) (taken at time t = 0): (a) ∥q∥ of LF-τ scheme (b) ∥q∥ of MAT-τ scheme Figure 5: Orthogonality preserving test case: evolution of the kernel part of the LF-θ-τ and MAT-θ-τ scheme with θ 1 = θ 2 = 1 2 and various values of (τ 1 , τ 2 ).Orthogonality preservation is only achieved when τ = θ.
But the role of the curl and gradient operators can be exchanged so that there also exist (µ h , λ h ) defined respectively by their values at the vertices and at the edges of the mesh such that In order to have a velocity field that is at the same time a non-conforming gradient and a non-conforming curl, a possibility is that one picks any random (periodic) α h defined at the edges of the mesh and first set v h = ∇ CR α h .Then, one computes the second Hodge decomposition of v h : Then, if we set u h = v h − ∇ T h ζ h , we obtain that u h is both a non-conforming gradient and a non-conforming curl: where ζh is defined at the edge midpoints as the half-sum of the values of ζ h at the vertices of the edges.Note that, in order to decompose v h , we have to solve a P 1 Lagrange finite element solution of the Laplace equation: Having constructed a spurious mode, we use it as an initial value for the various schemes developed in this article, as well as the "Inertial Oscillation" (IO scheme) which is simply the θ-scheme (with θ = 1  2 ) applied to (23), since this corresponds to the standard B-grid scheme applied to the spurious mode (r = 0 and ∇ D h • u h = 0) when no damping is used.We recall that these spurious modes belong to the orthogonal space ( 16)-( 17) since both their pressure components and the curls of their velocity components vanish.So their initial projections in the kernel space (15) also vanish.According to Figure 6a, we can see that the three non-orthogonality preserving schemes VL, PL-VJ and AT immediately create some waves in the kernel part.Moreover, concerning the long time behavior, the kernel part of the VL scheme is damped by the Laplacian diffusion while it seems to be a constant for the PL-VJ and AT schemes.On the other hand, the LF and MAT schemes do not create any wave in the kernel.Figure 6b shows that the orthogonal part is damped by all the numerical schemes with numerical diffusion.As a final remark, we can stress that Figure 6c illustrates the fact that damping spurious modes supported by the standard B-grid scheme (the IO scheme) can be achieved by schemes that preserve GEs and their orthogonal space (the LF and MAT schemes), which is a major improvement over the usual damping solution (the VL scheme).

Accuracy at low Froude number test case
We now consider an initial condition close to the discrete kernel, up to a perturbation of size M .This initial condition is simply given by , where q0 h stands for the kernel part given in Section 4.1 and q0 h is the orthogonal part considered in Section 4.2.In that case, the exact solution remains close to q0 h , in the sense that their difference remains of size M for all time (see Lemma 5).So we expect efficient schemes to follow this behaviour.However, Figures 7a and 7b indicate that the classical VL and PL-VJ schemes are not accurate at low Froude number because the norm of ∥q(t) − Pq 0 ∥ has a size which is much larger than M and independent of the value of M .These figures also illustrate that the solutions of the VL and PL-VJ schemes quickly move far from the kernel without regard to the size of parameter M .By contrast, the proposed LF, AT and MAT schemes are accurate at low Froude number because the norm of the total deviation remains of order O(M ).
with periodic boundary condition, a ⋆ = ω = 1 and the domain is the square 5].In this test case, we aim to investigate the long time behavior of the numerical solution when the initial condition is far from the GE. Figure 8a presents the projection into the kernel of the initial condition.This is the target solution for the simulation in long time since the numerical diffusion will damp the orthogonal part.Figure 8 clearly shows that in long time, the solution of the VL and PL-VJ schemes tend to a trivial steady state which consists of only a constant pressure field.On the contrary, the LF, AT and MAT schemes tend to the GE, since the final state is similar to projection Pq 0 of the initial condition shown in Figure 8a.
On the other hand, Figure 9 shows that the kernel components of the classical VL or PL-VJ schemes are damped during the simulation while it is a constant for the LF and MAT schemes, and nearly a constant for the classical AT scheme.Moreover, the total deviation (distance between the numerical solution and the initial projection) of VL or PL-VJ scheme decreases in early times, but eventually increases in long time, while it tends to zero with the other schemes.This is another evidence to conclude that the VL and PL-VJ schemes will tend to a wrong component in the kernel, while the other schemes have solutions that, as expected, tend to the correct GE.The fact that the classical AT scheme does not capture a projection in the kernel which is constant in time is due to the fact that it is not orthogonality preserving, and so there is an exchange of energy between the orthogonal space and the kernel.

Anticyclonic eddy propagation test case for Rossby waves
As a first extension of what was presented in the simplified case of the constant coefficient linearized system (2), we consider in this test case the Rossby waves generated by the variation of the Coriolis force with latitude.In particular, we focus on the β-plane approximation for the Coriolis source term given by ω = ω 0 + βy.We aim to investigate how the above schemes influence the numerical propagation of the Rossby waves.
This test case is motivated by the second test in [17].Here, all parameters are scaled according to our setting.Particularly, the computational domain is a rectangular [−1.0, 1.0] × [−0.6, 0.6] with periodic boundary conditions.The initial pressure is a Gaussian distribution at the center of the domain and the initial velocity field is chosen so as to satisfy the geostrophic balance a ⋆ ∇r = −(ω 0 + βy)u ⊥ .Therefore, the initial condition is given by We choose A = 0.6, B = 0.13, a ⋆ = 1 as well as the Coriolis parameters ω 0 = 4, β = 4/3.According to Figure 10, the anticyclonic vortex evolve purely westward.Qualitatively, the results with the MAT scheme presented here are in good accordance with those published in [17], and much less diffused than those obtained with the VL and PL-VJ schemes.

Extension to non-linear equations, non-flat topography and hybrid velocity dissipation
We propose the following extension of the presented schemes for the inclusion of the non-linear and topography terms; moreover, we discuss how to recover full Laplacian velocity damping, with a possible tuning of the respective weights of damping the normal and tangential components of the velocity jumps at the cell interfaces.We first define a further dissipation operator acting on the tangential components of the velocity jumps: where n ⊥ is a unit vector orthogonal to n.Then, using q = hu as main unknown, the scheme reads: Once (h n , q n ) is known, then compute in which, for each cell T i , the interpolated value (h c ) i is equal to one third of the sum of the values of h at the vertices of T i .There are several options for the definition of the ∇ T h • operator which we do not detail here.In practice, we chose to define this divergence on the primal cells as the mean value of the three divergences ∇ D h • that can be computed at the vertices of each cell.Note that now, the diffusion terms vanish on fields that verify a discrete GE taking into account the topography and the jump terms still vanish for discretely divergence free velocities when the additional tuning parameter ν ⊥ vanishes.On the contrary, we recover the full Laplacian velocity damping when ν ⊥ = ν u .Indeed in that case Scheme (32) is the fully discrete version of a semi-discrete scheme (continuous in time, discrete in space) devised so that, apart from the advection terms, it dissipates energy.Indeed, when we take the scalar product on primal cells of the first equation of (32) with q hc and then sum with the scalar product of the second equation with g(h + b), we obtain (omitting the time derivatives and advection terms) The sum of the first two terms vanishes like in the linear case.The third and fourth terms are always negative.In addition, when ν r = γ r and thanks to a discrete integration by part, the last two terms sum up to which is negative as long as all values of (h c ) remain positive.This energy dissipation ensures in practice the stability of the scheme.As a test case, we consider the stationary vortex of the nonlinear rotating shallow water equation studied in [18,19] on the domain [−0.5, 0.5] × [−0.5, 0.5] with periodic boundary condition.The initial velocity field is given by where r stands for the distance to the center of the domain and χ is the characteristic function of a domain.The corresponding initial water depth is the solution of the following ODE Let us note that when the coefficient ε is small, the nonlinear term in the definition of the water depth (34) is much smaller than the linear term and the stationary solution is very close to a geostrophic equilibrium.Note that we do not claim that the scheme presented here preserves this initial condition; however, we show that for small values of ε, it behaves much better than other schemes which do not preserve geostrophic equilibriums.In this test, in addition to the non-linear terms, we also consider the topography b(x) given by b(x, y) = 1 2 e −40(x 2 +y 2 ) (35) in order to show that the presented schemes also behaves correctly in presence of a non flat topography.Figure 11 (for ε = 0.01) shows the contours of (h + b) after some time of simulation and it is shown that the solution of the MAT scheme is closer to the initial condition than the PL-VJ scheme, which produces more diffused solutions.We used ν ⊥ = 0 in this test.

Conclusion
In this work, we propose new staggered type schemes on triangular cells that are accurate for the simulation of flows near the GE.The construction of these schemes is based on the adaptation of the Low Froude and Apparent Topography strategies to triangular grids.Theoretical analysis and numerical results in the linearized setting with 20 constant Coriolis parameter show that the new schemes are well-balanced.Additionally, unlike the classical VL, PL-VJ and Apparent Topography method, the Low Froude and the Modified Apparent Topography schemes are orthogonality preserving under an appropriate discretization in time.This seems particularly important to really separate the GE from the orthogonal perturbations around it.An important feature of the proposed schemes is that they are able to damp the spurious modes associated to the staggered triangular unknown arrangement without destroying the GE.Encouraging preliminary results show that the schemes can be used in and extended to more realistic and complex situations, like non-constant Coriolis parameters and non-linear dynamics; this needs to be further explored in future works.Coupling with more dissipative schemes was also shown possible via a tuning parameter that can be adjusted according to the expected physical regime.Additional effort could also be dedicated to the optimal time step choice for these schemes and to the incorporation of ideas drawn from [20,21] for the energy control in fully discrete extensions to the nonlinear case.Another possible extension would be to consider higher-order schemes like those in the P disc k -P k+1 family discussed in the introduction.

Appendix A. Analysis of the behavior of collocated schemes
Since collocated schemes have a long and successful history for the numerical approximation of systems of conservation laws such as (1) with ω = 0, one could have the natural idea to try to extend them to the treatment of the Coriolis source term (ω ̸ = 0).The aim of this Appendix is to show the difficulties one encounters in trying to devise collocated schemes that would verify a discrete version of the GEs on triangular meshes.
First, as far as Cartesian meshes are concerned, it is well-known that the classical Godunov scheme applied to (2) will fail to capture a discrete equivalent of the GE (3) [3].Without going into details, the reason is that the numerical diffusion of that scheme, that is necessary for stability reasons, also destroys GEs and, as a consequence, generates spurious modes at a rate which is proportional to h M , with h being the mesh spacing.In the 1D case, the only reason for the accuracy problem is linked to the numerical diffusion in the pressure equation [1].To remedy this issue, the authors of [3] introduce the so-called Apparent Topography method by adapting the hydrostatic reconstruction of [22]; this works well in the 1D case and extensions to the 2D Cartesian case were considered in [23,24,25,26].Other relevant works on this issue can be found in [27,28,29,30,31].
However, in the 2D Cartesian case, the numerical diffusion in the velocity equations also comes into play since it does not vanish on velocity fields verifying a discrete version of (4).This issue was discussed in [6], in which modified Godunov schemes on Cartesian meshes were investigated that correctly capture discrete GEs.Extensions to staggered schemes on Cartesian meshes are investigated in [32].
We now investigate what happens if the use of triangular grids is preferred to Cartesian grids, for example because triangular grids are much more flexible than Cartesian ones for accurate discretizations of complex detailed geometries and for local mesh refinement.An essential characteristic of triangular meshes is that the numerical diffusion of Godunov schemes on the velocity equation is different in nature (as compared to Cartesian grids) and does not affect velocity fields that verify a discrete version of (4).References [15,16] give a detailed explanation of this, which may be summarized by the fact that such discretely divergence free velocity fields can be characterized as the curls of P 1 Lagrange functions defined by their values at the vertices of the triangular cells.This has two major consequences: we will first show that this is hardly compatible with a discrete version of (3) when the pressure gradient is defined like in the Godunov scheme with cell-centered pressures; on the other hand, we shall secondly show that this is obviously compatible with a discrete version of (3) when the pressure gradient is defined with pressures located at the vertices.
To begin with, let us note that the semi-discrete collocated scheme on triangular mesh can be written as where In this scheme, κ r and κ u represent parameters that can be used to tune the diffusion terms.The classical Godunov scheme corresponds to the case κ r = κ u = 1.We refer to, e.g., [15] for the construction of the scheme applied to the homogeneous wave equation.Here, we take into account the effect of the Coriolis force in the right-hand side of the semi-discrete scheme and we analyze whether steady-states of (A.1)-(A.2),i.e., fields in the kernel of the operator L κ,h , are able to correctly represent discrete versions of the GE.
Proposition 7. On a triangular mesh, with the collocated Godunov scheme, we have: KerL κr̸ =0,h = q h := r h u h ∈ R 3N such that ∃a ∈ R : r i = a and u = 0 .(A.3) Moreover, we also have with the following definitions for C 1 h and C 2 h : where C 0 ♯ (Ω) denotes the space of periodic continuous functions over Ω and Proof.Denoting by ⟨• , •⟩ the scalar product weighted by the cell areas |T i | and following the properties of the collocated scheme [15] and using the energy conservation for the Coriolis force ⟨u ⊥ h , u h ⟩ = 0, we obtain We now suppose that κ u ̸ = 0 and we will consider the influence of κ r on the structure of the kernel KerL κ,h .From (A.5), a necessary condition for a velocity field to be in the kernel is This obviously implies that (using, for the second equality below, the fact that Aij ⊂∂Ti When κ r ̸ = 0, equation (A.5) also implies that ∀(i, j), r i = r j and it follows that r must be a constant: ∃a ∈ R such that r i = a for all i ∈ [1, N ].Then, the velocity equation of L κ,h q h = 0 leads to ωu ⊥ i = 0 ∀i.Therefore, we get (A.3).
In the other case , when κ r = 0, the second equality in (A.7) implies that the pressure equality in L κ,h q h = 0 is verified.Moreover, we recall from [15, Lemma 5.1] that (A.6) is equivalent to the fact that there exists a Lagrange piecewise P 1 conforming function, denoted by ϕ L h , defined by its values at the vertices of the mesh and periodic over the computational domain, and constants (b, c) ∈ R 2 such that on each triangle T i there holds Moreover, from the velocity equation of L κ,h q h = 0 and (A.6), we get that which means that q h ∈ C 1 h .In addition, by multiplying (A.9) by |T i | and summing over all i ∈ [1, N ], we obtain by periodicity that i=1,N |T i |u i = 0, which implies that (b, c) T in (A.8) vanishes.Therefore, (A.8) with b = c = 0 implies that q h ∈ C 2 h .In conclusion, we obtain and the converse inclusion is also true.
Remark 7. The kernel given by (A.3) clearly shows that the standard Godunov scheme (κ r ̸ = 0) is not able to preserve general geostrophic equilibriums.On the other hand, the situation is less clear when the diffusion in the pressure equation vanishes.Indeed, let us consider (A.4).An important property is that for all i ∈ [1, N ] provides 2N constraints while r has N degrees of freedom but ϕ L h only has a number of degrees of freedom equal to the number of vertices of the mesh.In a periodic triangular mesh with no holes, the number of vertices is half the number of cells, so that (A.10) is actually probably too constraining (2N constraints for only 3N   2   degrees of freedom) to admit non trivial solutions.
With the previous remark, it becomes obvious that one possibility is to define the pressure field at the vertices and replace the gradient of rCR h by the gradient of the conforming Lagrange P 1 function defined by these values at the vertices.Of course then, we have to change the pressure evolution equation accordingly.This staggered scheme, and several variations of it, is what we explore in the body of this article.

Figure 2 :Figure 3 :
Figure 2: Vortex test case: evolution of the kernel and orthogonal parts.

Figure 6 :
Figure 6: Energy evolution of the kernel part (a), of the orthogonal part (b) and of the full solution (c).

Figure 7 :
Figure 7: Evolution of the total deviation from the initial condition, when the initial condition is close to the discrete kernel.

Figure 8 :Figure 9 :
Figure 8: Pressure r(x, y, t) at time t = 150 obtained from staggered type schemes and the initial projection (top left).

Figure 10 :
Figure 10: The contours of h(x, y, t) at time t = 20 obtained from staggered type schemes for the case β = 4/3.

Figure 11 :
Figure 11: The contours of h(x, y, t) + b(x, y) at time t = 20 obtained from staggered type schemes for the case ε = 0.01.
r j 2 n ij is exactly the gradient on the cell T i of the non-conforming P 1 finite element function rCR h which is defined by the values ri+rj 2 at the midpoints of the mesh edges A ij .So C 1 h tells us that u h is, up to a multiplicative constant, the curl of the non-conforming P 1 function rCR h , while C 2 h tells us that u h is the curl of a conforming P 1 function ϕ L