Luca Arpaia

Keywords: Shallow Water Equations, Arbitrary Lagrangian Eulerian form, Finite Volume, Multidimensional Upwind Schemes, tsunami simulations, moving mesh method, r-adaptation Équation de Saint-Venant, Description Arbitraire Lagrangienne Eulérienne, Volume Finis, Méthode d'upwind multidimensionel, adaptation du

In this thesis we implement the Shallow Water equations (SWEs) on unstructured grids in order to simulate free surface ow over irregular Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

iii bathymetries, wetting/drying and other complex phenomena that typically occurs in hydrodynamic applications. In particular we would to accurately simulate tsunami events, from large scale wave propagation up to localized runup. To this aim we use two methods that are extensively compared along the manuscript: the Finite Volume method, which is very popular in the hydrodynamics and hydraulic community and a more recent technique called Residual Distribution which belongs to the class of multidimensional upwind schemes. To enhance the resolution of important ow feature such as bore development or small scale ooding, we use a dynamic mesh adaptation based on a redistribution of mesh nodes or r-adaptation (r stands for relocation).

The proper combination of this method with the ow solver is usually referred to as Moving Mesh Method. Among the many dierent moving mesh algorithms available we propose an Arbitrary Lagrangian Eulerian (ALE) form of the SWEs which elegantly permit to evolve the ow variables from one mesh to the updated one.

Résumé

Dans cette thèse, nous implémentons les équations de Saint Venant (SV), ou de Shallow Water, sur des grilles non structurées an de simuler des écoulements de surface libre sur des bathymétries irrégulières, incluant inondation et d'autres phénomènes complexes qui se produisent généralement dans des applications hydrodynamiques. En particulier, nous voudrions simuler avec précision les tsunamis, la propagation d'onde à grande échelle jusqu'à l'inondation très localisé. À cette n, nous utilisons deux méthodes qui sont comparées en profondeur le long du manuscrit: la méthode des volumesnis, très populaire dans la communauté hydrodynamique et hydraulique et une technique plus récente appelée Distribution du Résidu appartenant à la classe des schémas upwind multidimensionnels. Pour améliorer la résolution de certaines caractéristiques de l'écoulement telles que le développement du déferlement et les inondations à petite échelle, nous utilisons une adaptation de maillage dynamique basée sur une redistribution des n÷uds de maillage, aussi appelé adaptation de type r (r signiant "relocalisation"). La combinaison appropriée de cette méthode avec le solveur SV est généralement appelée Méthode de Maillage Mobile. Parmi les nombreux algorithmes de maillage mobile disponibles, nous proposons une forme Arbitrary Lagrangian Eulerian (ALE) des équation SV qui permettent de faire évoluer les variables de ux d'une maille à l'autre de manière élégante. Dans ce contexte, nous soulignons les principales contributions de la thèse:

1. Nous montrons l'importance de conserver toutes les propriétés standards d'un solveur Eulérien SWE tel que la préservation du lac au repos et la conservation de la masse également sur des maillages en mouvement.

2. Notre couplage ALE est comparé à l'approche de rezoning, avec une légère augmentation de la performance globale de l'algorithme en termes de précision et de temps CPU.

3. Nous étendons l'approche ALE sur la sphère an d'inclure l'eet de la courbure terrestre dans la dynamique de propagation des ondes à grande échelle du tsunami.

4. la simulation du tsunami 2011 de Tohoku-Honsu devrait prouver que la méthode de maillage mobile étudiée dans la thèse, bien que simple, pourrait être un bon candidat pour réduire le coût de calcul des simulations de tsunami.

Introduction

Context of the thesis: the TANDEM project

The 11 March 2011 an earthquake of M w 9.1 with epicenter o the coast of Tohoku triggered a tsunami on the eastern coast of the Honsu island. The devastating wave has hit the coastal prefectures of Iwate, Miyagi and Fukushima.

It was the tsunami attack to cause the major part of damages and fatalities (the Fire and Disaster Management Agency counted for 16,278 fatalities and 2,994 missing in Tohoku and Kanto regions). Among the damages, the accident of the Fukushima Daichi Nuclear Power Plant (NPP) was the one that mostly impressed the public opinion worldwide. In the Fukushima site, the designed tsunami height, whose computation was based on a standard procedure Yanagisawa et al. [2007], was 6.1 [m]. The 11 March 2011 a tsunami wave 16 [m] high hit the Fukushima Daichi NPP. The water ooded in the site, causing the loss of the cooling system which was followed by the hydrogen explosion of the reactors and by the subsequent propagation of radioactive material. In the days that followed Japan government declared nuclear emergency and thousands of people have been evacuated due to elevated radiation risk, evoking a scenario similar to the Chernobyl accident in 1986. This event mobilized the International Atomic Energy Agency (IAEA) which asked for new security standards concerning the NPP located along coastal regions. In France the design of new guidelines was particularly urgent: 19 NPP are active and 5 of them (Blaye, Flamanville, Paluel, Penly and Gravelines) are right on the Atlantic coast. In this context, the French National Research Agency (ANR) launched the TANDEM project (Tsunamis in the Atlantic and the English ChaNnel Denition of the Eects through numerical Modeling) which has funded the present thesis. The objective was to study the 2011 tsunami by numerical modeling and to dene the tsunami hazard for the French coastal regions. 10 French partners were involved: among them there were industrial companies (Electricté de France, Principia) and French research institutions (University of Pau, BRGM, SHOM, Ifremer, IRSN and INRIA). The project was structured in four parts or work-packages (WP):

• WP1 consists in a set of benchmarks to validate, compare and improve the dierent hydrodynamic softwares available among the partners. The idea was similar to the one of the Workshop on Long Wave Runup Model IIW-LWRM [1995]: collect a series of existing laboratory tests for long wave propagation and runup and check the ability of the hydrodynamic solver (both the physical model and numerical method) to well reproduce the experiments.

• WP2 studied the inuence of uncertain parameters in the computation of tsunami heights and ooded areas. Sensitivity analysis on the computational and physical parameters was conducted in order to provide a quantications of them. These includes mesh denitions as well as the seismic source.

• WP3 was dedicated to the simulation of the 2011 Tohoku tsunami which is widely recognized to be the largest instrumentally recorded earthquake in Japan history. This fact will permit us: 1) to recover realistic initial waveforms from seismic/tsunami data inversion 2) to compare the numerical simulation against high reliable data.

• WP4 provided an inventory of all the possible tsunamigenic sources for the French coast, the denition of possible tsunami scenario and their numerical simulation.

This manuscript contributes to WP1, WP2 and WP3. The main objective of the thesis consists in increasing the accuracy of hydrodynamic (and in particular tsunami) simulations through a mesh adaptation strategy, rather then employ very high order schemes. This choice is motivated by the fact that, for large scale tsunami simulations, static grids are inherently not ecient. The tsunami wave is localized in a certain area of the ocean, then crosses all the ocean basin, shoals, steepens near the coastlines and inundates the land. To enhance important features in tsunami simulations such as the resolution of wave patterns during propagation and shoaling or the runup/rundown stages in which ne scale ooding occurs, mesh adaptation could be the ideal candidate. In this thesis mesh adaptation will be achieved by redistributing (or relocating) the same number of mesh nodes within the domain.

State of the art

Tsunami science and tsunami modeling Introduction the derivation of analytical solutions. The physical model was quite simple: a long wave (a periodic wave, solitary waves or N-shaped waves) travels on a constant depth before reaching a sloping beach on which runup and reection occur. Even if simple, there were important unresolved questions such as the role of non-linearity and dispersion in the wave transformation. The simplest depth averaged model describing long wave propagation is the non linear Shallow Water equations (in this thesis abbreviated as SWEs). SWEs are derived from the conservation of mass and momentum under the hydrostatic assumption. It is assumed that the horizontal length scale of the ow is much larger then vertical one (the scale of the water depth); the consequence is that the vertical dynamics can be neglected. [START_REF] Carrier | Water waves of nite amplitude on a sloping beach[END_REF] resolved the SWEs for various waveforms climbing a sloping beach. Later [START_REF] Synolakis | The runup of solitary waves[END_REF] compared the analytical solution for a solitary wave with the experimental water level prole created in a wave tank: he found that the SWEs approximated very nicely the long (and non breaking) wave proles all along the runup process. It was also found that this set of equations correctly describes non linear wave transformation and wave breaking. The success of the SWEs made them the most commonly long wave approximation for runup calculation. It was also recognized that the lack of wave dispersion could give wrong wave propagation results in deep water and shoaling. To include these eects, the Boussinesqtype equations received considerable attention, see [START_REF] Peregrine | Long waves on a beach[END_REF]; [START_REF] Sorensen | A new form of the boussinesq equations with improved linear dispersion characteristics[END_REF]; [START_REF] Nwogu | An alternative form of the boussinesq equations for nearshore wave propagation[END_REF]. These eects are only relevant locally and for landslide tsunamis. Here we will consider the SWEs as our set of governing equations.

tested and analyzed against experimental data in order to obtain reliable numerical tools with precise range of applicability. It was also realized that tidal gauges where unt to record a tsunami wave. A new system of tsunamometers buoys to record tsunami height in deep ocean was implemented with an appropriate frequency response. In two thousands a tsunami forecasting system was operational in the Pacic Ocean. In 2003 there was the rst warning cancellation in Hawaii. After the megathrust earthquake of 2004 in Sumatra, numerical modeling was mature to robustly simulate the complete tsunami event, from generation to propagation and runup. Numerical models were able to explain the large runups and destruction observed in coastal Thailand Ioualalen et al. [2007]. When this thesis started in 2014, there were already many published works on the numerical modeling of the Tohoku-Honsu tsunami, see [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] Løvholt | Modeling propagation and inundation of the 11 march 2011 tohoku tsunami[END_REF]; [START_REF] Shimozono | Propagation and inundation characterstics of the 2011 tohoku tsunami on the central sanriku coast[END_REF].

Numerical methods in hydrodynamics

SWEs are of hyperbolic type, eventually with source terms to take into account dierent eects such as bathymetry and friction. The hyperbolic nature of the underlying equations have permitted to transfer earlier upwind Finite Dierence and Finite Volume (FV) methods, originally developed for the compressible Euler equations of gasdynamics in the aeronautical community, to hydrodynamics.

It was soon realized that there were specic issues. First and dierently from gasdynamics, there was the problem of embedding at a discrete level exact solutions more complex then the trivial uniform ow with constant depth. The lake at rest, which results from the physical balance between hydrostatic pressure and the pressure force played by the bottom on the uid, emerged as one of the most important state to be preserved numerically. This was a way to prevent the unexpected rise of numerical oscillation in regions where the ow is at rest. Such schemes are called Well Balanced (WB). The work of [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF] was pioneering in this area since it established the link between upwinding of the source term and Well Balanced property. This initial work in one dimension, has been led throughout the years to many dierent results allowing the construction of unstructured mesh discretizations verifying the WB via an appropriate coupling of the numerical ux and numerical source terms [START_REF] Hubbard | Flux dierence splitting and the balancing of source terms and ux gradients[END_REF], or based on dierent forms of the equations, as the well-balanced form of [START_REF] Russo | Central schemes for balance laws[END_REF]; [START_REF] Kurganov | Central-upwind schemes for the Saint-Venant system[END_REF]; [START_REF] Chertock | Well-balanced schemes for the shallow water equations with coriolis forces[END_REF], or the so-called pre-balanced form of [START_REF] Rogers | Mathemathical balancing of ux gradient and source terms prior to using roe's approximate riemann solver[END_REF]; [START_REF] Liang | Adaptive quadtree simulation of shallow ows with wet-dry fronts over complex topography[END_REF]; [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]. A review on Well Balanced can be found in [START_REF] Castro | Well-balanced schemes and path-conservative numerical methods[END_REF].

In coastal engineering, modeling the moving shoreline is a second important issue. At the beginning, the shoreline was treated as a boundary condition which is changing with time and it should be determined by the computations Introduction itself. An alternative way, that proved to be simpler, is to include the wet/dry interface in the computation, introducing the presence of cells of mixed type where both dry and wet nodes appears. This method was rst employed for Finite Dierence by [START_REF] Sielecki | The numerical integration of the nonlinear shallow-water equations with sloping boundaries[END_REF]. Many years after, a robust and eective treatment of wetting/drying was set in [START_REF] Brufau | A numerical model for the ooding and drying of irregular domains[END_REF]. In a Finite Volume context the situation with cells of mixed type can be treated as a Riemann problem for which one state is wet and one state is dry, the solution consisting of a rarefaction wave advancing in the dry region and causing ooding. However unphysical oscillation may appear at the front especially for the ow speed: the ow velocity in SWEs code is typically computed as the ratio between discharge and depth, both small quantities at the shore, and this lead to not accurate values. A depth threshold value is typically xed avoiding the computation of large velocities. Additional problems arise when bed slope is present. It was recognized the diculty of recovering the lake at rest for WB schemes on a general irregular topography, since oscillation appeared at the shore. Modifying the dierence in the bathymetry level at the interface to recover the equilibrium condition is a possible solution.

Another progress was the use of unstructured grids, which appears useful also in hydrodynamics. Their inherent advantages respect to Cartesian grids consist in the possibility to construct body tted mesh and the exibility when exploiting local mesh renement. Unstructured grids are characterized by the absence of preferential direction. We recall that the Godunov scheme, in two dimensions, is based on a direct extension of the one dimensional scheme.

However, while in one dimension the cell normals are always aligned to the wind direction, in two dimensions and unstructured grids this could be no longer true and upwinding is created also in the crosswind direction to the ow. For this reason, there was much interest to construct multidimensional upwind scheme with optimal diusion property, see [START_REF] Roe | Optimum positive linear schemes for advection in 2 and 3 dimensions[END_REF].

The idea was to construct schemes directly in two dimensions, embedding the multidimensional character of the physics into the numerical methods. This was possible after that [START_REF] Roe | Fluctuations and signals -a framework for numerical evolution problems[END_REF] had presented an alternative interpretation of the one dimensional FV for linear advection. This work shed a new light in the comprehension of upwinding and it was exploited to construct a class of genuinly multidimensional scheme that are nowdays referred to as Residual Distribution (RD). After gasdynamics Deconinck et al. [2000], RD were successfully extended to the SWEs by Paillere et al. [1998] and tested against standard runup benchmarks in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF].

Very high order of accuracy is the next challenge for the hydrodynamic codes. Although in this thesis we will stick to second order of accuracy, which is standard nowdays for commercial codes, we mention the increasing attention that the Discontinuous Galerkin and WENO methods have received in the last years. A review of very high order schemes obtained with these two methods can be found in [START_REF] Xing | A survey of high order schemes for the shallow water equations[END_REF].

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Moving Mesh Methods

h-renement is the most widely used technique to rene locally the mesh size.

The idea is very simple: based on some a posteriori error estimates, the mesh is locally rened or coarsened adding or deleting points from the mesh. Originally this technique was developed for structured Cartesian grids by [START_REF] Berger | Local adaptive mesh renement for shock hydrodynamics[END_REF] under the name of Adaptive Mesh Renement (AMR). Nowdays AMR tools are implemented in many geophysical codes such as GeoClaw Berger et al. [2011] and Gerris Flow [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF]. They have been already tested to improve 2011 Tohoku tsunami simulation accuracy in MacInnes et al. [2013] and [START_REF] Popinet | Adaptive modelling of long-distance wave propagation and ne-scale ooding during the tohoku tsunami[END_REF].

For unstructured grids, h-renement (node insertion/deletion) can be performed by appropriate local remeshing strategies, with impressive results shown in [START_REF] Alauzet | 3d transient xed point mesh adaptation for time dependent problems. application to cfd simulations[END_REF]; [START_REF] Isola | Arbitrary lagrangian eulerian formulation for two-dimensional ows using dynamic meshes with edge swapping[END_REF]. In this thesis we have focused on a more recent strategy of mesh adaptation for unstructured grids based on nodes redistribution (or relocation). These are known as r-adaptation techniques. We invite the interested reader to see the recent review of Budd et al. [2009a]. Roughly speaking, these methods, always based on a posteriori error estimates, cluster the points of a given reference mesh, keeping the mesh topology and number of mesh points unchanged. Now days they are widespread in the CFD community including aerodynamics [START_REF] Bibliography Li | Moving mesh methods with upwinding schemes for time dependent pdes[END_REF], magnetohydrodynamics [START_REF] Tan | Adaptive moving mesh methods for two-dimensional resistive magneto-hydrodynamic pde models[END_REF], multiphase ow Chen et al. [2008], meteorology Budd et al. [2009b] and recently they have been applied also to hydrodynamics Zhou et al. [2013b]. The central idea is the one dimensional equidistribution principle of de Boor [1973]. We search for a time dependent transformation x = M (X, t) from a computational domain described by coordinate X to a physical domain described by coordinate x. The map M should be computed such that it equidistributes a monitor function (error estimates) m on the reference domain, that is m∆x = ∆X. Given a uniform computational grid ∆X = const, if m is large (large error) the grid space ∆x is automatically rened.

The main reasons behind the choice of r-adaptation can be summarized:

• Although very powerful, remeshing techniques requires much higher overheads and a very complex data structures. This results in a higher computational cost, especially compared to a single step of an explicit discretization of the SWEs (as the one we will use in this thesis). On the contrary, in r-renement, nodal movement is obtained by solving an appropriate Moving Mesh Partial Dierential Equation (MMPDE) whose computational cost is of the same order of the CFD solver.

• the potential shown in the past for these techniques to capture shocks, boundary layers and singularities. The MMPDE of [START_REF] Ceniceros | An ecient dynamically adaptive mesh for potentially singular solutions[END_REF], which is the one implemented in this thesis, merged this advan-Introduction tages with a very ecient implementation consisting in the resolution of a sparse non-linear system.

• ecient conservative/accurate remaps of ow variables can be carried out from one grid to the updated one [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF] We have to mention that r-adaptation is less developed and less understood with respect to h-adaptation; The fact that mesh topology is xed poses several limitations: as stated in Budd et al. [2009a] it may never be possible to resolve all of the ne structures of a PDE as it evolves. Moreover great care has to be put in the denition of the monitor function and in the choice of the number of mesh points, since there is the serious risk that the nodes migrate to regions with high gradient, leaving an extreme low number of mesh points to approximate smooth structure of the ow. Mesh tangling is another issue that we will discuss and that deserves further studies.

In this thesis we will take an existing MMPDE (we said, the one of [START_REF] Ceniceros | An ecient dynamically adaptive mesh for potentially singular solutions[END_REF]) and we discuss extensively the coupling of the SWEs (or a general ow solver) with this MMPDE. In particular the SWEs and the MM-PDE can be either solved simultaneously or alternately. The latter has been successfully implemented by [START_REF] Huang | Moving mesh strategy based on a gradient ow equation for two-dimensional problems[END_REF]. Historically, depending on the framework in which we evolve the SWEs, two dierent alternate algorithms were developed and they will be tested in this thesis.

The rst is the rezoning method introduced in Tang and Tang [2003].

The SWEs are written in Eulerian framework and they are are solved on the last updated xed mesh at each time iteration. The advantage is that the CFD solver can be treated as a black box; the drawback is that it requires a remap/interpolation of the ow variables on the updated mesh from the previous mesh. This operation may be quite expensive as it needs to guarantee the same properties as the ow solver itself (high order accuracy, non-oscillatory character/positivity preservation, Well Balancedness, mass conservation).

At the opposite, one can evolve the ow with an Arbitrary-Lagrangian-Eulerian (ALE) formulation of the governing equations, as suggested e.g. in [START_REF] Huang | Moving mesh strategy based on a gradient ow equation for two-dimensional problems[END_REF]; [START_REF] Ni | Remapping-free ale-type kinetic method for ow computations[END_REF]. SWEs are rewritten in an arbitrary moving reference framework which follows mesh movement, see Donea [1983]. The resulting scheme is interpolation free and the properties of the ow solutions are only determined by the scheme. However, a proper ALE form of the numerical discretization has to be used. In particular, a well known requirement for ALE discretizations is the compatibility with a Geometric Conservation Law (GCL), which guarantees that no articial volume is produced in the computational domain due to mesh motion. The discrete counterpart of this property is known as the DGCL (cf. [START_REF] Thomas | Geometric conservation law and its application to ow computations on moving grids[END_REF], [START_REF] Etienne | Geometric conservation law and nite element methods for 3d unsteady simulations of incompressible ow[END_REF] for an overview). Ideally, in Shallow Water ows, we have to ensure the satisfaction of both the DGCL, and of the WB, while still Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. being able to conserve mass and momentum. Up to the our knowledge, this issue has not been addressed yet and it widely discussed in this manuscript.

Tsunami simulation and earth curvature eects

There is a last key ingredient for performing large scale tsunami simulations: earth curvature eects. In this thesis we have chosen to set properly the SWEs in standard latitude-longitude curvilinear coordinates. This is the standard model for numerical weather prediction and ocean modeling. Within this community, Finite Volume on the sphere has already received considerable attention, in the context of both Semi-Lagrangian methods and Eulerian methods, see the review of [START_REF] Machenhauer | Finite-volume methods in meteorology[END_REF]. In this thesis we will refer to the latter class of schemes: the Eulerian FV based on the ux form of the SWEs implemented in [START_REF] Ullrich | High-order nite-volume methods for the shallow water equations on the sphere[END_REF]; [START_REF] Pogorelov | A wave propagation algorithm for hyperbolic systems on curved manifold[END_REF].

Recently Residual Distribution have been extended to curvilinear coordinates only for steady problems [START_REF] Rossmanith | Residual distribution schemes for hyperbolic balance laws in generalized coodinates[END_REF]. They appear to us very interesting due to their simplicity to include the spherical geometry starting from a code in Cartesian coordinate. We stress that our objective is the inclusion of earth curvature in the tsunami dynamics. For this reason we don't need a global scale modeling which means that the complication of the poles, where the coordinate transformation become singulars, will be neglected. ALE-SWEs are not much discussed in the literature. Although ALE balance law on manifold were already presented in [START_REF] Savidis | Arbitrary lagrangianeulerian nite element formulation for geotechnical construction processes[END_REF], to the best of the our knowledge this thesis represent the rst attempt to set the SWEs in ALE form on the sphere.

In geophysics mesh adaptation is gaining increasing popularity. Atmospheric motion is a multiscale phenomenon and it is very hard, even with modern supercomputers to capture all the scale of interest in a single computation. Since the works of [START_REF] Skamarock | An adaptive multiblock highorder nite-volume method for solving the shallow-water equations on the sphere[END_REF]; [START_REF] Behrens | An adaptive semi-lagrangian advection scheme and its parallelization[END_REF], hadaptation on the sphere has strongly developed through AMR techniques, see for example the recent McCorquodale et al. [2015]. Recently also r-adaptation was reconsidered in order to increase the resolution of PDEs on the sphere. It was recognized that the problem is much harder with respect to r-adaptation in Cartesian coordinates. The main diculty is related to nd a unique map from the computational sphere to the physical sphere [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 mongeampère type equation[END_REF]. We will bypass the problem taking advantage of the fact that our domain is only a part of the whole sphere. We will discuss a naif approach consisting in mapping a portion of the sphere into a plane.

Introduction

Main results of the thesis

We highlight directly the main contribution of this thesis. In the next section we will detail each point with respect the structure of the manuscript.

1. We have developed an ALE Residual Distribution and Finite Volume scheme for the SWEs on moving mesh. These schemes are Well-Balanced and mass conserving.

2. The ALE approach has been tested (on standard long wave benchmarks) against the popular rezoning method, showing a slight increase of the overall performances in term of accuracy and CPU time. We have proposed an improved version of the rezoning algorithm which uses a simple linear advection scheme to interpolate the error estimate.

3. We have extended the ALE Residual Distribution and ALE Finite Volume for the SWEs to curvilinear coordinates. 4. We have run a moving mesh simulation of the 2011 Tohoku-Honsu tsunami.

In the following we provide also a list of the publications related to this manuscript:

1. L. A., M. Ricchiuto. ALE r-adaptive methods for the Shallow Water equations in curvilinear coordinates (in preparation).

2. L. A. and M. Ricchiuto. r-adaptation for Shallow Water ows: conservation, well balancedness, eciency. (accepted at Computers and Fluids).

3. L. A., M. Ricchiuto and R. Abgrall. An ALE formulation for explicit Runge-Kutta Residual Distribution. Journal of Scientic Computing, 190(34):14671482, 2014.

For completeness we mention another project in which the author was involved during the three years of PhD. Always in the context of the numerical simulation of long waves, we have investigated the large scale mechanism of tidal bore formation in convergent alluvial estuaries. The work was inspired by the experimental campaign in the Garonne river of Bonneton et al. [2015] in which it was shown that the dominant dimensionless parameter for bore formation appeared to be the dissipation one, which multiplies the frictional term of the SWEs. Although the basic conditions are well known (large tidal range, shallow and convergent channel, low freshwater discharge), a parametric estuarine classication in terms of bore occurrence did not exist in the literature. It was also generally accepted that tidal bores form in estuaries which amplify the incoming tidal wave. We have numerically investigated the estuarine parameter space and proposed a new scaling for the SWEs equations that ensures Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. a clear separation of the eect of non-linearity and friction dissipation. As a consequence, we ended up with a set of dimensionless parameters dening a space in which real alluvial estuaries developing bores are clearly divided from those where bores are not observed. The formation process of an undular tidal bore in such idealized convergent estuaries was reproduced by the use of Boussinesq-type equations and dispersive eects on the local scale of the bore were quantied. For the interested reader, the published papers related to this works are:

1. A.G. Filippini, L. A., P. Bonneton and M. Ricchiuto. Modeling analysis of tidal bore formation in convergent estuaries. (submitted to European Journal of Mechanics / B Fluids).

2. P. Bonneton, A.G. Filippini, L. A., N. Bonneton and M. Ricchiuto. Conditions for tidal bore formation in convergent alluvial estuaries. Estuarine, Coastal and Shelf Science. 172, 121-127, 2016.

Structure of the manuscript

Chapter 1 is devoted to present two existing explicit second order schemes for the Eulerian Shallow Water equations (SWEs) on xed unstructured grids.

They are a Residual Distribution scheme and a more classical Finite Volume scheme which is used for comparison.

In chapter 2 we propose an Arbitrary Lagrangian Eulerian (ALE) form of the SWEs allowing moving mesh simulations with wetting/drying fronts. In order to do this, we systematically review the forms of the SWEs which are best suited for preserving fundamental physical solutions on moving meshes such as the lake at rest or a uniform ow at constant depth. We use the resulting model equations to provide second order ALE Finite Volume and Residual Distribution discretizations. In particular we provide a simple recipe to marry the preservation of the lake at rest and mass conservation on moving meshes using a re-interpolation of the nodal bathymetry based on accurate quadrature of the given bathymetric data.

In chapter 3 the ALE moving mesh method is tested against the popular rezoning approach in which the Eulerian SWEs are solved on the last adapted mesh and an appropriate interpolation permits to transfer the solution from the previous mesh. We study the impact of cheaper and simpler interpolation algorithms which retains all the desired discrete properties. A novel simple ad-hoc error estimators allows to track the shorelines. Finally the ALE and the rezoning moving mesh algorithms are evaluated in terms of CPU time for a given resolution, using the standard benchmarks for near shore hydrodynamics of the WP1.

Introduction

In chapter 4 we will address the numerical approximation of the ALE-SWEs in curvilinear coordinates. The validation of the partners' codes on the sphere was not addressed in the TANDEM project. In particular the benchmarking work-package WP1 did not foresee any test-case to validate the accuracy and the Well-Balanced properties on the sphere. Firstly, we have chosen to set properly the Eulerian SWEs in standard latitude-longitude curvilinear coordinates. We will take cases from the standard test-suite of [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF] to validate our Finite Volume and Residual Distribution scheme for the SWEs in curvilinear coordinates on unstructured grids. Secondly we have set the ALE-SWEs on the sphere. A novel Finite Volume and Residual Distribution approximation of the resulting equations is proposed. Third, a simple moving mesh method will be established which allows point movements on the sphere. Numerical tests show that our resulting moving mesh algorithm, can improve the resolution of linear and nonlinear waves on the sphere, using a limited number of mesh points.

Finally in chapter 5 we address the issue of WP3. We present xed and moving mesh SWEs simulations of the 2011 Tohoku-Honsu event. The numerical tools developed along this manuscript will be compared with other partners involved in the WP3. This should provide evidence that the moving mesh method studied in the thesis, although simple, could be a good candidate to reduce the large computational cost of tsunami simulations.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Chapter 1 Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids

Shallow Water equations

The use of Shallow Water equations (SWEs) is ubiquitous in hydrodynamics and geophysics. SWEs is a nonlinear hyperbolic system of partial dierential equations that describes wave propagation in a uid that is shallow relative to the wavelength. SWEs can be derived from the conservation of mass and momentum of a vertical uid column under a crucial assumption, the so called hydrostatic approximation. We consider the vertical uid velocity negligible and the horizontal velocity about constant along the the vertical coordinate from the bottom to the free surface, which is true if the uid depth is shallow compared to the horizontal length scale. Moreover since the uid is at rest along the vertical coordinate the pressure is determined just by the weight of uid above, that is p = p(h), h here is the water depth. There are many circumstances under which these hypothesis lead to a good approximation of natural phenomena such as waves in the atmosphere, tides in ocean but also breaking waves on a beach and, more important for us, tsunamis. If we think of tsunamis generated by an oshore quake, they are very long waves with wavelength L in the order of ∼ 100 km, much larger then the water depth D ∼ 1 km. See e.g. [START_REF] Benoit | Long-distance propagation of tsunamis and surface waves: on the relative importance of dispersion and nonlinearity[END_REF]; [START_REF] Lannes | Pde modelling (parameters, notion of asymptotics)[END_REF] for the scaling analysis.

A set of Cartesian coordinates x = {x 1 , x 2 } ∈ Ω dened on the orthonormal basis {e 1 , e 2 } is given at this point. The vector basis dene a plane which we assume that the water free surface in undisturbed condition lies on. Then, the

Shallow Water equations

SWEs in engineering notation writes:

∂h ∂t + ∂hu j ∂x j = 0 ∂hu i ∂t + ∂T ij ∂x j + S i = 0
the repeated index j = 1, 2 is implied in summations and the index i = 1, 2

is the direction of the momentum equation. The momentum discharge vector is dened as hu = hu j e j and the momentum ux tensor is T = T ij e i e j with components T ij = huu ij + 1 2 gh 2 δ ij . The source term includes the bottom topography force and the frictional one

S i = h g ∂b ∂x i + c F u i b(x)
is the bottom topography, measured form the origin which is the undisturbed free surface level. c F is the friction coecient dened through Manning's empirical formula c F = gn 2 u h 4/3 and n is the Manning number. We also provide the vector form of the SWEs, which has the advantage of being more compact

∂u ∂t + ∂F j ∂x j + S(x, u) = 0, (1.1) u = h hu i , F j = hu j T ij , S = 0 gh ∂b ∂x i S b + 0 c F hu i S f (1.2)
We collect the ux components in the ux matrix

F = [ F 1 F 2 ] ∈ R 3 × R 2 .
Equations (1.1)(1.2) represents a non-homogeneous non-linear hyperbolic system of partial dierential equations. Given any vector ξ = ξ i e i , the ux Jacobian K(u, ξ) = ∂F j ∂u ξ j admits a full set of real eigenvalues and linearly independent eigenvectors, namely

K =   0 ξ 1 ξ 2 c 2 ξ 1 -uu • ξ u • ξ + uξ 1 uξ 2 c 2 ξ 2 -vu • ξ vξ 1 u • ξ + vξ 2   (1.3)
with eigenvalues λ(u, ξ) to the water depth, the speed of every point may dier due to the dierent height, leading to a distortion where the crest is moving faster respect to the trough. Two typical hyperbolic phenomena occurs: the wave front steepens (hydraulic jumps) and the tail attens (depression). This can be clearly seen for waves approaching the beach, at least before that the familiar breaking of the wave occurs.

λ 1,3 = u j ξ j ± c ξ , λ 2 = u j ξ j
Later in the text we will also make use of the Jacobian at rest A = K(u = 0)

A(h, ξ) =   0 ξ 1 ξ 2 c 2 ξ 1 0 0 c 2 ξ 2 0 0   (1.5)
that admits the non-linear eigenvalues λ = ±c ξ that characterize the acoustics of the system.

Conservation form

The SWEs (1.1)(1.2) constitute a balance law and allow to recover the integral balance of of fundamental physical quantities such as mass and momentum.

For a xed volume C with boundary ∂C and boundary normal n = n j e j we can state integral conservation for the vector of conservative variables u:

C u(x, t) dx = C u(x, 0) dx - t 0 ∂C F j n j ds + C S dx dt = 0 (1.6)
These are the most fundamental equations when the solution contains discontinuities. In this case, the dierential form does not hold anymore but

(1.6) still admits weak solutions corresponding to the Rankine-Hugoniot jump conditions, see [START_REF] Bernetti | Exact solution of the riemann problem for the shallow water equations with discontinuous bottom geometry[END_REF]:

s[h] = [hu] s[hu] = [hu] + 1 2 gh 2 - 1 2 h 2 L -(h L -(b L -b R ) 2 )
where we dened the jump of any scalar φ across the discontinuity as [φ] = φ R -φ L and the shock velocity as s.

Shallow Water entropy

Being conservative is a delicate issue for every numerical method approximating conservation laws. It means that the numerical solution satises the integral form of conservation laws, mimicking what the exact solution does.

For conservative methods, one can hope, as the mesh is rened, that the numerical solution converges to the weak solution of (1.6) and, even in presence of a discontinuity, we compute the correct wave speed. This was proved by Lax and Wendro [1960].

Shallow Water entropy

An hyperbolic equation should be supplied by some conditions, which ensure that the weak solution of the problem is physically correct. These conditions are called entropy conditions because in gasdynamics the correct weak solution is selected by the second principle of thermodynamics which states that the entropy of a system must be non increasing with time. For SWEs, the role of the entropy is played by the depth-integrated total energy per unit mass, see

E(u) = h gh 2 + gb + u 2 2
A conservation law for the SW energy can be obtained projecting mass and momentum equations (1.1), on the so called entropy vector

Ψ = ∂E ∂u = [ (h + b)g - u 2 2 u v ] T
We obtain the conservation of energy:

∂E ∂t + ∂ψ j ∂x j = 0 (1.7)
where the entropy ux is ψ = ψ j e j with components ψ j = u j E + 1 2 u j gh 2 .

However if we go back to the complete set of physical equation with viscosity, then we can argue that the admissible solution should be the one in the vanishing viscosity limit. [START_REF] Lax | Hyperbolic system of conservation laws and the mathematical theory of shock waves[END_REF] demonstrated that the vanishing viscosity limit for a system of conservation law with positive convex entropy ∂ 2 uu E > 0, is equivalent to an entropy inequality. For the SWEs this is

∂E ∂t + ∂ψ j ∂x j ≤ 0 (1.8)
(1.8) states that the energy, in absence of boundary conditions, cannot grow.

And this is a physical criteria that prevent from picking up unphysical weak solution, such as expansion shocks. Moreover (1.8) states that the energy is 1. Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids conserved whenever the solution is smooth, but it is dissipated across discontinuities. We will see that this fact is used in numerical methods to detect shock waves in the domain. Finally (1.8) can be seen as stability condition for the SWEs. Given a certain norm of the energy E V (Ω) , entropy inequality corresponds to:

E(t) V ≤ E 0 V
where E 0 V is the norm of the energy at the initial time.

Maximum principle, entropy stability and water depth positivity

Consider for a while, the following scalar hyperbolic equation that conserves the variable u(x, t), for a general non-linear ux function F (u):

∂u ∂t + ∂F j ∂x j = 0, x ∈ Ω, t ∈ [0, T ] (1.9)
L ∞ stability or maximum principle holds for the exact solution:

min x∈Ω u(x, 0) < u(x, t) < max x∈Ω u(x, 0), ∀t ∈ [0, T ] (1.10)
The above condition is very important when one searches for approximate solution u h ≈ u. In particular, maximum principle is used to have precise bounds on the numerical solution which in turn avoids the computations of oscillatory discontinuous solution, see the work of [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF] on maximumprinciple-satisfying high order schemes. This stability property cannot be automatically extended to the non-linear system case since the same maximum principle is not trivial for the general problem described by (1.1). However

(1.10) remains important. Historically, numerical methods have been studied

and tested in the scalar case, where a non-linear stability theory is available.

Only later, scalar methods were extended to the system case, more or less straightforwardly. Even if we do not have any guarantees that the stability property will work well as in the scalar case, this approach was successful.

To cope with the obstacle of L ∞ stability many researchers investigated other types of stability. An important class of methods are, for example, entropy stable schemes [START_REF] Osher | Riemann solvers, the entropy condition and dierence approximations[END_REF], [START_REF] Tadmor | The numerical viscosity of entropy stable schemes for systems of conservation laws[END_REF]), which aim to mimic the physical constraint that the energy cannot increase. By adding a certain amount of numerical viscosity, entropy stable schemes can verify a discrete inequality analogue to the entropy one (1.8).

For the SWEs, another stability condition is a physical constraint: the positivity of water depth h ≥ 0. Preserving the positivity of water depth is an important stability condition for numerical methods since the uncontrolled occurrence of negative depth leads soon to the code breakdown. We will discuss positivity preserving schemes [START_REF] Bouchut | Nonlinear stability of nite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]) later in this chapter.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

1.5. Lake at rest 1.5 Lake at rest

In the context of Shallow Water ows, an important role is played by the socalled "lake at rest" state which, denoting the free surface level η = h + b, is the particular steady solution characterized by the two invariants :

hu = 0, h + b = η 0 = const (1.11)
A numerical method approximating (1.1)(1.2) is said to enjoy the C-property or also to be is also an exact steady solution of the discrete equations. In other words, Well-balanced schemes provide a discrete analog of the balance between hydrostatic pressure and the pressure forces exerted by the bottom ∂F j ∂x j + S = 0 allowing to preserve (1.11) exactly at the discrete level, [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF].

1.6 Discrete approximation 1.6.1 Notation for mesh, geometry and unknowns Consider an unstructured discretization of the spatial domain Ω composed by non overlapping triangular elements. We will denote the grid (or mesh) by T h , h K being the local reference element length (mesh size, hereinafter). K is the generic triangle, |K| its area and n j = n m j e m is the normal to the j-th face of the triangle scaled by the corresponding edge length (see g.1.1). For every node i of the triangulation, D i denotes the subset of triangles containing i. With a little abuse in the notation j ∈ D i is the set of nodes j sharing an edge with node i. We then denote by C i the median dual cell obtained by joining the gravity centers of the triangles in D i with the midpoints of the edges meeting in i, its area is

|C i | = K∈D i |K| 3
In a Finite Volume context we dene also the boundary of the median dual cell as the interface ∂C i = j∈D i ∂C ij . The interface belonging to nodes i, j, denoted as ∂C ij , is the union of two segments connecting the baricenters of the neighboring triangles K i, j with the midpoint of the edge ij (cf. right picture in g.1.1). C ij is the area delimited by ∂C ij and by the two segments joining i with the gravity centers of the elements K i, j. Both C ij and ∂C ij 1. Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids can be split over the neighboring triangles. We dene the normal and the area associated to the interface ij

n ij = 1 2 K i,j n K ij , |C ij | = K i,j |C K ij | with |C K ij | = |K|

6

. We highlight the following geometrical relationship between nodal normals and interface normals: With the subscript i we distinguish discrete variables from continuous ones.

1 2 j∈K,j =i n K j = - n K i 2 = j∈K,j =i n K ij (1.12)
It is a widespread convention in literature, however we note that engineering notation with sub/superscripts to denote vector and tensor components, could confuse the reader. In the following, we will specify whenever standard summation does not apply. For the Finite Volume method we will evolve in time approximations of cell averages of conservative variables over the standard median dual cells and we denote the average as u i . On the contrary, the Residual Distribution method evolves values of the unknowns at mesh nodes. For simplicity, we shall still denote these values as u i .

Finite Volume for SWEs

Finite Volume (FV) is the most popular class of schemes for the discretization of hyperbolic PDEs in computational hydrodynamics. An exhaustive monograph on FV for hyperbolic problem is [START_REF] Leveque | Finite Volume Methods for Hyperbolic problems[END_REF] and a review of Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Discrete approximation

FV methods for the SWEs can be found in [START_REF] Toro | Godunov-type methods for free surface shallow water ows: a review[END_REF].

Here we present the standard well balanced node-centered unstructured scheme based on Roe's linearized Riemann solver as it has been developed by [START_REF] Nikolos | An unstructured node-centered nite volume scheme for shallow water ows with wet/dry fronts over complex topography[END_REF]. This scheme proved to successfully simulate complex problems with irregular topography and inundations. The FV idea is to break the domain into many volumes C i . The solution is then averaged within every cell with

u i (t) = 1 |C i | C i u(x, t) dx (1.13)
and it is updated by a discrete integral balance, equation (1.6)

u n+1 i = u n i - ∆t |C i | j∈D i R ij (u n i , u n j )
where we have

R ij = F ij + S ij (1.14)
with F ij a numerical approximation of the ux along ∂C ij and S ij an approximation of the integral of the source term on C ij .

F ij = 1 ∆t t n+1 t n ∂C ij F j n j ds dt, S ij = 1 ∆t t n+1 t n C ij

Sdx dt

After the averaging step, a Riemann problem is dened at each interface.

Numerical uxes can be computed if the complete solution of the Riemann problem is achieved, typically through the resolution of a non-linear system.

Note, however, that the exact solution of the Riemann solution seems to be redundant since the same solution is roughly averaged just at the end of the time step. For this reason one can think to replace the exact Riemann Problem with an approximate one. The seminal idea consists in linearizing the ux

∂F j ∂x k = ∂F j ∂u ∂u ∂x k .
In this work we use the original Roe linearization Roe [1982].

Recalling expression (1.3), K ij = K(u ij , n ij ) will be the ux Jacobian evaluated at the interface Roe state u ij . Within this hypothesis, the Riemann solution consists in the superposition of simple or p-waves. Numerical ux evolves according to the sum of p-waves entering the median dual cell through the interface, each one carrying a jump in the solution equal to W p :

F ij (u i , u j ) = F i • n ij + m p=1,αp<0 W p = F i • n ij + K - ij (u j -u i ) (1.15) = F j • n ij - m p=1,αp>0 W p = F j • n ij -K + ij (u j -u i ) (1.16)
Please note that the above scalar products are standard ones and the summation convention does not apply in the above formula and in the ones which 1. Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids follow. Matrices are computed through eigenvalues decomposition, e.g.

K + ij = R ij Λ + ij L ij .
To obtain expressions (1.15) and (1.16) we have decomposed the total jump in the solution as follows:

W p = α p λ p r p , α p = (u j -u i ) l p r p , l p are respectively the p-th right and left eingenvectors and are referred, as the eigenvalues λ p , to the Roe state at the interface, e.g. λ p = λ p (u ij , n ij ).

Averaging expression (1.15) and (1.16) we have the numerical ux:

F ij = F ij (u i , u j ) = F(u j ) + F(u i ) 2 • n ij - K ij 2 (u j -u i ) (1.17)
this scheme is known as Roe scheme. Second order of accuracy is achieved with a piecewise linear reconstruction of the vector of conservative variable over the dual cell. In practice we have replaced, in expression (1.17), nodal values u i and u j with linearly reconstructed ones at the cell interface, respectively ȗi and ȗj .

We have combined this with a Green-Gauss reconstruction [START_REF] Nikolos | An unstructured node-centered nite volume scheme for shallow water ows with wet/dry fronts over complex topography[END_REF]; [START_REF] Delis | Performance and comparison of cell-centered and node-centered unstructured nite volume discretizations for shallow water free surface ows[END_REF]. We will refer to this linear second order scheme as to the FROMM scheme, [START_REF] Leveque | Finite Volume Methods for Hyperbolic problems[END_REF]. If necessary in the reconstruction step, a slope limiter is used to enforce monotonicity

ȗi = u i + 1 2 l(u i , u j , ∇u i )
with l, the Van Albada limiter [START_REF] Van-Albada | A compartive study of computational methods for cosmeic gas dynamics[END_REF]. The resulting scheme is a two dimensional implementation of the MUSCL scheme of [START_REF] Van-Leer | Towards the ultimate conservative dierence scheme v. a second-order sequel to godunov's method[END_REF].

To march in time with second order of accuracy we use an explicit Runge-Kutta two (eRK2); the FV discrete evolution equations reads:

u * i = u n i - ∆t |C i | j∈D i R ij (u n ) (1.18) u n+1 i = u n i - ∆t |C i | j∈D i R ij (u n ) 2 + R ij (u * ) 2 (1.19)

Conservation

From the solution decomposition in simple waves (1.15) and (1.16), it is clear that Roe linearization ensures conservativeness. By denition the jump in the ux is:

(F i -F j ) • n ij = K ij (u ij , n ij ) (u j -u i )
and thus (1.15) and (1.16) coincide with the consequent exact conservation at the interface. An advantage of the numerical ux (1.17) is that conservation is true for any linearization, for example the simple arithmetic average u ij = 1 2 (ȗ i + ȗj ). This fact will guarantees the method is conservative even in more complex cases for which the Roe average is not directly available, as we will in chapter 4.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Discrete approximation

Positivity of water depth [START_REF] Perthame | On positivity preserving nite volume schemes for euler equations[END_REF] have shown that Godunov method applied to the Euler equations of gasdynamics, (with exact or approximate Riemann solver)

preserves density and pressure positivity, with a rigorous CFL condition. They also presented second order pressure/density positivity preserving schemes but some other requirement on the reconstruction should be assured. In this work to avoid the appearance of negative depths, following [START_REF] Nikolos | An unstructured node-centered nite volume scheme for shallow water ows with wet/dry fronts over complex topography[END_REF] we propose to switch to a rst order scheme when approaching regions dry regions with an appropriate CFL condition, see (1.38). The limiter is modied according to the denition given in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]:

l * = l(u i , u j )e - h K L ref h ref h min (1.20)
the exponential factor takes into account the occurrence of dry areas. L ref is the domain reference length; h ref is the maximum depth at t = 0; h min is the minimum depth among h i and h j . It is always of order one but it quickly tends toward zero where h is small.

Source terms

Non dierential source terms in (1.1) includes only friction. For second order of accuracy, we can use pointwise values S i = S(x i , u i ) and approximate the numerical source as

S ij = C ij S(x, u) dx = S i |C ij |, j∈D i S ij = S i |C i |
which corresponds to linearly reconstruct also the source term. More complex treatments can be constructed, especially when increasing the accuracy to more then second order.

Well-Balancedness

It is widely recognized that the numerical treatment of the bathymetric source term is as important as that of advective terms. In fact, this particular term contained in the numerical source S ij , in addition to not spoil the accuracy of the advective part, it should preserves of the Well-Balanced property discussed in section 1.5. The rst attempt to a well-balanced source term was due to Bermúdez and Vázquez [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF] who tried to incorporate the bathymetry source term into Roe method. A second important cornerstone was the extension from rst to second order accuracy proposed by [START_REF] Hubbard | Flux dierence splitting and the balancing of source terms and ux gradients[END_REF]. We present here the main conclusions of this last reference. They distinguished 1. Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids two contributions, the rst balancing the central part of the uxes, and the second the upwind dissipation term:

S b ij = S c ij + S * ij (1.21)
We dene the following average values

h - ij = hi + h i 2 , h ij = hj + hi 2
and the bathymetry variation vectors

∆b - ij =   bi -b i 0 0   , ∆b ij =   bj -bi 0 0   (1.22)
moreover we recover the denition (1.5) (Jacobian at rest) to introduce

A - ij = A(h - ij , n ij ), A ij = A(h ij , n ij ).
S c ij = A - ij ∆b - ij + 1 2 A ij ∆b ij (1.23)
Concerning the upwind balancing term, the original denition given in [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF]; [START_REF] Hubbard | Flux dierence splitting and the balancing of source terms and ux gradients[END_REF] leads to the following expression

S * ij = - sign(K ij ) 2 A ij ∆b ij (1.24)
The proof of Well-Balancedness rests on the property of the Roe average and the fact that, on the lake at rest state, we have

K ij = A ij . On each edge R ij = F(ȗ j ) -F(ȗ i ) 2 • n ij + (F(ȗ i ) -F(u i )) • n ij - |A ij | 2 (ȗ j -ȗi ) + 1 2 A ij ∆b ij + A - ij ∆b - ij - sign(A ij ) 2 A ij ∆b ij = 0
Using the fact that, on the selected equilibrium, (F(ȗ j ) - ȗi ) we can sum each term of the numerical ux with the corresponding bathymetric term. Note now that ȗjȗi + ∆b ij = 0 by hypothesis.

F(ȗ i ))•n ij = A ij (ȗ j -

An introduction to Residual Distribution

At the beginning of the eighties, Roe's research on the multidimensional extension of Godunov method brought him to formalize FV into a form called Fluctuation Splitting [START_REF] Roe | Fluctuations and signals -a framework for numerical evolution problems[END_REF][START_REF] Roe | Linear advection schemes on triangular meshes[END_REF]. This has paved the way to the development of multidimensional upwind schemes and to what today is referred

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

1.6. Discrete approximation as Residual Distribution (RD). The main ideas behind RD are recalled here, with the main purpose of introducing denitions that will be used extensively in the manuscript. The interested reader can nd the rigorous foundations of the method in [START_REF] Deconinck | Residual distribution schemes: Foundations and analysis[END_REF] and a recent review in [START_REF] Abgrall | High order methods for cfd[END_REF]. We start from the FV framework, which the reader is more condent with, and, as Roe actually did, we develop the concept of uctuation. We write a general ux (e.g. Roe uxes with MUSCL reconstruction

(1.17)) compactly as

F ij = H(ȗ i , ȗj )
and we decompose them in the contribution of each element, the FV update writes:

u n+1 i = u n i - ∆t |C i | j∈D i H(ȗ i , ȗj ) = u n i - ∆t |C i | K∈D i j∈K,j =i H(ȗ i , ȗj ) since j∈D i H(u i , u i ) = F i • j∈D i n ij = 0
we can add the following term:

u n+1 i = u n i - ∆t |C i | K∈D i j∈K,j =i (H(ȗ i , ȗj ) -H(u i , u i ))
(1.25)

This step seems arbitrary but it is not. In fact, in each element, we are implicitly making use of a very special ux. Using the property (1.12) allow to show that:

-

j∈K,j =i H(u i , u i ) = - j∈K,j =i F i • n ij = -F i • j∈K,j =i n ij = F i • n i 2 = - j∈K,j =i F i • n j 2
that is, we are taking into account also the numerical uxes along half the edges with vertex i in common. This contributions cancels out between the elements sharing the same edge and do not modify the original scheme. However, summing over the vertexes, the internal numerical ux contribution disappears for the conservation of an element and we are left with

Φ K = i∈K j∈K,j =i (H(ȗ i , ȗj ) -H(u i , u i )) = - i∈K j∈K,j =i H(u i , u i ) = - i∈K j∈K,j =i F i • n j 2 = 1 2 i∈K F i • n i ≈ ∂K F • n ds
1. Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids which is a second order quadrature formula approximating the ux balance over the element K. We call this quantity uctuation of the element K. The FV scheme could be seen as dening a proper splitting or distribution Φ K i of the uctuation, such that:

i∈K Φ K i = Φ K so that we don't loose consistency property F ij (u, u) = F(u) • n.
The split uctuation for FV writes (dropping superscript K in the residual notation, we take it for granted)

Φ F V i = j∈K,j =i (H(ȗ i , ȗj ) -H(u i , u i )) (1.26)
and veries the consistency condition. Substituting the denition (1.26) in

(1.25) we get a reformulation of the standard FV method in Godunov form

u n+1 i = u n i - ∆t |C i | K∈D i Φ F V i (1.27)
We call this compact update as the uctuation form of FV. The solution is evolved at each time step according to a balance of element's uctuations rather then interface's uxes. All this may seem just a trivial exercise but the promising aspect is that an abstract geometrical interpretation of upwinding arises, that is, the directions in which we distribute the uctuation with respect to the direction of simple waves celerities. We can in fact generalize the splitting step and create new schemes with improved properties. For example, this fact is used to introduce upwinding in a more clever way. Indeed early experiments demonstrated that, for multidimensional problems, rst order RD were less diusive then rst order FV [START_REF] Roe | Optimum positive linear schemes for advection in 2 and 3 dimensions[END_REF]; Paillere et al. [1998]. We also remark that FV, as they have presented in section 1.6.2, have been implemented in the RD form (1.27) and not in the classical ux form.

RD prototype scheme

Once we have presented FV in uctuation from (1.27), the last eort consists in generalize the splitting procedure and to consider a splitting operator dierent from (1.26). These lead us directly to the RD approximation of hyperbolic systems such as the SWEs (1.1) and which is constructed with the following three steps: given an approximation of conservative variables at the nodes u i (t) = u(x i , t), we introduce the following continuous numerical approximation

u h (x, t) = i∈T h ϕ i (x)u i (t)
(1.28) {ϕ i } i∈T h will be the standard P 1 continuous piecewise linear Lagrange kernel.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Discrete approximation

• For each element K, compute the ux balance, called uctuation:

Φ K (u h ) = ∂K F j (u h )n j ds + K S(u h , x) dx (1.29)
• Split the uctuation in contributions, one for each node of the element, through an appropriate distribution matrix β K j Φ K j :

j∈K

Φ K j = Φ K , or Φ K j = β K j Φ K (1.30)
for consistency reason it is clear that ⇒ j∈K β K j = I 3 .

• envoy at each node of the element the corresponding split uctuation, updating progressively the solution at vertexes. At the end, the RD method reads

|C i | du i dt + K∈D i Φ K i = 0, ∀i ∈ T h (1.31)
In gure (2.6) there is a sketch which summarize the three successive steps.

The key properties of the method are determined by the denition of the split uctuation, or, if you prefer, by the distribution matrix. In the following paragraph we relate the properties of the distribution matrix to stability and accuracy results. Leaving apart a rigorous treatment of the stability and accuracy analysis, we will try instead to highlight basic concepts. We hope that this could help the reader in a better comprehension of the RD results.

|C i | du i dt + K∈D i Φ K i = 0
We can develop the advective part with the consistency condition contained in (1.30)

i∈T h K∈D i Φ K i = K∈T h j∈K Φ K j = K∈T h ∂K F j n j ds = ∂Ω F j n j ds
we have assumed that numerical uxes are continuous along the triangle edges ∂K such as it happens in the continuous Finite Element method. Integrating in time with

U (t) = i∈T h |C i |u i (t) we get U (t) = U (0) - t 0 ∂Ω F j n j dsdt
which states that we have exact conservation over the full domain if the consistency condition is veried and if numerical uxes cancel out at the element's edges. Note that this is a continuity condition for F j that is satised by several denitions of the discrete ux. We will use reinterpolation of ux's nodal values.

Positivity

Consider again the scalar conservation law (2.65). For linear schemes the scalar uctuation φ K i can be put in a linear compact form

φ K i = j∈K, =i c K ij (u i -u j ) (1.
u n+1 i = 1 - ∆t |C i | j∈K,j =i c K ij u n i + j∈K,j =i ∆t |C i | c K ij u j = j∈K a K ij u j u n+1
i can be written as a convex combination of the solution at the previous time step, enforcing the positivity of the coecients a K ij . These schemes are said to be Positive (P). For scalar LED-RD under a proper CFL condition, L ∞ stability follows

min j∈T h u 0 j < u n i < max j∈T h u 0 j , ∀i ∈ T h , ∀n ∈ [0, M ]
The direct extension of Positive RD to the system case was successful in practice. However we mention that the extension of the stability result could be misleading. First of all, the existence of a maximum principle for a system of non-linear conservation laws is not a trivial task. This does not mean that a non oscillatory behavior near discontinuities is desired/expected. If we extend the positive analysis to systems, it is not easy to handle with the LED condition, this time applied to a matrix C ij < 0. To avoid this diculties, it is necessary to searches for some conditions relying on entropy consideration (we have seen that should always decrease) but this is beyond the scope of the brief paragraph, see [START_REF] Barth | An energy look at the n scheme[END_REF]; [START_REF] Ricchiuto | Contruction and analysis of compact residual discretizations for conservation laws on unstructured meshes[END_REF]. Instead, for the SWEs much attention should be put to preserve the physical constraint of water depth positivity. In section 1.6.5 we will proof that a positive RD scheme applied to

(1.1) preserves this property.

Linearity Preserving

We introduce now Linearity Preserving (LP) schemes which were developed for the accuracy analysis of RD. However the reader more familiar with FV will nd analogies with the k-th exactness schemes of Barth. We consider a steady scalar problem with a smooth exact solution v such that

∂F j (v) ∂x j + S(x, v) = 0
Keeping the notation of (1.28), v h is the continuous piecewise polynomial approximation of v obtained with a general Lagrangian basis {ϕ i } i∈T h of order q.

First we give a trivial denition. For steady scalar problem a numerical scheme is q-th order accurate if the truncation error is of order T E(v h ) = O(h q K ). The local truncation error in RD case writes

T E = i∈T h ϕ i K∈D i φ K i (v h )
1. Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids with φ K i being the split residual in scalar notation. Less trivial is the fact that the truncation error estimate T

E(v h ) = O(h q K ) is veried if the split uctuation is φ K i = O(h q+1 K ) (1.33)
the proof can be nd in [START_REF] Deconinck | Residual distribution schemes: Foundations and analysis[END_REF]. Given continuous and q-th order accurate ux and source terms approximations, [ F 1 h , F 2 h ] and S h , for a smooth exact solution one has also that

φ K = ∂K F j h n j ds + K S h dx = O(h q+1 )
Now we can take advantage of the RD formalism. The fundamental relation

φ K i = β K i φ K (
max K∈T h max j∈K β K j < C ∀φ K , v h , v 0 h
Using (1.33) we deduce that, for q = 2 (piecewise linear Lagrange basis), a scheme which is linearity preserving is second order accurate at steady state.

All these results extend directly to the system case, in particular to the SWEs.

Distributions

Through a proper choice of the distribution matrix/distributed uctuation, many classical schemes can be recovered. The FV scheme with Lax Friedrich (LxF) uxes is important for its ability to compute non oscillatory discontinuous solutions and to have, at least on some variables [START_REF] Perthame | On positivity preserving nite volume schemes for euler equations[END_REF],

precise bounds. The LxF uctuation writes

Φ LxF i = 1 3 Φ K + j∈K, =i α K 3 (u i -u j ) , α K = max j∈K K j
in the scalar case, the diusion parameter α K is chosen to verify the positivity (P) requirement of the coecients. For systems one can take some upper bound within the time step to the the largest absolute value of the ux Jacobians evaluated in the nodes of an element. Unfortunately LxF is only rst order accurate. Using the necessary condition of LP schemes, a root to the construction of non linear schemes which are LP and P consists in limiting the unbounded coecients of a P scheme. We have used the PSI limiter of [START_REF] Deconink | A multidimensional of roe ux dierence splitter for the euler equations[END_REF] Φ LLxF i = m=1,3

β m i l T m Φ LxF i r m , β m i = l T m Φ LxF + l T m Φ LxF j + (1.34)
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Discrete approximation

where the apex LLxF stands for Limited LxF scheme. Formally this scheme is linearity preserving and positive. However, things turned out not so simple.

Discontinuity are well handled while in smooth regions we encounter poor accuracy with wiggles. These oscillations can be seen as the rising of some unexpected destabilizing phenomena. Source of energy instability might be introduced by the PSI limiter. A possible solution to cure the problem is suggested by [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] is to add a SUPG term (which is LP but not P).

The SUPG scheme of Hughes and co-workers (see [START_REF] Hughes | Streamline upwind Petrov-Galerkin formulations for convection dominated ows with particular emphasis on the incompressible Navier-Stokes equations[END_REF];

Hughes et al. [2010] and references therein for details) can be written in terms of uctuations:

Φ SU P G i = β SU P G i Φ K , β SU P G i = 1 3 + K j T (1.35) with T = j∈K |K j | -1
. A limiter tune the streamline upwind diusion introduced. The nal distribution for such a scheme called LLxF stabilized or briey LLxF-SUPG reads

Φ LLxF -SU P G i = (1 -δ(u h )) Φ LLxF i + δ(u h )Φ SU P G i
The limiter is based on the energy considerations of section 1.3. [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] propose the following heuristic denition

δ = min 1, E ref u ref |K| |φ K E |L ref with φ K
E an approximation of the entropy/energy uctuation and is obtained projecting the uctuation on the entropy vector φ K E = Ψ T • Φ K . When the solution is smooth enough, we can scale the energy uctuation as φ

K E ∼ E ref u ref |K|L -1 ref .
On cells where the solution is discontinuous the energy uctuation becomes singular, φ K E → ∞ and δ → 0.

Source terms

The treatment of source term is very simple and eective in the RD framework.

In fact, one could think to include every source term in its wave propagation algorithm, see LeVeque [1998]. While, for FV, this could be laborious, for RD it results automatic. According to the rst step of RD algorithm (1.29), general source terms can be integrated with Gaussian quadrature formula with the prescribed order of accuracy (second for us):

Φ K S = K S(x, u) dx = |K| Nq q=1 ω q S q
Then it is added to the uctuation and split automatically with the distribution matrix of the advective scheme. For explicit schemes however, once should take into account these terms in the stability analysis and extends the CFL condition in presence of the source.

1. Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids 1.6.4 Explicit RD Runge-Kutta two scheme

The accurate simulation of unsteady problems, is crucial in many coastal and geophysical applications. We search for a time discretization that gives second order of accuracy and, in particular that leads to an explicit scheme in order to compare the results with the FV-eRK2 presented in section (1.6.2). For unsteady RD, the time derivative, once added to the uctuation, dene the element's residual

Φ K(t) = K du h dt dx + ∂K F j (u h )n j ds + K S(u h , x) dx
and we denote it with the same symbol of the uctuation, adding the superscript (t) to underline the presence of the temporal derivative. We hope that this does not generate confusion. The time part leads to the appearance of a mass matrix, (see [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for navier stokes simulations[END_REF] for a discussion on mass matrices) that must be inverted locally:

β K i K du h dt dx = j∈K m K ij du j dt , m K ij = K w i ϕ j dx
with w i a Petrov-Galerkin test function chosen such that a coherent RD splitting for the time part is recovered β K i = K w i dx.

In this thesis we have implemented the explicit Runge Kutta 2 or Predictor-Corrector (PC) scheme, described in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF]. Through an ecient mass-lumping strategy, the RD-eRK2 scheme allows to march in time explicitly and without the unnecessary additional cost of inverting the mass matrix. The algorithm follows:

1] Predictor step: for each element K ∈ T h

• Compute the residual Φ K(1) = Φ K (u n h ).
• Distribute the residual to the nodes of K such that j∈K Φ

K(1) j = Φ K(1)
• Compute the rst order prediction of the solution, denoted as u *

u * i = u n i - ∆t |C i | K∈D i Φ K(1) i (u n h ) (1.36) 2] Corrector step: for each element K ∈ T h • Compute the residual Φ K(2) = 1 2 Φ K (u n h ) + Φ K (u * h ) + K u * h -u n h ∆t dx
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

1.6. Discrete approximation

• Distribute the residual to the nodes of K such that j∈K Φ K(2) j = Φ K(2)
• Compute the second order correction from

u n+1 i = u * i - ∆t |C i | K∈D i Φ K(2) i (u * h , u n h )
(1.37)

Residual Distribution for SWEs

Trough this manuscript, we will compare the eRK2-FV method for the SWEs proposed by [START_REF] Nikolos | An unstructured node-centered nite volume scheme for shallow water ows with wet/dry fronts over complex topography[END_REF] to an RD discretization of the same system recently suggested by [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF]. The numerical resolution of the SWEs with RD schemes poses two major issues. The rst is numerical preservation of positive depths. The second is is the accurate discretization of source terms, bathymetry and friction up to now, recalling that we would like to embed at a discrete level the exact balance in the case of lake at rest.

Positivity of water depth

For SWEs simulations, the positivity of water depth is a physical constraint.

We study the LLxF uctuation, with a limitation performed as if we were in the scalar case, which means equation by equation, or R = L = I 3 in (1.34).

This variant is interesting because [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] showed that it preserves the positivity of water depth. We repeat the proof for sake of clarity only for the predictor step. The reader can nd the complete proof in the reference. The distributed uctuation for the mass equation writes

φ LLxF i = φ LLxF i φ K φ K φ LxF i φ LxF i = β LLxF i β LxF i γ i 1 3 φ K + α K 3 j∈K, =i (h i -h j ) = γ i 1 6 j∈K h j u j • n j + α K 3 j∈K, =i (h i -h j ) = γ i 1 6 u i • n i + 2α K 3 h i + j∈K,j =i γ i 1 6 u j • n j - α K 3 h j
we remark that we have put the LLxF uctuation into the form (1.32) which is particularly suited to nd local bounds of the numerical solution. Each 1. Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids element contributes with its residual to the update

h * i = h n i - ∆t |C i | φ LLxF i = 1 - ∆tγ i |C i | u i • n i 6 + 2α K 3 h i + j∈K,j =i ∆tγ i |C i | α K 3 - u j • n j 6 h n j = a K ii h n i + j∈K,j =i a K ij h n j h * i > 0 is positive if the extra-diagonal coecients a K ij > 0, which x a lower
bound to the LxF dissipation parameter

α K > 1 2 u j • n j , ∀j ∈ K, = i
We have neglected

∆tγ i |C i |
because it is always positive. Moreover it is necessary to enforce a CFL condition to have the positivity of the diagonal coecient [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF] shows that the PC scheme given by (1.36) and (1.37) with

a ii > 0 ∆t < 3|C i | K∈D i 1 2 u i • n i + 2α K
LLxF uctuation Φ LLxF i , veries the positivity of water depth h n+1 i ∀i ∈ T h , for the following choice of the LxF parameter and time step

α K = 1 2 max j∈K u j max j∈K n j , ∆t < min i∈T h |C i | K∈D i α K (1.38)
if the limitation procedure (1.34) is carried out equation by equation, that is R = L = I 3 in (1.34). [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF] uses this result to avoid negative depth and selects how to perform the limiting according to the water depth.

In particular when h 1 he switchs to a scalar version of LLxF scheme for each equation

R = L = I 3 if min j∈K h j < C H R, L else Well-Balancedness
For FV we have seen that bathymetric source terms are approximated in the same fashion of the uxes, in order to recover exactly at a discrete level the balance between topographic and hydrostatic terms. This idea have been incorporated into the residual approach almost naturally, as done for any other source terms, see section (1.6.3) and Ricchiuto et al. [2007]; [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] and references therein). In the lake at rest case u = [ η 0 -

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

1.6. Discrete approximation b 0 0 ] T , one should check that trivially Φ K = 0. That is, that the boundary integral of the hydrostatic term in the ux balances exactly the integral of the topographic source:

Φ K = ∂K F j (u h )n j ds + K S b (u h , b h ) = 0 ∂K 1 2 gh 2 h n j ds + 0 K gh h ∂b h ∂x i dx = 0 (1.39)
This is true under the hypotheses of exact integration w.r.t. the linear variation of depth and bathymetry. In practice we have transformed

∂K 1 2 gh 2 h n j ds = K gh h ∂h h ∂x i dx
and then we used a second order quadrature formula to approximate both the hydrostatic and the bathymetric term.

1.6.6 Wet/Dry cell treatment

The treatment of the wetting/drying phenomenon is crucial in many coastal applications. For example, in tsunami simulation, rundown and ooding stages are crucial and must be accurately reproduced to predict runup heights and inundated areas after the tsunami attack. Wetting and drying in Shallow

Water simulations is an active subject of research which is discussed thoroughly in Brufau et al. [2002,2004]; [START_REF] Delis | A robust highresolution nite volume scheme for the simulation of long waves over complex domains[END_REF] and [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] for the FV and RD methods respectively. The interested reader is referred to these references for all details. We limit ourselves to highlight the main diculties and to introduce the implementation details along with the reference. First, we dene a partially dry cell or wet/dry cell as an element of triangulation on which the water depth passes from a positive value to zero. We note that the treatment of these regions requires the introduction of two small quantities. The rst is a threshold value C H , such that a node is considered dry if h i ≤ C H . This artice avoids the computations of unphysical wet cells with 0 < h 1. The second, is a cut-o required to modify the mass uxes and velocities close to dry cells. This value will be denoted here by C U , and C H C U . C H is a small quantity compared to the scale of the phenomenon, for real scale simulation is in the order of C H ≤ 10 -4 m. The choice of C U is less trivial as it can depends on the particular problem, [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] relates this threshold coecient to the mesh size h K : the problem is cured by using, only for the dry nodes, a modied value of the bathymetry. In particular, in the computation of the element's residual (RD) or interface's ux (FV), we dene for wet-dry cells, the maximum water height at the wet nodes

C U = h 2 K L ref
η K max = max j∈K,h j >0 (h j + b j )
and based on this values we articially correct of the bathymetry value at dry nodes according to the following modication

b i = η K max , if b K i > η max
such that well-balanced is recovered.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Chapter 2

Well-Balanced conservative methods for the Shallow Water equations in ALE form

In the rst part of the chapter we introduce basic kinematics relationships that are involved in the formulation of the SWEs in Arbitrary Lagrangian Eulerian (ALE) framework. The reader can nd these formulas in almost every book on continuum mechanics, we refer for example to [START_REF] Chadwick | Continuum Mechanics: Concise Theory and Problems[END_REF].

In the second part we use the kinematics relationships to write the SWEs in a framework which is not Eulerian nor Lagrangian but it is arbitrary. We discuss dierent forms of the ALE-SWEs with respect to the accomplishment of Well-Balancedness and mass conservation at a discrete level. Finally it is presented a possible implementation of ALE Finite Volume and ALE Residual Distribution. A few test cases will conrm theoretical expectations.

Basic kinematics

Since the uid motion is columnar, we will use the word uid column or uid particle without distinction. A point Q occupies a certain position identied by the vector x which can be expressed in Cartesian coordinates with standard summation

x(Q) = x 1 e 1 + x 2 e 2 ≡ x i e i e 1 and e 2 dene the horizontal plane where the undisturbed uid lies. We dene also e 3 as the the vertical axis in the upward direction. Along e 3 , we measure the uid depth h(x) and the bathymetry b(x). The uid free surface follows as η = h + b, see gure 2.1. Within the Shallow Water context, we call material uid B, a set of uid columns that can be put in bijective correspondence with the points of a certain conguration B lying on a plane. We distinguish the particles P ∈ B from the points Q ∈ B. Please note 2.1. Basic kinematics Figure 2.1: Denition of axis x 1 /x 3 , water depth h, bathymetry b, free surface η with respect to an arbitrary origin. On the right, a water column/particle P .

that the bold font will be used in the following for vector and tensors. In the initial or reference conguration we use index (•) 0 to highlight that time is xed, t = 0. Although the notation changes from text to text, all along this manuscript, particles in the reference conguration are labeled by coordinates in capital letter X(Q 0 ), dierently from particles in a specic conguration that are labeled by low case letters x(Q). We assume the existence of a function such that

F 0 : B → B 0 , X = F 0 (P ) (2.1)
with ∃ F -1 0 . As sketched in gure (2.2), the body conguration changes with time. If, at each value of time is associated a unique conguration B(t), then the family of congurations is called motion of the uid B:

F : B → B, x = F (P, t) (2.2)
Combining relations (2.1) and (2.2) we obtain x = F (F -1 0 (X), t) or

B : B 0 → B x = B(X, t)
We assume the function B to be continuously dierentiable, thus a smooth Jacobian of transformation can be dened

J B = ∂x ∂X , J B = det J B
Particle collisions is avoided through the requirement ∃J -1 B , ∀P , or equivalently det J B = 0. Furthermore we add the physical condition that material volumes are positives, thus the determinant of the Jacobian will be strictly positive J B > 0. A sketch of these relations can be found in gure (2.2).

Well-Balanced conservative methods for the Shallow Water equations in ALE form

Lagrangian and Eulerian descriptions

Lagrangian and Eulerian descriptions

The uid B is endowed with various physical properties associated to the particles P . These properties, represented by scalar φ, vector v and tensor elds T , can be either dened with respect to the reference conguration or to the actual one. In the former case, called Lagrangian framework, the independent variable is X. For a scalar φ :

B → R or φ 0 : B 0 → R this corresponds to φ = φ • B = φ 0 (X, t) Lagrange (2.3)
The pullback operator is useful because it transfer the eld set from the actual description to the reference one, highlighting the transformations φ = φ • B : B 0 → B → R. The second part of the statement (2.3) is just a shorter notation for the same concept φ = φ 0 (X, t) : B 0 → R and will be used in the following. In presence of uid motion, the material derivative of φ is the time rate of change of φ measured by an observer which moves with the particle labeled with X, where X(Q 0 ) is the particle's position in the reference conguration.

We will shorten it

dφ dt = ∂φ ∂t X Lagrange (2.4)
Analogously, if we take x as independent variable

φ = φ 0 • B -1 = φ (x, t) Euler
and the time derivative with respect to a xed point labeled with x, is

∂φ ∂t = ∂φ ∂t x

Euler

The velocity of the particle at point Q, u = u j e j is calculated applying the aforementioned denition of Lagrangian derivative (2.4), to the position x(Q).

Its components write

u j = ∂x j ∂t X
Using the pullback one can write the scalar φ with respect to the actual coordinates and then derive:

dφ dt = ∂φ 0 ∂t X • B -1 = ∂φ ∂t x + ∂φ ∂x j ∂x j ∂t X = ∂φ ∂t x + ∂φ ∂x j u j
we get an important relationship between the Lagrangian time derivative and the Eulerian one. Similar arguments can be used to transport time derivatives

of a vector v = v i e i dv i dt = ∂v i ∂t X • B -1 = ∂v i ∂t x + ∂v i ∂x j u j

Well-Balanced conservative methods for the Shallow Water equations in ALE form

During uid motion the material volumes transforms according to dx m = J B dX m , with dx m = dx 1 dx 2 . From this fact stems the following relationship

∂J B ∂t X = J B ∂u j ∂x j (2.5) or, more compactly ∂ t | X J B -J B ∇ • u = 0.
The proof can be found in every textbook on continuum mechanics.

Arbitrary Lagrangian Eulerian kinematics

To embed adaptive mesh deformation in the numerical solution of (1.1)(1.2)

an appropriate Arbitrary Lagrangian Eulerian (ALE) formulation will be used.

The objective of the next section is to recall some basic aspects related to ALE.

For more details concerning the ALE formalism, the interested reader can refer to the original paper of [START_REF] Donea | Arbitrary lagrangian eulerian nite element methods[END_REF] or to a recent review Donea et al. [2004].

We introduce an arbitrary conguration A and we will refer to it as Arbitrary Lagrangian Eulerian (ALE). For us, A will be the conguration which the mesh is attached to. As for the body conguration B, it belongs to a Cartesian plane composed of points Q. The arbitrary conguration in the reference or initial state is denoted by A 0 and the points Q ∈ A 0 . If capital letters are used for the reference conguration, and lower case letters label the actual conguration, Greek letters are used to describe points in the ALE conguration, χ( Q). Moreover we are interested in arbitrary meshes attached to a conguration that changes with time A(t). If, to each value of time we associate a unique arbitrary conguration, we could assume the existence of a function

A : A 0 → A x = A(χ, t) (2.6)
As done for B, we require that A is continuously dierentiable, and consider the Jacobian matrix:

J A = ∂x ∂χ , J A = det J A > 0 ∀ Q (2.7)
which means that also A cannot admit fractures, collisions and negative volumes. From the point of view of a mesh cell, the occurrence of negative volumes is also known as tangling, an urgent problem in ALE simulations. For example, very often, when solving the motion of A(t), we cannot impose explicitly the positiveness of the J A , thus we cannot avoid mesh tangling a priori. This will be detailed in Chapter 3. For now, analogously to the particle/Lagrangian velocity, we introduce the mesh/ALE velocity, which is the velocity of the point Q which was lying, in the arbitrary reference domain, in χ( Q)

σ j = ∂x j ∂t χ

Arbitrary Lagrangian Eulerian kinematics

σ = σ j e j is arbitrary and independent of the ow (in general). In the previous paragraph we have transformed time derivative from the Lagrangian to the Eulerian description B 0 → B. Recalling briey that, for the scalar φ, the following relations holds

φ = φ • A = φ 0 (χ, t) ALE (2.8) φ = φ 0 • A -1 = φ (x, t) Euler (2.9)
We can do the transformation for the arbitrary conguration A 0 → A

∂φ ∂t χ = ∂φ 0 ∂t χ • A -1 = ∂φ ∂t x + ∂φ ∂x j ∂x j ∂t χ = ∂φ ∂t x + ∂φ ∂x j σ j (2.10)
and an interesting relationship emerge between Lagrangian and ALE time derivative:

dφ dt = ∂φ ∂t χ + ∂φ ∂x j u j -σ j (2.11)
The above relationship holds also for vector

dv i dt = ∂v i ∂t χ + ∂v i ∂x j u j -σ j (2.12)
The correct computation of ALE domain volumes in the transformation is assured by the relation dx m = J A dχ m or by the so-called Geometric Conservation Law (GCL), which generalize (2.5):

∂J A ∂t χ = J A ∂σ j ∂x j (2.13) More compactly ∂ t | X J A -J A ∇•σ = 0.
This can be interpreted as a constraint that ALE Jacobian and the domain velocity both must comply during the motion. As we will see in while, people involved in ALE simulations, are particularly interested to respect explicitly the GCL. If integrated over a close domain (2.13) represents a sort of volume/area conservation, which one must also preserve at the discrete level.

ALE remap

Lastly, we note that for a steady function in Eulerian framework φ = φ(x), we can write (2.10) imposing ∂ t φ = 0

∂φ ∂t χ = ∂φ ∂x j σ j (2.14)
This last equation represents the time variation of the function φ(x(t)) measured from an observer which is following the ALE domain motion x = A(χ, t).

Well-Balanced conservative methods for the Shallow Water equations in ALE form

Summing eq. (2.14) (pre-multiplied by J A ) to eq. (2.13) (pre-multiplied by φ) leads to the conservative ALE remap equation for the function φ

∂J A φ ∂t χ -J A ∂φσ j ∂x j = 0 (2.15)

Shallow Water equations in ALE form

We move away from a pure kinematic and inertial description of the ow, and we go back to hydrodynamics. Let us consider a shallow uid column P in B, with depth h and ow velocity u. Forces of dierent nature act on P .

In particular pressure term participates microscopically to momentum balance and they can be added to advective ux. The bottom topography exerts also contact forces on the uid, in particular pressure and frictional forces.

The transport formulas (2.11)(2.12) and the two volume conservation constraints, namely the uid volume (2.5) and the geometric conservation law (2.13), are used here to transform the eld equations

∂J B h ∂t X = 0 and ∂J B hu i ∂t X = J B F i
in a general arbitrary framework coincident with the domain motion. The SWEs write in this case (see also the Eulerian equations in (1.1)):

∂J A h ∂t χ + J A ∂ ∂x j hu j -hσ j = 0 (2.16) ∂J A hu i ∂t χ + J A ∂ ∂x j T ij -hu i σ j + J A S i = 0
(2.17) all vectors and tensors denitions hu i , T ij , S i were already given in section1.1 of chapter 1. Always as in chapter 1, we provide the vector form of the SWEs

∂J A u ∂t χ + J A ∂ ∂x j F j -σ j u + J A S(x, u) = 0, (2.18) u = h hu i , F j = hu j T ij , S = 0 gh ∂b ∂x i S b + 0 c F u i S f (2.19) Equations (2.18)(2.19
) is a non-homogeneous hyperbolic system of partial differential equations. In particular, given any vector ξ = ξ i e i the ux Jacobian

K ALE = ∂ ∂u (F j -σ j u)ξ j is K ALE (u, σ, ξ) = K -σ j ξ j I 3 (2.20)
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Discrete ALE equations

It admits a full set of real eigenvalues in the form λ(u, σ, ξ), namely

λ 1,3 = (u j -σ j )ξ j ± c ξ λ 2 = (u j -σ j )ξ j
and linearly independent eigenvectors that does not change with respect to the standard Eulerian SWEs ones.

Discrete ALE equations

The implementation of ALE Finite Element or Finite Volume poses a new issue: the preservation of time accuracy and non-linear stability properties on moving meshes. A naif extension of the xed grid algorithm (without worrying how to choose the mesh velocity, for example) does not assure that the GCL requirement is veried exactly and this may lead to spoil time accuracy and to the rise of numerical instabilities. [START_REF] Guillard | On the signicance of the geometric conservation law for ow computations on moving meshes[END_REF] proved that satisfying the GCL is a sucient condition for rst order time accurate simulation on moving meshes. We remark that, imposing a uniform ow in (2.18), gives back the GCL. [START_REF] Thomas | Geometric conservation law and its application to ow computations on moving grids[END_REF] proposed in fact to replace the GCL with the constraint of reproducing exactly uniform ows.

The compliance of the GCL discretely is referred to as Discrete GCL (DGCL).

In our case, a numerical method approximating solutions of (2.18)(2.19) on a moving mesh, is said to verify a DGCL if for S = 0, the state u = u 0 = const is an exact solution of the discrete equations. As we said, a numerical method veries the DGCL if it also embeds an exact discretization of the GCL (2.13).

Even if one could think that it is trivially important to conserve the total area of the computational domain, still the importance/benet of the GCL remains quite controversial, see the summery of [START_REF] Etienne | Geometric conservation law and nite element methods for 3d unsteady simulations of incompressible ow[END_REF]. This debate is important but, for the SWEs, the relevance of the DGCL lies in the fact that, while avoiding the appearance of spurious oscillations, it becomes a necessary condition to preserve some exact steady state on moving meshes, which we are interested to. DGCL schemes will be of key importance to retain

Well-Balancedness on moving meshes.

At a discrete level we work on integral equations, thus it is important to consider the integral form of the GCL. For example for a general volume V :

∂ ∂t χ V dx = V σ m n m ds (2.21)
where the left hand side is discretized with the time scheme. It is thus crucial to choose edge uxes (right hand side) such that the above expression is an identity. Of course at a discrete level, the closure of the problem will depend on the specic time scheme used. In general we can say that the unknown are the grid velocities at the edges and the conguration on which we integrate.

Regarding these quantities in particular, [START_REF] Mavriplis | Construction of the discrete geometric conservation law for high-order time-accurate simulations on dynamic meshes[END_REF] show that 2. Well-Balanced conservative methods for the Shallow Water equations in ALE form the satisfaction of the DGCL is not a sucient condition to have the prescribed high order time accuracy. In particular edge velocities and conguration should be chosen coherently with the time scheme. That is, using the quadrature points of the time integrator. Before going on we note that the grid velocity is not completely free. In particular besides the constraints on the positivity of the volume, we limit ourselves to a linear approximation of the mesh faces.

This means that only straight edges are allowed (no curved mesh). As a consequence we take the velocity space to be contained in the P 1 linear nite

element space σ = j∈T h ϕ i σ j (2.22)
with, as usual, {ϕ i } i∈T h being the standard piecewise linear P 1 Finite Element kernel. σ j is the mesh velocity vector of the node j of the triangulation. Still the grid velocity time history σ(t) is unknown and no hypothesis are made.

DGCL closure for eRK2

We consider explicit Runge Kutta type schemes, as the one presented in Chapter 1. For every RK-stage the time lapse is just one, we go from n to n + 1 and the time discrete approximation of (2.21) for a polygon V , namely the DGCL, writes

V dx n+1 - V dx n = ∆t ∂V σ m n m ds
We can further manipulate, decomposing the left-hand side in the contributions of each volume's edge

|V n+1 | -|V n | = ∆t j∈V υ j (2.23)
where |V | is the area of the polygon, j ∈ V is an index over the set of polygon's edges. The edges' velocities result:

υ j = ∂V j σ m n m j ds (2.24)
and they are denoted, in a Finite Volume context, as interface velocities. The question is how to verify (2.23) exactly, which means with an error driven by machine roundo ? [START_REF] Mavriplis | Construction of the discrete geometric conservation law for high-order time-accurate simulations on dynamic meshes[END_REF] and [START_REF] Isola | Arbitrary lagrangian eulerian formulation for two-dimensional ows using dynamic meshes with edge swapping[END_REF], propose to evaluate geometrically the edge velocities from the volume swept by the corresponding interface. We will use extensively their closure. We state the main conclusion in the form of the following preposition.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Combining DGCL and Well-Balancedness on moving meshes

Proposition (DGCL). A numerical method approximating eq. ( 2.18) veries the DGCL constraint eq. ( 2.23), if each edge velocity is computed through the (signed) area swept by the j-th edge of the polygon in one timestep ∆t

υ j = 1 ∆t Q j dx (2.25)
where the quadrangle Q j is dened by vertexes

x n i ,x n k ,x n+1 k ,x n+1 i
, i and k denotes the two vertexes of the j-th edge. We use the convention that, if the polygon is expanding, then the area swept is positive. [START_REF] Isola | An Interpolation Free Two-Dimensional Conservative ALE scheme over Adaptive Unstructured Grids for Rotorcraft Aerodynamics[END_REF] explicitly developed the formula (2.25), and found the popular midpoint closure proposed by [START_REF] Lesoinne | Geometric conservation laws for ow problems with moving boundaries and deformable meshs, and their impact on aeroelastic computations[END_REF]. We summarize the main result contained in the last reference: a numerical method approximating eq. ( 2.18) veries the DGCL constraint or eq. (2.23) if the midpoint conguration V n+1/2 and constant nodal grid velocities are used to computed the ALE ux, that is replacing in (2.24) the following expressions

n m j = n m,n j + n m,n+1 j 2 , σ m j = x m,n+1 j -x m,n j ∆t (2.26)
For eRK2, the closure of Mavripils and Fahart collapses. However formally, we prefer to maintains the notation of [START_REF] Mavriplis | Construction of the discrete geometric conservation law for high-order time-accurate simulations on dynamic meshes[END_REF] and (2.25) because, in chapter 4, it will prove to be more exible.

In the Finite Element case, the DGCL and the hypothesis (2.22) (σ ∈ P 1 ), dene a constant mesh velocity divergence on the element. We take our control polygon as the element itself V = K. We x a conguration K t = K(t); the DGCL writes

|K n+1 | -|K n | = ∆t K t ∂σ j ∂x j dx = ∆t K dx ∂σ j ∂x j K
and we obtain a rule to compute the elemental divergence of the mesh velocity vector:

∂σ j ∂x j K = ∆|K| ∆t|K t | (2.27)

Combining DGCL and Well-Balancedness on moving meshes

For a balance law, the source term S is dierent from zero, and the relevant state to be preserved may not be u = const but η = const. Note in particular, 2. Well-Balanced conservative methods for the Shallow Water equations in ALE form that we may write (2.18) as

J A ∂u ∂t χ -σ j ∂u ∂x j H 1 +u ∂J A ∂t χ - ∂σ j ∂x j H 2 +J A ∂F j ∂x j + S H 3 = 0
Existing Eulerian discretization methods do embed integral (or even local)

variants of H 3 = 0 (C-
H 1 + ∂b ∂t X -σ j ∂b ∂x j H 4 = ∂η ∂t χ -σ j ∂η ∂x j
which is of course null when η is the invariant associated to the equilibrium H 3 = 0.

This suggests that a better form of (2.17) for computations on moving meshes, is that obtained by summing Eq. (2.15) to Eq. (2.17). This leads to a

Well-Balanced ALE form of the problem reading

∂J A w ∂t χ + J A ∂ ∂x j F j -σ j w + J A S(x, u) = 0 (2.28) w = η hu i (2.29)
In this case one can do much better job in the approximation of the lake at rest solution. In particular, we can write (2.28) as

J A ∂w ∂t χ -σ j ∂w ∂x j H 1 +H 4 +w ∂J A ∂t χ - ∂σ j ∂x j H 2 +J A ∂F j ∂x j + S H 3 = 0
If η is constant any Eulerian method will be able to embed the condition H 1 + H 4 = 0 while, choosing appropriate schemes verifying both the DGCL and the WB, we will be able to satisfy all the compatibility requirements, and preserve steady equilibria independently on the mesh movement strategy.

As a particular case and for completeness, we recall the pre-balanced form of the Shallow Water equations of [START_REF] Rogers | Mathemathical balancing of ux gradient and source terms prior to using roe's approximate riemann solver[END_REF] which is obtained by

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Mass conservation on moving meshes

introducing in (2.28) the modied ux and source functions

Fj = hu j T ij , Sb = 0 gη ∂b ∂x i (2.30) with T ij = huu ij + 1 2 g (η 2 -2ηb
). These denitions lead to the pre-balanced (PB) ALE form of the shallow-water system reading

∂J A w ∂t χ + J A ∂ ∂x j Fj -σ j w + J A S(x, u) = 0 (2.31)
We recall that the this form satises the relation, see [START_REF] Rogers | Mathemathical balancing of ux gradient and source terms prior to using roe's approximate riemann solver[END_REF] ∂ F(w; b) ∂w = ∂F(u) ∂u so the pre-balanced system has the same eigenstructure of the standard one.

Mass conservation on moving meshes

We now consider the additional constraint of conserving the total water mass in the domain. We integrate, in space and in time, the mass conservation equation in Well-Balanced form

Ω(t) η(x(t), t) dx - Ω X η(X, 0) dx + ∂Ω(t)
(hu j -ησ j )n j ds dt = 0

Let H(t) be the total mass of water at time t, H(t) = Ω(t) h dx, and dene B(t) = Ω(t) b dx, we can rewrite mass conservation statement separating the terms in h and the terms in b

H(t)-H(0)+ ∂Ω(t) (hu j -hσ j )n j dsdt+B(t)-B(0)- ∂Ω(t)
bσ j n j ds dt = 0

(2.32) which states that, modulo the boundary conditions, we have conservation over the full domain if the ALE remap equation (2.15) is satised, namely if

B(t) -B(0) - ∂Ω(t)
bσ j n j ds dt = 0

So a scheme approximating (2.28) will be exactly mass conservative only if the bathymetry is evolved according to an integral form of an ALE remap (2.15). This is the strategy proposed in Zhou et al. [2013a]. However, as pointed out in the same paper, this approach leads to changes in the bathymetric altitudes which will depend on the scheme. For example, substantial smoothing of 2. Well-Balanced conservative methods for the Shallow Water equations in ALE form the bed slopes is observed. To deal with this issue, in Zhou et al. [2013a] the authors propose to regularly re-initialize the bathymetric data. This will however violate ( 14), and so a mass loss will be associated to each of these re-initialization steps. Here we propose an alternative solution, allowing to preserve mass down to almost machine accuracy. Assume for simplicity that the domain boundaries are not moving, or that σ • n = 0 is veried. We can write the mass error at time t as

E mass = H(t) -H(0) + ∂Ω hu j n j ds dt = B(0) -B(t)
We now remark that the two quantities on the right hand side are in principle equal, as they are both approximations of the integral of b(x) over the domain. If the domain boundaries are not moving, this quantity should remain constant in time. In practice however, these two integrals will be evaluated on a moving mesh. This means that, even if both the domain of integration and the data being integrated are constant, the quadrature points used will move, so the result will not be the same. To be more precise, with the numerical approximation, B(t) will be split in integrals over the set of median dual cell areas. In case of standard FV or P 1 RD discretizations

B(t) = i∈T h C i (t) b dx = i∈T h b i (t)|C i (t)| (2.33) with b i = b(x i (t)). Our idea is to compute dierent nodal values b i = b(x i (t))
such that the total mass error is reduced by simply increasing the accuracy which the elemental integrals are evaluated with. So we will set

C i (t) b(x(t)) dx ≈ |C i (t)| Nq f =1 ω q b(x q (t))
and we we will compute bathymetric nodal values

b i = Nq f =1 ω q b(x q (t))
(2.34) in gure 2.3 we provide a sketch of the treatment that we propose for the bathymetry. Here it is essential to underline that b(x q ), on the right hand side, is a given high accurate (analytical or reference one, interpolated on a ne mesh) representation of the bathymetry. In case we employ only one quadrature point that coincide with the node x q = x i we get back to the standard choice (2.33).

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. 2. Well-Balanced conservative methods for the Shallow Water equations in ALE form

Mass conservation and Well-Balancedness in presence of wet/dry cells

In presence of wetting/drying, there is a major complication related to the fact that the volume containing water mass is moving, and its movement is a-priori independent on the mapping dened by σ. Of course, in this case the ow equations are only solved in the wet region Ω W (t). So mass conservation now reads

Ω W (t) η(x(t), t) dx - Ω W,X η(X, 0) dx + ∂Ω W (t)
hu j -ησ j n j ds dt = 0

Water depth and discharge are both null at the shoreline ∂I W while the ALE ux is null at the domain boundaries. So we choose to dene two boundary

regions ∂Ω W = ∂Ω ∩ ∂Ω W + ∂I W and write H(t)-H(0)+ ∂Ω∩∂Ω W hu j n j ds dt+B W (t)-B W (0)- ∂I W (t)
bσ j n j ds dt = 0

The mass error, due to the deformation of the computational domain, becomes

E mass = -B W (t) -B W (0) - ∂I W (t)
σ j n j ds dt

As before, this quantity is not zero, as we do not use the ALE remap to evolve the (given) bathymetry. To link this residual error to the previous case, we add and remove the following quantity:

Q = B D (t) -B D (0) - ∂I D (t)
bσ j n j ds dt the sub-script • D denoting integrals over the dry area. We nally obtain

E mass = B(0) -B(t) + Q
The dierence between the rst two terms can be reduced as discussed before.

The reminder Q is a geometrical term associated to the deformation in dry areas. Unfortunately, we are not able to guarantee any a-priori control on this term, since, as we will see later, grid adaptation w.r.t. the shoreline benets from the possibility of exploiting points in the dry region. In this paper, this geometrical factor arising from deformation in dry areas will be accounted for by uniformly redistributing the mass excess/defect in the wet region.

A second non trivial issue related to wetting/drying and moving meshes is the Well-Balancedness. To guarantee that the mesh movement does not spoil the preservation of the lake at rest state close to partially wet cells, an ad-hoc treatment is introduced. This procedure impacts the way in which the new Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

η max,i = max K∈D i max j∈K h j >C H η j C H was introduced in section 1.6.6 as a small threshold (C H 1) to dene dry nodes; if h i ≤ C H , then h i = 0. 2. ∀i set ∆ b i = b * i -b n i (2.35) with b * i = η max,i if i is dry and b n i + C H > η max,i b n i otherwise 3.
Compute the new water depth as

h n+1 i = max(0, η n+1 i -b n+1 i + ∆ b i ) (2.36)
In the lake at rest, the value of η n+1 obtained from the discretization of the system 2.28 or equivalently 2.31 and to be substituted in (2.36) is η n+1 i = η n i ∀i ∈ T h (the ux/residuals are zero everywhere). As represented in gure 2.4

2. Well-Balanced conservative methods for the Shallow Water equations in ALE form there are basically two situations for the lake at rest case in presence of a mesh which is arbitrary moving. First a node can move from the wet side the dry side of the interface. From step 2, a dry node has ∆ b i = 0 (step 2 is not active).

So we have only step 3 that writes

h n+1 i = max(0, η n+1 i -b n+1 i + ∆ b i ) = max(0, η n i -b n+1 i <0 ) = 0
which is the correct value in the dry areas. Secondly there is the vice-versa:

a node moves from the dry to the wet side of the interface. This time, step 2 enforces the lake at rest:

h n+1 i = η n+1 i -b n+1 i + ∆ b i = b n i -b n+1 i + η max,i -b n i = η max,i -b n+1 i
which is the exact water depth. This guarantees that a constant at free surface level is exactly preserved also near shorelines. In gure 2.4 there is a graphical representation of the lake at rest preservation.

Notation for time dependent geometry

In the ALE case the ow unknowns have both explicit and implicit dependence on time u i (t) = u(x i (t), t). A discrete evaluation, for example in t n , will be denoted by u n i = u(x i (t n ), t n ). Also geometrical quantities and physical data will change in time according to the movement of the mesh. Keeping the same notation we have

C n i = C i (t n ), ∂C n i = ∂C i (t n ) and b n = b(x(t n )).

ALE Finite Volume for SWEs

ALE-FV are of great interest for aerodynamics problems involving moving boundaries such as aeroelastic simulations, see for example [START_REF] Lesoinne | Geometric conservation laws for ow problems with moving boundaries and deformable meshs, and their impact on aeroelastic computations[END_REF]. ALE methods become succeful also in the context of mesh adaptation since they can embed conservation very eciently for both r/h-adaptation,

see [START_REF] Isola | An Interpolation Free Two-Dimensional Conservative ALE scheme over Adaptive Unstructured Grids for Rotorcraft Aerodynamics[END_REF]. Instead, for the SWEs, ALE schemes are not very popular.

We believe that the main reason is that many scalar quantities involved in the SWEs source term are inherently Eulerian, that is "xed to the ground". propose an ALE extension to approximate equations in the form (2.28). We will show the equivalence between the scheme obtained using the well-balanced form (2.28), and the pre-balanced formulation (2.31) with the denition (2.30).

The FV discrete evolution equations reads

|C n+1 i |w * i = |C n i |w n i -∆t j∈D i R ij (w n , b n ) (2.37) |C n+1 i |w n+1 i = |C n i |w n i - ∆t 2 j∈D i R ij (u n , b n ) 2 + R ij (u * , b n+1 ) 2 (2.38)
where we have again

R i = F ij + S ij (2.39)
As in chapter 1, we use the Roe-type numerical ux, this time for the well balanced formulations (2.28). ALE uxes are evaluated at the interface:

F ALE ij = 1 ∆t t n+1 t n ∂C ij wσ m n m ds dt = 1 ∆t t n+1 t n
wυ ij dt and they can be incorporated automatically in the resolution of the Riemann problem:

F ij (u i , u j ) = F i • n ij -w i υ ij + m p=1,αp<0 W p = F j • n ij -w j υ ij - m p=1,αp>0
W p keeping in mind that, when decomposing the solution jump (1.6.2), one should consider the velocity of simple waves is modied according to the ALE eigenvalues, see for example [START_REF] Isola | An Interpolation Free Two-Dimensional Conservative ALE scheme over Adaptive Unstructured Grids for Rotorcraft Aerodynamics[END_REF]:

m p=1,αp<0 W p = (K ij -υ ij I 3 ) -(u j -u i ) m p=1,αp>0 W p = (K ij -υ ij I 3 ) + (u j -u i )
with I 3 the 3×3 identity matrix. At the end, the numerical ux (1.17) will result enhanced by:

F ij = F ij (ȗ i , ȗj ; bi , bj ) = F(ȗ j ) + F(ȗ i ) 2 •n ij -υ ij wj + wi 2 - K ij -υ ij I 3 2 (ȗ j -ȗ i ) (2.40)
2. Well-Balanced conservative methods for the Shallow Water equations in ALE form

∂V j = ∂C ij : |C n+1 i | -|C n i | = ∆t j∈D i υ ij (2.41)
The following interface interface velocity stems from the application of (2.25)

υ ij = ∂C ij σ m n m ds = K i,j ∂C K ij σ m n m ds = 1 ∆t K i,j Q K ij dx (2.42) with the quadrilateral Q K ij dened by vertexes x n G K ,x n ij ,x n+1 ij ,x n+1 G K , see gure 2.5.

Well Balancedness

The treatment of the bathymetric term does not change substantially with respect to the xed mesh case. The denitions (1.21),(1.22) and (1.23) are valid. We only have to modify the upwind part (1.24) as follows: Proof. For constant b, and constant u 0 , the discrete equations reduce to

S * ij = - sign(K ij -υ ij I 3 ) 2 (A ij -υ ij I 3 )∆b ij (2.
(2.41), which proves the rst part (DGCL)

|C n+1 i |u n+1 i -|C n i |u 0 = ∆t j∈D i υ ij u 0 ⇒ u n+1 i = u 0
For the second part, the proof rests on the property of the Roe average and on the fact that, on the lake at rest state, we have K ij = A ij . In particular, proceeding as in chapter 1

F ij = F(ȗ j ) -F(ȗ i ) 2 • n ij + (F(ȗ i ) -F(u i )) • n ij -υ ij w 0 - |A ij -υ ij I 3 | 2 (ȗ j -ȗi ) S b ij = A - ij ∆b - ij + 1 2 A ij ∆b ij - sign(A ij -υ ij I 3 ) 2 (A ij -υ ij I 3 )∆b ij
Note now that ȗjȗi + ∆b ij = wjwi which vanishes by hypothesis, so that the last two terms cancel each other. The rest of the proof is almost identical to the scalar case, and uses the fact that, on the selected equilibrium,

(F(ȗ j ) -F(ȗ i )) • n ij = A ij (ȗ j -ȗi )
and the constancy of w = w 0 .

Source term upwinding: a bridge between WB and PB form

Before moving on, it is interesting to note that the use of the FV discrete equations obtained by using the pre-balanced form of the shallow equations (2.31) are almost identical to those presented above which are instead derived from (2.28). We neglect friction that does not play any role in deriving the results that follows. First we dene

R ij = Fij + Sb ij (2.44)
Roe numerical ux for the PB formulation:

Fij = F(w j ; bj ) + F(w i ; bi ) 2 • n ij -υ ij wj + wi 2 - K ij -υ ij I 3 2 (w j -wi ) (2.45)
and for later purposes, we dene a non upwind numerical source term for the PB equations:

Sb ij = Ã- ij ∆b - ij + 1 2 Ãij ∆b ij (2.46)
2. Well-Balanced conservative methods for the Shallow Water equations in ALE form with (cf. equation (1.5)) Ãij = A(n ij , η ij ) where, in analogy with the notation used so far, η ij = (η j + ηi )/2. Similarly, we can also dene Ã-

ij = A(n ij , η - ij ) with η - ij = (η i + η i )/2.
We observe the following equivalence:

Proposition 2. The Pre-Balanced upwind FV discretization obtained from the pre-balanced form of the SWEs (2.31) with numerical uxes (2.45) and a non-upwind source term approximation (2.46) is equivalent to a minor modication of the Well-Balanced scheme which is always given by (2.37)-(2.38) (with all the denitions contained) but setting in (1.21) the upwind source term as:

S * ij = - sign(K ij -υ ij I 3 ) 2 (K ij -υ ij I 3 )∆b ij (2.47)
instead of (2.43).

Proof. We rewrite the well balanced FV spatial discretization incorporating the new denition of source upwinding: (2.47) reads:

R ij =F ij + S b ij = F(ȗ j ) -F(ȗ i ) 2 • n ij + (F(ȗ i ) -F(u i )) • n ij -υ ij wj + wi 2 - |K ij -υ ij I 3 | 2 (ȗ j -ȗi ) +A - ij ∆b - ij + 1 2 A ij ∆b ij + S * ij (2.48)
We now use the fact that denition (2.47) of S * ij is such that when added to dissipation of the numerical ux one gets

- |K ij -υ ij I 3 | 2 (ȗ j -ȗi ) + S * ij = - |K ij -υ ij I 3 | 2 (w j -wi ) (2.49)
We exploit the continuous equivalence

∂F j (u) ∂x j + S b = ∂ Fj (w, b) ∂x j + Sb
which is written here at discrete level:

(F(ȗ j ) -F(ȗ i )) • n ij + A ij ∆b ij = F(w j ; bj ) -F(w i ; bi ) + Ãij ∆b ij (2.50) (F(ȗ i ) -F(u i )) • n ij + A - ij ∆b - ij = F(w i ; bi ) -F(w i ; b i ) + Ã- ij ∆b - ij (2.
R ij = Fij + Sb ij (2.52)
with Fij dened as in (2.45) and Sij as in (2.46) which is exactly the pre- balanced FV discretization otained from (2.31) (cf. [START_REF] Rogers | Mathemathical balancing of ux gradient and source terms prior to using roe's approximate riemann solver[END_REF]; [START_REF] Liang | Adaptive quadtree simulation of shallow ows with wet-dry fronts over complex topography[END_REF]; [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]).

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

ALE Residual Distribution for SWEs

The last proposition shows that the well balanced discretization of [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF]; [START_REF] Hubbard | Flux dierence splitting and the balancing of source terms and ux gradients[END_REF] is equivalent to the use of the pre-balanced form of the equation for a particular choice of the upwind component of the source. The proposition also shows that another viable alternative would be for example

S * ij = - |A ij -υ ij I 3 | 2 ∆b ij (2.53)
which also leads to a well-balanced discretization (cf. proof of proposition 1).

Mass conservation on moving mesh

The bathymetric values b n+1 i (in the previous paragraph time dependency has been dropped for clarity) are computed as explained in the previous section, cf. eq. (2.34). For the FV method, the nodal values are obtained via Gaussian quadrature formula over each sub-triangle

C K ij b n+1 i = 1 |C n+1 i | j∈D i K i,j |K n+1 | 6 Nq q=1
ω q b n+1 q given b q = b(x q ) with x q = j∈K λ q j x j and y q = j∈K λ q j y j . The baricentric coordinates of the quadrature points λ q j are dened over the sub-triangles C K ij .

The one point quadrature with baricentric coordinate in i corresponds to a constant approximation of the bathymetry function over the median dual cell (zero order, r = 0) and coincides with the standard choice b n+1 i = b(x n+1 i ).

In the numerical experiments we will test the impact of rst and second order accurate formulas (denoted respectively r = 1, 2), in order to arbitrary decrease the mass error.

ALE Residual Distribution for SWEs

While rst ALE FV and FE methods dates back to many decades, RD scheme in ALE framework are quite recent. [START_REF] Michler | An arbitrary lagrangian eulerian formulation for residual distribution schemes on moving grids[END_REF] achieved rst order of accuracy with an Explicit Euler time integrator and later [START_REF] Dobes | Second order blended multidimensional upwind residual distribution scheme for steady and unsteady computations[END_REF] moved to high order time approximation such as BDF3 and Crank

Nicholson. In these references the common startegy to approximate the ALE term is to split it into two contributions ∂(σ j w) ∂x j = σ j ∂w ∂x j + w ∂σ j ∂x j 2. Well-Balanced conservative methods for the Shallow Water equations in ALE form [START_REF] Michler | An arbitrary lagrangian eulerian formulation for residual distribution schemes on moving grids[END_REF] referred to this term as Geometric Source Term and discretize it in such a way to verify the DGCL. We will take advantage of the above form and write the uctuation as:

Φ K (w h , b h ) = ∂K F j (u h )n j ds - K σ j ∂w ∂x j dx + K S(u h , x) dx (2.54)
As for FV, the uctuation is evaluated at the midpoint conguration. Thus, if not specied, for us K = 1 2 (K n + K n+1 ) in the following. Keeping in mind this, we provide directly the ALE extension of the eRK2-RD algorithm presented in section 1.6.4. We remark that it is in a form very close to the update (1.36)-(1.37). The only remarkable dierence is that the variable w appears instead of u; this is a consequence of the results of section 2.12.1 which are applicable to any scheme and thus holds also for RD. In order to verify WB a numerical method should approximate the Well-Balanced form of the SWEs (2.28):

1] Predictor step: for each element K ∈ T h

• Compute the residual Φ K(1) = Φ K (w n h , b n h ).
• Distribute the uctuation to the nodes of K such that j∈K Φ

K(1) j = Φ K(1)
• Compute the rst order prediction of the solution, denoted as w *

w * i = w n i - ∆t |C n+1 i | K∈D i Φ K(1) i (w n h , b n h ) (2.55) 2] Corrector step: for each element K ∈ T h • Compute the residual Φ K(2) = 1 2 Φ K (w n h , b n h ) + Φ K (w * h , b n+1 h ) + K w * h -w n h ∆t dx • Distribute the uctuation to the nodes of K such that j∈K Φ K(2) j = Φ K(2)
• Compute the second order correction from

w n+1 i = w * i - ∆t |C n+1 i | K∈D i Φ K(2) i (w * h , w n h , b n h , b n+1 h ) (2.56)
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

ALE Residual Distribution for SWEs

In Appendix D we provide the details that brought us to the above form which essentially comes from an analogy with stabilized Finite Element method.

With respect to the xed mesh algorithm, we are left with a new term into the uctuation (4.7) related to the ALE

Φ K,ALE = K σ j ∂w ∂x j dx
The problem is that we have the elements' edge velocities rather then nodal values, see section 2.5.1. As remarked by [START_REF] Etienne | Geometric conservation law and nite element methods for 3d unsteady simulations of incompressible ow[END_REF] this is a common problem for Finite Element methods since the DGCL provides a rule to compute the divergence of the mesh velocity but not the velocity itself. Still we can compute the above term as:

Φ K,ALE = K ∂σ j w ∂x j -w ∂σ j ∂x j dx = ∂K wσ j n j ds - K w ∂σ j ∂x j dx
We examine the rst part. Edge uxes are computed by quadrature formula.

At the same time, these uxes should verify the DGCL. Following the closure proposed in 2.5, we replace in (2.23) and (2.24), V = K and ∂V j = ∂K j and the DGCL becomes:

|K n+1 | -|K n | = ∆t j∈K υ j (2.57)
The velocity of the j-th edge (opposed to node j) comes from (2.25)

υ j = ∂K j σ m n m ds = 1 ∆t Q j dx (2.58)
with the quadrilateral Q j dened by vertexes

x n i ,x n k ,x n+1 k ,x n+1 i , see gure 2.6.
Once υ j are recovered, we share these velocities on the edge's quadrature points through Gaussian weights q = 1,

N Q ∂K wσ j n j ds = j∈K N Q q ω q w q υ j (2.59)
The second term is computed using the constant divergence statement (4.43)

K w ∂σ j ∂x j dx = K w dx ∂σ j ∂x j K = j∈K ϕ j dx w j ∂σ j ∂x j K = j∈K w j 3 ∆|K| ∆t (2.60)
where a weighted average solution on the cell appears.

Well-Balanced conservative methods for the Shallow Water equations in ALE form

Figure 2.6: DGCL for ALE-RD. Area swept by the edge ∂K j .

Well-Balanced on moving mesh

We give directly the following statement, which can be regarded as the extension of the Well-Balanced property for RD to moving meshes:

Proposition 3. The explicit predictor corrector residual distribution prototype (2.55), (2.56) with edge velocities (2.58) veries the DGCL for constant b, and the Well-Balanced property both on moving and xed meshes, provided that the same linear piecewise continuous approximation is used for w, b and consequently for u, and that all integrals involving these quantities are evaluated exactly w.r.t. this variation.

Proof. To prove the rst part (DGCL), we check that for constant bathymetry and no friction, the splitting terms on the right-hand sides in (2.55) and (2.56) are identically zero for a given constant state u 0 . We do the proof only for corrector, the predictor is a particular case. We use the hypothesis that quadrature formula for line integration are exact for P1 function, then the obvious ∂K n j ds = 0 and the property of Gaussian weights N Q q ω q = 1:

Φ K(2) = K u 0 -u 0 ∆t dx + F j 0 ∂K n j ds -u 0 j∈K N Q q ω q υ j + u 0 ∆|K| ∆t = u 0 ∆|K| ∆t - j∈K υ j = 0
This is immediately shown if we compute the edge velocities with (2.58). We have u n+1 i = u 0 .

To prove the second part of the proposition (Well-Balancedness), we proceed in an identical manner, except that now we assume that we have a con-Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

ALE Residual Distribution for SWEs

stant state w = w 0 . As before we consider only the corrector

Φ K(2) = K w 0 -w 0 ∆t dx + ∂K F j n j ds + ∂K S b dx -w 0 j∈K N Q q ω q υ j + w 0 ∆|K| ∆t = w 0 ∆|K| ∆t - j∈K υ j = 0

Mass conservation on moving mesh

In section 2.7 we have discussed the fact that mass conservation on moving meshes is related on the accuracy of the bathymetric values (2.34). For RD we used the following L 2 type projection

b n+1 i = 1 |C n+1 i | K∈D i |K n+1 | Nq q=1 ω q b n+1 q ϕ q
The three points quadrature with baricentric coordinates in the triangle's vertex corresponds to a piecewise linear approximation of the bathymetry function over the triangles (rst order, r = 1) and coincides with the standard choice b

n+1 i = b(x n+1 i ).
In the numerical experiments we will test the impact of second and third order accurate formulas (denoted respectively r = 2, 3).

Distributions

The distributions matrix can be extended directly from chapter 1. The SWEs Jacobian changed, see (4.8), and this must be taken into the distribution operator, especially for the SUPG matrix contained in (1.35), and the Lax Friederich parameter in (1.34).

Water depth positivity

The ALE-RD algorithm (2.55)(2.56) is in a form that closely resemble the Eulerian-xed mesh (1.36)(1.37). The result of water depth positivity preserving can be extended using the same arguments of the proof contained in [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF]. There are however two issues that must be taken into account. First, the uctuation is enhanced by the ALE term. Second, the mass equation is written in the variable η and not h. The predictor update for the mass equation from (2.55), followed by the WB correction (2.34) rewrites:

η * i = η n i - ∆t |C n+1 i | K∈D i φ K i (2.61) h * i = max(0, η * i -b n+1 i + ∆ b i ) (2.62)

Well-Balanced conservative methods for the Shallow Water equations in ALE form

where ∆ b i is dened in (2.35) and it is dierent from zero only if i is dry.

As a consequence of the max operator, introduced to avoid negative depths in the lake at rest, such an update is trivially positive. However let's try to understand which is exactly the role of the max operator for a general situation where the ow is not at rest. We rewrite the update (2.61) in h and we separate the ALE uctuation in h and b. Always for the predictor step we have:

ĥ * i = h n i - ∆t |C n+1 i | K∈D i φ h i I -b n+1 i + b n i - ∆t |C n+1 i | K∈D i φ b i II =0
(2.63)

h * i = max(0, ĥ * i + ∆ b i ) (2.64)
where in every cell K, φ h i is a splitting of the residual in h

φ h = ∂K hu j n j ds - K σ j ∂h ∂x j dx
b ALE i = b n i - ∆t |C n+1 i | K∈D i φ b i
We can dene II as:

II = b ALE i -b n+1 i = 0
The rst two terms on the right hand side of (2.63), named as I, have not been addressed yet. Basically they represent the update as if we were evolving the ALE mass equation in h, that is (2.16). With respect to solving the SWEs in the water depth, the term II can be seen as a perturbation related to the fact that we are not solving an equation for the bathymetry remap. Unfortunately, we do not control the sign of this perturbation. Very roughly, in (2.62)

h * i = max(0, I + II + ∆ b i )
if II + ∆ b i generates a negative water depth, the max operator nullies this contribution automatically.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

ALE Residual Distribution for SWEs

However, it is important to ensure that I, the part of the mass equation update in h, is still positivity preserving. This provides our ALE scheme of a CFL condition that avoids negative depths in case of at bathymetry and ensures that negative contributions, which are roughly nullied in (2.62), are only due to a non-exact ALE remap of the bathymetry. The uctuation φ h writes:

φ h = 1 2 j∈K h j u j • n j - j∈K h j υ i + υ k 2 - 1 3 ∆|K| ∆t = 1 2 j∈K h j u j • n j - 1 3 j∈K h j υ i + υ k 2 -υ j
in last development we have used the DGCL (2.57). With respect to the xed mesh case of section 1.6.5 the ALE uctuation is enhanced by the ALE term which have been computed according to (2.59) and (2.60). We set for simplicity:

υ i + υ k 2 -υ j = υ ik -υ j
As done in section 1.6.5 we write the LLxF uctuation as:

φ LLxF i = γ i 1 6 j∈K h j u j • n j - 1 9 j∈K h j (υ ik -υ j ) + α K 3 j∈K, =i (h i -h j ) = γ i 1 6 u i • n i - 1 9 (υ jk -υ i ) + 2α K 3 h i + j∈K, =i γ i 1 6 u j • n j - 1 9 (υ ik -υ j ) - α K 3 h j
for the denition of γ i > 0 see always 1.6.5. We remark that we have put the LLxF uctuation into the form (1.32) which is particularly suited to nd local bounds of the numerical solution. Considering the contribution of each element to the update (2.63) separately:

ĥ * i = h n i - ∆t |C n+1 i | φ LLxF i = 1 - ∆tγ i |C n+1 i | u i • n i 6 - υ jk -υ i 9 + 2α K 3 h n i + j∈K,j =i ∆tγ i |C n+1 i | α K 3 - u j • n j 6 + υ ik -υ j 9 h n j = a K ii h n i + j∈K,j =i a K ij h n
h * i > 0 is positive if the extra-diagonal coecients a K ij > 0, which x a lower
bound to the LxF dissipation parameter:

α K > 1 2 u j • n j - υ ik -υ j 3 , ∀j ∈ K, = i We have neglected ∆tγ i |C i | because it is always positive. Moreover it is necessary
to enforce a CFL condition to have the positivity of the diagonal coecient

a ii > 0 ∆t < 3|C n+1 i | K∈D i 1 2 u i • n i - υ ik -υ j 3 + 2α K
For the corrector step, the positivity analysis can be carried out exactly in the same fashion of [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF] to which we refer for details (adding of course the ALE uctuations for the predictor and corrector step). In practical calculation we have imposed:

α K = 1 2 max j∈K u j max j∈K n j + 2 max j∈K υ j , ∆t < min i∈T h |C n+1 i | K∈D i α K 2.
12 Numerical experiments

Well Balanced

To illustrate some concepts and to better highlight certain numerical eects, we reconsider the simplied model (2.65). We add a general source to mimick the SWEs. This model reads

∂u ∂t + ∂F j ∂x j + S(u, x) = 0 , x ∈ [0, 2] × [0, 1], t ∈ [0, 1] (2.65)
where, for a given ux vector F (u) = [ F 1 F 2 ], the source term is dened as S = a j ∂b ∂x j

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. To test the WB property on moving mesh for the ALE formulation, we use the simple case of linear advection of a smooth sinusoidal hill

η 0 (x) = 1 + cos 2 (2πr) if r ≤ 0.25, r = (x -0.5) 2 + (y -0.5) 2 η 0 (x) = 1 otherwise
A pseudo-bathimetry is dened by b(x) = 0.8e ψ(x,y) with ψ = -5 (y -0.9) 2 -50 (x -0.5) 2 . The following arbitrary mapping is used to move the mesh

x i (t) = X i + 0.1 sin (2πX 1 ) sin (πX 2 ) sin (2πt) i = 1, 2
We check the validity of the analysis of section on this scalar case by performing a grid convergence study (halving the mesh sizes h K in the reference domain), and by visually checking the preservation of the state η = 1. The computations are run with the RD scalar scheme, but the FV results are almost identical.

The results are summarized in gure 2.7. We can conrm that: when no perturbation is added, the well balanced ALE formulation (2.28) (ALE WB in the gures) preserves the constant state to machine accuracy (not shown in the gures), while the classical ALE form (2.19) (ALE NO WB in the gure)

does not, as the left and middle pictures clearly show.

For the smooth perturbation (and pseudo-bathymetry) considered here we observe second order of accuracy for both the formulations. However the presence of spurious oscillations in the at region increase substantially the absolute value of the error obtained with the unbalanced ALE form.

Accuracy

We consider the advection of a vortex problem proposed in Ricchiuto and

Bollerman [2009] to test the accuracy of the SWEs-ALE schemes. The spatial domain is the square [0, 1]. The vortex is initially centered in (0.5, 0.5) and is transported by a constant eld u = [ 6 0] until it has crossed the whole domain and get back to its initial position at T = 1/6. Periodic boundary conditions are imposed. The computations have been performed on 5 unstructured grids with the topology shown on gure 3. The coarsest has mesh size h = 1/56.

The other 4 meshes have been generated independently, halving the mesh size at each step. The following arbitrary mapping is used to move the mesh

x 1 (t) = χ 1 + 0.1 sin (2πχ 1 ) cos (2πχ 2 ) sin 2πt T x 2 (t) = χ 2 + 0.1 cos (2πχ 1 ) sin (2πχ 2 ) sin 2πt

Mass conservation

We perform the classical test benchmark, taken from [START_REF] Seaid | Non-oscillatory relaxation methods for the shallow water equations in one and two space dimensions[END_REF]. The set-up consists in a square domain [0 × 200] 2 m with a dam, placed at x = 95 m, separating an upper and a lower bassin which contain water at dierent levels, respectively at 10 m and 5 m. To check our mass conservation correction (cf.

section 2.7) we have added a bathymetry shaped as a circular hump centered in (x, y) = (0, 200) [m], and dened by an exponential law in the radial direction (cf. left picture on gure 2.9). We report on the right pictures on gure 2.9 the mass error measured without any correction, and with corrections based on dierent quadrature formulas (for the denition of E mass , see always section 2.7. We can clearly see that we are able to preserve the total mass in the domain practically up to machine accuracy.

Chapter 3

r-adaptation for hydrodynamics

In this chapter, to enhance important features in hydrodynamics simulations such as the resolution of wave patterns or the wetting/drying dynamics, we employ mesh adaptation. We point out the obvious fact that mesh adaptation introduces an additional level of complexity and must be used with great care to be eective. To improve the accuracy of the second order RD/FV-eRK2

Shallow Water solver presented in chapter 1 and 2, we employ mesh adaptation techniques based on nodes redistribution (or relocation) that are also known as r-adaptation techniques. Roughly speaking, these methods move the points of a given reference mesh, keeping the mesh topology and number of mesh points unchanged. In fact in this simple description, it is hidden the main advantage of r-renement with respect to h-renement: the data structure is unchanged and ecient conservative/accurate remaps of ow variables can be carried out from one grid to the updated one.

In the rst part of the chapter we introduce the reader to r-adaptation (a vast review can be found in Budd et al. [2009a]). In this framework we detail the r-adaptive technique of [START_REF] Ceniceros | An ecient dynamically adaptive mesh for potentially singular solutions[END_REF] which is implemented in this thesis. A small variant in the denition of the error estimate is considered to rene the mesh at the wet-dry interface. Once we have understood the tool to move the mesh, in the second part of the chapter we deal with the problem of resolving the SWEs (or general balance laws) on these moving grids. A popular method is the rezoning algorithm of [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF] which is based on conservative ALE remaps. ALE remap is presented in 3.5 and permits to interpolate the numerical solution among the dierent grids. We also recall that an elegant way to resolve the SWEs in an arbitrary reference framework

x = A(χ, t) has been presented in chapter 2. If we choose this arbitrary reference to coincide with the grid transformation, then the solution on the grid

x n , u(x n , t n ), is evolved directly on the new grid x n+1 , u(x n+1 , t n+1 ). Rezoning and ALE algorithms (section 3.6) are tested on scalar problems in section 3.6.4.

From this rst investigation we propose a third coupling algorithm which is a cheaper variant of the rezoning approach. Finally, in section 3.7 we present 3.1. Uniform mesh in metric space the coupled algorithms in terms of accuracy, and CPU time for some classical SWEs benchmarks.

Uniform mesh in metric space

The general idea in mesh generation is to exploit the fact that we can easily generate a uniform mesh in a Riemannian space, which will correspond to a non uniform physical mesh in the Euclidean space, which is used for the computations. With respect to chapter 2, we switch from the term "actual" to the term "physical" in order to highlight that this is the conguration in which our physical conservation law is settle. The physical space is opposed to the computational space, which exists only to compute the mesh. At the end, points in the physical mesh must be clustered where large gradients of the solution appears. With a brief example taken from [START_REF] Alauzet | Metric-based anisotropic mesh adaptation[END_REF], we would like to explain why a Riemannian metric space makes possible to obtain non uniform mesh in the physical space. A manifold M can be seen as a surface embedded in the Euclidean space E 3 and described with Cartesian coordinates:

M = x 1 , x 2 , x 3 (x 1 , x 2 ) (3.1)
with x 1 , x 1 belonging to Ω ⊆ R 2 an open region in the Euclidean plane . First we consider a pair of auxiliary variable to parametrize the surface

χ 1 , χ 2 ∈ Ω χ ⊆ R 2 such that M = x 1 (χ 1 , χ 2 ), x 2 (χ 1 , χ 2 ), x 3 (χ 1 , χ 2 )
We consider a straight curve γ ∈ Ω parameterized by the parameter t ∈ [0, 1] which connects two points P, Q with P (x 1 , x 2 ) and Q(x 1 , x 2 ):

γ(t) = (x 1 (t), x 2 (t))
Its image c(t) on the surface M is a curve of R 3 connecting two points lying on M, R and S. Also the parametric coordinates can be assumed to be function of the parameter t. The curve c(t) simply is:

c(t) = (x 1 (χ 1 (t), χ 2 (t)), x 2 (χ 1 (t), χ 2 (t)), x 3 (χ 1 (t), χ 2 (t))).
The length of the curve c(t) on M, is expressed by (A.14):

l M (P, Q) = 1 0 dχ i dt M ij (t) dχ j dt dt
We assume now that we can construct a mesh such that for every edge P Q, its arc length is constant, let's say unitary l M (P, Q) = 1. This mesh is said to be uniform in the metric space. But its image on Ω corresponds to straight 3. r-adaptation for hydrodynamics edges with non unitary distance. If now we think to the plane Ω as a plane with Euclidean metric, we have constructed a non uniform two dimensional mesh. This illustrates why a Riemannian metric curves the space: a uniform computational mesh corresponds to a physical mesh where points cluster if curvature appears.

Elliptic Moving Mesh PDE

The theory of elliptic grid generation provides a suitable theoretical framework for developing the set of PDEs that controls mesh movement, or moving mesh PDEs (MMPDEs). For clarity we repeat here the main results contained in [START_REF] Thompson | Handbook of Grid Generation[END_REF], and we provide the expression of the MMPDE that we chose among many, as a particular case of this general theory. Later we will see that, at least for the general case, there is an analogy with the theory of harmonic map, [START_REF] Ivanenko | Harmonic mappings[END_REF].

We introduce two (unstructured), simply connected grids: the computational and the physical mesh. We recall from chapter 2 the ALE transformation which transforms the reference mesh into the actual one, see (2.6) and

(2.7):

A : A 0 → A x = A(χ, t), J A = ∂x ∂χ
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Elliptic Moving Mesh PDE

the map should be invertible and we assume J A = det J A > 0 to avoid negative areas. Dierently from what stated in chapter 2, we will assume here that points Q ∈ A 0 and Q ∈ A do not lie in the Euclidean plane but that they both A and A 0 could belongs to a manifold in E 3 . As stated in [START_REF] Warsi | Mathematics of space and surface grid generation[END_REF] and [START_REF] Spekreijse | Elliptic generation systems[END_REF], in the context of mesh adaptation/generation we would like to map a region of manifold called M with local curvilinear coordinates χ into a region of a two dimensional Euclidean plane E 2 described by Cartesian coordinates x. With respect to (2.6) we consider A 0 ⊂ M and A ⊂ E 2 . We put this in evidence in the denition of the ALE map:

A : M → E 2 x = A(χ, t) (3.2)
and we specify the parametric domain to which curvilinear coordinates belongs:

χ ∈ Ω χ ⊂ R 2 and x ∈ Ω ⊂ R 2 see gure 3.1.
As for the rest, we refer to them respectively as computational and physical domain.

In the following, we will extensively make use of the concept of curvilinear coordinates. We give in appendix A a brief review of tensor analysis and a list of the formulas for divergence, gradient and Laplacian operator in generalized coordinates. Concerning the notation we use the one of appendix A but with dierent letters. We assume that M is equipped with metric tensor M of components M ij (

√ M = √ det M ), inverse metric tensor M -1 of components M ij and Christoel symbols Γ k ij .
We start by considering the Laplacian of a scalar φ in the curvilinear coordinates χ, see (A.13)

∇ 2 φ = M ij ∂ 2 φ ∂χ i ∂χ j -Γ k ij ∂φ ∂χ k (3.3)
We write the Laplacian of each physical component, separately. Setting in the above equation φ = x α , we realize that

∇ 2 x α = e i ∂x α ∂x i • e j ∂x α ∂x j = ∂ 2 x α ∂x j ∂x j = 0
and we have an equation for the physical coordinates

M ij ∂ 2 x α ∂χ i ∂χ j -Γ k ij ∂x α ∂χ k = 0 (3.4)
This equation can be further manipulated setting always in (3.3), φ = χ k :

∇ 2 χ k = -M ij Γ k ij ∂χ k ∂χ k = -M ij Γ k ij

r-adaptation for hydrodynamics

And we can write compactly a MMPDE for each physical coordinate

M ij ∂ 2 x α ∂χ i ∂χ j + ∇ 2 χ k ∂x α ∂χ k = 0
This is a set of uncoupled elliptic PDEs which gives back directly the grid in the physical space. This system can be solved numerically.

What is done typically at this point is to generate grids in the computational domain imposing Poisson equations ∇ 2 χ k = P k with control function P k that can be assigned arbitrary. For general

P k = -M ij Γ k ij ,
this is a generaliza- tion of (3.4) and it could still represent a moving mesh generator, although the transformation from a metric space has no more sense and the map x = A(χ, t) is meant A : E 2 → E 2 . We rewrite the MMPDE with control function:

M ij ∂ 2 x α ∂χ i ∂χ j + P k ∂x α ∂χ k = 0 (3.5)
3.3 Moving Mesh PDE from Harmonic Maps [START_REF] Dvinsky | Adaptive grid generation from harmonic maps on riemannian manifolds[END_REF] noted that elliptic meshes discussed in the previous paragraph could be generated through the theory of harmonic maps, as formulated by [START_REF] Eell | Harmonic mappings of riemannian manifolds[END_REF]. We refer to appendix B for a brief digression on the denition of harmonic map and the main statements. Keeping in mind that the map is A : M → E 2 , we gain that one metric is Euclidean H αβ = δ αβ (with δ αβ the Kronecker's delta) and, as a consequence, the energy density for the map (3.2) writes, see B.1:

e = M ij δ αβ ∂x α ∂χ i ∂x β ∂χ j = M ij ∂x α ∂χ i ∂x α ∂χ j
thus the energy functional writes

E(x) = 1 2 M ij ∂x α ∂χ i ∂x α ∂χ j √ M dχ (3.6)
which indeed admits the Euler-Lagrange equations, see (B.2)

M ij ∂ 2 x α ∂χ i ∂χ j -Γ k ij ∂x α ∂χ k = 0
so we have found again (3.4). The advantage of the approach of Dvinsky consists in the fact that it comes together with a prove of existence and uniqueness of the map. The Hamilton-Yau-Shoen theorem (theorem HYS, [START_REF] Hamilton | Harmonic Maps of Manifolds with Boundary[END_REF], [START_REF] Schoen | On univalent harmonic maps between surfaces[END_REF]) states that the map exists when the following two conditions are veried:

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. The rst condition is identically satised by the physical conguration because

Moving Mesh PDE from Harmonic Maps

A ⊂ E 2 . In addition, if the boundary of the physical domain is convex, then the maps exists. We point out that, in hydrodynamics simulations this last hypothesis seems to be a little too restrictive. Let's think for example to the physical domains in gure 3.2. For such cases we have no guarantees that a map into the logical domain exists. This was a reason why, at the beginning, much eort was spent in resolving the inverse map, that is

A -1 : M → E 2 χ = A -1 (x, t)
Following [START_REF] Dvinsky | Adaptive grid generation from harmonic maps on riemannian manifolds[END_REF] and later [START_REF] Huang | Moving mesh strategy based on a gradient ow equation for two-dimensional problems[END_REF], one can choose A 0 to have an Euclidean metric and A to be a region of a Riemannian manifold M. A 0 is at with zero curvature and the convexity hypothesis is no more urgent, since for the computational domain we can always take a convex domain.

Moreover, the HYS theorem guarantees the existence and the invertibility of the map for a general metric. In light of this and for completeness with respect to this small review on moving mesh methods, we report the (inverse) MMPDE of [START_REF] Huang | Moving mesh strategy based on a gradient ow equation for two-dimensional problems[END_REF] (given the complete Euler-Lagrange equations (B.3) you may set Γ α γδ = 0 and invert x α with χ α )

∂ ∂x i M ij √ M ∂χ α ∂x j = 0
which works also for non convex physical domain. Of course the drawback of taking the mapping in the opposite direction is that we need to transform it, somehow, back into the curvilinear coordinates x and this makes the resulting equation much more complex.

3. r-adaptation for hydrodynamics

MMPDE of Ceniceros and Hou

In the following we will stick to the direct map x = A(χ, t) and to the general MMPDE (3.5) that the interested reader will nd in [START_REF] Warsi | Mathematics of space and surface grid generation[END_REF]. Consider the following diagonal metric

M ij = 1 ω δ ij with √ M = ω -1 . The function ω is called monitor function. The inverse metric is M ij = ωδ ij .
We choose now the control function

P k = ∂ω ∂χ k
and substitute in the MMPDE (3.5)

δ ij ω ∂ 2 x α ∂χ i ∂χ j + ∂ω ∂χ k ∂x α ∂χ k = 0 ω ∂ 2 x α ∂χ i ∂χ i + ∂ω ∂χ k ∂x α ∂χ k = 0
We end up with the simple MMPDE rst obtained by [START_REF] Ceniceros | An ecient dynamically adaptive mesh for potentially singular solutions[END_REF]:

∂ ∂χ i ω ∂x α ∂χ i = 0 (3.7)
This mesh generator has many advantage, among them 1. it is a set of decoupled quasi-linear elliptic equations for which many ecient numerical methods are available 2. we obtain directly the physical coordinates of the mesh points.

Among the drawbacks there is fact that folded grids can occur. Another issue is the control of the mesh quality, for example the element's skewness and the mesh alignment with respect to the ow direction. With the MMPDE (3.7)

based on diagonal metric ω -1 δ ij , we have little control on these quantities. I will briey detail this issue later. Originally, [START_REF] Ceniceros | An ecient dynamically adaptive mesh for potentially singular solutions[END_REF] have tested an implementation of (3.7) to improve the resolution of blow-up problems and heat convection. The results were promising since the method prove to be computationally ecient and capable to follow complex ow evolutions up to small scale phenomena. This approach has been quite successful, and it has been used among the others, in [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF], Chen et al. [2008] and Zhou et al. [2013a].

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

MMPDE of Ceniceros and Hou

Ceniceros [START_REF] Ceniceros | An ecient dynamically adaptive mesh for potentially singular solutions[END_REF] derived and presented equation (3.7) as the Euler-Lagrange equation related to the following functional:

E(x) = 1 2 Ωχ ω ∂x α ∂χ i ∂x α ∂χ i dχ (3.8)
It is worth to note that for this particular mesh generator with P k = -M ij Γ k ij , we loose the analogy with harmonics maps, as we can see by checking the dierence with respect to the functional

(3.6) → M ij √ M = ωδ ij .
Finally we think it is interesting to write (3.7) with a sort of elastic analogy, replacing the position with the displacement respect to the computational mesh

x α = X α + δ α : ∂ ∂χ i ω ∂δ α ∂χ i = ∂ω ∂χ α (3.9)
where Σ αi = ωδ ij ∂δ α ∂χ j is a pseudo elastic stress tensor, and the right hand side act as force on the medium F α = ∂ω ∂χ α . We remark that the monitor function play a role in controlling both the stiness of the system and the force.

Boundary conditions

When solving problem (3.7) we will assume that the computational domain is a closed polygon whose boundary ∂Ω χ is composed by the union of m segments. ∂Ω χ is mapped into the boundary ∂Ω and we further assume that it is invariant to the transformation. For instance we will have no moving boundaries. In particular we consider free-slip boundary conditions δ • n = 0, χ ∈ ∂Ω χ (3.10) with δ = 0 at the polygon's vertexes. A standard method to impose boundary condition is contained in [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF] where it is introduced a second map A ∂ : ∂A 0 → ∂A which correspond to the trace of (3.2) on the boundary. This mapping is then used as Dirichlet conditions to solve the transformation for inner points. Alternatively as shown in [START_REF] Li | A moving mesh nite element algorithm for singular problems in two and three space dimensions[END_REF] the variational formulation could be complemented by a constraint equation to take into account (3.10). We will however stick to form (3.9), written in terms of displacement, which is suited to express directly the boundary conditions.

Monitor function

In general the metric tensor is related to the Jacobian of the transformation which we called J A , M = J T A J A (see appendix A). The Jacobian controls size and orientation of mesh elements. Following [START_REF] Huang | Mathematical principles of anisotropic mesh adaptation[END_REF] we write J A using the singular value decomposition (SVD):

J A = U ΛV T (3.11)

r-adaptation for hydrodynamics

where U and V are the orthogonal matrices associated with left and right singular vectors and Λ is the diagonal matrix consisting of the singular values. This decomposition has a geometrical meaning: U is a rotation matrix that species the orientation of physical mesh elements and Λ species the size and shape of these elements. Setting properly the SVD one can construct monitor matrix to achieve certain properties of the elements, such as anisotropy and skewness, see [START_REF] Alauzet | Metric-based anisotropic mesh adaptation[END_REF]. What we want to point out is the dierence with the diagonal metric ω -1 δ ij that decouples completely the two directions. That is, x 1 and x 2 are resolved independently with the scalar monitor function ω.

While this allows to achieve a quite good anisotropy, it does not guarantee the rened control of the decomposition (3.11). In practice this means that we will have limited control on mesh quality. Now we go back to the denition of ω. A classical denition for scalar problems, given by [START_REF] Winslow | Numerical solution of the quasi-linear poisson equation[END_REF], couples the mesh motion with the gradient of the solution of the underlying PDE on the physical mesh: ω = ω(∇u). As in Zhou et al. [2013b], we have selected the free surface η, in order to detect free surface wave patterns and bore development. We propose the following denition of the monitor function

ω = 1 + α (max (||∇η|| * , ||∇ 2 η|| * )) 2 + δφ 2 (3.12)
We see that when ||∇η|| * , ||∇ 2 η|| * , φ → 0, the metric becomes the Euclidean one M ij = δ ij and the mesh tends to the unperturbed state. • represent normalized L 2 -norms computed as

||∇η|| * = min 1, ||∇η|| β max||∇η|| , ||∇ 2 η|| * = min 1, ||∇ 2 η|| γ max||∇ 2 η||
The coecients α, β and γ are free parameters, allowing to optimize the mesh movement. Note that in all of the above formulas, the derivatives of η are computed on the physical (moving) mesh, making problem (3.7) nonlinear.

In (3.12) we have also tested the inuence of the Hessian of the free surface and a tracking of the wet/dry interface. There are in literature some examples of such front-tracking error functions. For example, in the context of phase change problems, J.A. Mackenzie and W.R.Mekwi [2007] dened ω = α/ β|x -x interf | + γ. This expression, like others, requires the knowledge of the distance function from the interface, whose computation may be quite costly. Here we propose a simpler approach explicitly exploiting the knowledge that h → 0 at the front. We have added a new term φ = max(∇f (x), ∇ 2 f (x)) where

C H < h < C U :      f (x) = 0, if h(x) ≤ C H f (x) = h-C H C U -C H , if C H < h(x) < C U f (x) = 1, if h(x) ≥ C U
In gure 3.3 we report an example to show how this technique may work. The mesh is rened over the wet/dry cells and, due to the presence of the Hessian norm in (3.13), it is also rened one cell after the interface in dry region. We have added the Hessian contribution after having observed that a more spread mesh renement at the interface, ensures robust ooding simulation.

Mesh Smoothing

The smoothness of the mesh is measured in terms of the variation of the local element's size h K over the domain Ω. More rigorously, given the map (2.6), a mesh has degree of regularity r if the Jacobian J ∈ C r (Ω) , see Budd et al.

[2009a]. Mesh regularity is important because it is widely recognized that abrupt variation of the mesh size leads to a deterioration of the numerical solution. Typically, meshes obtained by (3.7) are smooth, but, in presence of non convex boundaries ∂Ω loss of regularity can occur, ultimately leading to mesh tangling in extreme cases. [START_REF] Huang | Analysis of moving mesh partial dierential equations with spatial smoothing[END_REF] studied the property of the following one dimensional equation:

1 -ν ∂ 2 ∂X 2 x = x
and demonstrated that x satises the following smoothness condition:

|J| = ∂x ∂X x -1 < ν -2

r-adaptation for hydrodynamics

The two dimensional extension can be found in Budd et al. [2009a] and reads:

1 -ν ∂ 2 ∂X i ∂X i x α sm = x α (3.14)
This MMPDE is particularly ecient to resolve because the operator

L = 1 -ν ∂ 2 ∂X i ∂X i
can be inverted once on the computational mesh, and then used to compute a smoother variation of a given mesh (typically the solution of (3.7)) as:

x α sm = L -1 x α
To avoid mesh tangling at the corner of non convex domain, we have dened a variable diusion coecient based on the distance form the boundaries, d P,∂ = d(x P , ∂). Since the boundaries do not move, the distance is computed once and then successively interpolated.

ν i = max 0, 100 - d i,∂Ω h k 2
In gure 3.4 we can see a case where smoothing is not applied and the mesh tangles at the corner, while, applying a sequence of the two MMPDE (3.7) and

(3.14) at each time step, we observe a smoother mesh without tangling. We remark that the above denition of the diusion parameter is heuristic and there is still no guarantee that tangling never occur.

Numerical resolution of MMPDE

In practice, given an initial mesh in the computational domain, the weak form of (3.7) with boundary conditions (3.10) is discretized with a standard P 1 Galerkin nite element method. We search for an approximate solution

x h ∈ V h such that Ωχ ∂x α h ∂χ m ω ∂v h ∂χ m dχ = 0 ∀v h ∈ V h , α = 1, 2
If we use continuous piecewise polynomials, the FE solution

x α h = j∈T h ϕ j x α j {ϕ i (x 1 , x 2 )} N
1 is the Lagrangian basis and x α j coincides with the position of the node labeled by j. The standard development is:

j∈T h Kχ ∂ϕ j ∂χ m ω ∂ϕ i ∂χ m dχ x α j = 0 ∀i ∈ T h , α = 1, 2
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

MMPDE of Ceniceros and Hou

∂ϕ j ∂χ m = 1 2|K χ | n m j
We can write the algebraic system in terms of displacement δ α = x α -χ α .

j∈T h a ij δ α j = b α j ∀i ∈ T h , α = 1, 2 (3.15) with: b α j = j∈T h a ij χ α j a ij = K∈D i ∪D j Kχ ω dχ 1 4|K| 2 n m i n m j
Due to the dependence of ω on the derivatives of η on the new mesh, the weak form (3.15) denes a nonlinear system of algebraic equations which needs to be solved by means of some iterative procedure.

The choice of this procedure and its coupling with the ow evolution equations plays a crucial role in determining the balance between the gain brought by the adaptation procedure, and its cost overhead with respect to the evolution of the ow quantities with the explicit schemes discussed in chapter 2. In this regard it is worth noting that the eRK2-RD method is ten times faster compared to the implicit version, [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF]. For all these reasons, we have chosen a simple explicit Newton-Jacobi iteration method, as in Chen 3. r-adaptation for hydrodynamics et al. [2008]. In particular, at each time step, the displacement δ k = x kx n is computed from the following relaxed iteration

δk+1 i = δ k i - 1 a ii j∈D i a ij x k j (3.16) x k+1 i = x n i + µ i δk i (3.17)
Note that the update (3.16) is equal to the one of Chen et al. [2008], but recast in terms of displacements so to embed more naturally the boundary conditions.

As in the last references, to improve the control on the regularity of the mesh, we have introduced a relaxation phase in the iterations. In particular, the following denition of the relaxation parameter µ i has been used (cf. also [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF], Chen et al. [2008])

µ i = min (1, max (ϑ, τ ||∇η i ||))
To avoid nodes'depletion in regions with small solution variations, a threshold for the stiness is tuned by xing ϑ, if ϑ ∼ 0 the stiness in regions where ∇η ∼ 0 is strongly increased. Finally, we recall that the entries of the stiness matrix a ij depend on the value of the monitor ω, and thus on the value of the solution on the new grid. As a consequence an essential element of this method is a suciently accurate projection step allowing to remap the discrete solution on the moving mesh. This projection step has to be chosen very carefully, as it impacts the overall accuracy, monotonicity, and cost of the computation. This issue will be extensively covered in section 3.5.

Mesh Generation

We test the moving mesh algorithm dened by (3.16) on the benchmarks proposed in [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF]. The monitor function is computed according

to ω = √ 1 + αu 2 , with u assigned: u(x, y) = exp(-8(x 2 + 9y 2 -1) 2 ) (3.18) u(x, y) = exp(-100(y -x 2 + 0.5) 2 ) (3.19) u(x, y) = 50 exp(-2500(x 2 + y 2 ) (3.20) u(x, y) = 1 if y = x 0 if y = x (3.21) with x = x 1 and y = x 2 . The reference domain is a square [-1, 1] × [-1, 1].
The iteration is repeated in the pseudo-time loop until convergence is reached.

For the rst two smooth examples reported in gure 3.5 we have set α = 100.

In the left picture of gure 3.6 we show that also singularities are well handled

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

3.5. High order projections from ALE remaps with a proper choice of the parameter. Finally always in gure 3.6 on the right, there is a test for an oblique shock: the zoom shows that anisotropy in the direction of the shock is achieved without tangling occurrence.

High order projections from ALE remaps

As already said, the Newton-Jacobi iterations (3.16) and (3.17) require the projection of the solution values on the last updated mesh. The problem have already been formalized in section 2.3.1. We want to evaluate scalar and vector properties of the uid through the ALE transformation. For the vector of ow variables:

w = w • A = w 0 (χ)
3. r-adaptation for hydrodynamics Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

High order projections from ALE remaps

We recall here the remap equation (2.15):

∂J A w ∂t χ -J A ∂wσ j ∂x j = 0 (3.22)
and we remark that the interpolation has been formalized into a typical hyperbolic problem in which the solution is advected by the domain movement.

It is not a coincidence that it is usually referred to as advection remap.

The problem of computing updates for the solution values due only to the mesh movement can be elegantly solved by using remaps generated by the same schemes used to discretize the PDEs. Such a remap corresponds to the limit for ∆t → 0 of the schemes presented in sections 2.10 and 2.11 with a "frozen" the uid conguration B(t * ), but letting the mesh conguration A(t) continue to move. This provides an instantaneous approximation to the conservative ALE remap equation.

Practically, during the iterations (3.16) we generate a sequence of sub-grids

T k+1
h on which we want to interpolate the last available numerical solution. If we are at the k-th iteration, the last available solution is w k i = w(x k i ).

FV ALE remap

For the FV scheme, taking the limit for ∆t → 0 of (2.37) we obtain the one step projection over the sub-grid T k+1 h

|C k+1 i |w k+1 i = |C k i |w k i - j∈D i R ij (w k ) (3.23) and R ij (w) = j∈D i -∆x ij wj + wi 2 - |∆x ij 2 (w j -wi )
the interface velocity is replaced by an interface displacement (which coincides with the area swept by the same interface):

∆x ij = ∂C ij ∆x m n m ds ∆x m = x m,n+1 -x m,n
is the displacement of the interface during a single time step. The advantage of this approach is that it retains all the properties of the original method. A second order, non-oscillatory, well-balanced, mass conserving projection can be obtained by applying the limited high-resolution FV scheme, referred to as MUSCL in section 1.6.2. If the scalar, decoupled nature of the projection equations (all quantities independently are transported in the direction of the displacement) reduces the cost of these evaluations, it still means that the cost of one projection will be that of a single step of the FV scheme. As this may be repeated at every Newton-Jacobi iteration, this cost may lead to an important overhead.

r-adaptation for hydrodynamics RD ALE remap

In the ALE remap the RD uctuation denes

Φ K (w h ) = - K ∆x j ∂w ∂x j dx
We then take the limit ∆t → 0 in the RD predictor step, see (2.55):

w * i = w n i - 1 |C n+1 i | K∈D i Φ K(1) i (w k h ) with Φ K(1) i = Φ K (w k h ).
For the corrector, see (2.56):

w k+1 i = w * i - 1 |C n+1 i | K∈D i Φ K(2) i (w * h , w k h ) (3.24) with Φ K(2) = 1 2 Φ K (w k h ) + Φ K (w * h ) + K w * h -w k h ∆t dx
We remark a dierence with respect to FV. To obtain a second order remap the two steps projection must be carried, making the RD remap less ecient then the FV one. We believe that this is related to the presence of the mass matrix.

With such a two steps projection, the non linear splitting Φ K i = Φ LLxF -SU P G i of section 1.6.3 allows to project the ow variables retaining all the properties of the RD LLxF-SUPG scheme (second order, non oscillatory solutions, wellbalanced and mass conservation).

Adaptive algorithms

The coupling of the ow solver with the mesh at each time step is non-trivial, as the mesh equations depend on the solution on the (unknown) adapted mesh.

In particular the Shallow Water equations and the MMPDE can be either solved simultaneously or alternately. The latter algorithm, generally speaking, alternate at every time step the solution of the MMPDE and of the underlying PDEs as follows:

• compute the monitor function based on the current solution

• evolve the MMPDE and compute the new mesh

• evolve the PDEs and compute the new solution This has been successfully implemented by [START_REF] Huang | Moving mesh strategy based on a gradient ow equation for two-dimensional problems[END_REF] showing a signicant reduction of stiness problems with respect to the full-coupled Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

3.6. Adaptive algorithms approach. At the same time is it has been recognized that it can lead to a lag in the mesh movement with respect to the physical features.

In this thesis, depending on the framework in which we evolve the PDEs, two dierent alternate algorithms are tested. If the PDEs are written in Eulerian framework one gets the rezoning method suggested in [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF]. This approach is based on a sequence of mesh and ow iterations.

It uses the mesh solver as a black box, the ow equations being solved on a (dierent) xed mesh at each time iteration. Its drawback is that, at each time iteration, the ow solver requires a remap/interpolation on the new mesh which may be quite expensive as it needs to guarantee the same properties as the ow solver itself (high order accuracy, non-oscillatory character/positivity preservation, C-property, mass conservation). At the opposite, once the grid has been adapted, one can evolve the ow with an ALE formulation of the PDEs. In this case, the properties of the ow solutions are only determined by the scheme. In the moving mesh community this approach is called quasi Lagrangian, see [START_REF] Huang | A high dimensional moving mesh strategy[END_REF] and [START_REF] Cao | An r-adaptive nite element method based upon moving mesh pdes[END_REF]: time derivatives are transformed along mesh trajectories and the resulting PDEs are a non conservative form of our balance laws (2.18). Given a PDE in Eulerian form ∂u ∂t + L(u, x, t) = 0 it can be transformed using the relationship between Eulerian and ALE time derivative (2.10)

∂u ∂x χ - ∂x j ∂t χ ∂u ∂x j + L(u, x, t) = 0
We note that this approach hides the GCL. The consequences is that, approximation of the above ALE equations can represent exactly a constant uniform ow on moving meshes without the need to conserve the total mesh volume along the simulation; this latter statement in fact is not an issue in Huang

and Russell [1998]. A combination of the alternate algorithm and the ALE approach in conservation form can be found in Ni et al. [2009].

We have now all the basic blocks to perform adaptive mesh simulations.

These boil down to the ow evolution equations (section 2.10 and 2.11) and to the MMPDE, discussed in section 3.2. We propose hereafter 2 alternate techniques, which are extensively tested in the numerical results. A weakly coupled ALE method and a decoupled adaptation-evolution steps. Particular cases of these two implementations have already been considered in literature (see e.g. [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF] and [START_REF] Ni | Remapping-free ale-type kinetic method for ow computations[END_REF] for the ALE). Here we will compare their impact on the overall cost of the simulation, and on the quality of the results. The set of nodes' coordinates is called x = {x 1 , ..., x i , ..., x N }

with N the number of mesh nodes. We used only kmax = 5 iterations of the

r-adaptation for hydrodynamics

Newton-Jacobi method which, it is important to remark, do not ensure the convergence of the iterative method within each time step. In practice they are sucient to get a good node renement.

Moving Mesh ALE algorithm (ALE)

In this approach the balance law is written using the ALE formulation directly in a framework coincident with the moving domain. At every time step we get the solution on the adapted grid, independently on the interpolation scheme which is only needed now to evaluate the error monitor. The algorithms reads :

Step 1. Taken a triangular mesh T n h , compute the vectors of nodal coordi- nates x n , and the initial solution w n h . Set the initial conditions for the MMPDE, η 1 h = η n h and x 1 = x n .

DO k=1,kmax=5

Step 2. Compute the monitor function ω k = ω(η k h ) and, in turn, the moving mesh matrix a ij = a ij ω k . Move the mesh according to the Newton-Jacobi iteration (Eq. (3.16) and (3.17)). At each iteration we get x k+1 .

Step 3. Compute the interpolated free surface η k+1 Step 5. Let T n h = T n+1 h and w n h = w n+1 h .

IF (t > T) EXIT ELSE GO TO Step 1.

We see that the interpolated solution is only used to evaluate the error function. As a consequence, we thought that the interpolation step can be simplied a great deal without aecting the quality of the solution, as the numerical tests will conrm.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. Bottom) Rezoning (EUL2).

Adaptive algorithms

3. r-adaptation for hydrodynamics 3.6.2 Moving Mesh Eulerian algorithm/rezoning (EUL1)

In this case, the balance law is resolved numerically at every time step in a purely Eulerian framework, and on a xed mesh. The latter is then adapted to the new solution and an accurate guess for the values of the last solution on the new mesh is provided by the projection scheme. The algorithm reads:

Step 1. Taken a triangular mesh T n h , compute the vectors of nodal coordi- nates x n , and the initial solution w n h . Set the initial conditions for the MMPDE, w 1 h = w n h and x 1 = x n .

DO k=1,kmax=5

Step 2. Compute the monitor function ω k = ω(η k h ) and matrix a ij = a ij ω k .

Move the mesh according to the Newton-Jacobi iteration (Eq. (3.16) and

(3.17)). At each iteration we get x k+1 .

Step 3. Compute the full interpolated solution w k+1 Step 5. Let T n h = T n+1 h and w n h = w n+1 h .

IF (t > T) EXIT ELSE GO TO Step 1.

Since this time the interpolated solution will act as the initial condition for the new time iteration, great care has to be put in its computation. The interpolation step does not have to spoil the accuracy property of the numerical scheme. As a consequence, costly projections obtained from high resolution non-linear schemes have to be used to ensure that the quality of the results is not spoiled. Here ALE remaps based on high resolution schemes are used to this purpose.

Modied Moving Mesh Eulerian algorithm (EUL2)

In the previous algorithm, a double role emerges for the interpolation step.

Firstly we need an interpolated solution η k h at every Newton sub-step in order to evolve the mesh. Secondly we provide an interpolated solution w kmax+1

|C n+1 i |w kmax+1 i = |C n i |w n i - j∈D i R ij (w n i )

Adaptive algorithms: eciency

To test the eciency of the dierent coupling algorithms proposed in section 3.6 (ALE, EUL1, EUL2), we propose two scalar tests: linear rotation and Burgers'equation. We would like to conrm the feeling that, for ALE and EUL2, the numerical resolution of the remap equation into the MMPDE could be quite inaccurate (e.g. rst order) without destroying mesh quality. Can this result leads to more ecient algorithms with respect to the classical rezoning/EUL1?

The following MMPDE parameters are used in all the tests α = 10, β = γ = 0.15. We did not perform a systematic optimization relative to these parameters. The value used correspond to those that visually provided the best mesh.

Rotation

This is the classical rotation of a smooth sinusoidal hill, but with a source term. We recall our general scalar balance law (2.65):

∂u ∂t + ∂F j ∂x j + S(u, x) = 0 , x ∈ [-1, 1] × [-1, 1], t ∈ [0, π] (3.25)
where, for a given ux vector F (u) = [ F 1 F 2 ], the source term is dened as in section (2.12.1) S = a j ∂b ∂x j this time with b(x, y) = 0.8e ψ(x,y) and ψ = -5y 2 -5x 2 . The following denition of the uxes will be used, F (u) = au, with a = [ -2y 2x ]. The initial condition for the "free surface" error on the CPU time. We perform the same test for both the RD and FV scheme. In the left picture of gure 3.10 and 3.9 the convergence curves for the dierent combinations of moving mesh algorithms and interpolations schemes are reported. For the ALE algorithm and EUL2 (not shown in gure), we see that all the curves in blue color, corresponding to dierent remap schemes, are almost overlapped (dark blue is rst order LxF/upwind FV, electric blue is a linear centered approximation, light blue is a second order scheme SUPG/FROMM). On the contrary for the EUL1 algorithm there is only one interpolation scheme which guarantees stable and second order accurate results, actually the one which we evolve the PDE with, namely the SUPG/FROMM scheme. We can summarize saying that, even if the interpolation of the monitor function has a positive impact on the resolution of the MMPDE and reduces mesh delay, the specic scheme used, weakly inuences mesh conguration.

η = u + b is η 0 (x) = 1 + cos 2 (2πr) if r ≤ 0.25, r = x 2 + (y -0.5) 2 η 0 (x) = 1
In the right picture of gure 3.10 and 3.9 the performances of the dierent algorithms are compared in terms of error/time. For the RD method, the ALE algorithm shows the lowest CPU time for a xed error level (roughly 80% faster then a xed grid computation). The Eulerian algorithms are less ecient because the full two stage RK interpolation had to be implemented (60% gain for EUL2 and 35% for EUL1). For the FV scheme the eciency of the ALE and Eulerian algorithms is more similar (ALE and EUL2 80%, EUL1 70%). The reason is that, in this case, the second stage of the interpolation is not necessary, [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF]. Finally, for both RD and FV, the EUL2 represents a slight improvement respect to the EUL1 algorithm.

Burgers'equation

In this section we test test if the above conclusions are true when discontinuities develop. Solutions with discontinuities are obtained with a Burgers equation and discontinuous initial conditions. We set in (2.65) S = 0,

F (u) = [ u 2 2 u 2 2 ].
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. The initial condition is

Adaptive algorithms

u 0 = 1 if x ∈ [-0.6, -0.1] × [-0.35, 0
.15] u 0 = 0 otherwise A reference solution is computed on a xed unstructured mesh with h K = 1/100, which is referred to as the ne mesh. To test the eectiveness of the adaptation algorithm, computations are performed on a xed structured mesh with element reference size h K = 1/40, which is referred to as coarse. A good mesh renement is obtained in correspondence of the discontinuity. The comparison between EUL1, EUL2 and ALE strategies is made for both the scalar RD and FV.

From gure 3.13 to gure 3.20, the results for RD and FV are shown. This time we have found that the interpolation step slightly aect the mesh adaptation. We have reported the extreme cases where centered approximation and rst order remap scheme are used. We can observe that the rst order remaps produces slightly smoother mesh at discontinuities; of course this is as a consequence of the fact that we compute a smoother monitor function. The reader may refer to [START_REF] Huang | Analysis of moving mesh partial dierential equations with spatial smoothing[END_REF] for the consequence of a proper smoothing of the monitor function. For the ALE, the CPU time remains more then two times smaller respect the one obtained with the ne grid, while the two solutions are downright comparable. The EUL2 algorithm is obtained using a Galerkin remap into the MMPDE. The advantage respect to EUL1, in term of CPU time, is clear from table (3.1). In particular for RD, the full two stages eRK2 remap, makes for this test, the EUL1 not ecient.

Using these results we highlight that:

• in the following numerical test cases, the EUL2 and ALE algorithms will be used in their faster versions with inaccurate centered approximation remap into the MMPDE.

r-adaptation for Shallow Water ows

In the nal section of the chapter we present a thorough study of the coupling algorithms in terms of accuracy, and CPU time for both simple Shallow Water academic problems and some classical benchmarks involving the long wave Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. runup on complex bathymetries. In SWEs simulations, the MMPDE parameters used in the monitor function (see section 3.2) are α = 20, β = γ = 0.10 and δ = 3α, unless otherwise specied. For the relaxation parameters we used τ = 3 and ϑ = 0.7.

r-adaptation for Shallow Water ows

Asymmetric dam Break

This classical test benchmark, taken from [START_REF] Seaid | Non-oscillatory relaxation methods for the shallow water equations in one and two space dimensions[END_REF], is used to test the adaptive algorithm when bores develop. The set-up consists in a square domain [0 × 200] 2 m with a dam, placed at x = 95 m, separating an upper and a lower basin which contain water at dierent levels, respectively at 10 m and 5 m.

The sudden break of the dam leads to a depression wave advancing in the upper basin and a bore advancing in the lower basin. Two corners depression interact, forming a deep trough at the inlet of the dam.

The test is run with both the FV and RD scheme, on a coarse triangulation containing 14538 triangles and 7480 nodes, on a ne one, containing 77302 triangles and 39130 nodes, and on the coarse mesh with adaptive mesh deformation. The typical qualitative result obtained is provided in gures 3.21 and 3.22. The pictures show the potential of this adaptation procedure to provide with considerably fewer unknowns a better resolution of the breaking bore.

In gures 3.23,3.24 a comparison between the ALE algorithm and the EUL1 and EUL2 is shown. For both RD and FV, the ALE algorithm shows a well resolved bore and a correct computation of the trough with a signicant saving in CPU time. As shown on table 3.2, the savings obtained with the ALE algorithm go up to 60% for RD, and 50% for FV. For the RD scheme, the cost of a two-step interpolation, makes the EUL1 algorithm inecient, thus the EUL2 is a clear improvement. For FV both the interpolation based algorithms (EUL1 and EUL2) are not able of providing a considerable improvement in the resolution of the peaks and the trough upstream the dam (x 60 [m]), probably due to excessive numerical diusion in the interpolation. Some improvement is instead observed with the ALE algorithm, which also gives a much sharper capturing of the bore.

A second test is performed to provide qualitatively evidence that r-adaptation could handle complex ows. It consist in the simultaneous break of two sym-3. r-adaptation for hydrodynamics We observe the mesh rening to capture sharp bore interaction and vortex formation.

Small perturbation of a lake at rest

We consider the classical test of a small perturbation over an elliptic exponential hump (se e.g. [START_REF] Seaid | Non-oscillatory relaxation methods for the shallow water equations in one and two space dimensions[END_REF]; [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF] for details concerning the test setup). This test allows to check the ability of the algorithms proposed to catch relatively smooth wave patterns, and to conserve mass, and the lake at rest state in the unperturbed regions. To run the test, we use a coarse triangulation, containing 12142 nodes and 23852 triangles, and we compute reference solutions on a ner mesh, containing 50631 nodes and 100376 triangles.

The qualitative behavior of the methods proposed can be seen in gures 3.26 and 3.27 (same contour lines drawn in all the pictures). We can see that the mesh follows quite well the propagation and transformation of the waves, providing, on the coarse mesh, a resolution very close to the reference one. No numerical artifacts are observed in the unperturbed region, as a consequence of the exact preservation of the lake at rest state. To perform a more quantitative analysis we report in table 3.3 the CPU times of all the schemes, and the water height along the line at y=0.5 on gures 3.28 and 3.29. For clarity, only the 3. r-adaptation for hydrodynamics EUL2 method results are plotted in the latter gures, the EUL1 algorithm providing virtually identical solutions.

The cuts show how both the ALE and the rezoning algorithms provide solutions close to the reference one. The CPU time savings w.r.t. the reference are of the order of 70% for the ALE method, of 60% for the EUL2, and between 50% (for FV) and 40% (for RD) for the EUL1 algorithm.

Finally, gure 3.30 shows a study of mass conservation, providing additional proof that the corrections proposed allows to retain the physical mass in the domain virtually to machine accuracy. 3. r-adaptation for hydrodynamics

r-adaptation for hydrodynamics

Runup on a conical island

This is a very classical benchmark for tsunami simulation models. It aims at reproducing some of the experiments of [START_REF] Briggs | Laboratory experiments of tsunami runup on a circular island[END_REF]. We refer to the above reference, and to [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF]; [START_REF] Delis | A robust highresolution nite volume scheme for the simulation of long waves over complex domains[END_REF] for the test setup. The parameter in the MMPDE are α = δ = 20, β = γ = 0.2. This benchmark will allow to test the ability of the algorithms proposed to track dry fronts, as well as the mass conservation correction. We have run the test on two meshes, both progressively rened in the region of interaction between the wave and the conical island. The coarse one, contains 10401 nodes, and 20580 triangles, with mesh sizes h K going from 0.5 to 0.2 meters. The ne mesh contains 37982 nodes, and 75594 triangles, with mesh sizes going from 0.3 to 0.08 meters. The ne mesh results obtained with FV and RD are quite close (cf. gure 3.32), and similar to those typically shown in literature. They have been used as a reference for those obtained on the coarse mesh, with adaptive mesh deformation.

The qualitative behavior of the method is shown on gure 3.33. The pictures show the ability of the modied monitor function to track both the incoming and refracting waves, and the moving wet/dry interfaces. The gauge signals for the adaptive simulations are reported in gure 3.34 for the gauges g9 (upstream the island), g16 (lateral runup), and g22 (rear side runup). The results obtained on gauges 9 and 22 show that, for both FV and RD, the adaptive ALE algorithm provides results comparable to those obtained on the ne mesh. In particular, the interference between the two refracted waves that causes the peak and highest runup values on the back of the island, is well reproduced. This is also the case with the interpolation-based methods, which provide practically the same results (only EUL2 show in the plots). In the RD case, all the adaptive algorithms lead to a less impressive improvement in the lateral runup gauge 16.

CPU times are reported on table 3.4. We can see that the ALE adaptive computations allow still savings of the order of 71% w.r.t. the ne mesh computation. The percentages of CPU time reduction for the rezoning algorithms are close to 66% for the EUL2 method, and to 37% (for RD) and 44% (for FV)

for the EUL1 algorithm. Lastly, the tables also report the % of the total cost represented by the moving mesh algorithm alone including the re-computation of geometrical quantities. These show that, while for the ALE the overhead w.r.t. a xed mesh simulation is of 40%, the EUL2 and EUL1 algorithms counts for, respectively, 50% and 70% of the computation. This means that more time is spent adapting the mesh than in computing the ow. Clearly, this is a consequence of the costly projection steps on which the method relies.

Finally, gure 3.35 shows the study of mass conservation for this problem.

The pictures prove how a high accuracy correction of the nodal bathymetric heights, combined with the redistribution of the spurious geometric mass gen-Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. erated by the motion of dry nodes, allows to reduce the mass error practically to zero.

r-adaptation for Shallow Water ows

Monai valley benchmark

This test involves the tsunami runup over a complex 3D bathymetry, and is a standard test for tsunami simulation models Liu et al. [2008]. The experiment that it reproduces was carried out at Central Research Institute for Electric Power Industry (CRIEPI) in Abiko (Japan), and consisted of a 1/400 reproduction of the Hokkaido-Nansei-Oki tsunami of 1993 that struck Okushiri Island, with disastrous consequences especially in the region of the Monai village, on which the experiment itself focuses. For a full description of the setup, including all the necessary data to run the test, and with the results from the experiments, we refer to page of the center for tsunami research at NOAA [1993]. We have run this test on the grids reported on the right pictures of gure 3.36, statically adapted to the bathymetric variations [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF].

The coarse one contains 7000 nodes and 13720 triangles, with mesh sizes h K ranging from 0.1 to 0.025 meters; the ne mesh contains 36911 nodes and 18711 triangles, with sizes ranging from 0.05 to 0.01 meters. Note that the prescribed uniform mesh size for this test is usually of 1.4 cm NOAA [1993].

We have used the ne mesh results as a reference, to compare against the solutions obtained with adaptive mesh deformation on the coarse grid.

The qualitative impact of the adaptation algorithms has been visualized on gure 3.36, which reports plots relative to the instant of maximum runup. The top rows report the x grid results, while the bottom one shows the solution on the adaptive grid, and the mesh itself. The moving adaptive result shows a clear improvement in the reected bores, and, as we will see in more detail shortly, runup heights very close to those obtained on the ne mesh. Note that this is a dicult test for the r-adaptation, as the initial non-uniform mesh size distribution leads to strongly anisotropic triangles in the adaptive case, as clearly visible in the gure.

As already remarked in [START_REF] Ricchiuto | An explicit residual based approach for shallow water ows[END_REF], there is little inuence of the mesh size on the gauge signals. This is shown clearly by the water height signal in gauge 7, reported for completeness in gure 3.37. A much more interesting quantity to look at is the runup plot, which is provided in the top row of gure 3.38. In the pictures, the brown line represents the height of the maximum runup observed in the experimental setup in the narrow gulley with a cove at (x, y) ≈ (5.15, 1.875) [m] in the scaled down model. The gure shows that only with ner grids the correct runup height can be reached, and that both the ALE and rezoning methods allow to obtain the correct prediction on the coarser grid. To corroborate this result, we have placed an additional gauge (not present in the experiment). Its position is at (x g , y g ) = (5.05, 1.9)[m], very close to the maximum runup point. The water height time series in this gauge are reported in the bottom row of gure 3.38. These pictures conrm that the ALE algorithm is superior in allowing to retain the correct values of the maximum water heights, even though failing in reproducing the exact shape of the signal. The rezoning methods also provide a considerable improvement over the coarse mesh result, with water heights very close to the reference.

CPU times are given in table 3.5.

Lastly, the evolution of the mass conservation error is reported on gure 3.39. Again we can see the improvement brought by the corrections proposed here.

Solitary wave on a shelf with an island

Finally, as an application to a more complex ow, we consider a laboratory experiment, conducted in the wave tank of the Oregon State University, involv-Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. 3. r-adaptation for hydrodynamics Only in the latter case we compute the correct arrival time of the leading tsunami wave. On the contrary the same numerical discretization applied to the Cartesian SWEs gives completely a erratic result.

r-adaptation for Shallow Water ows

We begin this chapter by presenting the SWEs in curvilinear coordinates.

Then we extend the ALE transport formulas and the volume transformation statement (the GCL), developed in chapter 2 for Cartesian coordinates, to general curvilinear coordinates. Using these results, the SWEs can be written in a framework in which points moves arbitrary on a sphere. We discuss two novel implementations (FV and RD) of the resulting ALE-SWEs which are basically a development of the algorithms seen in chapter 2, see sections 2.10 and 2.11. A simplied moving mesh method will allow mesh adaptation on a portion of the sphere. Finally all the new code features (Eulerian RD/FV on the sphere, ALE-RD/FV on the sphere, moving mesh, ALE coupling with the moving mesh) are tested against standard benchmarks.

SWEs in curvilinear coordinates

SWEs in curvilinear coordinates

In the community of ocean modeling and numerical weather prediction the rst step is to set accurately the SWEs on a rotating sphere S 2 . There are many ways to do that, each one with its own pros and cons. One of the most popular approach [START_REF] Ritchie | Application of the semi-lagrangian method to a spectral model of the shallow water equations[END_REF]; [START_REF] Pogorelov | A wave propagation algorithm for hyperbolic systems on curved manifold[END_REF]; [START_REF] Lauter | A discontinuous galerkin method for the shallow water equations in spherical triangular coordinates[END_REF] is to take advantage of a parametric representation of the sphere in a covariant vector basis. To allow a more complete understanding of this approach we have included, in appendix A, basic concepts from tensor analysis and dierential geometry that are used here and in the following. At every point P ∈ S 2 we x an orthogonal but not orthonormal reference system with basis vectors {g 1 , g 2 , g 3 } together with local curvilinear coordinates {z 1 , z 2 , z 3 }. We use this notation to avoid confusion with standard Cartesian coordinates {x 1 , x 2 , x 3 } adopted in all the previous chapters. A point P ∈ R 3 which belongs to the sphere and is specied in Euclidean space by the vector of coordinates x, can be expressed in curvilinear coordinates:

x = g 1 z 1 + g 2 z 2 + g 3 R
g 1 and g 2 dene the spherical surface where the undisturbed uid lies. R is sphere radius and g 3 is the axis "going out from the sphere". Along it, we measure the uid depth h(z) and the bathymetry/topography b(z). The uid free surface follows as η = h + b. We straightforwardly dene the metric tensor G ij = g i • g j and its determinant G = det G. Transforming the divergence operator of the SWEs (1.1), in curvilinear coordinates (formulas (A.9) and 4. ALE r-adaptive methods for the Shallaw Water equations on the sphere (A.10)) we get:

∂h ∂t + 1 √ G ∂ ∂z j √ G hu j = 0 (4.1) ∂hu i ∂t + 1 √ G ∂ ∂z j √ G T ij + S i = 0 (4.2)
hu i are the components of the discharge vector hu = hu i g i and T = huu ij + 1 2 G ij gh 2 are the components of the ux tensor huu = huu ij g i g j . We remark that the velocity components are not dened in a unitary basis. This, however, this allows to write the conservation laws in the compact form above. The source term reads

S i = G ij gh ∂b ∂z j + c F hu i + S i γ -S i c
Besides friction and bathymetry (rst two addends), the source term includes also a geometrical force due to Earth curvature and the ctitious force associated to Earth rotation. The aforementioned geometrical force

S i γ = Γ i jk T jk
comes from the transformation of the dierential operator, see the appendix A for the denition of the Christoel symbol Γ i jk . The rotation of the Earth is an important element in the dynamics of the atmosphere and ocean. Without loss of generality we let the axis of rotation be the Cartesian z-axis e 3 . If we solve the SWEs in the reference frame of the rotating Earth, we introduce in the momentum equations a pseudo-force of the form

S i c = f ε ij hu j
with f = 2Ω sin λ being the Coriolis parameter, λ the latitude and ε the 2D

Levi-Civita symbol. In all the experiments the Earth rotation rate is taken as

Ω = 7.292 × 10 -5 [s -1 ] unless it is dierently specied.
At this point the SWEs have the following non autonomous vector form

∂u ∂t + 1 √ G ∂F j ∂z j + S(z, u) = 0, (4.3) u = h hu i , F j = √ G hu j T ij , (4.4) 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 ε ij hu j S c

√

G a capacity coecient appears in the time derivative. For these reason, some authors speak about a capacity non autonomous form of the SWEs on the sphere. The ux Jacobian K(ξ, u)

= ∂F j ∂u ξ j is K(ξ, u) = √ G   0 ξ 1 ξ 2 G 11 ghξ 1 -uu • ξ u • ξ + uξ 1 uξ 2 G 22 ghξ 2 -vu • ξ vξ 1 u • ξ + vξ 2   (4.6)
for any vector ξ = ξ j g j . Note that the scalar product that appears in the expression above is not the actual inner product on the curved manifold but simply u • ξ = u j ξ j . It become the same if the basis is orthonormal. The Jacobian (4.6) admits a full set of real eigenvalues and linearly independent eigenvectors. The eigenvalues are

λ 1,3 (u, ξ) = √ G u • ξ ± gh ξ c , λ 2 (u, ξ) = √ Gu • ξ
where the norm is computed with respect to the inverse metric ξ c = ξ i G ij ξ j .

At this point we remark that the eigenstructure changes with respect to the Cartesian case. Particularly important is the scaling factor G ij that multiplies the acoustic part of the Jacobian, giving it the correct dimension. The Jacobian at rest on the sphere reads (cf. expression ??):

A(ξ, h) =   1 √ Gξ 1 √ Gξ 2 √ GG 11 gh ξ 1 0 0 √ GG 22 gh ξ 2 0 0   (4.7)
We mention that, starting from the work of Pons et al. [1998], many authors that implement FV in covariant basis, prefer to go back to an orthonormal framework (e.g. when solving the Riemann problem), namely G ij = δ ij and √ G = 1, we refer to [START_REF] Pogorelov | A wave propagation algorithm for hyperbolic systems on curved manifold[END_REF] for the details of this approach which requires a double projection: rst to obtain the orthonormal Riemann problem, then to go back to covariant uxes. Here we propose to use directly (4.6) which does not require these projections. After the resolution of the standard orthonormal Riemann problem, the uxes are projected back in the covariant basis to update the solution of the SWEs in covariant formulation.

ALE in curvilinear coordinates

To extend the SWEs in curvilinear coordinates and in ALE framework, we upgrade the kinematics relationships seen in chapter 2. The main results, namely transport formulas between Eulerian, Lagrangian and ALE and the volume/area conservation statements are recovered here for general curvilinear

Basic kinematics

We start this section setting the problem in the general case: a material uid B is a set of uid columns ("particles") that can be put in bijective correspondence with the points of a region B belonging to a dierentiable Riemannian manifold. We keep the same notation of chapter 2, the particles are P ∈ B and the points Q ∈ B. We assign a local covariant vector basis {g 1 , g 2 } and local curvilinear coordinates {z 1 , z 2 } to the tangent plane T Q B dened at each point Q. A point position is identied by the vector z which can be expressed in the local curvilinear coordinate system

z(Q) = z 1 g 1 + z 2 g 2 ≡ z i g i
The covariant vector basis is sometimes expressed as

{g i = ∂ ∂z i } ∈ T Q B.

ALE in curvilinear coordinates

Q 0 in the reference conguration is denoted by capital letters Z(Q 0 ). For this conguration we assign a local covariant vector basis {G 1 , G 2 } and local chart {Z 1 , Z 2 } to the tangent plane

T Q 0 B 0 Z(Q 0 ) = Z 1 G 1 + Z 2 G 2
As B moves, we assume the existence of a function that transform the point positions from the reference conguration to the actual one. This function is always called B, as in chapter 2:

B : B 0 → B z = B(Z, t)
Even if in principle the reference and the actual congurations could belong to dierent manifolds, points for us will be always attached to the sphere, that is

B = S 2 B and B 0 = S 2
B 0 as sketched in gure 4.2 and

B : S 2 B 0 → S 2

B

The function B is assumed to be continuously dierentiable, thus a smooth Jacobian of transformation can be dened

J B = ∂z ∂Z (4.8)
In order to avoid particle collision we add the smoothness requirement ∃J -1

B , J B = det J B > 0.
The relationship between Lagrangian and Eulerian time derivatives for a scalar φ and a vector v modies as:

dφ dt = ∂φ ∂t Z • B -1 = ∂φ ∂t z + ∂φ i ∂z j u j dv i g i dt = ∂v i g i ∂t Z • B -1 = ∂v i ∂t z + ∂v i ∂z j u j + Γ i jk u j v k g i
due to the presence of Christoel symbols in the gradient denition, (see (A.11)). The kinematic relation (2.5) expressing the derivative of a material volume modies as follows in curvilinear coordinates

∂ ∂t Z √ GJ B = J B ∂ ∂z j √ Gu j = 0 (4.9) or, more compactly ∂ t | Z ( √ GJ B ) -J B √ G∇ • u = 0.
For us (4.9) will be the time derivative of a material volume on the sphere. This result appears here in a slightly dierent form with respect to Savidis et al. [2008]. Since we prefer the above form for future developments we prefer to give its proof.

Proof of (4.9). We extend a classical proof, see for example the lecture notes 4. ALE r-adaptive methods for the Shallaw Water equations on the sphere of [START_REF] Mantegazza | Balance laws in eulerian and lagrangian form[END_REF], to curvilinear coordinates. From the Jacobian denition (4.8) we inherit dz m = J B dZ m where we recall dz m = dz 1 dz 2 . Passing to the transformation of area dierentials on manifolds B 0 → B, the statement must be rewritten. We consider the denition of area dierentials given in appendix A. In the actual conguration the innitesimal area on the tangent plane writes, see eq. (A.15)(A.16)

da = |dz 1 × dz 2 | = √ Gdz 1 dz 2 (4.10)
Taking the material derivative of (4.10):

∂ ∂t Z √ Gdz m = ∂ √ G ∂t Z dz m + d ∂z 1 ∂t Z √ Gdz 2 + d ∂z 2 ∂t Z √ Gdz 1 = ∂ √ G ∂t X dz m + ∂u j ∂z j √ Gdz m √ G = G(z(t)
) and we can use chain rule

∂ ∂t Z √ Gdz m = ∂z j ∂t Z ∂ √ G ∂z j dz m + ∂u j ∂z j √ Gdz m = u j ∂ √ G ∂z j dz m + ∂u j ∂z j √ Gdz m = ∂ ∂z j √ Gu j dz m = ∂ ∂z j √ Gu j J B dZ m
Developing the left-hand side with

√ Gdz m = √ GJ B dZ m we get: ∂ ∂t Z √ GJ B = J B ∂ ∂z j √ Gu j
which is exactly (4.9).

ALE kinematics in curvilinear coordinates

Mesh movement on the sphere will be accomplished through an arbitrary conguration A called ALE which we assume it belongs to a dierentiable manifold. As for the material conguration, A is composed of points Q. Its initial or reference conguration is A 0 , composed instead of points Q. A local coordinate system {ζ 1 , ζ 2 } can be assigned in the neighborhood of every Q such that the point's position writes locally

ζ( Q) = ζ i γ i where {γ i = ∂ ∂ζ i } ∈ T QA

ALE in curvilinear coordinates

reference ALE conguration. If, at each time value is associated a unique conguration, then the arbitrary conguration will be in motion:

A : A 0 → A z = A(ζ, t)
Although the ALE conguration changes with time dragging unsteady mesh movement, it will always belong to the sphere, that is A = S 2 A and A 0 = S 2 A 0 .

As sketched in gure 4.2 we have:

A : S 2 A 0 → S 2 A Using the denition of ALE time derivative ∂ ∂t ζ
, the ALE/mesh velocity of the point Q can be introduced in the current vector basis σ( Q) = σ j g j

σ j = ∂z j ∂t ζ
One can then use the chain rule to recover the relationship between Eulerian and ALE time derivative for scalars φ and vectors v, see (2.10) in Chapter 2:

∂φ ∂t ζ = ∂φ 0 ∂t ζ • A -1 = ∂φ ∂t z + ∂φ ∂z j ∂z j ∂t ζ = ∂φ ∂t z + ∂φ ∂z j σ j (4.11) ∂v i ∂t ζ = ∂v i 0 ∂t ζ • A -1 = ∂v i ∂t z + ∂v i ∂z j σ j (4.12)
As in Chapter 2, a useful relationship emerges between the Lagrangian time derivative and the ALE one (4.11) and (4.12):

dφ dt = ∂φ ∂t ζ + ∂φ ∂z j u j -σ j (4.13) dv i g i dt = ∂v i ∂t ζ + ∂v i ∂z j u j -σ j + Γ i jk u j v k g i (4.14)
We take the ALE time derivative of an innitesimal area on the sphere. This is called geometric conservation law (GCL):

∂ ∂t ζ √ GJ A = J A ∂ ∂z j √ Gσ j (4.15) or, more compactly ∂ t | ζ ( √ GJ A ) -J A √ G∇ • σ = 0.
Proof. We don't repeat the proof that follows the one done in the previous paragraph. In fact, this is achieved replacing body reference quantities 4. ALE r-adaptive methods for the Shallaw Water equations on the sphere with arbitrary reference ones, √ G 0 with Ĝ, J B with J A and dZ with dζ.

The GCL (4.15) is a geometrical relation between the ALE Jacobian and the ALE velocity. One can check that the integral of (4.15) over the whole sphere states the conservation of the area of the ALE domain. In the context of the numerical solution of ALE-PDEs on the sphere, this corresponds to preserve the total area of the mesh during the simulation. For this reason in section 4.5 we will discuss how to preserve also its discrete counterpart.

ALE-SWEs in curvilinear coordinates

We consider now the eld equations or balance laws for the water column in curvilinear coordinates:

∂ √ GJ B h ∂t Z = 0 and ∂ √ GJ B hu i g i ∂t Z = √ GJ B F i (4.16)
where F i represents the forces acting on the water column. We transform them from a Lagrangian description to the ALE framework. To do this, we can employ the transport formulas (4.13),(4.14), the area transformation relationships, (4.9),(4.15), all discussed in the previous sections. In order:

1. we transform the time derivative in (4.16) using (4.13) and (4.14)

2. we substitute (4.9)

3. we multiply the resulting expression by J A 4. nally we add (4.15) multiplied in turn by h or hu i

As a result, the SWEs in curvilinear coordinates in the ALE framework write:

∂ √ GJ A h ∂t ζ + J A ∂ ∂z j hu j - √ Ghσ j = 0 (4.17) ∂ √ GJ A hu i ∂t ζ + J A ∂ ∂z j T ij - √ Ghu i σ j + √ GJ A S i = 0 (4.18)
As done in previous chapter we benet of the vector form to write them compactly:

∂ ∂t ζ √ GJ A u + J A ∂ ∂z j F j - √ Gσ j u + √ GJ A S(z, u) = 0, (4.19) u = h hu i , F j = √ G hu j T ij 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 ε ij hu j S c

Mesh and geometry in curvilinear coordinates

SWEs are written for the water depth h and depth-integrated momentum hu.

So far we have seen that this form is the appropriate one to conserve, at a discrete level, total mass and momentum which remain for us the most important quantities to be conserved. We will show that our numerical approximations of (4.17) and (4.18) allow mass and momentum conservation on the sphere.

We also mention that in geophysical applications it could convenient to manipulate the SWEs in order to have discrete conservation of more interesting quantities such as total energy or potential vorticity, see for example [START_REF] Ringler | A potential enstrophy and energy conserving numerical scheme for solution of the shallow-water equations on a geodesic grid[END_REF]. This is left for future work.

In chapter 2 we have seen that this form is not particularly suitable to preserves the lake at rest on moving mesh. This is related to the fact that solving the lake at rest in a moving reference framework corresponds to evolve the ALE remap equation for the water depth h, to which it is associated a discretization error. To cure this problem we introduced the ALE-WB and ALE-PB form of the SWEs. The same can be done here. We write equation

∂ √ GJ A b ∂t ζ -J A ∂ ∂z j √ Gbσ j = 0
If we add the ALE remap for b to the mass equation (4.17) we nally obtain the WB-ALE form

∂ ∂t ζ √ GJ A w + J A ∂ ∂z j F j - √ Gσ j w + √ GJ A S(z, u) = 0, (4.21) u = η hu i

Mesh and geometry in curvilinear coordinates

Points on the sphere are represented by a proper parametrization {z 1 , z 2 } through the use of geometrical mappings. From now on, we will use a standard latitude longitude (lat-lon) parametrization for which the denitions of geometrical quantities that describes the sphere, including Jacobian, metric tensor and Christoel symbols are given in appendix C. We recall that for the 4. ALE r-adaptive methods for the Shallaw Water equations on the sphere sphere the aforementioned quantities are available in simple analytical expression. The main drawback of this approach is that it introduces a singularity at the poles. Here we are not interested in the simulation of tsunamis that propagate all around the globe and we will simply circumvent this issue solving the SWEs far from the polar region:

(z 1 , z 2 ) ∈ [-180 • , 180 • ] × [-75 • , 75 • ]
where the metric coecients are well dened with lat-lon coordinates. We mention that this is a big limit if one is interested in global circulation model where it is urgent to accurately compute the solution at the poles. In this case the use of dierent grids such as the cubed sphere grid Ronchi et al. [1996], the Yin-Yang grid [START_REF] Kageyama | yin-yang grid: An overset grid in spherical geometry[END_REF] or rotated lat-lon grids, allows to resolve the pole problem. The implementation of rotated lat-lon grids, see e.g. [START_REF] Rossmanith | Residual distribution schemes for hyperbolic balance laws in generalized coodinates[END_REF], is left for future work.

We consider a discretization of the spatial domain S 2 composed by non overlapping triangular elements which cover the sphere. For the mesh quantities, the same notation of chapter 1 is employed here, we refer to section 1.6.1 for the denitions. We just introduce a new notation for areas computed on the sphere, for which we use capital italic font (e.g. A), in order to dierentiate these from areas dened in the parametric space, for which standard italic is used as in 1.6.1. For instance the area of a spherical triangle and the median dual cell on the sphere are denoted as

|K| = K √ Gdz, |C i | = K∈D i |K| 3 (4.22)
Geometrical quantities computed in the parametric space are dened as in Cartesian case. For clarity we repeat

|K| = K dz, |C i | = K∈D i |K| 3
Finally we specify n j = n m j g m as the normal to the j-th face of the triangle scaled by the corresponding edge length.

DGCL in curvilinear coordinates

A rst look to the DGCL and its relevance for SWEs-ALE simulation has been given in chapter 2. Here we come through the main denitions given in section 2.5 and extends them to general smooth manifolds (as the sphere). When approximating the SWEs the most fundamental equations are the integral ones. It is thus more useful to consider the integral GCL over a general control volume V . The time discrete approximation of (4.15), namely the Discrete Geometric Conservation Law (DGCL) writes:

V √ Gdz n+1 - V √ Gdz n = ∆t ∂V √ Gσ m n m ds (4.23)
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

DGCL in curvilinear coordinates

or more compactly, using spherical areas

|V n+1 | -|V n | = ∆t j∈K υ j (4.24)
we have decomposed the left-hand side in the contributions of the three edges.

The edge velocity results:

υ j = ∂V j √ Gσ m n m j ds (4.25)
In chapter 2 we have presented a class of DGCL numerical schemes for the approximation of plane SWEs based on previous work, see [START_REF] Mavriplis | Construction of the discrete geometric conservation law for high-order time-accurate simulations on dynamic meshes[END_REF]. We have also seen that the presented closure collapses to the one of [START_REF] Lesoinne | Geometric conservation laws for ow problems with moving boundaries and deformable meshs, and their impact on aeroelastic computations[END_REF]. A rst dierence when setting the problem in curvilinear coordinates is that, for general metric tensor √ G = 1, we cannot dene the DGCL through the original characterization of [START_REF] Thomas | Geometric conservation law and its application to ow computations on moving grids[END_REF], i.e. the preservation of a uniform ow. On the sphere one could think to replace uniform ows with the preservation of zonal/meridional ows. But, even general numerical schemes on xed grid do not have this property, so one could question the importance to have it on moving meshes. Let's say however that we would like to conserve the total area of the mesh. In this sense the proof of Farhat comes to be less general: elements' areas are computed through quadrature formulas and it is hard to nd proper quadrature points in time and in space to evaluate the term t n+1 tn σ j n j dsdt such that eq.(2.23) results an identity (e.g. in this case midpoint does not assure exactness). Instead the closure of Mavripils is the one that we prefer and, with a minor modication, allows to respect the GCL.

Proposition (DGCL in curvilinear coordinates) A numerical method approximating eq. ( 4.3) on a smooth manifold M, veries the DGCL constraint or eq. (2.23) if each edge velocity is computed through the (signed) area swept on M by the j-th edge of the polygon in one time step

υ j = 1 ∆t Q j √ Gdz (4.26)
where the quadrangle Q j is dened by vertexes

z n i ,z n k ,z n+1 k ,z n+1 i
, i and k denotes the two vertexes of the j-th edge. We use the convention that, if the polygon is expanding, then the area swept is positive.

A brief comment is necessary. Even if complicated analytical formulas exist to compute the area of a spherical quadrilateral, υ j will be computed with a quadrature formula to approximate the integral in (2.25). This means that, 4. ALE r-adaptive methods for the Shallaw Water equations on the sphere at a later stage, we have to enforce explicitly the DGCL (2.23) to compute the spherical area in the new conguration

|V n+1 | = |V n | + ∆t j∈K υ j
Up to the our knowledge this is the rst attempt to nd DGCL uxes in curvilinear coordinates for general metric.

Finite Volume for SWEs on the sphere

There are two main complications that one has to face when attempting to solve the SWEs on the sphere: the rst concerns the non-autonomous nature of the SWEs on a manifold, and the second the resolution of the Riemann problem in curvilinear coordinates. Although the literature on Finite Volume schemes on the sphere is large Machenhauer et al. [2009], we refer, with respect to the two aforementioned issues, to the Lax-Wendro scheme of Rossmanith et al. [2004], and more closely to the MUSCL-type scheme implemented in Ullrich et al. [2010].

We will compute approximations of solution averages over the standard median dual cells that we denote as u i

u i (t) = 1 C i C i u(z, t) √ Gdz m (4.27) A common choice is to set C i = |C i | √ G i .
Then we dene the numerical ux along the boundary of the dual cell ∂C ij and the numerical source on dual cell C ij , such as in chapter 1 and 2:

F ij = 1 ∆t t n+1 t n ∂C ij F j n j ds dt F ALE ij = 1 ∆t t n+1 t n ∂C ij wσ m n m √ Gds dt = 1 ∆t t n+1 t n wυ ij dt S ij = 1 ∆t t n+1 t n C ij S √

Gdx dt

The ALE-FV discrete evolution equations then reads:

|C n+1 i |w * i = |C n i |w n i -∆t j∈D i R ij (w n , b n ) (4.28) |C n+1 i |w n+1 i = |C n i |w n i - ∆t 2 j∈D i R ij (u n , b n ) 2 + R ij (u * , b n+1 ) 2 (4.29)
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

4.6. Finite Volume for SWEs on the sphere with R ij = F ij + S ij . Please note how the capacity coecient have been naturally considered by applying the average preserving property (4.27). We make use to the Roe-type numerical uxes, as in the planar case (see (2.40))

(to which we refer for the notation):

F ij = 1 2 ( Fi + Fj ) • n ij -υ ij wj + wi 2 - K ij -υ ij I 3 2 (ȗ j -ȗi ) (4.30)
Note that Roe uxes have been tested in [START_REF] Ullrich | High-order nite-volume methods for the shallow water equations on the sphere[END_REF] for the low Mach number tests of [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF] showing adequate accuracy. We refresh the notation of chapter 2:

K ij = K(n ij , z * ij ) and R ij = R(n ij , z * ij ) are
respectively the ux Jacobian and the matrix of right eigenvectors, evaluated with a Roe linearization z * ij . Due to the spatially varying ux function, the property of the Roe average in so no more available, however, through the numerical ux (4.30), we still compute a local conservative solution that assures

F ij = -F ji .
In our implementation we have used the following average

z * ij =      √ h j + √ h j √ h i v i √ G i + √ h i v j √ G j √ h j + √ h j      (4.31)
For completeness concerning FV for spatially varying ux function we mention the dierent perspective of Bale et al. [2002]. Their ux wave decomposition allows to truly takes into account the non autonomous nature of the ux, directly into the resolution of the Riemann problem.

To accomplish the DGCL, we follow the strategy described in section 4.5 which consist of computing spherical area swept by each interface of the dual cell. The resulting interface velocities υ ij , dened in (4.25), are prescribed by (4.26)

υ ij = ∂C ij σ m n m √ G ds = K i,j ∂C K ij σ m n m √ G ds = K i,j Q K ij √ G dz (4.32) using the quadrilateral Q K ij dened by vertexes x n G K ,x n ij ,x n+1 ij ,x n+1 G K .
Once this value is computed in the code, the dual cell area can be updated in the new conguration:

C n+1 i = C n i + ∆t j∈D i υ ij (4.33)
4. ALE r-adaptive methods for the Shallaw Water equations on the sphere

FV algorithm in curvilinear coordinates

Although one can use the algorithm described in the previous section, additional diculties arise in the resolution of the Riemann problem in generalized coordinates. In particular the question of the basis in which we solve the Riemann problem becomes urgent since the basis' vector is changing from point to point. Using straightforwardly the Cartesian solver makes no sense. Indeed the ux jump at the interface between two dual cells has no meaning, being the ux F i and F j , composed by vectors and tensors dened on dierent basis.

There is also second issue: the Finite Volume discretization of the geometric source term related to Earth curvature.

Concerning the resolution of the Riemann Problem, Rossmanith et al.

[2004] propose to transform vectors F i and F j to a common coordinate system, for example the coordinate system at the interface position z ij , through the action of parallel transport. We propose an extension of their algorithm and, if possible, their conclusions to the case of unstructured grids. We dene the parallel transport operator P on the vector u as the vector eld tangent to the curve parametrized by the geodesic parameter λ:

P (u) = du dλ = d dλ (u m g m ) = du m dλ g m + dz n dλ dg m dz n u m = du k dλ + dz n dλ Γ k nm u m g k
Parallel transport is an operation that takes a vector u at a point A and transports it to a point B along the curve z i (λ). Parallel transport accomplishes this in such a way that u remains parallel to z i (λ). This constraint is the parallel transport equation:

P (u) = 0 ⇒ du k dλ + Γ k nm dz n dλ u m = 0
If one evaluate the linearization of the above expression at the interface gets:

u k ij = u k i + ∆λ 2 du k dλ i + O(∆λ 2 ) = u k i - ∆λ 2 dz n dλ Γ k nm u m i + O(∆λ 2 )
where the indexes m, n and k loop with standard summation convention while indexes i and j represent the nodes; node i in particular is the node around which we have performed linearization. Index ij denotes values at the interface between node i and node j. We can transport the full vector of conservative variables:

h hu k ij = 1 0 0 I 2 - ∆z n ij 2 0 0 0 Γ k nm i h hu m i
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. 4.6. Finite Volume for SWEs on the sphere in our implementation we have simplied ∆λ dz n dλ = ∆z n ij = z n j -z n i . We can compactly write the discrete parallel transport operator which transports the vector of conservative variable from node i to the interface ij:

u ij = P ij u i P ij = 1 0 0 I 2 - ∆z n ij 2 0 0 0 Γ k nm i = I 3 + B ij
dening a parallel transport matrix P ij . Once conservative variables are transported at the common interface, then the resolution of the Riemann problem can be carried, through (1.15)(1.16) or (1.17) (chapter 1). Let us reconsider the Riemann problem (on a xed grid for simplicity) dened at the interface between nodes i and j, both lying on a certain geodesic parametrized with λ.

The uxes in these cells are F i and F j respectively. Let us apply the appropriate parallel transport matrices so that F i and F j both get (approximately) represented in the coordinate basis at z ij . Relations (1.15)(1.16) become:

F ij = P ij F i • n ij - m p=1,β p <0 W p = P ji F j • n ij + m p=1,β p >0 W p
Once the solution jump has been transported to the interface, it can be decomposed in p-waves and resolved for the coecient α p W p = α p λ p r p , α p = l p (P ji u j -P ij u i )

First order numerical uxes are written in a compact form

F ij = 1 2 (F j + F i ) • n ij - 1 2 |K ij |(u j -u i ) + S γ ij (4.34)
where the following terms can be isolated in a source

S γ ij = 1 2 (B ji F j • n ij + B ij F i • n ij ) - 1 2 |K ij |(B ji u j -B ij u i ) (4.35)
To obtain a formally second order update we can write numerical uxes in the following form

F ij = 1 2 ( Fj -Fi ) • n ij + ( Fi -F i ) • n ij - 1 2 |K ij |(ȗ j -ȗi ) + Sγ ij
This approximate source term result enhanced by a second order correction:

Sγ ij = S γ ij + S γ- ij
4. ALE r-adaptive methods for the Shallaw Water equations on the sphere

S γ- ij = B ij Fi -F i • n ij
What does this source term Sγ ij represent? On structured quadrilateral (Cartesian ) grids this corresponds to have directly an approximation for the geometric source term S γ in (4.3). In this case, the application of parallel transport within the wave propagation algorithm produces a natural numerical approximation of the geometric source term, which was the second issue reported.

Proof. Consider a Cartesian grid on the sphere with the following nodes:

z 1 i = z 1 l + i - 1 2 ∆z 1 z 2 j = z 2 l + j - 1 2 ∆z 2 with z 1 l , z 2
l is the lower left corner of the computational rectangular domain.

Let's consider a median dual cell area in gure 4.3 centered in A = (z 1 i , z 2 j ). Let's consider the vertical interface between node A and B = (z 1 i-1 , z 2 j ) where n 1 = -∆z 2 and n 2 = 0. Moreover

∆z 1 i-1,i = z 1 i -z 1 i-1 = ∆z 1 ∆z i,i-1 = z 1 i-1 -z 1 i = -∆z 1
We consider the central part of the source term (4.35) that, on a Cartesian grid reads

S γ i-1/2,j = 1 2 ∆z 1 (-∆z 2 ) 0 -a 2 Γ k m1 huu m1 i-1,j - 1 2 (-∆z 1 )(-∆z 2 ) 0 -a 2 Γ k m1 huu m1 i,j
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

4.6. Finite Volume for SWEs on the sphere Let's consider the horizontal interface between node A and C = (z 1 i , z 2 j-1 ) where n 1 = 0, n 2 = -∆z 1 and

∆z 2 j-1,j = z 2 j -z 2 j-1 = ∆z 2 ∆z 2 j,j-1 = z 2 j-1 -z 2 i = -∆z 2 and S γ i,j-1/2 = 1 2 (-∆z 1 )∆z 2 0 -a 2 Γ k m2 huu m2 i,j-1 - 1 2 (-∆z 1 )(-∆z 2 ) 0 -a 2 Γ k m2 huu m2 i,j
Summing all the contributions from the four interfaces, taking into account the source term appears with a minus sign in the global update (4.28):

S γ ij = -∆z 1 ∆z 2 0 q1 ω q1 a q1 Γ k 1n huu 1n q1 + q2 ω q2 a q2 Γ k 2n huu 2n q2 ≈ ∆z 1 ∆z 2 S γ (z) dz
we get a quadrature formula for the integral of the geometric source term. Different quadrature points are used for the the two terms containing Christoel symbols Γ k 1m and Γ k 1m . For the rst, quadrature points are placed on the mid- point of the upper and lower edge of the rectangular dual cell, for the latter the two quadrature points coincide with the midpoints of left and right edges, see again gure 4.3. Similar arguments hold for the upwind part.

For unstructured grids such an approximation of the geometric source is not obvious since the inclusion of parallel transport into the resolution of the Riemann Problem can become complicated. In fact, taking the central part of the source term (4.35) we have

B ij F i • n ij = 0 - a∆z n ij 2 Γ k mn huu mq n q ij i
which includes but not coincides with a discretization of the geometric source term. We will refer to this formulation as the one of Rossmanith et al. [2004] since it is an extension of their work to unstructured grids.

An alternative yet simpler approximation of the geometric source consists in approximating S γ directly on the underlying grid as done by Ullrich et al.

[2010] on a structured grid. An extension of this last work to unstructured grids consists in using the pointwise value:

C i S γ √ Gdz ≈ |C i | G i S γ i = |C i | G i 0 Γ k mn huu mn i (4.36)
such a discretization is second order accurate. We have tested both the formulations in gure 4.4 for a circular hump propagating on a sphere. Even if they provide similar accuracy, the complexity of (4.35), together with the additional cost to include high order reconstruction, drove us to the implement directly (4.36). This will represent the nal implementation of the FV.

=     h j -h i 0 0     , ∆b =     b j -b i 0 0    
and the following Jacobian at rest

A ij = A(n ij , h ij ) and A - ij = A(n ij , h - ij ) with
average values (see (4.7)):

h ij = hi √ G i G km i + hj √ G j G km j 2 , h - ij = √ G i G km i hi + √ G j G km j h i 2 The central part S c ij = 1 2 A ij ∆b + A - ij ∆b - ij (4.37)
and the upwind dissipation term :

S * ij = - K ij -υ ij I 3 2 ∆b ij (4.38)
However for the FV scheme presented so far we cannot prove Well Balanced.

The reason is that we cannot use the Roe average property in the demonstration of (1.6.2), spoiling Well Balancedness. To cure this main problem we have isolated the hydrostatic part of the ux and the geometric source term.

Setting Π = 1 2 gh 2 we have to approximate

C ij P G ij n j √ Gdz + C ij Γ i jk P G jk √ Gdz (4.39)
We now write the above quantity into a non conservative form. Putting the ux (rst term) in divergence form and deriving term by term we get

1 √ G ∂ ∂z j √ GΠG ij = Γ m mj ΠG ij + ∂G ij ∂z j Π + ∂Π ∂z j G ij
where we have used relation (A.6). We add the metric term and using Ricci's Lemma (A.7) the rst, second and fourth term cancel out. We can simplify expression (4.39):

Π ∂G ij ∂z j + G ij Γ m mj + ∂Π ∂z j G ij + ΠΓ i jk G jk = gh ∂h ∂z j G ij
4. ALE r-adaptive methods for the Shallaw Water equations on the sphere It will be this non conservative form of the hydrostatic term that we approximate as

C ij √ GG ij gh ∂ ∂z j 0 h dz = 1 2 A ij ∆h ij + A - ij ∆h - ij
As in the Cartesian case two contributions appear, one related to rst-order piece-wise constant solution approximation, the second one comes from the integration of the solution's variations due to the gradient-reconstruction. We remark that the presence of variable-coecients in the matrix A makes the non-conservative discretization no longer equivalent to the conservative one.

Strictly speaking, we are loosing exact conservation to the extent that Rankine-Hugoniot condition are no more satised exactly at a discrete level. We remark however that the variation of the metric coecients are very smooth. This problem of Well-Balanced scheme on the sphere is already known and we note that our WB approximation on the sphere resembles to the work of Ullrich et al. [2010]. For dierent WB scheme for the SWEs on the sphere, we refer to Rossmanith et al. [2004] and [START_REF] Yang | A parallel well-balanced nite volume method for shallow water equations with topography on the cubed-sphere[END_REF].

Now we have all the formulas to prove Well Balancedness with moving meshes.

Proposition 1. The nite volume discrete equations (4.28)-(4.29) with denitions (4.30), (4.32),(4.37) and (4.38) veries the Well-Balanced property both on moving and xed mesh, provided that the same reconstruction procedure is used for u and b.

Proof. For the constant lake at rest solution w 0 = [const 0 0 ] T , numerical ux and the source term reduce to

F ij = 1 2 A ij ∆h ij + A - ij ∆h - ij -υ ij w 0 - |A ij -υ ij I 3 | 2 (ȗ j -ȗi ) S ij = 1 2 A ij ∆b ij + A - ij ∆b - ij - |A ij -υ ij I 3 | 2 ∆b ij
Note now that ȗjȗi + ∆b ij = wjwi which vanishes by hypothesis. We are left with the DGCL

|C n+1 i |w n+1 i = |C n i |w n 0 -w 0 ∆t j∈D i υ ij
which is veried by hypothesis (4.33). This complete the proof:

w n+1 i = |C n i | -∆t j∈D i υ ij w 0 |C n+1 i | = w 0
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Residual Distribution for SWEs on the sphere

We remark that the proof does not rely, as usually done for the Cartesian coordinates, on the property of the Roe average, which again is not available in this case.

4.7 Residual Distribution for SWEs on the sphere Steady Residual Distribution have been extended to the rotating sphere in [START_REF] Rossmanith | Residual distribution schemes for hyperbolic balance laws in generalized coodinates[END_REF]. In this section we propose a high order extension for time dependent problem on moving meshes. This is the rst such generalization of RD to our knowledge. To obtain second order of accuracy in time, we extend the eRK2 ALE-RD of [START_REF] Arpaia | An ale formulation for explicit runge-kutta residual distribution[END_REF] detailed in chapter 2 and appendix D. Compared to FV, the extension to the sphere appears almost straightforward. In the following we provide directly the algorithm that we comment extensively later:

Given the uctuation, see (4.7):

Φ K (w h , b h ) = ∂K F j (u h )n j ds + K σ j ∂w ∂z j √ Gdz + K S(u h , z) √ Gdz
If not specied, for us K = 1 2 (K n + K n+1 ) in the following. Keeping in mind this, the two steps writes 1] Predictor step: for each element K ∈ T h

• Compute the residual Φ K(1) = Φ K (w n h , b n h ).
• Distribute the uctuation to the nodes of K such that j∈K Φ

K(1) j = Φ K(1)
• Compute the rst order prediction of the solution, denoted as w *

w * i = w n i - ∆t |C n+1 i | K∈D i Φ K(1) i (w n h , b n h ) (4.40) 2] Corrector step: for each element K ∈ T h • Compute the residual Φ K(2) = 1 2 Φ K (w n h , b n h ) + Φ K (w * h , b n+1 h ) + K w * h -w n h ∆t √ Gdz (4.41)
• Distribute the uctuation to the nodes of K such that j∈K Φ

K(2) j = Φ K(2)
4. ALE r-adaptive methods for the Shallaw Water equations on the sphere • Compute the second order correction from

w n+1 i = w * i - ∆t |C n+1 i | K∈D i Φ K(2) i (w * h , w n h , b n h , b n+1 h ) (4.42)
All the integral are still approximated via Gaussian quadrature. We believe it is worth noting both the simplicity and the rigor of the RD extension from

Cartesian coordinates to curvilinear ones, specially if compared to FV. In particular:

1. the heterogeneity of the medium is taken into account integrating the coecient √ G (hidden in the uxes or in the area denition) with the precision degree of the quadrature formula.

2. the geometric source term comes from the fact that the ux is written in curvilinear coordinates. As for any other source term, it is included automatically in the wave propagation algorithm through the splitting approach, see section 1.6.3. We recall that for FV an analogous approach was much more complicated.

RD mass matrix computation on manifold

In this paragraph we detail the practical computation of the split residual (4.41). In particular it is worth to report how we compute the unsteady term in (4.41) since this could impact consistency and (we will see) conservation

properties on the sphere. First, we compute the following important quantities as:

m GAL ij = K ϕ i ϕ j √ G dz ≈ K 12 (1 + δ ij ) |K j | = K ϕ j √ G dz ≈ K 3 (4.43)
where the |K| has been dened in (4.22). We refer to m GAL ij as the Galerkin mass matrix. On Cartesian plane the above denitions collapse respectively to the Galerkin mass matrix and to |K j | = |K|

3

. In chapter 1 we have already mentioned that the time discretization in RD methods is obtained from an analogy with the theory of Finite Element. From the weak form it is well known that an element mass matrix of the form

m K ij = K w i ϕ j √ G dz
β K i = 1 K √ G dz K w i √ G dz (4.44)
which is veried, for example, by:

w i = ϕ i + β K i -ℵ K i , m K ij = m GAL ij + |K j |(β K i -ℵ K i ) (4.45)
where ℵ K i is:

ℵ K i = |K i ||K| -1 (4.46)
With this in mind, we can treat separately the splitting of time residual (4.41).

Invoking the analogy RD/Finite Elements (4.44) we let a mass matrix appears, mass matrix that can be computed with the consistent w i of (4.45). At the end, the corrector residual (4.41) will be distributed according to:

Φ K(2) i = β K i 1 2 Φ K (w * ) + Φ K (w n ) + j∈K K j w * j -w n j ∆t + j∈K δm K ij w * j -w n j ∆t (4.47)
We see that the last two terms in the above splitting (4.47) are related to the the three contributions of the mass matrix (4.45). In fact, two of them have been grouped in the so-called dissipation matrix, see [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF]:

δm K ij = m GAL ij -|K j |ℵ K i (4.48)

ALE part

We treat the ALE residual as done in chapter 2:

Φ K,ALE = K √ Gσ j ∂w ∂z j dz = K ∂ √ Gσ j w ∂z j -w ∂ √ Gσ j ∂z j dz = ∂K √ Gwσ j n j ds - K w ∂ √ Gσ j ∂z j dz
We examine the rst part. Edge uxes are computed by quadrature formula.

Moreover these uxes should verify the DGCL which means that every edge should sweep the area ∆tυ j . We can split these areas on the edge quadrature points through the weights

∂K

√

Gwσ j n j ds = j∈K N Q q ω q w q υ j 4. ALE r-adaptive methods for the Shallaw Water equations on the sphere while the second term using the usual passages

K w ∂ √ Gσ j ∂x j dz = K w √ Gdz (∇ • σ) K = j∈K ϕ j √ Gdz w j (∇ • σ) K = j∈K |K j |w j ∆|K| ∆t|K| = j∈K |K j |w j |K| ∆|K| ∆t
where a weighted average solution on the cell appears. Finally using denition (4.43) for this term we get form

K w ∂ √ Gσ j ∂x j dz = j∈K w j 3 ∆|K| ∆t (4.49)

Properties of the RD scheme

In the following we analyze some of the properties of the scheme proposed.

Conservation

Flows at geophysical scale are typically very smooth and for numerical methods approximating very smooth phenomena, the verication of Rankine-Hugoniot conditions could be not as urgent as it was for coastal and hydraulic applications. However global conservation of the ow variables on the sphere remains an important issue. It is interesting to check the global conservation statement for RD, in light of all the denitions given. For sake of clarity we take the source term equal to zero, since it will not aect the analysis. We sum the corrector (4.42) over all the nodes of the triangulation i ∈ T h :

i∈T h C n+1 i w n+1 -w * ∆t = - K∈T h j∈K Φ K(2) j
As seen in previous chapters, uxes in divergence form cancel in the sums. On the right hand side, (4.47) reduce to the sum of the unsteady terms and the ALE term (4.49)

i∈T h C n+1 i w n+1 i -w * i ∆t = - K∈T h j∈K β K j m∈K     w * m + w n m 2 ∆|K| 3∆t

Residual Distribution for SWEs on the sphere

Recalling the denition of the dissipation mass matrix (4.48), we observe that:

i∈K j∈K δm K ij (w * j -w n j ) = 0
if the following condition is veried:

j∈K m GAL ij = |K i | (4.50)
which is true for the denitions given in (4.43). Now we use consistency

j∈K β K j = I 3 : i∈T h C n+1 i w n+1 i -C n+1 i w * i = - i∈T h K∈D i w * i + w n i 2 ∆|K| 3 + |K| 3 (w * i -w n i )
where we have also used again K∈T h j∈K = i∈T h K∈D i . We sum and subtract i∈T h K n i w n i to the left hand side and, using denition (4.22) we get:

i∈T h C n+1 i w n+1 i -C n i w n i = - i∈T h K∈D i w * j + w n j 2 ∆|K| 3 - |∆K| 3 w * j -w n j 2
which proves that we conserves the ow variables over the sphere:

i∈T h C n+1 i w n+1 i = i∈T h C n i w n i
Well Balancedness for moving mesh on the sphere

The discretization of the hydrostatic term have to ensure the discrete balance between hydrostatic and bottom forces when approximating a lake at rest ow.

We can check that ( 1 2 gh 2 h = Π for shortness):

Φ K = ∂K F j (u h )n j ds + K S(u h , b h ) √ Gdz (4.51) = 0 ∂K Πn j G ij √ Gds + 0 K gh h ∂b h ∂z j G ij √ G + 0 K Π Γ i jk G jk √ G = 0
the residual is not zero. Instead, for Cartesian coordinates ( √ G = 1) under the assumption of P 1 and exact integration, the above statement is an equality. Unfortunately, in generalized curvilinear coordinates we loose these discrete equivalences because of the non linear metric coecients which spoils the exactness of the integrals. As done for FV, it is convenient to rewrite the hydrostatic pressure term in a non-conservative form Proof. We do the proof only for corrector, the predictor is a particular case.

1 √ G ∂ ∂z j √ G ΠG ij + Γ i jk ΠG jk = gh ∂h ∂z j G ij
We impose w 0 = [const 0 0 ] T into the residual (4.41) (∆w 0 = 0):

Φ K(2) = ∂K F j 0 n j ds -w 0 j∈K N Q q ω q υ j + w 0 ∆|K| ∆t + K S 0 √ Gdz = w 0 ∆|K| ∆t - j∈K υ j + ∂K F j 0 n j ds + K S 0 √ Gdz = 0
where we have used the property of Gaussian weights N Q q ω q = 1 to develop the second term. The term in brackets is zero if we compute the edge velocities with (2.25). We are left with the balance between hydrostatic ux and bathymetry terms. The ux term is now transformed to a non conservative form (4.52) which is more suited to accomplish Well Balancedness. We integrate the ux term in the divergence form and we can write the balance (4.51)

as 0 K gh h ∂h h ∂x j G ij √ G + 0 K gh h ∂b h ∂x j G ij √ G = 0
which is zero if the same linear piecewise continuous approximation is used for h, b. To have Well Balancedness it is sucient to employ the same quadrature formula for the two integrals. Under these hypothesis Φ K(2) = 0 and we have w n+1 i = w 0 .

We remark that the integrals, even under the piecewise linear approximation, cannot be computed exactly, spoiling exact conservation. This loss of conservation is the sense that Rankine-Hugoniot condition are no more satised exactly at a discrete level. However the integral that has to be computed is in the following form:

K ϕ i G ij √ Gdz
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Moving mesh on the sphere

The metric coecients, which are innitely smooth on the sphere, could be expanded in Taylor series and the whole term could be integrated with an arbitrary high degree of precision making the error much smaller then the discretization error.

Moving mesh on the sphere

Although h-renement on the sphere has already received considerable attention (see the Adaptive Mesh Renement methods of [START_REF] Behrens | An adaptive semi-lagrangian advection scheme and its parallelization[END_REF][START_REF] Behrens | Atmospheric and ocean modeling with an adaptive nite element solver for the shallow-water equations[END_REF]; [START_REF] Hubbard | A three-dimensional, adaptive, godunov-type model for global atmospheric ows[END_REF]), the moving mesh method on the sphere is a quite recent subject of research. Up to the author's knowledge there are only two published works: the harmonic map of [START_REF] Di | Moving mesh methods for singular problems on a sphere using perturbed harmonic mappings[END_REF] and the Monge-Ampere equation of [START_REF] Weller | Mesh adaptation on the sphere using optimal transport and the numerical solution of a mongeampère type equation[END_REF]. In chapter 3 we have seen that a powerful and general tool for mesh adaptation is the theory of harmonic map (see also appendix B). Immediately one can think to use this general framework to map from the computational manifold with its own metric to the physical mesh which lie on a sphere (in the same fashion as we where considering, in two dimensional case, the Euclidean plane):

A : M → S 2 A z = A(ζ, t), J A = ∂z ∂ζ
Unfortunately such a mapping violates the hypothesis of the HYS theorem:

the positive curvature of a sphere makes harmonic map not unique. Keeping in mind that we will deal with problems at large scale but not at global scale we can hope to simplify the problem. [START_REF] Weller | Mesh adaptation on the sphere using optimal transport and the numerical solution of a mongeampère type equation[END_REF] for example suggest that a naif idea could be to project the mesh on a plane, adapt on the plane with a certain MMPDE from the ones seen in chapter 3 and nally project back to the sphere. The same authors, which were always interested in global circulation models, were then forced to move away from this because the problem was moved on how to impose proper boundary conditions to the mesh node position, for examples at the poles. However, if one is not interested in moving the mesh around the poles, this idea remains valid.

We assume that the computational domain is described by a certain parametrization of the sphere ζ with metric components γ ij and square root of the metric tensor √ γ. We consider then following direct transformation from the sphere to a two dimensional Euclidean plane:

A : S 2 A 0 → E 2 z = A(ζ, t), J A = ∂z ∂ζ
We can extend the functional of [START_REF] Ceniceros | An ecient dynamically adaptive mesh for potentially singular solutions[END_REF] with minor modi-4. ALE r-adaptive methods for the Shallaw Water equations on the sphere cations (recall the ALE notation ζ = ζ i γ i )

E(z) = 1 2 Ω ζ ωγ i ∂z α ∂ζ i • γ j ∂z α ∂ζ j √ γdζ = 1 2 Ω ζ ωγ ij ∂z α ∂ζ i ∂z α ∂ζ j √ γdζ (4.52)
to which the following Euler-Lagrange equations or MMPDE is related

1 √ γ ∂ ∂ζ i ω √ γγ ij ∂z α ∂ζ i = 0 (4.53)
This will be our MMPDE on the sphere. It consist, as in the two dimensional case, of a decoupled system of PDEs. Minor modication in the algorithm must be done to implement the Laplace-Beltrami operator, see (A.13), instead of the Laplacian. We recall also that √ γ and γ ij are dened in the computational domain (greek letters): they don't depend from that physical conguration.

The only non-linearity is associated to the monitor function ω = ω(z) which depends on the evolving underlying solution.

Mesh generation

The moving mesh algorithm dened by the sequence of iterations (3.16) can be used to approximate our new MMPDE (4.53). We remark that, as done for the CFD, we do not solve the MMPDE in the polar regions. Roughly speaking, points will not move at the poles.

We extend the benchmarks proposed in [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF] and tested in section 3.4.5 to the sphere. The monitor function is computed according to ω = √ 1 + αu 2 , with u assigned:

u(ψ, λ) = exp(-8(4ψ 2 + 9λ 2 -1) 2 ) schemes in spherical coordinates and in presence of Earth rotation. It consists of a zonal ow with the corresponding geostrophic height eld. Typically several rotation orientations are specied. In our case, in order to avoid ow entering/leaving the pole regions, only the case where the axis of rotation coincides with the z-axis will be simulated. The velocity and height elds are initially (and for all time) given by:

h(z, 0) = h 0 - 1 g ΩR u 0 + u 2 0 2 sin z 2 sin z 2 u 1 (z, 0) = u 0 cos z 2 u 2 (z, 0) = 0 with gh 0 = 2.94 × 10 4 [m] 2 [s] -2 and u 0 = 2πR 12days
. Simulation are performed on 5 unstructured grids, halving the mesh size. The meshes are generated always with the meshing software Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensionalnite element mesh generator with built-in pre-and post-processing facilities[END_REF][START_REF] Geuzaine | Gmsh: a three-dimensionalnite element mesh generator with built-in pre-and post-processing facilities[END_REF],

with the frontal algorithm (Mesh.Algorithm=6). In gure 4.7 we can see the mesh topology for the rst three levels of mesh renement.

The error norm

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

4.9. Numerical experiments is computed following [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF]:

||e|| L2 = 2π 0 π/2 -π/2 (h(z) -h ex ) 2 √ G dz
In gure 4.11 we report the dimensionless error against the mesh size, after 5 days of computations, as specied by always in [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF] l

2 = ||e|| L2 ||h ex || L2
For this steady case we have run only xed mesh computations. We observe that the error is of the same order of magnitude for FV and RD. For both schemes it decreases with more then second order of accuracy, this being particularly evident for FV.

Advection of cosine bell

We revisit test case #1 from the Williamson test suite [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF] in order to assess numerical accuracy for unsteady problems and in presence of mesh movement. A cosine bell is transported once along the equator (R = 6371220, g = 9.80616 and Ω = 0):

h(z, 0) = h 0 2 1 + cos πr π/3 + h 0 if r < π 3 h 0 otherwise u 1 (z, 0) = u 0 cos z 2 (4.56) u 2 (z, 0) = 0
where h 0 = 1000 and u 0 = 2πR 12days . r is the great circle distance between z and the center of the cosine bell z C :

r = R arccos sin z 2 C sin z 2 + cos z 2 C cos z 2 cos(z 1 -z 1 C )
In general several orientations for the wind are specied but so far we have tried only the zonal conguration to avoid interaction with polar regions. For this test, Williamson propose to test the advective components in isolation.

We prefer to test all the SWEs using the technique of manifactured solution.

See appendix E for the derivation of the source term to be added to the SWEs in order that (E.1) is an exact solution. Moreover, with respect to Williamson original solution we have translated the initial height of h 0 to avoid the appearance of dry regions. Simulation are performed on 5 unstructured grids generated with the meshing software Gmsh, halving the mesh size. In gure 4.9 there is an example of the grid topology for the rst three levels of mesh renement. With respect to the steady case, the order of accuracy fall below 4. ALE r-adaptive methods for the Shallaw Water equations on the sphere the second order. RD seems to be slightly more accurate in terms of global error and order of accuracy (1.88 against 1.78), see gure 4.11.

To test the accuracy of the ALE scheme we have added an unsteady perturbation to the position of mesh nodes, see gure 4.10. They will be mapped according to:

z 1 (t) = ζ 1 + 0.3 cos (ζ 1 ) cos (ζ 2 ) sin 4πt T z 2 (t) = ζ 2 + 0.3 cos (ζ 1 ) cos (ζ 2 ) sin 4πt
T Always in gure 4.11, with dashed lines are reported the convergence lines for the ALE scheme. We observe that, for both RD and FV, the unsteady and arbitrary mesh perturbation does not spoil the order of accuracy.

Circular hump on a non rotating sphere

We consider the test contained in Rossmanith et al. [2004]. For the initial conditions we place a circular depth disturbance at the equator: semi-spherical grids: a coarse one with 7122 grid points and 14034 elements; a ne one with 39699 and 78900 elements. Isolines (only for FV) are plotted in gure 4.12 for the coarse and the ne simulation (rst two rows), showing a good shock capturing and good symmetry in all radial directions.

h(z, 0) = 2 if arccos(cos(z 1 ) cos(z 2 )) ≤ 0.2 0.2 otherwise u(z, 0) = 0 This initial condition is symmetric about the point (z 1 , z 2 ) = (0 • , 0 • ).
Then the coarse mesh is used as the computational domain for a moving mesh ALE simulation. The moving mesh parameter are α = 50 and β = γ = 0.15 for RD and α = 50 and β = γ = 0.3 for FV. In gure 4.12 (last row) we can observe that the meshes exhibit a strong adaptation on the shocks and also a good symmetry in the radial direction. As a consequence the ALE simulation provides a good resolution of both the inner and the outer shock, of the order of simulation run on the xed ne mesh.

In gure 4.13 we show some cuts of the uid depth along the equator. In the same gure our FV and RD implementations are also compared against the numerical solution of [START_REF] Pogorelov | A wave propagation algorithm for hyperbolic systems on curved manifold[END_REF] which consists in a high-resolution FV (Lax-Wendro ux with MC limiter) on a Cartesian grid composed of 34680 points. From the cuts and from the table 4.1 reporting the CPU times, we can conclude that the moving mesh method is eective. We can observe that the shock resolution improves with respect to the xed coarse mesh and we also see a consistent CPU time reduction with respect to the xed ne mesh (58% for RD and 65% for FV).

Circular hump on a rotating sphere

This test case is also taken form [START_REF] Pogorelov | A wave propagation algorithm for hyperbolic systems on curved manifold[END_REF]. It is an extension of the previous one in presence of the sphere rotation. The dimensionless rotation rate is Ω = 5. We have run this test with two semi-spherical grids: a coarse one with 8361 nodes and 16488 elements and a ne one with 46521 nodes and 92488 4. ALE r-adaptive methods for the Shallaw Water equations on the sphere 4. ALE r-adaptive methods for the Shallaw Water equations on the sphere elements. As before the coarse mesh represents the computational domain for an moving mesh ALE simulation. The moving mesh parameter are α = 50 and β = γ = 0.1. In gure 4.14 we report snapshots of the depth isolines and of the adaptive meshes. This times mesh nodes are clustered both around the shock waves and also near some smoother feature where high gradients appear. We remark that, for this case, FV gives better and smoother result with respect to RD where some oscillations appear (visible in gure 4.14). This is probably related to a non proper activation of the limiter and it subject of current investigation.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. for example in [START_REF] Mori | The 2011 tohoku earthquake tsunami joint survey group (2012), nationalwide post event survey and analysis of the 2011 tohoku earthquake tsunami[END_REF]. This big amount of data helped researchers in better understanding the phenomenon. For instance, data were used to recover initial waveform; this inverse problem is explained in section 5.1 with all the physical hypotheses that are introduced. Both buoy data and and eld survey are also very useful to test the capability of hydrodynamics codes in modeling the main physical features of a tsunami event. In section 5.2 and 5.3 we compare the accuracy of the numerical simulation against observations.

Initial waveform

The selection of the initial waveform is crucial to obtain realistic wave heights and runups. First, the waveform depends on the fault's rupture dynamics.

Then the rupture pulse produces seismic waves which propagate to the seaoor.

Finally the seaoor displacement generates wave motion in the ocean up to the free surface. The resulting waveform represents the initial condition for the SWEs code. As we see the phenomenon is complex, involving non linear wave propagation in dierent media. Typically, in the context of tsunami modeling, it can be simplied using many physical assumptions.

Fault rupture model

First of all we need to understand how to properly account for the rupture mechanism. The rupture geometry is described in the sketch in gure 5.1. We Tsunami height: the height, relative to a stated reference level, to which a particular location at certain time instant is covered by water.

Tohoku-Honsu tsunami simulation

can see two plates moving w.r.t one another having the fault surface Ω F as interface; the movement is characterized by

• two angles: the strike angle formed between the nord and the alongstrike direction and the dip angle formed by the horizontal and the downdip direction • an heterogeneous eld of displacement vectors contained in the fault plane. The relevant quantity is the module of the slip vector which is simply called the slip and for us is d(x, t) ∈ R.

Fault's rupture is a complex phenomenon which depends from pre-rupture stress condition, geometrical setting and frictional property of the fault which are largely unknowns. To simplify the problem we approximate the scalar slip eld: the fault plane is subdivided in N F rectangular subfaults, each one with its slip dened at the center d i , i = 1, N F , see the right picture of gure 5.1.

If the dynamics of the rupture is taken into account, the time domain is also approximated and subdivided in multiple windows allowing a time dependent activation of the subfaults. The approximated slip eld d i is also called source model. The computation of the slip eld d i from available time dependent data (both seismometeres and tsunamometers) is referred to as the inversion problem. This inversion problem could be very hard to solve and strongly depends on the algorithm used and on the data, leading to quite dierent slip elds for the same earthquake/tsunami event. For the 2011 Tohoku tsunami we have assisted to a proliferation of source models among which is quite dicult to select a better one. Moreover the choice of the specic inversion is determinant to reproduce accurately wave heights and runups. This issue was addressed recently in MacInnes et al. [2013] and in the TANDEM project. The source model of [START_REF] Satake | Time and space distribution of coseismic slip of the 2011 tohoku earthquake as inferred from tsunami waveform data[END_REF] emerged as one of the best for matching both nearshore and oshore height records as well as inundations. Satake used 55 subfaults, each 50 [km] long and 50 [km] wide, only the shallower column of subfaults having an halved width of 25 [km]. The strike and slip angles are adopted from USGS W-phase inversion result (strike 193 • , slip angle 81 • ) while the dip angle varies along the depth from 8 • for the shallowest subfaults to 16 • for the deeper ones. It has been realized that, for great earthquakes, the rupture propagates over large distances in a nite time, in the order of a few minutes. To consider these time dependent eects in the inversion, Satake et al. [2013] used a delayed rupture mechanism in such a way that, at the end, a temporal variation of the slip eld was included. 5 timewindows were considered: every subfault activates at time m∆t + t d . The time interval was set to ∆t = 30 [s] based on physical observations; m = 0, 1, 2, 3, 4 is an integer that represents each time-window and t d = |x i |/V r is a delay which depends on the rupture velocity V r and on the distance |x i | of the subafault's center from the hypocenter. Satake computed 275 unknowns values of slip ( 55Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Initial waveform

subfaults to be computed at each of the ve time windows) based on tsunami waveform data (GPS buoys, tidal gauges and DART buoys). Practically, the slip distribution obtained with this inversion algorithm is sampled every 30 [s], from 0 [s] (the time origin is 14:46:18 Japan Standard Time) until 300 [s] when the rupture mechanism stops. The result of the inversion can be found in the original paper Satake et al. [2013] and we don't report them here. For clarity, at each sampling time, the slip looks like the right picture of gure 5.1.

Okada dislocation theory

Seismic wave propagation in principle can be performed with modern Finite Elements techniques. In the tsunami community what is commonly done is to work within the following simplifying hypotheses:

1. the earth is an elastic material 2. wave propagation is neglected and one computes an instantaneous equilibrium conguration.

The problem reduces to a classical dislocation problem van Zwieten et al. [2013]:

an elastic body contains a dislocation plane (fault plane) on which a discontinuous displacement instantaneously occurs. The question is which is the equilibrium conguration of the body after this discontinuous initial condition has been imposed. An analytical solution exists for a nite rectangular dislocation plane immersed in a homogeneous half space (the domain boundary consists of a at surface and a far eld boundary which is at rest; isotropy and homogeneity of the material are considered). This is the celebrated Okada solution [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF]. For earthquakes, the interesting quantity is the displacement at the at boundary which represents the seaoor. The great advantage of this approach is that we can easily implement the subfaults conguration of Satake. Due to linearity, we superpose at each point of the seaoor the solution related to each subfault displacement, as if it was acting independently from the others. In this way, the slip distribution interpolated every 30 [s] is used as initial condition for the Okada analytical expression to compute the seaoor displacement.

Free surface motion

The free surface displacement is assumed to be identical to the static vertical displacement of the seaoor, basically as if the water would translate instantaneously. This rough assumption is classically justied by several observations (see the review in [START_REF] Dutykh | The contribution of horizontal sea-bed displacements into tsunami generation processes[END_REF]): since tsunamis are long waves, the vertical dynamics can be neglected 2) elastic waves are much faster then gravity waves (in the order of ten times faster), as a consequence the seaoor 5. Tohoku-Honsu tsunami simulation deformation is assumed to be instantaneous 3) the bathymetry has in general mild slope, which means that the eect of horizontal seaoor deformation is negligible for tsunami generation.

For the 2011 Tohoku-Honsu earthquake the horizontal seaoor displacement was large near the trench axis. If large horizontal seaoor displacements occur where the bottom present a strong slope, such as in correspondence of the trench, the tsunami wave can be amplied by 30% [START_REF] Tanioka | Tsunami generation by horizontal displacement of ocean bottom[END_REF]. For this reason we have accounted also for the horizontal seaoor displacement into the computation of the waveform. A partner of the TANDEM project, the Bureau de Recherche en Géologiques et Miniérs (BRGM) provided us the free surface adding the horizontal displacement to the waveform.

Then the initial condition for the SWEs code is generated as follows. Starting from the time origin, we impose the computed seaoor displacement at t = 0 [s] as initial condition for the SWEs code. This initial condition is propagated and at every sampling time (30 [s]) the waveform related to the interpolated slip is superposed to the propagating wave. Progressing this sequence of wave propagation and superposition of free surface displacement, we obtain the waveform at t = 300 [s] which is then propagated up to the coast. The initial waveforms at t = 300 [s] are reported for completeness in gure 5.2 (with and without horizontal displacement).

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Results with tted mesh

Results with tted mesh

We present the numerical simulation of the Tohoku-Honsu tsunami with the Residual Distribution and the Finite Volume schemes developed in previous chapters. The initial waveform does not include horizontal displacement, thus it is the one of the left picture in gure 5.2. A body tted mesh was created by Electricité de France (EDF) and provided to us within the partnership of the TANDEM project. EDF made available also the results of the numerical SWEs simulation that they run with the software Telemac2D. This is a Finite Element hydrodynamics software developed at EDF and freely available online, the reader may refer to [START_REF] Hervouet | Hydrodynamics of Free Surface Flows: Modelling with the nite element method[END_REF]. The details concerning this Telemac2D tsunami simulation can be found in Legal [2017]. This will allow a comparison, for the same set-up of the Tohoku-Honsu event, of dierent SWEs codes which was one of the main objective of the TANDEM project.

The mesh has 1080181 elements and 547469 nodes with a strong variation of the local mesh size that ranges from 100000 [m] to 120 [m] in the coastal region, where the maximum resolution is necessary to correctly reproduce the complex Japan's boundary. In gure 5.3 the reader can appreciate the mesh topology. The mesh is adapted to the bathymetry variation and to the distance to boundaries. Reective boundary conditions are imposed everywhere along the coastline.

In gure 5.4 the simulated sea level displacement is compared against experimental data recorded by the GPS buoys, see gure 5.3 to check the position of the GPS buoys. The three numerical schemes provide similar free surface levels. We can see in particular that both our RD and FV implementations match well with the results produced by Telemac2D. We can underline that 5. Tohoku-Honsu tsunami simulation RD scheme is slightly less dissipative then FV, providing better resolved secondary peaks. The tsunami simulation agrees quite well with experimental data, with an accurate computation of the arrival time of the leading wave and well reproducing the complex interaction between incoming and trapped waves. In gure 5.5 and as observed by Chen et al. [2014], we see that the elevation wave, once reected, undergoes a complex transformation: reected waves travel parallel to the coast from north to south and from south to north or focus to form higher waves such as it happens in the Sendai Bay (see the red area in the last picture of 5.4). Always from gure 5.4 we clearly see that all the numerical simulations largely underestimate the rst elevation wave peak, up to 3.5 [m] error with respect to the observations at the GPS802. Aside form this drawback we can conclude that the Shallow Water model reproduce a realistic scenario for the Tohoku event. The tsunami wave was a weakly non-Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. One possible explanation for the leading wave underestimation could be the absence of the contribution of horizontal seaoor displacement in the waveform computation considered so far. Horizontal displacement has been included and its eect can be analyzed by inspecting the GPS buoys in gure 5.7. Actually the computation of the rst wave peak improves for all the buoys in the case in which horizontal displacement is retained. In particular the wave height recorded at the buoy G802 (where the highest tsunami height can be observed) is one meter higher. Despite this improvement, maximum wave heights at the buoys are still underestimated with respect to observation (up to 2.5 [m] at the GPS802).

Results with tted mesh

We have tried to improve the results with the moving mesh algorithm discussed in chapter 3 but, unfortunately, this was not immediate. First of all, the topology of the mesh in gure 5.3 is not suited for moving mesh computations due to the strong variation of the mesh size. Here, to perform a moving mesh tsunami simulation, we have used a computational grid which is rened smoothly as the boundary is approached. The minimum mesh size at the coast is h k = 1000 [m] compared to a minimal size of 120 [m] of the mesh provided by EDF, while the maximum is h k = 50000 [m]. The grid has 151650 nodes and 299207 elements. In the initial waveform it is not considered the horizontal displacement. The moving mesh parameters, see section 3.4.2, are α = 50, β = 0.1 and δ = 200.

The initial propagation of the tsunami wave and the corresponding mesh movement are visualized in gure 5.8. The images show snapshots of the initial state, after 20 and 40 . We can observe that mesh adaptation increases the resolution of wave shoaling and allows to well capture the smooth wave pattern.

However as the wave interact with the boundary mesh tangling has revealed to be more important then we expected. In gure 5.10 we show typical tangling situations at the boundaries that we have faced during the simulation. The smoothing technique discussed in section 3.4.3 become fundamental to solve this issue and avoid tangling in most situations. However to be able to perform a full simulation without mesh tangling, an a posteriori automatic untangling procedure has been set up. This procedure consist in an a-posteriori reduction of the nodal displacement by locally increasing the value of the relaxation parameter µ i in equation (3.17

Results with embedded mesh

when all the cells have been untangled:

1. we set a counter n = 0 and a constant relaxation factor µ i = 1 ∀i ∈ T h .

Then we compute the relaxation step with (3.17)

x k+1 i = x n i + δk with d G,i the distance between the center of the tangled element and node i.

3. we perform a new relaxation step (3.17):

x k+1 i = x n i + µ i δk i and we go back to step 2.

The simulated free surface level for the adaptive ALE simulation is reported in gure 5.11 where it is compared with the xed grid simulation. We observe that the moving mesh is slightly curing the excessive diusion of the leading tsunami wave encountered for the xed mesh simulation. This should suggest that the xed computation obtained with the mesh of EDF is locally not converged yet. As a nal remark we note that the overall cost of the moving mesh algorithm with the untangling procedure grows making it ineective.

For this reason we have moved to an embedded mesh that, in the benchmarks of sections 3.7.3, 3.7.4 and 3.7.5, provided an accurate wave reection when coupled with a strong mesh adaptation at the shoreline. In that cases tangling was not an issue.

Results with embedded mesh

We have tested an embedded approach in which the whole coast has simply been treated as a wet-dry boundary. We will in particular present simulations preformed with our RD approach on both a xed and a moving grid.

5. Tohoku-Honsu tsunami simulation

Results with embedded mesh

In the above section we got the evidence that, for the Tohoku-Honsu event, the horizontal displacements must be taken into account in the computation of the initial waveform. For this reason, the initial condition considered in all the tests of this section will include this eect and will coincide with the one in the right picture of gure 5.2.

With the embedded approach, friction becomes important during the inundation stage. It is well known that having accurate maps of the Manning coecient (roughness maps) has a considerable inuence on inundations, specially in densely populated areas. Here we will employ a simpler, yet widely used approach of considering a constant Manning of n = 0.03125 [sm -1/3 ]. This value represents land covered with coarse sand and is typically used when no detailed roughness map is available [START_REF] Gayer | Unconditionally stable spacetime discontinuous residual distribution for shallow-water ows[END_REF].

First we have run a xed grid computation with an unstructured embedded mesh composed of 689861 nodes and 1379054 triangles. The mesh is adapted to the bathymetry: the minimum local mesh size is at the coast (120 [m]) and the maximum mesh size is 5 [km]. We note that this grid is ner with respect to the tted one detailed in the previous section and it is used here to run a reference simulation. In gure 5.12 the reader can nd the simulated free surface level at the buoys. As suspected from the previous paragraph, we realize that the more rened mesh reduce numerical diusion and the wave prole of the leading wave appear much sharper. More importantly, the height of the rst wave is considerably closer to the experimental value with respect to the result obtained with the tted mesh of the previous section. A side eect is that, due to the friction in the wetting/drying phase, the reected wave are slightly damped.

With the embedded reference mesh we model also the ooding of the Japan coastal prefectures. We focus the attention on three bays in the south of the Iwate prefecture: in gure 5.13, as they appear from south to north, they are the Kesennuma, Hirota and Ofunato bay. In gure 5.14 we provide a comparison between the inundated areas predicted by the RD scheme and the surveyed areas for these three densely populated bays. With a mesh size of 120 meter we can say to have a good prediction of the runup process, except very ne scale ooding such as the inundation of two narrow channels at the end of Hirota bay. In the simulation they remain dry even if water was observed.

We test now if the moving mesh methods could provide the same resolution of the reference simulation, using a coarser computational mesh. The computational mesh is adapted to the bathymetry but the minimum local mesh size (in the proximity o the coast) is 360 [m], three times larger with respect to the reference mesh. It results a computational mesh composed by 364864 nodes and 728874 elements. The moving mesh parameters are α = 50, β = 0.1 and δ = 200. of the free surface level and the corresponding adapted meshes are plotted every twenty minutes at t = 00 (g. 5.15 and 5.16), t = 20 ( (g. 5.17 and 5.18) and t = 40 (g. 5.19 and 5.20). We observe that the mesh points move and cluster around the complex tsunami wave pattern. This becomes visible especially when the wave shoals and the mesh points gather in correspondence of the largest gradients of the incoming waves train. The shoreline is also detected by the mesh. We recall that tangling at the coast was the main limitation of the tted method. In this embedded moving mesh computation tangling has not occurred even without employing the smoothing step and the a-posteriori relaxation implemented in the previous section. In gure 5.21 we can observe a closer zoom of a very complex coastline detected by the moving mesh method without tangling. Finally from gure 5.22 we can see that the simulated sea level computed with the moving mesh at the wave gauges is similar to the one computed with the reference. Mesh adaptation at the shoreline allows to reproduce sharply the complex reection of the tsunami wave.

Results with embedded mesh

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Results with embedded mesh

Conclusions

Summery of the thesis results

In this thesis we have analyzed the use of r -adaptation for the shallow water simulation of complex wave interactions and wave runup on irregular bathymetries.

In chapter 2 we have provided a thorough theoretical setting to construct Well Balanced Finite Volume and Residual Distribution schemes on moving grids with an ALE and a rezoning approach. This schemes includes a mass conserving correction of the nodal bathymetric heights, based on a quadrature of the given bathymetric data. Being based on the actual data, this correction requires no re-initialization, as e.g. the ALE remap used in Zhou et al. [2013a] which leads to a numerical deviation from the real bathymetric data. Both simple academic test cases and long wave benchmarks have conrmed our theoretical expectations in terms of the conservation of steady equilibria and mass conservation.

In chapter 3 we have coupled these schemes with the Laplacian-type radaptation method of [START_REF] Ceniceros | An ecient dynamically adaptive mesh for potentially singular solutions[END_REF] and we have investigated different coupling strategies between the ow solver and the adaptation in terms of accuracy and cost (CPU time). The delicate point is here the overhead of the mesh adaptation method when the ow solver is based on fully explicit multi-stage methods. Our results show that, as long as possible, one should stick to the use of a fully ALE method coupled with the mesh deformation solver, used with a simplied solution remapping for the error sensor. This turns out to be the most ecient in terms of accuracy for a given CPU time, as well as the most robust in providing substantial improvements both for smooth and non-smooth features, including an improved prediction of runup.

We have also proposed a simplied rezoning method which allows to run the ow solver on a xed mesh. The method proposed allows to save signicant CPU time and can be used in situations where local remeshing is necessary, and a full ALE method with nite time step values cannot be used. These results improve on, and complete the studies done in the past in e.g [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF]; [START_REF] Ceniceros | An ecient dynamically adaptive mesh for potentially singular solutions[END_REF]; [START_REF] Ni | Remapping-free ale-type kinetic method for ow computations[END_REF]; Zhou et al. [2013a,b]; [START_REF] Chen | Second-order accurate godunov scheme for multicomponent ows on moving triangular meshes[END_REF] providing quantitative as well as qualitative elements.

In chapter 4 we have extended the Shallow Water solver (both Residual Distribution and Finite Volume) and the moving mesh method to curvilinear coordinates. The Residual Distribution approximation enhances to unsteady problems the scheme of [START_REF] Rossmanith | Residual distribution schemes for hyperbolic balance laws in generalized coodinates[END_REF] and it is appealing due to its simplicity and rigor in treating all the issues that the spherical geometry poses such as the medium heterogeneity and the geometric source term. Numerical results show that the resulting ALE r-adaptive method can improve the resolution of steep solutions of the SWEs on the sphere.

In chapter 5 we have employed r-adaptation to the simulation of the Tohoku-Honsu tsunami. This includes the large scale wave propagation and coastal ooding. We show that with a combined use of r-adaptation an embedded meshes it is possible to perform accurate tsunami simulations with a limited number of mesh points.

Limitations and future perspective

Future developments will involve the extension of our analysis to both multistep and higher order methods, as well as the addition of dispersive eects, based on the approach of Filippini et al. [2016]. The reduction of the adaptation overhead obtained with the full ALE approach (and with the simplied rezoning) also opens the door to new developments. In particular, both the underlying mesh adaptation method, as well as its discretization, and iterative solution will be object of future work. Although successful we have observed that the implementation of the simple Laplacian poses some strong limitation on the control of mesh quality. A possible avenue is the combination of the simple Laplacian approach used here, providing a very sharp approximation of discontinuous features, with elastic deformation [START_REF] Stein | Automatic mesh update with the solid-extension mesh moving technique[END_REF] allowing greater control on mesh quality. Concerning the discretization and resolution of the MMPDE, improved iterative methods can certainly be benecial to produce grids with improved quality. In [START_REF] Tang | Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws[END_REF] it is mentioned that the algebraic system is solved with Gauss-Seidel iterations, but no quantitative informations whatsoever are given w.r.t. the cost overhead, or of the total cost of the adaptive simulations compared to a xed ne mesh one. Future work will be also dedicated to overcome the pole singularity and to propose a Residual Distribution method for global circulation model. The pole problem for the moving mesh method is more complicated and necessitates further study. However the work of [START_REF] Di | Moving mesh methods for singular problems on a sphere using perturbed harmonic mappings[END_REF] could be implemented to then test the ALE-SWEs moving mesh method against benchmarks on the whole sphere.

Appendix A Tensors, vectors and nabla operators in curvilinear coordinates Notation for curvilinear coordinates Consider a transformation from a Cartesian coordinate system x = {x 1 , x 2 , x 3 } to a general curvilinear coordinates system z = {z 1 , z 2 , z 3 } x = G(z)

The orthonormal basis in Cartesian coordinates is {e i } = {e 1 , e 2 , e 3 }. The non-orthonormal covariant basis vectors is dened by {g j } = {g 1 , g 2 , g 3 }. A vector x in Euclidean space can be expressed in curvilinear coordinates as x = g 1 z 1 + g 2 z 2 + g 3 z 3 = 3 i=1 g i z i ≡ g i z i

where the symbol ≡ means that we have used the standard summation convention. The basis vectors of the curvilinear framework can be calculated from the orthonormal basis associated to the Cartesian coordinate system The components ∂x j ∂z i fully describe the transformation from cartesian to curvi- linear coordinate system and they can be collected in a transformation tensor, or jacobian

J =     ∂x 1 ∂z 1 ∂x 1 ∂z 2 ∂x 1 ∂z 3 ∂x 2 ∂z 1 ∂x 2 ∂z 2 ∂x 2 ∂z 3 ∂x 3 ∂z 1 ∂x 3 ∂z 2 ∂x 3 ∂z 3     (A.2)
which admits a determinant det(J ) = 0. The contravariant basis is useful to describe tangent plane to a given surface. Given a point P embedded in cartesian space, the contravariant basis g k is perpendicular to the plane dened by the covariant basis g i and g j in P . It is dened as

g k = g i × g j (g i × g j ) • g k
The metric tensor G is composed of the coecients given by the scalar product between the covariant basis vector G ij = g i • g j with determinant det(G) = 0. The inverse metric tensor G -1 is:

G ij = g i • g j
An important relation between the metric tensor and the jacobian G = J T J .

In the following we will dene the square root of the determinant of the metric tensor with √ G which will mean:

√ G = √ det G = det J (A.3)

Physical components of tensors

Now, let v be a vector in the curvilinear framework, it can be written in the covariant basis as v = v j g j , always using standard summation convention. For our purposes it is useful also to dene physical components of a vector in a unitary covariant basis. This means that the curvilinear basis vectors have to be normalized

g * i = g i |g i |
The vector reads v = (v j |g j |)g * j = v * j g * A. Tensors, vectors and nabla operators in curvilinear coordinates

Christoel symbol

The 3 × 3 second-kind Christoel symbol reads

Γ k ij = g k • ∂g i ∂z j
We recall two important property of the Christoel symbol which we will use later in the text. The rst one reads as follows

Γ i ij = 1 √ G ∂ √ G ∂z j (A.6)
the second important relationship is known as Ricci's Lemma and concerns the derivative of the coecients of the metric covariant coecient G ij with respect to z k . For us it will be useful the derivation with respect to z j . In this particular case:

∂G ij ∂z j = -G mj Γ i jm -G im Γ j jm (A.7)
for the demonstration, see [START_REF] Nguyen-Shafer | Tensor Analysis and Elementary Dierential Geometry for Physisict and Engineers[END_REF].

Nabla operator in curvilinear coordinates

Nabla operator in a general curvilinear coordinates reads

∇ = g i ∂ ∂z i (A.8)
Given a scalar φ, a vector v and a second order tensor T , we list the nabla operator in curvilinear coordinates that we use in this thesis. The divergence of a vector v in curvilinear coordinates writes

∇ • v = g i ∂ ∂z i • v = 1 √ G ∂ ∂z i √ G v i , i = 1, 2, 3 (A.9)
The divergence of a second order tensor is

∇ • T = g k ∂ ∂z k • T = 1 √ G ∂ ∂z i √ G T ij + T ik Γ j ik g j , i, j = 1, 1, 2 (A.10)
The gradient of a vector is ∇v = g i ∂ ∂z i v = g j g i ∂v j ∂z i , i = 1, 2, 3

(A.11)
Finally the Laplacian of a scalar (also called Laplace-Beltrami operator) is reported here with two dierent formulas, both useful:

∇ 2 φ = g i ∂ ∂z i • g j ∂φ ∂z j = G ij ∂ 2 φ ∂z i ∂z j -Γ k ij ∂φ ∂z k (A.12) = 1 √ G ∂ ∂z i √ GG ij ∂φ ∂z j (A.13)
Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

and Christofell symbols are

Γ 1 = g 1 • ∂g i ∂z j = 0 -tan λ -tan λ 0 , Γ 2 = g 2 • ∂g i ∂z j =
sin λ cos λ 0 0 0 of the stabilization bubble: for this part we use a non conservative form, see [START_REF] Arpaia | An ale formulation for explicit runge-kutta residual distribution[END_REF], and secondly we use the step shifting operator ∆ k instead of the original ∆ k . We mention that this approximation is necessary to end up with a fully explicit scheme and does not spoil the formal order of accuracy.

The details can be found in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF] k

= 1 ⇒ ∆ 1 = [•] n -[•] n k = 2 ⇒ ∆ 2 = [•] * -[•] n
The stabilization part is

∆t K∈D i      K γ i ∆ k w ∆t dx V K + K γ i ∂F j ∂x j -σ j ∂w ∂x j k dx V I K = K γ i S k dx V II K     
Now we can consider the sum of I and III:

I + III = ∆ k Ω i ϕ i w dx -∆t Ω i ϕ i w k ∂σ j ∂x j dx = K∈D i (I K + III K )
where we have decomposed the integral over neighbor elements. Using the fact that w k = ∆k w/2 + w n , we can write on each element

I K + III K = ∆ k j∈K K ϕ i ϕ j dxw j -∆t j∈K K ϕ i ϕ j ∂σ j ∂x j dx ∆ k w j 2 + w n j
Moreover we can use the denition of constant divergence over the element, eq.(4.43)

K ϕ i ϕ j ∂σ j ∂x j dx = K ϕ i ϕ j ∇ • σ dx = K ϕ i ϕ j dx ∆|K| ∆t|K|
Wee see the Galerkin mass matrix m G ij = K ϕ i ϕ j dx.

I K + III K = ∆ k j∈K m G ij w j - j∈K m G ij ∆ k w j 2 + w n j ∆|K| |K|

D. Proof of the eRK2-RD-ALE update

We follow exactly [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF] and we perform mass lumping that is j∈K m G ij = K ϕ i dx = |K|/3 so that

I K + III K = 1 3 ∆ k (|K|w i ) - 1 3 ∆|K|w n i - j∈K 1 2 m G ij ∆|K| |K| ∆ k w j V III K Remembering that |C i | = K∈D i |K| 3 I + III = ∆ k (|C i |w i ) -∆ k |C i |w n j - K∈D i j∈K 1 2 m G ij ∆|K| |K| ∆ k w j = |C n+1 i |∆ k w i - K∈D i j∈K 1 2 m G ij ∆|K| |K| ∆ k w j
for each element we rewrite

V K = IX K + X K K γ i ∆ k w dx = K w i ∆ k w dx IX K - K ϕ i ∆ k w dx X K
and then we sum

II K + IV K + V K + V I K + V II K into the residual of element K Φ K(k) i = w i ∆ k w ∆t dx + w i ∂F j ∂x j -σ j ∂w ∂x j k dx + K S k dx
we remark that we have computed the uctuation at the k as it has been dened in (4.7). We miss X K and V III K but they can be summed up together:

-

j∈K m G ij ∆ k w j - K∈D i j∈K 1 2 m G ij ∆|K| |K| ∆ k w j = - j∈K 1 + 1 2 ∆|K| |K| m G ij ∆ k w j
Lumping again this matrix:

-

j∈K 1 + 1 2 ∆|K| |K| |K| 3 ∆ k w j = -|C n+1 i |∆ k w i
We dene ∆ k -∆ k = ∆ k , the whole algorithm can be recast in compact form as

∆ k w i = - ∆t |C n+1 i | K∈D i Φ K(k) i (D.1)
We obtain the two stage, for k = 1 the predictor (2.55) and for k = 2 the corrector (2.56).

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Appendix E

Manifactured source term

We construct a source term to run test case #1 from [START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF].

We impose in equation (4.2) the following advection bell travelling over a background irrotatational velocity eld a = a j g j h(z, t) = f (z -at) a 1 = 1

(E.1)
a 2 = 0 called "manifactured" solution. We get a residual:

∂f ∂t + 1 √ G ∂ ∂z j √ G f a j = S h ∂f a i ∂t + 1 √ G ∂ ∂z j √ G f aa ij + G ij gf ∂f ∂z j + Γ i jk f aa jk = S hu
The continuity equation is satied exactly:

∂f ∂t + ∂ ∂z j f a j + f a j √ G ∂ ∂z j √ G = 0
Since ∂G ∂z 1 = 0 we have ∂f ∂t + ∂ ∂z j f a j = S h = 0

which is an advection equation exactly satied by the linear wave solution (E.1). For momentum equation

a i ∂f ∂t + ∂ ∂z j f aa ij + f aa ij √ G ∂ ∂z j √ G + G ii gf ∂f ∂z i + Γ i km f aa km = S i a i ∂f ∂t + ∂ ∂z j f a j + G ii gf ∂f ∂z i + Γ i km f aa km = S i
The terms which are left represents a source term for the Shallow Water system S i = G ii gf ∂f ∂z i + Γ i km f aa km

These extra terms provide that (E.1) is an exact solution of the system. The term gf ∇f balances the nonlinearity nullifying the hydrostatic part and the term Γ km f aa km balances the geometrical term. The modied shallow water system writes

∂h ∂t + 1 √ G ∂ ∂z j √ G hu j = 0 ∂hu i ∂t + 1 √ G ∂ ∂z j √ G huu ij + G ij gh ∂η ∂z j + Γ i jk huu jk = S i

√

  gh being the local wave celerity. The nonlinear eigenstructure gives more insight about the physics of waves in the shallow water regime. Shallow water waves are propagating at speed ± √ gh relative to the uid, thus wave in deeper water moves faster. If the amplitude of the wave is not small compared 1. Shallow Water equations: continuous properties and Well-Balanced approximations on unstructured grids

  As we can see, the bathymetric term disappears when the bottom is at least piecewise continuous b R = b L . On the contrary, in presence of a bottom jump b R = b L , this term play an important role in Rankine-Hugoniot relationship and in the resolution of the Riemann problem. It has been shown that a standing wave in correspondence of the bottom jump appears, see Bernetti et al. [2008], Alcrudo and Benkhaldoun [2001]. Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Figure 1

 1 Figure 1.1: Left: Nodal normal denition. Center: Finite Volume dual cell. Right: dual cell interface with normals.

Figure 1

 1 Figure 1.2: The abstract steps to construct a Residual Distribution approximation.

Figure 2

 2 Figure 2.2: A sketch of Eulerian, Lagrangian and ALE congurations.

Figure 2 . 3 :

 23 Figure 2.3: Computing the bathymetry nodal values with the standard FV method. pointwise value(left) and the method proposed here (right). Top) quadrature points. Middle) 1D sketch of bathymetry values that are computed from quadrature rule (2.34). Bottom) 1D sketch of analytical and computed bathymetry integral.

2. 8 .

 8 Figure 2.4: Preservation of the lake at rest on wet/dry cells. Left: a node is passing from wet to dry. Right. a node is passing from dry to wet.

  For example, the bathymetry b(x) or the friction coecient c F (x), that also could vary in space. In the moving mesh case one should transport all this information on the new meshes. This means a time consuming computation of b n+1 /c n+1 F at each timestep, trough an accurate interpolation or through the remap equation (2.15). For an implementation of the ALE-FV for in hydrodynamics, the reader can refer to the cell centered FV on moving mesh developed Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. 2.10. ALE Finite Volume for SWEs in Zhou et al. [2013a], which is based on a xed grid ow evolution but with an ALE remap for the interpolation step. Here we consider the standard well balanced node-centered FV-MUSCL scheme introduced in chapter 1 and we

Figure 2

 2 Figure 2.5: DGCL for ALE-FV. Area swept by the interface C ij .

  as we see it is enhanced by the ALE part. Similarly, in every cell K, φ b i denes a splitting of the ALE residual in b:φ b = -K σ j ∂b ∂x j dxIn 2.63, two parts can be put in evidence (within brackets). The last three terms underlined, named with II, should correspond to the approximation of the ALE remap equation for the bathymetry. They do not sum up to zero because we have deliberately choose to not compute b n+1 i with the ALE remap equation which means by:

j 2 .Figure 2

 22 Figure 2.7: Linear Advection. Left: Lake at rest for the NO WB ALE formulation and failing in verifying Well Balanced. Middle: comparison between the numerical solution and exact one on the symmetry line x = 0.5. Right: convergence order for the L 2 -norm of the error.

  2.12. Numerical experiments with b = b(x, y) a given function, and with the ux Jacobian a = ∂F (u)/∂u. The following denition of the uxes will be used, F (u) = au, with a = [ 1 0 ]. Introducing the variable η = u + b, equation (2.65) admits a non-trivial steady state given by u(x, y) + b(x, y) = η 0 = const (2.66)

T 2 .

 2 Well-Balanced conservative methods for the Shallow Water equations in ALE form

Figure 2

 2 Figure 2.8: Vortex advection. 2nd level of mesh renement: original mesh and transformed mesh. Mesh convergence for: left) RD. right) FV.

Figure 2

 2 Figure 2.9: Dambreak with circular hump. Left: bathymetry. Right: dimensionless mass error for dierent quadrature formula of the bathymetry integral.

Figure 3 . 1 :

 31 Figure 3.1: An example taken from Alauzet [2010]: a uniform mesh in a metric space and a non uniform mesh in Euclidean space.

Figure 3 . 2 :

 32 Figure 3.2: Examples of non convex physical domain in hydrodynamic simulation.

  (3.13) and f is a function which is constant everywhere except in the narrow region Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

  Figure 3.3: Example of mesh adaptation to the wet/dry interface.

Figure 3

 3 Figure 3.4: Tangling for non convex domain. right) non smoothed mesh left) smoothed mesh

Figure 3

 3 Figure 3.5: Top) test (3.18). Bottom) test (3.19)

Figure 3

 3 Figure 3.6: Top) test (3.20). Bottom) test (3.21).

h

  according to the scalar version of FV/RD projections, (3.23) or (3.24). ENDDO Step 4. Let x n+1 = x kmax+1 and T n+1 h = T kmax+1 h . Evolve the underlying balance law in ALE framework with the FV/RD-eRK2 scheme, see Eq. (2.38) or Eq. (2.56) on the grid T

Figure 3

 3 Figure 3.7: Moving mesh algorithms. Top) ALE. Middle) Rezoning (EUL1).

  scheme, see Eq. (2.38) and Eq. (2.56) with σ = 0, on the grid T n+1 h .

h

  Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.3.6. Adaptive algorithmson the nal updated mesh in order to give a proper initial condition for the ow solver. If what we thought about the simplication of the ALE remap would be conrmed by numerical experiments, then we could use a simplied version of the remap algorithm to project only the free surface variable for then estimating the error. A full high resolution remap is used only after reaching k = kmax in the adaptation loop, to perform the interpolation. For example for FV:

otherwise

  Figure 3.8: Rotation test. left) adapted mesh for the second level mesh. center) adaptive FV. right) adaptive RD.

Figure 3 . 9 :Figure 3

 393 Figure 3.9: Rotation test with RD. left) order of convergence. right) Error vs CPU time.

Figure 3 .Figure 3

 33 Figure 3.12: Burgers Equation computed with MUSCL scheme on xed grid. Left: coarse mesh h K = 1/40. Right: ne mesh h K = 1/100.

Figure 3 Figure 3 Figure 3 Figure 3 Figure 3 Figure 3 Figure 3

 3333333 Figure 3.14: Burgers equation computed with ALE-RD with rst order projection scheme (LxF) to compute the error estimate in the MMPDE. Top) adapted grid at nal time, 20 equispaced solution isolines between 0 and 1. Bottom) comparison of the solution along the symmetry lines and the lines at y = 0.4

Figure 3 .

 3 Figure 3.21: Asymmetric dam-break computed with RD scheme. 30 equispaced contour lines for h and adapted mesh.

Figure 3 .Figure 3 Figure 3

 333 Figure 3.22: Asymmetric dam-break computed with FV scheme. 30 equispaced contour lines for h and adapted mesh.

Figure 3 .

 3 Figure 3.25: Double symmetric dam break, t = 10, t = 35, t = 50. Solution iso-lines and mesh computed with adaptive RD scheme.

  Figure 3.26: Small perturbation of a lake at rest (RD scheme). Solution isolines at t = 0.24, t = 0.48 are shown for xed grid and adaptive computations. Top: xed coarse grid. Middle: xed ne grid. Bottom: adaptive ALE scheme.

Figure 3 .Figure 3

 33 Figure 3.27: Small perturbation of a lake at rest (FV scheme). Solution isolines for t = 0.24, t = 0.48 are shown for xed grid and adaptive computations. Top: xed coarse grid. Middle: xed ne grid. Bottom: adaptive ALE scheme.

Figure 3 Figure 3

 33 Figure 3.29: Small perturbation of a lake at rest (FV scheme). Solution at t = 0.48 along line y = 0.5. Left: ALE. Right: EUL2.

Figure 3 .

 3 Figure 3.31: Conical island. Left: sketch of the computational domain with gauges. Middle: static coarse mesh topology. Right: static ne mesh topology.

Figure 3 .

 3 Figure 3.32: Conical Island computed with xed ne mesh: total water height η signal registered at the gauges g9, g16, g22 and comparison with experimental data.

Figure 3 .

 3 Figure 3.33: Conical Island: contour lines for total water height η and adapted mesh at dierent time instants, t = 6.0, 8.0, 10.0 [s]

Figure 3 .Figure 3

 33 Figure 3.36: Okushiri experiment. Contour lines for h and mesh at t = 16.5 [s]. Top: xed coarse grid. Middle: xed ne grid. Bottom: adaptive ALE scheme.

Figure 3 Figure 3

 33 Figure 3.38: Okushiri experiment computed with adaptive ALE schemes. Top row: maximum runup for RD (left) and FV (right). Bottom row: total water height η signal registered at the gauge placed in the Monai valley. RD (left) and FV (right).

Figure 3 .

 3 Figure3.41 shows visualizations of the wave patterns arising from this complex interaction. In the gure, the top row shows the results obtained on the coarse grid. The second row reports the results on the ne grid. The ALE results, and the corresponding grids, are reported in the third and fourth row. Figure3.42 shows visualizations comparing the ALE results (top half of the pictures, with snapshots of the video of the experiment, available online [https : //www.youtube.com/watch?v = I4uT HW BpaZg]. The results are those obtained with the RD scheme, but very close ones are obtained with the FV method, not reported here due to shorten the presentation. The ALE results on the coarse mesh provide a ow description which is even clearer of the one obtained on the ne mesh, and clearly allows to resolve wave and vortical structures otherwise absent on the xed coarse grid simulations. The comparison with the experimental snapshots shows a very satisfactory qualitative agreement with the patterns observed in the wave tank. Finally, gure 3.43 provides the time series in gauge ADV3. We can see that the adaptive simulation computes better resolved proles of the waves reected from the bar. The gain in time is between 40-50% with respect to using a reference mesh.

Figure 3 .

 3 Figure 3.40: Solitary wave on a shelf: 3D visulization of the bathymetry, the scale of the z-axis is 5:1 with respect to x and y-axis

Figure 3

 3 Figure 3.41: Solitary wave on a shelf (RD scheme). Solution isolines at t = 12, 17 and 25 [s] are shown for xed grid and adaptive computations. First row: xed coarse grid. Second row: xed ne grid. Third and fourth: adaptive ALE scheme.

Figure 4

 4 Figure 4.1: Example of a SWEs tsunami simulation of the 2011 Tohoku-Honsu event. In red) Cartesian SWEs with Mercator coordinates. In blue) spherical SWEs with lat-lon coordinates.

  3)(4.4)(4.5) is a non-homogeneous hyperbolic system of PDEs.Following[START_REF] Leveque | Wave propagation algorithm for multidimensional hyperbolic systems[END_REF], we remark that the ux vector depends explicitly Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami. 123 4.2. ALE in curvilinear coordinates on the position making the system non-autonomous. If one multiplies (4.3) by

  Figure 4.2: A sketch of Eulerian, Lagrangian and ALE congurations on a sphere.

  The initial or reference conguration is denoted by B 0 and the position of the point Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

  0 is the covariant vector basis. As in chapter 2 lower case letters label the actual conguration and Greek letters are used for the Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

  (4.20) Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

  (4.11) for the bathymetric function (does not depends explicitly on time b = b(z)) ∂b ∂t ζ = ∂b ∂z j σ j multiplying by J A √ G and summing it to the (4.15) (multiplyed by b) we have the ALE remap in curvilinear coordinates

Figure 4

 4 Figure 4.3: Cartesian grid of Rossmanith et al. [2004] and quadrature points for numerical approximation of the geometric source term

  4. ALE r-adaptive methods for the Shallaw Water equations on the sphere

Figure 4 . 4 :

 44 Figure 4.4: Circular hump propagating on a sphere. Comparison between two dierent treatments of the geometric source term on the xed unstructured grid in the bottom-right picture. Dashed-blue) extension to unstructured grids of the centered approximation of Ullrich et al. [2010]. Continuous-black) extension to unstructured grids of the wave propagation algorithm Rossmanith et al.

  free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.4.6. Finite Volume for SWEs on the sphereWell Balanced for moving mesh on the sphere The discretization of the topography/bathymetry term follows what done for plane Cartesian coordinates (see sections 1.6.2 and 2.10.1) which in turns consisted in the upwind discretization of[START_REF] Hubbard | Flux dierence splitting and the balancing of source terms and ux gradients[END_REF];[START_REF] Nikolos | An unstructured node-centered nite volume scheme for shallow water ows with wet/dry fronts over complex topography[END_REF]. We distinguish two contributions S b ij = S c ij + S * ij . We introduce the nodal dierences ∆h

  free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

  4. ALE r-adaptive methods for the Shallaw Water equations on the sphere and then integrate exactly with respect to linear variation of h and b. One can show the following result: Proposition 3. The explicit predictor corrector residual distribution prototype (4.40), (4.42) with edge velocities (2.25) veries the Well-Balanced property both on moving and xed meshes, provided that the same linear piecewise continuous approximation is used for w, b and consequently for u.

  λ) = 50 exp(-250(ψ 2 + λ 2 ) (4.55) with ψ = z 1 and λ = z 2 . The iteration is repeated in the pseudo-time loop until convergence is reached. For the rst smooth example reported on the left of gure 3.5 we have set α = 100. In the left picture of gure 4.6 we show that also singularities are well handled with a proper choice of the parameter. 4.9 Numerical experiments 4.9.1 Global steady state zonal geostrophic ow Case #2 from Williamson test suite[START_REF] Williamson | A standard test set for numerical approximations to the shallow water equations in spherical geometry[END_REF] is a steady state solution of (4.2) which allow to measure the order of accuracy of numerical Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Figure 4

 4 Figure 4.5: Top) test (4.54). Bottom) test (4.55)

Figure 4 . 9 :

 49 Figure 4.9: Advection of cosine bell. Left: mesh topology. Right: uid depth isolines.

Figure 4 .

 4 Figure 4.10: Advection of cosine bell. 2nd level of mesh renement: original mesh and transformed meshes at two dierent instants.

  Therefore, the solution should remain symmetric for all time. We will use two Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

Figure 4 .

 4 Figure 4.11: Advection of cosine bell. Left: convergence for RD. Right: convergence for FV.

Figure 4 .

 4 Figure4.12: Circular hump on a sphere. Snapshots of FV depth isolines (20 levels between 0 and 0.55) at dierent times, t = 0.3, 0.6, 0.9. Top) coarse mesh. Middle: ne mesh. Bottom: adaptive ALE with adaptive mesh.

Figure 4

 4 Figure 4.13: Circular hump on a sphere. Solution along a symmetry line.Comparison between x mesh simulations and adaptive ALE.

Figure 4 .

 4 Figure 4.14: Circular hump on a sphere. Snapshots of FV depth isolines (20 levels between 0 and 0.55) at dierent times, t = 0.4, 0.8, 1.2. top) coarse mesh. Middle: ne mesh. Bottom: adaptive ALE.

  Figure 5.1: Left) Fault geometry and nomenclature. Right) Approximation of the fault plane in subfaults and approximated slip eld, (to generate this image we have used the software SlipReal Mai and Beroza [2002])

Figure 5

 5 Figure 5.2: Waveforms in [m] at t = 300 [s] used as initial condition for the SWEs code and computed accounting for: left) only vertical seaoor displacement Right) both horizontal and vertical seaoor displacement.

Figure 5

 5 Figure 5.3: Left) GPS buoy position. Right) Mesh topology. zoom of the Japan coast and zoom of the Sendai bay.

Figure 5 . 4 :

 54 Figure 5.4: Comparisons of the simulated and observed sea level displacement recorded at GPS buoys. The simulated curves refer to FV, RD and Telemac2D.

Figure 5 . 5 :

 55 Figure 5.5: Snapshot of the simulated oshore (RD) free surface level in meters every one hour starting from time 14:51:18: t = 0:00, 1:00, 2:00 and 3:00.

Figure 5

 5 Figure 5.7: Simulated (RD) free surface level recorded at GPS buoys with and without seaoor horizontal displacement in the waveform computation.

i 2 .

 2 ∀K ∈ T h we check tangling. If there is no tangling the loops ends. Otherwise, if one cell tangles |K| < 0 we ag it and we compute its center position x K G , then we set n = n + 1 τ = 2n 2 and we compute for every node ∀i ∈ T h µ i = min µ i , 1 -e -

Figure 5 . 8 :

 58 Figure 5.8: Snapshot of free surface level (in meters) and corresponding moving mesh renement every twenty minutes starting from time 14:51:18: t = 0:00, 0:20. and 0:40.

Figure 5 . 9 :

 59 Figure 5.9: Snapshot of the simulated (ALE-RD) free surface level (in meters)and corresponding moving mesh renement every twenty minutes starting from time 14:51:18: t = 0:40, 1:00 and 1:20.

Figure 5 .

 5 Figure 5.10: Examples of mesh tangling at the boundaries.

Figure 5 .

 5 Figure 5.11: Simulated (RD) free surface level recorded at GPS buoys with xed tted mesh and ALE moving tted mesh.

From

  gure 5.15 to 5.20 we show the solution and the mesh, in an alternate fashion, for a zoom of the physical domain close to Japan coast. The snapshots 5. Tohoku-Honsu tsunami simulation

Figure 5 .

 5 Figure5.12: Simulated (RD) sea level displacement recorded at GPS buoys with an embedded reference mesh.

Figure 5 .Figure 5

 55 Figure 5.13: Simulated ooding in the southern Iwate prefecture. Snapshot of the numerical simulation (RD) at t = 0:00, 0:40 and 0:50. The picture in bottom-right position is the embedded mesh with the coastline drawn in green.

Figure 5 .

 5 Figure 5.15: Snapshot of the free surface level (same color legend of previous gures) starting from time 14:51:18 at t = 0:00.

Figure 5 .

 5 Figure 5.22: Simulated sea level displacement recorded at GPS buoys with xed embedded mesh and moving embedded mesh.

j

  using the following components v * j = |g j |v j (A.4) the star denotes physical contravariant vector components dened on the unitary basis. For tensors T = T ij g i g j the physical contravariant tensor components in the unitary contravariant basis writeT * ij = |g i ||g j |T ij (A.5)

  here in the scalar notation) leads to the concept of Linearity

	Preserving schemes. RD schemes are LP if the distribution coecients are
	uniformly bounded with respect to the solution and data of the problem, hence
	exists a constant C such that

  property) and ALE discretizations can provide exact integral variants of H 2 = 0 (DGCL). However, unfortunately Eulerian methods are unable to embed exact integral (or local) forms of the advection equation

H 1 . So, in correspondence of steady equilibria, these methods will always have a truncation associated to the term H 1 . On the other hand, at the continuous level we ca use (2.14) to deduce that

  .10. ALE Finite Volume for SWEs(1.22), (1.23), and (2.43) verify the DGCL for constant b, and the Well-Balanced both on moving and xed meshes, provided that the same reconstruction procedure is used for u and b.

	43)
	since the ALE Jacobian has an extra term due to mesh motion. With these
	denitions we have now the following characterization.

Proposition 1. The nite volume discrete equations (2.37)-(2.38) with numerical ux and source respectively given by (2.40)(1.21) and with (4.32), Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

2

Table 3

 3

.4: Conical island. CPU times.

Table 3 .

 3 5: Okushiri experiment. CPU times.

Table 3 .

 3 6: Solitary wave on a shelf. CPU times.

Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.

  emerges, where w i is the Petrov-Galerkin test function. Within a RD context, it is important to build mass matrices which are consistent with the space approximation, thus with a given set of distribution coecients. Assigned a Adaptive techniques for free surface ow simulations. Application to the study of the 2011 Tohoku tsunami.4.7. Residual Distribution for SWEs on the sphere certain distribution {β K i } i∈K , this consistency property is translated in the following constraint on w i Ricchiuto and Abgrall [2010]:

Luca Arpaia

1.6. Discrete approximation

3.7. r-adaptation for Shallow Water ows

h i G i +h j G j 2 √ h i u i √ G i + √ h i u j √ G j

4.9. Numerical experiments

Acknowledgments v vi Luca Arpaia

Chapter 5

Tohoku-Honsu tsunami simulation

The 11 March 2011 a catastrophic earthquake of magnitude M w 9.0 ruptured the plate o the coast of the Tohoku region. This is usually referred to as the 2011 Tohoku-Honsu earthquake. Many authors that have studied the fault rupture mechanism agree that the earthquake was compact with huge plate slips concentrated near the trench axis, see [START_REF] Fujii | Tsunami source of the 2011 o the pacic coast of tohoku earthquake[END_REF]; [START_REF] Shao | Focal mechanism and slip history of the 2011 m9.1 o the pacic coast of tohoku earthquake, constrained with teleseismic body and surface waves[END_REF]; [START_REF] Ammon | A rupture model of the 2011 o the pacic coast of tohoku earthquake[END_REF]; [START_REF] Iinuma | Coseismic slip distribution of the 2011 o the pacic coast of tohoku earthquake (m9.0) estimated based on gps data -was the asperity in miyagi-oki ruptured?[END_REF]; Satake et al. [2013]. We will give a denition of the slip vector in 5.1.1, for now see gure 5.1 for an illustration of its physical meaning. Maximum slips may vary from one author to the other but they stay in the range 30 -70 [m]. The fact that large slips occurred near the trench, at rather low depths, produced a large seaoor uplift that in turn triggered a devastating tsunami wave. At the end, signicant damages occurred in the coastal prefectures of Fukushima, Myiagy and Iwate and they were mostly caused not by the earthquake itself, but by abnormal wave heights (up to 40 [m] in northern Japan) that overwhelmed protection barriers, breakwaters and other sea defenses.

Very interestingly for us, the Tohoku-Honsu earthquake was the largest instrumentally recorded earthquake in Japan history. A part from extensive ground motion dataset, the tsunami was recorder by 10 nearshore GPS gauges, 32 tidal gauges and by the oshore DART buoy system. All these data are available online: for the GPS wave gauges data one can access the website of the Nationwide Ocean Wave information network for Ports and Harbours (NOWPHAS) NOWPHAS [2017]; the DART buoy data are available from the website of the National Oceanic and Atmospheric Administration (NOAA) NOAA [2017a]. Equally important was the post-tsunami eld survey which have been conducted meticulously along the coast of Japan. Inundation areas as long as runup and tsunami height 1 from in situ measurement, are reported 1 Hereafter for the tsunami terminology we refer to the NOAA [2017b]: Inundation area: an area that is ooded with water. Runup: Maximum height of the water onshore observed above a reference sea level. Usually measured at the horizontal inundation limit.

Elementay dierential geometry

We consider a 2-dimensional Riemannian manifold M. At a point R dened by curvilinear coordinates {z 1 , z 2 }, R(z 1 , z 2 ) let {g 1 , g 2 } be the vector basis.

In the nighbourhood of this point we construct a tangent plane.

Arc lenght between points in curvilinear coordinates

Consider now two points P (z 1 , z 2 ) and Q(z 1 , z 2 ) connected by a curve c parametrized with t ∈ [0, 1]. The arc length of ds between the points Q and Q is given by

And the arc lenght of P Q is computed with a line integral from t = 0 to t = 1:

Surface area in curvilinear coordinates

The area dierential of the tangent plane can be calculated using contravariant components of an innitesimal displacement in the tangent plane

Appendix B

Harmonic maps

We provide the denition of harmonic maps between Riemannian manifold, as emerged from the work of [START_REF] Eell | Harmonic mappings of riemannian manifolds[END_REF]. Given two Riemannian manifold M and N , the former is described by curvilinear coordinates χ with metric tensor G while the latter is described by curvilinear coordinates x, being H its metric. Let A : M → N be a smooth map. Let Ω be a domain of M. We dene an energy functional in local coordinates:

where the energy density is e = ||∇x|| 2 . From (A.11), the gradient of a vector is the following tensor: ∇x = g i h α ∂x α ∂χ i computing its norm, the energy functional reads:

where standard summation convention is assumed. The energy functional reads:

The map x = A(χ) is called harmonic if it is a critical point of the energy functional. Moreover, the harmonic map is the solution of the following Euler-Lagrange equation:

Here Γ k ij and Γ α γδ denote the Christoel symbols on M and N respectively.

Using (A.13) we get:

Lat-Lon coordinates

We consider a sphere S 2 with radius R described by orthogonal curvilinear coordinates {z 1 , z 2 } = {ψ, λ} whose basis are mutually perpendicular but not unitary. z 2 = λ is the latitude and z 1 = ψ the longitude. A vector

We got the transformation {x i } = G({z j })

For lat-lon coordinates, the transformation tensor, according to equation (A.2), is

The covariant basis results

The metric tensor reads

The innitesimal area of a sphere is computed as

Proof of the eRK2-RD-ALE update

We extend the same steps contained in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF] to the ALE framework. As the authors of the same paper underline, the procedure to construct a second order explicit RD scheme is based on the following three steps:

• analogy between RD and Stabilized Finite Element, valid for the P1 case

• step shifting operator for the time part in the stabilization bubble

• high order mass-lumping

The reader will nd these steps in the following. The test function is w i = ϕ i + γ i . Well-Balanced SWEs in ALE framework (2.28) are multiplied by the the test function and integrated over the domain Ω. Using dx m = J A dχ m and the hypothesis that ∂ t | χ ϕ i = 0, one obtain the weak form (Galerkin)

We have split the ALE term in a pure advective term and a Geometric Source Term, as is commonly done in other RD-ALE discretizations [START_REF] Dobes | Second order blended multidimensional upwind residual distribution scheme for steady and unsteady computations[END_REF]; [START_REF] Michler | An arbitrary lagrangian eulerian formulation for residual distribution schemes on moving grids[END_REF]. The next step is to add the contribution