DMA D Lannes
email: david.lannes@ens.fr

A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations

We introduce a new class of two-dimensional fully nonlinear and weakly dispersive Green-Naghdi equations over varying topography. These new Green-Naghdi systems share the same order of precision as the standard one but have a mathematical structure which makes them much more suitable for the numerical resolution, in particular in the demanding case of two dimensional surfaces. For these new models, we develop a high order, well balanced, and robust numerical code relying on an hybrid finite volume and finite difference splitting approach. The hyperbolic part of the equations is handled with a high-order finite volume scheme allowing for breaking waves and dry areas. The dispersive part is treated with a finite difference approach. Higher order accuracy in space and time is achieved through WENO reconstruction methods and through a SSP-RK time stepping. Particular effort is made to ensure positivity of the water depth. Numerical validations are then performed, involving one and two dimensional cases and showing the ability of the resulting numerical model to handle waves propagation and transformation, wetting and drying; some simple treatments of wave breaking are also included. The resulting numerical code is particularly efficient from a computational point of view and very robust; it can therefore be used to handle complex two dimensional configurations.

Introduction

Propagation of surface water waves in an incompressible, homogeneous, inviscid fluid, is governed by the so called free surface Euler (or water waves) equations. Despite recent advances on the analysis of these equations (see for instance the recent monograph [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF] and references therein), this problem remains mathematically and numerically challenging. In particular several configurations of interest for applications such as wetting and drying and wave breaking still remain out of reach.

For coastal applications very important unsteady processes, such as wave run-up in the swash zone, coastal flooding due to storm waves, tsunami and tidal bore propagation, must be handled by phase-resolving models (as opposed to the so called phase averaged approach [START_REF] Phillips | The Dynamics of the Upper Ocean[END_REF]). For coastal applications, these models are based on the nonlinear shallow water (NSW) equations that give a good description of the nonlinear transformation of waves, including wavebreaking with broken waves represented as shocks. The NSW (or Saint-Venant) equations can be derived from the full water waves equations by neglecting all the terms of order O(µ), where the shallowness parameter µ is defined as

µ := h 2 0 L 2 ,
with h 0 the typical depth and L the typical horizontal scale. However, the NSW equations are not fully satisfactory since they neglect all the dispersive effects that play a very important role in many situations, and in particular in wave shoaling. In order to take dispersive effects into account, one must keep the O(µ) terms neglected by the NSW equations and neglect only the O(µ 2 ) terms; this leads to a more accurate -but mathematically and numerically more involved -set of equations known under several names: Serre [START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF][START_REF] Su | Korteweg-de Vries equation and generalizations. III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF][START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF], or Green-Naghdi [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF], or fully nonlinear Boussinesq [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF][START_REF] Miles | Weakly dispersive nonlinear gravity waves[END_REF] equations. We call them here Green-Naghdi (GN) equations. Simpler models such as the weakly nonlinear Boussinesq equations can be further derived from these equations upon making additional smallness assumptions on the amplitude of the waves and/or the topography variations. On these aspects, we refer the reader to several reviews [START_REF] Barthelemy | Nonlinear shallow water theories for coastal waves[END_REF][START_REF] Brocchini | Nonlinear shallow water equation modeling for coastal engineering[END_REF][START_REF] Lannes | Derivation of asymptotic two-dimensional timedependent equations for surface water wave propagation[END_REF][START_REF] Bonneton | Recent advances in Serre-Green naghdi modelling for wave transformation, breaking and runup processes[END_REF][START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF]. Among recent numerical studies of the Green-Naghdi equations, let us mention [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF][START_REF] Li | High order well-balanced CDG-FE methods for shallow water waves by a Green-Naghdi model[END_REF]; note also that the Green-Naghdi equations have recently be generalized in order to allow general nonzero vorticity [START_REF] Castro | Fully nonlinear long-waves models in presence of vorticity[END_REF].

Unfortunately, and despite the fact that they include the dispersive effects neglected by the NSW equations, the GN equations are not satisfactory neither. If they provide a very good description of the waves before wave breaking (see [START_REF] Alvarez-Samaniego | A Nash-Moser theorem for singular evolution equations. Application to the Serre and Green-Naghdi equations[END_REF] and also [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF] for a mathematical justification of this fact), they cannot be used to handle broken waves. Several approaches have been proposed to get rid of this important limitations. The most common one is to add an artificial viscous term to the momentum equation, whose role is to account for the energy dissipation that occurs during wave breaking (see for instance [START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking, and runup. I: 1D[END_REF][START_REF] Chen | Boussinesq modeling of wave transformation, breaking, and runup. II: 2d[END_REF], and [START_REF] Cienfuegos | A wave-breaking model for Boussinesq-type equations including mass-induced effects[END_REF] for a related approach).

Another approach to handle wave breaking is based on the fact mentioned above that the NSW equations, as a particular example of hyperbolic system, exhibit shocks that turn out to provide a good description of broken waves [START_REF] Bonneton | Modelling of periodic wave transformation in the inner surf zone[END_REF][START_REF] Brocchini | Nonlinear shallow water equation modeling for coastal engineering[END_REF], which can therefore be computed using finite volume schemes (we also refer to recent improvements in the "hyperbolic treatment" of hydraulic jumps [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF]). In [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF] for instance, the GN equations are written in a (pseudo)conservative form; since the system is no longer hyperbolic due to the dispersive terms, the Riemann problem associated to this system can only be approximately solved. A similar approach had also been used for non breaking waves, possibly including the dispersive terms as source terms for the NSW equations [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF][START_REF] Erduran | Hybrid finite-volume finite-difference scheme for the solution of Boussinesq equations[END_REF][START_REF] Erduran | Further application of hybrid solution to another form of Boussinesq equations and comparisons[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] Shiach | A temporally second-order accurate Godunovtype scheme for solving the extended Boussinesq equations[END_REF][START_REF] Soares-Frazao | A second-order semi-implicit hybrid scheme for one-dimensional Boussinesq-type waves in rectangular channels[END_REF]. The ability of the NSW equations to handle wave breaking can also be exploited in the numerical simulation of the GN equations by switching the dispersive terms off in the vicinity of wave breaking. This approach has been used for simpler weakly nonlinear Boussinesq models, like the Madsen and Sorensen equations [START_REF] Orszaghova | From the paddle to the beach : a Boussinesq shallow water numerical wave tank based on Madsen and S Ãÿrensen's equations[END_REF][START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Tonelli | Finite volume scheme for the solution of 2D extended Boussinesq equations in the surf zone[END_REF] or the Nwogu equations [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesqtype equations[END_REF][START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF]. It has also been recently implemented in the FunWave model [START_REF] Wei | A time-dependent numerical code for extended Boussinesq equations[END_REF] based on some fully nonlinear Boussinesq equations [START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF]; note also that the shock capturing approach has also been used recently for the non-hydrostatic models of [START_REF] Smit | Depth-induced wave breaking in a nonhydrostatic, near-shore wave model[END_REF]. Note that this latter reference, as well as [START_REF] Ma | Shock-capturing non-hydrostatic model for fully dispersive surface wave processes[END_REF] for instance, are based on the so the called non-hydrostatic shallow water equations that are promising and whose detailed analysis and comparison with Boussinesq type models remains to be done. This "switching" strategy was also adopted in [START_REF] Tissier | Nearshore Dynamics of Tsunami-like Undular Bores using a Fully Nonlinear Boussinesq Model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF], albeit in a slightly different form: the GN equations are written in a form particularly suited to the implementation of a splitting scheme between its hyperbolic (NSW) and dispersive components. More precisely, a second order splitting scheme is used. We compute the approximation U n+1 = (ζ n+1 , V n+1 ) at time (n + 1)δ t in terms of the approximation U n at time nδ t by solving

U n+1 = S 1 (δ t /2)S 2 (δ t )S 1 (δ t /2)U n ,
where S 1 (•) is the solution operator associated to the NSW equations, and S 2 (•) the solution operator associated to the dispersive part of the GN equations. For the numerical computation of S 1 (•), we use a high order, robust and wellbalanced finite volume method, based on a relaxation approach [START_REF] Berthon | A Positive Preserving High Order VFRoe Scheme for Shallow Water Equations: A Class of Relaxation Schemes[END_REF]. In the vicinity of shocks detected during the computation of S 1 (δ t /2), we "skip" the computation of S 2 (•) [START_REF] Tissier | Nearshore Dynamics of Tsunami-like Undular Bores using a Fully Nonlinear Boussinesq Model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] . Elsewhere, S 2 (•) is computed using a finite difference scheme. Note that another motivation for this scheme was to get a numerical benefit of the mathematical stability properties that allow one to establish the well posedness of the GN equations [START_REF] Alvarez-Samaniego | A Nash-Moser theorem for singular evolution equations. Application to the Serre and Green-Naghdi equations[END_REF].

The treatment of vanishing depth (shoreline) for nonlinear dispersive systems such as the GN or weakly nonlinear Boussinesq equations is another mathematical and numerical difficulty. For the (hyperbolic) NSW equations, several techniques have been developed to handle such configurations, such as [START_REF] Brocchini | A comparison of two different types of shoreline boundary conditions[END_REF][START_REF] Van't Hof | Modelling of wetting and drying of shallow water using artificial porosity[END_REF][START_REF] Lynett | Modeling wave runup with depth-integrated equations[END_REF][START_REF] Marche | Evaluation of well-balanced bore-capturing schemes for 2D wetting and drying processes[END_REF] for instance. Consequently, several of the numerical schemes developed for the GN or Boussinesq equations and based on an extension of a finite volume code for the NSW equations are able to handle the shoreline [START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF][START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF][START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF]; this is in particular the case with the splitting scheme described above [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Tissier | Nearshore Dynamics of Tsunami-like Undular Bores using a Fully Nonlinear Boussinesq Model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF]. However, the issue of robustness for higher order schemes is usually not addressed.

Additionally, one of the major challenging issues is that solving the GN equations is particularly computationally demanding, especially for practical large scale applications. The splitting approach developed in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF] has been validated with many experimental data [START_REF] Tissier | Nearshore Dynamics of Tsunami-like Undular Bores using a Fully Nonlinear Boussinesq Model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] for one dimensional surface waves. While the extension of this approach to two-dimensional surface waves does not raise theoretical difficulty, it should require a lot of computational time in complex configurations. Indeed, the computation of S 2 (δ t ) in the above splitting scheme (dispersive part) requires the inversion of a second order differential operator (I + µαT) (see [START_REF] Berthon | A Positive Preserving High Order VFRoe Scheme for Shallow Water Equations: A Class of Relaxation Schemes[END_REF] below); this is computationally demanding in the 2D case for several reasons: ♯ 1 It is a matricial operator of size 2 × 2 which is moreover not in diagonal form; once discretized on a cartesian mesh, this leads to the resolution of an 2(N x N y ) × 2(N x N y ) linear system with many nonzero diagonals, with N x and N y the number of cells respectively in the first and second direction.

♯ 2 It is a time dependent operator depending on the the water height h and the corresponding matrix has therefore to be modified at each time step. These computational difficulties are actually one of the reasons why most of the 2D numerical codes are based on weakly nonlinear Boussinesq models (see for instance the recent studies [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF][START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesqtype equations[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF]) rather than the fully nonlinear Green-Naghdi equations. Of course, one could treat the dispersive components as source terms and therefore avoid the inversion of the operator (I + µαT), but one would loose the high robustness of the code. In this work, in order to keep the stabilizing effects of the inversion step without losing too much in terms of computational time, we derive a new family of physical models that are equivalent to the standard GN equations (in the sense that they differ only by O(µ 2 ) terms). It is possible to implement our splitting strategy on these new systems with the benefit of removing the two numerical obstructions mentioned above. This is done in two steps:

♯ 1 For the system derived in §2.2, the matricial operator T is replaced by a new differential operator T diag , which is still matricial but has a diagonal structure. Its inversion is therefore equivalent to the inversion of two scalar operators; numerically, this leads to the resolution of two (N x N y )×(N x N y ) linear systems with only a few nonzero diagonals.

♯ 2 For the system derived in §2.3, the matricial operator T is replaced by an operator T b diag which, as T diag , has a diagonal structure, but for which the dependence on h has been removed; it is therefore time independent.

Additionally, like the GN formulation introduced in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], these new models do not involve any third order derivatives on the flow variables.

The simple structure of the systems thus obtained leads to considerable improvements in terms of computational time, without any loss compared to the standard GN system in terms of accuracy; we can therefore consider complex 2D configurations including wetting and drying processes. This computational efficiency also makes possible the construction of dispersion enhanced models for configurations that involve deep water waves for instance.

In Section 2, the derivation of these new physical models is presented, and their stability and dispersive properties are discussed. Section 3 is then devoted to the presentation of the numerical scheme. The hyperbolic/dispersive splitting is recalled in §3.1 before we describe the vertex-centered finite-volume spacial discretization of the hyperbolic part and the finite-difference discretization of the dispersive parts in §3.2 and §3.3 respectively. In particular, we recall and apply the strategy of [START_REF] Xing | Positivity-preserving high order wellbalanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF] to combine high-order accuracy and positivity preserving properties for the water depth, leading to a robust fourth order model. The boundary conditions and the time discretization are described in §3.4 and §?? respectively We then explain in §3.5, how we handle wetting and drying processes, while a brief description of the way we treat wave breaking is provided in §3. [START_REF] Berkhoff | Verification of numerical wave propagation models for simple harmonic linear water waves[END_REF]. Finally, we present in Section 4 several numerical validations of our model. We first consider the case of solitary waves in §4.1 and use it as a validation tool for our numerical scheme, studying accuracy and convergence. Then, we study the transformation of a solitary wave propagating over a composite beach, the propagation and overtopping of a solitary wave over a sloping seawall, and the propagation of periodic waves over a submerged bar, and compare our results with experimental data. We then turn on considering two-dimensional cases with the shoaling of a regular wave trains over a semi-circular shoal and an elliptic shoal. Lastly, we highlight the robustness of our approach through the computation of run-up of solitary waves over varying topography in §4.8 .

The physical model

Throughout this paper, we denote by ζ(t, X) the elevation of the surface with respect to its rest state, and by -h 0 + b(X) a parametrization of the bottom, where h 0 is a reference depth (see Figure 1). Here X stands for the horizontal variables X = (x, y), and t is the time variable; we also denote by z the vertical variable. If U hor denotes the horizontal component of the velocity field in the fluid domain, we then define the averaged horizontal velocity V as Denoting by a the typical amplitude of the waves, by a bott the typical amplitude of the bottom variations, and by L the order of the wavelength of the wave, it is possible to define dimensionless variables and unknowns as

V (t, X) = 1 h ζ -h0+b U hor (t, X, z)dz, z 0 72-h 0 z = ζ(t, X) z = -h 0 + b(X) X = (x, y)
X = X L , t = √ gh 0 L t and ζ = ζ a , b = b a bott , V = V √ gh 0 .
We also define three dimensionless parameters as

ε = a h 0 , µ = h 2 0 L 2 , β = a bott h 0 ;
here ε denotes the nonlinearity parameter, µ is the shallowness parameter while β accounts for the topography variations.

The non-dimensionalized Green-Naghdi equations

According to [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF], the classical Green-Naghdi equations can be written under the following non-dimensionalized form (we omit the tildes for dimensionless quantities for the sake of clarity):

∂ t ζ + ∇ • (hV ) = 0, (I + µT [h, b]) ∂ t V + ε(V • ∇)V + ∇ζ + εµQ 1 (V ) = 0 (1)
(the momentum equation being actually satisfied up to O(µ 2 ) terms), where we still denote by h the non-dimensionalized water depth,

h = 1 + εζ -βb,
and the linear operator

T [h, b]• and the quadratic form Q 1 [h, b](•) are defined for all smooth enough R 2 -valued function W by T [h, b]W = R 1 [h, b](∇ • W ) + βR 2 [h, b](∇b • W ) (2) Q 1 [h, b](W ) =-2R 1 (∂ 1 W • ∂ 2 W ⊥ + (∇ • W ) 2 ) + βR 2 (W • (W • ∇)∇b) (3)
(here ∂ 1 and ∂ 2 denote space derivatives in the two horizontal directions) with, for all smooth enough scalar-valued function w,

R 1 [h, b]w = - 1 3h ∇(h 3 w) -β h 2 w∇b, (4)
R 2 [h, b]w = 1 2h ∇(h 2 w) + βw∇b. (5)
Notation 2.1. For the sake of clarity, and when no confusion is possible, we often write

T , Q, R 1 and R 2 instead of T [h, b], Q[h, b], etc.
It is also shown in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] following classical ideas of [START_REF] Witting | A unified model for the evolution of nonlinear water waves[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF] or [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary conditions and validations[END_REF] that the frequency dispersion of (1) can be improved by adding some terms of order O(µ 2 ) to the momentum equation. Since this equation is already precise up to terms of order O(µ 2 ), this manipulation does not affect the accuracy of the model. In this perspective, (1) is a particular case (corresponding to α = 1) of the following family of asymptotically equivalent models parametrized by α > 0,

∂ t ζ + ∇ • (hV ) = 0, (I + µαT ) ∂ t V + ε(V • ∇)V + α -1 α ∇ζ + 1 α ∇ζ + εµQ 1 (V ) = 0. (6)
It is of course possible to rewrite these equations in (h, hV ) variables, which yields, following [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], the following equivalent formulation of (6),

∂ t h + ε∇ • (hV ) = 0, (I +µαT) ∂ t (hV )+ε∇•(hV ⊗V )+ α-1 α h∇ζ + 1 α h∇ζ +εµhQ 1 (V ) = 0, (7)
with

TW = hT ( 1 h W ).
This formulation has two main advantages:

1. It does not require the computation of third order derivatives 1 , while this is necessary in the standard formulation of the GN equations.

2. The presence of the operator (I + µαT) -1 in the second equation makes the model very stable with respect to high frequency perturbations, which is of course of highest interest for numerical computations. Note that this property is also used to prove that these equations are (nonlinearly) well posed (see [START_REF] Li | A shallow-water approximation to the full water wave problem[END_REF][START_REF] Israwi | Large Time existence For 1D Green-Naghdi equations[END_REF] for the 1D case and [START_REF] Alvarez-Samaniego | A Nash-Moser theorem for singular evolution equations. Application to the Serre and Green-Naghdi equations[END_REF] for 2D surface waves).

1 Of course, third order derivatives are qualitatively present in the model, but the fact that they can be factored out by (I +µαT) makes it possible not to compute them. Indeed, this equation can be written under the equivalent form

∂t(hV ) + ε∇ • (hV ⊗ V ) + α -1 α h∇ζ + (I +µαT) -1 1 α h∇ζ + εµhQ 1 (V ) = 0.
The main drawback of this method is that inverting2 the operator (I + µαT) is computationally demanding in the 2D case. Indeed, it is a matricial differential operator of order 2 acting on two-dimensional vectors; the operator T can be written in matricial form as

T = T 11 T 12 T 21 T 22 ,
where T ij (1 ≤ i, j ≤ 2) are second order scalar differential operators. The fact that the antidiagonal terms are nonzero complicates considerably the numerical computation of (1 + µαT) -1 . This is the reason why we derive in the following section a new family of models for which T is replaced by an operator T diag of the form

T diag = T 11 0 0 T 22 , (8)
where T 11 and T 22 are second order scalar differential operators. Since for T diag the antidiagonal coefficients are zero, the computation of (1 + µαT diag ) -1 reduces to inverting two (namely, 1 + µα T 11 and 1 + µα T 22 ) second order scalar differential operators. This is of course much simpler than inverting 1 + µαT.

The "diagonal" 2D Green-Naghdi equations

We derive here a new family of systems asymptotically equivalent to [START_REF] Berthon | A Positive Preserving High Order VFRoe Scheme for Shallow Water Equations: A Class of Relaxation Schemes[END_REF] but with a structure more adapted to the computation of 2D surface waves, as explained in the previous section. More precisely, we want to rewrite the second equation of (1), namely,

(I + µT ) ∂ t V + ε(V • ∇)V + ∇ζ + εµQ 1 (V ) = 0,
under the form

(I + µT diag ) ∂ t V + ε(V • ∇)V + ∇ζ + εµQ 1 (V ) + µQ 2 (ζ) = 0,
where T diag has a "diagonal" structure (as discussed in the previous section), while Q 2 (ζ) involves only second order space derivatives of ζ, and where O(µ 2 ) terms have been discarded. Denoting

W := [∂ t V + ε(V • ∇)V , (9)
we can invoke the definition (2) of T to write

T W = - 1 3h ∇(h 3 ∇ • W ) :=A1 + β 2h ∇(h 2 ∇b • W ) -h 2 ∇b∇ • W + β 2 ∇b∇b • W :=A2
we therefore turn to transform the two components A j (j = 1, 2) of this expression.

• Transformation of A 1 . We can remark that

A 1 = - 1 3h ∂ 1 (h 3 ∂ 1 W )- 1 3h ∂ 2 (h 3 ∂ 2 W )+ 1 3h ∇ ⊥ (h 3 ∇ ⊥ •W )+h(∇ ⊥ h•∇)W ⊥ ,
with ∇ ⊥ = (-∂ 2 , ∂ 1 ) T , with ∂ 1 and ∂ 2 the partial derivatives in the first and second spatial directions. Now, using the second equation of ( 1) and recalling the definition (9) of W , we have that W = -∇ζ + O(µ) and therefore

A 1 = - 1 3h ∂ 1 (h 3 ∂ 1 W ) - 1 3h ∂ 2 (h 3 ∂ 2 W ) -h(∇ ⊥ h • ∇)∇ ⊥ ζ + O(µ).
• Transformation of A 2 . Substituting as above W = -∇ζ + O(µ) in the expression for A 2 , we get directly

A 2 = - β 2h ∇(h 2 ∇b • ∇ζ) -h 2 ∇b∆ζ -β 2 ∇b∇b • ∇ζ + O(µ).
It follows from this analysis that

T W = - 1 3h ∂ 1 (h 3 ∂ 1 W ) - 1 3h ∂ 2 (h 3 ∂ 2 W ) -h(∇ ⊥ h • ∇)∇ ⊥ ζ - β 2h ∇(h 2 ∇b • ∇ζ) -h 2 ∇b∆ζ -β 2 ∇b∇b • ∇ζ + O(µ),
with W as in [START_REF] Berthon | Asymptotic preserving scheme for the shallow water equations with source terms on unstructured meshes[END_REF]. Replacing in (1) and neglecting O(µ 2 ) terms (which does not affect the precision of the model), we obtain

∂ t ζ + ∇ • (hV ) = 0, I + µT diag ∂ t V + ε(V • ∇)V + ∇ζ + εµQ 1 (V ) + µQ 2 (ζ) = 0, (10)
where

T diag = - 1 3h ∂ 1 (h 3 ∂ 1 •) - 1 3h ∂ 2 (h 3 ∂ 2 •),
and

Q 2 (ζ) = -h(∇ ⊥ h • ∇)∇ ⊥ ζ - β 2h ∇(h 2 ∇b • ∇ζ) + β h 2 ∆ζ -β(∇b • ∇ζ) ∇b. (11)
The system (10) forms the "diagonal" 2D Green-Naghdi equations. For our numerical simulations, we rewrite it in terms of the (h, hV ) variables. Proceed as in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], and using the identities

∂ t h = ε∂ t ζ, h∂ t V = ∂ t (hV ) + ε∇ • (hV )V, ∇ • (hV ⊗ V ) = ∇ • (hV )V + h(V • ∇)V, it is easy to deduce from (10) that    ∂ t h + ε∇ • (hV ) = 0, I + µT diag ∂ t (hV ) + ε∇ • (hV ⊗ V ) +h∇ζ + εµhQ 1 (V ) + µhQ 2 (ζ) = 0, (12)
with

T diag W = hT diag 1 h W ,
and where we recall that h = 1 + εζ -βb, while Q 1 and Q 2 are defined in ( 3) and ( 11) respectively.

Remark 2.2. The interest of working with (12) rather than [START_REF] Berthon | A Positive Preserving High Order VFRoe Scheme for Shallow Water Equations: A Class of Relaxation Schemes[END_REF] is that instead of inverting the operator (I + µT), we have to invert (I + µT diag ) that has the diagonal structure (8), with

T 11 f = T 22 f = - 1 3 ∂ 1 (h 3 ∂ 1 f h ) - 1 3 ∂ 2 (h 3 ∂ 2 f h ).
The other difference with [START_REF] Berthon | A Positive Preserving High Order VFRoe Scheme for Shallow Water Equations: A Class of Relaxation Schemes[END_REF] is the presence of the term µhQ 2 (ζ) in the last term, which shares with [START_REF] Berthon | A Positive Preserving High Order VFRoe Scheme for Shallow Water Equations: A Class of Relaxation Schemes[END_REF] the useful properties that no computation of third order derivative is needed (see Footnote 1).

2.3

The "constant-diagonal" 2D Green-Naghdi equations

The system [START_REF] Bonneton | Modelling of periodic wave transformation in the inner surf zone[END_REF] is much easier to solve numerically in two surface dimensions than the standard version [START_REF] Berthon | A Positive Preserving High Order VFRoe Scheme for Shallow Water Equations: A Class of Relaxation Schemes[END_REF] due to the diagonal structure of the operator I + µT diag that must be inverted. However, this system is time dependent and a great gain of computational time could be achieved if we were able to replace this diagonal operator by a diagonal operator with time independent coefficients. We show here that it is indeed possible to replace the inversion of I + µT diag by the inversion of I + µT b diag , where T b diag corresponds to T diag when the fluid is at rest (i.e. (ζ = 0),

T b diag = h b T b diag 1 h b , with T b diag = - 1 3h b ∂ 1 (h 3 b ∂ 1 •) - 1 3h b ∂ 2 (h 3 b ∂ 2 •),
and where h b is the water depth at rest,

h b = 1 -βb = h -εζ.
Remarking that for all scalar function f , one has

hT diag 1 h f = h b T b diag 1 h b f + S[h 2 -h 2 b ]f with S[h 2 -h 2 b ]f = - 1 6 ∇(h 2 -h 2 b ) • ∇f - h 2 -h 2 b 3 ∆f + 1 6 ∆(h 2 -h 2 b )f, (13)
we can write the second equation of ( 12) under the form

I + µT b diag ∂ t (hV ) + ε∇ • (hV ⊗ V ) + h∇ζ + εµhQ 1 (V ) + µhQ 2 (ζ) +µS[h 2 -h 2 b ] ∂ t (hV ) + ε∇ • (hV ⊗ V ) = 0.
We infer in particular from this equation that

∂ t (hV ) + ε∇ • (hV ⊗ V ) = -I + µT b diag -1 (h∇ζ) + O(µ), (14)
and therefore

I + µT b diag ∂ t (hV ) + ε∇ • (hV ⊗ V ) + h∇ζ + εµhQ 1 (V ) + µhQ 2 (ζ) -µS[h 2 -h 2 b ] I + µT b diag -1 (h∇ζ) = O(µ 2 ).
Dropping O(µ 2 ) terms, the "diagonal" Green-Naghdi equations ( 12) are equivalent to the following "constant-diagonal" system,   

∂ t h + ε∇ • (hV ) = 0, I + µT b diag ∂ t (hV ) + ε∇ • (hV ⊗ V ) + h∇ζ +εµhQ 1 (V ) + µhQ 2 (ζ) + µQ 3 (ζ) = 0. with, recalling that S[h 2 -h 2 b ] is defined in (13), Q 3 (ζ) = -S[h 2 -h 2 b ] I + µT b diag -1 (h∇ζ). (15)
More generally, we can introduce a parameter α > 0 to find as in [START_REF] Berthon | A Positive Preserving High Order VFRoe Scheme for Shallow Water Equations: A Class of Relaxation Schemes[END_REF] the following family of asymptotically equivalent systems,

       ∂ t h + ε∇ • (hV ) = 0, I + µαT b diag ∂ t (hV ) + ε∇ • (hV ⊗ V ) + α -1 α h∇ζ + 1 α h∇ζ +εµhQ 1 (V ) + µhQ 2 (ζ) + µQ 3 (ζ) = 0 (16)
(note that the Q 3 term involves the inversion of a linear system but, as shown below, this extra cost is off-set by the gain obtained by using a time independent operator). Moreover, it allows to avoid the computation of third order derivatives on the flow variables.

Remark 2.3. Using ( 14) one more time, and dropping as usual O(µ 2 ) terms, one can replace the second equation of ( 16) by (we take here α = 1 for the sake of simplicity)

∂ t (hV ) + ε∇ • (hV ⊗ V ) + h∇ζ + εµhQ 1 (V ) + µhQ 2 (ζ) + µQ 3 (ζ) + µQ 4 (ζ) = 0 with Q 4 (ζ) = -T b diag (h∇ζ)
. With this model, there is no operator to invert and the gain in computational time is important, but we lose the stabilizing effects associated to this inverse operator. This model is consequently too instable to be of interest.

Remark 2.4. Since I + µT b diag -1 (h∇ζ) = h∇ζ + O(µ) one could replace Q 3 (ζ) by Q 3 (ζ) = -S[h 2 -h 2 b ](h∇ζ)
in the second equation of ( 16), keeping the same O(µ 2 ) precision. This would avoid the resolution of an extra linear system but leads to instabilities. This is due to the fact that Q 3 (ζ) contains third order derivatives in ζ that can create high frequencies instabilities. This problem does not arise with Q 1 (V ) and Q 2 (ζ) which contain at most second order derivatives. We refer to §2.4 for more comments on this point.

Remark 2.5. Going back to variables with dimension, the system of equations ( 16) reads

       ∂ t h + ∇ • (hV ) = 0, I + αT b diag ∂ t (hV ) + ∇ • (hV ⊗ V ) + α -1 α gh∇ζ + 1 α gh∇ζ +h Q 1 (V ) + gQ 2 (ζ) + gQ 3 (ζ) = 0, ( 17
)
where the dimensional versions of the operators Q 1 , Q 2 and Q 3 correspond to (3), ( 11) and ( 15) with β = 1, and where h (resp. h b ) now stands for the water height with dimensions (resp. at rest),

h = h 0 + ζ -b, and
h b = h 0 -b. (18)
Remark 2.6. It is classical, for some practical applications, to add a quadratic friction term to the right-hand side of the momentum equation. A commonly used possibility is to use the Manning-Chezy formulation of the friction dissipation: we add -n 2 h 10/3 hV hV to the right hand side of [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF], where n is the Manning friction coefficient.

Choice of α and stability issues

This section is devoted to a qualitative discussion on the stability of our new model [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF]. Considering the case of flat bottoms and in dimension 1 for the sake of simplicity, we investigate the linear behavior of small perturbation ( ζ, V ) to a constant state solution (ζ, V ). The linear equations that govern these perturbations are

         ∂ t ḣ + V • ∇ ḣ + h∇ • V = 0, I - α 3 h 2 0 ∆ ∂ t V + V • ∇ V + α -1 α g∇ ζ + 1 α g∇ ζ +g h 2 -h 2 0 3 ∆ I - α 3 h 2 0 ∆ -1 ∇ ζ = 0, (19)
with h = h 0 + ζ. Looking for plane wave solution of the form (ζ 0 , V 0 )e i(k•X-ωt) to this system, one finds the following dispersion relation

1 ghk 2 (ω -V • k) 2 = 1 + α -1 3 h 2 0 k 2 - h 2 -h 2 0 3 k 2 1 + α 3 h 2 0 k 2 1 + α 3 h 2 0 k 2 , ( 20
)
with k = |k|. This dispersion relation differs from the one found when applying the same procedure to the standard GN equations with improved dispersion (7) α , namely,

1 ghk 2 (ω GN,α -V • k) 2 = 1 + α -1 3 h 2 k 2 1 + α 3 h 2 k 2 ; ( 21
)
in this case the perturbations are always stable (i.e. ω GN,α is always a real number) if α ≥ 1. The choice of α is then classically made to obtain a good matching with the dispersion relation of the full Euler equations around the rest state (ζ, V ) = (0, 0) (and therefore h = h 0 ). We refer to Section 2.6 of [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] for a more complete discussion on the choice of α.

For our new model ( 17), the dispersion relation ( 20) is the same as ( 21) derived for the standard GN equations ( 7) α , when we consider the linearization around the rest state (ζ, V ) = (0, 0). The discussion concerning the choice of α in order to improve the dispersive properties of the model therefore follows the usual procedure and, as in [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary conditions and validations[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] we take throughout this article α = 1.159. However, when the linearization is performed around a non trivial constant state (ζ, V ) = (0, 0), ( 20) and ( 21) are obviously different; there is indeed an additional term, namely, - [START_REF] Chen | Boussinesq modeling of wave transformation, breaking, and runup. II: 2d[END_REF]. When h > h 0 , this new term has a negative sign and plays a destabilizing role; it is however compensated by the other terms and ω, as provided by [START_REF] Chen | Boussinesq modeling of wave transformation, breaking, and runup. II: 2d[END_REF] is always real at high frequencies, thus ensuring the local in time mathematical well-posedness and numerical stability of the linearization [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF] of our new model [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF]. A quick functional study of the right-hand-side of [START_REF] Chen | Boussinesq modeling of wave transformation, breaking, and runup. II: 2d[END_REF] shows moreover that all modes are stable under the condition h 2 h 2 0 < 2α, which, with our choice of α, leads to h/h 0 < 1.52; in most of the situations considered for applications, the overall height h does not exceed by more than 52% the depth at rest, and this condition is satisfied. When this condition is not satisfied, there is a local amplification of some modes, but it is harmless. In the extreme situation h = 1.8h 0 for instance, the right-hand-side of (20) becomes negative for 5 h 0 k 2 10 but remains very small (the minimum is roughly -0.01). The resulting amplification has therefore a very small rate and occurs for a very small time period since after wave breaking the condition h/h 0 < 1.52 is rapidly restored.

h 2 -h 2 0 3 k 2 1+ α 3 h 2 0 k 2 in the numerator of the right-hand- side of
Remark 2.7. The above discussion allows us to make more precise the comments made in Remark 2.4.

Replacing Q 3 (ζ) by Q 3 (ζ) = -S[h 2 -h 2 b ](h∇ζ)
in the second equation of ( 16), would modify the dispersion relation [START_REF] Chen | Boussinesq modeling of wave transformation, breaking, and runup. II: 2d[END_REF] and one would obtain instead

(ω -V • k) 2 = ghk 2 1 + α -1 3 h 2 0 k 2 -k 2 h 2 -h 2 0 3 1 + α 3 h 2 0 k 2 ; quite obviously, if h 2 -h 2 0 h 2 0 > α -1
, the numerator will become negative for sufficiently large values of k 2 , and the root ω will become complex inducing a high frequency instability of the model (see Fig. 2). Choosing α large enough can partially controls this instability, but at the cost of damaging the dispersive properties of the model; this why we work is the model (17) which, as seen above, is always stable.

Further optimization of the model

We have obtained above a family of "constant-diagonal" Green-Naghdi systems depending on one paramater α. It is possible to improve the dispersive properties of these models by adding two additional parameters θ and γ using a change of variables for the velocity. In order to derive such a three parameters family of "constant diagonal" Green-Naghdi systems, let us first go back to [START_REF] Bingham | Velocity potential formulations of highly accurate Boussinesq-type models[END_REF],

∂ t ζ + ∇ • (hV ) = 0, I + µT diag ∂ t V + ε(V • ∇)V + ∇ζ + εµQ 1 (V ) + µQ 2 (ζ) = 0;
following an original idea of [START_REF] Nwogu | An alternative form of the Boussinesq equations for nearshore wave propagation[END_REF], generalized to the fully nonlinear case in [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF] (we actually use the slightly different approach of [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF] and [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF], §5.2), we introduce a new velocity V θ as follows

∀θ ≥ 0, V = (1 + µθT diag )V θ .
We therefore have

V • ∇V = V θ • ∇V θ + µθ (T diag V θ ) • ∇ V θ + (V θ • ∇)T diag V θ + O(µ 2 );
in order to relate this quantity to (1

+ µθT diag ) (V θ • ∇)V θ , let us write T diag = a∆ + (A • ∇) with a = - h 2 3 and A = -h∇h, so that T diag (V θ • ∇)V θ = (T diag V θ ) • ∇ V θ + (V θ • ∇)(T V θ ) + k=1,2 2a(∇V θ,k • ∇)∂ k V θ -V θ,k (∂ k A) • ∇ V θ .
We therefore have

V • ∇V = (1 + µθT b diag ) (V θ • ∇)V θ ) + µθB, with B = k=1,2 2 3 h 2 (∇V θ,k • ∇)∂ k V θ + V θ,k ∇(h∂ k h) • ∇ V θ .
Substituting this expression in the velocity equation, this yields

(I + µ(1 + θ)T diag ) ∂ t V θ + ε(V θ • ∇)V θ + ∇ζ + εµQ 1 (V ) + µθ(B + (∂ t T diag )V θ ) + µQ 2 (ζ) = O(µ 2 ),
with the notation

∂ t T diag = (∂ t a)∆+(∂ t A•∇).
Remarking also that the equation on ζ written in terms of V θ is given by

∂ t ζ + ∇ • (hV θ ) + µθ∇ • (hT diag V θ ) = 0,
we can write

h ∂ t V θ + ε(V θ • ∇)V θ = ∂ t (hV θ ) -(∂ t h)V θ + ε(hV θ • ∇)V θ = ∂ t (hV θ ) + ε∇ • (hV θ ⊗ V θ ) + µεθ∇ • (hT diag V θ )V θ .
The velocity equation then becomes

(I + µ(1 + θ)T diag ) ∂ t (hV θ ) + ε∇ • (hV θ ⊗ V θ ) + h∇ζ + εµhQ 1 (V θ ) + µhQ 2 (ζ) + µθh(B + (∂ t T diag )V θ ) + µεθ∇ • (hT V θ )V θ = O(µ 2 ).
Expressing T diag in terms of T b diag as in §2.3, we get

(I + µ(1 + θ)T b diag ) ∂ t (hV θ ) + ε∇ • (hV θ ⊗ V θ ) + ∇ζ + εµhQ 1 (V θ ) + µhQ 2 (ζ) + µ(1 + θ)Q 3 (ζ) + µθQ 4 (V θ ) = O(µ 2 ), with Q 4 (V θ ) given by Q 4 (V θ ) = ε∇ • (hT diag V θ )V θ + h(B + (∂ t T diag )V θ ).
Up to O(µ 2 ) terms in the resulting equations, we can replace ∂ t h by -ε∇ • (hV θ ) in the expression for ∂ t T diag , leading to

∂ t T diag = ε 2 3 h∇ • (hV θ )∆ + ε∇ h∇ • (hV θ ) • ∇
and therefore

Q 4 (V θ ) =ε∇ • (hT diag V θ )V θ + ε 2 3 h 2 ∇ • (hV θ )∆V θ + εh∇ h∇ • (hV θ ) • ∇V θ + k=1,2 2 3 h 3 (∇V θ,k • ∇)∂ k V θ + hV θ,k ∇(h∂ k h) • ∇ V θ (22)
As for [START_REF] Brocchini | Nonlinear shallow water equation modeling for coastal engineering[END_REF], we can introduce another parameter α > 0 and obtain

(I+µα(1 + θ)T b diag ) ∂ t (hV θ ) + ε∇ • (hV θ ⊗ V θ ) + α -1 α ∇ζ + 1 α ∇ζ + εµhQ 1 (V θ ) + µhQ 2 (ζ) + µ(1 + θ)Q 3 (ζ) + µθQ 4 (V θ ) = O(µ 2 ).
The third parameter γ ≥ 0 can be introduced by applying (1 + µγT b diag ) to the equation on the surface elevation,

(1 + µγT b diag ) ∂ t ζ + ∇ • (hV θ ) + µθ∇ • (hT diag V θ ) = O(µ 2 ).
Neglecting all the O(µ 2 ) terms, the three parameters family of constant diagonal Green-Naghdi equations is given by

     (I + µγT b diag ) ∂ t ζ + ∇ • (hV θ ) + µθ∇ • (hT diag V θ ) = 0 (I + µα(1 + θ)T b diag ) ∂ t (hV θ ) + ε∇•(hV θ ⊗ V θ ) + α -1 α ∇ζ + 1 α ∇ζ +εµhQ 1 (V θ ) + µhQ 2 (ζ) + µ(1 + θ)Q 3 (ζ) + µθQ 4 (V θ ) = 0, (23)
with Q j (j = 1, . . . , 4) are given by ( 3), ( 11), ( 15) and ( 22) respectively.

Remark 2.8. Taking θ = γ = 0, ( 23) coincides of course with the one parameter family (16) indexed by α. An optimization of the phase and group velocities, as in [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF], leads to the following coefficients: α = 1.024, θ = 0.188 and γ = 0.112. In some stiff configurations where higher harmonics are released, working with the three parameters optimized model [START_REF] Davis | An unsymmetric-pattern multifrontal method for sparse LU factorization[END_REF] leads to considerable improvements, as shown in the numerical computations of §4.5.

Remark 2.9. The three-parameters Green-Naghdi systems derived in [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF] could of course be generalized to the 2D case, but the computational cost would be prohibitive. It is the computational gain of our new constant diagonal GN systems that make such dispersive improvements possible in 2D.

Numerical methods

As the standard GN equations ( 1), our new model ( 17) is well-adapted to the implementation of a splitting scheme separating the hyperbolic and the dispersive parts of the equations. We present in §3.1 this splitting scheme inspired by [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF]; We explain in §3.2 and §3.3 how we treat respectively the hyperbolic and dispersive parts of the equations, insisting on the particularities of the twodimensional case which was not considered in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF]; note however that even for 1D configurations, the numerical schemes proposed here for the two parts of the splitting bring significant improvement compared to [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF]. Indeed, the use of the constant diagonal formulation leads to a considerable decrease of the computational time. This gain allows us to use a more accurate and robust (but more computationally expensive) high-order reconstruction, namely the 5 th order WENO reconstruction together with the water height positivity preserving strategy recently introduced in [START_REF] Xing | Positivity-preserving high order wellbalanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF]. The discretization of the 3 parameters constant diagonal models is investigated in §3.7.

The splitting scheme

We decompose the solution operator S(•) associated to ( 17) -with a possible friction term as commented in Remark 2.6-at each time step by the second order Strang splitting scheme

S(δ t ) = S 1 (δ t /2)S 2 (δ t )S 1 (δ t /2), (24)
where S 1 and S 2 are respectively associated to the hyperbolic and dispersive parts of the Green-Naghdi equations [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF]. We have:

• S 1 (t) is the solution operator associated to NSWE    ∂ t h + ∇ • (hV ) = 0, ∂ t (hV ) + ∇( 1 2 gh 2 ) + ∇ • (hV ⊗ V ) = -gh∇b -n 2 h 10/3 hV hV. (25)
• S 2 (t) is the solution operator associated to the remaining (dispersive) part of the equations,

       ∂ t h = 0, ∂ t (hV ) - 1 α gh∇ζ + I + αT b diag -1 [ 1 α gh∇ζ + h Q 1 (V )+gQ 2 (ζ) + gQ 3 (ζ) = 0. ( 26
)
In this paper, S 1 (t) is computed using a vertex centered finite-volume approach and S 2 (t) using a finite-difference approach.

In order to discretize system [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF], the numerical two-dimensional domain Ω is discretized with a cartesian grid made of rectangular cells. We denote (S i ) 1≤i≤Nv the vertices of the rectangular cells. The dual cell C i is the volume control associated with the vertex S i . It is delimited in joining the barycenters of the 4 rectangles surrounding S i (see Fig 3). We obtain the associated dual mesh, made of N c rectangular cells (C i ) 1≤i≤Nc . In this vertex-centered approach, we associate any averaged finite-volume value with the center of the corresponding dual cell. Remark 3.1. Considering an initial mesh with N x + 2 and N y + 2 nodes respectively in the first and second direction, we have N x × N y interior vertices and therefore N c = N x × N y cells in the dual mesh (the introduction of the 2 extra nodes is for the sake of boundary conditions).

S i C i C j(3) C j(1) C j(2)
C j(4)

x y y

Figure 3: Initial cartesian mesh ( ) and associated dual cells ( ) for a regular mesh of constant steps δ x and δ y . Each interior nodes (•) is the center of a computational cell. The boundary nodes (•) are associated with boundary "ghost" cells.

Considering a cell C i , and neglecting the issue of boundary conditions for a while, we set:

• j(k) : the index of the cell C j(k) neighboring C i , for each k ∈ {1, . . . , 4},

• |C i | : area of C i .
• Γ ij : boundary interface defined by C i and C j .

• ℓ ij : length of Γ ij .

• n ij : unit normal to Γ ij pointing to C j .

• wn

i : the averaged value of any arbitrary quantity w on cell C i at time t n = nδ t .

Spatial discretization of the S 1 (•) component

Let us consider the following system, neglecting for a while the friction source term:

   ∂ t h + ∇ • (hV ) = 0, ∂ t (hV ) + ∇( 1 2 gh 2 ) + ∇ • (hV ⊗ V ) = -gh∇b (27)
It is well-known that a numerical imbalance problem is induced when the surface gradient term is splitted into an artificial flux gradient and a source term that includes the effect of the bed slope for a non-uniform bed. This can be eliminated for first order or second order schemes using for instance the celebrated hydrostatic reconstruction method [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. For higher-order schemes, such as the WENO schemes used in §3.2.3, it is not trivial to remove the artificial source. This imbalance problem can be solved by a reformulation relying on a deviation from an unforced but separately specified equilibrium state, as done for instance in [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF][START_REF] Rogers | Mathematical balancing of flux gradient and source terms prior to using Roe's approximate Riemann solver[END_REF]. This is the strategy adopted in this work, and recalled in the next section.

Pre-balanced formulation

Let us denote η = h 0 + ζ and q = (q x , q y ) = hV the horizontal discharge.

Following [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF], we remark that:

- g 2 ∇h 2 -gh∇b = - g 2 ∇(η 2 -2ηb) -gη∇b, (28)
and reformulate Eqs. ( 27) in an alternative way:

u t + ∇ • F(u, b) = S(u, b), (29) with u
=   η q x q y   F(u, b) =   q x q y uq x + 1 2 g(η 2 -2ηb) vq x uq y vq y + 1 2 g(η 2 -2ηb)   , (30)
and

S(u, b) =   0 -gηb x -gηb y   . (31)

First order finite volume discretization

On each cell C i , we denote ūn i = (η n i , ( qx ) n i , ( qy ) n i ) a piecewise constant approximation of the exact solution vector at time t n . Let us also denote bi a piecewise constant approximation of the topography. This approximate state is evolved at time t n + δt with the use of a finite volume scheme. Considering a first order forward Euler scheme with a δt time step for the timemarching algorithm, the fully discrete finite volume method to approximate the weak solutions of ( 29) is given under the following form:

|C i | ūn+1 i -ūn i δt + 4 k=1 ℓ ij(k) F s (ū * ij(k) , ū * j(k)i , b i , b j , n ij(k) ) = 0, (32)
with

F s (ū * ij , ū * ji , b i , b j , n ij ) = F (ū * ij , ū * ji , b ij , b ij , n ij ) -S c,ij .
In this last expression:

• F (ū * ij , ū * ji , b ij , b ij , n ij )
is the numerical flux through the interface between cells C i and C j . In the numerical validations shown in §3, we use a positive preserving relaxation-VFRoe scheme [START_REF] Berthon | A Positive Preserving High Order VFRoe Scheme for Shallow Water Equations: A Class of Relaxation Schemes[END_REF], consistent with the exact flux and satisfying the conservation property.

• ū * ij = (η * ij , h * ij u i ) and ū * ji = (η * ji , h * ji u j
) are faces values computed from the i ("left") and j ("right") Riemann states obtained with the following reconstructions:

b ij = max(b i , b j ), ∆ ij = max(0, b ij -η i ), b ij = b ij -∆ ij , (33)
h * ij = max(0, η i -b ij ), h * ji = max(0, η j -b ij ), ( 34
)
η * ij = h * ij + b ij , η * ji = h * ji + b ij . (35)
• S c,ij is a centered discretization of the topography source term (31), defined as follows:

S c,ij =   0 g η * ij + η * ji 2 (b i -b ij ) n ij   . ( 36
)
Remark 3.2. The resulting first order scheme is shown in [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF] to preserve motionless steady-states, even with the occurrence of dry states.

Higher order discretization

To reach higher order space discretization, we have implemented 3 rd and 5 th order accuracy WENO reconstructions (WENO3 and WENO5 in the following), following [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF], together with the weight splitting method for the 5 th order method [START_REF] Shi | A technique of treating negative weights in WENO schemes[END_REF] to treat the negative weights. Classically, given a quantity w and a cell C i , this approach provides, for all time t n , interpolated quantities at each of the 4 faces of the cell by performing dimension-wise reconstructions. The new values are then injected into the Riemann solver instead of the piecewise constant values.

As far as robustness is concerned, we aim at keeping the water height nonnegative. Such a feature can be achieved, while preserving the high order accuracy, under a more restrictive CFL-like condition following the recent approach introduced and applied in [START_REF] Xing | Positivity-preserving high order wellbalanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF] for discontinuous-Galerkin and WENO approaches. Let us briefly recall the flowchart of this strategy in the one dimensional case with constant step δ x , for the WENO5 reconstruction :

1. for a given cell C i and time step n, let us denote h n i-1 2 + and h n i+ 1

2

-the high order reconstructed values for the water height obtained from the WENO reconstruction at each face of the cell. We build a 4 th order polynomial h n i (x) , using the Hermite type reconstruction, such that

h n i (x i-1 2 ) = h n i-1 2 + , h n i (x i+ 1 2 ) = h n i+ 1 2 - 20 and
1 δx Cj h n i (s)ds = h n i , j = i -1, i, i + 1.
2. let S i = {x i k , k = 1, . . . , p} be the set of p Gauss-Lobatto quadrature points on the segment C i , and {ω k , k = 1, . . . , p} the quadrature weights on the segment [- 1 2 , 1 2 ]. For the 5 th order WENO reconstruction, we choose p = 4, as the quadrature has to be exact for integrals of 4 th order polynomial. We compute point values {h n i (x i k ), k = 1, . . . , p} to get m n i = min

x i k ∈Si h n i (x i k ).
3. we modify the polynomial h n i using a conservative linear scaling around the cell average :

ĥn i (x) = θ n i (h n i (x) -h n i ) + h n i , (37)
where

θ n i = min e -h n i m n i -h n i , 1
,
and e is a threshold used to define what is numerically called a dry cell.

In practice, we set e = 10 -7 .

we compute updated high order values ĥn

i-

1 2 + = ĥn i (x i- 1
2 ) and ĥn

i+ 1 2 -= ĥn i (x i+ 1
2 ) and use them within the finite volume scheme instead of

h n i-1 2 + and h n i+ 1 2 -.
Then, as shown in [START_REF] Xing | Positivity-preserving high order wellbalanced discontinuous Galerkin methods for the shallow water equations[END_REF], the resulting scheme is robust under the CFL condition δt δx (|u| + c) ∞ ≤ ω 1 , where (|u| + c) ∞ is some approximated maximum faces wave speed, computed within the Riemann solver. For the 4 points Gauss-Lobatto rule, we have ω 1 = 1/12. We do not give further details concerning this procedure. Note that the two dimensional implementation is slightly more complicated, even considering direction by direction WENO reconstruction to get the 4 faces high order reconstructed values. Indeed, we need to consider tensor products of Gauss and Gauss-Lobatto quadrature points and construct more polynomials in each cell. The reader is referred to the [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF] for the detailed implementation in the 2D case, straightforwardly applied here. Remark 3.3. Following [START_REF] Xing | Positivity-preserving high order wellbalanced discontinuous Galerkin methods for the shallow water equations[END_REF], and for regular grids, the resulting two-dimensional high order scheme is robust (i.e. the water height remains positive if it is positive initially) under the following CFL condition:

δt δx (|u| + c) ∞ + δt δy (|v| + c) ∞ ≤ ω 1 , (38)
where (|u|+c) ∞ and (|v|+c) ∞ are some approximated maximum faces wave speeds, computed within the Riemann solver respectively in the first and second direction. We have ω 1 = 1/6 for the WENO3 method and ω 1 = 1/12 for the WENO5.

We emphasize that the strategy detailed above aims at preserving the high-order accuracy of the WENO schemes, and this is the reason why such stringent restrictions on the time step are applied. These restrictions can lead to a significant increase of the computational time. In practice, to keep a high level of computational efficiency, the computations are performed with a relaxed CFL condition and we use the CFL condition [START_REF] Israwi | Large Time existence For 1D Green-Naghdi equations[END_REF] only when a preliminary calculation to the next time step produces negative water height. Note that another efficient approach is the recent one from [START_REF] Hu | Positivity-preserving method for highorder conservative schemes solving compressible Euler equations[END_REF], which ensures the preservation of the water height positivity without additional time step restriction, by considering linear combinations of the high-order and first-order fluxes in some particular cells, identified as about to generate negative values of h. The price to pay may be a local loss of accuracy near the wet/dry interfaces.

Remark 3.4. As shown in [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF], the source term discretization [START_REF] Hsiao | Tsunami-like solitary waves impinging and overtopping an impermeable seawall: Experiment and RANS modeling[END_REF] does not need to be modified to preserve well-balancing and consistency when high-order reconstructions are used.

Note that the additional friction term introduced in Eqs. [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF] (see Remark 2.6) can be embedded into our finite volume discretization in a natural way, following the modified Riemann solver approach introduced recently in [START_REF] Berthon | A numerical scheme for Shallow Water Equations with friction[END_REF]. The resulting scheme is shown to preserve the balance properties, the overall accuracy and the robustness under the CFL-like condition used for system [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF]. The reader is referred to [START_REF] Berthon | A numerical scheme for Shallow Water Equations with friction[END_REF] and [START_REF] Berthon | Asymptotic preserving scheme for the shallow water equations with source terms on unstructured meshes[END_REF] for the practical implementation.

Spatial discretization of the S 2 (•) component

The system ( 40) is solved at each time step using a classical finite-difference technique. The spatial derivatives are discretized using fourth-order formulae.

In the one-dimensional case, considering a regular grid of constant discretization steps δ 1 , we obtain for an arbitrary quantity w, at the node i:

(∂ 1 w) i = 1 12δ 1 (-w i+2 + 8w i+1 -8w i-1 + w i-2 ), (∂ 2 1 w) i = 1 12δ 2 1 (-w i+2 + 16w i+1 -30w i + 16w i-1 -w i-2 ). ( 39
)
These formulae are applied straightforwardly direction-wise and eventually combined to compute cross-direction derivatives. The computation of the terms I +

αT b diag -1 1 α gh∇ζ + h(Q 1 (V ) + gQ 2 (ζ)) + gQ 3 (ζ) and I + αT b diag -1 (gh∇ζ)
results in the resolution of sparse unsymmetric linear systems. For each term, the diagonal structure of the operator I + αT b diag involved in system (17) allows to consider two sparse linear systems of size (N x N y ) 2 instead of one system of size (2N x N y ) 2 . Additionally, as I + αT b diag is a constant operator, the associated matrix is build and LU -factorized in the compressed sparse row format (CSR) [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] in a preprocessing step at the beginning of each computation, and used through the whole computation at each time step (and substep with higher order time Runge-Kutta discretizations) without being modified. It results in a dramatic decrease of the computational cost when compared with the diagonal model [START_REF] Bonneton | Modelling of periodic wave transformation in the inner surf zone[END_REF], see §4.1 and Fig. 6. For the validations of §4, the factorization and the resolution of the resulting triangular linear systems are performed using the unsymmetric multifrontal method [START_REF] Davis | An unsymmetric-pattern multifrontal method for sparse LU factorization[END_REF].

Remark 3.5. As the solution operator S 2 (•) does not modify the value of h, the whole splitting method ( 24) is obviously robust under the CFL condition [START_REF] Israwi | Large Time existence For 1D Green-Naghdi equations[END_REF], using the additional robustness limitation (37).

Boundary conditions and time discretization

The boundary conditions for the hyperbolic and dispersive parts of the splitting are treated as in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]: we introduce some ghosts cells all around the computational domain, and suitable reflecting relations are imposed on the computed quantities. For the dispersive part, the boundary conditions are imposed by reflecting the coefficients associated to stencil points that are located outside of the domain (periodically for periodic conditions, evenly for Neumann conditions and oddly for Dirichlet conditions). Outgoing/absorbing boundary conditions are obtained with combining Neumann boundary conditions for both steps of the splitting method with a classical sponge layer approach. We have to slightly increase the computational domain to include the sponge layers, which length has to be calibrated from the incoming waves (generally 2 wavelengths). Periodic wave generation is achieved by using a generation/relaxation zone method [START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF]. We implemented the relaxation functions proposed in [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]. At the inlet boundary, we progressively impose on a one-wavelength long generation layer the targeted wave train using stream functions [START_REF] Fenton | A fifth order stokes theory for steady waves[END_REF]. This method is easy to implement and give very satisfying results, without spending too much time for calibration. Combined with a classical sponge layer approach it allows to generate incoming waves while allowing the possible reflected waves to be damped, and therefore to converge towards stationary states.

As far as time discretization is concerned, we choose to use explicit methods. The systems corresponding to S 1 and S 2 are integrated in time using a third or fourth-order SSP-RK scheme [START_REF] Gottlieb | Strong stability preserving high order time discretization methods[END_REF].

Wetting and drying

No particular treatment is needed to allow the occurrence of wetting and drying. The high order Finite-Volume approach detailed in §3.2 is robust under the suitable CFL-like condition and consequently preserves the positivity of the water height. The dispersive part of the computation does not modify the values of the water height. However we emphasize that the constant diagonal model ( 17) is derived under the implicit assumption that h b ≥ 0. To obtain a robust and consistent model allowing the occurrence of dry state, we slightly modify the definition [START_REF] Castro | Fully nonlinear long-waves models in presence of vorticity[END_REF] of h b as follows:

h b = max(h 0 -b, ǫ b ),
where ǫ b is a freely chosen threshold that only aims at preserving the regularizing effects of the operator I + αT b diag -1 when dry states occurs. It is important to insist on the fact that this manipulation does not create any additional consistency error in the model, since it is compensated by a similar change in the definition (15) of Q 3 (ζ). In the numerical validation of §4, we set ǫ b = 10 -1 . Additionally, our numerical investigations have shown that there is no need to suppress the dispersive effects near dry areas, the resulting numerical model being robust enough.

Wave breaking

Several criteria for the initiation of breaking can be found in the literature, see for instance [START_REF] Bjorkavag | Wave breaking in Boussinesq models for undular bores[END_REF][START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking, and runup. I: 1D[END_REF][START_REF] Okamoto | The relative trough Froude number for initiation of wave breaking: theory, experiments and numerical model confirmation[END_REF][START_REF] Schaeffer | A Boussinesq model for waves breaking in shallow water[END_REF]. Recently, a new approach specifically developed for our splitting numerical model has been introduced in [START_REF] Tissier | Nearshore Dynamics of Tsunami-like Undular Bores using a Fully Nonlinear Boussinesq Model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF].

In this approach, we identify wave fronts, and distinguish broken fronts and fronts likely to break from others, using the dissipation. Indeed, the energy dissipation forms peaks at the steepest parts of the wave fronts when shocks are forming. Then, an efficient criterion for initiation and termination of breaking, based both on the critical front slope [START_REF] Schaeffer | A Boussinesq model for waves breaking in shallow water[END_REF] and the Froude number across the fronts, is applied to determine dynamically the area where to locally switch in space and time between Green-Naghdi and NSW equations, using the first Finite-Volume sub-step as a predictor. After breaking, the wave fronts are handled as shocks by the NSW equations. Indeed, we suppress the S 2 (δt) step in the splitting approach, and only the hyperbolic part of the equations is solved for the wave fronts. The breaking wave dissipation is represented by shock energy dissipation [START_REF] Bonneton | Modelling of periodic wave transformation in the inner surf zone[END_REF]. The switching between the two systems is performed abruptly, without any smooth transition zone. In this way, wave propagation is governed by one given set of equations in each cell, and not by a non-physical mix of both sets. The transition does not generate any disturbances and therefore no numerical filtering is applied. The reader is referred to [START_REF] Tissier | Nearshore Dynamics of Tsunami-like Undular Bores using a Fully Nonlinear Boussinesq Model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] for implementation details and extensive validations. This method is implemented in our numerical model and used for the one dimensional test cases 4.2 and 4.4.

Considering the 2d case, we emphasize that an efficient 2d extension of the previous approach is far from trivial and is currently under study. A simple initiation criterion is introduced in [START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2DH improved Boussinesq equations[END_REF] to characterize breaking wave fronts can be used for 2d simulations. This criterion is based on the ratio of the surface elevation to the local water depth, and the solution operator S 2 (t) is suppressed at grid points where this ratio is greater than the threshold value 0.8. We emphasize that this wave-breaking criterion is particularly unsophisticated. Indeed, there is no efficient method for the switch back to S-GN equations, and it does not allow to follow the waves during their propagation, inducing limitations for the study of irregular waves trains over uneven bottoms. However our numerical investigations have shown that it was efficient enough for the breaking of solitary waves. This breaking wave criterion is also used in [START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF]. Some two-dimensional extensions of the more sophisticated ideas of [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] are left for future works, and we do not investigate here 2d simulations with breaking waves.

The 3 parameters constant diagonal model

For the numerical discretization of the 3 parameters constant diagonal model ( 23), the second time step S 2 (δt) is modified as follows:

         ∂ t h + I + γT b diag -1 θ∇ • (hT diag V ) = 0, ∂ t (hV ) - 1 α gh∇ζ + I + α(1 + θ)T b diag -1 1 α gh∇ζ +hQ 1 (V ) + ghQ 2 (ζ) + g(1 + θ)Q 3 (ζ) + θQ 4 (V ) = 0. ( 40
)
As h is now modified by the S 2 (δt) step, the computation of the terms I +

γT b diag -1 θ∇ • (hT diag V
) results in the resolution of an additional sparse unsymmetric linear system, at each time step (the fact that this linear system is time independent in this new model makes 2D computations possible, while a 2D generalization of [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF] for instance is not realistic). Let us remark that this additional computation obviously does not disturb the well-balanced property of the whole numerical method, as ∇ • (hT diag V ) vanishes whenever V = 0. Therefore, the motionless steady states at still preserved. However, we do not provide here a robust discretization of the 3 parameters constant diagonal model. Our numerical simulations have shown that the 3 parameters model can be successfully used for test cases involving wet/dry interfaces and run-up phenomena, providing that the dispersive step S 2 (δt) is suppressed in the vicinity of dry areas.

Numerical validations

In this section, we assess the accuracy and the dispersive properties of the new constant diagonal [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF]. We recall that the second order accuracy in time of the time splitting method is numerically confirmed in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], in which a time convergence analysis is performed. We also highlight that the 1d cases performed in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] are successfully reproduced with the new constant diagonal model (with a considerable gain of computational time). We do not include these results in the present study. Instead, we perform a convergence analysis of the space accuracy in the first test case. We then study the propagation of solitary waves over a composite beach, a solitary wave overtopping a seawall and the propagation of highly dispersive periodic waves over a submerged bar. The following cases are devoted to true validations in 2D configurations, through the comparison with experimental data for the shoaling of regular waves over uneven bottoms. Then, we focus on robustness, with a test case involving wetting and drying processes. Note that some numerical results produced by the 3-parameters constant diagonal model with optimized coefficients are shown in §4.5. Test cases 4.6, 4.7 and 4.8 have also been reproduced with the 3-parameters constant diagonal model without any noticeable improvements of the results.

In all the test cases, we use the WENO5 reconstruction, together with the 4 th order SSP-RK time marching algorithm, allowing to use a CFL number of 1.5. Our numerical investigations highlight that the time step restriction from the CFL condition of the finite-volume step is enough to ensure stability for the whole numerical method. For computational efficiency measurement, note that all the computations are performed on a single core of a 1.7 Gh Intel i5 dual-core processor.

Propagation of a solitary wave

It is known that for horizontal bottoms, the Green-Naghdi model (1) in the one-dimensional framework has exact solitary wave solutions given, in variables with dimensions, by

             h(x, t) = h 0 + a sech 2 (κ(x -ct)), u(x, t) = c 1 - h 0 h(x, t) , κ = √ 3a 2h 0 √ h 0 + a , c = g(h + H), (41)
Such solitary waves are also solution of the "constant diagonal" model up to a O(µ 2 ) remainder. Therefore, this family of solutions can be used as a validation tool for our present numerical scheme. In this test, we consider the following system :

   ∂ t h + ∇ • (hV ) = 0, I + T b diag ∂ t (hV ) + ∇ • (hV ⊗ V ) + gh∇ζ +h(Q 1 (ζ) + Q 2 (V )) + Q 3 (ζ) = R(ζ, V ), where R(ζ, V ) is the O(µ 2
) remainder obtained when one injects the solitary wave (41) into system [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF]; this residual is therefore given by an explicit formula deduced from [START_REF] Kazolea | Mathematical and Computational Modeling for the Generation and Propagation of Waves in marine and Coastal Environments[END_REF]. We consider the propagation of a solitary wave of relative amplitude a/h 0 = 0.2, in the x-direction, in a 2000 m long domain with a constant depth h 0 = 10 m. The initial surface and velocity profiles are centered at x 0 = 200 m; the solitary wave is therefore far from the boundaries of the domain in the x-direction (see Fig. 4). The width of the domain is arbitrary and does not influence the computation. The boundary conditions used here are therefore not important, and we choose for the sake of simplicity periodic boundary conditions in the x-direction and solid wall boundary conditions in the y-direction. The first computation is performed with a space discretization We show in Fig. 4 the free surface elevation of the propagating wave at t = 0 s, t = 30 s and t = 60 s. We observe a propagation without noticeable deformation, even after 1800 m. In order to quantify the numerical accuracy and assess the convergence of the numerical scheme, the numerical solution is computed for several cell sizes δ x , over a duration T = 1 s. Starting with δ x = 5 m and, we successively divide the space step by two, while keeping the time step small enough to ensure that the leading error order is provided by the spatial discretization. At t = 1 s " the relative errors E L2 (h) and E L2 (q x ) on the free surface elevation and the averaged discharge are computed using the discrete L 2 norm ||.|| 2 :

Number of cells

E L 2 (ζ) E L 2 (
E L2 (ζ) = ||ζ num -ζ sol || 2 ||ζ sol || 2 ; E L2 (hu) = ||(hu) num -(hu) sol || 2 ||(hu) sol || 2
where (ζ num , (hu) num ) are the numerical solutions and (ζ sol , (hu) sol ) denotes the analytical ones coming from [START_REF] Kazolea | Mathematical and Computational Modeling for the Generation and Propagation of Waves in marine and Coastal Environments[END_REF].

Results are gathered in Table 1 and Fig. 5, where E L2 (ζ) and E L2 (hu) are plotted against δ x in log scales, for the considered relative amplitude a/h 0 = 0.2. The convergence of our numerical scheme is clearly demonstrated. Furthermore, computing a linear regression on all points yields a slope equal to 3.9 for ζ and 3.8 for the discharge hu. These results are consistent with the numerical scheme, for which the leading order of the numerical error in space is given by the fourth order finite difference schemes used to discretize the dispersive part of the equations. To conclude this section, and to highlight and quantify the computational improvements obtained with the constant diagonal formulation, we compare the computational time (cpu-time) per time-step between the diagonal (12) and the constant diagonal (17) models, for an increasing number of cells and for the first order accuracy in space and time scheme. Two methods are successfully used to compute the solutions of the linear systems involved in the diagonal model:

1. a direct solver: the unsymmetric multifrontal method [START_REF] Davis | An unsymmetric-pattern multifrontal method for sparse LU factorization[END_REF]. We recall that for the diagonal formulation, as the discretization matrix of the operator I + αT diag is time-dependent, the LU factorization has to be computed at every time time step, whereas with the diagonal-constant model, the discretization matrix of I + αT b diag is build and LU -factorized as a preprocessing step at the beginning of the computation.

2. an iterative solver: the Bi-Conjugate Gradient Stabilized (BiCGStab) method, together with the ILUT preconditioner of SPARSKIT [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Again, for the the diagonal formulation, the discretization matrix has to be build at every time step but, thanks to the iterative nature of the solver, no LU factorizations are needed.

For both solvers, the matrix are stored in CSR format. Starting with N x = N y = 50, we progressively increase the number of cells to reach N x = N y = 600.

Resulting cpu-time per iteration are shown on Fig. 6. We clearly observe that the use of the direct solver for the diagonal formulation rapidly leads to very large computational costs, due to the LU-factorization. That's the reason why we turn on using the iterative solver. Even with the iterative solver, we can rapidly observe the gain obtained using the diagonal-constant model [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF]. 12) and the diagonalconstant model [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF] for an increasing number of cells.

Solitary wave propagation over a composite beach

We study in this case the propagation of a solitary wave over a composite beach, simulating geometrical dimensions of the Revere Beach. The original experiment was performed in a tank by the U.S. Army Corps of Engineers at the Coastal Engineering Research Center in Vicksburg, Mississippi. The constructed beach consists of three piece-wise linear segments, terminated with a vertical wall on the left. The slope s of the topography is define as follows:

s(x, y) =        0 if x ≤ 15.
) 42
The schematics of the beach is shown on Fig. 7. In the following, our numerical results are compared with experimental data measured from gages 5, 7 and 9, respectively located at x 5 = 15.04 m, x 7 = 19.4 m and x 9 = 22.33 m (exactly at the locations corresponding the slope variations). Three wave configurations were studied in the original experiment, with targeted solitary waves of relative amplitudes 0.05 (♯ A), 0.3 (♯ B) and 0.7 (♯ C). Only ♯ B and ♯ C are simulated in the following. Indeed ♯ A, which involves a very small relative amplitude wave can be quite accurately reproduced with a classical non-dispersive shallow water model, while the two other cases clearly involve some dispersive effects. We provide the solitary wave of targeted height, centered at x = 0 as the initial condition, using formula [START_REF] Kazolea | Mathematical and Computational Modeling for the Generation and Propagation of Waves in marine and Coastal Environments[END_REF] and we observe the propagation over the beach, reflection on the vertical wall before traveling back to the inlet boundary. We place a reflective wall boundary on the right and lateral boundaries of the domain, simulating the vertical wall, and a sponge layer to absorb the possibly small amplitudes reflected waves at the left boundary. Experimental data are provided as time series of the wave elevation at several gages located along the wave flume. We show on Fig 10 the comparison between data and computed results for ♯ B and ♯ C, at gages 5, 7 and 9. Note that we obtain similar qualitative results for the remaining gages. Results are obtained with 500 cells in the x direction. Note that the solitary wave encounters breaking during propagation.

Shoaling of solitary waves

Let us now investigate the dispersive properties of the new model through the study of the nonlinear shoaling, using the data issued from a laboratory study performed at the LEGI (Grenoble, France). We consider in this test a 36m channel with constant bed slope and solitary waves generated at the left inflow boundary. Measurements of the free surface are available in the vicinity of the breaking point during the simulation, at several wave gages. We consider 4 series of experiments, involving increasing wave's relative amplitudes, from a/h 0 = 0.096 to a/h 0 = 0.534. We report the numerical results obtained with the diagonal constant model on Fig. 8, together with the experimental data taken from [START_REF] Guibourg | Modélisation numérique et expérimentale des houles bidimensionnelles en zone cotière[END_REF]. A very good matching is observable, even for large wave's amplitude. These results assess the good accuracy of our model in the reproduction of the shoaling process.

4.4 Wave overtopping a seawall. In this case we highlight the robustness of our model, through the simulation of a solitary wave overtopping a seawall, following the experiments performed in the Tainan Hydraulic Laboratory, National Cheug Kung University, Taiwan [START_REF] Hsiao | Tsunami-like solitary waves impinging and overtopping an impermeable seawall: Experiment and RANS modeling[END_REF]. We consider a 22 m long wave flume, with a water depth of 0.2 m and a propagating solitary wave of 0.07 m height. The topography including the seawall, is illustrated on Fig. 11. In this case, the propagating solitary wave breaks on the sloping beach before reaching the wall. In the experiment, several wave gages were placed along the basin (see [START_REF] Hsiao | Tsunami-like solitary waves impinging and overtopping an impermeable seawall: Experiment and RANS modeling[END_REF] for exact locations) and during the simulation we extract time series of the water free surface at corresponding locations. The comparison between experimental data and numerical results are gathered on Fig. 11. We use 400 cells in the propagation direction, giving δ x = 0.03 m. A small friction term, corresponding to a Manning coefficient n = 0.012 is added. We obtain very satisfying results, especially when compared with those introduced in [START_REF] Hsiao | Tsunami-like solitary waves impinging and overtopping an impermeable seawall: Experiment and RANS modeling[END_REF] and obtained with a Volume Of Fluid (VOF) method for the Reynolds averaged Navier-Stokes equations. We highlight that the use of the WENO5 reconstruction, together with the positivity preserving limitation described in section §3.2.3 leads to similar results than those shown in [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF], but using fewer cells. In this test, we highlight the ability of the constant-diagonal and 3 parameters constant-diagonal models to describe the propagation and the interaction of highly dispersive waves. We simulate the propagation of regular periodic waves over a submerged bar, using the set-up introduced in [START_REF] Beji | Experimental investigation of wave propagation over a bar[END_REF], and first used as a test in [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF], see Fig. 12. Periodic waves are generated at the left boundary, with an amplitude of 0.01 m, a time period of 2.02 s and mean water depth h 0 = 0.4 m.

Periodic waves propagation over a submerged bar

When the incident wave encounters the upward part of the bar, it shoals and steepens, which generates higher-harmonics as the nonlinearity increases. These higher-harmonics are then freely released on the downward slope, and become deep-water waves behind the bar. Comparisons are performed between the one-parameter optimized constantdiagonal model (α = 1.159), the three-parameters constant-diagonal model (α = 1.024, θ = 0.188, γ = 0.112) and the data taken from the experiment, for the last four wave gauges. Time series of the free surface elevation at the four last wave gauges of the experiment are plotted on Fig. 13, using δ x = 0.03 m. Note that the constant-diagonal model is not able to provide an accurate free surface evolution at the wave gauge 11. This last gauge is the most discriminating one as the higher-harmonics are completely released and can be regarded as highly dispersive waves. We observe the improvements obtained with the 3 parameters constant-diagonal model at this last wave gauge. Note that the use of the 3-parameters model induces an additional computational cost of approximately 20%.

Periodic waves propagation over a semi-circular shoal

This test reproduces one 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 effect induced by a semicircular shoal on wave trains of different periods. The experiments were performed in a basin of 25.6 m long and 6.096 m wide. Its middle portion consists in a semi-circular shoal leading the water depth to decrease from h 0 = 0.4572 m (at the inlet boundary) to 0.1524 m at the end of the tank. The bottom topography is described as follows: (43) with Γ(y) = y(6.096 -y). This has become a standard benchmark for the validation of dispersive numerical models, to test nonlinear refraction and diffrac-tion. lateral boundaries. For test ♯ B and ♯ C, we use a relatively coarse mesh of δ x = δ y = 0.12 m. For test ♯ A we use δ x = 0.05 m and δy = 0.1 m, as it is a more demanding test case. We perform the propagation during 100 s and we extract numerical data during the last 25 s. To highlight the efficiency of the constant diagonal numerical model, approximately 25 minutes are needed to perform the whole 100 s of propagation for test ♯ A and 6 minutes for tests ♯ B and ♯ C. The water surface elevation is measured along the centerline of the basin and harmonic analysis (DFT) was performed to obtain the amplitude of frequency components.The time series of the surface elevation, at each cell centre along the midline, is analyzed in the frequency domain to obtain the first, second and third harmonic amplitudes and we consider their spatial evolution. The numerical results are compared to the experimental data on Fig 18 . We can observe some discrepancies between the numerical results and the experimental data for test ♯ C. However, similar results are observable on previous studies (see for instance [START_REF] Bingham | Velocity potential formulations of highly accurate Boussinesq-type models[END_REF][START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF] ) and are usually attributed to the shorter evolution distance for this test or to the presence of free reflected waves.

z(x, y) =    0 if 0 ≤ x ≤ 10.67 -Γ(y), ( 10

Periodic waves propagation over an elliptic shoal

In this test, we reproduce the experiment carried out in [START_REF] Berkhoff | Verification of numerical wave propagation models for simple harmonic linear water waves[END_REF] to study the diffraction of 2D monochromatic wave train over a varying bottom. The wave tank is 20 m wide and 22 m long. The bathymetry consists of an elliptic shoal built on a ramp of constant slope, forming a 20 • angle with the y axis (see Fig. 19). Introducing the rotated coordinates These sections provide a good coverage of the waves propagation in the experiment. To obtain these numerical results, 50s of propagation are computed, with δ x = 0.065 m and δ y = 0.1 m. Time series of the surface elevation are obtained at several locations along the respective sections, between t = 40 s and t = 50 s. The signal is then analyzed with the zero up-crossing method to isolate single waves and compute mean wave elevation, normalized by the incoming wave amplitude. We obtain satisfying results, comparable with previous and recent studies (see [START_REF] Kazolea | Mathematical and Computational Modeling for the Generation and Propagation of Waves in marine and Coastal Environments[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF][START_REF] Tonelli | Hybrid finite volume -finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]). Approximatively 38 minutes are needed to perform the 50 s of simulation. Note that his test case has also be runned with the 3 parameters model, without any real improvements in the quality of the results.

Conclusion

We proposed here a new strongly nonlinear shallow water model with the same accuracy as the standard Green-Naghdi (or Serre, or fully nonlinear Boussinesq) equations. This model has a remarkable structure that allowed us to implement a very computationally efficient splitting method, based on an hybrid finitevolumes and finite-differences approach. A 5 th order WENO reconstruction is used within the finite volume scheme, whereas fourth order finite-differences are used for the finite-difference part, together with a fourth-order SSP-RK time marching algorithm. An high order accuracy preserving additional limitation is applied to ensure the preservation of the water height positivity. We validated this approach with several experimental data and showed its robustness by using it on complex 2D configurations, including waves propagating and shoaling, diffraction and refraction or run-up and overtopping. The three parameters model appears to be useful in challenging situations, when high-order harmonics are released but does not provide any significant improvements otherwise. The computational efficiency and robustness of this code opens many perspectives for applications in coastal oceanography, and could also be improved in various directions: -The 2D treatment of wave breaking we use here allows us to handle stiff configurations, but it is very rough. The approach used in the 1D case in [START_REF] Tissier | Nearshore Dynamics of Tsunami-like Undular Bores using a Fully Nonlinear Boussinesq Model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] should be generalized to the 2D case. Another interesting perspective -especially for the description of hydraulic jumps -would be to implement the new strategy of [START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF] which provides a description of the size of the bore (while in the classical description, it is a shock of zero width).

-The modeling of rip currents and more generally of other wave-current interactions requires that one takes into account rotational effects. A generalization of the physical model proposed here should be developed to handle rotational effects, and our splitting approach should be implemented on this model. Finally, we insist that our new models open perspectives beyond the numerical approach used here; other numerical methods can indeed take advantage of their constant diagonal structure. Discontinuous-Galerkin [START_REF] Eskilsson | Spectral/hp discontinuous Galerkin methods for modelling 2D Boussinesq equations[END_REF], Finite-Volume [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF] or Finite-element [START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] approaches, allowing the discretization on unstructuredmeshes, have for instance been implemented on weakly nonlinear Boussinesq models but the gain in computational time made possible by the constant diagonal Green-Naghdi systems should allow their generalization to the fully nonlinear case. This is of great interest in situations where weakly nonlinear models fail to capture correctly the evolution of the waves (see for instance [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]). The main difficulty is related to the stable discretization of higher order derivatives but the fact that the new constant diagonal models do not involve any third order derivatives on the flow variables allows interesting simplifications of the associated weak formulation. Also, it may lead to considerable computational savings in the use of higher-order polynomial interpolation. This is investigated, in the 1d case, in [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF].

Figure 1 :

 1 Figure 1: Sketch of the domain

Figure 2 :

 2 Figure 2: Instabilities in the case of uncontrolled third order derivatives

Figure 4 :

 4 Figure 4: Propagation of a solitary wave over a flat bottom: snapshots of the 2D free surface at t = 0 s and of 1D profiles at several times, a/h 0 = 0.2 and h 0 = 10 m

Figure 5 :

 5 Figure 5: Propagation of a solitary wave over a flat bottom: L 2 -error on the water height and the discharge for a/h 0 = 0.2, h 0 = 10 m.

Figure 6 :

 6 Figure 6: Propagation of a solitary wave over a flat bottom: comparison between the cpu-time per time step between the diagonal model (12) and the diagonalconstant model[START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF] for an increasing number of cells.

Figure 7 :

 7 Figure 7: Solitary wave on a composite beach -Set up of the composite beach with vertical wall on the right boundary, and initial condition for the case a/h 0 = 0.7.

Figure 8 :

 8 Figure 8: Guibourg test case : comparison between computed (solid lines) and experimental (dots) time series of total free surface at several gauges before the breaking point.

Figure 9 :Figure 10 :

 910 Figure 9: Wave overtopping over a seawall: sketch of the topography and locations of the wave gages.

Figure 11 :

 11 Figure 11: Wave overtopping over a seawall: time series of the free surface. Comparison with experimental data at wave gages for type 1 solitary wave.

Figure 12 :

 12 Figure 12: Experimental set-up and locations of the wave gauges.

Figure 13 :

 13 Figure 13: Periodic waves propagation over a submerged bar: time series of the free surface. Comparison between numerical results (-) with experimental data at wave gauges (o).

  .67 -G(y) -x)/25 if 10.67 -Γ(y) ≤ x ≤ 18.29 -Γ(y), 0.30480 if 18.297 -Γ(y) ≤ x.

Figure 14 :Figure 15 :

 1415 Figure 14: Periodic waves propagation over a semi-circular shoal -3D view of the topography

  x r = x cos(20 • ) -y sin(20 • ) , y r = x sin(20 • ) + y cos(20 • ), the topography is given by the formula z = z b + z s , where:z b (x, y) =(5.82 + x r )/50 if x r ≥ -5.82 , 0 elsewhere , z s (x, y) =

Figure 18 :

 18 Figure18: Propagation of periodic waves over a semi-circular shoal, case A, B, C: perspective view of the surface elevation along the centerline at t = 100 s (on the left) and comparison of the computed and experimental results for the wave amplitudes for the 1 st , 2 nd and 3 rd harmonics along the centerline (on the right).

Figure 19 :

 19 Figure 19: Periodic wave propagation over an elliptic shoal -3D view of the topography

Figure 20 :

 20 Figure 20: Periodic wave propagation over an elliptic shoal -Comparison of the computed mean wave height with the experimental data along sections 2, 3, 5, and 7

Figure 21 :

 21 Figure 21: Solitary wave over a conical island -Time series of the surface elevation at 4 gauges, for the case a = 0.1h 0 (on the left) and a = 0.2h 0 (on the right), obtained with both constant-diagonal GN and NSW equations.

Figure 22 :

 22 Figure 22: Solitary wave over a conical island -Run-up around the conical island. Measured maximum run-up (red dots) and computed inundated area (blue circles) for Case B (left) and Case C (right). The solid line stands for the initial shoreline.

Figure 23 :

 23 Figure 23: Solitary wave over a conical island -Case C: free surface at times t=5, 8, 9, 11, 11.5 and 13s.

Table 1 :

 1 Propagation of a solitary wave over a flat bottom-Convergence study: relative L 2 -error table for the conservative variables

	hu)

Or more precisely solving linear systems arising from the discretization of (1+µαT)F = G.

Acknowledgments

The authors would like to thank the referees for their suggestions that helped us improve the manuscript, as well as Prof P. Bonneton, Dr. M. Ricchiuto, Dr. M. Tissier and Dr. F. Chazel for fruitful discussions and Prof. A. Delis, S.-C. Hsiao and T.-C. Lin for providing some experimental data. The authors would also like to acknowledge additional financial and scientific support of the French INSU-CNRS (Institut National des Sciences de l'Univers-Centre National de la Recherche Scientifique) program LEFE-MANU (Méthodes Mathématiques et Numériques) and from the ANR BLANC Bond. D. L. also acknowledges support from the ANR BLANC Dyficolti.

The periodic incoming wave train is generated using a generation zone of length 5 m, and we applied a sponge layer upstream and downstream to absorb the outgoing and reflected wave. We apply solid wall boundary conditions at the

Tsunami wave on a conical island

We now study the propagation of a solitary wave over a conical island, based on the laboratory experiments performed in [START_REF] Liu | Runup of solitary waves on a circular island[END_REF]. Several authors have used this test as a benchmark to study the run up phenomena, see for instance [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF]. The basin dimensions are [0, 25] m × [0, 30] m, and we use a regular mesh with δ x = 0.1 m and δ y = 0.12 m. We compute 30 s of propagation of the wave for 2 cases, corresponding to ♯ B and ♯ C of the experiment. In these 2 cases, we generate a solitary wave, using formulae (41) at the left inlet boundary, with relative amplitudes respectively of a/h 0 = 0.1 and a/h 0 = 0.2, with a mean water depth h 0 = 0.32 m.. Denoting r the distance from the point (x 0 , y 0 ) = (12.96, 13.80), we consider an ideal island defined as follows:

We follow the evolution of the solitary wave after the impact at the front side of the island, and observe that all the features of the propagation are properly reproduced, according to the numerical results found in the references above. We can clearly observe on Fig. 23 the run up issuing from the reflection at the front of the island and then the separation in two secondary waves surrounding the island, finally colliding at the rear side. Time series of the surface elevation are available for the experiment at several gauges around the island. We show on Fig. 21 22. We obtain satisfying results, especially concerning the magnitude of the run-up and backwash at both side of the island. Note that the more important discrepancies are observable at gauge 9.

For ♯ C for instance, the run-up and run-down are slightly over-estimated and the subsequent free surface oscillations are not captured. This phenomena is consistent with other numerical results shown in the literature [START_REF] Fuhrman | Simulation of nonlinear wave run-up with a high-order Boussinesq model[END_REF][START_REF] Yamazaki | Depth-integrated, nonhydrostatic model for wave breaking and run-up[END_REF][START_REF] Tonelli | Finite volume scheme for the solution of 2D extended Boussinesq equations in the surf zone[END_REF][START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF]. These discrepancies after the initial run-up may be largely explained by the fact that the front sides of the experimental solitary waves were generated with more accuracy than the rear side, which included a spurious tail [START_REF] Liu | Runup of solitary waves on a circular island[END_REF][START_REF] Fuhrman | Simulation of nonlinear wave run-up with a high-order Boussinesq model[END_REF]. To highlight the efficiency of the constant diagonal formulation, note that it takes approximatively 8.5 minutes to compute the 30 s of propagation, for ♯ C, on a single core of an Intel i5 processor. With the regular CFL and the chosen mesh, the time step varies around 0.05 s but can decrease to slightly less than 0.005 s when the additional limitation of §3.2.3 is activated to ensure robustness. Running the same test case and configuration for the NLW equations leads to 360 s of computation, allowing us to estimate the cost of the dispersive step to approximately 30 % of the whole computational time. Note finally that we also investigated this case with the 3 parameters model without noticing any improvements.