Stéphane Clain

C Reis

R Costa

J Figueiredo

M A Baptista

J M Miranda

tsunami simulation with a finite volume solver and the hydrostatic reconstruction technique

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The 28 February 1969 event was a submarine earthquake Ms7.9, with epicentre located on the Horseshoe Abyssal Plain, south-east of the Gorringe bank, approximately the same location as the noticeable Lisbon earthquake and tsunami of 1st November 1755. The earthquake struck western Portugal and Morocco at 02:40:32.5 UTC and a small tsunami was registered in tide stations of mainland Portugal, Azores archipelago, Spain and Morocco [START_REF] Baptista | Maximum entropy analysis of portuguese tsunami data[END_REF]). [START_REF] Baptista | Maximum entropy analysis of portuguese tsunami data[END_REF] [START_REF] Baptista | Maximum entropy analysis of portuguese tsunami data[END_REF]) analysed tsunami data recorded at the tide stations of Portugal mainland, Azores and Cadiz in Spain. Later, [START_REF] Heinrich | Numerical simulation of the 1969 tsunami along the Portuguese coasts[END_REF] [START_REF] Heinrich | Numerical simulation of the 1969 tsunami along the Portuguese coasts[END_REF]) and [START_REF] Gjevik | Modeling tsunamis from earthquake sources near Gorringe Bank southwest of Portugal[END_REF] [START_REF] Gjevik | Modeling tsunamis from earthquake sources near Gorringe Bank southwest of Portugal[END_REF]), presented numerical simulations of the tsunami propagation along the Portuguese coastline and in the Tagus estuary.

Tsunami simulation codes are mainly based on two numerical techniques: the finite difference and the finite volume methods and, to a lesser extend, the finite element method. These solvers tackle the shallow-water equations, in Cartesian or spherical coordinates, equipped with additional source terms to take the bathymetry variation, the friction or the turbulence into account. Finite difference schemes for simulating tsunami propagation were first implemented due to its simplicity in the structured grids context ( MOST Titov and Gonzalez [1997], TUNAMI [START_REF] Imamura | method of tsunami simulation with the leap-frog scheme[END_REF], [START_REF] Imamura | Tsunami modelling manual, UNESCO IOC international training course on Tsunami Numerical Modelling[END_REF] or COMCOT Wang [2009]). Nevertheless, the finite difference schemes suffer of three drawbacks: 1) they are not entirely conservative (the total mass is not preserved, even if the losses can be very small), 2) they are not well-balanced (the configuration ocean at rest is not maintained) and, at last, 3) to reduce the numerical effects of the bathymetry source term, a formulation based on the primitive variables (height and velocity) is often used however does not provide the correct solution when discontinuities are involved. The nonconservative system using the velocity as an unknown function must not be used since only the conservative formulation provides the correct Rankine-Hugoniot conditions. Moreover, second-order versions create some spurious oscillations in the vicinity of discontinuities and a large amount of artificial viscosity is added to stabilise the scheme leading to a dramatic reduction of the accuracy [START_REF] Zhou | Numerical solutions of the shallow water equations with discontinuous bed topography[END_REF]; [START_REF] Gallouët | Some approximate Godunov scheme to compute shallow-water equations with topography[END_REF]; [START_REF] Nikolos | An unstructured node-centered finite volume scheme for shallow water flows with wet/dry fronts over complex topography[END_REF]).

Since the nineties, the finite volume method turns to be a more efficient technique to solve hyperbolic non-linear sys-1 tems for complex modelling in atmospheric sciences and geophysical fluid dynamics (Zhou et al. [2015]) or hydrological models [START_REF] Salah | Tools and Algorithms to Link Horizontal Hydrologic and Vertical Hydrodynamic Models and Provide a Stochastic Modeling Framework[END_REF]). Software for tsunami simulation using the finite volume scheme were developed and implemented such as ANUGA [START_REF] Roberts | ANUGA User Manual[END_REF]) and GEOCLAW [START_REF] Leveque | CLAWPACK software, available from netlib. att. com in netlib[END_REF]), among others. They preserve the mass and the use of conservative variables guarantees the validity of the shock propagation. Non-conservative flux is a key to provide a wellbalanced scheme. GEOCLAW is based on the f-wave decomposition which integrates the bathymetry jump at the cell's interface [START_REF] George | Finite volume methods and adaptive refinement for global tsunami propagation and local inundation[END_REF]) while ANUGA use the Divergence Form for Bed slope source term (tagged DFB) developed in [START_REF] Valiani | Divergence Form for Bed Slope Source Term in Shallow Water Equations[END_REF] to take the hydrological jumps into account.

In this study, we propose to introduce and implement a different finite volume technology which enables to deal with varying bathymetry and provides high-accurate solutions. The hydrostatic reconstruction method [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]; [START_REF] Berthon | Efficient well-balanced hydrostatic upwind schemes for shallow-water equations[END_REF]) became a very popular and simple numerical technique for handling with both the topography and the dry/wet interface. One critical aspect of tsunami simulations is the variation of the bathymetry where the non-conservative term has to be discretized with caution. Indeed, some steady-state situations have to be preserved at the discrete level leading to the socalled C-property [START_REF] Bermúdez | Upwind methods for hyperbolic conservation laws with source terms[END_REF]; [START_REF] Bermúdez | Upwind schemes for the two-dimensional shallow water equations with variable depth using unstructured meshes[END_REF]) and the well-balanced concept [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF]; [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF]). If not preserved, the scheme is disqualified since it produces additional non-physical forces and the solution is not acceptable. At last, in the on-land tsunami propagation the dry/wet situation is a critical point. Here, both the water height h and the mass flow q = hu converge to zero while the velocity u does not vanish. Accurate approximations where h << ∆x (let us say h < ∆x/100), ∆x being the characteristic length of the cell, turn to be difficult since the numerical diffusion takes a larger contribution. It results that the velocity may be unbounded close to the dry zone (and diverges to infinity in some case) leading to a too small time step deriving from the CFL condition. Moreover, the finite volume method is powered with the MUSCL technique (Monotonic Upwind conservative schemes Van Leer [1974] and see the book of [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] for an overview of the method)to improve the scheme accuracy with a very low additional computational cost. The aim of this paper is threefold.

1. The presentation of a fully detailed numerical scheme, based on recent methods developed in the mathematical community for the shallow-water system, including the nonconservative term with a special dedication to the tsunami propagation problem. All the scheme stages are presented: the generic finite volume framework, the MUSCL technique, the C-property leading to the construction of the non-conservative fluxes and, at last, the hydrostatic reconstruction to deal with the dry/wet situations;

2. The verification and validation of the numerical code, where the numerical benchmarks used by [START_REF] Synolakis | Validation and verification of tsunami numerical models[END_REF] and Tinti and Tonini [2013], adopted from the Long-Wave Run-up Models Workshops, edition 1, 2 and 3 [START_REF] Liu | A report on the international workshop on long wave runup[END_REF]; [START_REF] Yeh | Long-wave runup models[END_REF]; [START_REF] Synolakis | Tsunami science before and beyond Boxing Day 2004[END_REF]) were carried out to assess the robustness and the accuracy of the method.

3. The study of an historical real scenario, the 28 February 1969 tsunami. As result, we show that even with the complex bathymetry and reflection characteristics of the Tagus estuary, the code reproduces the registered waveforms recorded by some of the tide stations. We also show the code performance in open-sea propagation, using Lagos (Algarve, Portugal) and Casablanca (Morocco) stations.

Model and numerical scheme

Tsunami modelling involves the shallow-water system equipped with non-conservative terms such as the bathymetry effect, the friction or the turbulence. We skip the two last terms since we focus on the bathymetry term and its numerical treatment. The system writes ∂th + ∂x(hu) + ∂x(hv) = 0

(1)

∂t(hu) + ∂x hu 2 + g 2 h 2 + ∂y(huv) = -gh∂xb (2) ∂t(hv) + ∂x(huv) + ∂y hu 2 + g 2 h 2 = -gh∂yb (3)
where h denotes the water height, (u, v) the fluid velocity vector gathering the components along the x and y axis, b the bathymetry, g = 9.81 ms -2 the gravity acceleration. Moreover, η = h + b stands the free surface, and vector W = (h, hu, hv) gathers the conservative quantities while V = (h, u, v, b) corresponds to the physical variables vector.

Mesh and notations

Let Ωx = [0, L] with L a positive real number. We define the cells ci

• = [x i-1/2 , x i+1/2 ] with interfaces x i+1/2 = i∆x, x i-1/2 = x i+1/2 -∆x and centre xi = x i-1/2 +x i+1/2 2 , i = 1, • • • , I,
where we have set ∆x = L/I. In the same way, the cells c•j = [y j-1/2 , y j+1/2 ] represent a partition of Ωy = [0, H], H > 0 with interfaces y j+1/2 = j∆y, y j-1/2 = y j+1/2 -∆y and centre yj =

y j-1/2 +y j+1/2 2 , j = 1, • • • , J with ∆y = H/J. Domain Ω = [0, L] × [0, H] is decomposed into non- overlapping I × J cells cij = ci• × c•j, i = 1, • • • , I, j = 1, • • • , J, with centroid (xi, yj) and interfaces e i+1/2,j = {x i+1/2 } × c•j, e i-1/2,j = {x i-1/2 } × c•j, e i,j+1/2 = ci• × {y j+1/2 }, e i,j-1/2 = ci• × {y j-1/2 }. For a prescribed time final time T , 0 = t 0 < t 1 < • • • < t n < • • • < t N =
T is a subdivision with non constant time step ∆t n = t n+1 -t n that will be controlled by the CFL condition.

Real number φ n ij represents an approximation of the mean value over cell cij for any function φ = h, η, b, u, v at time t n . We recall that for regular functions over the cell cij (say C 2 ), the point-wise value at (xi, yj) is a second-order approximation of the mean value. We denote by W n ij = (h n ij , (hu) n ij , (hv) n ij ) the vector of the conservative variables and by

V n ij = (h n ij , u n ij , v n ij , b n ij
) the vector of the physical variables where we have set

u n ij = (hu) n ij h n ij , v n ij = (hv) n ij h n ij
when h n ij > 0 and 0 otherwise. Real numbers φ n i+1/2j,L and φ n i+1/2j,R represent approximations on the left and right side of interface e i+1/2j while φ n ij+1/2,L and φ n ij+1/2,R stand for approximations on the left (lower) and right (upper) side of interface e ij+1/2 (see Figure 1).

Non-conservative problem and well-balanced scheme

Due to the change of topography in space, the nonconservative term -gh∇b is required to preserve some steady-state configurations such as the water at rest, where the velocity is null and η is constant in the wet area. Indeed, writing the mass flow equation with null velocity, steadystate assumption yields

∂t(hu) + ∂x hu 2 + g 2 h 2 + ∂y(huv) = gh∂xh = gh∂x(η -b) = -gh∂xb,
and the same property holds for the v component. Consequently, the non-conservative term -gh∇b on the right-side compensates the hydrostatic pressure variation due to the bathymetry variations and the well-balanced property has to be mimic at the discrete level. Discretization of the nonconservative term is still an important challenge to provide some nice properties such as stability and well-balancedness. For example, the simple discretization

-gh∂xb ≈ ghij bi+1,j -bi-1,j 2∆x (4)
is not eligible. Indeed, such an approximation coupled with a classical numerical flux (Rusanov, HLL) will induce nonphysical motions, and after some steps the water moves alone, leading to a wrong simulation. Consequently, preserving the water at rest configuration (the C-property) is the minimal requirement that the numerical scheme has to respect. A scheme is then said wellbalanced (for the water at rest case) if it preserves this specific situation. Other types of steady-state solution can also be preserved (the so-called moving water stationary solutions) and there exists a large literature on this subject in the mathematical community (we refer to [START_REF] Delestre | SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies[END_REF] for an overview of the up-to-date techniques). Nevertheless, the water at rest is the main steady-state one has to preserve in the applications involving lakes or oceans.

To highlight the importance of such a property, we prove hereafter why expression (4) is not correct. Indeed, assume that the lake is at rest at time t n , then we have

u n ij = v n ij = 0 while h n ij +b n ij = η n ij = η ∈ R.
Computing the next time step for the velocity using the first-order finite volume scheme equipped with a viscous flux such as the Rusanov or the HLL one [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics, 3 rd revision[END_REF] leads to

h n+1 ij u n+1 ij = - g∆t 4∆x (h n i+1,j ) 2 -(h n i-1,j ) 2 - -h n ij b n i+1,j -b n i-1,j 2∆x ∆t = - g∆t 4∆x (h n i+1,j -h n i-1,j )(h n i+1,j + h n i-1,j ) - -h n ij b n i+1,j -b n i-1,j 2∆x ∆t = g∆t 4∆x (b n i+1,j -b n i-1,j )(h n i+1,j + h n i-1,j ) - -h n ij b n i+1,j -b n i-1,j 2∆x ∆t = -g∆t b n i+1,j -b n i-1,j 2∆x h n ij - h n i+1,j + h n i-1,j 2 = g∆t b n i+1,j -b n i-1,j 2∆x b n ij - b n i+1,j + b n i-1,j 2 = g∆t b n i+1,j -b n i-1,j 4∆x (2b n ij -b n i+1,j + b n i-1,j ).
Clearly, the right-hand side does not vanish except for the particular case of a constant or linear bathymetry (affine strictly speaking). Hence h n+1 ij u n+1 ij = 0. The steady-state flow is no longer preserved and such a scheme leads to an erroneous evaluation of numerical approximations. On the other hand, we can check that if one employs the following discretization

-gh∂xb ≈ g h n i+1,j + h n i-1,j 2 • b n i+1,j -b n i-1,j 2∆x ,
we obtain a well-balanced scheme for the water at rest situation since the velocity remains null. A major difficulty is that the choice of the discretization of the non-conservative term strongly depends on the discretization of the conserva-tive part and the last choice turns to be inadequate if, for example, we use the Roe numerical flux [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics, 3 rd revision[END_REF]. Therefore the non-conservative contribution has to be adapted in function of the numerical flux used in the conservative term.

MUSCL Technique

Finite volume is an excellent framework due to its built-in conservation property but it suffers of an important numerical diffusion if one uses first-order approximations. Monotonic Upstream-Centred Scheme for Conservation Laws (MUSCL) technique, initially developed by [START_REF] Van Leer | Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second-order scheme[END_REF], is a popular method to easily increase the accuracy while preserving the robustness (see [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] for an overview of the MUSCL method). It is based on two steps: a local linear reconstruction to achieve the secondorder and a limiting procedure for preventing the solution from generating non-physical oscillations.

To compute the first derivatives for any function φ = h, η, u, v, b, we define the slopes for the x and y direction respectively with

p n i+1/2j (φ) = φ n i+1j -φ n ij ∆x , p n ij+1/2 (φ) = φ n ij+1 -φ n ij ∆y ,
and one can achieve a more accurate approximation on edge e i+1/2j taking for example (see Figure 2)

φ n i-1/2j,R = φ n ij - p n i+1/2j + p n i+1/2j 2 ∆x 2 , φ n i+1/2j,L = φ n ij + p n i+1/2j + p n i+1/2j 2 ∆x 2 ,
where we have skipped the reference to φ for the sake of simplicity. Such a reconstruction will give rise to oscillations in the vicinity of a discontinuity due to the Gibbs phenomenon and a non-linear limiting procedure has to be implemented to preserve the monotonicity in each direction. The traditional MUSCL way consists in substituting the original slopes by a limited version π(α, β) depending on the left and right unlimited slope (α, β respectively), such that some stability criterion is fulfilled. Classical limiter operators such as the minmod or the van-Leer limiters are involved in the limiting process [START_REF] Van Leer | Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second-order scheme[END_REF].

Then we define the stabilized reconstructed values on the left and right side of the vertical interfaces with

φ n i+1/2j,L = φ n ij + π p n i-1/2j , p n i+1/2j ∆x/2, φ n i-1/2j,R = φ n ij -π p n i-1/2j , p n i+1/2j
∆x/2, while we set for the horizontal interfaces

φ n ij+1/2,L = φ n ij + π p n ij-1/2 , p n ij+1/2 ∆y/2, φ n ij-1/2,R = φ n ij -π p n ij-1/2 , p n ij+1/2 ∆y/2. Notice that φ n i+1/2j,L , φ n i-1/2j,R
are different and that we recover the constant piecewise representation when π = 0 with

φ n i+1/2j,L = φ n i-1/2j,R = φ n i .
Therefore, the first-order case/method is a possible outcome of the second-order one. In the following W n i+1/2j,L stands for the reconstructed conservative vector on the left side of interface e i+1/2,j and V n i+1/2j,L the corresponding vector using the physical variables.

An important point is that the reconstruction cannot be performed with h, η and b at the same time for compatibility reasons. It has been proved in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] that the good choice is to first carry out the MUSCL procedure on h and η. Then we deduce the values for b setting (notice that the approximations of b now depend on the time):

b n i+1/2j,L = η n i+1/2j,L -h n i+1/2j,L , b n i-1/2j,R = η n i-1/2j,R -h n i-1/2j,R . We evaluate approximations b n ij+1/2,L an b n ij-1/2,R
for the horizontal interfaces in the same way.

Hydrostatic reconstruction

To design a discretization of the source term that preserves the C-property, we use the hydrostatic reconstruction proposed by [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] which enables to deal with complex flows with the dry/wet situation. In particular, such a technique preserves the positivity of the water height providing a very good robustness and accuracy.

Consider a generic interface e and denote bL, hL and bR, hR the bathymetry and water height for the left and right sides of the interface. We set b * = max(bL, bR) and perform the hydrostatic reconstruction setting

h * L = max(0, hL -b * + bL), h * R = max(0, hR -b * + bR).
Therefore h * L and h * R correspond to the water heights which are really involved in the pressure at the interface e. Figure 3 left panel shows that h * R < hR due to the step, the middle panel presents the dry/wet case where h * L = h * R , while the right panel gives an example of hydrostatic reconstruction with a piecewise linear bathymetry.

We now adapt the generic principle to a vertical interface e i+1/2,j for instance. We then denote by b

* ,n i+1/2j = max(b n i+1/2j,L , b n i+1/2j,R
) where b n i+1/2j,L and b n i+1/2j,R are approximations of topography on interface e i+1/2,j and set the new hydrostatic reconstruction variables:

h * ,n i+1/2j,L = max(0, h n i+1/2j,L -b * ,n i+1/2j + b n i+1/2j,L ), η * ,n i+1/2j,L = h * ,n i+1/2j,L + b * ,n i+1/2j , with h n i+1/2j,L , b n i+1/2j
,L an approximations of the water height on the left side of interface e i+1/2,j . We proceed in the same way to compute h * ,n i+1/2j,R and η * ,n i+1/2j,R . For the sake of consistency, we also use the notation

u * ,n i+1/2j,L = u n i+1/2j,L , v * ,n i+1/2j,L = v n i+1/2j,L and u * ,n i+1/2j,R = u n i+1/2j,R , v * ,n
i+1/2j,R = v n i+1/2j,R for the velocity. We apply an identical algorithm for the horizontal interfaces of cell cij to evaluate φ

* ,n i-1/2j,R , φ * ,n i-1/2j,R , φ * ,n ij+1/2,L , φ * ,n ij+1/2j,R , φ * ,n ij-1/2,L , φ * ,n ij-1/2j,R for φ = h, η, u, v, b.
At last, we shall denote by W * ,n i-1/2,L and V * ,n i-1/2,L the conservative and physical vectors after applying the hydrostatic reconstruction.

Second-order in space finite volume scheme

The explicit finite volume scheme writes

U n+1 ij = U n ij -∆t ∆x F n i+1/2j + ε n i+1/2j,L -F n i-1/2j -ε n i-1/2j,R -∆t ∆y F n ij+1/2 + ε n ij+1/2,L -F n ij-1/2 -ε n ij-1/2,R + ∆tS n ij (5)
with

F n i-1/2j = F(W * ,n i-1/2,L , W * ,n i-1/2,R ,
ν) a numerical flux for the conservative contribution that casts under the viscous form

F(WL, WR, ν) = F (WL, ν) + F (WR, ν) 2 -λ(WR -WL),
where

F (W, ν) =   hU • ν huU • ν hvU • ν   + gh 2 2   0 νx νy  
is the physical flux in the normal direction ν and λ > 0 is the scheme viscosity. In the following, the HLL numerical flux will be used in all the numerical simulations. Function

ε n i-1/2j,L = ε h n i-1/2j,L , h * ,n i-1/2j,L
is the nonconservative flux at edge e i-1/2j at the left side of the interface. We recall that, in non-trivial cases such as a flat bottom, the non-conservativity yields

ε n i-1/2j,L + ε n i-1/2j,R = 0,
where this quantity represents the impulsion variation due to the brutal change of topography at the interface. Of course,

ε n i-1/2j,L = ε n i-1/2j,R = 0 if one has b n i-1/2j,L = b n i-1/2j,R
, which corresponds to a local continuous bathymetry. We shall use the following discretization to deal with the discontinuous part of the non-conservative term

ε n i+1/2j,L = g 2   0 ((h n i+1/2j,L ) 2 -h * ,n i+1/2j,L ) 2 0  
and a similar expression for the three other non-conservative contributions, namely ε n i-1/2j,R and ε n ij-1/2,L , ε n ij+1/2,R for the three other interfaces.

At last, when dealing with the second-order approximation, the gradient of the bathymetry has to be substituted with its numerical representation corresponding to the regular part of the non-conservative term. Consequently, the contribution of the source term over the cell writes

S n ij = -g    0 h n i+1/2j,L +h n i-1/2j,R 2 × b n i+1/2j,L -b n i-1/2j,R ∆x h n ij+1/2,L +h n ij-1/2,R 2 × b n ij+1/2,L -b n ij-1/2,R ∆y    .
Notice that the source term for the cell does not involve the hydrostatic reconstructed variable but the original ones. Thus, if one adopts a first-order scheme, we have b n i+1/2j,L = b n i-1/2j,R = bij. Hence the source term in the cell vanishes and the contributions of the bathymetry variations are computed with the non-conservative flux. On the contrary, if the bathymetry is continuous, the non-conservative flux deriving from the bathymetry discontinuity vanishes and the change of topography contribution is exclusively computed with the source term Sij in the cell. To sum up, ε and S aim at computing the non-conservative term but ε is dedicated to the brutal bathymetry variation at the interfaces while S treats of the smooth topography variations inside the cells.

Full second-order finite volume scheme

A second-order method in time is required to guarantee a global second-order method for smooth solutions. The usual RK2 (Heun method) is usually employed to be effective and robust. Assuming that we know all the vectors

W n ij , i = 1, • • • , I, j = 1, • • • , J at time t n ,
we proceed in two sub-steps. We first compute a first-order candidate solution W n,1 ij for time t n + ∆t applying successively the MUSCL procedure, the hydrostatic reconstruction and the finite volume scheme (4). With W n,1 ij in hand, we proceed in the same way computing a second approximation W n,2 ij for time t n + 2∆t. Then the Heun method consists in defining

W n+1 ij = 1 2 W n ij + 1 2 W n,2 ij
which provides a second-order approximation in time at t n+1 . At last, W n+1 ij corresponds to a full second-order approximation both in space and time as long as we respect the CFL condition for the stability.

Numerical tests and validation

Several benchmarks were carried out to perform the verification and validation of the numerical scheme, through comparisons of the code predictions with analytical solutions, laboratory experiments and field measurements.

The validation stage guarantees that the numerical method correctly solves the equations where consistency, accuracy and stability are assessed with representative situations such as bathymetry change, dry/wet interface and run-up on an inclined beach.

To this end, two types of benchmarks are considered: the synthetic benchmark, where the numerical solution is compared to an analytical one, and the laboratory benchmark involving a confrontation of the numerical approximation with data deriving from laboratory experiences.

The benchmarking process used to accomplish the geophysical component is an assumption of the benchmark problems introduced by the Long-Wave Run-up Models Workshops 1990 (Catalina Island in California), 1996 (Friday Harbor in Washington) and 2004 (Catalina Island) [START_REF] Liu | A report on the international workshop on long wave runup[END_REF]; [START_REF] Yeh | Long-wave runup models[END_REF]; [START_REF] Synolakis | Tsunami science before and beyond Boxing Day 2004[END_REF]). Later, as a product of these workshops, a technical memorandum was compiled by the National Oceanic and Atmospheric Administration (Synolakis et al. [2007]) and organized into four different categories: 1) Basic hydrodynamic considerations, including the mass conservation and convergence; 2) Analytical benchmarking, including the single wave on a simple beach problem; 3) Laboratory benchmarking, including tests for solitary wave on a simple beach, solitary wave on a composite beach, solitary wave on a conical island, tsunami run-up onto a complex 3D beach (Monai Valley) and tsunami generation and run-up due to 3D landslide; and 4) Field benchmarking with the Rat Islands and the Okushiri tsunamis.

In our study, the basic hydrodynamic considerations were performed to validate the numerical code. Mass conservation, convergence and stability are the major issues a software has to address. For non-conservative problems, the question of the well-balanced (or C-property) is also important, namely, some steady-state situations have to be preserved. The most popular one is the so-called lake at rest configuration where a constant free surface with varying bathymetry and null velocity at t = 0 will be maintained for t > 0.

For the numerical code verification, three problems were selected to test the code run-up and inundation numerical solutions of the non-linear shallow-water equations, based on [START_REF] Synolakis | Validation and verification of tsunami numerical models[END_REF] and Tinti and Tonini [2013] papers. Two synthetic sanity check benchmarks, namely: tsunami run-up onto slope plane beach benchmark (Carrier and Greenspan [1958]; [START_REF] Carrier | Tsunami run-up and draw-down on a plane beach[END_REF]), and tsunami run-up of a planar surface oscillating in a paraboloidal basin benchmark [START_REF] Thacker | Some exact solutions to the nonlinear shallow-water wave equations[END_REF]). One laboratory benchmark due to the Monai, Okushiri tsunami was carried out to compare the simulation with real experimental data in order to assess the validity of both the model and the numerical method.

Synthetic Benchmarks

Analytical benchmarking is a procedure that is used to identify the dependence of the results on the problem parameters [START_REF] Synolakis | Validation and verification of tsunami numerical models[END_REF]; Tinti and Tonini [2013]).

tsunami run-up onto a sloping plane beach

The uniformly sloping beach one-dimensional benchmark consists in a comparison of a leading-depression N wave approximation with a given initial profile travelling across a 1/10 slope numerical solution with the one obtained by the analytical integral formula given in [START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF]; [START_REF] Carrier | Tsunami run-up and draw-down on a plane beach[END_REF]. The goal of the benchmark is to compute and present the snapshots of the numerical free surface and velocity profiles at t = 160 s, t = 175 s and t = 220 s, and to compare them with the exact solution. We plot in Figure 4 the free surface using the first-order (left panel) and the second-order method (right panel) for different times. In the same way, Figure 5 displays the velocity at the same time instants. We qualitatively observe the convergence of the approximation as ∆x decreases and conclude that the second-order clearly provides the best solution.

To quantify the impact of the second-order method, we report in Table 1 the error of the free surface minimum at t = 175 (exact minimum is -21.34) and the convergence orders. The second-order scheme provides the smallest errors and the convergence rate is around 1.6, whereas the firstorder method presents more limitations to approximate the exact solution. Note that the dry-wet interface is not well collocated since its position depends on whether a cell is wet or dry. We identify this issue as the major problem to provide an accurate solution and it is the main limiting factor of the numerical scheme. Effective second-order error can not be achieved due to this important restriction and a better location of the dry/wet interface inside the cell will be a crucial challenge.

tsunami run-up of a planar surface oscillating in a paraboloidal basin

The second benchmark derives from an analytical solution proposed by [START_REF] Thacker | Some exact solutions to the nonlinear shallow-water wave equations[END_REF] for a parabolic basin and has been used by Tinti & al. Tinti and Tonini [2013] to check the numerical schemes.

The bathymetry is a paraboloid of equation

b(x, y) = b0 1 - x 2 L 2 - y 2 l 2
with L and l the lengths of the ellipse semi-axis. The free surface is a plane that oscillates around the horizontal axis Oy and the component of the velocity along the Ox axis is null. The exact solution has an analytical expression given by:

η(x, y, t) = 2A b0 L cos(ωt) x L - A 2L cos(ωt) ,
for the free surface and

u(x, y, t) = Aω sin(ωt), v(x, y, t) = 0, ω = 2gb0 L 2
for the velocity, where A is free parameter we shall set with the initial conditions. In numerical simulations, we use L = 4700 m, l = 1300 m, b0 = 201.42 m, g = 9.80 ms -2 and A = 235 m leading to ω = 1.3368 10 -2 s -1 and T = 470 s. By construction, η is a constant value at time t = T /4 and t = 3T /4. Then, an easy way to check the code is assessing the flatness defaults by computing the difference between the maximum and the minimum values of η at t = 3T /4. We carry out three simulations using meshes M1 = 250 × 75, M2 = 500 × 150 and M3 = 1000 × 300 and we report in Table 2 the differential which corresponds to the flatness default ∆ f . We get a first-order scheme for the free surface while we report a second-order of accuracy with the MUSCL technique. Moreover, the MUSCL method does not produce any spurious oscillations and the free surface flateness is preserved.

Laboratory Benchmark: Monai

The laboratory benchmark is an 1/400 scale laboratory experiment of the extreme Monai run-up, consequence of the 1993 Okushiri tsunami (village of Monai in Okushiri Island 1993 [START_REF] Matsuyama | An experimental study of the highest run-up height in the 1993 Hokkaido Nansei-oki earthquake tsunami[END_REF]). Laboratory measurements were performed in a tank 205 m long, 6 m deep and 3.4 m wide and three points of interest (PoI) were considered, simulating virtual tide gauges. This PoI's are located at: 1) tide gauge 1: 4.521 m; 1.196 m; 2) tide gauge 2: 4.521 m; 1.696 m and 3) tide gauge 3: 4.521 m; 2.196 m.

We carry out two simulations using the first-order and the second-order version of the code. We display in Figure 7 two elevations corresponding to time t = 14 s (left panel) and t = 16 s (right panel) corresponding to the first-order (top row) and the second-order (bottom row) methods. Clearly, the second-order method provides qualitative better approximations with a sharper front and a larger dry zone which result from a less diffusive scheme.

To better quantify the error between the laboratorial data and the numerical approximations, we compare the water height measurement at three locations and draw the associated numerical curves in Figure 8. We observe the good correspondence until 50 s. From then on, the shallow-water model is not enough representative for larger time due to the lack of dispersive terms. The second-order scheme provides a better approximation, in particular smaller structures (local variations) are captured whereas they do not appear with the first-order scheme.

Case study application Simulation of the 28th February 1969 tsunami

The 28th February 1969 submarine earthquake had a magnitude Ms=7.9, and its epicentre located on the Horseshoe Abyssal Plain, 36.01 o N and 10.57 o W. Its source was interpreted as a thrust fault with a small strike slip component [START_REF] Fukao | Thrust faulting at a lithospheric plate boundary the Portugal earthquake of 1969[END_REF]) (see Figure 9). The earthquake generated a small tsunami that affected the coasts of Portugal, Spain and Morocco [START_REF] Heinrich | Numerical simulation of the 1969 tsunami along the Portuguese coasts[END_REF]; [START_REF] Gjevik | Modeling tsunamis from earthquake sources near Gorringe Bank southwest of Portugal[END_REF]). The tsunami was recorded by several tide gauges, including the Cascais, Pedrouços, Terreiro do Paço, Cabo Ruivo and Lagos, in Portugal, and Casablanca, in Morocco. Figure 9 shows the epicenter of the 1969 earthquake and the loaction of several tide gauges where records are available.

Two alternative fault planes were suggeste by [START_REF] Fukao | Thrust faulting at a lithospheric plate boundary the Portugal earthquake of 1969[END_REF] and [START_REF] Matias | Tsunamigenic earthquakes in the Gulf of Cadiz: fault model and recurrence[END_REF]. The strike angle is 55 o for the case of the NE-SW dipping fault plane and 235 o for the SW-NE dipping fault plane. The dip angle is 52 o and the rake 90 o . Considering a seismic moment of M0 = 6.10 20 N.m, the dimensions of the fault plane were fixed as 80 km long by 50 km width Heinrich, [START_REF] Heinrich | Numerical simulation of the 1969 tsunami along the Portuguese coasts[END_REF].

We carried out several simulations of the tsunami generation and propagation, and compared them with the recorded waveforms. The initial time of the simulation, t=0, is the time of the earthquake. We ran the simuation for 180 minutes. All tide records were digitized, linearly interpolated and detided to isolate the tsunami signal. To remove the tide we used a polynomial fitting of degree 9. The simulations were carried out for Mean Sea Level conditions. The bathymetric grid is an uniform spacing of 200 m, and was generated from the General Bathymetric Chart of the Oceans (GEBCO, http://www.gebco.net/).

To prescribe the initial conditions, we assumed the earthquake rupture to be instantaneous and the water incompressible, so that the initial sea surface displacement mimics the coseismic deformation of the sea bottom. We computed the co-seismic deformation using the half-space elastic approach (Okada [1985]), implemented in Mirone suite [START_REF] Luis | Mirone: A multi-purpose tool for exploring grid data[END_REF]. The initial sea surface displacement is shown in Figure 10 for the SW-NE dipping fault plane (left panel) and for the NE-SW dipping fault plane (right panel). The maximum vertical displacement is 1.65 m. Null velocity is assumed at the initial time. Transmission conditions are prescribed on the whole boundary to make the waves freely flow out of the domain without spurious oscillations.

Numerical simulations were carried out with the firstand the second-order method and the two fault plane solutions. The comparison between the synthetic results and the recorded data focused on four aspects: 1) amplitude of the first wave, 2) tsunami travel time, 3) frequency of the first waves, and 4) first wave polarity.

Wave amplitude

Figure 12 displays the comparison of the water elevation at the different locations with respect to the time, while Table 4 provides the maximum water height values for the recorded and simulated waves.

The main difficulty to fit the simulation with the recorded data is the grid size, which prevents to take into account structures of size lower than 200 m, and to prescribe an accurate location of the sensors. We performed the simulations using as Points of Interest (PoI) virtual sensors situated in the neighbor cells of the real sensor localization, ensuring that synthetic waveforms were representative of the area. We focused on the three first tsunami waves since local geographical configurations produce new reflection waves that are latter on superposed on the waveform, and are not well reproduced by the numerical simulation.

• The Cascais PoI synthetic waveform shows an amplitude about 11 cm, 26 cm, 29 cm and 24 cm larger than the observed waveform (46.7 cm), for the different method orders and fault orientation, respectively. Such differences may derive from the large diffusivity inherent to the low order scheme which dramaticaly reduces the accuracy (a third of the water height). Since the gravitational wave is compressed and focused as tsunami enters shallower waters close to the coasts, the high diffusion scheme is not able to correctly reproduce the energy concentration and provides erroneous estimates of the tsunami impact. Second-order scheme is, at least, accurate enough to give a relevant shape of the wave when reaching the coast.

• The Pedrouços tide gauge simulation shows an amplitude 10 cm smaller in the first-order scheme and strike of 55 o . For the remaining synthetic solutions, we observe an over-estimation of about 19 cm, 20 cm and 16 cm.

• The Terreiro do Paço comparison shows that the first wave peak amplitude is 19 cm and 46 cm, respectively for the first-order, depending on the fault orientation. For the second order method, the values are 44 cm and 40 cm. The recorded waveform shows an amplitude of 42 cm.

• In Lagos there is good agreement between the amplitude of the recorded waveform and the modelling results. The first-order simulation presents differences of 7 mm and 4.8 cm. The second-order simulation shows differences of 8.4 cm and 4.8 cm.

• The Casablanca second-order simulation, used mainly to analyse the polarity, shows differences of about 15 cm and 4 cm.

Second-order technique achieves a better approximation in terms of values and form of the wave. As a pattern, the first-order simulation results are smaller than the secondorder results. In what concerns the choice of the fault plane, the SW-NE solution behaves sligthly better than the NE-SW.

Tsunami travel time

Tsunami travel times are generally well reproduced by the model, taking into account the tsunami propagation between the PoI and the real tide gauge, which is not reproduced by the numerical model.

• The Cascais tide gauge recorded the peak of the first wave about 36 min after the earthquake while the first-order simulation gives a travel time of about 40 min and 41 min and the second-order simulation 41 min and 39 min, depending on the polarity of the initial condition.

• The Pedrouços travel time is 66 min (±4.4 min) while the numerical simulations give travel times of 51 min, 53 min 53 min and 50 min, respectively for the first-and secondorders (NE-SW and SW-NE fault planes).

• The Terreiro do Paço observed travel time is 72.6 min (±4.6 min) while the simulation gives a travel time of 63 min, 50 min for the first-order and 65 min, 63 min for the secondorder.

• In Lagos PoI, the synthetic travel time is 40.7 min, while it is about 35 min, 33 min , 32 min and 35 min, respectively for the first-and second-orders (NE-SW and SW-NE fault planes).

• In Casablanca, the arrival time is not coherent with the recorded waveform. This situation was already noted by previous studies [START_REF] Guesmia | Numerical simulation of the 1969 Portuguese tsunami by a finite element method[END_REF]).

In all cases, the numerical simulation provides smaller arrival times than the ones obtained from the recorded data. This is consistent with the bathymetric resolution and the very shallow area where the sensor are located, with a depth shallower than 2-4m. The situation is particularly imporant for the tide gauges located on the shallow margins of the Tejo river. Nevertheless, second-order simulations reproduce better the recorded waveforms, as expected. Figure 11 presents the wave propagation over time..

Period

Tsunami propagation results into a series of waves due to the reflections and the refraction with the coast leading to a complexe waveform constituted of the superposition of several frequencies. Nevertheless, one can observe that the lower frequency content is almost constant (at least for the three first periods) in all tide gauges, and mainly a function of the tsunami source size. To do a quantitative comparison we estimated the largest period by the delay between the first and the second peak and from the average of the three peaks for both observed and synthetic waveforms (see Figure 12).

• The Cascais observed waveform is probably affected by the proximity with marine strutures, resulting in three main superposed frequencies. It is possible to identify the high and low frequencies and perform the comparison with the simulated waveforms.

• Pedouços recorded waveform also shows a complex frequency content (see Figure 12). The numerical code reproduces fairly the lowest frequency.

• In Terreiro do Paço, the period of the first wave is better reproduced by SW-NE fault simulation, while the frequency estimated from the average of the two initial waves is better reproduced by NE-SW fault simulation.

• In Lagos the recorded waveform has higher frequencies than all simulated waveforms.

• In Casablanca both observed and simulated periods are similar.

Choice of the Fault Plane

Due to the location of most of the tide gauges with respect to the tsunami sources, it is difficult to judge which is the best fault plane from the polarity of the first tsunami wave and we got inconlusive results. To get a better insigth, we focus the polarity analysis on the tide gauge of Casablanca which, due to its location relatively to the earthquake source, constitutes the best PoI to evaluate the polarity. Results (see Figure 12) show that the SW-NE dipping fault plane fits better the tsunami waveform.

Conclusions

From the numerical point of view, the finite volume technique coupling the hydrostatic reconstruction and the second-order scheme MUSCL technique provides an efficient code with mass conservation and ocean at rest preservation. Shocks are well-evaluated and the MUSCL technique reduces the numerical diffusion and increases the accuracy without generating non-physical oscillations. The code was submitted to 3 different benchmarks to assess the accuracy, the robustess and the C-property. The benchmarks were performed with the first-and second-order schemes for the run-up and inundation numerical solutions of the non-linear shallow water equations. Comparisons between the synthetic and the real waveforms show an very good agreement.

Finally, we test the code with a real case, the propagation of the small tsunami generated by the 1969 tsunamigenic earthquake (Ms7.9) in the Tagus estuary. Although the complex bathymetry, the numerical code had a favourable performance and the comparison synthetic/recorded waveforms shows a good agreement. We also use the numerical simulation to ascertain the earthquake focal mechanism proposed by [START_REF] Fukao | Thrust faulting at a lithospheric plate boundary the Portugal earthquake of 1969[END_REF] and challenged by [START_REF] Matias | Tsunamigenic earthquakes in the Gulf of Cadiz: fault model and recurrence[END_REF] in what concerns the choice of the fault plane. The comparison between mareographic measurements and synthetic tsunami waveforns favours the fault plane dipping SW-NE as proposed by [START_REF] Fukao | Thrust faulting at a lithospheric plate boundary the Portugal earthquake of 1969[END_REF]. As expected,the second order scheme code generates waveforms with a better approximation to the real waveforms recorded by the considered tide stations.

The main conclusions of this paper are:

• Volume finite method is adequate for the shallow water equations resolution;

• The implemented techniques increase the performance of the method;

• Mathematical and geophysical benchmarks (analytical solutions, laboratory experiments and field measures) certifies the numerical code capability to simulate tsunami runup and propagation with accuracy and robustness;

• Tsunami data favours the focal plane of the 1969, 28th February earthquake, proposed by [START_REF] Fukao | Thrust faulting at a lithospheric plate boundary the Portugal earthquake of 1969[END_REF].

Future developments concerns the introduction of a more recent limiting technique, namely the MOOD method, in substitution to the MUSCL one [START_REF] Figueiredo | Second-order finite volume MOOD method for the shallow water with dry/wet interface SYM-COMP[END_REF], which provides effective second-order of accuracy and preserves some essential properties such as the positivity of the water height. x i x i+1/2 p i+1/2

φ i+1/2,L φ i φ i+1 x i+1 x i-1/2 x i-1 p i-1/2 φ i-1 φ i-1/2,R

Figure 1 .

 1 Figure 1. Notations and grid.

Figure 2 .

 2 Figure2. Slope and interface value reconstructions for the Ox direction. We skip the second index j for the sake of clarity.

Figure 4 .

 4 Figure 4. Free surface at t = 160 (top), t = 175 (middle) e t = 220 seconds(bottom) for the run-up on linear ramp. the left panel corresponds to the first-order scheme and the right panel to the second-order scheme.

Figure 5 .

 5 Figure 5. Water velocity at t = 160 (top), t = 175 (middle) e t = 220 seconds(bottom) for the run-up on linear ramp. the left panel corresponds to the first-order scheme and the right panel to the second-order scheme.

Figure 6 .

 6 Figure 6. Free surface of the paraboloid basin at time t = 3T /4 corresponding to the exact solution η = 0. The free surface is situated between -0.003 (blue) and 0.003 (red).

Figure 7 .

 7 Figure 7. Comparison of the first-(top) and secondorder (bottom) scheme at time t = 14 s (left panel) and t = 16 s (right panel).

Figure 8 .

 8 Figure 8. Comparison of laboratorial and numerical data, for three tide gauges located at points (4.521, 1.196) (tide gage 1), 4.521, 1.696 (tide gage 2) and 4.521, 2.196 (tide gage 3).

Figure 9 .

 9 Figure 9. Location of Gorringe Bank. Epicenter of the 1969 earthquake and respective fault (parallel to the Horseshoe fault), represented by a star. Location of the tide gauge stations in Casablanca, Lagos, Cascais and in Tagus estuary, represented by triangles (left). Zoom of the Cascais and Tagus estuary tide stations location (right).

Figure 10 .

 10 Figure 10. Initial deformation, considering a NE-SW (left) and SW-NE (right) fault orientation. Location of the simulation PoI's.

Figure 11 .

 11 Figure 11. Tsunami wave propagation near and in the Tagus estuary, for time t = 40min, t = 60min, and the restant with time intervals of 10min.

Figure 12 .

 12 Figure 12. Tsunami waves amplitude (m), for the Cascais, Pedrouços, Terreiro do Paço, Lagos and Casablanca tide gauges. Comparison of the first and second order numerical simulation with the tide gauge waveform.

Table 1 .

 1 Errors and convergence rates for the free surface minimum at time t = 175 s.

		∆x error 1st order error 2nd order
		50	7.30	-	4.00	-
	min t = 175	20 10	3.22 1.25	0.9 1.3	0.98 0.31	1.5 1.7
		5	0.58	1.1	0.10	1.6

Table 2 .

 2 Flatness defaults at time t = 175 for several meshes and orders.

	Mesh	∆ f : 1st-order ∆ f : 2nd-order
	250 × 75	0.82	0.072
	500 × 150	0.44	0.021
	1000 × 300	0.21	0.006

Table 3 .

 3 PoIs used in the 28 February 1969 event with respective location.

	Station	Location
		Latitude Longitude
	Cascais	38.693	9.411
	Pedrouços	38.690	9.259
	Terreiro do Paço	38.704	9.136
	Lagos	37.066	8.667
	Casablanca	33.610	7.386

Table 4 .

 4 Amplitude: comparison between the recorded waveform and the synthetic results obtained from a first and second order simulation. Station Recorded [m] 1st-order, strike 55 o [m] 1st-order, 235 o [m] 2nd-order, 55 o [m] 2nd-order, 235 o

	Cascais	0.467	0.577	0.732	0.759	0.703
	Pedrouços	0.433	0.334	0.627	0.633	0.594
	Terreiro do Paço	0.417	0.192	0.456	0.436	0.403
	Lagos	0.429	0.436	0.477	0.513	0.477
	Casablanca	0.622	-	-	0.476	0.582

Table 5 .

 5 Arrival time: comparison between the recorded waveform (Red.), the synthetic results obtained from a first and second order simulation and a tsunami travel time estimative from Mirone code (Mir.) at the northeast and southwestern points of the fault.

	y j+1/2	e ij+1/2			φ ij+1/2,R
	y j-1/2	y j	e i+1/2j	φ ij	φ ij+1/2,L	φ i+1/2j,L φ i+1/2j,R
		x i
		x i-1/2	x i+1/2

Acknowledgments. Historical data for Cascais and Lagos (1969 Lisbon Tsunami) are available at http://www.dgterritorio.pt/ca

The tagus estuary data (typewriter document) are available at the dom Luiz institute library http://idl.ul.pt/node/33. This work is funded by the Portugal-France research agreement, through the research project GEONUM FCT-ANR/MAT-NAN/0122/2012.