Luca Arpaia

Mario Ricchiuto

Well balanced residual distribution for the ALE spherical shallow water equations on moving adaptive meshes

Keywords: Shallow Water equations, curvilinear coordinates, Arbitrary Lagrangian Eulerian formulation, Moving Mesh, Residual Distribution, Well-Balanced

We consider the numerical approximation of the Shallow Water Equations (SWEs) in spherical geometry for oceanographic applications. To provide enhanced resolution of moving fronts present in the flow we consider adaptive discrete approximations on moving triangulations of the sphere. To this end, we restate all Arbitrary Lagrangian Eulerian (ALE) transport formulas, as well as the volume transformation laws, for a 2D manifold. Using these results, we write the set of ALE-SWEs on the sphere. We then propose a Residual Distribution discrete approximation of the governing equations. Classical properties as the DGCL and the C-property (well balancedness) are reformulated in this more general context. An adaptive mesh movement strategy is proposed. The discrete framework obtained is thoroughly tested on standard benchmarks in large scale oceanography to prove their potential as well as the advantage brought by the adaptive mesh movement.

Introduction sec:intro

The Earth's curvature has a strong impact on wave propagation at large scales. For this reason, efficient and accurate techniques to solve flow equations on the sphere are of interest in many fields, such as ocean modeling, weather prediction as well as coastal hazards. Our objective is to explore the possibility of improving the discrete resolution of flows relevant to large-scale oceanographic applications by adaptively deforming an initial mesh, possibly already adapted to some of the data (e.g. initial depth or depth at rest). Indeed, mesh adaptation is useful to handle not only complex wave patterns, but also complex bathymetries and coastlines. In this respect, the possibility of working on triangular grids may be useful especially when approaching shallow areas. On these triangulations, in the spirit of the work done in ArR:17 [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF], we aim to assess the potential of adaptive mesh deformation, as well as its limitations for complex multi-scale computations.

To achieve this objective, there are several aspects to be dealt with. The first one is the discrete formulation of the Shallow Water Equations (SWEs) on the sphere. The literature on the topic is quite large. On structured quadrilateral grids, it is common to formulate the SWEs in covariant curvilinear coordinates (with the set of covariant coordinates identified by a patchwise parametrization of the sphere).

In the finite volume context, a thorough analysis is discussed in Ros:04,Ros:06 [START_REF] Rossmanith | A wave propagation algorithm for hyperbolic systems on curved manifold[END_REF][START_REF] Rossmanith | A wave propagation algorithm for hyperbolic systems the sphere[END_REF] and other formulations are proposed e.g. in UlJ:10,Cas:18 [START_REF] Ullrich | High-order finite-volume methods for the shallow water equations on the sphere[END_REF][START_REF] Castro | Well-balanced methods for the shallow water equations in spherical coordinates[END_REF]. On triangular meshes most works focus on high order local finite element-like mappings as in e.g.

Comb:09,Bern:09 [START_REF] Comblen | A finite element method for solving the shallow water equations on the sphere[END_REF][START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF], or on a fully 3D formulation with Lagrange multipliers to locally re-project vectors on the appropriate basis. A compromise is the initial work of Ros:13 [START_REF] Rossmanith | Residual distribution schemes for hyperbolic balance laws in generalized coordinates[END_REF], in which covariant coordinates are combined with triangular grids and with a Residual Distribution (RD) discretization. An ad-hoc patching at the poles guarantees the C 0 continuity of the solution. This provides us with an initial framework to extend our previous work. In particular, we have shown in RiB:09,r11,Ric:15 [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | On the C-property and generalized C-property of residual distribution for the shallow water equations[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] that RD method allows to handle very naturally source terms with a simple design of schemes preserving many known steady equilibria on unstructured triangulations. This is relevant when considering the SWEs in curvilinear coordinates which embed several geometrical and physical forcing terms. The RD method also allows to combine high order with monotonicity and positivity preservation principles. The work of Ros:13 [START_REF] Rossmanith | Residual distribution schemes for hyperbolic balance laws in generalized coordinates[END_REF] however is limited to steady state problems, on fixed grids, with constant bathymetry, and linear schemes. Here we present a more general development of nonlinear second order discrete approximations, for steady and time dependent problems on complex bathymetries.

Another important aspect is the design of the adaptation strategy on the sphere. The geophysics community has focused on quad-tree AMR techniques based on cell subdivision, see for example SkK:92 [START_REF] Skamarock | Adaptive grid refinement for two-dimensional and three-dimensional nonhydrostatic atmospheric flow[END_REF] or McU:15,Fer:16 [START_REF] Mccorquodale | An adaptive multi-block high-order finitevolume method for soliving the shallow water equations on the sphere[END_REF][START_REF] Ferguson | Analyzing the adaptive mesh refinement (AMR) characteristics of a high-order 2D cubed-sphere shallow-water model[END_REF], and

LGB:11,Pop:11 [START_REF] Leveque | Tsunami modelling with adaptively refined finite volume methods[END_REF][START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF] for tsunami propagation. Here we consider a different approach, consisting in generating an initial grid adapted to the problem's data (e.g. topography) and then using mesh deformation to track more accurately the flow features. An elegant way to deal with the underlying moving mesh is to formulate balance laws directly in an Arbitrary Lagrangian Eulerian (ALE) moving reference framework attached to the grid, see the original work of HIRT1974227 [START_REF] Hirt | An Arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF] and PS:03,KSD:12,Zho:13,ArR:17 [START_REF] Prusa | An all-scale anelastic model for geophysical flows: dynamic grid derformation[END_REF][START_REF] Kuhnlein | Modelling atmospheric flows with adaptive moving meshes[END_REF][START_REF] Zhou | An adaptive moving finite volume scheme for modeling flood inundation over dry and complex topography[END_REF][START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF] for application in geophysics.

Given an ALE map x = A(X, t) between a reference (fixed) domain and an actual (moving) one, there are two ways to formulate the PDEs in ALE form: either with respect to the reference coordinates X or the actual ones x. In the former case, the derivation of the ALE equations in the reference domain becomes immediate in a spacetime setting: the transformation of the 4-divergence differential operator under the ALE map introduces an additional ALE transport term. An overview of the fundamental relations necessary to achieve this in a spacetime setting are discussed in e.g. Ba08 [START_REF] Bazilevs | Isogeometric fluid-structure interaction: theory, algorithms, and computations[END_REF] (or chapter 2 in the book CFSI12 [START_REF] Bazilevs | Computational Fluid-Structure Interaction[END_REF]), see KSD:12 [START_REF] Kuhnlein | Modelling atmospheric flows with adaptive moving meshes[END_REF] for applications in geophysics, and D12,BOSCHERI201648,BOSCHERI2014484,VILAR2014188 [START_REF] Dobrev | High-order curvilinear finite element methods for Lagrangian hydrodynamics[END_REF][START_REF] Boscheri | High order accurate direct Arbitrary-Lagrangian-Eulerian ADER-WENO finite volume schemes on moving curvilinear unstructured meshes[END_REF][START_REF] Boscheri | A direct Arbitrary Lagrangian Eulerian ADER-WENO finite volume scheme on unstructured tetrahedral meshes for conservative and non-conservative hyperbolic systems in 3D[END_REF][START_REF] Vilar | A discontinuous Galerkin discretization for solving the twodimensional gas dynamics equations written under total Lagrangian formulation on general unstructured grids[END_REF] for finite volume/element formulations on curved meshes based on similar principles. There is certainly a similarity in the use of the reference coordinates X to represent the mesh movement and the aforementioned use of curvilinear coordinates to represent a geometrical manifold. Such analogy between the ALE form of balance laws and e.g. the governing equations of relativistic hydrodynamics, which are written in a time-dependent curved spacetime, is discussed in t16 [START_REF] Teukolsky | Formulation of discontinuous Galerkin methods for relativistic astrophysics[END_REF].

In moving mesh applications, another strategy used very often is to express the differential operators in the actual configuration, and thus discretize them on the actual (moving) mesh. Within such approach, a special role is played by the Geometric Conservation Law (GCL), which enforces conservation of the total volume of the moving domain ThL:79,LeF:96,GuF:00 [START_REF] Thomas | Geometric conservation law and its application to flow computations on moving grids[END_REF][START_REF] Lesoinne | Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations[END_REF][START_REF] Guillard | On the significance of the geometric conservation law for flow computations on moving meshes[END_REF], In ThL:79 [START_REF] Thomas | Geometric conservation law and its application to flow computations on moving grids[END_REF] it was proposed to replace the GCL with the constraint of reproducing exactly uniform flows, and in GuF:00 [START_REF] Guillard | On the significance of the geometric conservation law for flow computations on moving meshes[END_REF] it has been shown that satisfying the Geometric Conservation Law (GCL) is a sufficient condition for first order time accuracy on moving meshes. The compliance of the GCL at a discrete level is referred to as Discrete GCL (DGCL). Although rooted in fundamental relations of continuous mechanics Ba08,CFSI12 [START_REF] Bazilevs | Isogeometric fluid-structure interaction: theory, algorithms, and computations[END_REF][START_REF] Bazilevs | Computational Fluid-Structure Interaction[END_REF], the practical importance/benefit of the GCL and of its discrete counterpart is sometimes controversial, see e.g. Eti:09 [START_REF] Etienne | Geometric conservation law and finite element methods for 3D unsteady simulations of incompressible flow[END_REF].

The approach that we use is of mixed type SaA:08 [START_REF] Savidis | Arbitrary Lagrangian Eulerian finite element formulation for geotechnical construction processes[END_REF]. The spherical map in written with respect to the curvilinear covariant coordinates. This means that the parametrization of the sphere is embedded in the PDE system with the appropriate metric and source terms as in e.g. Ros:04,Ros:06,UlJ:10,Ros:13,Cas:18 [START_REF] Rossmanith | A wave propagation algorithm for hyperbolic systems on curved manifold[END_REF][START_REF] Rossmanith | A wave propagation algorithm for hyperbolic systems the sphere[END_REF][START_REF] Ullrich | High-order finite-volume methods for the shallow water equations on the sphere[END_REF][START_REF] Rossmanith | Residual distribution schemes for hyperbolic balance laws in generalized coordinates[END_REF][START_REF] Castro | Well-balanced methods for the shallow water equations in spherical coordinates[END_REF]. The ALE map, on the contrary is written in the actual configuration, which does not require the computation of the ALE Jacobian matrix. In this setting, the ALE-SWEs on a 2D manifold will be derived using classical transport formulas and volume conservation statement. A derivation of the governing equations in spacetime is also given in the Appendix. It will be shown that, for the SWEs on the sphere, the satisfaction of a DGCL is a necessary condition to retain well-balancedness. To this end, a new ALE closure for RD schemes is proposed, since the one of ArR:14 [START_REF] Arpaia | An ALE formulation for explicit Runge-Kutta residual distribution[END_REF] does not extend to the sphere, and more in general to manifolds.

The schemes obtained are thoroughly tested both on fixed and moving adaptive meshes on classical applications for the shallow water on the sphere, and on a realistic computation of relevance in oceanography: the regional propagation and coastal impact of the 2011 Tohoku-Honsu tsunami.

The paper is organized as follows. The standard SWEs in spherical geometry are recalled in section sec:SW-sphere sec:SW-sphere In this section we recall the Shallow Water Equations (SWEs) on the sphere in covariant curvilinear coordinates. At every point of the sphere S 2 we fix an orthogonal but not orthonormal reference system with covariant basis vectors {g 1 , g 2 } together with local curvilinear coordinates {x 1 , x 2 }. The vectors g 1 and g 2 define the tangent plane to the sphere and are used to construct all quantities allowing to locally describe the manifold. For completeness, in app:a1.1 app:a1.1

Sections

Appendix A.1 we recall some basic concepts from differential geometry, as well as many definitions used throughout the following sections. We will consider the evolution of the fluid depth h(x, t) measured along the radial axis. The bathymetry/topography height with respect to a given reference level is denoted by b(x). The fluid free surface is defined here as η = h + b. We denote hu i the components of the discharge vector hu = hu i g i . We remark that the velocity components are not defined in a unitary basis. The SWEs in covariant curvilinear coordinates

write 1 : ∂h ∂t + 1 √ G ∂ ∂x j √ G hu j = 0 ∂hu i ∂t + 1 √ G ∂ ∂x j √ G T ij + S i = 0
(1) eq:SW-s G denotes the determinant of the metric tensor defined in ( eq:Gij_def eq:Gij_def

A.3) in

app:a1 app:a1

Appendix A, and

T ij = huu ij + 1 2 G ij gh 2
are the components of the momentum flux T = T ij g i g j , involving the components of the inverse metric tensor, also defined in the appendix.

The source term in the momentum balance reads

S i = G ij gh ∂b ∂x j + c F hu i + S i γ -S i c ( 2
)
1 standard summation convention is used, for the notation we refer again to app:a1 app:a1 Appendix A

Besides the effects of bathymetry and friction (first two terms), the source includes also a geometrical force due to Earth curvature S i γ , and Coriolis term S i c . The geometric source arises from the space variation of the covariant basis vectors and reads

S i γ = Γ i jk T jk (3)
with Γ i jk being the Christoffel symbols. The Coriolis term is associated to the rotation of the Earth, which plays an important role in ocean dynamics. It reads

S i c = f G ij ε jk hu k (4)
with the parameter f = 2Ω sin λ with λ the latitude and Ω the Earth rotation rate, always taken as Ω = 7.292 × 10 -5 s -1 unless differently specified, and having denoted by ε the 2D Levi-Civita symbol.

The SWEs system can also be expressed in a compact vector form

∂u ∂t + 1 √ G ∂F j ∂x j + S(x, u) = 0, (5)
u =   h hu i   , F j = √ G   hu j T ij   , (6)
S =   0 G ij gh ∂b ∂x j   S b +   0 c F hu i   S f +   0 Γ i jk T jk   S γ -   0 f G ij ε jk hu k   S c (7)
Equations ( eq:SW-sph eq:SW-sph

5)-(

eq:SW2-sph eq:SW2-sph 7) constitute a non-homogeneous hyperbolic system of PDEs. As in Lev:97 [START_REF] Leveque | Wave propagation algorithm for multidimensional hyperbolic systems[END_REF] we remark that the flux vector depends explicitly on the position making the system non-autonomous.

Arbitrary Lagrangian Eulerian framework on the sphere sec:ale

We now introduce an independent configuration A, the Arbitrary Lagrangian Eulerian (ALE) configuration, which we assume to belong to a differentiable manifold (the sphere). For us A will be the configuration the mesh is attached to. A is composed of points Q identified by the vector x which can be expressed in the local curvilinear coordinate system x(Q) = x i g i , with i = 1, 2. The initial or reference configuration of the continuous mesh is denoted by A 0 and the reference position of the points is denoted by greek letters χ( Q). For this configuration we assign a local covariant vector basis {γ 1 , γ 2 } and local chart {χ 1 , χ 2 } of the tangent plane so that χ( Q) = χ i γ i . As the mesh moves, also A does. We assume to be able to uniquely map the initial configuration to the reference one via the mapping:

A : A 0 → A x = A(χ, t) (8)
The function A is assumed to be continuously differentiable, with a smooth Jacobian J A = ∂x ∂X . Classical requirements on the Jacobian J A are invertibility ∃J -1 A , and the positivity of its determinant J A = det J A > 0. For a spacetime definition of the ALE transformation see (

app:a1.2 app:a1.2 Appendix A.2).
The fluid is endowed with various physical properties which may be represented by scalar, vector and tensor fields. These fields and their rate of change in time and space may be defined in the reference configuration or in the actual one. Here we choose the actual description: that is, in the following, all spatial derivatives are expressed with respect to the actual coordinates. The ALE time derivative of a scalar φ is its time rate of change measured by an observer which moves with the mesh point labeled with χ, where χ( Q) is the particle's position in the reference configuration. We will shorten it ∂φ ∂t χ

. The ALE/mesh velocity of a point Q in the actual vector basis, σ(Q) = σ j g j is calculated by applying the aforementioned definition of ALE derivative, to the position x(Q). Its components write

σ j = ∂x j ∂t χ (9)
Using chain rule, we can define relations between Eulerian and ALE time derivatives for these properties:

∂φ ∂t χ = ∂φ ∂t + ∂φ ∂x j σ j (10)
Finally, taking the ALE time derivative of an infinitesimal area on the sphere we obtain the Geometric Conservation Law on the manifold Tru:66 [START_REF] Truesdell | The mechanical foundations of elasticity and fluid dynamics[END_REF]: [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] eq:gcl-s or, more compactly

∂ ∂t χ √ GJ A = J A ∂ ∂x j √ Gσ j (
∂ t | χ ( √ GJ A ) -J A √ G∇ • σ = 0.
The GCL ( eq:gcl-s eq:gcl-s 11) is a geometrical relation between the ALE Jacobian and the ALE velocity. One can check that the integral of ( eq:gcl-s eq:gcl-s 11) over the whole computational domain states the conservation of the total area. At a discrete level this corresponds to preserve the total area of the mesh during the simulation.

For this reason in section sec:dgcl-s sec:dgcl-s 5.3 we will discuss how to preserve also its discrete counterpart.

ALE-SWEs on the sphere ale-SW-sphere

Using the relations derived in the previous sections, and more particularly ( eq:tra-a_C eq:tra-a_C 10) and ( eq:gcl-s eq:gcl-s 11), we can easily recast the SWEs ( eq:SW-s eq:SW-s 1) in an ALE form, obtaining the ALE-SWEs on a 2D manifold with respect to the actual parametric coordinates:

∂ √ GJ A h ∂t χ + J A ∂ ∂x j √ Ghu j - √ Ghσ j = 0 (12) ∂ √ GJ A hu i ∂t χ + J A ∂ ∂x j √ GT ij - √ Ghu i σ j + √ GJ A S i = 0 ( 13
)
As before, the system can be also cast in compact vector form as

∂ ∂t χ √ GJ A u + J A ∂ ∂x j F j - √ Gσ j u + √ GJ A S(x, u) = 0 (14)
Finally, it is useful for the following to recall that the bathymetry, which does not depends explicitly on time b = b(x), verifies the following identity:

∂b ∂t χ = ∂b ∂x j σ j (15)
multiplying by J A √ G and summing it to the GCL ( eq:gcl-s eq:gcl-s 11) (multiplied by b) we obtain the ALE remap on the 2D manifold

∂ √ GJ A b ∂t χ -J A ∂ ∂x j √ Gbσ j = 0 ( 16
) eq:b_ALE_remap
Please note that scalar, vector, tensor fields and spatial derivatives are written w.r.t. the actual coordinates

x, in all the expressions above. For completeness in ( app:a2 app:a2

Appendix B) we report an alternative derivation of the governing equations ( eq:SW1-as eq:SW1-as 14) with a spacetime approach.

Discretization: parametrization and constraints

eq:not-as

In this section we discuss the basic aspects of the discretization. These involve the main geometrical notations as well as the form of the equation used for the discretization. We will introduce an appropriate form of the ALE-SWEs allowing to combine well-balancedness with geometric conservation. We then show how to embed discrete analogs of these constraints in a residual distribution setting.

Unstructured grid representations: main notation

Points on a manifold are represented by a proper parametrization through the use of geometrical mappings. In particular, when working in ALE formalism, a point i will be expressed in the parametrization of the initial configuration by a vector χ i = {χ 1 , χ 2 } i , while x i = {x 1 , x 2 } i will refer to the configuration attached to the mesh. Note that the discussion here is general, even though the numerical applications will only consider a latitude longitude parametrization of the sphere. This particular choice is known to introduce a singularity at the poles, which is nowadays well handled by meteorologists Wil:07 [START_REF] Williamson | The evolution of dynamical cores for global atmospheric models[END_REF]. Since we are interested in oceanic applications we circumvent this issue by adopting a simple fix. A patch-wise parametrization of the sphere (see e.g.

KaS:04 [START_REF] Kageyama | Yin-Yang grid": An overset grid in spherical geometry[END_REF] RoI:96 [START_REF] Ronchi | The "cubed sphere": A new method for the solution of partial differential equations in spherical geometry[END_REF] and references therein) is deployed, namely a polar cap defined by a limiting latitude λ lim , in which we use a rotated reference system (x , y , z ) = (x, z, -y). To account for this "patching", a rotation operator needs to be applied when updating vectors in elements crossing the patch boundary. In our case, this applies both to physical velocities, as well as mesh nodes velocities and position. In all right hand side evaluations this is accounted for following Ros:13 [START_REF] Rossmanith | Residual distribution schemes for hyperbolic balance laws in generalized coordinates[END_REF] to which we refer for details.

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 each element, we define the spherical area and median dual cell of a node as

|K(t)| = K(t) √ Gdx, |C i | = K∈Di |K| 3 ( 17
) eq:sp-tr
For moving meshes we also have:

|K(t)| = K0 J A (t) √ G dχ . ( 18
) eq:sp-tr-ALE
denoting by K 0 the element in the initial parametric space. Please note that, after the change of coordinates, we assume in the above expression that G(x(χ, t)). This duality of integral formulas via the ALE mapping applies to all integrals and will not be recalled explicitly unless necessary.

For a given mesh, we label by h the largest element diameter in the parameter space M 2 . For a smooth enough initial mesh, and on a smooth manifold, we unduly assume that the number of elements [START_REF] Kuhnlein | Modelling atmospheric flows with adaptive moving meshes[END_REF] eq:mesh-smooth for some bounded constant. We also denote by u h the P 1 approximation of a quantity u [START_REF] Zhou | An adaptive moving finite volume scheme for modeling flood inundation over dry and complex topography[END_REF] eq:fem having denoted by ϕ i the standard P 1 bases. We now define elemental contributions to the Galerkin mass matrix (m GAL ij [START_REF] Bazilevs | Isogeometric fluid-structure interaction: theory, algorithms, and computations[END_REF] eq:gal-mass as well as to the lumped mass matrix ) eq:gal-mass-L 5.2. Well-Balancedness and GCL on the sphere sec:wb-dgcl

N K verifies N K |K| ≤ c|M 2 | (
u h (x) = i u i ϕ i (x) (
) K := K ϕ i ϕ j √ G(x) dx , (
(m GAL ij ) L,K := K i δ ij , K i := K √ G(x) ϕ i dx . ( 22
For shallow water flows, a very important role is played by the lake-at-rest steady state: [START_REF] Dobrev | High-order curvilinear finite element methods for Lagrangian hydrodynamics[END_REF] eq:lake_at_res

hu = 0, h + b = η 0 = const (
The exact preservation of this (initial) steady state at the discrete level is referred to as the C-property,

or Well-Balancedness

BeV:94 [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF]. In other words, well-balanced schemes provide a discrete analog of the balance (cf. equation ( eq:SW-sph eq:SW-sph 5))

∂F j ∂x j + √ GS = 0. ( 24
)
Written in the form ( eq:SW1-as eq:SW1-as [START_REF] Ferguson | Analyzing the adaptive mesh refinement (AMR) characteristics of a high-order 2D cubed-sphere shallow-water model[END_REF], the ALE-SWEs require the simultaneous satisfaction of several compatibility constraints in order to devise well-balanced schemes. These constraints include:

1. well-balancedness for the Cartesian SWEs on fixed grids: Appendix A);

∂ x i (gh 2 /2) + gh∂ x i b = 0, ∀ i ; 2.
We will show in the following of the paper how we can embed exactly the first two constraints at a discrete level. Concerning the third, replacing u = 0 in the first of ( eq:SW1-as eq:SW1-as [START_REF] Ferguson | Analyzing the adaptive mesh refinement (AMR) characteristics of a high-order 2D cubed-sphere shallow-water model[END_REF] we get

0 = ∂ t ( √ G J A h) χ -J A ∂ x j ( √ G σ j h) 0 = ∂ t ( √ G J A η) χ -J A ∂ x j ( √ G σ j η) -∂ t ( √ G J A b) χ -J A ∂ x j ( √ G σ j b) (25)
If η = η 0 we see that a necessary condition for the lake-at-rest to be a solution for system ( eq:SW1-as eq:SW1-as [START_REF] Ferguson | Analyzing the adaptive mesh refinement (AMR) characteristics of a high-order 2D cubed-sphere shallow-water model[END_REF] is to verify both the GCL, and the ALE remap of the bathymetry. Some authors ZhC:13 [START_REF] Zhou | A well-balanced stable GRP scheme for shallow water equations for adaptive unstructured triangular meshes[END_REF] have proposed to directly use a discretization of ( eq:b_ALE_remap eq:b_ALE_remap [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF]) to evolve the topography on a moving mesh. The drawback of this approach is that it significantly modifies the bathymetry during the computation, which cannot be accepted in practical applications. For this reason, here we follow ArR:17 [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF] and modify ( eq:SW1-as eq:SW1-as 14) by solving instead the system obtained by adding to the first equation the ALE remap of the bathymetry: [START_REF] Vilar | A discontinuous Galerkin discretization for solving the twodimensional gas dynamics equations written under total Lagrangian formulation on general unstructured grids[END_REF] eq:SW-aswb

∂ ∂t χ √ GJ A w + J A ∂ ∂x j F j - √ Gσ j w + √ GJ A S(x, u) = 0 , w =   η hu   (
Note that, with this notation, we have F j (u) = F j (w; b), as the bathymetry is required to evaluate the flux starting from the w array. Concerning now the fourth requirement, we consider the momentum balance obtained when replacing the lake-at-rest solution in ( eq:SW-aswb eq:SW-aswb 26). Setting Π = gh 2 /2 for short, we have [START_REF] Teukolsky | Formulation of discontinuous Galerkin methods for relativistic astrophysics[END_REF] eq:wbs

0 = J A ∂ x j ( √ GG ij Π) + √ GΓ i jk G jk Π + gh √ GG ij ∂ x j b 0 = J A √ GG ij (gh∂ x j b + ∂ x j Π) + Π ∂ x j ( √ GG ij ) + √ GΓ i jk G jk (
The first group in parentheses vanishes if the first property (well-balancedness for Cartesian SWEs on fixed grids) is satisfied. The remainder is nothing else that the combination of the two metric relations.

Indeed just by replacing the second relation of ( eq:Ric eq:Ric

A.6) into Ricci's Lemma we obtain 0 = √ G ∂G ij ∂x j + G ij Γ m mj + Γ i jm G mj 0 = √ G ∂G ij ∂x j + G ij ∂ √ G ∂x j + √ GΓ i jm G mj = ∂ √ GG ij ∂x j + √ GΓ i jm G mj (28)
The symmetry of the Christoffel symbols Ngu:14 [START_REF] Nguyen-Shafer | Tensor Analysis and Elementary Differential Geometry for Physicist and Engineers[END_REF] allows to recover exactly the last expression in parentheses in ( eq:wbs eq:wbs 27). The above relation is thus necessary for the lake-at-rest solution to be verified. At the discrete level, the easiest way to make sure this condition is embedded is to combine high order quadrature with the direct use of the analytical formulas. Even though low quadrature errors can be obtained with relatively few quadrature points, due to the smoothness of the functions involved, in this way the C-property would be satisfied within this quadrature error. Also, compared to the Cartesian case, the necessity of using higher order quadrature would still give a considerable overhead. For this reason we have chosen to embed the metric relations analytically. Whenever possible we discretize the momentum balance in the following form:

∂ ∂t χ √ GJ A hu i + J A ∂ ∂x j √ Ghu i (u j -σ j ) + J A √ G G ij ∂ ∂x j (gh 2 /2) + gh ∂ ∂x j b +J A √ GS i = 0 . ( 29
) eq:qdm-wb
For later use in the paper, we introduce the short notation

∂(J A √ Gw) ∂t χ + J A R(w, b) = 0
(30) eq:ut_lu

for the system obtained with the first in ( eq:SW-aswb eq:SW-aswb

26), plus (

eq:qdm-wb eq:qdm-wb 29).

DGCL on the sphere sec:dgcl-s

We consider now the satisfaction of the discrete analog of the GCL an a 2D manifold ( eq:gcl-s eq:gcl-s 11). As already said, differently from the Cartesian case, a uniform flow is not an exact solution of ( eq:SW1-as eq:SW1-as

14) embedded with the GCL (

eq:gcl-s eq:gcl-s 11), even for constant bathymetries and in absence of friction and Coriolis force. As a consequence we cannot define the DGCL through the original characterization of ThL:79 [START_REF] Thomas | Geometric conservation law and its application to flow computations on moving grids[END_REF], i.e. the preservation of a uniform flow. On the sphere, we could replace the preservation of uniform flow with some other particular cases, as e.g. the preservation of zonal or meridional flows. However, the most relevant characterization of the GCL for the applications of this paper is the one we have introduced in the previous section, namely the satisfaction of the lake-at-rest state.

It is thus in the perspective of obtaining well balanced schemes on moving grids that we consider the issue of mimicking ( eq:gcl-s eq:gcl-s 11) at the discrete level. As we will shortly see, the main issue is properly defining the components of the mesh velocity σ j along certain directions, as well as an approximation of the divergence operator ∂ x j ( √ Gσ j ). The DGCL closure used in previous works by the authors ArR:14,ArR:17 [START_REF] Arpaia | An ALE formulation for explicit Runge-Kutta residual distribution[END_REF][START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF], and inspired by LeF:96,GuF:00 [START_REF] Lesoinne | Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations[END_REF][START_REF] Guillard | On the significance of the geometric conservation law for flow computations on moving meshes[END_REF], is based on a piecewise linear P 1 interpolation of the mesh displacements. This definition allows to recover naturally an element-by-element DGCL when integrating the divergence of the displacements over the half time step average element

K n+1/2 = (K n+1 + K n )/2. Due to the presence of the coefficient √ G,
the extension of this approach is not at all straightforward. We have thus chosen to use a generalization of the closure originally proposed in MaY:06 [START_REF] Mavriplis | Construction of the discrete geometric conservation law for high-order time-accurate simulations on dynamic meshes[END_REF]. To describe this closure, we consider an element K of the mesh, or equivalently its projection in the parametric space K. We study the discrete approximation of the area variation

|K n+1 | -|K n | = K0 J A √ Gdχ n+1 - K0 J A √ Gdχ n
We decompose this variation in components associated to the area swept by each face, as illustrated in 1. We have in particular

|K n+1 | -|K n | = ∆t j∈K υ fj . ( 31
) eq:delta-K Denoting by Q j the quadrilateral ({x i , x k } n ; {x i , x k } n+1 ) (cf. figure fig:swept-area fig:swept-area 1
), υ fj is the average edge velocity tangent to the manifold and normal to the edge defined by [START_REF] Savidis | Arbitrary Lagrangian Eulerian finite element formulation for geotechnical construction processes[END_REF] eq:ivc-s

υ fj := 1 ∆t Qj √ Gdx (
In practice, the integral is computed by a high order quadrature formula

Qj √ Gdx = |Q j | xq ω q √ G(x q ) ( 33
) eq:ivc-s1
where the quadrature grid is built simply by tensorizing 1D quadrature points over each edge, and the signed volume |Q j | is evaluated using the surveyor's (or shoelace) area formula shoelace [START_REF] Braden | The surveyor's area formula[END_REF].

To discretize the ALE-SWEs we now need to provide a discrete approximation of the curvilinear divergence appearing in the GCL. For the purposes of the paper we would like this definition to be compatible with the following relations

K n+1 -K n = t n+1 t n K0 J A ∂ ∂x j ( √ Gσ j ) K dχ dt (34)
K n+1 i -K n i = t n+1 t n K0 J A ϕ i ∂ ∂x j ( √ Gσ j ) K dχ dt ( 35
) K0 J n+1 A √ Gϕ i w h - K0 J n A √ Gϕ i w h = t n+1 t n K0 J A ϕ i ∂ ∂x j ( √ Gσ j ) K w h dχ dt (36)
with w h the approximation of the same data on the different meshes associated to the mappings with Jacobians J n+1 A and J n A . The first relation is the spacetime integral of the GCL while the second and the third relations are useful to derive a consistent finite element method, and mass lumping. The practical evaluation of each of the integrals on the left hand side plays here a fundamental role in our ability to satisfy exactly the three relations above. We have the following characterization.

Proposition 1. (Approximate quadrature on the manifold and DGCL) Consider a one point quadrature approximating the integral over a time step as

t n+1 t n K0 J A f (χ) dχ dt = ∆t K0 JA f (χ) dχ dt = ∆t K f (x) dx dt
and the constant-by-element approximation of the divergence

∂ ∂x j ( √ Gσ j ) K := 1 | K| j∈ K υ fj .
(37) eq:div-dgcl Consider also the approximate quadrature defining full and lumped mass matrices ( eq:gal-mass-L eq:gal-mass-L

22) and (

eq:gal-mass eq:gal-mass [START_REF] Bazilevs | Isogeometric fluid-structure interaction: theory, algorithms, and computations[END_REF] given by

(m GAL ij ) K ≈ |K| 3 1 + δ ij 4 , (m GAL ij ) L,K = K i δ ij ≈ |K| 3 δ ij . ( 38
) eq:gal-mass-d
We have the following properties:

1. Relations ( eq:dgcl-rels1 eq:dgcl-rels1

34), (

eq:dgcl-rels2 eq:dgcl-rels2

35) and (

eq:dgcl-rels3 eq:dgcl-rels3 Appendix C).

Remark 1. (Higher order quadrature in mass matrices). The approximate quadrature ( eq:gal-mass-d eq:gal-mass-d 38) is only required to embed exactly the GCL in the scheme. In absence of mesh movement, higher order (or exact) quadrature can be used to define the mass matrix entries, which will differ from ( eq:gal-mass-d eq:gal-mass-d 38) within the quadrature error just proved. In practice we have seen no advantage in using more accurate quadrature in these formulas. For completeness, the modifications necessary to account (on fixed meshes) for the exact evaluation of the mass matrices for Residual Distribution and stabilized finite elements are provided in ( app:a4 app:a4

Appendix D).

Note that in practice, when using the sphere with lat-lon parametrization, the gradient of √ G appearing in the error estimate in the appendix is given by R 2 sin λ with λ the latitude, which is clearly bounded.

Finally, note that the choice of JA is arbitrary and has no impact on the satisfaction of the DGCL. For accuracy reasons, here we have kept a mid-point approximation JA = (J n A + J n+1 A )/2. This means that, most geometrical quantities are evaluated as averages on the elements K n+1 and K n . This is implicitly assumed in all the formulas which will follow. For simplicity, in the text we will drop the superscript altogether, unless necessary in the formula.

ALE Residual Distribution on the sphere sec:RD

A first extension of Residual Distribution (RD) to general covariant coordinates has been proposed by Ros:13 [START_REF] Rossmanith | Residual distribution schemes for hyperbolic balance laws in generalized coordinates[END_REF]. The work of the reference is limited to steady state computations, with linear schemes and without bathymetry. This paper presents a general second order and well balanced formulation allowing time dependent computations on moving meshes. The construction proposed here is a generalization of the work done in Cartesian coordinates in ArR:17 [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF]. We refer the interested reader to Dec:17,AR:17 [START_REF] Deconinck | Residual distribution schemes: Foundations and analysis[END_REF][START_REF] Abgrall | High order methods for CFD[END_REF] and references therein for a thorough discussion of the properties of residual distribution.

The RD discretization of unsteady problems can be described starting with an analogy with stabilized continuous finite elements. Given a triangulation of the parametric domain S 2 , to solve L(w h ) = 0, one can formally recast a RD scheme as

S 2
w i L(w h ) dχ = 0 (39) eq:pg where w i = ϕ i + γ i is a Petrov-Galerkin test function, with γ i an appropriately defined bubble function.

For consistency, on each element

K j∈K w j = 1 , (40) eq
:rd-consiist which implies j γ j = 0. We define the element residual associated to ( eq:SW-aswb eq:SW-aswb 26) as the area integral over the element of the full differential operator:

Φ K = K0 L(w h ) dχ . (41)
The consistency ( eq:rd-consiistency0 eq:rd-consiistency0 40) allows to define the distributed residuals

Φ K i := K0 w i L(w h ) dχ (42) eq:rd-split0
which verify the consistency condition

j∈K Φ K j = Φ K . ( 43
) eq:rd-consiist
In actual computations, the main component of the distributed residuals is evaluated by means of distribution matrices β K i which also verify consistency

j∈K β K j = I . ( 44
) eq:rd-consiist
The bridge with finite elements is completed by the requirement

β K i = 1 |K 0 | K0 w i dχ = wK i . (45)
which leads to form of the distribution used in practice: [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF] eq:rd-split1

Φ K i = β K i Φ K + K0 (w i -wK i )L(w h ) dχ (
The formula above shows the key design elements of the method which are the integration strategy used to evaluate Φ K , and the definition of the distribution matrices and of the associated test functions. Note that the first term in the above formula is enough to have a consistent and conservative method (in the sense of ( eq:rd-consiistency1 eq:rd-consiistency1

44)

). The second term has an impact on the accuracy and stability, but none on consistency as clearly j (w j -wK j ) = 0.

Two-step explicit method

Concerning the time marching technique, we consider here a fully explicit multi-stage strategy, relying on the ideas discussed in RiA:10,Ric:15 [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF]. For recent implicit space-time strategies we refer to hr11,shr13,hrs20 [START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF][START_REF] Sármány | Unconditionally stable space-time discontinuous residual distribution for shallow-water flows[END_REF][START_REF] Hubbard | Space-time residual distribution on moving meshes[END_REF].

Concerning the scheme used here, we recall for completeness in ( app:a4 app:a4

Appendix D) its construction which, starting from ( eq:pg eq:pg [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF]), relies on the following steps:

1. replacing L(w h ) by ( eq:ut_lu eq:ut_lu 30);

2. mapping integrals on the reference domain and their evaluation on the current configuration (mesh);

3. use of a simplified residual based stabilization involving an extrapolation of the time derivative;

4. use of Proposition 1, as well as of ( eq:dgcl-rels1 eq:dgcl-rels1

34), (

eq:dgcl-rels2 eq:dgcl-rels2

35), and (

eq:dgcl-rels3 eq:dgcl-rels3 36);

5. lumping of the Galerkin mass-matrix.

In its final form, the scheme evolves nodal values of the state vector w with a two-step procedure, in which the new solution value is obtained from:

|C n+1 i |(w n+1 i -w * i ) = -∆t K i Φ K i (47) RC-corr
where the nodal residuals Φ K i are a splitting of the element residual

j∈K Φ K j = Φ K = 1 ∆t K0 (J n+1 A √ Gw * h -J n A √ Gw n h )dχ + 1 2 Φ K (w * h , b n+1 h ) + 1 2 Φ K (w n h , b n h ) (48) phi-corr
with Φ K (w, b) the fluctuation (the area integral over the element K of only advective and source terms)

Φ K = ∂K √ G   (hu j -ησ j ) hu i (u j -σ j )   h n j ds+ K √ G   0 G ij (∂ x j (gh 2 /2) + gh∂ x j b)   h dx + K √ G   0 S f + S γ + S c   h dx (49) phi
Please note that the fluctuation integral is evaluated with respect to the actual parametric coordinates after using dx = JA dχ. The predicted nodal values w * i are computed from

|C n+1 i |(w * i -w n i ) = -∆t K i ΦK i (w n h , b n h ) (50) RC-pred
with the ΦK i defining a splitting of the geometrically non-conservative fluctuation (cf. (

app:a4 app:a4 Appendix D)) j∈K ΦK j = ΦK = ∂K √ G   hu j hu i u j   h n j ds + K √ G   0 G ij (∂ x j (gh 2 /2) + gh∂ x j b)   h dx - K √ Gσ j ∂w h ∂x j dx + K √ G   0 S f + S γ + S c   h dx (51) phi-t

Evaluation of the integrals

To begin with, unless explicitly specified, in all quadrature formulas we have used analytical expressions for both the Jacobian √ G, and the metric tensor G ij . At the quadrature points, the mass flux components (hu i ) h are evaluated simply using the P 1 expansion, while velocity components are obtained from the interpolated values of momentum (hu i ) h and depth h h . The wet/dry region is handled as discussed in

RiB:09,Ric:15 [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF]. Consistently with the approximation ( eq:gal-mass-d eq:gal-mass-d [START_REF] Ronchi | The "cubed sphere": A new method for the solution of partial differential equations in spherical geometry[END_REF] of Proposition 1, volume integrals of the unknown are computed as

K0 J n+1 A √ Gw h dχ = j∈K |K n+1 j |w j = |K n+1 | 3 j∈K w j ( 52
) eq:vol-int
In all boundary integrals, scaled constant-by-face values of the normal mesh velocity σ j n j are computed consistently with the DGCL closure ( eq:ivc-s eq:ivc-s

32) and (

eq:div-dgcl eq:div-dgcl 37), namely

( √ Gσ j n j ) f := υ f ( 53
) eq:sigma-n
Concerning the non-conservative ALE term in ΦK , this has been evaluated as

K √ Gσ j ∂w h ∂x j dx = ∂K √ Gw h σ j n j ds - K w h ∂ √ Gσ j ∂x j dx (54)
where the face values ( eq:sigma-n eq:sigma-n 53) are used in the boundary integral, while the element integral is evaluated by means of the constant divergence of the DGCL ( eq:div-dgcl eq:div-dgcl 37):

K √ Gσ j ∂w h ∂x j dx = i∈K υ fi w fi -w K |K n+1 | -|K n | ∆t = i∈K υ fi (w fi -w K ) ( 55
)
where the average values w f and w K are simple arithmetic averages over the faces and over the element.

The integrals of the source terms in ( phi phi

49) and (

phi-t phi-t 51) are computed as

K √ GS dx ≈ |K| 3 j∈K S j (56)
Finally, the hydrostatic term is computed as

K √ GG ij (∂ x j (gh 2 h /2) + gh h ∂ x j b h )dx = K √ GG ij gh h ( ∂ x j η h ∂ x j h h + ∂ x j b h )dx = K √ GG ij gh h dx ∂ x j η h K ≈ |K| q ω q √ GG ij (x q )gh h (x q ) ∂ x j η h K (
57) eq:hyd Several (second or higher order) quadrature formulas have been tested for the term in parentheses without any visible impact on the results.

Distribution of the residual and main properties

In RD schemes stability and accuracy properties depends on how the residual is distributed. The computation of the split residuals for the two RK steps Φ K i and ΦK i is performed following exactly the steps discussed in Ric:15 [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF]. In particular, following the construction in the reference, the residual splitting ( eq:rd-split1 eq:rd-split1

46) can be more explicitly written as [START_REF] Weller | Mesh adaptation on the sphere using optimal transport and the numerical solution of a Monge-Ampere type equation[END_REF] eq:RD

Φ K i = β K i Φ K + δ(w h )β SD i Φ K + δ(w h ) j∈K m * ij w * j -w n j ∆t (
In the above expression, the first term is the consistent one, namely the one satisfying j β K j = 1. The last two terms do not contribute to the consistency of the method, as indeed one can show that Ric:15,RiA:10 [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF][START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF] j∈K

β SD j∈K Φ K + j∈K k∈K m * jk w * k -w n k ∆t = 0 (59)
These additional terms are essentially stabilization terms allowing to avoid spurious modes appearing in smooth regions (cf.

Abg:06,RiB:09 [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] for details), and can be embedded in the last integral in ( eq:rd-split1 eq:rd-split1

46

). The first one is a streamline upwind type stabilization, the second is a correction to a mass-lumping operator. These terms are in practice evaluated exactly as discussed in Abg:06,RiB:09,Ric:15 [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] to which we refer for details. Concerning the first term, as in Ric:15 [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] it is a distribution based on a blending, via the smoothness sensor δ(w h ), between a central splitting and a limited nonlinear splitting. The latter is obtained using a standard procedure in RD schemes involving the application of a bounded and sign-preserving mapping to a first order and positivity preserving Lax-Friedrich's like distribution (see for example Abg:06 [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] or the chapters AR:17,Dec:17 [START_REF] Abgrall | High order methods for CFD[END_REF][START_REF] Deconinck | Residual distribution schemes: Foundations and analysis[END_REF]). This allows to combine formal second order of accuracy with depth non-negativity preservation.

Finally note that the distribution functions for the predictor step ( RC-pred RC-pred 50) are easily deduced from ( eq:RD eq:RD

58) setting

w * = w n , and replacing Φ K by ΦK .

Remark 2. (Metric terms treatment, accuracy, stability). It is worth noting that the scheme is formally identical to the one obtained in the planar case. The treatment of the geometric source term boils down to the modification of the cell residual. This is a substantial simplification when conceiving and coding the scheme compared to other methods (e.g Finite Volume).

Concerning the properties of the final discretization, the two main tools of analysis are a consistency property w.r.t. a variational form evaluated for smooth solutions, and a positivity preserving property rooted in the positive coefficient schemes theory (see AR:17,Dec:17 [START_REF] Abgrall | High order methods for CFD[END_REF][START_REF] Deconinck | Residual distribution schemes: Foundations and analysis[END_REF] and references therein). For the schemes used here, high order of accuracy is formally ensured by the boundedness of the distribution matrices, and more generally of the test functions {w j } j∈K in ( eq:rd-split1 eq:rd-split1 [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF]. For the analysis of the depth positivity-preservation property of this Lax-Friedrich's-type scheme (as well as of the one obtained applying the limiter), the interested reader can refer to Ric:15 [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] for the Eulerian case. In the ALE case, an additional limiting step is required, as proposed in ArR:17 [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF]. This is to control the mismatch between the update in the mass equation, performed in terms of η, and the one of the bathymetry based on ( eq:batale eq:batale [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF]) and not on the ALE remap.

Well-balancedness, DGCL and mass conservation

The first property we can easily show is the following.

Proposition 2. (DGCL and well-balanced on the sphere) For any given mesh displacement on the sphere, the residual distribution prototype (

RC-corr RC-corr 47)-( phi-t phi-t 51
), with the quadrature strategy defined in Proposition 1

and by ( eq:aleterm eq:aleterm [START_REF] Chen | Second-order accurate Godunov scheme for multicomponent flows on moving triangular meshes[END_REF]) and ( eq:hyd eq:hyd 57), and with the distribution functions ( eq:RD eq:RD

58

) verifies the C-property on moving and fixed meshes, as well as for Cartesian and curvilinear coordinates.

Proof. The proof is obtained by applying the scheme to the initial solution w 0 = [η 0 , 0, 0] t . With the hypotheses made, and by virtue of ( eq:aleterm eq:aleterm [START_REF] Chen | Second-order accurate Godunov scheme for multicomponent flows on moving triangular meshes[END_REF]) combined with ( eq:delta-K eq:delta-K 31), and using ( eq:hyd eq:hyd 57) with η constant over the mesh, we deduce immediately that ΦK = 0, thus ΦK i = 0, and w * i = w 0 i .

In the corrector step the element residual, evaluated with constant η and invoking ( eq:hyd eq:hyd 57), reduces to

Φ K = |K| n+1 -|K| n ∆t w 0 - j∈K υ fj w 0 = 0 ( 60
)
due to the definition of υ f . Combined with the previous results and ( eq:RD eq:RD

58

) this implies w n+1 = w 0 .

This proposition corroborates the choices made in the previous sections concerning the form of the SWEs. As already mentioned, however, the choice of w as a main unknown has a drawback concerning mass conservation. Indeed if we sum both ( RC-corr RC-corr

47) and (

RC-pred RC-pred 50) over all the mesh nodes we obtain

i |C n+1 i |w n+1 i - i |C n+1 i |w * i = -∆t i K i Φ K i = -∆t K i∈K Φ K i = -∆t K Φ K (61)
Using ( phi-corr phi-corr

48) and (

eq:vol-int eq:vol-int

52) we can further modify the mass balance to

i |C n+1 i |w n+1 i - i |C n i |w n i = - ∆t 2 K Φ K (w n h , b n h ) - ∆t 2 K Φ K (w * h , b n+1 h ) (62)
If we only consider the mass equation, the right hand side gives

- ∆t 2 K ∂K √ G(hu j -ησ j ) n h n j + ∂K √ G(hu j -ησ j ) * h n j = 0 ( 63
)
face integrals of neighbor elements canceling each other. As a consequence we end with

i |C n+1 i |h n+1 i = i |C n i |h n i - i |C n+1 i |b n+1 i - i |C n i |b n i (64) mass-error

As remarked in

ArR:17 [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF] the error in total mas conservation is related to the error in satisfying a discrete ALE remap of the bed. To minimize the mass error, we use here the same technique proposed in ArR:17 [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF]. We set: [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF] eq:batale where b(χ) is a high resolution (analytical, or reconstructed using a finer mesh) description of the topography. We then evaluate the last integral using a higher order quadrature as done in ArR:17 [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF]. In particular, if the quadrature is of order κ we have, considering a smooth enough manifold:

|C n+1 i |b n+1 i := C 0 i J n+1 A √ G b(χ) dχ (
i |C n+1 i |b n+1 i = i Ci J n+1 A √ G b(χ) dχ + i O(h κ+3 ) = M 2 √ Gb(x) dx + O(h κ+1 ) ( 66
)
where the O(h κ+3 ) is the local error modeled as a polynomial of degree (at least) κ + 1, readily obtained using e.g. a truncated Taylor series, and accounting for the fact that |K| = O(h 2 ). The second equality uses ( eq:mesh-smooth eq:mesh-smooth [START_REF] Kuhnlein | Modelling atmospheric flows with adaptive moving meshes[END_REF]. Note now that the first term in the final expression is actually time independent as long as

M 2 = M 2 0
, which is in practice the case. Since also the integrand is a time independent quantity, we have

i |C n+1 i |b n+1 i - i |C n i |b n i = O(h κ+1 ) (67)
which combined with ( mass-error mass-error 64) shows that we can conserve h within an arbitrary accuracy.

Adaptive mesh movement strategy sec:mmpde

The mesh adaptation technique used here is based on an error dependent mesh movement. The main principle is to define a mapping A(χ, t) from a reference/initial configuration, described by the coordinate χ, to a physical mesh described by coordinate x, such that the error associated to a certain density is equi-distributed. The error density is often referred to as a monitor function m(x)

DeB:73 [START_REF] De Boor | Good Approximations by Splines with Variable Knots II[END_REF]. The differential form of this principle is referred to as Moving Mesh Partial Differential Equation or MMPDE for short.

The simplest method in the multidimensional case is the elliptic generator based on the MMPDE:

∂ ∂χ i m ∂x α ∂χ i = 0 α = 1, 2 (68) eq:mmpde-ch
This PDE is the Euler-Lagrange equation corresponding to the minimization of the energy functional CeH:01 [START_REF] Ceniceros | An efficient dynamically adaptive mesh for potentially singular solutions[END_REF]:

E(x) = 1 2 Ωχ m||∇x|| 2 dχ = 1 2 Ωχ m ∂x α ∂χ i ∂x α ∂χ i dχ . ( 69
) eq:efun-ch
It consists in a set of decoupled quasi-linear elliptic equations mapping each coordinate. Many efficient numerical methods are available to solve ( eq:mmpde-ch eq:mmpde-ch 68). The counterpart of its simplicity is that special care must be taken to avoid the degeneration of the quality of the mesh, or mesh tangling. If carefully implemented, the MMPDE ( eq:mmpde-ch eq:mmpde-ch 68) can be quite effective (see, e.g.

TaT:03,JIN200668,ChT:08,XuN:13,Zho:13,ZhC:13,PAN2015207 [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF][START_REF] Jin | An adaptive grid method for two-dimensional viscous flows[END_REF][START_REF] Chen | Second-order accurate Godunov scheme for multicomponent flows on moving triangular meshes[END_REF][START_REF] Ni | Remapping-free ALE-type kinetic method for flow computations[END_REF][START_REF] Zhou | An adaptive moving finite volume scheme for modeling flood inundation over dry and complex topography[END_REF][START_REF] Zhou | A well-balanced stable GRP scheme for shallow water equations for adaptive unstructured triangular meshes[END_REF][START_REF] Pan | Generalized coordinate transformation and gas-kinetic scheme[END_REF] and references therein).

The design of moving mesh methods on the sphere is a relatively recent subject of research. Optimal transport method have been discussed in WeB:16,McR:18 [START_REF] Weller | Mesh adaptation on the sphere using optimal transport and the numerical solution of a Monge-Ampere type equation[END_REF][START_REF] Mcrae | Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements[END_REF] using finite differences or finite-elements to solve the Monge-Ampere type nonlinear PDEs involved. The results shown in terms of mesh quality and mesh refinement are very interesting. It is however unclear whether the error reduction compensates the CPU time overhead of mesh adaptation. Based on our previous experience, here we have opted for a simpler version of mesh movement on the sphere, trying to remain as much as possible close to the simplicity of ( eq:mmpde-ch eq:mmpde-ch 68). It is known that functional ( eq:efun-ch eq:efun-ch 69) can be generalized to model a mapping A :

S 2 A0 → S 2 A x = A(χ, t),
from a computational sphere with metric tensor G ij and local coordinates χ to a physical sphere with metric H αβ and local coordinates x:

E(x) = 1 2 Ωχ mG ij H αβ ∂x α ∂χ i ∂x β ∂χ j √ Gdχ . ( 70
) eq:hm-e
Unfortunately, the above functional is similar to the harmonic map energy of DiL:06 [START_REF] Di | Moving mesh methods for singular problems on a sphere using perturbed harmonic mappings[END_REF] where it is shown that the minimization problem associated may lack a unique solution. This is what has led several authors to consider more complex techniques DiL:06,WeB:16,McR:18 [START_REF] Di | Moving mesh methods for singular problems on a sphere using perturbed harmonic mappings[END_REF][START_REF] Weller | Mesh adaptation on the sphere using optimal transport and the numerical solution of a Monge-Ampere type equation[END_REF][START_REF] Mcrae | Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements[END_REF] To obtain a simpler formulation, we propose to use a mapping form the reference sphere to the parametrization of the adapted one. So we map the sphere to a plane A : S A0 → E 2 x = A(χ, t). This corresponds to a flat target metric is H αβ = δ αβ and to an energy functional: [START_REF] Loevholt | Modeling propagation and inundation of the 11 March 2011 Tohoku tsunami[END_REF] eq:hm-e2 having used the fact that G ij is diagonal. In practice, all goes similarly to when mapping two planes.

E(x) = 1 2 Ωχ mG ij ∂x α ∂χ i ∂x α ∂χ j √ Gdχ = 1 2 Ωχ mG ii ∂x α ∂χ i ∂x α ∂χ i √ Gdχ (
The catch is that when minimizing ( eq:hm-e2 eq:hm-e2

71), we will obtain an adapted mesh which will depend on the parametrization. In our opinion this is not a problem as long as the same parametrization is used for the flow equations. The MMPDE obtained in this case is

1 √ G ∂ ∂χ i m √ GG ii ∂x α ∂χ i = 0 ( 72
) eq:mmpde-chs which, as in the Cartesian case, constitutes a decoupled system of PDEs in which the only non-linearity is associated to the monitor function m = m(x) which will depend on the solution (as in the Cartesian case). We use a straightforward P 1 Finite Element method to discretize the MMPDE, and a Newton-Jacobi iteration to update nodal displacements reading

δ k+1 i = δ k i - K i a K ii -1 K i j∈K a K ij (x k j ) K , x k+1 i = x n i + δ k+1 i , (73)
with a K ij the entries of the stiffness matrix arising from a straightforward finite element approximation of ( eq:mmpde-chs eq:mmpde-chs 72), with the only source of nonlinearity being the monitor function m = m(x). In some of the computations a standard Laplacian mesh smoothing step is also added HuR:97,BuH:09 [START_REF] Huang | Analysis of moving mesh partial differential equations with spatial smoothing[END_REF][START_REF] Budd | Adaptivity with moving grids[END_REF].

To test the moving mesh algorithm we use a benchmark from

McR:18 [START_REF] Mcrae | Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements[END_REF]. The monitor function is: In this paper we use the arc-length type monitor function of Win:67 [START_REF] Winslow | Numerical solution of the quasi-linear Poisson equation[END_REF]:

m = 1 + 10 sech 2 5 ||x -x 1 || 2 - π 2 4 + 10 sech 2 5 ||x -x 2 || 2 - π 2 4 (74) with x 1,2 = (± √ 3 2 , 0, - 1
m = 1 + α i f 2 i * , f i * = min 1, f i β max x∈Ω f i (75) eq:omega
In the tests of section sec:exp sec:exp 8, the physical quantities f i have been chosen depending on the phenomenon under study. For gravity waves, we have selected, as in ZhC:13 [START_REF] Zhou | A well-balanced stable GRP scheme for shallow water equations for adaptive unstructured triangular meshes[END_REF], only one term f 1 = ||∇η|| the free surface gradient, in order to detect free surface wave patterns and bore development. In presence of wet/dry transitions as in ArR:17

Moving Mesh PDE-ALE algorithm

For adaptive computations, the MMPDE-ALE algorithm starts with a triangulation at time n and the solution w n h . Setting the initial conditions for the MMPDE m 1 h = m n h and x 1 = x n , where m n h is the monitor function evaluated with w n h , we proceed as follows.

DO k=1,kmax

Step 1. Move the mesh according to the Newton-Jacobi iteration (69).

Step 2. Use the ALE schemes to get a remap on the new mesh as discussed e.g. in ArR:17 [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF] and re-evaluate the monitor function using the new solution

ENDDO

Step 3 Let x n+1 = x kmax+1 define the new triangulation, evolve the underlying balance law in ALE framework on the half time step mesh as discussed in section sec:RD sec:RD 6

Step 5. IF (t = T final ) EXIT ELSE GO TO Step 1.

As remarked in

ArR:17 [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF], the remap from one mesh to the new one here is only needed to evaluate the error function. This allows to simplify the mapping if necessary. The interested reader can refer to ArR:17

[1] for a discussion.

Numerical experiments

sec:zon

In this paragraph we test only the Eulerian (fixed grid) RD on some classical benchmarks taken from the test suite of Wil:92 [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF]. Grid convergence studies gave been conducted on 5 meshes generated with Gmsh and with reference sizes starting from h K R -1 = 0.16 and halved from one mesh to the other. The meshes thus obtained are composed of 1306, 4990, 19568, 75478, and 300426 triangles respectively.

Relative errors are evaluated following

Wil:92 [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF] as:

l 2 = ||e|| L2 ||h ex || L2 (76)
||e|| L2 = 2π 0 π/2 -π/2 (h(x) -h ex ) 2 √ G dx (77)
Case #2 is an exact steady geostrophic equilibrium allowing to measure the order of accuracy in presence of Earth rotation. The velocity ( given directly in covariant basis) and height fields are initially given by:

h(x, 0) = h 0 - 1 g ΩR u 0 + u 2 0 2 (-cos λ cos ϕ sin α + sin λ cos α) 2 u 1 (x, 0) = u 0 (cos α + cos ϕ tan λ sin α) (78)

fig:zon-err

ig:zon-icoerr

Case #5 is a perturbation of the previous one. The initial velocity and height fields are given by ( eq:wil-2 eq:wil-2

78)

with α = 0, h 0 = 5960 m and u 0 = 20 m/s. An isolated mountain is now added with the profile:

b = b 0 1 - r R (79) 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 mesh of the gmsh suite (h K R -1 = 0.02). 6 provides the contour lines of the water height after 5, 10, and 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 , (80)
Note that mass error is only related to the change of the coordinate system across the edges that delimit the polar cap. Regarding energy, none of the schemes proposed have 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 a linear approximation. The resulting error evolution is reported on figure

Advection of cosine bell

To assess the numerical accuracy of our algorithm for unsteady problems and in presence of mesh movement, we revisit test #1 of the Williamson suite Wil:92 [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF]. For this test, the advective component is tested in isolation: only the mass equation is a prognostic equation while the wind (always expressed in

u 1 (x, 0) = u 0 (cos α + cos ϕ tan λ sin α) (81)
u 2 (x, 0) = -u 0 sin ϕ sin α (82)
We compute the transport of a cosine bell given by:

h(x, 0) =      h 0 2 (1 + cos(3r)) + h 0 if r < π/3 h 0 otherwise (83)
where h 0 = 1000, u 0 = 2πR 12days , and r is the great circle distance from the center of the bell (ϕ C , λ C ):

r = R arccos (sin λ C sin λ + cos λ C cos λ cos(ϕ -ϕ C )) (84)
We have tested two different wind directions α = 0, and α = π/2 for which the bell is transported through the north and south poles, crossing four times a polar cap interface. To test the accuracy of the ALE formulations we have added an independent unsteady motion of the mesh nodes according to the transformation:

ϕ(t) = Φ + 0.5 cos (Φ) cos (Λ) sin 4πt T (85) λ(t) = Λ + 0.5 cos (Φ) cos (Λ) sin 4πt T ( 86
)
with Φ, Λ the computational or reference coordinates. The typical deformations obtained can be seen on 8. This motion will act as a background perturbation to the cosine bell advection.

The mesh convergence for fixed and moving meshes is reported on figure

fig:cb-conv fig:cb-conv 9.
The expected second order of accuracy is recovered. Moreover, we see roughly the same level of error for the two configurations (zonal and meridional advection). This suggests that the polar patching does not affect significantly the error constant, nor the slope. We also observe that the unsteady mesh distortion does not spoil neither the order of accuracy nor the error level. This validates the ALE formulation, at least in terms of accuracy

Circular hump

This is a test proposed in

Ros:04 [START_REF] Rossmanith | A wave propagation algorithm for hyperbolic systems on curved manifold[END_REF]. The initial condition is a circular depth disturbance at the equator:

h(x, 0) =    2 if arccos(cos(x 1 ) cos(x 2 )) ≤ 0.2 0.2 otherwise , u(x, 0) = 0 (87)
This initial condition is symmetric about the point (x 1 , x 2 ) = (0 • , 0 • ), and should remain symmetric in absence of rotation. We run computations on two half sphere meshes generated with gmsh: a coarse one with 7122 points and 14034 elements; a fine one with 39699 and 78900 elements. Adaptive computations on the coarse grid are performed using the MMPDE of section sec:mmpde sec:mmpde 10. We can see on all meshes good shock capturing and good symmetry in all radial directions. We can observe that the mesh adaptation itself also preserves a good symmetry of the grid, and allows an extremely sharp resolution of the inner and outer shocks. We can see the sharp capturing of the discontinuities on all grids, and the improvement brought by the MMPDE adaptation on the coarse one. In these plots the adaptive mesh solution is indistinguishable from the fine mesh one.

Finally on table

tab:ch-cpu tab:ch-cpu 2 we report the CPU times for all computations. We compute the overall time and also isolate the time associated to the MMPDE part (both mesh equation, and solution projections). We can see that the moving adaptive mesh computations are twice less expensive than the fine mesh ones, and also that even with our simplified formulation adaptation counts for non-negligible component of the simulation, mode than half for the RD scheme.

We repeat the same test adding the effects of the rotation. We set the dimensionless rotation rate to Ω = 5. In figure

fig:chc-iso fig:chc-iso
12 we report snapshots of the depth contour lines on fixed (coarse and fine) and adaptive meshes. As in the previous case, we see the higher resolution obtained with the adaptive moving mesh and the ALE formulation. This time mesh nodes are clustered both around the shock waves and also in correspondence of smoother features where high gradients appear. On table tab:chc-cpu tab:chc-cpu 3 we report the CPU times. tab:ch-cpu

Barotropic instability

This test, contained in Gal:04 [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 an instability. Here we test the RD scheme on a fine mesh with reference size of h K R -1 = 0.01, and compare with adaptive simulations on coarse grids with h K R -1 ≈ 0.02.

The fine mesh has roughly 300k triangles, while the coarse ones have approximately 80k triangles. For the adaptive simulations we compare results on the icosphere as well as on a gmsh mesh. In particular, adaptation has been performed using in the monitor function relative vorticity and its gradient through the formula of section subsec:mon subsec:mon [4]) and it is considered here as a reference solution. The last row on the bottom show the same fixed mesh computation on the coarse grid. We can clearly see that the onset of the instability as well as its evolution are wrongly reproduced on this mesh level. Note that the result is the same independently on the mesh topology, the icosahedral mesh providing very similar results to the one generated with gmsh.

Second and third rows report the adaptive simulations on the coarse meshes. We can clearly see that the adaptive mesh movement allows to correctly reproduce the global instability process, with a slight smoothing of last part of the instability (ϕ ≤∼ -180 • ) for the gmsh mesh, while the result on the coarse The complex interaction between incoming and trapped waves that follows is also well reproduced, see

others published results

Che:14,LoK:12,ShS:12 [START_REF] Chen | The March 11, 2011 tohoku M9.0 earthquake-induced tsunami and coastal inundation along the Japanese coast: A model assessment[END_REF][START_REF] Loevholt | Modeling propagation and inundation of the 11 March 2011 Tohoku tsunami[END_REF][START_REF] Shimozono | Propagation and inundation characteristics of the 2011 Tohoku tsunami on the central Sanriku coast[END_REF]. We remark that the adaptive computations allow to reproduce the fine mesh results with half of the mesh elements.

We now consider the approximation of the flooding of the Japan coastal prefectures. We focus on three 22 showing the meshes adapted to the initial coastline and to the coastline after 40 minutes. We can see here the limitations of the moving mesh approach which is stretched here to its maximum. For some of the features, as for example the thinner lands close to the Osimaseto Strait and those around Hirota bay, clearly the density In both pictures the flooded computational areas are represented in green, while red symbols denote the surveyed inundation in the real event. The fixed fine mesh computation gives a good prediction of the runup process, with exception of some fine detail as the inundation of two narrow channels at the end of Hirota bay, which remain dry throughout the simulation.

The result of the MMPDE-ALE scheme is acceptable as it does cover most of the surveyed inundation areas. However, there is no improvement w.r.t. the fine mesh, and we clearly see the main limitation of the approach: the effect of the mesh stretching introduced in the dry region close to the wet/dry areas (cf. right picture on figure

fig:floo-2 fig:floo-2 22
) is that as soon as a point becomes wet, the inundated region covered geographically is unrealistically large. To cope with this phenomenon one could either look for a recipe to a-priori tune the density of mesh nodes, so that all features may be resolved. One may also try to move more points in the dry region to improve the resolution, however for such multi scale problems, a more realistic path is to combine the MMPDE with a re-meshing step. We think the result obtained here is a good example of the how much MMPDE can be stretched for this type of problems.

Conclusions and outlook

We have discussed approximations of the rotating shallow water equations on the sphere using fixed and moving meshes. The main principles of the ALE form of the equations have been derived, and the relations between the GCL and well balancedness have been shown. A numerical approximation based on residual distribution has been proposed, with a numerical closure allowing to combine DGCL and Cproperty. A simplified elliptic PDE for mesh adaptation has been proposed and tested as a basis for a time dependent MMPDE-ALE method. The numerical results have shown that on fixed unstructured meshes the RD scheme proposed achieves accuracy close to te one usually presented in literature on structured grids. For flows with discontinuities and unstable fronts the MMPDE-ALE adaptation allows to enhance this accuracy, allowing in principle a considerable reduction in CPU time. For the cases considered here this reduction goes from 30% to 50%. We have also shown that for genuinely multi-scale problems, as the Tohoku-Honsu tsunami computation of section sec:toho sec:toho 8.5, this technique reaches its limitations and a re-meshing phase is clearly necessary.

This work is a first step toward several generalizations which will involve both higher order approximations based on stabilized continuous polynomial approximations, an ad hoc combination of mesh movement and local remeshing, and possibly improved treatments of the pole problem, following e.g. the finite element mixed 3D/2D approximation proposed in Bern:09 [START_REF] Bernard | High-order discontinous Galerkin schemes on general 2D manifolds applied to the shallow water eqautions[END_REF].

σ i ∂x i ∂χ j    , J -1 A =    1 ∂t ∂x j ∂χ i ∂t ∂χ i ∂x j    =    1 0
-Σ i ∂χ i ∂x j    (A.9) eq:ale_jac

In the above expression we have used the definitions of covariant mesh velocity, in actual covariant components σ i and in reference ones Σ i . The relationships between the two is given by Σ i = ∂χ i ∂x j σ j . One can show that the determinant of the ALE Jacobian reduces to the determinant of the spatial Jacobian Finally we distinguish tensor quantities in actual covariant components T = T ij g i g j from tensor quantities in reference covariant components T = T ij γ i γ j .

J A = det( ∂χ i ∂x j ).
e.g. T i0 = hu i . Then the divergence in curvilinear reference/ALE coordinates writes

∇ • T = 1 J ∂ ∂t χ J hu i + 1 J ∂ ∂χ j J T ij + Γi j0 hu j + Γi jk T jk (B.1) eq:div2
Note that the total Jacobian J comes from the composition of the spherical and of the ALE transformations, J = J A J F with determinant J = J A √ G. Moreover we have separated the Christoffel symbols related to time and space variation of the basis vectors.

The last step consists in transforming each { hu i , T i1 , T i2 } into actual components {hu i , T i1 , T i2 } through the inverse Jacobian J -1 A ( eq:ale_jac eq:ale_jac A.9)

hu i = hu i (B.2) T ij = ∂χ j ∂x k T ik -Σ j hu i (B.3)
If we replace the transformed variables ( eq:vartrans eq:vartrans B.3) into ( eq:div2 eq:div2

B.1) we note that ALE terms related to Christoffel symbols cancel out. This leads to the ALE-SWEs in the reference parametric coordinates: which appear in the same form of PS:03 [START_REF] Prusa | An all-scale anelastic model for geophysical flows: dynamic grid derformation[END_REF] except that momentum components are covariant and not physical.

1 J A √ G ∂ ∂t χ J A √ G hu i + 1 J A √ G ∂ ∂χ j J A √ G ∂χ j ∂x k T ik -J A √ G Σ j hu i + Γi
To get the counterpart of ( eq:ALE-SWEs-ref eq:ALE-SWEs-ref B.4) into actual components it suffices to apply chain rule to each term.

The advective term then reads

∂ ∂χ j J A √ G ∂χ j ∂x k T ik = J A ∂χ j ∂x k ∂ ∂χ j √ GT ik + √ GT ik ∂ ∂χ j J A ∂χ j ∂x k = J A ∂ ∂x k √ G T ik (B.5)
The second term on the right hand side is set to zero thanks to the metric identity. The ALE term goes similarly:

∂ ∂χ j J A √ G Σ j hu i = J A ∂χ j ∂x k ∂ ∂χ j √ G σ k hu i + √ Gσ k hu i ∂ ∂χ j J A ∂χ j ∂x k = J A ∂ ∂x k √ G σ k hu i (B.
ψ h + ψ j 6 (f h (x) -f j ) + (ψ h -ψ j ) f h (x) + f j 6 √ G(x) - √ G(x * ) dx
We can now use standard arguments to bound this quantity as

≤ K j∈K K |ψ h + ψ j | 6 |f h (x) -f j | + |ψ h -ψ j | |f h (x) + f j | 6 √ G(x) - √ G(x * ) dx ≤ c 1 |S 2 | ψ L ∞ sup j ∂ x j f L ∞ + sup j ∂ x j ψ L ∞ f L ∞ sup j ∂ x j √ G L ∞ h 2
To recover ( RC-pred RC-pred 50) we use the second in ( eq:gal-mass-d eq:gal-mass-d [START_REF] Ronchi | The "cubed sphere": A new method for the solution of partial differential equations in spherical geometry[END_REF], to replace the lumped mass matrix with |C i |.

The corrector step is derived in a similar manner. We start from the second stage Runge-Kutta iteration for the unstabilized Galerkin scheme: As proposed in RiA:10,ArR:14 [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Arpaia | An ALE formulation for explicit Runge-Kutta residual distribution[END_REF] , we stabilize the scheme using the approximate shifted residual in which w n+1 is replaced by w * :

S 2 ϕ i √ GJ n+1 A w n+1 dχ - S 2 ϕ i √ GJ n A w n dχ+ ∆t 2 S 2 ϕ i JA ∂ x j (F j - √ Gσ j w) n + √ GS n dχ + ∆t 2 S 2 ϕ i JA ∂ x j (F j - √ Gσ j w) * + √ GS * dχ + K i K0 γ i √ GJ n+1 A w * - √ GJ n A w n + ∆t 2 JA ∂ x j (F j - √ Gσ j w) + √ GS n + ∆t 2 JA ∂ x j (F j - √ Gσ j w) + √ GS * dχ = 0
(D.6) eq:a6

We now add and remove the term For piecewise constant (per element) values of γ i , we can see that in general the form of the full mass matrix on a manifold is

m K ij = (m GAL ij ) K + K j γ i
If proposition 1 is valid, the last expression becomes

m K ij = |K| δ ij + 1 12 + |K| 3 γ i
The standard RD mass matrix corrections are obtained by setting γ i = β K i -I/3, while the final form ( RC-corr RC-corr 47) of the scheme is obtained by recalling the hypothesis RiA:10

[46]:

Φ K i = K0 (ϕ i + γ i ) √ GJ n+1 A w * - √ GJ n A w n ∆t + 1 2 JA ∂ x j (F j - √ Gσ j w) + √ GS n + 1 2 JA ∂ x j (F j - √ Gσ j w) + √ GS * dχ

  sec:ale sec:ale 3 and sec:ale-SW-sphere sec:ale-SW-sphere 4 are then devoted to the ALE-SWEs formulation on the sphere. The main notation and ingredients required for the discrete approximation are the introduced in section eq:not-as eq:not-as 5. These are used in section sec:RD sec:RD 6 to construct well balanced second order residual distribution. The mesh movement strategy is discussed in section sec:mmpde sec:mmpde 7, and section sec:exp sec:exp 8 discusses an extensive evaluation of the schemes. The paper is ended by a summary and an outlook on future developments.2. SWEs on the sphere sec:SW-sphere

Figure 1 :

 1 Figure 1: Swept area decomposition of area variations for mesh element K. fig:swept-area

  figure fig:swept-area fig:swept-area

2 ) 18 [ 59 ]

 21859 fig:t5-t6 fig:t5-t6 2. The mesh is highly compressed around the bands. The result is quite close to those presented in McR:18 [59], although clearly the optimal transport method of the reference gives a smoother mesh. Somewhat nicer results are obtained when including Laplacian smoothing, as shown on the right picture on the same figure.

Figure 2 :

 2 Figure 2: Adaptation to the monitor function ( eq:test_mmpde eq:test_mmpde 74). Left: computational (icosahedral) mesh. Center: adaptation without smoothing. Right: adaptation with Laplacian smoothing with smoothing coefficient 1/2.

[ 1 ] 1 = 17 [ 1 ]

 11171 we have added a second term f 2 = F (h) where the function F is a regularized Heaviside function going from 1 to zero in vicinity of dry regions. When considering barotropic instability we have taken f ||∇ × u|| and f 2 = ||∇∇ × u|| in order to follow roll-up of the vorticity field. For all other details on our implementation of the MMPDE method please refer to ArR:and references therein. The values of the tuning coefficients α i , β are given for each test in section sec:exp sec:exp 8.

8 . 1 .

 81 fig:sphere-mesh fig:sphere-mesh 3. Adaptation is performed with 5 iteration of the Newton-Jacobi

Figure 3 :

 3 Figure 3: Computational grids on the sphere. Top row: mesh #1 generated with gmsh. Bottom row: third level of refinement for icosahedral mesh.

u 2

 2 fig:zon-err fig:zon-err4 the error maps after 5 days of simulation on the fourth grid (h K R -1 = 0.02). Particular attention should be put along the interface line (it is a discontinuity for the coordinates). Error iso-countours are very smooth for the RD scheme proposed, with seemingly no effect of the patching treatment used for the poles. For comparison we also report on figure fig:zon-icoerr fig:zon-icoerr

  tab:ch-cpu relative to the computation on the fourth-level icosahedral mesh (81920 nodes) which has roughly the same reference size. The error (||e|| L2 = 4.084761e -05 for α = 0, ||e|| L2 = 2.634469e -04 for α = π 4 )result is very close, showing no impact of the mesh topology.

Figure 4 :

 4 Figure 4: Steady zonal flow. Absolute value of local error on gmsh grid h K R -1 = 0.02. Black areas correspond to 2 m error, white areas to zero error. Left: α = 0. Right: α = π/4.

Figure 5 :

 5 Figure 5: Steady zonal flow. Absolute value of local error on icosahedral grid (forth-level refinement) h K R -1 ≈ 0.02. Black areas correspond to 2 m error, white areas to zero error. Left: α = 0. Right: α = π/4.

Figure

  Figure fig:mou-1 fig:mou-1

2 7.

 2 fig:mou-2 fig:mou-2 7. The evolution obtained, as well as errors obtained, for the mass and energy are comparable to or better than those typically reported in literature for computations on structured grids, and with the high resolution results by the German Weather Service (DWD) http://icon.enes.org/swm/stswm/node5.html.

Figure 6 :

 6 Figure 6: Zonal flow over an isolated mountain. Snapshot of the numerical solution for fixed grid simulation on gmsh mesh h K R -1 = 0.02 at day 5,10 and 15. Contour levels are from 5050 m to 5950 m in intervals of 50 m.

Figure 7 :

 7 Figure 7: Zonal flow over an isolated mountain. Relative conservation error for mass (left) and energy (right) on gmsh mesh h K R -1 = 0.02.

  figure fig:cb-mesh fig:cb-mesh

Figure 8 :

 8 Figure 8: Advection of cosine bell. Transformed meshes number 2 at three different instants. At day 12 (final time) the transformed mesh coincides with the computational mesh.

Figure 9 :

 9 Figure 9: Advection of cosine bell: grid convergence. Left: α = 0. Right: α = π/2 fig:cb-conv

7

  fig:ch-iso fig:ch-iso

Figure 10 :

 10 Figure10: Circular hump on a sphere: depth iso-contours (RD solutions, 20 levels between 0 and 0.55) at times t = 0.6 (top) and t = 0.9 (bottom). From left to right: fixed coarse mesh solution, fixed fine mesh solution, ALE adaptive mesh solution, adaptive mesh.

Figure 11 :

 11 Figure 11: Circular hump on a sphere. Solution along the equator.

Figure 12 :Figure 13 :

 1213 Figure 12: Circular hump on a rotating sphere: depth iso-contours (RD solutions, 20 levels between 0 and 0.55) at times t = 0.8 (top) and t = 1.2 (bottom). From left to right: fixed coarse mesh solution, fixed fine mesh solution, ALE adaptive mesh solution, adaptive mesh.

17 [ 67 ] with the approach of SaF: 13 [ 68 ]

 17671368 fig:gal1 fig:gal1 15, showing a nice

Figure 14 :

 14 Figure 14: Barotropic instability. Vorticity field at day 6. Contour lines are from -1.1e -4 s -1 to -1.5e -4 s -1 in intervals of 2e -5 s -1 . Top: fixed fine mesh. Second and third rows: adaptive (respectively icosahedral and gmsh) coarse mesh. Bottom: fixed coarse.

Figure 15 :

 15 Figure 15: Barotropic instability. Relative vorticity field with corresponding adapted mesh at day 5, 5.5 and 6.

Figure 16 :

 16 Figure 16: Tohoku-Honsu tsunami: domain with position of the GPS probes (left) and initial wave height (right). Colour legend in meters.

Figures

  Figures fig:ml-10 fig:ml-10 17, fig:ml-20 fig:ml-20 18, and fig:ml-30 fig:ml-30 19 provide visualizations of the first 10, 20, and 30 minutes of propagation, showing water height and the adaptive meshes obtained in vicinity of the Japanese coast. We can see how the

Figure 17 :

 17 Figure 17: Tohoku tsunami. Snapshot of the free surface level (close up of the Iwate and north of the Miyagi prefectures, same colour legend of previous figures) at t = 0:10' (time origin is 11 March 2011 14:51:18 LT).

Figure 18 :

 18 Figure 18: Tohoku tsunami. Snapshot of moving mesh (close up of the Iwate and north of the Miyagi prefectures, same colour legend of previous figures) at t = 0:20' (time origin is 11 March 2011 14:51:18 LT).

Figure 19 :

 19 Figure 19: Tohoku tsunami. Snapshot of the free surface level (close up of the Iwate and north of the Miyagi prefectures, same colour legend of previous figures) at t = 0:30' (time origin is 11 March 2011 14:51:18 LT).

  fig:floo-1 fig:floo-1 21, from south to north. The figure shows on the left the fine mesh with the initial coastline superimposed, a visualization of the initial coastline in the middle and the position of the inundation front after 40 minutes from the beginning of the event. The figures give an idea of the geometrical complexity of the front line, with inundation scales of the order of magnitude of 10 meters, which can be compared with the propagation scales of several hundreds of kilometers. The handling of this complexity with the MMPDE is represented on figure fig:floo-2 fig:floo-2

Figure 20 :

 20 Figure 20: Simulated sea level displacement recorded at GPS buoys with fixed embedded mesh and moving embedded mesh. fig:ssh

Figure 21 :

 21 Figure 21: Tohoku tsunami. Flooding of the Iwate prefecture computed with RD scheme on fine fixed grid mesh. Left: zoom of the mesh, Right: snapshot of the inundated areas at at t = 0:00' and t = 0:40' (time origin is 11 March 2011 14:51:18 LT). fig:floo-1

Figure 22 :

 22 Figure 22: Tohoku tsunami. Flooding of the Iwate prefecture computed with RD scheme. Snapshot of the mesh at t = 0:00' and t = 0:40' (time origin is 11 March 2011 14:51:18 LT).

Figure 23 :

 23 Figure 23: Tohoku tsunami. Green zones are flooded areas in Iwate prefecture. Left: fix coarse simulation. Middle: fix fine simulation. Right: adaptive ALE simulation on coarse mesh. Red dots denotes surveyed points, from Che:14 [70].

  considered an ALE transformation in the parametric domain from the actual curvilinear coordinates x = {x 1 , x 2 } with covariant basis {g 1 , g 2 } to a set of reference/ALE curvilinear coordinates χ = {χ 1 , χ 2 } with covariant basis {γ 1 , γ 2 }:x = A(χ, t) (A.8)Two different time derivatives have been introduced: the Eulerian time derivative ∂ ∂t and the ALE one ∂ ∂t χ . Then, the spacetime Jacobian matrix and its inverse for the ALE transformation are

For= g i • ∂g j ∂x m ∂x m ∂t χ = Γ i jm σ m = Γi jk ∂χ k ∂x m σ

  the time-dependent ALE transformation, Christoffel symbols change with time. We can express them with respect to the actual components such as in ( eq:a1-christ eq:a1-christ A.5) or with respect to the reference/ALE components:Γi jk = g i • ∂g j (A(χ, t)) ∂χ k (A.10) the relationship between the two beingΓ i jm = g i • ∂g j ∂x m = g i • ∂g j (A(χ, t)) ∂χ k relativeto the time derivative (denoted by the index 0) can be expressed by chain rule Γi j0 = g i • ∂g j (A(χ, t)) ∂t χ m = Γi jk Σ k (A.13)

  jk ∂χ k ∂x m T jm = 0 (B.4)

  * ) dx We now use the relation (ab) 2 -(ab) 1 = (a 2 + a 1 )(b 2 -b 1 )/2 + (a 2 -a 1 )(b 2 + b 1 )/2 and recast the error as = K j∈K K

  ϕ i JA ∂ x j (F j -√ Gσ j w) n + √ GS n dχ + ∆t 2 S 2 ϕ i JA ∂ x j (F j -√ Gσ j w) * + √ GS * dχ = 0(D.5) eq:a5

S 2 ϕ i √ GJ n+1 Aw

 n+1 * dχ to obtain the following error correction form of the stabilized scheme w * )dχ = -

  a discrete analog of the GCL on moving grids (

				eq:gcl-s eq:gcl-s 11);
	3. the ALE remap of the bathymetry (	eq:b_ALE_remap eq:b_ALE_remap 16);
	4. the metric relations (	eq:Ric eq:Ric A.6) given in (	app:a1 app:a1

Table 1 :

 1 Steady zonal flow. l 2 error.

	5 the error plot

Table 2 :

 2 Circular hump on a non rotating sphere. CPU times.

Table 3 :

 3 Circular hump on a rotating sphere. CPU times.

tab:chc-cpu

[START_REF] Lesoinne | Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations[END_REF]

Table 4 :

 4 Barotropic instability. CPU times.

	tab:jet-cpu

Email address: l.arpaia@brgm.fr, mario.ricchiuto@inria.fr (Luca Arpaia and Mario Ricchiuto)Preprint submitted to Journal of Computational Physics December 4, 2019

Acknowledgments

Work partially funded by the TANDEM contract, reference ANR-11-RSNR-0023-01 of the French Programme Investissements d'Avenir. Some of the computations presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine. We are grateful to our colleagues from BRGM Orléans (and in particular R. Pedreros) for useful discussions and for the data used in the Tohoku simulations.

Appendix A. Notation and concepts from differential geometry app:a1

This appendix recalls some basic concepts from differential geometry, and used in the main body of the paper. We adopt a spacetime formalism which is better suited to time-deforming frame of reference.

In particular, we work with spacetime coordinates {x 0 , x 1 , x 2 }, where the first is understood to be time x 0 = t, while {x j } j=1,2 are purely spatial ones. In all the formulas that will follow, standard summation convention is implied, unless otherwise stated. We consider a steady transformation from the set of Cartesian coordinates r = {r 1 , r 2 , r 3 } to a set of curvilinear coordinates x = {x 1 , x 2 }. The curvilinear coordinates describe a sphere embedded in R 3 with covariant vector basis which are tangent to the sphere {g 1 , g 2 }. The transformation reads:

We can define the Jacobian matrix of this mapping, and of the inverse mapping, as

An important role is played by the metric tensor G, which is a symmetric tensor that allows to define a measure of length on the sphere ds = G ij dx i dx j , with

3) eq:Gij_def

The metric tensor can also be evaluated directly from the mapping Jacobian using the relation

The determinant of G will be denoted by G = det G. The determinant of the transformation then reads

We use the inverse metric tensor G -1 whose entries are denoted by G ij . We also introduce for later use the second-kind Christoffel symbols:

with g i the contravariant vector basis. We recall two important properties of the Christoffel symbols:

for the demonstration, see Ngu:14 [START_REF] Nguyen-Shafer | Tensor Analysis and Elementary Differential Geometry for Physicist and Engineers[END_REF].

Even though the above relations are valid on general 2D manifolds, all the numerical results in this paper are obtained for a latitude longitude (lat-lon or λ-ϕ) parametrization of the sphere for which the metric tensor and Christoffel symbols of second kind write respectively as:

with R the sphere radius.

Appendix B. ALE-SWEs on a 2D manifold derived with a spacetime approach app:a2

With the definitions of ( app:a1 app:a1

Appendix A), the derivation of the ALE-SWEs in curvilinear coordinates is immediate. We consider only momentum conservation since the mass equation is a particular case.

Momentum conservation is expressed in terms of the 4-divergence of a second order tensor T . We denote, as it is custom, the first row of the tensor as the momentum vector which represents the conserved quantity, Replacing ( eq:adv1 eq:adv1 B.5), ( eq:adv2 eq:adv2

B.6) and (

eq:chris eq:chris A.12) in ( eq:ALE-SWEs-ref eq:ALE-SWEs-ref B.4) we get the ALE-SWEs on a 2D manifold with respect to the actual parametric coordinates:

Except for the physical source terms, the above set of equations coincides with ( eq:mom-as eq:mom-as 13) which is the momentum equation discretized in this article. In section sec:ale-SW-sphere sec:ale-SW-sphere

4 it has been obtained by a more classical procedure combining transport formulas and the GCL.

Appendix C. Proof of proposition 1 app:a3

The validity of ( eq:dgcl-rels1 eq:dgcl-rels1

34), (

eq:dgcl-rels2 eq:dgcl-rels2

35) and (

eq:dgcl-rels3 eq:dgcl-rels3

36) as a consequence of the constant definition ( eq:div-dgcl eq:div-dgcl 37) is trivially shown using ( eq:delta-K eq:delta-K

31) and the definition of the integrated edge velocities (

eq:ivc-s eq:ivc-s

32)-(

eq:ivc-s1 eq:ivc-s1 33). Concerning the accuracy of the projection, we remark that by virtue of the mean value theorem

) eq:relG for some x * within the element. Estimating the error of ( eq:gal-mass-d eq:gal-mass-d

38

) is equivalent to estimating the error of

To this end, we follow the consistency analysis used e.g. in AR:17,Dec:17 [START_REF] Abgrall | High order methods for CFD[END_REF][START_REF] Deconinck | Residual distribution schemes: Foundations and analysis[END_REF]: we take the difference between the two expressions, sum over all the elements, multiply by the sampled nodal values of an arbitrary test function ψ ∈ C r (r ≥ 1), and sum over all the nodes. We end up with the following global error

The trick is now to use ( eq:relG eq:relG C.1) to remove nodal values ψ j f j , leading to

Stabilized FEM and RD in moving curvilinear coordinates app:a4

We sketch here the derivation of the predictor corrector scheme ( RC-pred RC-pred 50)-( eq:RD eq:RD

58

), and provide some details concerning the computation of the mass matrices with and without the simplified quadrature of proposition 1.

We start with the following explicit unstabilized predictor:

GS n dχ = 0 (D.1) eq:a1

which we recast as

(D.2) eq:a2

For a fixed mesh, J A = 1 and σ j = 0. For a moving mesh (and only for a moving mesh) we assume that the mass matrices are computed following proposition 1. Otherwise, on fixed meshes we can either use the approximate quadrature of the proposition, or exact (or higher order) quadrature of ( eq:sp-tr eq:sp-tr

17), (

eq:gal-mass eq:gal-mass

21)

and ( eq:gal-mass-L eq:gal-mass-L

). In any case, either by cancellation, or by virtue of ( eq:dgcl-rels3 eq:dgcl-rels3

36), we can recast the scheme as

3) eq:a3

As in

ArR:14 [START_REF] Arpaia | An ALE formulation for explicit Runge-Kutta residual distribution[END_REF] we 1) stabilize the predictor step using an incomplete residual which does not involve the time increment, and 2) lump the Galerkin mass matrix to obtain:

(D.4) eq:a4

where γ i is a properly defined stabilization bubble. In particular, following RiA:10,ArR:14 [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Arpaia | An ALE formulation for explicit Runge-Kutta residual distribution[END_REF] it is assumed that for RD schemes the definition of the right hand side is such that