A G Filippini
email: andrea.filippini@inria.fr

M Ricchiuto
email: mario.ricchiuto@inria.fr

M Kazolea
email: maria.kazolea@inria.fr

A Flexible 2D Nonlinear Approach for Nonlinear Wave Propagation, Breaking and Run up

Keywords: Green-Nagdhi equations, unstructured, Finite Volume, Finite Element, breaking equations Green-Nagdhi, maillages non-structurée, méthode des volumes finis, méthode des éléments finis, déférlement

We present a hybrid solution strategy for the numerical solution of the twodimensional (2D) partial diferential equations of Green-Nagdhi (GN), which simulates fully nonlinear, weakly dispersive free surface waves. We re-write the standard form of the equations by splitting the original system in its elliptic and hyperbolic parts, through the definition of a new variable, accounting for the dispersive effects and having the role of a non-hydrostatic pressure gradient in the shallow water equations. We consider a two-step solution procedure. In the first step we compute a source term by inverting the elliptic coercive operator associated to the dispersive effects; then in a hyperbolic step we evolve the flow variables by using the non-linear shallow water equations, with all non-hydrostatic effects accounted by the source computed in the elliptic phase. The advantages of this procedure are firstly that the GN equations are used for propagation and shoaling, while locally reverting to the nonlinear shallow water equations to model energy dissipation in breaking regions. Secondly and from the numerical point of view, this strategy allows each step to be solved with an appropriate numerical method on arbitrary unstructured meshes. We propose a hybrid finite element (FE) finite volume (FV) scheme, where the elliptic part of the system is discretized by means of the continuous Galerkin FE method and the hyperbolic part is discretized using a third-order node-centred finite volume (FV) technique. The performance of the numerical model obtained is extensively validated against experimental measurements from a series of relevant benchmark problems.

Introduction

Accurate simulations of water wave's propagation and non-linear wave transformations is of fundamental importance to marine and coastal engineering. Over the last decades, significant efforts in the development of depth averaged models have been made in order to provide the means of accurately predicting near-shore wave processes such as shoaling and runup, diffraction, refraction and harmonic interaction. One of the most applied depth averaged models is the Non-linear Shallow Water Equations (NSWE) which are able to model important aspects of wave propagation phenomena, the general characteristics of the ru-nup process, and the wave breaking with broken waves represented as shocks, but they are not appropriate for deeper waters and shoaling since they neglect all the dispersive effects that play a very important role. In order to take dispersive effects in to account we must keep the O(µ 2 ) terms from the full water waves equations, which where neglected in the derivation of the NSWE. µ is the shallowness parameter defined as water depth to wavelength ratio h 0 /L. This leads to the Green-Naghdi (GN) equations [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF] known also as Serre equations. The range of validity of the model may vary as much as far the non-linearity parameter (defined as the ratio of wave amplitude to water depth A/h 0 ) is concerned, but it requires the shallowness parameter µ to be small.

In this work we use the improved GN system of equations in the form proposed in Bonetton et al. (2011). This formulation has been recovered by adding some terms of O(µ 2 ) to the momentum equation in order to improve the frequency dispersion description of the original GN model, using a tuning parameter α. The two dimensional form of the system can be written in the following form:

h t + ∇ • (hu) = 0
(1)

(I + αT ) q t + ∇ • q ⊗ q h + gh∇η + g α h∇η + hQ 1 (u) = 0
where the operators T (•) and Q(•) are:

T (•) = - 1 3 ∇ h 3 ∇ • (•) h - h 2 2 ∇ • (•) h ∇b + 1 2 ∇ h 2 ∇b • (•) h + h ∇b • (•) h ∇b, (2)
Q(•) = 2 3h ∇ h 3 ∇(•) 1 • ∇ ⊥ (•) 2 + (∇ • (•)) 2 + h 2 ∇(•) 1 • ∇ ⊥ (•) 2 + (∇ • (•)) 2 ∇b + + 1 2h ∇ h 2 (•) • (•) • ∇ ∇b + (•) • (•) • ∇ ∇b ∇b. (3)
We denote h(x, t) = h 0 + η(x, t) -b(x) the total water depth, where η(x, t) the free surface elevation with respect to the water rest state h 0 , b(x) the topography variation and u(x, t) the flow velocity. (•) 1 and (•) 2 indicates respectively the first and second component of the vector (•) and ∇ ⊥ states for the normal gradient operator.

The operator T (•) plays a key role, as its inversion is necessary to be able to obtain evolution equations for the physical variables. T (•) can be written in compact form involving two operators S 1 (•), S 2 (•) and their adjoints S 1 (•), S 2 (•), as:

T (•) = S 1 hS 1 (•) h + S 2 hS 2 (•) h (4)
where

S 1 (•) = h √ 3 ∇ • (•) - √ 3 2 ∇b • (•), S 2 = 1 2 ∇ • (•). (5)
Note that this formulation is essential to show the coercivity of the operator (I +αT ) see [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF] and referenced therein, for further details. Note also that when α = 1 we retrieve the original GN equations.

Discretization strategy

To numerically solve (1), we rewrite the system of the two dimensional enhanced GN equations as:

(I + αT ) φ = W -R, (6)
h t + ∇ • q = 0, (7)
q t + ∇ • q ⊗ q h + gh∇η = φ (8)
by splitting the original system in its elliptic and hyperbolic parts, through the definition of the new variable φ = [φ x , φ y ] T . Φ accounts for the dispersive effects and has the role of a nonhydrostatic pressure gradient in the Shallow water equations. W = T (gh∇η) and R = hQ(u).

In this work we solve ( 6)-(8) using a hybrid Finite Element (FE)-Finite volume (FV) scheme where the elliptic part of the system is discretized by means of the continuous Galerkin FE method. The hyperbolic part of the system is discretized by the two dimensional node centred FV formulation used in [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF].

Inria

Elliptic Phase

We can prove the coercivity of the operator (I + αT ) [START_REF] Filippini | Free surface flow simulation in estuarine and coastal environments: numerical development and application on unstructured mehes[END_REF] via the variational form of T . This property descends on the self-adjoint nature of the operator T and is of primary importance to insure the invertibility of the matrix derived from the discretization (I + αT ), which is a necessary condition for the discrete equation to be solved.

The result of the elliptic equation will not be used inside the hyperbolic part in nodes where h i ≤ wd h ( wd h is a very small threshold value), and simply φ i = (0 , 0) T is set, so we can define a specific (per unit depth) dispersive correction. ψ such that: φ = hψ.

We can, thus, formulate the elliptic step using the unknown ψ. The definitions of the operator T and of S 1 (•) and S 2 (•), lead to the following variational form of the elliptic equation:

Ω h ν • ψ + α Ω S 1 (ν) h S 1 (ψ) + α Ω S 2 (ν) h S 2 (ψ) = RHS , (9)
being RHS a compact form to write the variational formulation of the right hand side of the equation, and being ν the vector of components respectively (ϕ i , 0) T or (0 , ϕ i ) T (with ϕ i the standard linear Lagrange basis functions), if the equation in the x or y direction is taken into account. ψ i = (0 , 0) T is imposed when h i ≤ wd h . The use of ψ as an unknown in the FE phase leads to a symmetric positive definite bilinear form on the left hand side of the variational equation ( 9), namely for the matrix of the differential operator (h + αT ). The elliptic part of the system is, thus, solved for ψ. The value of φ, to use in the conservative hyperbolic equations, are then recovered nodally a posteriori by reverting the definition φ = hψ

Let Ω h denote, an unstructured triangulation of the spatial domain Ω, with h the reference element size and with K denoting the generic element of area |K|. For a node i ∈ Ω h , let K i is the set of elements containing the node i. We consider piecewise linear continuous approximations η h and q h with standard piecewise linear continuous Lagrange bases.

Denoting with (•) x and (•) y the component of the vectorial quantity (•) along the x-axis and y-axis respectively, we define Ψ = [Ψ x , Ψ y ] T , with Ψ x = [ψ x 1 (t), ψ x 2 (t), ..., ψ x N (t)] T and Ψ y = [ψ y 1 (t), ψ y 2 (t), ..., ψ y N (t)] T , and

U = [U, V ] T , being U = [u 1 (t), u 2 (t), ..., u N (t)] T and V = [v 1 (t), v 2 (t), ..., v N (t)] T .
The discrete form of the elliptic equation (6) will thus read:

(M G H + αT)Ψ = W -R , (10)
W = T δ , (11)
R = Q (h, U) , (12)
with δ an approximation of g∇η. The matrix M G H is the four block matrix with zero off-diagonal terms and two equal blocks M G H on the main diagonal.

M G H = M G H 0 0 M G H ,
where M G H is the symmetric positive definite matrix whose entries are represented by M G H i,j =

Ω h h ϕ i ϕ j .
The operator T(h h , b h ) is a matrix differential operator of order two, acting on two dimensional vectors, which can be written as:

T = T 1 1 T 1 2 T 2 1 T 2 2 .
We call T m,n i,j

the element (i, j) belonging to the block (m, n) of the matrix T. Its fully discrete expression is obtained evaluating, with approximate numerical quadrature over each mesh element and the hypothesis of piecewise linear variations of all the quantities involved, the following series of integrals (obtained by the assumption of periodic boundary conditions):

T m,n i,j (h h , b h ) = 1 3 Ω h ∂ Xm (ϕ i ) h 3 h ∂ Xn (ϕ j ) - 1 2 Ω h ϕ i h h ∂ Xn (ϕ j ) ∂ Xm (b h ) + - 1 2 Ω h ∂ Xm (ϕ i ) h 2 h ∂ Xn (b h ) ϕ j + Ω h ϕ i h h ∂ Xn (b h ) ϕ j ∂ Xm (b h ) ,
where {X} = {x, y} and m = 1, 2 and n = 1, 2.

Proceeding similarly, the full discrete form for the operator Q[h h , u h , b h , δ bh ] is found by evaluating, with a quadrature approximation over each triangle of the mesh and with the hypothesis of piecewise linear variation of the quantities involved, the following series of integrals:

Q[h h , u h , b h , δ bh ] = - 2 3 Ω h ∇ϕ i h 3 h ∇u h • ∇ ⊥ v h + ∇ • u h 2 + + Ω h ϕ i h 2 h ∇u h • ∇ ⊥ v h + ∇ • u h 2 ∇b h + - 1 2 Ω h ∇ϕ i h 2 h u 2 h ∂ x δ x bh + v 2 h ∂ y δ y b h + 2 u h v h ∂ x δ y bh + + Ω h ϕ i h u 2 h ∂ x δ x bh + v 2 h ∂ y δ y b h + 2 u h v h ∂ x δ y bh ∇b h . (13)
Following the linear dispersion analysis performed on the one-dimensional scheme in [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF] and of the resulting optimum configuration that have been found to minimize the dispersion error, we use non-lumped mass matrix in equation ( 10) and for the reconstruction of the auxiliary variables δ h and δ bh .

M G δ h i = Ω h ϕ i g h h ∇η h and M G δ bh i = Ω h ϕ i ∇b h .
Note that the additional linear system is solved just once at the beginning of the simulation. Finally, the complete discrete linear system is solved by making use of the standard functionalities of the MUMPS algebraic library to factorize the (M G H + αT) matrix in any time step of the simulation.

Hyperbolic Phase

For the discretization of the hyperbolic part of the system ( 6)-( 8) we use the node-centered FV technique developed and validated in [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF]. In this finite volume framework, the triangulation of Ω used for the finite element discretization of the elliptic phase provides what we refer to as mesh. In the node-centered discretization, a median-dual partition is used to generate non-overlapping control volumes, covering the entire computational domain. We Inria denote by C i the median dual cell obtained by joining the gravity centers of the triangles in K i with the midpoints of the edges meeting in the node i, being K i the set of triangles of the mesh containing i . We can thus write:

|C i | = K∈Ki |K| 3
. The interface belonging to adjacent nodes i, j is named ∂C ij and is composed by the union of two segments connecting the barycenters of the two triangles satisfying K i, j with the midpoint of the edge ij . The boundary of the median dual cell of i can thus be defined as: ∂C i = j∈Ki ∂C ij . Moreover, we define r ij the vector connecting nodes i and j, while the normal and the area associated to the interface ij are:

n ij = 1 2 K i,j n K ij , |C ij | = K i,j |C K ij | with |C K ij | = |K| 6
For simplicity we rewrite the system of conservation laws ( 7)-( 8) as

U t + ∇ • H(U) = S b + Φ. ( 14
)
The FV integration over each computational cell C i leads to the semi-discrete form of the scheme as:

∂U i ∂t + 1 |C i | j∈Ki ∂Cij Fn x + Gn y = 1 |C i | j∈Ki Cij S b + Φ , (15)
where Φ = Cij Φ, is computed by exact integration over C i of the piecewise linear polynomial φ h , obtained from the elliptic phase. n = [n x , ny ] T is the unitary outward vector normal to the boundary of the computational domain and U i is the volume averaged value of U at a given time. F and G are the numerical flux vectors across each internal and boundary face. The numerical fluxes are evaluated solving a Riemann problem at cell interfaces using the approximate Riemann solver of Roe (1981). To reach higher-order spatial accuracy an extension of the MUSCL methodology of [START_REF] Van Leer | Towards the ultimate conservative difference scheme III. Upstream Centered Finite Difference Schemes for Ideal Compressible Flow[END_REF] is used. Each component of the physical variables and bed topography, b, is extrapolated using solution gradients obtained using a combination of centered and upwind gradients. In this way a third-order well-balanced scheme is obtained.

To obtain a well-balanced FV scheme, an upwind discretization approach for the bed topography source term is adopted to satisfy the so-called C-property in hydrostatic (flow at rest) conditions To this end, the topography source term, S b must be linearized in the same way and evaluated in the same Roe-average states as the flux terms. More details on this and on the wet/dry fronts special considerations that are needed t to accurately model transition between wet and dry areas and maintain the high-order spatial accuracy and mass conservation can be found in [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF] and Filippini (2016).

Interface with FE

The FV scheme evolves the average solution U i over the dual cells C i but the FE elliptic phase needs to know the values of the variables in the nodes os the mesh. When the mesh is symmetric the two values coincide. This is not the case when unstructured meshes are used. The nodal values are thus reconstructed at every time step as:

w h i = w i -(∇w) i • r i (16)
where r i is the distance between node i and the center of gravity of the dual cell C i .

RR n°9013

3 Time integration and boundary conditions

Concerning the time dicretization the third order strong stability preserving (SSP) Runge-Kutta (RK) scheme was adopted under the usual CFL stability restiction. In proximity of dry areas and in general when the water depth is low, the friction effect must be taken in to account. The friction term τ = -C f uu/2 is added to the momentum equation as a source term. C f can be a constant number or it can be define by a low which account for the local bottom steepness and water depth. The friction term dicretized in a semi-implict approach, since the explicit treatment of this term may impose a severe time step restriction in the presence of dry areas. We use the semi-implicit approach of Chertock et al. ( 2015) which conserves the high-order accuracy of the initial scheme when friction becomes dominating.

In this work only fully reflective and absorbing boundary conditions have been used. φ n = 0 , ∂ n φ t = 0 . While the first condition simply seems reasonable and can be derived from the momentum equation, considering that u • n = 0, the second one, which forces the tangent component of φ to be constant in the normal direction with respect to the wall, seems artificial. Note that this condition can be obtained by requiring that ∇ ∧ Φ on the wall boundaries (with ∇∧ indicating the curl operator), and using the fact that ∂ t (φ • n) = 0 along these boundaries. The imposition is made by directly modifying the right hand side of (10) and the related entries of the matrix operator (M G H + αT). A strong imposition of the reflective wall boundary is also applied in the hyperbolic phase through the condition: u • n = 0. Absorbing boundaries are also applied in order to dissipate completely the energy of the incoming waves, trying at the same time to eliminate any non-physical reflection. This kind of boundaries requires the definition of a sponge layer in which the surface elevation and the momentum are damped by multiplying their values by a coefficient (see [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF] and [START_REF] Filippini | Free surface flow simulation in estuarine and coastal environments: numerical development and application on unstructured mehes[END_REF] for details).

Wave breaking

A hybrid strategy for wave breaking treatment is implemented in the scheme. We first estimate the location of breaking waves using explicit criteria, applying the NLSW equations to solve the flow in the flagged cells and the GN ones elsewhere. Following the work of [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq type equations[END_REF] we use the combination of the two above phase-resolving criteria for the triggering mechanism The values of γ and φ c are depending on the type of the breaker. The first criterion flags for breaking when η t is positive, since breaking starts on the front face of the wave, while the second criterion, acting complementary to the first, is useful for the detection of hydraulic jumps. In this work the value of φ c = 30 o is used. Moreover, the estimation of the Froude number of the wave is used to established when to switch of the breaking and to detect non-breaking bores. A practical implementation of the breaking mechanism can be found in [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq type equations[END_REF] and [START_REF] Filippini | Free surface flow simulation in estuarine and coastal environments: numerical development and application on unstructured mehes[END_REF].

Numerical Results

Inria

2D solitary wave propagation

The accuracy of the two-dimensional scheme proposed for the GN equations is verified by performing a convergence analysis on a solitary wave propagation problem. A solitary wave of amplitude 0.2m is left propagating for t = 1s inside a computational domain of [0, 70] × [0, 0.8]m, characterized by a value of still water depth of h 0 = 1m In order to measure the rate of convergence of the numerical solution to the exact one, we compute the relative error on the total water depth

E L2 (h) = ||h num -h ex || 2 /||h ex || 2
, where h num is the numerical solution and h ex is the analytical one. The numerical model used triangular grids consisting of triangles with side length of h K = [0.2, 0.1, 0.05, 0.025]m. Figure 1 shows the results obtained, together with the slopes 2.5 and 3 as references, confirming the expected order of accuracy for the numerical scheme.

log 10 (h K )

Wave diffraction over a semi-circular shoal

We consider here the reproduction of the tests carried out in [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF] involving the study of the focusing eeffect induced by a semicircular shoal on wave trains of different periods. The experiments carried out in a wave tank 6.096m wide and 25.6m long, its middle portion consisted in a semicircular shoal leading the water depth to decrease from h 0 = 0.4572m (at the wave maker) to 0.1524m at the end of the tank. The bottom topography is described by the equation:

z(x, y) =      0 x ≤ 10
.67 -Λ(y), (10.67 -Λ(y) -x)/25 10.67 -Λ(y) < x < 18.29 -Λ(y), 0.3048

x ≥ 18.29 -Λ(y),

where Λ(y) = 6.096y -y 2 . The depth h is obtained as h = h 0 -z. This test has become a standard test case for 2D dispersive numerical models to test nonlinear refraction and diffraction, we refer for example [START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF], [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF], [START_REF] Ricchiuto | Upwind residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF], [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF] among others.

Two cases are considered here with wave trains characterised by [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF] (a) T = 1s, A = 0.0195m, h 0 /λ = 0.306 ;

(b) T = 2s, A = 0.0075m, h 0 /λ = 0.117 ;

where T denotes the period, A the wave's amplitude and λ the wavelength. For both cases, the harmonic analysis of free surface elevation measurements taken along the tank centreline are available, and are used to verify the capabilities of a model to reproduce nonlinear refraction and diffraction. The computational domain is [-10, 36] × [0, 6.096]m . Periodic waves are generated using the internal wave generator of [START_REF] Wei | Generation of waves in Boussinesq models using a source function method[END_REF], centred at x = 4m. Sponge layers of 6m length are placed at both ends of the domain, while reflective boundary conditions are imposed along the top and bottom boundaries. Following [START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF] and [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF] for case (b) a relatively sparse triangular grid was used, consisting of equilateral triangles with side length of 0.1m, leading to a mesh of N=28,151 while a more refined mesh of side length 0.05m and N=56,211 nodes had to be used for case (a). The CFL number used was set equal to 0.5 for all cases.

Figures 2 and5 give an illustration of the fully developed 3D free surface elevation for case (a) and (b) respectively, while figures 3 and 6 depicts the free surface elevations at the centreline at the final time of the simulation. The incoming waves are linear in the deeper portion of the tank, but as they propagate onto the topography they become steeper due to shoaling. After the focusing, wave energy gradually spreads out due to diffraction. Figure 4 compares the spatial evolution of the first and second harmonics with the experimental data. In order to make sure that a steady periodic state is obtained, we started sampling the solution after at least 15 periods of oscillations. Then we performed a Discret Fourier Transform on the time series of the free surface elevation along the centreline, measured over one period of the main incoming wave. It can be observed that both the first and second harmonics increase in magnitude in the focal zone and the numerical harmonics are consistent with the laboratory data but slightly underestimate it. Case (b) has a weaker dispersive degree but a higher relative nonlinearity. The incoming waves are linear, but after the focusing on the shoal, higher harmonics become significant due to non-linear effects and the energy transfer in to the second and third harmonics can de seen in figure 7.The results are consistent and in very good agreement with the experimental data.

Solitary wave run up on a plane beach

The next test case is one of the most intensively studied problems in long-wave modelling. The solitary wave-run-up on a plane beach. [START_REF] Synolakis | The run-up of solitary waves[END_REF] carried out laboratory experiments for incident solitary waves of multiple relative amplitudes over a planar beach with a slope 1 : 19.85. Performing this test, we want to asses the ability of our model to describe propagation, shoreline motions, breaking and run-up. Detailed description of the test case, along with the initial conditions, can be found in [START_REF] Synolakis | The run-up of solitary waves[END_REF]), [START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2DH improved Boussinesq equations[END_REF] and [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF] among many others. The incident wave height used in this work is A = 0.28m and h 0 = 1m. This wave breaks strongly both in the run-up and run-down phases of the motion.

The computational domain used is [0, 70] × [0, 0.08]m. Wall boundary conditions were used at the top and bottom of the domain and a sponge layer is applied offshore with length Ls = 5m The CFL number was set equal to 0.2. The hybrid wave breaking model was used with γ = 0.6 Inria and tan(φ c ) = 0.37. Finally, a Manning coefficient of N m = 0.01 was used to define the glass surface roughness used in the experiments. Figures 8 and9 compares the measured surface profiles and the numerical model's results on different non-dimensional times.

The solitary wave, util time t g/h = 15 , propagates to the shore and shoaling over the inclined bathymetry. The breaking procedure starts around t g/h = 15 and the experimental wave, breaks around t g/h = 20. The numerical solution is represented like a bore storing the water spilled from the breaking wave behind the front. The numerical model reduce to NSWE during the breaking event, which approximate the turbulent breaking wave by a flow discontinuity hence, the computed front face becomes steeper. The actual wave is not discontinuous but contains air bubbles and turbulence. The bore collapse on the slope at time t g/h = 25. Because of volume conservation in all models, the computed solution fully recovers until the water reaches the maximum run-up point around t g/h = 45. A good agreement with the experimental data is observed for times t g/h = 60 and t g/h = 60.

Solitary wave propagation over a three dimensional reef

Swigler et Lynett (2011) performed laboratory experiments at the O.H. Hinsdale Wave Research Laboratory of Oregon State University to study the specific phenomena which occurs when a tsunami like wave approaches the coast, including shoaling, refraction, breaking and run-up. The basinis 48.8m long 26.5m wide and 2.1m deep. An extensive description of the set-up of the problem can be found in [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq type equations[END_REF] and in [START_REF] Filippini | Free surface flow simulation in estuarine and coastal environments: numerical development and application on unstructured mehes[END_REF]. Nine wave gauges were placed into the basin in order to measure the variation of the free surface elevation: gauges 1, 2, 3, 7 were located at y = 0m and x = 7.5, 13, 21, 25m gauges 4, 5, 6, 8 were located at y = 5m and x = 7.5, 13, 21, 25m, while gauge 10 have been set at y = 10m and x = 25m.

Compared to the experimental case, the computational domain has been extended from x = 0m to x = 5m in order to be able to completely contain the initial solitary wave. It has been discretized by means of a non-uniform unstructured grid, adapted to the bed curvature, as shown in 10, and characterized by reference maximum and minimum size respectively: max(h K ) = 0.3m and min(h K ) = 0.125m. A solitary wave of amplitude A = 0.39m, corresponding to = 0.5, is initially placed in x = 0 and wall reflecting boundary conditions are imposed in each boundary of the domain. We used a Manning coefficient N m = 0.0014 for representing bed roughness. A CFL number of 0.5 was used, together with γ = 0.6 for the breaking detection criterion.

Figure 11 shows the computed free water surface at different time instants. With red colour we denote the time evolution of the breaking regions detected by the criteria of the breaking mechanism. As the solitary wave propagates towards the beach it shoal, increases its steepness and nonlinearity, up to reaching a breaking point at t = 5 on the center line of the domain, when it reaches the apex of the triangular shelf. At t = 6.5s the central part of the wave has completely overtopped the concrete cone, while on the two sides, the surge continues to shoal, diffracting around the base of the cone. By t = 8.5s, the refracted and diffracted waves collide on the lee side of the shelf. After t = 9s, the water starts to withdraw from the con top and a bore-front forms, from the combined waves after the diffraction, and propagates on the shelf behind the cone and then onshore. After t = 15s, a new bore is creates from the the drawn-down of the water and collides with the refracted waves.

The next figure 12 plots the computed free surface time series on the wave gauges 2, 3, 7 and 9 against the experimental data. The arrival of the first incoming wave is correctly captured in gauge 2, as it is for the refracted and diffracted waves at the lee side of the cone, as can be seen from gauge 3, except for the minimum of water height registered at t = 8.5s. A slight deviation from the measurements is displayed after t = 35s, maybe due to late arrival of the reflected waves from the extended left wall boundary. The signal at the gauges located at the north side of the cone indicates that wave shoaling, breaking and propagation on the shelf is accurately predicted, together with the complex nonlinear interaction between diffracted and refracted waves.

Conclusions

In this work a fully nonlinear-weakly dispersive unstructured grid wave model is presented. The model is a first unstructured grid generalization of the fully nonlinear weakly dispersive models based on the seminal work of Wei and Kirby (1995), and Wei et al. (1995). The underlying model of partial differential equations if the enhanced fully nonlinear -weakly dispersive Green-Naghdi system proposed bu [START_REF] Bonneton | Recent advances in Serre-Green Naghdi modelling for wave transformation, breaking and runup processes[END_REF]. Following [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF], the system is decomposed in the hyperbolic shallow water equations plus an algebraic correction satisfying an auxiliary system of partial differential equations. These equations are elliptic and defined by a self-adjoint operator, as observed in Alvarez-Samaniego and Lannes (2008). We exploit this property to construct a hybrid finite element-finite volume scheme. In particular, the self adjoint nature of the elliptic operator is used to construct a coercive continuous finite element approximation providing the algebraic source term to add to a classical hyperbolic finite volume solver. Following [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF], and in the spirit of Wei and Kirby (1995), we have used a second order approximation of this elliptic high order correction. For the hyperbolic solver, we have instead used a non-dispersive, higher order MUSCL approach. The dispersion properties of the resulting scheme are shown in [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF] to be comparable to those of a fully fourth order finite difference approximation. Wave breaking is embedded in the model simply locally neglecting the dispersive correction, and following the detection criteria of [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq type equations[END_REF]. The results show an accurate capturing of highly dispersive wave transformation, as well as an accurate resolution of complex interactions involving wave breaking and runup. This work opens the door to a large number of model improvements such as e.g.

• the study of optimized strategy to solve the assemble and solve the block symmetric linear systems obtained in the elliptic phase;

• the use of simplified variants of the Green-Naghdi equations, such as those proposed by [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF] which would allow to obtain constat in time matrices for the dispersive terms;

• the study of the benefits of higher orders for the elliptic phase in terms of error reduction at constant CPU time;

• the use of other schemes in the hyperbolic phase, such as e.g. the successful Discontinuous Galerkin method;

• the combination of our model with static and dynamic mesh adaptation techniques tailored to these flows;

• the extension to spherical coordinate systems;

Initial application and benchmarking of this model are being carried out in the framework of the French research program TANDEM (http://www-tandem.cea.fr), and will be pushed further in the near future.

Inria

•

  the surface variation criterion: |η t | ≥ γ √ gh with γ ∈ [0.35, 0.65] • the local slope angle criterion: ||∇η|| ≥ tan φ c with φ c the critical angle value.

Figure 1 :

 1 Figure 1: Grid convergence for the propagation of a 2D solitary wave

Figure 2 :Figure 3 :

 23 Figure 2: Shoaling of regular waves, case (a): perspective view of the free surface

Figure 4 :

 4 Figure 4: Shoaling of regular waves, case(a): comparison of the computed and experimental results for the wave amplitudes for the first and second harmonics along the centreline

Figure 7 :

 7 Figure 7: Shoaling of regular waves, case(b): comparison of the computed and experimental results for the wave amplitudes for the first and second harmonics along the centreline Inria

Figure 8 :

 8 Figure 8: Solitary wave run-up on a plane beach: snapshots of the free surface elevation for dimensionless times t g/h = 15, 20, 25, 45 (from left to right and top to bottom ).

Figure 9 :

 9 Figure 9: Solitary wave run-up on a plane beach: snapshots of the free surface elevation for dimensionless times t g/h = 60, 80.

Figure 10 :

 10 Figure 10: Solitary wave propagation over a three dimensional reef: Close up view of the adapted mesh used for the computation(left) and computed free surface solution with friction at times t = 3.5s.

Figure 11 :

 11 Figure 11: Solitary wave propagation over a three dimensional reef: computed free surface solution with friction at times t = 3.5, 5.5, 6.5, 8.5, 9.5, 11.5s (from top to bottom). The red area represents the region where wave breaking is detected and the NLSW equations are solved.

Figure 12 :

 12 Figure 12: Solitary wave propagation over a three dimensional reef: computed time series of the free surface elevation on gauges positions

RR n°9013

avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

Work partially funded by the TANDEM contract, reference ANR-11-RSNR-0023-01 of the French Programme Investissements d' Avenir.