Numerical approximation of the shallow water equations with Coriolis source term

. We investigate in this work a class of numerical schemes dedicated to the non-linear Shallow Water equations with topography and Coriolis force. The proposed algorithms rely on Finite Volume approximations formulated on collocated and staggered meshes, involving appropriate diﬀusion terms in the numerical ﬂuxes, expressed as discrete versions of the linear geostrophic balance. It follows that, contrary to standard Finite-Volume approaches, the linear versions of the proposed schemes provide a relevant approximation of the geostrophic equilibrium. We also show that the resulting methods ensure semi-discrete energy estimates. Numerical experiments exhibit the eﬃciency of the approach in the presence of Coriolis force close to the geostrophic balance, especially at low Froude number regimes.


Introduction
The question of the accuracy of numerical schemes around stationary solutions or/and in asymptotic regimes has been a subject of great interest over the last two decades, see the seminal works [4,10,12] in late nineties and the books [5,9] ten years later.In the context of geophysical flows and more particularly for finite-volume methods applied to shallow water equations, a lot of works have been devoted to the accuracy around the lake-at-rest equilibrium and more recently extended to nonzero velocity one dimensional stationary states.The accuracy of numerical schemes around geostrophic equilibria which are of fundamental importance for large scales atmospheric or oceanographic flows was less investigated, see [16] for a general introduction to geophysical rotating fluid dynamics.To our knowledge, the first work in this field is due to Bouchut, Le Sommer and Zeitlin [6] but was fully accurate only for one-dimensional flows, as exhibited in [2] where accurate and stable Godunov type schemes were designed for linear two-dimensional shallow water equations.Recently two independent works [13,15] proposed an IMEX type scheme for fully nonlinear equations.
Here we aim at designing explicit schemes that are proved to be accurate and stable in the nonlinear framework.Our work is mostly based on the ideas developed in [2] for the linear rotating case and in [7,8] for the nonlinear non-rotating case.We investigate schemes applied to collocated and staggered schemes.
Let Ω be an open bounded domain of R2 and let T > 0. The nonlinear Shallow Water equations with Coriolis force formulated on Ω × (0, T ) read: where h is the water height and u = (u x , u y ) the horizontal velocity 1 .The Coriolis force is accounted for in the momentum equations through the angular speed ω.Following [7,14], the pressure forces appear under a non conservative form through the scalar potential φ = gh, where g is the standard gravity constant.For the sake of simplicity, a flat topography is considered in the present work, but the proposed approaches naturally extended to varying bottoms. 2  It is well-known that the total energy of the system decomposes as E = E + K where stand respectively for potential and kinetic energies.We recall that the energy E plays the role of a mathematical entropy associated to the hyperbolic system (1) and regular solutions satisfy the following conservation law: When developing numerical methods, main objectives are accuracy and stability.To get stability, a crucial objective is to build numerical approximations satisfying a discrete counterpart of (2) that ensures that the discrete energy is nonincreasing.To achieve this, a general strategy is to consider a sufficient amount of numerical diffusion in the scheme.But in some physical contexts such as low Froude number regimes or near specific stationary states, these diffusive terms may considerably degrade the accuracy of the approximations and specific schemes are needed.Here we are interested in flows around the geostrophic balance To address such an issue, based on the linear case [2], and on the non-linear case without Coriolis force [7,8], we propose a numerical approach involving discrete versions of these equilibria in the numerical fluxes.As a preliminary step, the strategy can be understood at the continuous level by investigating how the model (1) behaves with respect to some generic perturbations (q, π): Hence q and π must be respectively seen as (small) perturbations with respect to the flow rate and to the hydrostatic pressure.
The solutions to the modified equations (4) satisfy the following energy balance: which motivates a choice for q and π involving resp.the quantities ∇φ + ω u ⊥ and div u.Let us remark that these quantities govern the geostrophic equilibrium (3) associated to System (1) linearized around the steady state ( h, u) = (h 0 , 0) for a constant h 0 : From a numerical point of view, diffusion terms are thus expected to have regularizing effects in the sense that they allow to recover a discrete counterpart of (5).Moreover, such terms are intended to vanish close to the geostrophic equilibrium, which must improve the quality of the approximations in this regime.Hence, the present approach can be seen as a non-linear extension of the linear case [2].
The paper is split into two parts: in the first section, we present the two geometric frameworks (collocated vs. staggered), we build some discrete operators with mimetic properties and we propose two well-balanced numerical schemes whose semi-discrete (in time) versions are proved to be stable.The second section is dedicated to the assessment of those schemes by means of the stationary vortex test case.
1. Design of the numerical schemes

. Mesh and notations
To discretize System (1), we consider a tessellation T of the computational domain Ω ⊂ R 2 made of nonoverlapping rectangular cells of sizes ∆x, ∆y.The set of all edges of the mesh is denoted by E and the set of vertices by V.
• A generic cell of T is denoted by K, its area by m K and its boundary by ∂K.
• The length of any edge e is denoted by m e .
• Given a cell K and an edge e ∈ ∂K, K e is the neighbouring cell to e (other than K) and n e,K the outward normal pointing to K e .• Given a vertex x, a dual cell is built by joining the centers of the neighbouring cells.Its area is denoted by m x .Notations are pictured on Figs. 1.In particular, the discrete Green formula yields In what follows, for any piecewise constant function ϕ we will denote and introduce the vector quantity In the same way, for a piecewise constant vector function Φ, we will note Finally, the quantity will refer to the the positive and negative parts of any scalar function ϕ.Before describing the schemes, we also need to define the following discrete divergence and gradient operators: x if (Φ K ) K∈T is rather considered, then Φ e is replaced by Φ e in (7) and the corresponding operator is denoted by div c K Φ; We can also define a divergence operator at the vertices when (Φ K ) K∈T is considered by applying Formula (7) to the dual cell (see Fig. 1(b)): • In the same way, considering (Φ e ) e∈E and (Ψ K ) K∈T , the upwind divergence operator applied to the formal product Ψ ⊗ Φ is defined by Likewise, if (Φ K ) K∈T is rather considered, then Φ e is replaced by Φ e in (9) and the corresponding operator is denoted by div c,up K (ΨΦ); • Given a scalar function (φ e ) e∈E , the discrete gradient reads Likewise, if (φ K ) K∈T is rather considered, then φ e is replaced by φ e in (10) and the corresponding operator is denoted by ∇ c K φ.If (φ x ) x∈V is considered, then φ e is replaced 1 2 x∈∂e φ x and the corresponding operator is denoted by ∇ x K φ.In particular, we have the mimetic properties (see [2])

.2. Numerical scheme
In this collocated version, all unknowns are located at the center of the cells except the perturbations (q located at the middle of the interfaces and π at the vertices).
With the discrete operators (7-8-9-10) introduced above, the scheme reads, in the semi-discrete case: where the interface fluxes are defined at the level of the edge e as: with and where γ and ν are positive dimensionless constants, and Λ is a positive characteristic velocity. 3  Remark 1.1.Some remarks must be done on the formulation of the scheme: (1) The origin of the expression (14) for the perturbation q will be explained in the proof below.We first mention that it is consistent with the geostrophic equilibrium (it is aimed to vanish at the equilibrium).Second, we only consider the normal component of u e in order to have the correct number of constraints.Moreover, as the operator δφ e is along the normal vector, q e = 0 would imply (without the normal vector) a too restrictive kernel (characterised by u K,x = 0 or u K,y = 0 depending on the edges).
(2) As m K = ∆x ∆y in this case, q e does not actually depend on K.
(3) Let us now explain why there is a 1  2 coefficient in the correction term of the right hand side of (12).The previous remark combined to the cartesian structure induces that there are only two terms by component in the mean e q e . 4

Semi-discrete energy
We now show this scheme ensures a discrete counterpart of (5) through semi-discrete mechanical energy estimates.To obtain such a result, we first need the two following lemmas, describing the evolution of potential and kinetic energies. 3Typically, we take Λ = max K u 0 K + gh 0 K . 4To extend the scheme to a non-cartesian grid, the correction term e c e,K qe must be consistent with q K with coefficients c e,K such that c e,K m K independent from K so that K m K A K K = 0, and Lemma 1.2 (Semi-discrete potential energy identity).We set E K = 1 2 g (h K ) 2 for K ∈ T. Then: where the interface value for φF is defined as: and the discrete conservative term: Proof.Using the definition of E K and the mass conservation law (1), we get: Then, invoking the decomposition φ K = φ e + 1 2 (φ K − φ Ke ) and the expression of the numerical fluxes (13): Lastly, using the relation hu e = h K u K + 1 2 (h Ke u Ke − h K u K ) and the discrete Green formula (6), basic computations yield: which allows to conclude.
Lemma 1.3 (Semi-discrete kinetic energy identity).We set with anti-symmetric contributions given by: Proof.We write: First, we remark that: Then, after some basic computations, we get the relation The second term of the right hand side being non-positive, we get the announced result inserting (21) and ( 22) into (20).
Gathering relations ( 16) and (18), we obtain the following estimate for the total energy With the choices ( 14) and ( 15), we finally obtain a discrete counterpart of (5) thanks to (11).
Remark 1.4.The contributions A E K (17) and A K K (19) involve anti-symmetric terms and therefore have no impact on the total energy budget: up to boundary terms.In fact, these quantities can be interpreted as a bias to the conservative parts of potential and kinetic energies respectively. Then: Theorem 1.5 (Semi-discrete mechanical energy balance).The scheme (12,13,14,15) satisfies the following semi-discrete inequality:

Steady states
Let us first assess the ability of the scheme to capture the lake at rest, characterized by φ = cte and u = 0.One observes that the stabilization terms q e and π x vanish in such a configuration, leading directly to the exact preservation of the steady state in (12).
Then focus on the geostrophic equilibrium.In this case the nonlinear study is particularly complex and we choose to investigate the linear case.Although it leaves aside nonlinear interactions, this strategy may not seem irrelevant since most of low Froude number regimes encountered in practical contexts are governed by linear dynamics.Linearizing the scheme (12) around ũ = 0 and h = h 0 for some constant h 0 > 0, we get: where )n e,K and πx = −νΛ max{∆x, ∆y}div c x u.
One should observe here that the discrete geostrophic equilibrium, characterized by: =⇒ div c K u K = 0 does not imply q e = 0 in the general case.Conversely, q e = 0 =⇒ π x = 0 according to [2] but does not imply ∇ c K φ + ωu ⊥ K = 0.As exhibited in the numerical validations, although it does not exactly preserve the geostrophic balance, the present scheme behaves much better than standard Godunov-type approaches.However, as detailed in the next section, it is possible to correct such a failure by means of a staggered mesh.

Mesh and notations
On the basis of the primal mesh T, a dual grid T * can be built with respect to the Rannacher-Turek (RT) or Crouzeix-Raviart finite-element spaces, according to the formalism described in [1].To set ideas, we focus here on the RT case considering a primal mesh T made of regular quadrilateral elements like in the previous section, but the following framework also applies to simplicial meshes or mixed triangular/quadrilateral elements.
• A dual cell D ∈ T * is associated to each internal primal edge e as the quadrilateral admitting e as diagonal and the mass centers of the two neighbouring primal cells K and K e as vertices (see Fig. 2). 5 • This dual element D is the union of the two triangles D K , D Ke separated by the edge e, whose areas will be denoted m D K and m D Ke , following the notations of the primal mesh.• Given a dual cell D and an edge e ⊂ ∂K, D f is the neighbouring cell to f (other than D) and n f,D the outward normal pointing to D f .In this context, the water height evolves on the primal mesh T as for the previous scheme while the velocity field is located on the dual elements D ∈ T * .
By analogy with the collocated frame, the upwind discrete divergence is defined on a dual element as: 5 The edge e and the corresponding diamond cell D will be used indifferently.
while the dual gradient is defined through the following formula: Note that this definition ensures the grad/div duality with respect to the L 2 inner product in the present staggered environment (for instance, see [11]):

Numerical scheme
The scheme we consider is the following: where the primal fluxes are defined by interface by: π K = νΛ max{∆x, ∆y}h K div K u.
(28) As in the collocated frame, γ and ν are positive constants and Λ is consistent with a velocity.Considering the momentum equations (25) and the diffusion term (27), it is necessary to define a dual water height h D .Following [1], we compute It is then possible to choose dual fluxes F f (involved in the divergence term of the momentum equations ( 25)) so that a discrete law of mass conservation also holds on diamond cells: which can be rewritten in the more compact form: As detailed in [1], these dual fluxes are expressed as linear combinations of the primal ones.These relations are purely geometrical and formulated locally, depending on the nature of the staggered mesh.Such a relation is mandatory to derive a kinetic energy budget as shown in the next paragraph.
Remark 1.6.Let us notice that in (27), both components are taken into account in the definition of q D unlike (14) where only the normal component was involved.In the present case, q D = 0 implies that the normal component of the velocity field vanishes (since the gradient is normal to the interface) but the value of interest is the reconstructed mean velocity Here u D is an unknown in itself (unlike in the collocated case where u e is reconstructed from the unknowns u K ).

Semi-discrete energy
Indeed, defining the dual kinetic energy by K D = 1 2 h D u D 2 and using the relation we have, with (29) and applying the relation ( 22) on the dual mesh: Then, looking at the potential energy, basic computations give: We are now left with the estimates (30) and (31), which are not expressed on the same mesh.A convenient way to understand the global impact of these quantities on the energy budget is to integrate them over their respective meshes, setting: and focus on the residual contributions at the level of dual elements.We first remark that, thanks to the duality relation (24), the second term of the right hand side of (30) furnish a first stabilizing term: according to the definition of π K (28).Then, considering a generic element D ∈ T * , the non conservative terms in the kinetic energy balance (30) result in the following quantities on the left and right hand sides: On the other hand, the term issuing from (31) needs to be counted twice (contributions of the elements K and K e sharing the edge e), which yields: Using the definition (23), we easily verify that the terms of the left hand side are in exact balance and the remaining term reads: Replacing q D by its expression (27), we have thus established the following result: Theorem 1.7 (Semi-discrete mechanical energy balance).The scheme (25-26-27) satisfies the following semidiscrete inequality: As a remark, defining the quantity K K = 1 m K e⊂∂K m D K D , which can be understood as a kinetic energy associated to the cell K, the previous result can be reformulated under local energy estimates.Note that the present formalism embeds the implicit and explicit approaches recently proposed in [8] without Coriolis force.

Steady states
The linearized version of the scheme (25) is the following: where From this, we immediately see that discrete kernel of the system corresponds to the geostrophic balance expressed on the dual grid T * : ∇ D φ + ω u ⊥ D = 0 , and is consequently exactly preserved by the linearized scheme.We may indeed check that div K (∇ ⊥ D φ) = 0 by construction which implies div K u = 0 when the geostrophic equilibrium holds.

Time discretisation
For the discretisation in time, fluxes are taken explicit.We only remind the reader of the standard discretisation of the right hand side (Coriolis force).The ODE is solved by means of the θ scheme Here we choose θ u = 1 and θ v = 0 so that the system is solved explicitly.The time step is chosen following [2] such that 2. Numerical assessments of the schemes

Presentation of the test cases
To assess the numerical schemes designed in this work, we consider the stationary vortex described in [3].Set in the domain Ω = [−0.5,0.5] × [−0.5, 0.5], the equations are supplemented with the initial velocity field given in polar coordinates The initial water height is then computed so that the steady version of the nonlinear equations is satisfied, i.e.such that As mentioned in [2], this case induces a Froude number of order O(ε).

Numerical results
Given the initial condition prescribed in § 2.2, we compare the robustness of several numerical schemes.The data are a steady solution of the nonlinear equations but their approximation on the grid is not necessarily preserved by any numerical scheme.We remind that the present paper was aimed at designing a scheme that is robust in the low Froude number regime by means of diffusion terms based on linear equilibria.We observe on Figure 3 that the smaller ε, the more accurate the solution obtained with the collocated scheme (12) while the accuracy of the standard Godunov scheme is barely improved.This is confirmed by Fig. 4 where a kind of relative error is plotted: where N is the index of the final time and h = max K h 0 K .The Godunov scheme broadly provides results of which the accuracy is independent from ε while the modifications proposed in this paper leads to a second order scheme until the accuracy is damaged by the truncation error driven by ∆x.
These results tend to show that the target is reached insofar as the resulting scheme is robust in the low Froude number regime where the standard Godunov scheme was inaccurate.We mention that both corrections terms are necessary to obtain stable results which implies γ > 0 and ν > 0.

Conclusion
We managed to construct some corrections to the standard Godunov scheme so that the resulting collocated and staggered schemes are accurate and energy-decreasing around the geostrophic equilibrium in the low Froude number regime.The stationary vortex test case go to show the improvement due to these new terms.
The proofs are given in the semi-discrete case.In particular, the semi-discrete energy inequality requires γ ≥ 0 and ν ≥ 0 so that it decreases.But numerical experiments showed that it is necessary to have them nonzero.Future works will thus deal with the stability analysis of the fully discrete version.Moreover, further numerical investigations are needed to ensure the schemes behave well in a large range of applications.
∩ ∂Ke ∆x ∆y (a) Primal cell K; edge e = ∂K ∩ ∂Ke; outward normal n e,K to e pointing to Ke

Figure 2 .
Figure 2. Geometric settings: for the primal mesh with cells K and K e , refer to Figs. 1 -for the dual mesh with diamond cells D and D f separated by the interface f = ∂D ∩ ∂D f , n f,D is the outward normal of f pointing to D f .

Figure 3 .
Figure 3. Stationary vortex: cross-section of the water height at final time for different values of ε; comparisons between the exact solution, the numerical solution obtained with the classic Godunov scheme, the collocated scheme (12) (A-B-C) and the staggered scheme (25) (D-E-F)

Figure 4 .
Figure 4. Error in L 2 norm between the numerical solution and the exact solution at final time