domain: sde chim sdv shs

text:
© A-M Guihard-Costa

Jean-Philippe Bedell

Eduardo-A Ferreira

D A Silva

M Guihard-Costa

Coordonnées

SOCIAL-ECOLOGICAL FRAMEWORK

The development of an industrial chemical complex (for ammonium sulphate, nitric acid, ammonium nitrate, synthetic resins) has caused extensive pollution of the environment of the Aveiro "ria" (a major lagoon located in the northern Portuguese coast ) for half a century.

DISRUPTING EVENT

Implementation of improved practices by the industry, from the 1990s.

KEY TOPICS

• Ecotoxicology • Dynamics of pollutants, remediation and resilience • Environment, health and population dynamics • Lagoon ecosystem and industrial activity • Social perception of risks related to industrial activity

TERRITORY

The study area is the Estarreja municipality located at Central Region of Portugal, at the north coast of Aveiro District. This region was affected by an intense industrial pollution since the early 50's, because it hosts the second major Portuguese complex of chemical industries, producing mainly ammonium sulphate, nitric acid and ammonium nitrate, but also synthetic resins (PVC). The industrial area is located next to the largest brackish water estuary of Portugal, which extend their channels to the area where is the Chemical Complex of Estarreja. The "Aveiro ria" is an ecologically rich environment, and the natural habitat of numerous species of waterfowl and terrestrial. This industrial activity has, until 1986, produced a large volume of toxic waste solids which were disposed in areas without properly preparation for this purpose. Furthermore, until 1975, liquid waste effluents, with potential toxic elements to the ecosystems and human health, were rejected in several sewage outlets.

SOCIAL-ECOLOGICAL ISSUES

During the 90's, the technologic progresses made in this sector allowed a significant reduction of the pollution from the local chemical industries. Also during this decade, several rehabilitation interventions resulted in an important reduction of the environmental legacy (solid waste and liquid effluents). Despite the major advances, there are still gaps of knowledge which justify the development of multidisciplinary research studies, using innovative approaches.

ÉVÉNEMENT FONDATEUR

Mise en oeuvre de pratiques vertueuses par les industriels à partir de 1990.

THÉMATIQUES CLÉS

PROBLÉMATIQUE SOCIO-ÉCOLOGIQUE

Au cours des années 1990, d'importantes innovations technologiques ont permis de réduire considérablement l'émission de polluants par les usines chimiques locales. De même, au cours de cette décennie, un certain nombre d'interventions de réhabilitation ont permis une réduction significative du passif environnemental, tant au niveau des résidus solides que des effluents liquides. Cependant, de nombreuses lacunes subsistent dans la connaissance de l'évolution de l'environnement ; elles justifient la nouvelle approche pluridisciplinaire qui est celle de l'OHMi Estarreja. L'OHMi Estarreja a pour objectif d'étudier les effets multiples de l'activité industrielle sur la relation homme-milieu, en prenant en compte la dimension temporelle, marquée par une amélioration significative des processus technologiques industriels, ainsi que par des changements dans la structure sociale de la population et dans les modes de vie. La perception et l'implication des populations locales est également pour l'OHMi un enjeu primordial.

  The OHMi Estarreja aims to study the effects of industrial activity humanenvironment relationship, taking into account the temporal dimension, marked by a significant change in industrial technological processes, as well as in the social structure and life habits of the population. The perception and incorporation of the local populations is very relevant to this Observatory. Observatoire Hommes-Milieux International (OHMi) Estarreja District d'Aveiro, Portugal CRÉATION Mars 2010 CADRE SOCIO-ÉCOLOGIQUE (OU « FAIT STRUCTURANT ») Développement d'un complexe d'industries chimiques (sulfate d'ammonium, acide nitrique, nitrate d'ammonium, résines synthétiques) ayant généré une intense pollution du milieu depuis les années cinquantes et pendant près d'un demi-siècle, à proximité de la lagune d'Aveiro.

  polluants, résilience et remédiation • Environnement, santé et dynamique des populations • Ecosystème lagunaire et activité industrielle • Perception sociale des risques liés à l'activité industrielle TERRITOIRE Le territoire étudié correspond à la commune d'Estarreja, située dans la Région Centre du Portugal (district d'Aveiro). La commune a été soumise à une pollution industrielle intense depuis le début des années 50, en raison de la présence du deuxième plus important complexe d'industries chimiques du pays, produisant principalement du sulfate d'ammonium, de l'acide nitrique et du nitrate d'ammonium, mais également des résines synthétiques (PVC). La zone industrielle se situe à proximité de la plus grande lagune d'eau saumâtre du Portugal, milieu écologiquement riche, habitat naturel de nombreuses espèces aquatiques et terrestres. L'activité industrielle a produit de grandes quantités de déchets toxiques solides, déchets qui ont été stockés, jusqu'en 1986, dans des lieux inadaptés. De plus, jusqu'en 1975, des effluents liquides chargés en éléments potentiellement toxiques pour les écosystèmes et pour la santé humaine étaient également rejetés dans les nombreux canaux de drainage de cette zone.


Luca Arpaia
email: l.arpaia@brgm.fr

Mario Ricchiuto
email: mario.ricchiuto@inria.fr

Andrea Gilberto Filippini
email: a.filippini@brgm.fr

Rodrigo Pedreros
email: r.pedreros@brgm.fr

An e cient covariant frame for the spherical shallow water equations: well balanced DG approximation and application to tsunami and storm surge

Keywords: shallow water equations, spherical geometry, discontinuous galerkin, well balanced schemes, tsunami, storm surge

In this work we consider an e cient discretization of the Shallow Water Equations in spherical geometry for oceanographic applications. Instead of the classical 2d-covariant or 3d-Cartesian approaches, we focus on the mixed 3d/2d form of [Bernard et al., JCP 2009] which evolves the 2d momentum tangential to the sphere by projecting the 3d-Cartesian right-hand side on the sphere surface. Di↵erently from the last reference we consider the exact representation of the sphere instead of the finite element one, mixed with a covariant basis projection of the momentum. This leads to several simplifications of the Discontinuous Galerkin scheme: the local mass matrix goes back to the standard block-diagonal form; the Riemann Problem does not require any tensor or vector rotations to align the bases on the two sides of an edge.

Second we consider well balancing corrections related to relevant equilibrium states for tsunami and storm surge simulations. These corrections allow to compensate for the inherent non-exactness of the quadrature induced by the non-polynomial nature of both the geometrical mapping and of the covariant basis. In other words, these corrections are the order of the cubature error. We show that their addition makes the scheme exactly well balanced, and is is equivalent to recasting the integral of the hydrostatic pressure term in strong form. The proposed method is validated on academic benchmarks involving both smooth and discontinuous solutions, and applied to realistic tsunami and an historical storm surge simulation.

Introduction

Coastal flooding is mainly caused by meteorological and telluric phenomena which may occur at the scale of the oceanic basin. The numerical simulation of such events requires to solve the governing equations in a spherical geometry. In this work we consider the Shallow Water Equations (SWEs) which provide satisfactory results for long wave phenomena such as tide/storm surge and tsunami.

The literature devoted to discrete formulation of the SWEs on the sphere is quite large. Two methods have emerged in the past to deal with the Earth curvature. A common approach among contemporary ocean models consists in a two-dimensional parametrization of the sphere through a proper curvilinear coordinate system. Then the primitive or shallow water equations are written directly on the curved manifold meaning that all the di↵erential operators are transformed in curvilinear coordinates. A metric source term arises from the space variation of the curvilinear vectors basis. The underlying equations in curvilinear coordinates are typically discretized on regular grids with finite di↵erences/finite volume schemes. Regional models uses the traditional latitude-longitude parametrization. Unfortunately for global ocean models, the latitude-longitude coordinate system has a singular point in the Arctic Ocean at the North Pole, where the meridians converge on water. The Jacobian of the coordinate transformation is not defined and moreover the singularity imposes a severe restriction on the maximum time step allowed for stability. Nowadays the pole problem is considered merely as an economical one [START_REF] Williamson | The evolution of dynamical cores for global atmospheric models[END_REF]. In particular, meteorologists have used a variety of fixes to circumvent the singularity and all of such fixes employ multiple curvilinear coordinates defined on di↵erent patches covering the whole sphere, [START_REF] Kageyama | Yin-Yang grid": An overset grid in spherical geometry[END_REF][START_REF] Ronchi | The "cubed sphere": A new method for the solution of partial di↵erential equations in spherical geometry[END_REF]. Then a specific treatment of the flux term at the edges dividing the patches, where the tangent basis changes, is deployed: a rotation of momentum vector [START_REF] Nair | A discontinuous galerkin transport scheme on the cubed sphere[END_REF] or an orthonormalization of the Riemann Problem [START_REF] Ullrich | High-order finite-volume methods for the shallow water equations on the sphere[END_REF][START_REF] Rossmanith | A wave propagation method for hyperbolic systems on the sphere[END_REF] assure a common reference system at the edge and local conservation. To bypass the pole problem ocean modelers have used instead ad hoc curvilinear coordinates conceived to place the singularity point on lands, see for example the tripolar grid of [START_REF] Madec | A global ocean mesh to overcome the north pole singularity[END_REF]. Finally, for Finite Volume methods, some authors [START_REF] Rossmanith | A wave propagation method for hyperbolic systems on the sphere[END_REF] claimed that extra care should be taken in the resolution of the Riemann Problem in curvilinear coordinates.

A second approach consists instead in resolving the governing PDEs in a three-dimensional Cartesian framework and then adding a constraint to force the currents to remain tangent to the sphere [START_REF] Coté | A Lagrange multiplier approach for the metric terms of semi-Lagrangian models on the sphere[END_REF] [START_REF] Giraldo | Nodal high-order discontinuous galerkin methods for the spherical shallow water equations[END_REF]. With this approach, it is more common to represent the sphere approximately through a local finite element map for each element. Thanks to the local map, no special treatment is required in polar regions to preserve accuracy and to conserve global mass.

A third approach seems to combine the advantages of both methods [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF]: momentum time derivative is written with respect to 2d components while the right-hand side is first expressed in 3d and then projected 2 back onto the sphere surface by a simple scalar product with respect to the tangent basis. The advantage is that the number of unknown is kept to a minimum (water depth and two momentum components) and, at the same time, the right-hand side maintains a Cartesian form, thus it is independent from the parametrization of the sphere, and moreover there is no need to transform the di↵erential operators. For example, in tsunami applications this could be laborious for depth averaged non-hydrostatic models with dispersive terms that involve mixed high order derivatives. If requested, Riemann solvers are formulated easily in 3d Cartesian framework, and then projected on the sphere surface along with the right-hand side.

In the last reference the geometrical error related to the approximation of the curvature of the sphere is controlled by using higher order super-parametric finite element approximations of the spherical geometry.

Here we present some improvements of the mixed 3d/2d formulation of [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF] applied to spherical geometry. First we use an exact representation of the spherical elements, which allows to avoid the pole problem altogether. A similar idea is used in the recent work by [START_REF] Baldauf | Discontinuous galerkin solver for the shallow-water equations in covariant form on the sphere and the ellipsoid[END_REF]. In our case, we then introduce a 2d covariant representation of the momentum vector. This combination completely cancels the geometrical error related to the sphere curvature. In practice, we can use any local reference mapping we want to represent the 3d spatial coordinates, including the classical latitude longitude transformation. Due to the orthogonality of covariant vector basis, the momentum equations remain decoupled as in the planar case, which simplifies the structure of the update. Moreover, the Riemann solvers also simplify thanks to continuity of the geometrical quantities as well as of the covariant basis at the elements edges. More precisely, no rotations of the velocities/momentum components, nor averaging of the normals are required.

The discontinuous Galerkin (DG) method is becoming increasingly popular in geophysical fluid dynamics, proving also in this field its potential in terms of high order accuracy, conservation properties and ease of massively parallel implementation. It has been applied to the spherical SWEs in di↵erent forms: at the opposite ends of the spectrum we mention the purely two-dimensional form on the cubed sphere [START_REF] Nair | A discontinuous galerkin transport scheme on the cubed sphere[END_REF] and the fully 3d Cartesian form on a finite elements quadrilateral grids [START_REF] Giraldo | Nodal high-order discontinuous galerkin methods for the spherical shallow water equations[END_REF]. Here we discretize the spherical SWEs in the proposed 3d/2d-covariant form. Then we look to the well balanced property, or the discrete preservation of some steady states which are of particular interest for tsunami and storm surge computations. These states include the lake at rest, which results from the steady balance between hydrostatic pressure and the force played by the bottom and the inverse barometer in which the atmospheric pressure e↵ect is also added to the balance. In general well balanced methods are more robust and allow to work with coarser mesh. A well balanced DG method typically relies on the fact that the so-called strong form (with a second integration by parts of the volume flux term) embeds more easily then the weak form the well balanced property at a discrete level. Based on the fact that the strong form is more suited to retain well balancing, some authors ( [START_REF] Ern | A well-balanced runge-kutta discontinuous galerkin method for the shallowwater equations with flooding and drying[END_REF][START_REF] Xing | Positivity-preserving well-balanced discontinuous galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF]), while being stuck to the weak DG form (with one integration by parts of the flux term), compute exactly the integrals in order to maintain a discrete equivalence to the strong form. Other authors (see e.g. [START_REF] Vater | A limiter-based well-balanced discontinuous galerkin method for shallow-water flows with wetting and drying: Triangular grids[END_REF]) write directly the scheme in the strong form. This has close similarities with fully residual based methods (see e.g. [START_REF] Abgrall | High order methods for CFD[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF]), and allows to obtain well balancing independently on the quadrature formula used. Such flexibility is used to retain general well balanced states on the sphere [START_REF] Bonev | Discontinuous galerkin scheme for the shallow water equations with application to tsunami modeling and prediction[END_REF]. Here we start from the weak DG form and we show that, when applied to the mixed 3d/2d formulation of the SWEs, it is unable to preserve equilibria due to inexactness of the quadrature formula with non-polynomial functions as the ones contained in the components of the tangent vectors. We add a well balancing correction related to relevant equilibrium states for tsunami and storm surge simulations. These corrections involve the quadrature error. There are thus small by definition, and compensate for the lack of summation by parts property making the scheme exactly well balanced. We show that this idea is equivalent to recasting the integral of the hydrostatic pressure term in the strong form.

The paper is organized as follows. In section 2 we introduce the notation for the spherical geometry.

In section 3 we introduce the SWEs in the 3d-2d covariant formulation. In section 4 we discuss their well balanced DG approximation. In section 5 we detail computational aspects such as shock-capturing and time-stepping. In section 6 we evaluate the proposed methods on standard global benchmarks for the shallow water equations on the sphere and on regional realistic tsunami/storm surge computations.

Tangent maps for the sphere

We consider a sphere S 2 with radius R described by generic curvilinear coordinates {X 1 , X 2 } and tangent basis {t 1 , t 2 }. The coordinate vector x then writes:

x = x i e i ⌘ x 1 e 1 + x 2 e 2 + x 3 e 3 (1) = X ↵ t ↵ ⌘ X 1 t 1 + X 2 t 2 (2)
Note that in the above definitions and in the following ones repeated indices imply summation. In particular, Latin indices are used for the components of vectors in Cartesian coordinates (from 1 to 3), while Greek indices are used for vectors in local curvilinear coordinates (from 1 to 2). The transformation from the curvilinear space to the physical space is defined here as x = T (X). It has a Jacobian or transformation tensor:

J T = @x @X (3)
with components J i↵ T and an inverse Jacobian:

J 1 T = @X @x (4)
with components (J 1 T ) ↵i . The tangent vectors can be obtained by di↵erentiation as the columns of the Jacobian. It is useful to normalize the basis:

t ↵ = @x @X ↵ , t ⇤ ↵ = t ↵ ||t ↵ || (5) in general t ⇤ ↵ • t ⇤ 6 = ↵ .
Finally we define the metric tensor constructed from the Jacobian matrix as:

G = J t T J T (6)
and G = det G, having denoted by t the matrix transpose. If the inverse Jacobian cannot be computed directly from the definition (4), it can be replaced by the pseudo-inverse:

J 1 T = G 1 J t T (7)
Depending on the transformation, di↵erent representations of the spherical geometry can be obtained.

Finite Element map

T can be a local finite element transformation [START_REF] Giraldo | Nodal high-order discontinuous galerkin methods for the spherical shallow water equations[END_REF][START_REF] Bonev | Discontinuous galerkin scheme for the shallow water equations with application to tsunami modeling and prediction[END_REF] that maps a reference element to a spherical element, denoted by

E K : K ! K, x = X j2DOF s ' j (X)x j (8)
where X 1 , X 2 are the reference coordinates and x j are the physical coordinates of the degrees of freedom.

The Jacobian and the tangent vectors follow from the application of the definitions (3) and ( 5), while the inverse Jacobian (4) has to be intended as the pseudo-inverse. A high order accurate geometry is usually attained by increasing the polynomial degree of the basis functions ' j .

Spherical map

In this work we will use the exact spherical transformation, defined here as:

G : [ ⇡, +⇡] ⇥ h ⇡ 2 , + ⇡ 2 i ! S 2 , x = G(X) (9)
where

x 1 = R cos X 2 cos X 1 x 2 = R cos X 2 sin X 1 (10)
x 3 = R sin X 2
with X 1 the longitude and X 2 the latitude. Explicit expressions for the Jacobian, inverse Jacobian and tangent vectors follow from (3), (4), and (5). For example the Jacobian reads:

J G = 2 6 6 6 4 R cos X 2 sin X 1 R sin X 2 cos X 1 R cos X 2 cos X 1 R sin X 2 sin X 1 0 R cos X 2 3 7 7 7 5 (11)
The tangent vectors are also called covariant vectors and they are denoted as g ↵ = @x @X ↵ . Please note that the normalized covariant basis is orthonormal g ⇤ ↵ • g ⇤ = ↵ . As already mentioned, the spherical transformation is singular at the poles.

An alternative and equivalent exact map that avoids the pole singularity but does not keep orthogonality of the tangent basis is the following local transformation that maps a reference element to a spherical element:

A K : K ! K, x = R F (X) ||F (X)|| (12)
where F = P 3 j=1 ' j (X)x j is the P1 finite element transformation from the reference element to the 3d element. The above transformation has a smooth Jacobian (3) everywhere:

@x k @X ↵ = R ||F || 3 X j=1 @' j @X ↵ x k j R P 3 j=1 ' j (X)x k j ||F || 2 @||F || @X ↵ (13)
from which we can compute the normalized tangent vectors with (5).

Shallow Water Equations in spherical geometry

We consider the SWEs in spherical geometry with uneven bathymetry and in presence of earth rotation.

We express them in a generic vector form:

@h @t + r • (hu) = 0 (14) @hu @t + r • T = S (15)
with h the fluid depth, hu the discharge/momentum vector, and having denoted the momentum flux tensor by T . The source S includes the e↵ects of bathymetry, bottom friction, Coriolis force, and meteorological forcing:

S = ghrb c F hu 2⌦k ⇥ hu gh ⇢ 0 rp atm + f w (16)
In the formula above b(x) denotes the bathymetry/topography, c F the friction coe cient given by the Manning formula, ⌦ the Earth rotation rate and k the Earth rotation axis. p atm is the atmospheric pressure, ⇢ 0 is the water density and f w represents the wind forcing. Turbulent eddy viscosity term has been neglected.

The expressions of nabla r•, of the momentum vector hu and of the tensor T depend on the coordinate system in which they are expressed. The most common representations used for the shallow water equations are either the 3d Cartesian coordinates or the curvilinear coordinates defined in section 2. The momentum vector hu can be expressed in both systems (see ( 1) and ( 2)) as:

hu = hu i e i ⌘ hu 1 e 1 + hu 2 e 2 + hu 3 e 3 ( 17
) = hu ↵ t ↵ ⌘ hu ↵1 t 1 + hu ↵2 t 2 (18)
We remark that, for curvilinear coordinates, momentum components are not defined in a unitary basis.

It may be convenient expressing momentum in a unitary basis:

hu = hu ↵ ||t ↵ || t ↵ ||t ↵ || (19) = hu ⇤↵ t ⇤ ↵ (20)
The relationship between the velocity components written in the two basis is:

u i = J i↵ T u ↵ , u i = J ⇤,i↵ T u ⇤↵ (21)
where J ⇤,i↵ T = J i↵ T /||t ↵ || (repeated index summation does not apply here) is the Jacobian (3) normalized by columns. Similarly the second order flux tensor T can be expressed in both coordinate systems as:

T = T ij e i e j (22)
= T ↵ t ↵ t (23)
with Cartesian components T ij = hu i u j + P ij and curvilinear components T ↵ = hu ↵ u + GP ↵ . The hydrostatic pressure is defined as P = gh 2 /2 and G is the determinant of the metric tensor (6).

Equilibrium states

For coastal flows, an important role is played by the so-called lake at rest state which, denoting the free surface level ⌘ = h + b, is the particular steady solution characterized by the two invariants [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF]:

hu = 0, Z = ⌘ = const = Z 0 (24)
In presence of atmospheric pressure forcing, a relevant state is the inverted barometer balance that is an exact solution of the SWEs in case of full adjustment of sea level to changes in barometric pressure [START_REF] Gill | Atmosphere-Ocean Dynamics[END_REF]:

hu = 0, Z = ⌘ + patm g⇢0 = const = Z 1 ( 25
)
Full adjustment is never reached since the pressure field is unsteady p atm (x, t), thus one may question if this is really a relevant equilibrium, but it is interesting to be able to handle this possibility numerically.

In [START_REF] Mercer | Barotropic waves generated by storms moving rapidly over shallow water[END_REF] it has been shown that, at least in deep water, the free surface response is almost isostatic and characterized by small perturbations around the inverted barometer balance. In light of these findings and for completeness, we have investigated the more general balance [START_REF] Pasquetti | Viscous stabilizations for high order approximations of Saint-Venant and Boussinesq flows[END_REF].

A numerical method approximating ( 14)-( 15) is said to be well balanced if [START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF] or [START_REF] Pasquetti | Viscous stabilizations for high order approximations of Saint-Venant and Boussinesq flows[END_REF] are also exact solutions of the discrete equations. In other words, well balanced schemes provide a discrete analog of the balance between hydrostatic pressure and the forces exerted by the bottom (and by atmospheric pressure)

r • P I = ghr ✓ b + p atm g⇢ 0 ◆ (26)

Fully diagonal mixed 3d/2d-covariant formulation

Following [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF] we proceed as follows to obtain the fully discrete equations.

1. we project the momentum equation [START_REF] Abgrall | High order methods for CFD[END_REF], on the tangent plane:

@ @t (hu • t ⇤ ↵ ) + (r • T ) • t ⇤ ↵ = S • t ⇤ ↵ 2.
the scalar product contained in the time derivative is performed with momentum expressed in curvilinear coordinates while the scalar product involving the rest of the equation is performed in 3d Cartesian coordinates:

@ @t (hu ⇤ t ⇤j t ⇤j ↵ ) + @T jk @x k t ⇤j ↵ = S j t ⇤j ↵ (27)
e.g. (↵ = 1) @ @t

(hu ⇤1 t ⇤ 1 • t ⇤ 1 + hu ⇤2 t ⇤ 2 • t ⇤ 1 ) + ✓ @T 11 @x 1 + @T 12 @x 2 + @T 13 @x 3 ◆ t ⇤1 1 aaaaaaaaaaaaaaaaaaaaaaa + ✓ @T 21 @x 1 + @T 22 @x 2 + @T 23 @x 3 ◆ t ⇤2 1 aaaaaaaaaaaaaaaaaaaaaaa + ✓ @T 31 @x 1 + @T 32 @x 2 + @T 33 @x 3 ◆ t ⇤3 1 = S j t ⇤j 1
One advantage of this formulation with respect to full 3d equations is that the number of unknowns is kept minimum (h, hu 1 , hu 2 ). Another attractive feature is that the flux function stays in 3d form and does not depend on a particular transformation. This means that line integrals are defined intrinsically and mass is easily conserved circumventing implementation issue related to the use of composite meshes in the 2d approach. The drawback of this formulation compared to the fully two-dimensional one is that for each momentum equation 9 components of the fluxes have to be evaluated instead of three, and that in general the curvilinear momentum components are coupled.

However, di↵erently from [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF] we propose to use the exact spherical representation, thus we have to project momentum on the exact tangent basis. We can choose either the covariant framework [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF] or the local transformation [START_REF] Ern | A well-balanced runge-kutta discontinuous galerkin method for the shallowwater equations with flooding and drying[END_REF]. Here we have used locally in each element the covariant framework guaranteeing the decoupling of the momentum components, and for which simple analytical expressions exist for the tangent vectors g ↵ = @x @X ↵ . We project onto the normalized tangent basis (5) to get:

@ @t (hu • g ⇤ ↵ ) + (r • T ) • g ⇤ ↵ = S • g ⇤ ↵ ( 28
)
or by components (e.g. ↵ = 1):

@ @t (hu ⇤1 g ⇤ 1 • g ⇤ 1 | {z } =1 +hu ⇤2 g ⇤ 2 • g ⇤ 1 | {z } =0 ) + ✓ @T 11 @x 1 + @T 12 @x 2 + @T 13 @x 3 ◆ g ⇤1 1 aaaaaaaaaaaaaaaaaaaaaaa + ✓ @T 21 @x 1 + @T 22 @x 2 + @T 23 @x 3 ◆ g ⇤2 1 aaaaaaaaaaaaaaaaaaaaaaa + ✓ @T 31 @x 1 + @T 32 @x 2 + @T 33 @x 3 ◆ g ⇤3 1 = S j g ⇤j 1
So in practice we work with:

@ @t (hu ⇤↵ ) + @T jk @x k g ⇤j ↵ = S j g ⇤j ↵ ( 29
)
This choice allows to take advantage of the orthogonality of the covariant basis and work with a fully diagonal system which is the main advantage of the formulation proposed here. Please note also that this algorithm evolves physical components of momentum hu ⇤↵ in m 2 /s (eastward for ↵ = 1, and northward for ↵ = 2).

Well balanced discontinuous finite element discretization

We consider a discretization of S 2 composed by non overlapping triangular elements K. To each element K we can associate an element defined in the parametric space, which we denote by K. For a given mesh, we label by the roman h the characteristic mesh size measured in physical space. Following a high order DG procedure, we project the vector of conserved variables and the bathymetry into the finite dimensional space:

V p h = {' : '(x) 2 L p K , x 2 K} ( 30
)
where L p K are Lagrange polynomials of degree p defined on each parametric element K. We obtain the weak formulation multiplying the mass equation ( 14) and the momentum equation ( 28) by a smooth test function ' i that belongs to the same space [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF] and integrating it on a spherical triangle K:

@ @t Z K h h ' i dx + Z @K hu g h ' i • n ds Z K hu h • r' i dx = 0 (31) @ @t Z K hu h • g ⇤ ↵ ' i dx + Z @K T g h g ⇤ ↵ ' i • n ds Z K T h : r (g ⇤ ↵ ' i ) dx = Z K S h • g ⇤ ↵ ' i dx (32)
where T g h denotes the numerical flux evaluated at the elements boundaries, and the symbol : is the scalar product between second order tensors A :

B = A ij B ji . As noted in [10], if we regroup g ⇤ ↵ ' i = v i,↵ , then (32)
corresponds to consider a vectorial test function for each momentum equation.

Several aspects need to be clarified. First of all the mapping required to define the quadrature formulas, which will involve non-algebraic transformations. Secondly, the well balancedness w.r.t the equilibria presented earlier is not trivial in this formulation as exact quadrature is impossible precisely due to the non-algebraic nature of the mappings adopted. These aspects are covered in the next sections.

Reference element mapping and quadrature

All the integrals are computed on the standard reference triangle K. For the transformation which map a reference triangle to a spherical triangle we use [START_REF] Ern | A well-balanced runge-kutta discontinuous galerkin method for the shallowwater equations with flooding and drying[END_REF], which allows to bypass the pole singularity in the evaluation of the integral and the di↵erential operators. Consider for example the volume integral in the momentum equation which we decompose into

Z K T h : r (g ⇤ ↵ ' i ) dx = Z K T kj h @ @x k g ⇤j ↵ ' i dx = Z K T kj h ' i @g ⇤j ↵ @x k dx + Z K T kj h g ⇤j ↵ @' i @x k dx (33)
Both the resulting integrals are mapped using [START_REF] Ern | A well-balanced runge-kutta discontinuous galerkin method for the shallowwater equations with flooding and drying[END_REF], so the di↵erential is rescaled based on the mapping Jacobian

dx = J A dX. ( 34
)
with [START_REF] Rossmanith | A wave propagation method for hyperbolic systems on the sphere[END_REF]. The derivatives of the test functions are evaluated classically via the gradient in the mapped element as

J A = p G and G = det(J t A J A ), see
@ @x k = @X @x k @ @X = (J 1 A ) k @ @X (35)
with the inverse Jacobian J 1 A = G 1 J t A , see [START_REF] Madec | A global ocean mesh to overcome the north pole singularity[END_REF]. All the remaining quantities are evaluated analytically via the mapping [START_REF] Ern | A well-balanced runge-kutta discontinuous galerkin method for the shallowwater equations with flooding and drying[END_REF], including the gradients of the covariant basis vectors, which can be easily obtained analytically starting from [START_REF] Baldauf | Discontinuous galerkin solver for the shallow-water equations in covariant form on the sphere and the ellipsoid[END_REF] Line integrals are evaluated similarly via a mapping to a reference 1d element. As already remarked, for the the time derivative we can write

@ @t Z K hu h • g ⇤ ↵ ' i dx = X j2DOF s Z K ' i ' j dx @hu ⇤↵ j @t (36)
which shows that a fully block diagonal structure is retained.

As a final remark, we note that di↵erently from two-dimensional formulations, the right-hand side appears in Cartesian form and it is well-posed on the whole sphere. However the tangent basis in [START_REF] Arpaia | Well-balanced residual distibution for the ALE spherical shallow water equations on moving adaptive meshes[END_REF] is not well defined at the poles, for the latitude-longitude parametrization that we have chosen. For this reason, only to evaluate the tangent basis and its derivatives, a polar cap defined by a limiting latitude z lim , is deployed. Elements are flagged according to the baricenter position: if z K G > z lim we use a rotated spherical transformation replacing [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF] with:

x 1 = R cos X 2 cos X 1 x 2 = R sin X 2 (37)
x 3 = R cos X 2 sin X 1
The Jacobian (3), inverse Jacobian (4) and basis vectors (5) are modified accordingly. Note that this is only a↵ecting the local representation of the 3d basis which remains perfectly continuous across element edges. We found that the DG framework is well suited to handle such discontinuity at the elements edges with no modification necessary to the algorithm. In particular:

• each element updates the momentum components hu ⇤↵ in its own reference system, thus zonal/meridional components in standard elements and rotated components in polar cap elements.

• in a DG method, communication between elements is needed only for evaluating numerical fluxes.

Since they are computed in a 3d framework (and independent from the local coordinate system) no specific treatment is necessary at the edges separating the polar cap.

In practice we have set z lim = 0.9R, however we will show that this value has no impact on the numerical results.

Integration by parts correction

We are interested in the well -balanced nature of the scheme. We will first consider the lake at rest case.

The e↵ect of the atmospheric pressure is added introducing an apparent topography [START_REF] Bouchut | Frontal geostrophic adjustment and nonlinear-wave phenomena in one dimensional rotating shallow water[END_REF] b ! B = b+ patm g⇢0 .

Imposing the lake at rest, the momentum components of the DG method reduce to:

Z @K P g h Ig ⇤ ↵ ' i • n ds Z K P h I : r (g ⇤ ↵ ' i ) dx = Z K gh h rb • g ⇤ ↵ ' i dx (38)
One drawback of the mixed formulation emerges: well balancing is lost because of the projection on the tangent basis, an operation that involves the evaluation of non-algebraic functions in all of the three integrals of [START_REF] Knapp | The international best track archive for climate stewardship (IBTrACS): unifying tropical cyclone best track data[END_REF]. Indeed even with constant bathymetry, and constant pressure, in general we cannot remove the geometrical integration error

✏ K P0 = P 0 Z K I : r (g ⇤ ↵ ' i ) dx Z @K Ig ⇤ ↵ ' i • n ds (39)
This is an issue already encoutered on curved finite element (polynomial) meshes [START_REF] Kopriva | Metric identities and the discontinuous spectral element method on curvilinear meshes[END_REF], and is even more critical for the non-polynomial approximation used here. We remark that this was also highlighted in [START_REF] Bonev | Discontinuous galerkin scheme for the shallow water equations with application to tsunami modeling and prediction[END_REF] in the context of the study of the well balancing of a fully 3d DG shallow water approximation on the sphere. As observed in the last reference, a more suitable form to guarantee the well balanced character of the method is the strong form. We obtain a somewhat similar result here, but we proceed in a di↵erent manner, and finally only the pressure term is converted into the strong form. Our idea is to balance the scheme with a small correction which is within the quadrature error, and thus much smaller than the truncation error of the scheme. To achieve this we re-inject into the scheme (32) the integration by parts error associated to the pressure terms, given by

✏ K P = Z K P h I : r (g ⇤ ↵ ' i ) dx Z @K P h Ig ⇤ ↵ ' i • n ds + Z K (r • P h I) • g ⇤ ↵ ' i dx (40)
Note that this term is identically zero for exact quadrature, and reduces to [START_REF] Holland | An analytic model of the wind and pressure profiles in hurricanes[END_REF] if the pressure is constant. Note also that the only anzatz made is that P h is the pressure associated to the finite element polynomial projection, and in our specific case P h it is assumed to be a polynomial of degree 2p, thus infinitely di↵erentiable. We can characterize the correction term with the following property.

Proposition (IBP correction bounds) Let (40) be evaluated by means of quadrature formulas which are exact for polynomials of degree k K for the volume integrals and k @K for the boundary integral. The correction can be bounded as

k✏ K P k  C K h min(k @K +1,k K +1) (41)
The proof of the above property is reported in the appendix. In practice, we have computed all integrals with formulas exact for k K = k @K = 3p. A numerical verification of the proposition for these choices is also reported in the appendix. When ( 40) is added to the discrete momentum equation, the variational form becomes

@ @t Z K h h ' i dx + Z @K hu g h ' i • n ds Z K hu h • r' i dx = 0 (42) @ @t Z K hu h • g ⇤ ↵ ' i dx + Z @K T g h g ⇤ ↵ ' i • n ds Z K T h : r (g ⇤ ↵ ' i ) dx = Z K S h • g ⇤ ↵ ' i dx + ✏ K P ( 43
)
4.3. Numerical fluxes: mass conservation and well balanced with discontinuous data

For the mixed 3d/2d formulation, all geometrical quantities are uniquely and objectively defined in 3d, so we can use practically any Riemann solver in 3d Cartesian coordinates. The line integral in the momentum equation ( 32) reads:

Z @K T g h • g ⇤ ↵ ' i • n ds = Z @K T g h • n | {z } H •g ⇤ ↵ ' i ds (44)
where H 2 R 3 is a specified numerical flux H(u L , u R ) and u L and u R are the left (internal) and the right (external) states u = {h, hu i }. In practice we use the Lax-Friedrich flux.

We remark that the problem of using a common framework between the left and right state in the original [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF] is no longer necessary since the normals are continuous at the edges. Only a change of basis is needed at this point to compute 3d velocities because they are only available from momentum update, in 2d. This is done with the right relationship in [START_REF] Bouchut | Frontal geostrophic adjustment and nonlinear-wave phenomena in one dimensional rotating shallow water[END_REF] and the definition of the sphere Jacobian (3). The important consequence of this is that the method proposed conserves mass exactly without any special treatment of the normals on element boundaries.

To recover also the well balanced property in presence of discontinuous data, besides the addition of the hydrostatic correction in [START_REF] Michel | Spectral analysis of continuous FEM for hyperbolic PDEs: influence of approximation, stabilization, and time-stepping[END_REF], as in [START_REF] Xing | Positivity-preserving well-balanced discontinuous galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF] we evaluate the numerical flux at a reconstructed state u >

[23] and add a high order correction term. The resulting well balanced mass/momentum numerical flux (distinguished by •) reads:

H = H(u > L , u > R ) (45)
Ĥ = H(u > L , u > R ) + P (h L ) P (h > L ) n (46)
for the definition of the reconstructed state u = {h > , h > u ⇤i } please refer to [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. A strategy for recovering a well balanced scheme in presence of discontinuous data for the bathymetry on non-conforming elements is also discussed in [START_REF] Bonev | Discontinuous galerkin scheme for the shallow water equations with application to tsunami modeling and prediction[END_REF]. In real applications, discontinuous data b, p atm arise, for example, when the correspondent fields are computed at the degrees of freedom from a local minimization problem. For the bathymetry this can be beneficial to avoid Gibbs phonomena when data is nearly discontinuous on one cell and the polynomial order is p > 1. For the atmospheric pressure field we assume that it is continuously interpolated at the degrees of freedom. In this case the well balanced correction related to the inverted barometer (first line is mass correction, second line is momentum correction):

+ 1 2 1 g⇢ 0 0 @ ↵ LF g h ⇤ R +h ⇤ L 2 n 1 A (p atm,L p atm,R ) (47)
is identically zero, being p atm,L = p atm,R .

Proposition. (Well balanced) The mixed 3d/2d-covariant DG method (42)(43) with numerical fluxes (46) is well balanced in the sense that it preserves to machine accuracy initial data verifying [START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF] or [START_REF] Pasquetti | Viscous stabilizations for high order approximations of Saint-Venant and Boussinesq flows[END_REF].

Proof. The main assumption is that either

h h +b h = ⌘ 0 =const, or h h +b h +(p atm /(g⇢ 0 )) h = ⌘ 0 =const,
within each cell. For the lake at rest state the flux reduces to T h = P h I. Then the discrete momentum equation ( 43) reduces to:

Z @K ⇣ Ĥ P h I • n ⌘ • g ⇤ ↵ ' i ds + Z K (r • P h I) • g ⇤ ↵ ' i dx = Z K gh h rb • g ⇤ ↵ ' i dx (48)
where we have canceled the first term of the correction [START_REF] Saha | The NCEP climate forecast system reanalysis[END_REF] with the third term in [START_REF] Michel | Spectral analysis of continuous FEM for hyperbolic PDEs: influence of approximation, stabilization, and time-stepping[END_REF] since they are equal and with opposite sign. Similarly to [START_REF] Xing | Positivity-preserving well-balanced discontinuous galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF], the proof is based on the fact that, in the lake at rest case, the Lax-Friedrich flux evaluated at the reconstructed state is

H(u > L , u > R ) = P h (h > L )
and one can check that the well balanced numerical flux thus simplifies as Ĥ(u > L , u > R ) = P h (h L )I • n. Since the left state is the internal state we can further write P h (h L )I • n = P h I • n. Then the integrand of the first term in 48 cancels out at each quadrature point and the whole penalization term vanishes. Moreover, for any order

h h , b h 2 V p
h on can show that:

r • P h I + gh h rb h = r g( P j ' j h j ) 2 2 + gh h rb h = gh h r⌘ h (49)

Computational details

Shock capturing

The applications considered in this work may involve the formation of discontinuous solutions during the computation. For tsunami simulations, wave breaking occurring in shallow areas can cause the linear DG method to develop unphysical oscillations that eventually grow in time. We cope with the possibility of non-smooth solutions with an entropy viscosity approach [START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF][START_REF] Pasquetti | Viscous stabilizations for high order approximations of Saint-Venant and Boussinesq flows[END_REF] which consists in adding a di↵usionlike term to the right-hand side of ( 14), [START_REF] Abgrall | High order methods for CFD[END_REF]. Here we discuss the details related to the implementation of the di↵usive term in a spherical geometry. We consider only momentum equation since mass equation is a particular case. Given a di↵usive flux tensor G = µrhu, the DG discretization in the mixed formulation of section 3.2 writes:

I µ := Z K r • G h • g ⇤ ↵ ' i dx = Z @K G h • n • ' i g ⇤ ↵ ds Z K G h : r (g ⇤ ↵ ' i ) dx (50)
To stabilize the shocks only the volume integral is retained and moreover we neglect the variation of the tangent basis so that integral (50) simplifies to:

I µ = Z K G kj h g ⇤j ↵ @' i @x k dx (51)
Note that G kj g ⇤j ↵ = G ↵k is a di↵usive flux tensor projected on the tangent plane. As for the advective flux, since we evolve momentum in the tangent basis, we do not have access directly to G kj in Cartesian basis. However the di↵usive tensor components are easily recovered from relationships ( 21) and ( 35):

G kj = µ @hu j @x k = µ @ @x k ⇣ J ⇤,j G hu ⇤ ⌘ = µ hu ⇤ @g ⇤j @x k + µ J ⇤,j G @hu ⇤ @x k (52)
Near the discontinuities the viscosity coe cient µ should be of order h K in order to kill oscillations while in smooth regions it is expected to be below the approximation order. To this end, we follow the principles originally proposed e.g. in [START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF][START_REF] Pasquetti | Viscous stabilizations for high order approximations of Saint-Venant and Boussinesq flows[END_REF] and define this coe cient starting from the total energy residual. We refer to the last references for the formulation, and to [START_REF] Filippini | UHAINA: A parallel high performance unstructured near-shore wave model[END_REF] and references therein for its parametrization.

Time stepping

High order temporal accuracy is achieved, as it is standard, with a SSPRK method [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF]. The time step is computed according to:

t = CF L min K ⇢ K 2 max f ↵ LF f ( 53
)
with ⇢ K the elemental inscribed circle radius and ↵ LF f the Lax-Friedrich parameter at face f . For linear stability reason one should use CF L = 1 2p+1 (being p the polynomial order of the approximation), but for practical use in this work we always set CF L = 0.5. For stability reasons, when included the Manning friction term is treated via a semi-implicit splitting as e.g. in [START_REF] Hou | A stable 2d unstructured shallow flow model for simulations of wetting and drying over rough terrains[END_REF][START_REF] Delis | A robust high-resolution finite volume scheme for the simulation of long waves over complex domains[END_REF].

Results

In all computations presented hereafter we have set, unless otherwise stated, R = 6371.22 ⇥ 10 3 m, g = 9.80616 m 2 s 1 , ⌦ = 7.295 ⇥ 10 5 s 1 and ⇢ 0 = 1025 kg/m 3 . Simulations have been carried out with linear (P 1 ), quadratic (P 2 ) and cubic polynomials (P 3 ).

Global atmospheric tests

In this paragraph we test the DG scheme on some classical benchmarks taken from the test suite of Williamson et al. (W92) [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF] and [START_REF] Galewsky | An initial-value problem for testing numerical models of the global shallow-water equations[END_REF].

Global Steady-State Nonlinear Zonal Geostrophic Flow. W92 Case #2 is an exact steady geostrophic equilibrium allowing to measure the order of accuracy in presence of Earth rotation. The velocity and height fields are initially given by:

h(', ) = h 0 1 g ✓ ⌦R u 0 + u 2 0 2 ◆ ( cos cos ' sin ↵ + sin cos ↵) 2 u ⇤1 (', ) = u 0 (cos cos ↵ + cos ' sin sin ↵) (54)
u ⇤2 (', ) = u 0 sin ' sin ↵ with gh 0 = 2.94⇥10 4 m 2 s 2 and u 0 = 2⇡R 12days and ↵ = 0. Grid convergence studies gave been conducted on 5 icosahedral grids, from level 1 (h K = 3571 km) to level 5 (h K = 223 km). Relative errors are evaluated following [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF]:

e p = ||e|| L p ||h ex || L p (55)
with e = h h ex and h ex the exact solution. The errors obtained at day 5 are reported on table 1 below. Convergence curves are compared in figure 1 against the results of the original [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF]. To make the comparison e↵ective, only for this experiment, we have tried to stick to their same scheme, DG in the weak-form [START_REF] Arpaia | Well-balanced residual distibution for the ALE spherical shallow water equations on moving adaptive meshes[END_REF] with Lax-Friedrich flux. We have remarked that the use of the well balanced momentum form [START_REF] Michel | Spectral analysis of continuous FEM for hyperbolic PDEs: influence of approximation, stabilization, and time-stepping[END_REF] has almost no impact on the errors. The theoretical order of convergence has been obtained, with smaller errors with respect to the reference. We explain this improvement by the cancellation of the geometrical error in the sphere representation. As it is customary we report on figure 2

with b 0 = 2000 m, R = ⇡/9, r 2 = min(R 2 , (' ' 0 ) 2 + ( 0 )
2 ) and with the center of the cone in (' 0 , 0 ) = (3/2⇡, ⇡/6). Simulations are performed on the fourth icosahedral grid (h K = 446 km). Figure 3 provides height contour lines for P 1 ,P 2 and P 3 computations after 15 days of physical time. For this case one typically checks the conservation of total mass and energy

E = h ✓ 1 2 ||u|| 2 + g h 2 + gb ◆ (57)
Mass error should be only related to round-o↵ errors. Regarding the energy, the proposed scheme have no provable energy conservation properties, which makes the definition of the actual discrete evaluation of the energy somewhat arbitrary. In practice, here the energy is computed directly from nodal quantities, and then integrated using the high order finite element approximation. After 15 days mass is conserved within a relative error of 1e 13. Relative energy conservation error is reported on the right column in figure 3: it is comparable to or better than those typically reported in literature. In the same figure we compare with a second order residual distribution scheme proposed by the authors [START_REF] Arpaia | Well-balanced residual distibution for the ALE spherical shallow water equations on moving adaptive meshes[END_REF], with a third order FV scheme on a cubed sphere [START_REF] Ullrich | High-order finite-volume methods for the shallow water equations on the sphere[END_REF] and with the high resolution results by the German Weather Service (DWD) http://icon.enes.org/swm/stswm/node5.html.

Unstable Jet. The last global test taken from [START_REF] Galewsky | An initial-value problem for testing numerical models of the global shallow-water equations[END_REF] consists of a geostrophically balanced mid-latitude jet, to which a small perturbation is added to initiate the instability. We use this test to confirm that the polar treatment has no impact on the numerical solution and on its derivatives. We look to the vorticity field for P 1 , P 2 and P 3 solutions but with two di↵erent z lim , the one used in other experiments z lim = 0.9R and z lim = 0.75R for which the polar cap interface falls exactly into the jet. Mesh convergence has been carried on two successive icosahedral grids (level 5 and 6, respectively 223 km and 111 km), see figures 4, 5. At both mesh level, P 1 solution is insu cient to support the vortex evolutions. Increasing the order, P 2 , P 3 computations at mesh level 5 provide the correct dynamics but the vorticity contours are a little bit jagged. At level 6 much of the jaggedness disappears.

Small amplitude barotropic waves

With this test we look to the impact of the well balancing correction for classical oceanic applications such as tidal and surge waves. We have taken a basin spanning from 0 to 25 E and from 35 N to 55 N .

The bathymetry is a smooth Gaussian hump mimicking a continental shelf and an isolated mountain:

b = max b 0 exp{ r 2 0 /R}, b 1 (1 r 1 /R) with b 0 = 2000 m, R 0 = ⇡/12, r 2 0 = 8.5(' ' 0 ) 2 + ( 0 ) 2 , (' 0 , 0 ) = (0, ⇡/4) and b 1 = 4000 m, R 1 = 5⇡/180, r 2 1 = (' ' 1 ) 2 + ( 1 )
2 , (' 1 , 1 ) = (10⇡/180, 48⇡/180). The mesh is an unstructured triangulation with topology shown in the top of figure [START_REF] Rossmanith | A wave propagation method for hyperbolic systems on the sphere[END_REF] and constant mesh size h K = 100 km. For a z lim R lake at rest state with initial level ⌘ 0 = 5000 m the well balanced spherical formulation (43) exhibits a maximum perturbation of the velocity norm in the range 1e 12 to 1e 13 depending on the polynomial order. After this well balancing test, a small amplitude tidal signal ⌘ 0 + A sin( 2⇡ T t) is imposed at the north boundary (A = 0.01 m, T = 12.42 hrs). The eastern boundary is a wall to let a southward propagating barotropic Kelvin wave develop. We compare [START_REF] Michel | Spectral analysis of continuous FEM for hyperbolic PDEs: influence of approximation, stabilization, and time-stepping[END_REF] against the non well balanced scheme based on the weak formulation [START_REF] Arpaia | Well-balanced residual distibution for the ALE spherical shallow water equations on moving adaptive meshes[END_REF]. In the left block of figure [START_REF] Rossmanith | A wave propagation method for hyperbolic systems on the sphere[END_REF] two snapshots of upwelling/downwelling at the lateral wall show that the P1 solution of the non well balanced scheme is a↵ected by small wiggles. The perturbation tend to disappear, as expected, with mesh resolution (see right block) and with polynomial order (not shown). A second test has been carried by forcing a surge wave with a small amplitude atmospheric depression:

p atm (', , t) = p 0 p exp ✓ r(', , t) 2 ◆
with r the distance from the center of the storm which is translating southward at constant speed U . The parameters are set to p 0 = 101 kP a, p = 0.3 kP a, = 350 km and U = p 100g. In figure [START_REF] Madec | A global ocean mesh to overcome the north pole singularity[END_REF] we show the impact of the well balanced correction. Again, at the coarse resolution (left block), the small-amplitude gravity waves that propagates faster then the storm are much wiggled in the non well balanced case.

Increasing mesh resolution (see right block) and polynomial order (not shown) makes again such wiggles negligible.

Circular hump centered at the equator

This is a test proposed in [START_REF] Rossmanith | A wave propagation method for hyperbolic systems on the sphere[END_REF]. The initial condition is a circular depth disturbance at the equator:

h(', ) = 8 < : 2 if arccos(cos ' cos )  0.2 0.2 otherwise , u(', ) = 0 (58)
This initial condition is symmetric about the point (', ) = (0 , 0 ), and should remain symmetric in absence of rotation. Shock capturing is necessary to well represent the discontinuity. We run computations on two half sphere meshes generated with gmsh: a coarse one with 2436 points; a fine one with 9187 points.

In figure 8 we compare the solution contours obtained with P1/P2 on the coarse mesh against a P1 solution on the fine mesh. We observe mesh/order convergence for this non-smooth case. At t = 0.9 the P2 solution on the coarse mesh is of the same order of accuracy of the P1 solution on the fine mesh.

We repeat the same test adding rotation e↵ects. We set the dimensionless rotation rate to ⌦ = 5. In figure 9 we report depth contour lines at t = 1.2. As in the previous case, we observe the higher resolution obtained with P2 scheme both in smooth and discontinuous regions. Finally, to show the potential of this method, we perform two realistic simulations in the Caribbean Sea.

The computational domain, the bathymetry and the mesh were specified by the consortium of the project Interreg Carib-Coast (https://www.interreg-caraibes.fr/carib-coast). The domain in figure 10 spans from 50 W to 84.7 W and from 8.4 N to 22.1 N . The bathymetry at the degrees of freedom is computed from GEBCO (www.gebco.net) and from SHOM for the French west indies (https://data.shom.fr/donnees).

Land boundaries are treated as reflective walls. To avoid wetting/drying of cells near the coastline, we have filtered the bathymetry at land masses contained in the computational domain with b filt = min(b, 1).

The filter assures that all the cells are always wet, avoiding positivity preserving issues and well balancing of partly dry cells, all aspects not addressed in the present study. The mesh has variable mesh size ranging from h K,max in deep ocean to h K,min at the coast in correspondence of Martinique, Guadaloupe and Saint Martin islands, which represent some interest sites of the Carib-Coast project. A fine and a coarse mesh have been generated, respectively with h K,min = 300 m, h K,max = 10 km and h K,min = 600 m, h K,max = 20 km. The Manning friction term is set to n = 0.025 for this case.

Tsunami in Carribean Sea

The initial tsunami waveform is associated with a random rupture of the Lesser Antilles subduction zone. We use the fault configuration proposed by Feuillet et al. characterized by a magnitude M w = 8.4

(for the fault parameters see [START_REF] Feuillet | Arc parallel extension and localization of volcanic complexes in Guadeloupe, Lesser Antilles[END_REF]). Randomness is associated to the heterogeneity of the fault slip which is computed by a Karuhnen-Loeve expansion [START_REF] Leveque | Generating random earthquake events for probabilistic tsunami hazard assessment[END_REF]. Then, the initial waveform is the solution of the Okada linear model. We have equipped the domain with several gauges placed in the near and far-field. In the near-field we distinguish coastal gauges in bays or harbours (Point-a-Pitre, Fort-de-France) and nearshore gauges on the continental shelf (Martinique NE, Martinique W). At far-field we have placed the gauges in the Carribean sea (Kingston) and close to east boundary (Atlantic) in order to look both to tsunami propagation after the wave has passed the Antillean Arch and free propagation in deep ocean. In figure 12 we report the sea surface height recorded at the gauge stations for P1/P2 simulations on the coarse mesh and we compare them with the P1 simulation on the fine mesh. The coarse P2 solution has roughly the same accuracy of the fine P1 solution in the near-field and it is superior for far-field propagation. For the near-field gauges the P1 solution captures su ciently well the incoming, trapped and reflected wave train, even if the leading wave peak is slightly underestimated with respect to P2 (40 cm at Point-a-Pitre).

For far-field gauges P1 solution is very di↵usive, probably due to the simple Lax-Friedrich flux employed in this study. Similar results emerge from the snapshot of the sea surface height at t = 5400 s 11. In figure 12, beyond mesh and order convergence we have validated the simulation against the structured high order FV code FUNWAVE-TVD [START_REF] Shi | A high-order adaptive time-stepping TVD solver for boussinesq modeling of breaking waves and coastal inundation[END_REF] on a regular grid with x = 800 m (which has roughly the same degrees of freedom of the P2 coarse simulation). The comparison is valid for the first two hours (after FUNWAVE-TVD solution is polluted by wave reflection from the boundary, see the gauge Atlantic Ocean close the western boundary).

Irma hurricane

We perform the simulation of the atmospheric surge produced by the passage of hurricane Irma in the Caribbean sea (from 30 th August 2017 to 14 th September 2017). As known, such phenomena link up great depressions and strong winds, leading to sea level rise and flooding risk at the coast. In order to account for atmospheric forcing, the source term of the shallow water equations has to be modified as shown in

(3), with the wind force f w defined by:

f w = c D ⇢ air ⇢ 0 |u 10 |u 10 (59)
where ⇢ air is the air density, u 10 is the wind velocity field ten meters above the ground and c D is the wind stress coe cient as proposed by Charnock [START_REF] Charnock | Wind stress on a water surface[END_REF] with Charnock parameter set to 0.022 [START_REF] Muller | Atmospheric storm surge modeling methodology along the french (Atlantic and English Channel) coast[END_REF]. To construct the hurricane pressure and wind fields, we used information issued from the NOAA IBTRACKS [START_REF] Knapp | The international best track archive for climate stewardship (IBTrACS): unifying tropical cyclone best track data[END_REF] database and we applied the Holland parametric model [START_REF] Holland | An analytic model of the wind and pressure profiles in hurricanes[END_REF] Computed hurricane winds are, then, converted to surface winds and averaged over 10 minutes. After that, cyclonic wind and pressure fields are finally merged with data from CFSR model [START_REF] Saha | The NCEP climate forecast system reanalysis[END_REF] in order to represent large-scale atmospheric conditions in the whole computational domain. Our simulations run from 05 th September 2017 to 08 th September 2017, simulating the hurricane strike on Caribbean French islands, Virgin islands and Puerto Rico. Figure 13 shows the hurricane surge propagation, simulated by the P1 coarse mesh. Indistinguishable results, at this scale, are obtained using the P2 coarse mesh, thus they are not reported here. Several NOAA/NOS and SHOM tide-gauges are available in the area, allowing to validate our simulations against real observations. Measured sea level is the result of several components: the tide, a seasonal adjustment due to baroclinicity [START_REF] Pringle | Baroclinic coupling improves depth-integrated modeling of coastal sea level variations around Puerto Rico and the U.S. Virgin Islands[END_REF], the setup due to breaking waves and, finally, the atmospheric surge due to wind and pressure. The first three components cannot be reproduced by our model. Therefore, before comparison, we processed raw signals such as to remove the tide and we applied a mean sea level o↵set to account for baroclinicity [START_REF] Joyce | IOOS coastal and ocean modeling testbed: Hurricane-induced winds, waves, and surge for deep ocean, reef-fringed islands in the Caribbean[END_REF].

In this way, we were able to produce the pictures of figure 14, comparing simulations and observations at four gauges located close to the hurricane track. We can appreciate a global agreement between computed water levels by both P1 and P2 simulations and historical data. The use of P2 approximations, in this case, leads to very few improvements in the results, only noticeable on the pick representation at the Further investigations could improved the results, but this is beyond the purposes of this work and they are not performed here.

Conclusion

In this work we have discussed an e cient and accurate formulation of the spherical Shallow Water equations in the form originally proposed in [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF] which implies a mixed 3d/2d treatment: the time derivative is expressed in 2d form while the right hand side is formulated in 3d Cartesian form. We have shown that using, for the 2d part, covariant vector basis, the implementation simplifies and the accuracy increases. In this work we have opted for a Discontinuous Galerkin discretization but the mixed 3d/2d-covariant form is in principle also more amenable to a continuous Finite Elements discretization w.r.t to the original [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF]. In particular local mass matrix goes back to the standard block-diagonal form, line integrals or Riemann Problem do not imply tensor/vector rotation thanks to the continuity of the normals at the edges. We have also discussed a well balanced form of the resulting DG scheme resorting the hydrostatic pressure term to the strong form. Finally with a realistic tsunami and an historical storm surge simulation we have validated the method for two di↵erent oceanographic applications.

Ongoing work is devoted to extension of the discretization of the 3d/2d-covariant Shallow Water equations to stabilized continuous Finite Element (SUPG or Continuous Interior Penalty see e.g. [START_REF] Abgrall | High order methods for CFD[END_REF][START_REF] Michel | Spectral analysis of continuous FEM for hyperbolic PDEs: influence of approximation, stabilization, and time-stepping[END_REF] and references therein), and to perform a one to one comparison with the current implementation.

We can use now [START_REF] Madec | A global ocean mesh to overcome the north pole singularity[END_REF] to write

@kF k = h kF k 3 X m=1 3 X j=1 F m (X) m j @' j (X)
where note that the derivatives of the linear basis functions are of O(1), in particular, they assume the values 1, 0 or -1 on the standard unit triangle containing the origin of the parametric space. We can write

@x l (X) = kF k 2 @p 1 l (hX) hp 1 l (hX) 3 P m=1 3 P j=1 F m (X) m j @' j (X) kF k 3
All it remains is to note that the derivative of p 1 l (hX) can be written as a 0l h for some O(1) constant a 0l , so that

@x l (X) = h kF k 2 a 0l p 1 l (hX) 3 P l=1 3 P j=1 F l (X) l j @' j (X)
kF k 3 in which the numerator is indeed a quadratic polynomial of the components of F , which combined with (64) leads to @x l (X) = h p 2 l (hX) kF k 3 . We need now to prove that, if true for k 1, the statement is also true for k, so we assume that

@ k 1 x l (X) = h k 1 p k l (hX) kF k 2k 1
and compute its di↵erential:

@ k x l (X) =h k 1 @p k l (hX)kF k 2k 1 p k l (hX)(2k 1)kF k 2k 2 @kF k kF k 4k 2 = h k 1 @p k l (hX)kF k 2 p k l (hX)(2k 1)kF k@kF k kF k 4k 2 (2k 3) =h k 1 @p k l (hX)kF k 2 hp k l (hX)(2k 1) 3 P l=1 3 P j=1 F l (X) l j @' j kF k 2k+1
As before, we remark that @p k l (hX) = hp k 1 l (hX), which allows to write

@ k x l (X) = h k pk 1 l (hX)kF k 2 hp k l (hX)(k + 1) 3 P l=1 3 P j=1 F l (X) l j @' j ) kF k 2k+1
which indeed by virtue of (64) can be written as (65), and concludes the induction.

392

To conclude this preliminary proof, we note that given an upper bound to the mesh size H > h, we can easily find a uniformly bounded constant C H such that kF k is bounded in between R and R C H H and 1/kF k is bounded as well. Using this remark, and the boundedness and the coe cients of p k l , the first of the properties sought is obtained for k = 1, the second for k = m + 1. The k = 1 property allows readily to prove (62) by using the definition

J A = p G , G = det(J t A J A )
Since the determinant of a 2⇥2 matrix is a quadratic form of the matrix entries, this means that G = O(h 4 ) by virtue of (60) and p G = O(h 2 ), which implies (62) by the definition of J A . Lastly the properties (60)-( 61)-( 62) can be used to show that r X J A = O(h 3 ), hence 63. ⇤

To complete the proof of the estimate, we now rewrite each of the integrals in the parametric space, and perform a truncated Taylor series expansion with the appropriate degree in the parametric coordinates to obtain an estimate of the quadrature error. We will show the procedure for one of the integrals of (40), and provide comments concerning possible di↵erences when evaluating others. So, let us consider the evaluation of

I 1 = Z K (P h I) : ' i rg ⇤ ↵ dx (66)
As we said in the main body of the paper g ⇤ ↵ is available in a uniformly bounded analytical form, for example in latitude-longitude coordinates. Note that, with the usual assumption of smoothness, removing from the variational form its analog evaluated on an exact regular solution would lead to estimating the last integral via expressions of the type Z K ((P ex h P ex )I) : ' i rg ⇤ ↵ dx which provide at best an estimate of an O(h p+3 ), which is related to the approximation error, and it is independent on the quadrature. However the theoretical and numerical evidence indicates a stricter bound, as we will show shortly.

First we write the integral in the parametric space:

I 1 = Z K (P h (x(X))I)' i (X)rg ⇤ ↵ (x(X))J A (X)dX
The above integral can be written regrouping the terms as

I 1 = Z K f (X)J iA (X)dX
having separated for convenience the terms that are defined directly in the reference space J iA (X) = ' i (X)J A (X) from the terms defined in the physical space f (X) = f (x(X)). We remark now that we can write, within each parametric element, J iA as a 0 th order truncated Taylor series:

J iA = J iA (X ⇤ ) + r X J iA (X ⇤⇤ X ) • (X X ⇤ )
where, following the Lagrange representation of the Taylor series remainder, for a given X, the point X ⇤⇤ X lies along the segment (X ⇤ , X). We can easily estimate the remainder from

r X J iA (X ⇤⇤ ) • (X X ⇤ ) = ' i (X ⇤⇤ X )r X J A (X ⇤⇤ X ) • (X X ⇤ ) + r X ' i • J A (X ⇤⇤ X )(X X ⇤ )
with r X ' i constant and of O(1) in the P 1 case as already remarked. From (63) r X J A (X ⇤⇤ ) = O(h 3 ) and J A (X ⇤⇤ ) = O(h 2 ). It follows that the second term of the above expression is not of higher order so the rest of the truncated series is of O(h 3 ) .

We now also develop f in a Taylor series truncated at the degree k K which is the assumed polynomial exactness of the volume quadrature formula. Combining with the previous development, we obtain using standard multi-index notation and again a Lagrange representation of the Taylor series remainder:

I 1 = Z K k K X |↵|=0 1 ↵! D ↵ X f ⇤ (X X ⇤ ) ↵ J iA (X ⇤ )dX + Z K k K 1 X |↵|=0 1 ↵! D ↵ X f ⇤ (X X ⇤ ) ↵ r X J ⇤⇤ iA • (X X ⇤ )dX + Z K X |↵|=k K +1 1 ↵! D ↵ X f ⇤⇤ (X X ⇤ ) ↵ J iA (X ⇤ )dX + Z K X |↵|=k K 1 ↵! D ↵ X f ⇤⇤ (X X ⇤ ) ↵ r X J ⇤⇤ iA • (X X ⇤ )dX
where we recall that in multi-index notation D ↵ X = @ |↵| @X 1 ↵ 1 X 2 ↵ 2 . For brevity we have omitted the explicit dependence of the remainder terms on the X ⇤⇤ X point (which in principle is di↵erent for the two development considered), since it has no relevance in the following. By classical arguments, we claim now that the first two integrals, involving degree k K polynomials, represent the value obtained with the approximate formula used in our method, so that the remaining terms provide an estimate of the quadrature error:

E I1 := Z K X |↵|=k K +1 1 ↵! D ↵ X f ⇤⇤ (X X ⇤ ) ↵ J iA (X ⇤ )dX + Z K X |↵|=k K 1 ↵! D ↵ X f ⇤⇤ (X X ⇤ ) ↵ r X J ⇤⇤ iA •(X X ⇤ )dX
Since the volume of the reference element K is of an O(1), we need to estimate the derivatives of f , and combine with the available estimates J iA (X ⇤ ) = O(h 2 ), and r X J ⇤⇤ iA • (X X ⇤ ) = O(h 2 ). To estimate the derivatives of f we now use the fact that |r ↵ f | = |D ↵ f | and (60)-( 61) and obtain:

|r X f (x(X))| =|J A r x f |  C 1 h|D 1 x f | |r 2 X f (x(X))| =|r X J A r x f + J 2 A r 2 x f |  C 2 h 2 max m2 |D m x f | |r 3 X f (x(X))| =|r 2 X J A r x f + (r X J A )J A r 2 x f + 2J A r X J A r 2 x f + J 3 A r 3 x f |  C 2 h 3 max m3 |D m x f | etc
The kth derivative can be shown by induction (omitted for brevity) to be bounded by

|D k X f |  C k h k max mk |D m x f |
so that the quadrature error for the term under consideration can be estimated to

|E I1 |  C 1 I max kk K |D k x f |h k K +2 + C 1 I max kk K +1 |D k x f |h k K +3  C I1 max kk K |D k x f |h k K +2 (67)
We can proceed in a similar manner for the other integrals in [START_REF] Saha | The NCEP climate forecast system reanalysis[END_REF] to obtain:

I 2 = Z K (P h I) : r x ' i g ⇤ ↵ dx ) |E I2 |  C I2 max kk K |D k x f |h k K +1 I 3 = Z K r x P h • ' i g ⇤ ↵ dx ) |E I3 |  C I3 max kk K |D k x f |h k K +1
Note that for these terms, the presence of the derivative in the physical space r x let the estimate one order of magnitude larger in h with respect to (67). Similarly for the line integral

I 4 = Z @K P h Ig ⇤ ↵ ' i • n ds ) |E I4 |  C I4 max k@k K |D k x f |h k @K +1
By summing up the di↵erent integrals we have obtained for the Well Balanced correction (40) the following bound

|✏ K P | = |I 1 + E I1 + I 2 + E I2 + I 3 + E I3 + I 4 + E I4 | = |E I1 + E I2 + E I3 + E I4 |  Ch min(k K +1,k @K +1)
having used the fact that for exact quadrature the sum of all the integrals I j vanishes. This proves the result of the proposition. Note that the constant on the left actually depends on the integrand, so in particular on the pressure field and its derivatives, so it is bounded provided these are.

We perform a numerical verification with the values used for all the numerical simulations of the paper, namely k K = k @K = 3p. We consider a smooth pressure field P (', ) = 2 + sin(2⇡') cos(2⇡ ) defined on a spherical triangle which is successively refined halving the mesh size. In Figure [START_REF] Abgrall | High order methods for CFD[END_REF] we check the order of accuracy for |✏ K P | computed as |✏ K P | 1 = max K |✏ K P |. For P 1 , P 2 we recover numerically the estimate of Proposition (41), while for P 3 we measure one order less w.r.t such theoretical estimate.

  With h h + b h = ⌘ 0 =const, also the volume terms in 48 cancel exactly, recovering the well balanced property. The same arguments holds for the inverted barometer case but replacing b with B = b + patm g⇢0 ⇤.

  the error maps after 5 days of simulation on the fourth grid (h K = 446 km). Zonal Flow over an Isolated Mountain. W92 Case #5 is a perturbation of the previous test. The initial velocity and height fields are equal to (54) with ↵ = 0, h 0 = 5960 m and u 0 = 20 m/s. An isolated mountain is added: b = b 0 L2 [Bernard,2008] P 2 -L2 [Bernard,2008] P 3 -L2 [Bernard,2008]

Figure 1 :

 1 Figure 1: Global Steady-State Nonlinear Zonal Geostrophic Flow. Left: mesh convergence for L1 and L2 error norm. Right: mesh convergence for L2 error and comparison with [10].

Figure 2 :

 2 Figure 2: Global Steady-State Nonlinear Zonal Geostrophic Flow. Error plot. Top: P 1 . Bottom: P 3 . Note that color scale change of one order of magnitude for each plot.

Figure 3 :Figure 4 : 1 z lim R 2 (

 3412 Figure 3: Zonal flow over an isolated mountain. Left: Snapshot at day 15 (depth contour levels from 5050 m to 5950 mintervals of 50 m). Right: relative energy error (results from [5, 32] and from the DWD reported for comparison where relevant)

Figure 5 :

 5 Figure 5: Unstable Jet. Fine resolution h K = 111 km. Dark grey area represents a circle with radius R q 1

Figure 6 :

 6 Figure 6: Small amplitude tidal wave. Coarse h K = 100 km (left block) and fine h K = 50 km (right block) simulation. Top: mesh and bathymetry. Middle and bottom: snapshot of the free surface at t = 18h and 23h. Left: well balanced scheme. Right: non well balanced scheme. Contour levels are from 0.02 m to 0.02 m in intervals of 0.002 m.

Figure 7 :

 7 Figure 7: Small amplitude surge wave. Coarse h K = 100 km (left block) and fine h K = 50 km (right block) simulation. Top: mesh and bathymetry. Middle and bottom: snapshot of the free surface at t = 4h and 16h. Left: well balanced scheme. Right: non well balanced scheme. Contour levels are from 0.02 m to 0.02 m in intervals of 0.002 m.

Figure 8 :Figure 9 :

 89 Figure 8: Circular Hump. Contour levels at t = 0.9 from 0.04 m to 0.5 m in intervals of 0.045 m. Bottom: cut along the equator

Figure 10 :

 10 Figure 10: Realistic simulations in Caribbean Sea. Top: Full fine mesh. Bottom: Zoom of the fine mesh in correspondence of Martinique and Guadalupe island.

Figure 11 :

 11 Figure 11: Tsunami in the Carribean Sea. Color map of sea surface height ⌘ after 5400 s. Numbers (1) to (6) in the top picture show the position of the six stations considered in Figure 12, which corresponds to: (1) Pointe-a-Pitre, (2) Fort-de-France, (3) Martinique W, (4) Kingston, (5) Atlantic Ocean, (6) Martinique NE.

Figure 12 :Figure 13 :

 1213 Figure 12: Tsunami in the Carribean Sea. Water height recorded at the gauge stations located at Point-a-Pitre (topleft), Fort-de-France (top-right), Martinique W (middle-left), Kingston (middle-right), Atlantic Ocean (bottom-left) and Martinique NE (bottom-right)

Figure 14 :

 14 Figure 14: Hurricane Irma. Comparison of simulated sea surface height signals (scale in meters) at a selection of gauge stations, located at Saint Martin (top-left), Lamashur Bay (top-right), Charlotte Amalie (bottom-left) and Culebra (bottomright).

Figure 15 :

 15 Figure 15: Numerical validation of the well balanced correction estimate for a smooth pressure field. Left: Domain and pressure field. Right: mesh convergence for |✏ K P |1. The dashed and continuous lines represent the x and y components of the correction.

Table 1 :

 1 Global Steady-State Nonlinear Zonal Geostrophic Flow. L 1 ,L 2 error norm and convergence rate. The order is computed from a linear fit of the 5 points.

	grid # (↵ = 0)	e 1 P 1	e 2 P 1	e 1 P 2	e 2 P 2	e 1 P 3	e 2 P 3
	1	4.6261e-02 5.1549e-02 2.7469e-03 3.8020e-03 3.4459e-04 4.8972e-04
	2	8.5532e-03 9.8352e-03 2.2773e-04 2.9368e-04 1.4836e-05 2.1296e-05
	3	1.4574e-03 1.7889e-03 2.2788e-05 2.8277e-05 8.1059e-07 1.0731e-06
	4	2.8337e-04 3.6591e-04 2.5750e-06 3.1282e-06 4.7899e-08 6.3057e-08
	5	6.2480e-05 8.2598e-05 3.0803e-07 3.7606e-07 2.9741e-09 3.9206e-09
	Order	2.3980	2.3319	3.2711	3.3160	4.1919	4.2260

(polar cap region). Top: z lim = 0.9R. Bottom: z lim = 0.75R. Contour levels from 1.1e 4 to 1.5e + 4.

Acknowledgments

Work partially funded by the Carib-Coast INTERREG project (https://www.interreg-caraibes.fr/caribcoast). Some of the computations presented in this paper were carried out using the BRGM experimental testbed.

Appendix: Bounds on the integration by parts correction

To prove the bounds [START_REF] Pringle | Baroclinic coupling improves depth-integrated modeling of coastal sea level variations around Puerto Rico and the U.S. Virgin Islands[END_REF] we start by characterizing the mapping [START_REF] Ern | A well-balanced runge-kutta discontinuous galerkin method for the shallowwater equations with flooding and drying[END_REF] used for the quadrature with a classical property (see e.g. [START_REF] Ciarlet | Interpolation theory over curved elements, with application to finite element methods[END_REF][START_REF] Engvall | Mesh quality metrics for isogeometric Bernstein-Bezier discretizations[END_REF] for curved finite element examples).

Property 1. (Bounds on the mapping derivatives) The mapping [START_REF] Ern | A well-balanced runge-kutta discontinuous galerkin method for the shallowwater equations with flooding and drying[END_REF] and its derivatives verify the following estimates:

Proof. To prove the estimates we start by noting that, denoting by

with l j = (x l j x l G )/h the rescaled components of the distance of the nodal positions from the average, which is of O(1) for shape regular meshes, and with ' j the linear finite element basis. The above expression can be written as a linear multivariate polynomial in terms of the components of hX, with O(1) coe cients.

We write succinctly

and similarly

x l (X) = p 1 l (hX) kF k where now with p k (hX) we denote a multivariate polynomial of degree k in the components of hX, with O(1) coe cients. To obtain the first two estimates we prove that any k th order di↵erential of the mapping (w.r.t. any component or combination of the components of X) verifies

To prove by induction we check the initial condition: @x l (X) = kF k@p 1 l (hX) p 1 l (hX)@kF k kF k 2