Kévin Pons
email: kevin.pons@principia.fr

Richard Marcer
email: richard.marcer@principia.fr

Mehmet Ersoy
email: mehmet.ersoy@univ-tln.fr

Frederic Golay
email: frederic.golay@univ-tln.fr

Richard Marcer Principia

Frédéric Golay

Frédéric Golay

Adaptive mesh refinement method. Part 2: Application to tsunamis propagation

come

Introduction

Tsunamis are generally referred to any impulsed gravity wave. It can be generated by many sources. The most common tsunamis are a consequence of earthquakes and landslides. Such events can displace a very large volume of water from its equilibrium. The displaced water mass, under the gravity action, attempts to come back to its equilibrium generating gravity waves. Depending of the tsunami source, the order of magnitude of a tsunami wave length is generally larger than ten kilometers. Since this horizontal length is much larger than the oceans depth, vertically integrated model are generally used for tsunamis propagation. One of the most simple and one of the most relevant model among depth-averaged equations is the Saint-Venant model (also called Non Linear Shallow Water -NLSW), see for instance [START_REF] Leveque | High-resolution finite volume methods for the shallow water equations with bathymetry and dry states[END_REF][START_REF] Leveque | Tsunami modelling with adaptively refined finite volume methods[END_REF][START_REF] Popinet | Adaptive modelling of long-distance wave propagation and fine-scale flooding during the tohoku tsunami[END_REF]. Nevertheless, due to the non-dispersive behaviour of this set of equations, the numerical results for some tsunamis problem may lead to under forecasting, for instance, premature breaking-waves. To be more accurate in that case, the use of dispersive and fully nonlinear equations are necessary, such as Green-Naghdi equations ( [START_REF] Shi | A high-order adaptive timestepping tvd solver for boussinesq modeling of breaking waves and coastal inundation[END_REF][START_REF] Popinet | A quadtree-adaptive multigrid solver for the serre-green-naghdi equations[END_REF][START_REF] Duran | A discontinuous galerkin method for a new class of green-naghdi equations on simplicial unstructured meshes[END_REF]) or some non-hydrostatic models ( [START_REF] Zijlema | Swash: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters[END_REF][START_REF] Ma | Shock-capturing non-hydrostatic model for fully dispersive surface wave processes[END_REF]). However, almost all dispersive models are based on the hyperbolic structure of the Saint-Venant equations. Therefore, in this work, we restrict ourselves to the Saint-Venant model and the generalisation to dispersive models can be easily adapted.

Solving accurately this model over very large ranges in spatial scale with high resolution inexorably leads to heavy computational time [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF][START_REF] Popinet | Adaptive modelling of long-distance wave propagation and fine-scale flooding during the tohoku tsunami[END_REF][START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF]. In principle, Adaptive Mesh Refinement (AMR) [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF] or Adaptive Moving Mesh (AMM) [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF] methods allow to solve in a reasonable CPU time these equations adjusting the computational effort locally to maintain a quiet uniform level of accuracy. The first one relies on macro cells which can be refined (and then possibly coarsened) while the other method allows all mesh nodes to move to generate big or small cells (with a fixed nodes number for all the simulation). The zones where the mesh needs to be moved or refined is determined thanks to a mesh refinement criterion. Depending on the application, equations and numerical methods, a variety of different criteria might be used based on an error estimation procedure or a feature detection technique. Basically, the zones where a given threshold is exceeded are moved or refined. In particular, for AMR methods, the more the threshold is small, the more accurate is the result to the expense of the computational time. We focus in the rest of the paper to AMR method.

• The overall efficiency and accuracy of the AMR method can be improved using an Automatic Mesh Refinement Threshold (AMRT) [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF] based on the decreasing rearrangement of the mesh refinement criterion function. The decreasing rearrangement provides a description of the criterion in terms of local maxima which are sorted canonically from the smallest to the largest. Thus, this method allows in practice to set the threshold automatically small enough to catch the most relevant scales to refine with a simple mesh refinement criterion. The method does not require accurate error estimate generally hard to obtain and provided case by case. The method needs only an error indicator as a mesh refinement criterion to provide relevant numerical results having a good balance between the accuracy and the computational cost of the numerical method. Moreover, the AMRT method is almost parameterless while the most of the AMR methods are based on a "user calibrated" threshold, case-dependent, which requires time and effort to set it correctly (see for instance [START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF]). One can easily adapt the parameterless method to AMM methods. • Next, the overall performance of the method can be improved in a parallel framework and more specifically with a Block-Based AMR method [START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF]. Each block corresponds to the initial unstructured mesh and according to the mesh refinement procedure each block can be sub-divided. The cells of each block are redistributed in a fixed number of domains. The number of domain being fixed, each domain are loaded in a given MPI process. The BB-AMR technique provides an efficient control of the shared memory leading to well-balanced computational time between cores by domain-like decomposition. Nevertheless, data locality is critical to obtain good performance since the memory access times are not uniform and may become expensive. Consequently, the re-meshing step has to be carefully managed. Taking advantage of such a block-based structure (as in domain decomposition), we can define two different time-step: the first one is based on the (micro-)CFL δt n condition (i.e computed through the finest cells) which ensures the stability of the numerical scheme. The second one is a (macro-)CFL ∆t n δt n condition based on the finest block and it corresponds to the re-meshing time step. Thus, the global stability of the scheme is preserved for any time and the re-meshing cost is minimised.

In this paper, we present the AMRT method and we focus on its efficiency for the two-dimensional Saint-Venant equations in the context of tsunami propagation. The method is validated through several test cases.

The paper is organized as follows: the first section of the paper is devoted to the presentation of the tsunami model (Sect. 2.1) and its finite volume approximation (Sect. 2.2). The second section is dedicated to the presentation of the multi dimensional finite volume solver on unstructured meshes with AMRT (Sect. 3.1) and the Block-Based Adaptive Mesh Refinement technique (Sect. 3.2). The computational results are compared to experimental data in Sect. 4. The first test case (Sect. 4.1) concerns the flume experiments conducted at Oregon State University [START_REF] Roeber | Shock-capturing boussinesq-type model for nearshore wave processes[END_REF][START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF]. These experiments involve the propagation, run-up and reflection of high amplitude solitary waves on a reef. The second test case (Sect. 4.2) is the propagation of a solitary wave over a complex three dimensional shallow shelf [START_REF] Lynett | Experimental study of solitary wave evolution over a 3d shallow shelf[END_REF]. In this experiment, the propagation, run-up, drying and flooding phenomena are involved. The last test case (Sect. 4.3) concerns the Monai-Walley tsunami run-up onto a complex three dimensional beach [START_REF] Liu | Advanced numerical models for simulating tsunami waves and runup[END_REF].

Unless otherwise indicated, bold characters are used for vector notation and w T denotes the transpose of the vector w.

Equations, properties and finite volume approximation

In this section, we present the Saint-Venant equations which are the standard and one of most the simple governing model used for tsunami propagation. We first recall its mathematical and physical properties which are useful for numerical purposes, such as the entropy inequality. Then, we present the multi dimensional finite volume approximation.

Saint-Venant equations and entropy inequality

The two dimensional non linear Saint-Venant system is

∂ t h + div(hu u u) = 0 ∂ t (hu u u) + div hu u u ⊗ u u u + g 2 h 2 I = -gh∇Z (1)
where the unknowns h(x x x,t) and u u u(x x x,t) = (u 1 (x x x,t), u 2 (x x x,t)) are respectively the height of the water and the depth-averaged velocity of the water at a space-time point (x x x,t), x x x = (x 1 , x 2 ) ∈ R 2 is the space coordinate, t > 0 is the time, g is the gravitational constant set to 9.81m/s 2 , Z(x x x) is the bottom surface elevation and I is the identity matrix. The three equations above express respectively, the conservation laws of mass and momentum driven by the fluxes

f f f 1 (w w w) =   hu 1 hu 2 1 + g 2 h 2 hu 1 u 2   and f f f 2 (w w w) =   hu 2 hu 1 u 2 hu 2 2 + g 2 h 2  
and forced through the source G G G(x x x, w w w) = (0, -gh∂ x 1 Z(x x x), -gh∂ x 2 Z(x x x)) T where the conservative variable w w w = (h,

hu 1 , hu 2 ) T . Noting f f f (w w w) = ( f f f 1 (w w w), f f f 2 (w w w)) T , System (1)
reads in its vectorial form

∂ t w w w + div f f f (w w w) = G G G(x x x, w w w) . (2)
System ( 2) is strictly hyperbolic on the set {h(x x x,t) > 0}. It is well-known that solving equations (2) with high accuracy is a challenging problem since solutions can and will breakdown at a finite time, even if the initial data are smooth, and develop complex structure (shock wave interactions). In such a situation, the uniqueness of the (weak) solution is lost and is recovered by completing the system (2) with an entropy inequality of the form:

S (x x x,t) = ∂ s(x x x, w w w) ∂t + divψ ψ ψ(s(x x x, w w w)) ≤ 0 , (3)
where (s, ψ ψ ψ) stands for a convex entropy-entropy flux pair [START_REF] Serre | Systems of Conservation Laws 1: Hyperbolicity, entropies, shock waves[END_REF][START_REF] Benzoni-Gavage | Multi-dimensional hyperbolic partial differential equations: First-order Systems and Applications[END_REF]. The quantity S is called the entropy production. This inequality in equation ( 3) is used to select the physical relevant solution. Moreover, the entropy satisfies a conservation equation only in regions where the solution is smooth and an inequality when the solution develops shocks. In simple case, it can be proved that the missing term in (3) to make it an equality is a Dirac mass. Therefore, the entropy production can be seen as a "smoothness indicator" at the numerical level, see for instance [START_REF] Puppo | Numerical entropy and adaptivity for finite volume schemes[END_REF][START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF].

For the Saint-Venant system (2), a convex entropy-entropy flux pair is given by

s(x x x, w w w) = h |u| 2 2 + gh 2 2 + ghZ(x x x), ψ ψ ψ(s(x x x, w w w)) = s(x x x, w w w) + gh 2 2 u
where

|u| 2 = u 2 1 + u 2 2 .

Finite volume approximation

The computational domain Ω ⊂ R d , d = 1 or 2, is split into a set of control volumes, also referred as cells,

Ω = ∪ k C k of mesh size |C k |.
The source term is upwinded at the cell interface through the hydrostatic reconstruction [START_REF] Audusse | A fast and stable wellbalanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. As a consequence the obtained numerical scheme is well-balanced on fixed grid [START_REF] Audusse | A fast and stable wellbalanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Simeoni | Remarks on the consistency of upwind source at interface schemes on nonuniform grids[END_REF] that is to say, it preserves the still water steady states

h + Z = constant (4)
and the positiveness of the water height. In the framework of adaptive mesh, we refer to [START_REF] Leveque | Tsunami modelling with adaptively refined finite volume methods[END_REF][START_REF] Donat | Well-balanced adaptive mesh refinement for shallow water flows[END_REF][START_REF] Pons | Adaptive Mesh Refinement Method Applied to Shallow Water Model: A Mass Conservative Projection[END_REF] and the reference therein.

Semi-discrete approximation

In order to get the semi-discrete approximation we proceed in two steps:

Step 1: we first construct an approximation for the homogeneous system. On a given cell C k , noting w k (t)

w w w k (t) 1 |C k | C k w w w (x x x,t) dx 1 dx 2
the approximation of the mean value of the unknown w w w(x x x,t) on C k at time t, and integrating (2) over each cell, we obtain:

C k ∂ w w w(t) ∂t dx 1 dx 2 + ∑ a ∂C k/a f f f (t, w w w) • n n n k/a ds = 0 (5)
where n n n k/a denotes the unit normal vector on the boundary ∂C k/a between cells k and a.

Next, F w w w k (t), w w w a (t), n n n k/a the flux approximation being written

F F F w w w k (t), w w w a (t), n n n k/a ≈ 1 |∂C k/a | ∂C k/a f f f (t, w w w) • n n n k/a ds ,
the semi-discrete finite volume approximation of equations (2) (see for instance [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF][START_REF] Eymard | Finite volume methods[END_REF]) is obtained:

∂ w w w k (t) ∂t + 1 |C k | ∑ a |∂C k/a |F F F w w w k (t), w w w a (t), n n n k/a = 0 ( 6
)
where F F F w w w k (t), w w w a (t), n n n k/a is defined via the Godunov solver, i.e. computed with the exact solution of the 1D Riemann problem at the interface k/a with the states w w w k (t) and w w w a (t). A second order approximation in space, the MUSCL reconstruction is also implemented (see for instance [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF][START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF]). To be consistent with the above discretisation of the equations (2), note that the entropy inequality (3) is approximated as follows

S k (t) = ∂ s k (t) ∂t + 1 |C k | ∑ a |∂C k/a |φ φ φ w w w k (t), w w w a (t), n n n k/a (7)
where φ φ φ (w w w k (t), w w w a (t), n n n k/a ) represents the entropy flux calculated from the resolution of the Riemann problem at the interface of cells k and a as done above.

Step 2: To take into account the source term G G G, we use the hydrostatic reconstruction. For each interface, the numerical flux F F F w w w k (t), w w w a (t), n n n k/a is replaced by the flux at left hand side F F F k w w w k (t), w w w a (t), n n n k/a , ∆ Z k/a and the right hand side F F F a w w w k (t), w w w a (t), n n n k/a , ∆ Z k/a of the interface k/a. These new fluxes are given by

F F F k w w w k (t), w w w a (t), n n n k/a , ∆ Z k/a = F F F w w w * k , w w w * a , n n n k/a +   0 g 2 (h 2 k -(h * k ) 2 )n 1 g 2 (h 2 k -(h * k ) 2 )n 2   F F F a w w w k (t), w w w a (t), n n n k/a , ∆ Z k/a = F F F w w w * k , w w w * a , n n n k/a +   0 g 2 (h 2 a -(h * a ) 2 )n 1 g 2 (h 2 a -(h * a ) 2 )n 2  
by means of reconstructed states (satisfying the still water steady state equation ( 4))

w * k = (h * k , u k ) , w * a = (h * a , u a ) , h * k = max(0, h k -max(0, ∆ Z k/a )) , h * a = max(0, h a -max(0, -∆ Z k/a )) .
In theses formula, ∆ Z k/a stands for the jump of Z across the interface k/a. The scheme is by construction well-balanced, i.e., the still water steady states are ex-actly satisfied (see [START_REF] Pons | Adaptive Mesh Refinement Method Applied to Shallow Water Model: A Mass Conservative Projection[END_REF] for further details). Finally, the numerical entropy fluxes are replaced with the reconstructed one as done above for the numerical fluxes F F F.

Time integration

The time integration of eqs. ( 6) and (3) can be achieved in a classical way either by a Runge-Kutta or Adams-Bashforth scheme. For instance, by integrating eq. ( 6) (and (3)) during the time step ]t n ,t n+1 [ of length δt n and by evaluating the numerical fluxes at time t n , the well-known second order Runge-Kutta method is

w w w k (t n+1 ) = w w w k (t n ) - δt n |C k | ∑ a |∂C k/a |F F F w w w k (t n+1/2 ), w w w a (t n+1/2 ), n n n k/a (8)
where

w w w k (t n+1/2 ) = w w w k (t n ) - δt n 2|C k | ∑ a |∂C k/a |F F F w w w k (t n ), w w w a (t n ), n n n k/a .
Remark 1. Note that, even if the Adams-Bashforth scheme is known to be less stable and less accurate, it can be easily handled in the framework of local time stepping to save computational time as done by Ersoy et al. [START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF]. Therefore, the code can support sub-cycling in time.

To be consistent with the discretisation of equation ( 8), a discrete version of the entropy production, called numerical density of entropy production S n k , is defined by integrating the equation [START_REF] Cockburn | A priori error estimates for numerical methods for scalar conservation laws. Part II: Flux-splitting monotone schemes on irregular Cartesian grids[END_REF] with a second order Runge-Kutta scheme.

BB-AMR method

Defining a robust mesh refinement criterion for multi dimensional configurations is not enough to design a suitable numerical solver. The mesh refinement threshold and the treatment of data are also crucial points and in particular the way to distribute refined cells and to share the memory in a parallel process. This point is handled within an Automatic Mesh Refinement Threshold framework in a hierarchical Block-Based way called BB-AMR. The global strategy to adapt the mesh and to manage data are presented.

AMR and Automatic thresholding methods

Numerical approximation of equations [START_REF] Cockburn | A priori error estimates for numerical methods for scalar conservation laws. Part II: Flux-splitting monotone schemes on irregular Cartesian grids[END_REF] leads to the so-called numerical density of entropy production which is a measure of the amount of violation of the entropy equation (as a measure of the local residual as in [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Houston | A posteriori error analysis for numerical approximations of Friedrichs systems[END_REF][START_REF] Karni | A smoothness indicator for adaptive algorithms for hyperbolic systems[END_REF][START_REF] Karni | Local error analysis for approximate solutions of hyperbolic conservation laws[END_REF]). Therefore, the numerical density of entropy production provides information on the need to locally refine the mesh (e.g. if the solution develops discontinuities) or to coarsen the mesh (e.g. if the solution is smooth and well-approximated), see for instance [START_REF] Puppo | Numerical entropy and adaptivity for finite volume schemes[END_REF][START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF]. As shown in [START_REF] Puppo | Numerical entropy and adaptivity for finite volume schemes[END_REF], if the numerical entropy flux is well-chosen, the entropy production S built on first order monotone schemes is almost everywhere negative definite. Some positive overshoots of the entropy production S are possible for non monotone data but their amplitude decreases fast under grid refinement for smooth flow. Moreover, as numerically shown for one dimensional hyperbolic equations [START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF][START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF] and for multi dimensional flows [START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF], the support of the relative error almost coincides with the support of the numerical density of entropy production for first and second order schemes but the absolute magnitude max k S n k is not representative of the error. Therefore, it is important to emphasize that the numerical density of entropy production is not an error estimate, or derived from, but only an error indicator. Thus, the numerical density of entropy production does not allow to relate to the equidistribution of the error as for a priori or a posteriori based error criteria which are generally hard to obtain.

In what follows, we define the quantity S n k = S n k also referred as the local numerical entropy production. The local numerical entropy production is compared to

S m = 1 |Ω | ∑ k S n k .
In practice, two coefficients 0 β min β max 1 are thus defined to determine the ratio of numerical production of entropy leading to mesh refinement or mesh coarsening.

For each cell C k :

• if S n k > S m β max , the mesh is refined and, • if S n
k < S m β min the mesh is coarsened whenever it is possible following the rule defined later. This approach is certainly the most simple and the most relevant for smooth flows. For discontinuous flows, the comparison of the error indicator to the mean value S m may lead to some problems to catch relevant small scales to be refined. In Fig. 1, we have illustrated formally a discontinuous flow, see Fig. 1(a), and a smooth one, see Fig. 1(b). In Fig. 1(a), the peak localised at x = x B = 3.75 corresponds either to a sharp gradient or a discontinuity in the solution while the other part is the region where the solution is smooth. One can see in Fig. 1(a) that if the parameters β min is not set small enough then the smooth part is not refined and only a small region around the sharp gradient or the discontinuity of the solution is refined. This is exactly what we mean by "some problems to catch relevant small scales to be refined". As illustrated in Fig. 1(b) and in general for smooth flows, the mean value S m allows to catch relevant scales to be refined.

Therefore, the parameters β min and β max have to be well-chosen (see for instance [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF]Section 4] for a review of some thresholding methods). The threshold parameters β min and β max are thus user-calibrated for each problem to reach a good compromise between computational cost and accuracy. More precisely, β min and β max allow to set a percentage of mesh refinement and mesh coarsening with respect to the quantity S m . It is not surprising that these settings will deteriorate or improve the accuracy of the numerical solution. For instance, the more β min and β max is small, the more accurate are the results to the expense of the computational time. Calibrating such parameters is a computational waste of time especially for real life numerical simulations which can take several hours to several weeks of computations. Unless to have an error estimate, the adaptive method with an error indicator has a major drawback concerning the selection of the threshold parameters β min and β max . For the sake of simplicity, assuming β min = β max = 1, one way to overcome this problem is to set the threshold α automatically. Since α = S m is in general too high for discontinuous flows and large enough for smooth flows (as illustrated in Fig. 1), it is then natural to look for α in the interval (0,

S m ].
In what follows, we focus on the one dimensional problem because it is easier to visualise. Keeping in mind that we want to catch relevant small scales to be refined, we look for an automatically set threshold α ∈ (0, S m ] such that for each cell C k , if S n k > α the mesh is refined otherwise the mesh is coarsened. Remark 2. As it is, such a strategy is time consuming because the mesh would flipping back and forth between fine and coarse grids. Therefore in practice, once the threshold parameter α is set "correctly", one can define the two quantities β min and β max to define a buffer region where the mesh is left alone. For instance, the top 10% is refined and the bottom 10% is coarsened while the rest of the mesh remains at the same level of refinement. Let us also add that even if β min is set to 1, the method will detect automatically pertinent scales.

Following [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF], a relevant way to define automatically the threshold α, referred as α PE in the sequel, is to construct it such that

f (α PE ) = max 0<α S m f (α) with f (α) = αd(α) (9)
In this expression, d : R → R + , α → meas{S(x) > α} is the decreasing rearrangement, also known as the distribution function, of the mesh refinement criterion S. meas{X} stands for the Lebesgue measure of the set X. For instance in Fig. 1(b), for α = S m , d(α) = meas{S(x) > α} = l1 + l2. The decreasing rearrangement d provides a complete description of the criterion in terms of local maxima which are sorted canonically from the smallest to the largest. Following [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF], assuming that the function S admits a finite number of strict local maxima, then d is only continuous and each singularity of its derivative corresponds to a local maximum as displayed in Fig. 2 for discontinuous flows and smooth flows. For instance in Fig. However in these examples with the distribution function only, it is hard to compute numerically a threshold α close to the first local maximum α A which represents the so-called "relevant scales to be refined". Therefore, we use the "weighted distribution function" f (α) = αd(α) which is represented in Fig. 3 and4(b) for discontinuous flows and smooth flows. As proved in [29, Corollary 1], the function f transforms almost the inflexion points of d in local maximum. As a consequence, the relevant threshold α PE is in general smaller than α A for discontinuous flows or equal to S m for smooth flows. Moreover, the threshold α PE is easier to compute numerically from f than d.

We have illustrated the computed threshold for the previous example for discontinuous flows in Fig. 5(a) (see also a zoom in displayed in Fig. 4(b)) and for smooth flows in Fig. 5(b). With this method in practice, we can compute the relevant threshold to refine for both flows. Namely, the threshold α PE verifies the desired behaviour α PE < S m for discontinuous flows and α PE ≈ S m for smooth flows. Thus, the presented AMR process, in comparison with existing methods, is almost parameterless.

Following [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF], from a numerical viewpoint, the distribution function is defined as follows. For each time step, let us consider a given discrete mesh refinement criterion

S(x) = ∑ N k=1 S n k 1 C k (x)
where N is the number of cells at time t n and 1 C k (x) is the indicator function of the mesh C k . Let (α j ) 0 j M be an increasing sequence of M + 1 threshold parameter such that

(α j ) 0 j M = S m j M β 0 j M (10)
Then, the distribution function d j = d(α j ) is given by where # is the number of elements in the set {S n k > α j }. In [29, Section 4], it is shown that the method is quite insensitive to the parameter β , set to 2 here. For all applications M is fixed to M = 1000 points. Interested readers can found details in [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF] and a review of existing threshold methods.

d j = #{k ; S n k > α j } (11)

Mesh refinement process and data management

For the one dimensional case, the local mesh refinement procedure is constructed following dyadic tree applied at each time step. "Macro cells" are used to be easily refined by generating hierarchical grids. Each cell can be split in two. For stability reasons, the mesh refinement level cannot exceed 2 between two adjacent cells as detailed in [START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF] and recalled hereafter. By convention, if the maximum the level of refinement is l max = 1 then the mesh is not adapted.

The multi dimensional extension of the mesh refinement procedure is a difficult task. Interesting works have been presented for 2D Cartesian grid or quad-tree [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Yiu | Quadtree grid generation: Information handling, boundary fitting and cfd applications[END_REF][START_REF] Min | A second order accurate level set method on non-graded adaptive cartesian grids[END_REF][START_REF] Zhang | A two dimensional hydrodynamic and sediment transport model for dam break based on finite volume method with quadtree grid[END_REF], octree for 3D simulations [START_REF] Losasso | Simulating water and smoke with an octree data structure[END_REF][START_REF] Fuster | Numerical simulation of droplets, bubbles and waves: state of the art[END_REF], and anisotropic AMR [START_REF] Coupez | Solution of high-reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing[END_REF][START_REF] Hachem | Immersed stress method for fluid structure interaction using anisotropic mesh adaptation[END_REF]. For the Saint-Venant system (1) the extension from 1D to 2D leads naturally to quad-tree meshing. But, the presence of a complex moving interface (composed of rarefaction and shocks) implies to re-mesh at each time step, which is obviously a costly process. Guided by the need to reach a relevant compromise between the contradictory aims of solution accuracy and computing speed, a regular block-based mesh approach is introduced, somehow like in [START_REF] Williamschen | Parallel anisotropic block-based adaptive mesh refinement algorithm for three-dimensional flows[END_REF][START_REF] Zheng | Block-based adaptive mesh refinement finite-volume scheme for hybrid multi-block meshes[END_REF]. If the mesh is not refine at each time step, the patch where the grid should be refine must be enlarged. The grid generation in the framework of BB-AMR can be developed as follow.

• Firstly, the initial computational domain is divided in several quadrilateral elements, called blocks as represented in Fig. 6(a).

• Secondly, we define for each block a discretization n x in x-direction, n y in ydirection and a level of refinement l r such that the total number of cell in xdirection of the block will be 2 l r n x . As the mesh refinement level cannot exceed 2 between two adjacent blocks (see [START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF]), the level of mesh refinement is automatically adapted if necessary as illustrated in Fig. 6(b) for the blocks B5 and B2. • Thirdly, in order to balance the CPUs load, the cells of each block are almost equi-distributed in a fixed number of domains according to the Cuthill-McKee numbering, see Fig. 6(c). • The unstructured mesh1 is built for each domain. At the interface between two domains, ghost cells are created in order to transfer the necessary data to the explicit time integration scheme at each time step. The number of domain being fixed, each domain are loaded in a given MPI process. These MPI processes are then load on a fixed number of cores (not necessary the same). For each refined cells (or blocks), averaged values are projected onto each sub-cell and fluxes are computed as simply as possible to avoid heavy computation following Osher and Sanders [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF]. Even if the Osher and Sanders projection leads to non conservative and non consistent scheme at the space-time grid between two levels of refinement, it has the advantage to be less time consuming than complex approach such as [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Tang | A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space grids[END_REF][START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF]. Let us also note that the apparent loss of consistency on the local truncation error seems not to affect the actual error of the scheme (see for instance [START_REF] Cockburn | A priori error estimates for numerical methods for scalar conservation laws. Part II: Flux-splitting monotone schemes on irregular Cartesian grids[END_REF][START_REF] Simeoni | Remarks on the consistency of upwind source at interface schemes on nonuniform grids[END_REF]). Therefore, a simple issue is proposed in [START_REF] Ersoy | Adaptive multiscale scheme based on numerical density of entropy production for conservation laws[END_REF]. As said above, the loss of the consistency and conservativity error can be minimised using smoothing grid technique, that is to say, we prevent two adjacent cells from having a level difference greater than two. More Details on the BB-AMR can be found in [START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF][START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF]. Finally, to compute the threshold automatically in any dimensional framework, the decreasing rearrangement d of the mesh refinement criterion function S is defined (cf. eq. 11). Then, the threshold is automatically defined as explained in the previous section (cf. eq. 9). In the context of Block Based AMR, if one cell or more of any block has a criterion below the threshold, then the entire block is refined.

Remark 3 (A comment on the Osher and Sander projection for vertically averaged Saint-Venant system). To be rigorous with respect to the underlying physical properties of the system, the projection have to be mass, momentum and energetically conservative. The smoothing grid technique works well for the most of the system but may fails for vertically averaged models like Saint-Venant equations. The main reason is related to the topography variation which set some difficulties for this simple projection method. As it is for any three dimensional models, the bathymetry does not act on the flow as a boundary condition but acts directly into the model as a source term. As a consequence, we have to pay attention to the new definition of the bathymetry onto the new refined or coarsened mesh as done in, for instance, [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Leveque | Tsunami modelling with adaptively refined finite volume methods[END_REF][START_REF] Pons | Adaptive Mesh Refinement Method Applied to Shallow Water Model: A Mass Conservative Projection[END_REF] which propose to interpolate the free surface η = h + Z instead of the conservative variable h. This interpolation strategy prevents the introduction of new extrema in the water surface elevation and preserves a flat sea surface. In that paper, the simple approach developed in [START_REF] Pons | Adaptive Mesh Refinement Method Applied to Shallow Water Model: A Mass Conservative Projection[END_REF] is used. This method is mass conservative and allows the computation of the solution close to or at dry states near the shore, being well-balanced by construction if the initial scheme is. For the ease of the reader and for the simplicity of the reader we describe the method in the one dimensional

Keeping in mind that h = η -Z and we impose η = cst for steady states, this relation yields to

Z 0 = η 1 - |C 00 | |C 0 | - |C 01 | |C 0 | + |C 00 | |C 0 | Z 00 + |C 01 | |C 0 | Z 01 Straightforwardly, if |C 00 | = |C 00 | = |C 0 | 2 , the mass conservation is satisfied if Z 0 = Z 00 +Z 01 2 .
step 2 (mesh refinement): Once the bathymetry is defined, the water height is reconstructed from h i = max(0, η i -Z i ) and the free surface conservation η = cste yielding to h 0i = max(0, h 0 + z 0z 0i ) where the max plays a crucial role near shore for dry states. The water speed u 0i is reconstructed form the momentum conservation relation similar to the previous one. step 2 (mesh coarsening): In that case, to go from the fine grid values to a coarsegrid value we average by setting h 0 = |C 00 |h 00 +|C 01 |h 01

|C 0 |
and define the water speed u 0 = |C 00 |h 00 u 00 +|C 01 |h 01 u 01 |C 0 |h 0 .

Numerical validations

We now focus on the overall performance of the BB-AMR scheme with automatic thresholding confronted to experimental and state-of-the-art numerical results 2 . For each test case, the numerical simulations are carried out on adaptive grids (using the BB-AMR method previously presented) and on uniform grids. The minimum adaptive cell size δ x adaptive is equal to the uniform one δ x uni f orm . The comparisons of both simulations allow to show the reliability and the efficiency of the BB-AMR method. To simplify the analysis, we have considered only 1 domain, i.e., 1 processor. Interested readers can found a detailed numerical study with several domains in [START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF]. For all computations,

• the MUSCL reconstruction and the Runge Kutta method are used to achieve a second order scheme in space and time. • to simplify the analysis and to focus on the automatic threshold, we have chosen β min = β max = 1. Indeed, one can improve the above results by setting 0 < β min < β max < 1.

• the CFL number is set to 0.9, • For the sake of simplicity, the re-meshing time step is fixed constant to 0.25 s.

This choice roughly verifies the macro CFL condition. • the code is written in Fortran 90, compiled with gfortran on Linux machine Intel(R) Core(TM) i5-2500 CPU @ 3.30GHz.

Solitary wave propagation over a two dimensional reef

This benchmark aims to reproduce a set of laboratory experiments carried out at the O.H. Hinsdale Wave Research Laboratory, Oregon State University (OSU, see Roeber et al. [START_REF] Roeber | Shock-capturing boussinesq-type model for nearshore wave processes[END_REF] and Roeber and Cheung [START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF]). These experiments involve the propagation, run-up, splash-up and reflection of high amplitude solitary waves on two-dimensional reefs. Their purpose is on one hand to investigate processes related to breaking, bore formation, dispersion, and passage from sub-to super-critical flows, and on the other hand, to provide data for the validation of near-shore wave models in fringing reef. Even if the Saint-Venant model ( 2) is not relevant for dispersive flows, we show, through this test case, that globally experimental data are well-reproduced and the cells are well-refined thanks to the automatic thresholding technique.

Experimental settings

The geometry of the test considered here is shown in Fig. 7. The length of the basin is of 104 m, however the computational domain is delimited by a reflecting wall placed at x = 83.7 m. The reef starts at x = 25.9 m with a nominal slope of 1/12. The actual slope is such that the height of 2.36 m is reached after x = 28.25 m. At this station a 0.2 m height crest is mounted. The offshore slope of the crest is the same of the reef, and the length of its plateau is of 1.25 m. The on-shore side has a slope of 1/15 giving a nominal length for the crest basis of 6.65 m (using the actual offshore slope, a crest basis of 6.644 m is obtained). This gives an offshore length of the crest slope (starting at x = 28.25 m) of 2.4 m. Thus, the bathymetry is The initial depth at still water is h 0 = 2.5 m, yielding to a partially submerged crest, and a depth behind it (on-shore side) of 0.14 m. The initial solution consists of a solitary wave of amplitude A = 0.75 m providing a non linearity ratio of A h 0 = 0.3. For numerical purpose, the initial data consist of a solitary wave centered at x 0 = 17.6 m of amplitude [START_REF] Eymard | Finite volume methods[END_REF] and velocity u(x, 0) = g(h 0 + A)η(0, x)/h 0 [START_REF] Fuster | Numerical simulation of droplets, bubbles and waves: state of the art[END_REF] which corresponds to the experimental data at the dimensionless time t = t g h 0 = 55.03 as displayed in Fig. 8(a).

Z(x x x) =            -2.5 if x < 25.
η(x, 0) := h(x, 0) + Z(x x x) = A cosh 3A (4h 3 0 ) (x -x 0 ) 2

Numerical results vs experimental results

Figures 8(a) and 8(b) shows the propagation of the solitary wave over the slope starting at x = 25.9 m. Figure 8(c) represents the surface elevation at the dimen-Fig. 7 Experimental settings and wave gauges locations sionless time t ≈ 69 where the wave develops into a plunging breaker onto the reef crest 3 . From an experimental viewpoint, t ≈ 69 corresponds to the subsequent time of overturning of the free surface. In Fig. 8(d), the wave hitting the reef creates a downstream and upstream propagating bore as shown in Fig. 8(e). The downstream propagating bore is reflected at the end wall at x = 83.7 m and propagates in opposite direction as displayed in Fig. 8(f).

In Figs. 8(a)-8(f), the experimental surface elevation (represented with blue circles) for several dimensionless time t are compared to the obtained numerical results (green and red lines-circles). The numerical solution is computed on a uniform grid (green lines-circles) composed of 1000 cells and on an adaptive grid (red lines-circles) initially composed of 200 blocks, composed initially with 1 cell, with l max = 3 where l max is the maximum level of the mesh refinement. Free boundary condition at x = 0 m and reflecting boundary condition at x = 83.7 m are imposed.

The obtained numerical results, as shown in Figs. 8(a)-8(f), are rather wellcomputed when compared to experimental data. However, the wave height is slightly underestimated (see Fig. 8(c)) while the hydraulic bore height is overestimated (see Fig. 8(d)).

We now compare in Figs. 9(a)-9(f), the computed and the experimental water elevation at the wave gauges displayed in Fig. 7. The oscillating shape of the experimental measures almost corresponds to dispersive effects which cannot be detected (f) t = 125.03

Fig. 8 Surface profiles of solitary wave propagation over an exposed reef crest. Confrontation of experimental data (blue circles) to numerical data computed on a uniform grid (solid green line) and on an adaptive grid (solid red lines). The solid cyan line represents the mesh level and the black one the bathymetry with the Saint-Venant model [START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, efficiency[END_REF]. However, the shape of the numerical results almost coincides with the experimental data for wave gauges located before x = 54.4 m as shown in Figs. 9(a)-9(c). The point x = 54.4 m is located just before the crest and in particular it corresponds to the region where the overturning of the free surface is experimentally observed. Numerically a small shift is observed in Figs. 9(d)-9(f).

AMR performance

We now compare the performance of the adaptive and uniform method. The adaptive scheme uses an average of 281 (varying up to 506) cells against 1000 cells (see Fig. 10(b)) for the uniform one and the CPU-time is almost around 7 s against 30 s respectively. This computational time can be improved if the top 10% is refined and the bottom 10% is coarsened for instance.

We observe that the results on adaptive grid are similar to the uniform grid ones and computed about 4 times faster. We see also in Figs. 8(a)-8(f) that for each time step, the threshold parameter is automatically well-set to capture efficiently the region to refine (see also Fig. 10(a)). During the first 50 non-dimensional time, we have already pointed out that the solitary wave propagates and splash-up. This phenomenon is characterized by a steep gradient flow followed by a discontinuous flow for the Saint-Venant system, see Fig 8. In view of the remarks in Sect. 3.1, we observe that for almost all t ∈ (55, 100), the threshold is set to α PE < S m . For t > 100, the flow is almost smooth and therefore α PE ≈ S m . It is also interesting to highlight that only the region of interest are automatically refined as shown in Fig. 8(e). Therefore, we recover the desired behavior predicted by Pons and Ersoy [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF]Section 4]. Moreover, in view of the confrontation with the computation on uniform grid, this test shows the reliability/efficiency of the automatic selection of the threshold for multiple scale flows.

Solitary wave propagation over an irregular three-dimensional shallow shelf

In this test case, the numerical simulation of a solitary wave propagation over a complex three dimensional bathymetry is performed. This test case was experimentally introduced in [START_REF] Lynett | Experimental study of solitary wave evolution over a 3d shallow shelf[END_REF] to understand the turbulence and kinematic properties associated with a breaking solitary wave. In this experiment, the propagation, run-up, drying and flooding phenomena are involved.

Experimental settings

A laboratory experiment was conducted in a large wave basin which was 48. and fluid velocities were recorded at several locations by wave gauges (WG) and velocity captors (ADV). A single solitary wave of height of 39 cm at x = 5 m is produced with a piston-type wave-maker. The bathymetry begins with a flat part allowing to generate a solitary wave. Then, starting at x = 10.2 m of the wavemaker and extending to x = 25.5 m, a complex three dimensional shelf is built as illustrated in Fig. 11. The associated slope to this zone is variable but less and less steep. Beyond that, a very small one dimensional slope finishing on a flat onshore zone is built. The experimental basin is closed by walls. Water height and velocities are recorded during 45 s. Finally, the material used for the walls and the bathymetry is made of smooth concrete so that one can assume almost a frictionless flow. For further informations on this experiment, interested readers can found more details in [START_REF] Lynett | Experimental study of solitary wave evolution over a 3d shallow shelf[END_REF].

Numerical vs experimental results

The domain is numerically extended to x = -5 m with a water depth of 78 cm to impose a solitary wave instead of reproducing the wave maker movement. The simulated solitary wave is a first order solution of the Boussinesq equation (see equations ( 13) and ( 14)) with x 0 = 5 m and A h 0 = 0.5. For computational purpose, we have considered 128 initial blocks composed of 7 500 (varying up to 25 000) cells for the initialization of the adaptive computation and almost 33 000 cells for uniform mesh computation. The simulation time is 30s. Reflecting boundary conditions are prescribed to walls.

To illustrate the main propagation phenomena, we have displayed the numerical solutions calculated at time t = 0.5 s (see Fig. 12), t = 2.5 s (see Fig. 13), t = 5.75 s (see Fig. 14) and t = 23.75 s (see Fig. 15). The Fig. 12 represents the profile of the solitary wave which propagates towards the coast. Figures 13 and14 show respectively the wave hitting the conic island and then the beginning of the flooding. The global numerical results are now confronted to the experimental data recorded at the waves gauges 2, 4, 6 and 7 (WG) for the water height, see Fig. 11 and16, and to the velocity captors 2 and 3 (ADV) for the first u component of the velocity vector, see Fig. 11 and17. For the wave gauge 2 to 7 (see Figs. [START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF](a)-16(d)), the global shape of the free surface is in a good agreement with the experimental data but some differences are observed. The numerical solution computed with the Saint-Venant system (2) slightly underestimates the water level at the wave gauge 2 (see Fig. 16(a)) and overestimates at the wave gauge 4 (see Fig. 16(b)). Moreover, the numerical solution is rather accurate up to time t = 15 s at the wave gauges 6 and 7 (see Fig. 16(c) and 16(d)). For t 15 s , the dispersive effects, not reproduced with the Saint-Venant system, begin significant and induce slightly different results between experimental and numerical data. Similar conclusion can be drawn up for the first component of the vector velocity on the captor 2 (see Fig. 17(a)) and 3 (see Fig. 17(b)). Globally, the numerical results on this test case are in a good agreement and can be improved with a dispersive model.

AMR performance

To highlight the efficiency of the mesh refinement parameter and the automatic selection of the threshold, the numerical density of entropy production (without absolute value) and the associated mesh are displayed at several times in Figs. 18-21, (time t = 0.5 s in Fig. 18, for t = 2.5 s in Fig. 19, for t = 5.75 s in Fig. 20 and for t = 23.75 s in Fig. 21). The time evolution of the mesh refinement threshold selecting the region to be refined is shown in Fig. 22(a). It can be noticed that according to the complex flow zones the numerical density of entropy production indicates which areas need to be refined. Secondly, the mesh density is well distributed on these zones which illustrate the efficiency of the automatic thresholding method.

For instance, in Fig. 14, the region which needs to be refined is clearly the zone around the "conic island". According to Fig. 20, the numerical density of entropy production in this area shows that the local maxima (lowest value without the absolute value) are of almost of order 0.005 and clearly follows the wave front. We see in Fig. 22(a), approximately at time t = 5.75 s, that the threshold is fixed to almost 0.001 and therefore allows to perform a suitable refinement in the region of interests as observed in Fig. 14 or 20.

Compared to the computation on a uniform grid (see Figs. 16 and 17), the adaptive mesh refinement scheme requires only an average of 13000 cells against 33000 cells (see Fig. 22(b) for the number of cells for all the simulation). For the same accuracy, the adaptive scheme allows to save almost 60% of cells with respect to the uniform mesh. During the adaptive simulation, the maximum refined area is reached around the time t = 12 s which corresponds to the time where the free surface oscillates behind the island and the wave run-up on the "coast". This result makes sense because the mesh has to follow a lot of oscillation and not only a solitary wave. Globally, all these cells saved allow to speed up the computation by 2.5 time for this test case. This test case concerns the Hokkaido-Nansei-Oki tsunami, in 1993, that struck the Okushiri Island in Japan. The tsunami run-up records was about 30 m height and the currents speed of order 10-18 m/s for which relevant high-quality data were saved. 4 The largest value run-up (32 m) was recorded near the Monai-Walley beach.

Experimental settings

To understand this complex run-up, a 1/400 scale laboratory model of Monai was realized in a large-scale tank of 205 m long, 6 m deep, and 3.5 m wide at the Central Research Institute for Electric Power Industry (CRIEPI) in Abiko (Japan). The offshore incident wave on a water depth d = 13.5 cm is prescribed. There are reflective side-walls at y = 0 and 3.5 m as for the x boundaries. The bathymetry as well as the coastal topography reproduced in the laboratory experiment are represented 5 in Fig.

23(b) and 23(c) .

The input wave at x = 0 m is a leading-depression height of -2.5 mm with a crest of 1.6 cm, as displayed in Fig. 23(a). In the experiment, the waves are measured at thirteen locations and complete time histories are given at three locations, (x, y) = (4.521, 1.196), (4.521, 1.696), and (4.521, 2.196) in meters, see also [START_REF] Liu | Advanced numerical models for simulating tsunami waves and runup[END_REF]. In contrast with the previous experimental test cases, the dispersive phenomena can be fully neglected here. The initial condition for this test case is a lake at rest and an imposed water height on the left side wall, see Fig. 23(a). Wave gauges data 6 experimental results are given up to time t = 100 s. However, the information concerning the input wave profile (see Fig. 23(a)) at x = 0 m is only available up to time t = 22 s. For numerical purpose, we prescribe for t > 22 s, a constant free surface level η = 0 m on the left wall. In practice, due to this assumption, some artificial reflecting waves can occur on the wave-maker beyond t > 22 s and can reach the free surface gauges around t = 30 s. Finally, to avoid any wrong comparisons due to these artificial reflections, the comparisons between the experimental gauges and the numerical results will be only done over the first 30 s of the experiment. For computational purpose, we have considered 240 initial blocks composed of 8 000 (varying up to 40 000) cells for the initialisation of the adaptive simulation and 62 000 cells for uniform mesh simulation. The time simulation is 30s. Reflecting boundary conditions are prescribed to walls.

To illustrate this test case we show in Figs. 24-27, the propagation of the wave at time t = 11.25 s (Fig. 24), t = 13.25 s (Fig. 25), t = 16 s (Fig. 26) and t = 17.5 s (Fig. 27). The Fig. 24 represents the tsunami wave entering inside the computational domain. The Fig. 25 shows the submersion of the conical island. Figures 26 and27 show respectively the large flooding and run-up of the wave and its reflection on the cliff. The numerical solutions are now quantitatively confronted to the experimental measures. The comparisons are shown in Fig. 28 at the three waves gauges WG 1 (see Fig. 28(a)), 2 (see Fig. 28(b)) and 3 (see Fig. 28(c)). These figures show that the wave amplitude and the wave shape are accurately predicted at three different location in the computed domain. As emphasized before, dispersive effects being negligible, the numerical results are accurate. However, it can be noticed in Fig. 28 a small amplitude errors occurring at the beginning of the simulations. This error can be attributed to the lack of accuracy of the prescribed initial condition and boundary condition at x = 0. This problem was already encountered by Popinet [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF].

AMR performance

To highlight the efficiency of the mesh refinement parameter and the automatic selection of the threshold, the numerical density of entropy production (without absolute value) and the associated mesh are displayed at several times in Figs. 29303132, (time t = 11.25 s in Fig. 29, for t = 13.25 s in Fig. 30, for t = 16 s in Fig. 31 and for t = 17.5 s in Fig. 32). We also display in Fig. 33(a) the evolution of the threshold parameter which selects the region to be refined. We observe mainly that the wave front and subsequent waves are very wellcaptured by the criterion and the mesh refinement threshold is well chosen to refine in those areas (see for instance Fig. 27 and32). Between 3 and 7 seconds, the tsunami propagates in the domain involving a grow-up of the number of cells. Then the wave passes through the island and reaches the reef around 15 seconds. After this, the wave is reflected and goes back to the wave maker with a second interaction with the island. The maximum number of cells is reached when the wave is passing on the island and begins to run-up the coast (see Figs. [START_REF] Pons | Adaptive Mesh Refinement Method Applied to Shallow Water Model: A Mass Conservative Projection[END_REF] and 33(a)), and when the wave goes back to the wave maker. These results are coherent with the physical process involved since it corresponds to the both times where multiple active flow regions are observed. The overall scales of the tsunami propagation with its reflections are very well captured as displayed in Fig. 29, Fig. 31 and Fig. 32.

Finally, the time evolution of the number of cells is represented in Fig. 33(b). Compared to the computation on the uniform grid (see Fig. 28), the adaptive mesh refinement method requires an average of almost 25000 cells against 62000 cells. For the same accuracy, the AMR method allows to save almost 60% of cells with respect to the uniform simulation. Thus, the AMR method allows to speed up the computation almost 3 times.

Conclusion

In this paper, we have considered the one and two dimensional non linear non dispersive Saint-Venant system for tsunamis propagation problems. Solving accurately this model over very large ranges in spatial scale leads to heavy computational time. Therefore, we have considered a well-balanced finite volume scheme in adaptive mesh framework. The BB-AMR method is improved using an Automatic Mesh Refinement Threshold based on the decreasing rearrangement of the mesh refinement criterion function [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF] yielding to a robust AMR method.

To test the efficiency of the overall method, and in particular of the automatic threshold technique, we have considered three classical large scale tsunami propagation test cases involving run-up, reflection and/or flooding and drying phenomena. As expected for the non dispersive test cases, we have obtained results in a very good agreement. In the case of weakly dispersive flows, we have obtained some small errors as pointed out before by several authors. However, the presented results can be improved using depth-averaged dispersive model such as the Boussinesq (weakly non linear) or the Green-Naghdi (fully non linear) one for instance. As emphasized in [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF], the use of the automatic threshold allows to localise almost all regions to be refined balancing the computational cost and the accuracy. The Automatic Mesh Refinement Threshold technique makes the AMR method almost parameterless and robust while the existing AMR methods are based on a fixed threshold and therefore are test case-dependent.

The overall efficiency and performance of the method has been shown through these three different test cases. Finally, let us highlight that the automatic thresholding method is independent of the numerical method used and the equations solved. Therefore it can be applied for a wide range of numerical applications.

10 Fig. 1

 101 Fig. 1 Comparison of the mean value S m for two given criterion profiles S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ) where x A = 1.25 and x B = 3.75. S represents either a discontinuous flow (left) or a smooth flow (right).

  2(b), the first singularity of d is localised at α = α A = 1 and it corresponds to the first local maximum α A = S(x A ) of S at x A = 1.25 in Fig. 1(b). The second singularity at α = α B = 2 (see Fig. 2(b)) represents the maximum α B = S(x B ) of S at x = x B = 3.75 (see Fig. 1(b)). For discontinuous flows, due to the large peak, a zoom in is necessary to display the local maximum in x A , from the distribution function d, as shown in Fig. 4(a).

10 Fig. 2

 102 Fig. 2 Decreasing rearrangement d(α) of S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ), where x A = 1.25 and x B = 3.75, for discontinuous flow (left) and smooth flow (right).

10 Fig. 3

 103 Fig. 3 Function f of S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ), where x A = 1.25 and x B = 3.75, for discontinuous flow (left) and smooth flow (right).

Fig. 4 10 Fig. 5

 4105 Fig. 4 Zoom in on d(α) and f (α), of S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ), where x A = 1.25 and x B = 3.75, for discontinuous flow.

Fig. 6

 6 Fig. 6 Illustration of the BB-AMR process

  framework: let us denote C 0 the initial cell of length |C 0 | and, C 00 and C 01 the two sub-cells of length |C 00 | and |C 01 |, then step 1 (mass conservation property): For a constant density of water ρ, the mass is conserved if and only if old mesh ρdV = new mesh ρdV ⇐⇒ ρ|C 0 |h 0 = ρ (|C 00 |h 00 + |C 01 |h 01 )

8 m long and 26 . 5 mFig. 9 Fig. 10

 265910 Fig.9Surface profiles of solitary wave propagation in time at wave gauges 1 to 6. Confrontation of experimental data (solid blue line) to numerical data computed on a uniform grid (solid green line) and on an adaptive grid (solid red line)

Fig. 11

 11 Fig. 11 Experimental settings

Fig. 12

 12 Fig. 12 Numerical result for the water height at time t = 0.5s . The colours correspond to the magnitude of the kinetic energy

Fig. 13

 13 Fig. 13 Numerical result for the water height at time t = 2.5s. The colours correspond to the magnitude of the kinetic energy

Fig. 14

 14 Fig. 14 Numerical result for the water height at time t = 5.75s. The colours correspond to the magnitude of the kinetic energy

Fig. 15

 15 Fig. 15 Numerical result for the water height at time t = 23.75s. The colours correspond to the magnitude of the kinetic energy

Fig. 16 Fig. 17

 1617 Fig.[START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF] Free surface results at different positions : experimental data versus numerical simulation with and without mesh adaptivity

Fig. 18 Fig. 19 Fig. 21 Fig. 22

 18192122 Fig. 18 Time evolution of the adaptive mesh : t = 0.5s. The colours correspond to the numerical density of entropy production S

Fig. 23

 23 Fig. 23 Experimental settings

Fig. 24

 24 Fig. 24 Numerical result for the water height at time t = 11.25s. The colours correspond to the magnitude of the kinetic energy

Fig. 25

 25 Fig. 25 Numerical result for the water height at time t = 13.25s. The colours correspond to the magnitude of the kinetic energy

Fig. 26

 26 Fig. 26 Numerical result for the water height at time t = 16s. The colours correspond to the magnitude of the kinetic energy

Fig. 27

 27 Fig. 27 Numerical result for the water height at time t = 17.5s. The colours correspond to the magnitude of the kinetic energy

Fig. 28

 28 Fig.[START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF] Free surface results at different positions : experimental data versus numerical simulation with and without mesh adaptivity

Fig. 29 Fig. 30 Fig. 32 Fig. 33

 29303233 Fig. 29 Adaptive mesh : time t = 11.25s. The colours correspond to the numerical density of entropy production S

•

  The re-numbering and re-meshing being expensive, the mesh is unchanged during a "macro-time" step ∆t n δt n , called AMR time-step, given by the smallest size block |B k | (rather by the smallest size cell |C k | |B k |) and the maximum velocity. The (micro-)time step δt n is the usual time step, computed with the smallest size cell, required to the time advancement to ensure the stability of the scheme.

  9 ,

	(x-25.9) 30.8 2.565 -2.5 if 25.9 < x < 56.7 ,
	0.065	if 56.7 < x < 57.65 ,

(57.95-x) 3 0.2 + 0.0650 if < 56.7 < x < 60.95 , -0.1360 if x > 60.95 .

The Cartesian grids used in this chapter are treacherous because there are often associated to structured mesh. This choice has been done to get a simple meshing tool but is unstructured by construction.

The numerical soft used for the following test cases is the EOLENS code developed by the Institut de Mathématiques de Toulon (IMATH) and Principia.

The authors would thanks V. Roeber who provide the experimental data for this test case.

Several sources and data can be found, see for instance http://nctr.pmel.noaa.gov/ benchmark/Laboratory/Laboratory\_MonaiValley/ or http://isec.nacse. org/workshop/2004_cornell/bmark2.html

sources and pictures are available at http://isec.nacse.org/workshop/2004\ _cornell/bmark2.html

Wave gauges data are available at the address http://nctr.pmel.noaa.gov/ benchmark/Laboratory/Laboratory_MonaiValley/

Acknowledgments

This work is partially supported the Project MTM2011-29306-C01-01 from the MICINN (Spain) and the French national research project TANDEM (Tsunamis in the Atlantic and the English channel : definition of the effects through numerical modeling), the French government (Projets Investissement d'Avenir, agreement reference number ANR-11-RSNR-0023-01). Moreover, the authors wish to thank the referees for their remarks and the careful reading of the numerical method presented in this paper.