Arnaud Duran
email: arnaud.duran@math.univ-toulouse.fr

Fabien Marche
email: fabien.marche@math.univ-montp2.fr

A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes

Keywords: Green-Naghdi equations, discontinuous Galerkin, high-order schemes, free surface flows, Shallow Water Equations, dispersive equations

In this paper, we introduce a discontinuous Finite Element formulation on simplicial unstructured meshes for the study of free surface flows based on the fully nonlinear and weakly dispersive Green-Naghdi equations. Working with a new class of asymptotically equivalent equations, which have a simplified analytical structure, we consider a decoupling strategy: we approximate the solutions of the classical shallow water equations supplemented with a source term globally accounting for the non-hydrostatic effects and we show that this source term can be computed through the resolution of scalar elliptic second-order sub-problems. The assets of the proposed discrete formulation are: (i) the handling of arbitrary unstructured simplicial meshes, (ii) an arbitrary order of approximation in space, (iii) the exact preservation of the motionless steady states, (iv) the preservation of the water height positivity, (v) a simple way to enhance any numerical code based on the nonlinear shallow water equations. The resulting numerical model is validated through several benchmarks involving nonlinear wave transformations and run-up over complex topographies.

Introduction

The mathematical modelling and numerical approximations of free surface water waves propagation and transformations in near-shore areas has received a lot of interest for the last decades, motivated by the perspective of acquiring a better understanding of important physical processes associated with the nonlinear and non-hydrostatic propagation over uneven bottoms. Great improvements have been obtained in the derivation and mathematical understanding of particular asymptotic models able to describe the behavior of the solution in some physical specific regimes. A recent review of the different models that may be derived can be found in [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF]. In this work, we focus on the shallow water and large amplitude regime: the water depth h 0 is assumed to be small compared to the typical wave length λ:

(shallow water regime) µ := h 2 0 λ 2 ≪ 1.
while there is no assumption on the size of the wave's amplitude a:

(large amplitude regime) ε := a h 0 = O(1).

Under this regime, the classical Nonlinear Shallow Water (NSW) equations can be derived from the full water waves equations by neglecting all the terms of order O(µ), see for instance [START_REF] Lannes | The water waves problem: mathematical analysis and asymptotics[END_REF]. This model is able to provide an accurate description of important unsteady processes in the surf and swash zones, such as nonlinear wave transformations, run-up and flooding due to storm waves, see for instance [START_REF] Bonneton | Modelling of periodic wave transformation in the inner surf zone[END_REF],

but it neglects the dispersive effects which are fundamental for the study of wave transformations in the shoaling area and possibly slightly deeper water areas. Keeping the O(µ) terms in the analysis, the corresponding equations have been derived first by Serre [START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF] in the horizontal surface dimension d = 1, by Green and Naghdi [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF] for the d = 2 case, and have been recently mathematically justified in [START_REF] Alvarez-Samaniego | Large time existence for 3d water-waves and asymptotics[END_REF].

As far as numerical approximation is concerned, it is only recently that the Green-Naghdi (GN) equations really have received attention and several methods have been proposed, mainly in the d = 1 case. Let us mention the Finite-Difference (FD) approach of [START_REF] Antunes | Boussinesq and Serre type models with improved linear dispersion characteristics: Applications[END_REF], Finite-Volume (FV) approaches of [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. I: Model development and analysis[END_REF][START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. II: Boundary conditions and validation[END_REF][START_REF] Favrie | A robust and rapid numerical method for dispersive models admitting a lagrangian : application to Serre-Green-Naghdi equations for long free surface gravity waves[END_REF], the pseudospectral approach of [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for fully-nonlinear 1d Serre equations[END_REF], continuous Finite Elements (FE) methods in [START_REF] Mitsotakis | On the Galerkin/Finite-Element Method for the Serre equations[END_REF][START_REF] Mitsotakis | A modified galerkin/finite element method for the numerical solution of the Serre-Green-Naghdi system[END_REF], and discontinuous Galerkin (DG) approaches (possibly mixed with FEM) in [START_REF] Panda | Discontinuous Galerkin methods for solving Boussinesq-Green-Naghdi equations in resolving non-linear and dispersive surface water waves[END_REF], [START_REF] Maojin | High order well-balanced CDG-FE methods for shallow water waves by a Green-Naghdi model[END_REF], or more recently [START_REF] Dong | A reconstructed central discontinuous galerkin-finite element method for the fully nonlinear weakly dispersive Green-Naghdi model[END_REF] for a reconstructed-DG method. In the d = 2 case, we find pioneering FD approaches [START_REF] Antunes Do Carmo | Numerical solution of the generalized Serre equations with the Mac-Cormack finite-difference scheme[END_REF], FV methods [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF], hybrid FV-FD methods [START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF][START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF] and a pseudo-spectral method in the rotating case [START_REF] Pearce | A pseudo-spectral algorithm and test cases for the numerical solution of the two-dimensional rotating Green-Naghdi shallow water equations[END_REF]. Let us also mention the recent non-hydrostatic models [START_REF] Choi | An efficient curvilinear non-hydrostatic model for simulating surface water waves[END_REF] and [START_REF] Aissiouene | A combined finite volume -finite element scheme for a dispersive shallow water system, Networks and heterogeneous media[END_REF] or the so called high-level GN equations [START_REF] Zhao | Application of higher-level GN theory to some wave transformation problems[END_REF] After having developed hybrid FV-FD formulations based on a temporal splitting approach (see [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][START_REF] Bonneton | Recent advances in Serre-Green-Naghdi modelling for wave transformation, breaking and runup processes[END_REF] for the d = 1 case and [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] for the d = 2 case on structured cartesian grids), we have recently introduced in [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF] an approach based on a full DG method in the d = 1 case. In this approach, we decouple the hyperbolic and elliptic parts of the model by computing the solutions of the NSW equations supplemented by an additional algebraic source term, which fully accounts for the whole dispersive correction, and which is itself obtained from the resolution of auxiliary elliptic and coercive linear second order problems. With this approach we have obtained:

♯ more flexibility, as the proposed discrete formulation handles an arbitrary order of accuracy in space and, although limited in the d = 1 case in [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF], it conceptually generalizes to the dimension d = 2 with unstructured meshes. It can moreover be straightforwardly used to enhance any numerical code based on the NSW equations, ♯ more efficiency, as we built this approach on the new set of GN equations issued from [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] which allows to perform the corresponding matrix assembling and factorization in a preprocessing step, leading to considerable computational time savings.

Note that similar decoupling strategies have recently been used to implement some hybrid approaches based on the classical GN equations (see [START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF] for an hybrid FV-FD scheme in the d = 2 case on cartesian meshes and [START_REF] Filippini | A flexible genuinely nonlinear approach for wave propagation, breaking and runup[END_REF] for a FV-FEM scheme in the d = 1 case. Although full DG or Spectral-Element methods have been developed respectively for high order Boussinesq equations in [START_REF] Engsig-Karup | Nodal DG-FEM solution of high-order Boussinesq-type equations[END_REF] and fully nonlinear water waves in [START_REF] Engsig-Karup | A stabilised nodal spectral element method for fully nonlinear water waves[END_REF], there is still no studies, up to our knowledge, that allows to approximate the solutions of GN equations on arbitrary unstructured meshes in the multidimensional case: this is the main purpose of this work.

Starting from the new model recently introduced in [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF], we consider a Discontinuous Galerkin (DG) formulation and generalize to the d = 2 case the decoupling strategy of [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF]:

♯1 we compute the flow variables using a RK-DG approach [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems[END_REF] by approximating the solutions of the hyperbolic set of NSW equations supplemented with an additional source term that accounts for the whole dispersive correction, ♯2 the dispersive source term is obtained from the computation of auxiliary scalar linear elliptic problems of second order, approximated through a mixed formulation and Local Discontinuous Galerkin (LDG) stabilizing interface fluxes [START_REF] Cockburn | The Local Discontinuous Galerkin method for time-dependent convection-diffusion systems[END_REF][START_REF] Castillo | An a priori error analysis of the Local Discontinuous Galerkin method for elliptic problems[END_REF].

This results in a global formulation of arbitrary order of accuracy in space, which is shown to preserve both the non-negativity of the water height and the motionless steady states up to the machine accuracy. These properties are particularly important in the study of propagating waves reaching the shore. Note that while this multidimensional strategy is clearly promising in terms of flexibility, the emphasize is put here on efficiency. First, the use of the diagonal constant model issued from [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] allows to considerably simplify the computations associated with the elliptic sub-problems. Indeed, the time evolutions of the velocity vector's components are here fully decoupled, thanks to the simplified analytical structure of the model. The problem therefore reduces to the approximation of scalar problems without any third order derivatives, instead of the more complex vectorial ones obtained with the classical GN equations. The corresponding matrix can be assembled and the associated LU factorization stored in a pre-processing step. Secondly, the use of a nodal approach [START_REF] Hesthaven | Nodal high-order methods on unstructured grids -i. time-domain solution of maxwell's equations[END_REF] together with the pre-balanced formulation of the hyperbolic part of the model (see for instance [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF][START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF]), allows to combine (i) an efficient quadrature free treatment for the integrals which are not involved into the equilibrium balance, (ii) a quadrature-based treatment with a lower computational cost, needed to exactly compute the surface and face integrals involved in the preservation of the steady states at rest, (iii) a direct nodal product method, in the spirit of pseudo-spectral methods, for the strongly nonlinear terms occurring in the source terms of the elliptic sub-problems. This remainder of this work is organized as follows: we describe the mathematical model and the notations in the next section. Section 3 is devoted to the introduction of the discrete settings and the DG formulations for both hyperbolic and elliptic sub-problems. This approach is then validated in the last section through convergence analysis and comparisons with data taken from experiments with several discriminating benchmark problems.

The physical model

Let us denote by x = (x, y) the horizontal variables, z the vertical variable and t the time variable. In the following, ζ(t, x) describes the free surface elevation with respect to its rest state, h 0 is a reference depth, -h 0 + b(x) is a parametrization of the bottom and h := h 0 + ζ -b is the water depth, as shown on Figure 1. Denoting by U hor the horizontal component of the velocity field in the fluid domain, we define the vertically averaged horizontal velocity

v = (u, v) ∈ R 2 as v(t, x) = 1 h ζ -h0+b
U hor (t, x, z)dz, and we denote by q = hv the corresponding horizontal momentum.

z 0 72-h 0 z = ζ(t, x) z = -h 0 + b(x) x = (x, y)
Figure 1: Sketch of the domain

The Green-Naghdi equations

According to [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], the Green-Naghdi equations may be written as follows:

∂ t ζ + ∇ • (hv) = 0, (I + T [h, b]) ∂ t v + (v • ∇)v + g∇ζ + Q 1 [h, b](v) = 0, (2)
where 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), (3)
Q 1 [h, b](w) = -2R 1 [h, b](∂ 1 w • ∂ 2 w ⊥ + (∇ • w) 2 ) + R 2 [h, b]((w • ∇)∇b), (4)
(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, (5)
R 2 [h, b]w = 1 2h ∇(h 2 w) + w∇b, (6)
and I stands for the 2 × 2 identity matrix. We also recall in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] that the dispersion properties of (2) can be improved by adding some terms of order O(µ 2 ) to the momentum equation, which consequently does not affect the accuracy of the model. An asymptotically equivalent enhanced family of models parametrized by α > 0 is given by

   ∂ t ζ + ∇ • (hv) = 0, (I + αT [h, b]) ∂ t v + (v • ∇)v + α -1 α g∇ζ + 1 α g∇ζ + Q 1 [h, b](v) = 0. ( 7
)
or alternatively written in (h, hv) variables:

   ∂ t h + ∇ • (hv) = 0, (I + αT[h, b]) ∂ t (hv)+∇•(hv ⊗ v) + α-1 α gh∇ζ + 1 α gh∇ζ + hQ 1 [h, b](v) = 0, (8)
where we have introduced the operator T[h, b] defined as follows

T[h, b]w = hT [h, b]( w h ).
As shown in [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF], (8) can be recast as

       ∂ t h + ∇ • (hv) = 0, (9a) ∂ t (hv) + ∇ • (hv ⊗ v) + gh∇ζ + D o = 0, (9b)
(I +αT[h, b])(D o + 1 α gh∇ζ) = 1 α gh∇ζ + hQ 1 [h, b](v). ( 9c
)
which highlights the fact that the dispersive correction of order O(µ) only acts as a source term D o in (9b), and is obtained as the solution of an auxiliary second-order elliptic sub-problem (9c).

Remark 1. These formulations have two main advantages:

1. They do not require the computation of third order derivatives, while this is necessary in the standard formulation of the GN equations, 2. The presence of the operator (I +αT[h, b]) -1 makes the models robust with respect to high frequency perturbations, which is an interesting property for numerical computations.

Remark 2. These formulations have two main drawbacks:

1. Solving linear systems arising from discrete formulations of (8) may be computationally expensive as I + αT[h, b] is a matricial second order differential operator acting on two-dimensional vectors. This structure entails a coupling of the time evolutions of the two components of hv through (9b).

I + αT[h, b

] is a time dependent operator (through the dependence on h) and this is of course the same for any associated discrete formulations: the corresponding matrices have to be assembled at each time step or sub-steps.

Remark 3. It is known that equations (2) exactly conserve the energy

E(t) = R 2 (ζ(t) 2 + h(t)|v(t)| 2 + h(t)T [h(t), b]v(t) • v(t)).
This cannot be expected for [START_REF] Antunes | Boussinesq and Serre type models with improved linear dispersion characteristics: Applications[END_REF] when α > 1: the conservation of energy is satisfied only up to O(µ 2 ) terms, which is the precision of the model.

A Green-Naghdi model with a simplified analytical structure

In [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF], some new families of models are introduced to overcome these drawbacks, without loosing the benefits listed in Remark 1. In particular, defining a modified water depth at rest (which therefore does not depend on time):

h b = max(h 0 -b, ε 0 ) = max(h -ζ, ε 0 ),
where ε 0 is a numerical threshold introduced to account for possible dry areas in a consistent way, the one-parameter optimized constant-diagonal GN equations read as follows:

   ∂ t h + ∇ • (hv) = 0, 1 + αT[h b ] ∂ t (hv) + ∇ • (hv ⊗ v) + α -1 α gh∇ζ + 1 α gh∇ζ + Q[h, b](ζ, v) = 0, (10)
where for all smooth enough R-valued function w:

T[h b ]w = - 1 3 ∇ • (h 3 b ∇ w h b ), (11)
and

Q[h, b](ζ, v) = h(Q 1 [h, b](v) + gQ 2 [h, b](ζ)) + gQ 3 [h, h b ] 1 + αT[h b ] -1 (h∇ζ) , (12)
is a second order nonlinear operator with

Q 2 [h, b](ζ) = -h(∇ ⊥ h • ∇)∇ ⊥ ζ - 1 2h ∇(h 2 ∇b • ∇ζ) + h 2 ∆ζ -(∇b • ∇ζ) ∇b, (13)
and for all smooth enough R-valued function w

Q 3 [h, h b ]w = 1 6 ∇(h 2 -h 2 b ) • ∇w + h 2 -h 2 b 3 ∆w - 1 6 ∆(h 2 -h 2 b )w. ( 14
)
Remark 4. We actually show in [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] that it is indeed possible to replace the inversion of

I + αT[h, b] by the inversion of 1 + αT[h b ],
where T[h b ] depends only on the fluid at rest (i.e. ζ = 0), while keeping the asymptotic O(µ 2 ) order of the expansion. The interest of working with (10) rather than ( 8) is that T[h b ] has a simplified scalar structure, i.e. it can be written in matricial form as

-1 3 ∇ • (h 3 b ∇ 1 h b •) 0 0 -1 3 ∇ • (h 3 b ∇ 1 h b •) . ( 15
)
From a numerical viewpoint, this simplified analytical structure allows to compute each component of the discharge hv separately, and to alleviate the computational cost associated with the dispersive correction of the model, as the discrete version of 1 + αT[h b ] may be assembled and factorized once and for all, in a preprocessing step.

Remark 5. The other difference with [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. I: Model development and analysis[END_REF] is the presence of the modified quadratic term Q[h, b], which shares with (8) the nice property that no computation of third order derivative is needed. The price to pay

is the inversion of an extra linear system, through the computation of

Q 3 [h, h b ] 1 + αT[h b ] -1 (h∇ζ) .
However, this extra cost is largely off-set by the gain obtained by using the time independent scalar operator

T[h b ]. Remark 6. One could replace Q 3 [h, h b ] 1 + αT[h b ] -1 (h∇ζ) by Q 3 [h, h b ] (h∇ζ)
in the second equation of [START_REF] Favrie | A robust and rapid numerical method for dispersive models admitting a lagrangian : application to Serre-Green-Naghdi equations for long free surface gravity waves[END_REF], keeping the same asymptotic O(µ 2 ) order. This would avoid the resolution of this extra linear system but leads to strong instabilities. We refer to [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] for more comments on this point.

Remark 7. As shown in [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF], the model (10) can be recast as follows:

         ∂ t h + ∇ • (hv) = 0, ∂ t (hv) + ∇ • (hv ⊗ v) + gh∇ζ + D c = 0, 1+αT[h b ] (D c + 1 α gh∇ζ) = h( 1 α g∇ζ + Q 1 [h, b](v) + gQ 2 [h, b](ζ)) + Q 3 [h, h b ]K, 1+αT[h b ] K = gh∇ζ, (16)
and we see that the dispersive correction D c acting as a source term in the second equation, is obtained as the solution of auxiliary scalar second-order elliptic sub-problems. This formulation will be referred to as (CG α ) in the following. Remark 8. We highlight that (CG α ) still does not involve third order derivative computation, which is a very interesting property. It is indeed shown in [START_REF] Filippini | A flexible genuinely nonlinear approach for wave propagation, breaking and runup[END_REF] that when third order derivatives on the free surface occur, it is important to introduce some sophisticated approximation strategies for the free surface gradient to reduce the dispersion error.

Remark 9. Considering the case of flat bottoms and d = 1 for the sake of simplicity, the linear dispersion relation associated with [START_REF] Dong | A reconstructed central discontinuous galerkin-finite element method for the fully nonlinear weakly dispersive Green-Naghdi model[END_REF] is classically obtained through the study of the linear behavior of small perturbation ( ζ, v) to a constant state solution (ζ, v) of ( 16) and looking for plane wave solutions of the form (ζ 0 , v 0 )e i(k•x-ωt) to the linearized system. Skipping computing details, one finds the following dispersion relation

1 gh|k| 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 , ( 17
)
where h = h 0 + ζ. Note that the dispersion relation obtained for the non-enhanced model is obtained choosing α = 1. For the applications considered here, we are interested in obtaining a model with the widest possible range of validity: in addition to the shallow water zone, we aim at providing a correct description of the waves in intermediate to moderately deep water, that is for increased values of |k|h 0 . The discussion concerning the choice of α to improve the dispersive properties of the model follows the usual procedure (see for instance [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. II: Boundary conditions and validation[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]) and we take α = 1.159 throughout this article. Such enhancement allows to achieve an improved matching between the phase and group velocities associated with [START_REF] Antunes Do Carmo | Numerical solution of the generalized Serre equations with the Mac-Cormack finite-difference scheme[END_REF] and the reference phase and group velocities coming from Stokes linear theory up to |k|h 0 = 4.

A pre-balanced Green-Naghdi formulation

Paving the way towards the construction of a well-balanced and efficient discrete formulation in §3, we now adjust the pre-balanced approach of [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF] for the case d = 1 and [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF] for the case d = 2 to the Green-Naghdi equations (CG α ). We introduce the total free surface elevation η = h + b, denote W = t (η, q) and use the following splitting of the hydrostatic pressure term:

gh∇ζ = 1 2 g∇(η 2 -2ηb) + gη∇b, (18)
to obtain a one-parameter pre-balanced constant-diagonal Green-Naghdi model (more simply referred to as (PCG α ) model in the following):

             ∂ t η + ∇ • q = 0, ∂ t q + ∇ • q ⊗ q η -b + 1 2 g∇(η 2 -2ηb)I + D c = -gη∇b, 1+αT[h b ] (D c + 1 α gh∇η) = h( 1 α g∇η + Q 1 [h, b](v) + gQ 2 [h, b](η)) + Q 3 [h, h b ]K, 1+αT[h b ] K = gh∇η. ( 19
)
Remark 10. The use of the pre-balanced formulation allows to straightforwardly construct a discrete formulation that exactly preserves motionless steady states, as detailed in §3.2.2 and Proposition 17. It also helps to slightly reduce the number of quadrature nodes needed to exactly achieve the preservation of motionless steady states, see Remark 20 for details.

3 Discrete formulation

Settings and notations

Let Ω ⊂ R d , d = 2, denote an open bounded connected polygonal domain with boundary ∂Ω. We consider a geometrically conforming mesh T h defined as a finite collection of |T h | nonempty disjoint open triangular elements T of boundary ∂T such that Ω = T ∈T h T . The meshsize is defined as h = max

T ∈T h h T
with h T standing for the diameter of the element T and we denote |T | the area of T , p T its perimeter and n T its unit outward normal. Mesh faces are collected in the set F h and the length of a face F ∈ F h is denoted by |F |. We denote by |F h | the number of faces. A mesh face F is such that either there exist

T 1 , T 2 ∈ T h such that F ⊂ ∂T 1 ∩ ∂T 2 (F is called an interface and F ∈ F i h ) or there exists T ∈ T h such that F ⊂ ∂T ∩ ∂Ω (F is called a boundary face and F ∈ F b h ). For all T ∈ T h , F T := {F ∈ F h | F ⊂ ∂T } denotes the set of faces belonging to ∂T and, for all F ∈ F T , n T F is the unit normal to F pointing out of T .
In what follows, we consider P k (T h ) the broken bivariate polynomial space defined as follows:

P k (T h ) := {v ∈ L 2 (Ω) | v |T ∈ P k (T ) ∀ T ∈ T h }, (20)
where P k (T ) denotes the space of bivariates polynomials in T of degree at most k, and we define

X h = P k (T h ) × P k (T h ) 2 . We denote N k = dim(P k (T )) = (k + 1)(k + 2)/2.
To discretize in time, for a given final computational time t max , we consider a partition (t n ) 0≤n≤N of the time interval [0, t max ] with t 0 = 0, t N = t max and t n+1 -t n = ∆t n . For any sufficiently regular function of time w, we denote by w n its value at discrete time t n . More details on the computations of ∆t n and the time marching algorithms are given in §3.4.1.

Discrete formulation for the advection-dominated equations

Postponing to §3.3 the computation of the projection of the dispersive correction D(W, b) on the approximation space, we focus here on the discrete formulation associated with the 2 first equations of [START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF] written in a more compact way:

∂ t W + ∇ • F(W, b) + D(W, b) = B(W, b), (21)
with

B(W, b) = t 0, -gη t ∇b , D(W, b) = t 0, t D c , and
F(W, b) =   t q q ⊗ q η -b + 1 2 g∇(η 2 -2ηb)I 2   . (22)

The discrete problem

We seek an approximate solution W h = t (η h , t q h ) of ( 21) in X h . Requiring the associated residual to be orthogonal to P k (T h ), the semi-discrete formulation reduces to the local statement: find (η h , t q h ) in

X h such that T d dt W h π h dx - T F(W h , b h ) • ∇π h dx + ∂T (F(W h , b h ) • n T ) π h ds + T D(W h , b h )π h dx = T B(W h , b h )π h dx, (23)
for all π h ∈ P k (T h ) and all element T ∈ T h , where b h refers to the L 2 projection of b on P k (T h ) and

D(W h , b h ) stands for a polynomial description of D(W, b) in X h , to be obtained in §3.3.
Remark 11. Considering a local basis {φ i } N k i=1 for a given element T ∈ T h , the local discrete solution W h|T may be expanded as:

W h|T (x, t) = N k i=1 Wi (t)φ i (x), ∀x ∈ T, ∀t ∈ [0, t max ], (24)
where { Wi (t)} N k i=1 are the local expansion coefficients, defined as Wi (t) = t (η i , t qi ) and qi = t ((q x ) i , (q y ) i ). Many choices are of course possible for the basis. In the following, {φ i } N k i=1 refers to the interpolant (nodal) basis on the element T and we choose the Fekete nodes [START_REF] Taylor | An algorithm for computing Fekete points in the triangle[END_REF] as approximation points.

Equipped with such local basis, following Remark 11, and splitting the boundary integral into faces integrals, the local statement ( 23) is now equivalent to

N k i=1 T φ i φ j dx d dt Wi (t) - T F(W h , b h ) • ∇φ j dx + F ∈FT F F T F φ j ds + T D(W h , b h )φ j dx = T B(W h , b h )φ j dx, 1 ≤ j ≤ N k , (25)
where F T F is a stabilizing numerical approximation of the normal interface flux F(W h , b h ) • n T F , to be defined in the following.

Interface fluxes and well-balancing

We recall in the following a simple choice to approximate the interface fluxes [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF], leading to a wellbalanced scheme that preserves motionless steady states. This modified flux can also be seen as the adaptation of the ideas of [START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF] to the pre-balanced formulation [START_REF] Pearce | A pseudo-spectral algorithm and test cases for the numerical solution of the two-dimensional rotating Green-Naghdi shallow water equations[END_REF]. Consider a face F ∈ F T (for the sake of simplicity, we only focus on the case F ∈ F i h and do not detail the case F ∈ F b h ). Let us denote denote W -and W + respectively the interior and exterior traces on F , with respect to the elements T . Similarly, b -and b + stand for the interior and exterior traces of b h on F . We define:

b * = max(b -, b + ), b = b * -max(0, b * -η -) (26)
and

ȟ-= max(0, η --b * ), ȟ+ = max(0, η + -b * ), (27) η
-= ȟ-+ b, η+ = ȟ+ + b, (28)
leading to the new interior and exterior values:

W-= t (η -, ȟ- η --b -q -), W+ = t (η + , ȟ+ η + -b + q + ). ( 29
)
Now we set

F T F = F h ( W-, W+ , b, b, n T F ) + F T F , (30)
as the numerical flux function through the interface F , where:

1. the numerical flux function F h is the global Lax-Friedrichs flux:

F h (W -, W + , b -, b + , n T F ) = 1 2 ((F(W -, b -) + F(W + , b + )) • n T F -a(W + -W -)), (31)
with a = max

T ∈T h λ T and λ T = max ∂T q h|T η h|T -b h|T • n T + g(η h|T -b h|T ) . (32)
2. F T F is a correction term defined as follows:

F T F =   0 0 g η-( b -b -) 0 0 g η-( b -b -)   • n T F . (33)
Note that the modified interface flux [START_REF] Filippini | A flexible genuinely nonlinear approach for wave propagation, breaking and runup[END_REF] only induces perturbations of order k + 1 when compared to the traditional interface fluxes.

Remark 12. Let denote in the following by w T the averaged value of the discrete approximation w h|T on the element T , for any R or R 2 -valued function w. Let consider the corresponding first order FV scheme for the averaged free-surface:

η n+1 T = η n T - ∆t n |T | F ∈FT F F h ( W-, W+ , b, b, n T F ) ds, (34)
with F h defined following [START_REF] Engsig-Karup | Nodal DG-FEM solution of high-order Boussinesq-type equations[END_REF] and W-, W+ , b obtained from ( 26)-( 29), starting from W -and W + defined respectively, for each face F , as the first-order piecewise constant interior and exterior values W n T1 and W n T2 , with

T F = {T 1 , T 2 }. Then, assuming that h n T ≥ 0, ∀T ∈ T h , we have h n+1 T ≥ 0 ∀T ∈ T h , under the condition λ T p T |T | ∆t n ≤ 1, ∀T ∈ T h . ( 35
)
This is a positivity property of the FV scheme based on Lax-Friedrichs flux within the pre-balanced formulation, which is shown in [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF] and is mandatory to obtain a positive high-order DG scheme, see [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF].

Positivity

The enforcement strategy of the water height non-negativity within the DG formulation [START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF] requires positivity of the water height at carefully chosen quadrature nodes at the beginning of each time step. Then, the positivity of the water height is ensured providing the use of the LF flux [START_REF] Engsig-Karup | Nodal DG-FEM solution of high-order Boussinesq-type equations[END_REF] and an additional time-step restriction, see Proposition 17.

We adjust these ideas to the (PCG α ) formulation [START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF] to enforce the mandatory property that the mean value of the water height h T = η T -b T on any given element T ∈ T h remains positive during the time marching procedure. The main ideas are summarized for an explicit first order Euler scheme in time for the sake of simplicity:

♯ 1 considering the broken polynomial space P k (T ), we assume that the face integrals are computed using (k + 1)-points Gauss quadrature (see Remark 20). The special quadrature rules introduced in [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF] are obtained by a transformation of the tensor product of a β-points Gauss-Lobatto quadrature (with β the smallest integer such that 2β -3 ≥ k) and the (k + 1)-point Gauss quadrature. This new quadrature includes all (k + 1)-point Gauss quadrature nodes for each face F ∈ F T , involves positive weights and it is exact for the integration of η n h over T . In the following, let us denote S k T the set of points of this quadrature rule. We show on Fig. 2 the resulting quadrature nodes used for k = 2 and k = 3 orders of approximation on a reference element.

♯ 2 for each element T ∈ T h , h n h|T is computed from η n h|T and b h|T and we need to ensure that h n h|T (x) ≥ 0, ∀x ∈ S k T , which is a sufficient condition to ensure the non-negativity property for the DG scheme [START_REF] Aissiouene | A combined finite volume -finite element scheme for a dispersive shallow water system, Networks and heterogeneous media[END_REF], under the CFL-like condition [START_REF] Zhang | Error estimates to smooth solutions of runge-kutta discontinuous galerkin methods for scalar conservation laws[END_REF]. This condition is enforced using the accuracy preserving limiter of [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF]. Assuming h n T ≥ 0, we replace h n h|T by a conservative linear scaling around this element average:

ĥn h|T = θ T (h n h|T -h n T ) + h n T , (36)
with

θ T = min h n T h n T -m T
, 1 , and m T = min Remark 13. Note that this approach also ensures that the water height remains positive at the k + 1 Gauss quadrature nodes used to compute the faces integrals in (25).

x∈S k T h n h|T (x).

Discrete formulations for the second order elliptic sub-problems

We are now left with the computation of the projection D(W h , b h ) of the dispersive correction on the approximation space X h . As its first component identically vanishes, the computation of the discrete version of D c (denoted by D c,h ∈ P k (T h ) 2 in the following) is obtained as the solution of the discrete problems associated with

1+αT[h b ] K = gh∇η, (37) 1+αT
[h b ] (D c + 1 α gh∇η) = h( 1 α g∇η + Q 1 [h, b](v) + gQ 2 [h, b](η)) + Q 3 [h, h b ]K. ( 38
)
Note that both problems ultimately reduce to the construction of a discrete formulation associated with the generic scalar problem 1+αT

[h b ] w = f,
where the source term f is successively defined as the directional scalar components of the source terms occurring in (37)- [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF], that is to say

f =    gh∇η for (37) 1 α gh∇η + h(Q 1 (v) + gQ 2 (η)) + Q 3 [h, h b ]K for (38), (39)
where we keep a vectorial notation for the sake of simplicity. Remarking now that the following identity holds for any smooth enough scalar-valued function:

T[h b ]w = - 1 3 ∇ • H b ∇w + 1 6 ∇ • w∇H b , (40)
with the notation H b = h 2 b , we consider the following mixed formulation, introducing a diffusive flux p:

w -α∇ • 1 3 H b p -1 6 w∇H b = f, p -∇w = 0, ( 41
)
and the following associated discrete problem: find (w h , t p h ) in X h such that

T w h π h dx + α 3 T H b h (p h • ∇π h )dx - α 3 ∂T H b h (p h • n T )π h ds (42) - α 6 T w h (∇H b h • ∇π h ) dx + α 6 ∂T w h (∇H b h • n T )π h ds = T f h π h dx, T p h • φ h dx + T w h ∇ • φ h dx - ∂T w h (φ h • n T ) ds = 0, (43)
for all (π h , φ h ) ∈ X h , where H b h and ∇H b h respectively denote the L 2 projection of h 2 b and ∇h 2 b on P k (T h ) and (P k (T h )) 2 , and f h is a polynomial representation of f in P k (T h ), with f defined according to [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]. We choose to use the Local Discontinuous Galerkin (LDG) approach [START_REF] Cockburn | The Local Discontinuous Galerkin method for time-dependent convection-diffusion systems[END_REF] to compute the associated stabilizing fluxes. For any given element T ∈ T h , we have:

∂T H b h (p h • n T )π h ds = F ∈FT F H b h ( p T F • n T F )π h ds, (44)
∂T

w h (φ h • n T ) ds = F ∈FT F w T F (φ h • n T F ) ds, (45)
∂T

w h (∇H b h • n T )π h ds = F ∈FT F w T F (∇H b h • n T F )π h ds, (46)
where the interface fluxes are computed as follows:

w T F = { {w h } } + (β • n T F ) , a + b w h , (47)
p T F = { {p h } } -β( p h • n T F ) + ξ |F | w h n T F , (48)
with the classical notations for the face average { {w} } = (w + + w -)/2 and face jump w = w + -w -, and w -and w + stand respectively for the interior and exterior traces of w h with respect to the considered face F (and similar notations for the R 2 -valued quantity p). The penalization parameter ξ and the upwinding parameter β are both set to 1.

Remark 14. The choice of the LDG fluxes ( 47)-( 48) allows to eliminate locally the auxiliary discrete flux p h , and to globally assemble the matrix corresponding to the discrete formulation of

1 + αT[h b ] .
The auxiliary flux does not appear in the final formulation. Note also that, for the sake of simplicity, we also use the upwind flux [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] to approximate the faces contributions associated with the purely advective part of [START_REF] Taylor | An algorithm for computing Fekete points in the triangle[END_REF].

Remark 15. The relative simplicity of the discrete formulation ( 42)-( 43) is a direct consequence of the simplified analytical structure of the (PCG α ) formulation. One can of course adapt this approach to the initial formulation [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. I: Model development and analysis[END_REF], leading in practice to a more costly and complicated global assembling process, as the operator T [h, b] directly acts on a R 2 -valued function.

For the sake of efficiency, the computations of the source terms integrals T f h π h dx occurring in [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF] are performed in a collocation way, in the spirit of pseudo-spectral methods. More precisely, we build approximated integrands with polynomial representations (Q j,h ) j=1,2,3 in P k (T h ) 2 of the operators (Q j ) j=1,2,3 using direct products of the discrete flow variables (η h , v h ) and their derivatives at the Fekete nodes. These derivatives are however weakly computed with the DG approach, using LDG fluxes. In practice, this is simply achieved by adjusting the approach of [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF] to the d = 2 case to easily compute the required first and second order derivatives of η h and v h in each space direction. Again, the second order spatial derivatives are written in mixed form and we use the stabilizing fluxes ( 47)- [START_REF] Tonelli | Hybrid finite volume-finite difference scheme for 2dh improved Boussinesq equations[END_REF]. The corresponding discrete formulations are simplified through local elimination of the diffusive fluxes, leading to the assembling of global matrices for first and second order weak derivatives in each direction. For instance, considering derivatives in the first direction, the mixed form

v + ∂ x w = 0 , u + ∂ x v = 0,
leads to the following local discrete formulations for 1 ≤ j ≤ N k :

N k i=1 ṽi M ij = N k i=1 wi S x ij - F ∈FT F w T F φ j n x T F ds , N k i=1 ũi M ij = N k i=1 ṽi S x ij - F ∈FT F v T F φ j n x T F ds, (49)
where n x T F is the first component of n T F and with

M ij := T φ i φ j dx , S x ij := T φ i ∂ x φ j dx.
These systems can be globally rewritten as follows:

M Ṽ = S x W -(E x -F x ) W , M Ũ = S x Ṽ -(E x + F x ) Ṽ - ξ h F x W , (50)
where: (i) Ũ , Ṽ and W are N k × card(T h ) vectors gathering the expansion coefficients ũ, ṽ and w for all mesh elements, (ii) M and S x are the square card(T h ) × N k global mass and stiffness matrices, with a block-diagonal structure, (iii) E x and F x are square card(T h ) × N k matrices, globally assembled by gathering all the mesh faces contributions, accounting respectively for the average { {•} } and jump • operators occurring in the definition ( 47)-( 48) of the LDG fluxes. Note that each interface F ∈ F i h contributes to four blocks of size N k in the global matrices and each boundary face F ∈ F b h contributes to one block. Equipped with these global structures, first and second order global differentiation matrices can be straightforwardly assembled, leading to the following compact notations:

Ṽ = D x W , with D x = M -1 (S x -E x + F x ) , (51)
Ũ = D xx W , with D xx = M -1 (S x -E x -F x ) D x - ξ h F x . (52)
Similar constructions are performed for D y and D yy .

Remark 16. The use of such a direct nodal products method helps to reduce the computational cost by reducing the accuracy of quadrature but may threaten the numerical stability for strongly nonlinear or marginally resolved problems with the possible introduction of what is known in the field of spectral methods as aliasing driven instabilities. Several well-known methods may help to alleviate this issue, like the use of stabilization filtering methods. In this work however, and considering the test cases studied in §4, we did not need to use any additional stabilization mechanism.

3.4 Time-marching, boundary conditions and wave-breaking

Time discretization

The time stepping is carried out using the explicit SSP-RK scheme [START_REF] Gottlieb | Strong stability preserving high order time discretization methods[END_REF]. Up to k = 3, we consider RK-SSP schemes of order k + 1. A fourth order SSP-RK scheme is used for k ≥ 3. For instance, writing the semi-discrete equations as d dt W h + A h (W h ) = 0, advancing from time level n to n + 1 is computed as follows with the third-order scheme:

       W n,1 h = W n h -∆t n Ãh (W n h ) , W n,2 h = 1 4 (3W n h + W n,1 h ) -1 4 ∆t n Ãh (W n,1 h ) , W n+1 h = 1 3 (W n h + 2w n,2 h ) -2 3 ∆t n Ãh (W n,2 h ) . ( 53
)
with the formal notation Ãh = A h • ΛΠ h , W h ← ΛΠ h W h being the limitation operator possibly acting on the approximated vector solution (see §3.4.3), and ∆t n is obtained from the CFL condition (57).

Boundary conditions

For the test cases studied in the next section, boundary conditions are imposed weakly, by enforcing suitable reflecting relations at virtual exterior nodes, at each boundaries, allowing to compute the corresponding interface stabilizing fluxes. Solid-wall (reflective) but also periodic conditions (as the computational domains geometries of the cases studied in §4 are rectangular) can be enforced following this simple process. These simple boundary conditions possibly have to be complemented with ad-hoc absorbing boundary conditions, allowing the dissipation of the incoming waves energy together with an efficient damping of possibly non-physical reflections, and generating boundary conditions that mimic a wave generator of free surface waves. We have implemented relaxation techniques and we enforce periodic waves combined with generation/absorption by mean of a generation/relaxation zone, following the ideas of [START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF], using the relaxation functions described in [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF], and the computational domain is locally extended to include sponge layers which may also include a generating layer. More explicitely, the relaxed solution along the domain will take the form :

W relax = F a W h + (1 -F a )F g R(t)W imp , (54)
where F a , F g stand for the absorption and generation profiles and R governs the time evolution of the generation process. Above, W imp contains the imposed profile. Concerning the relaxation functions, we follow [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF], taking the exponential forms :

F a (x) = 1 - exp ((x r ) n -1)) exp(1) -1 , F g (x) = 1 - exp ((1 -x r ) n -1)) exp(1) -1 , where x r = x -x R ∆ R
, n is a positive parameter, and x R , ∆ R are respectively the beginning and the width of the relaxation zone. In agreement with other works, the length of the sponge layers ∆ R is calibrated for each test case (generally 2 or 3 wavelengths); the parameter n is fixed to 3. The reader is referred to the above references for more details.

Wave-breaking and limiting

Obviously, vertically averaged models cannot reproduce the surface wave overturning and are therefore inherently unable to fully model wave breaking. Moreover, if the GN equations can accurately reproduce most phenomena exhibited by non-breaking waves in finite depth, including the steepening process occurring just before breaking, they do not intrinsically account for the energy dissipation mechanism associated with the conversion into turbulent kinetic energy observed during broken waves propagation. Several methods have been proposed to embed wave breaking in depth averaged models. Many of them focus on the inclusion of an energy dissipation mechanism through the activation of extra terms in the governing equations when wave breaking is likely to occur, and the reader is referred to [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] for a recent review of these approaches. More recently, hybrid strategies have been elaborated for weakly nonlinear models [START_REF] Tonelli | Hybrid finite volume-finite difference scheme for 2dh improved Boussinesq equations[END_REF][START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF] and for fully nonlinear models [START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF][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][START_REF] Tissier | Nearshore dynamics of tsunami-like undular bores using a fully nonlinear Boussinesq model[END_REF]. Roughly speaking, the idea is to switch from GN to NSW equations when the wave is ready to break by locally suppressing the dispersive correction and the various approaches may differ by the level of sophistication of the detection criteria and the switching strategies. Denoting that, from a numerical point of view, such an approach allows to avoid numerical instabilities by turning off the computation of higher-order derivatives and nonlinear and non-conservative terms in the vicinity of appearing singularities, we use in [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF] a purely numerical smoothness detector to identify the potential instability areas, switch to NSW equations in such areas by simply locally neglecting the D c source term and use a limiter strategy to stabilize the computation, letting breaking fronts propagate as moving bores.

As the aim of the present work is only to introduce and validate our fully discontinuous discrete formulation in the multidimensional case, we do not focus on the development of new wave breaking strategies and we simply adjust this simple approach to the d = 2 case to possibly stabilize the computations performed in §4.

We detect the troubled elements (following the terminology of [START_REF] Qiu | A comparison of troubled-cell indicators for Runge-Kutta discontinuous Galerkin methods using weighted essentially nonoscillatory limiters[END_REF]) using the criterion proposed in [START_REF] Krivodonova | Shock detection and limiting with discontinuous galerkin methods for hyperbolic conservation laws[END_REF], and based on the strong superconvergence phenomena exhibited at element's outflow boundaries. More precisely, for any T ∈ T h , denoting by ∂T in the inflow part of ∂T , we use the following criterion:

I T = F ∈∂Tin F (h --h + ) ds h (k+1)/2 T |∂T in | h h|T , (55)
which has already successfully been used as a troubled cells detector in purely hyperbolic shallow water models, see among others [START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF][START_REF] Kesserwani | Locally limited and fully conserved RKDG2 shallow water solutions with wetting and drying[END_REF][START_REF] Zhu | Runge-Kutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshes[END_REF].

If I T ≥ 1 then we apply a slope limiter on each scalar component of W h|T , based on the maxmod function, see [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF]. This limiting strategy is not recalled here, as we straightforwardly reproduce the implementation described in [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF] for the NSW equations in pre-balanced form. Of course, waves about to break do not embed any free surface singularities yet, but our numerical investigations have shown that when the wave has steepened enough, the free surface gradient becomes large enough to activate the criteria. This simple approach, although quite rough and far less sophisticated than recent strategies introduced for instance in [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF][START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF], allows us to obtain good results in the various cases of §4.

Main properties

We have the following result:

Proposition 17. The discrete formulation [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] together with the interface fluxes discretization (30) and a first order Euler time-marching algorithm has the following properties:

1. it preserves the motionless steady states, providing that the integrals of (25) are exactly computed for the motionless steady states. In other terms, we have for all n ∈ N:

η n h ≡ η e q n h ≡ 0 ⇒ η n+1 h ≡ η e q n+1 h ≡ 0 , (56)
with η e constant,

assuming moreover that h

n T ≥ 0, ∀T ∈ T h and h n h|T (x) ≥ 0, ∀x ∈ S k T , ∀T ∈ T h , then we have h n+1 T ≥ 0, ∀T ∈ T h under the condition λ T p T |T | ∆t n ≤ 2 3 ωβ 1 , (57)
where ωβ 1 is the first quadrature weight of the β-point Gauss-Lobatto rule used in §3.2.3.

Proof.

1. assuming that the following equilibriumW h = W e h = (η e , 0) holds, we have to show that ∀T ∈ T h and 1 ≤ j ≤ N k :

T F(W e h , b h ) • ∇φ j dx - F ∈FT F F e T F φ j ds - T D(W e h , b h )φ j dx + T B(W e h , b h )φ j dx = 0, ( 58
)
where F e T F is the interface numerical flux obtained at equilibrium. Looking at [START_REF] Filippini | A flexible genuinely nonlinear approach for wave propagation, breaking and runup[END_REF], and highlighting that for each interface F we have η-= η+ = η e and therefore W-= W+ , it is easy to check that F e T F = F(W -, b -)• n T F , thanks to the consistency of the numerical flux function F h . Consequently, we have

T F(W e h , b h ) • ∇φ j dx - F ∈FT F F e T F φ j ds = - T ∇ • F(W e h , b h )φ j dx, (59)
= - T B(W e h , b h )φ j dx. ( 60
)
assuming that the integrals are computed exactly (see Remark 20), and observing that we have

∇ h • F(W e h , b h ) = B(W e h , b h ).
For the integrals associated with the dispersive source term, it is straightforward to check from the definitions (4) and ( 13) of (Q j ) j=1,2 that Q 1,h = 0 whenever q h = 0 and that Q 2,h = 0 whenever η h = η e . Moreover, using again the fact that ∇ h η e = 0, the polynomial representation K h of K, obtained as the solution of the discrete problem ( 42)-( 43) associated with [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF], taking f h = 0 as source term and suitable boundary conditions, identically vanishes leading to Q 3,h = 0 using [START_REF] Panda | Discontinuous Galerkin methods for solving Boussinesq-Green-Naghdi equations in resolving non-linear and dispersive surface water waves[END_REF]. Note that this result holds no matter what method we use to compute polynomials representations (Q j,h ) j=1,2,3 . [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF], adjusted to the pre-balanced framework in [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF], can be straightforwardly reproduced, as the dispersive source term D(W h , b h ) has no influence on the mass conservation equation. Note that this is the condition (57) (more restricting than [START_REF] Castillo | An a priori error analysis of the Local Discontinuous Galerkin method for elliptic problems[END_REF]) which is used in practice. Condition [START_REF] Castillo | An a priori error analysis of the Local Discontinuous Galerkin method for elliptic problems[END_REF] only ensures the positivity of a FV scheme based on the LF flux.

the proof of

Remark 18. The previous analysis can be extended to the third order SSP scheme of §3.4.1 (or any TVD high order Runge-Kutta and multi-step time discretizations) at the cost or more technicality, see for instance [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF] and [START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF].

Computational cost

In order to provide some insights concerning the computational cost per iteration of the resulting discrete formulation, let us first consider the advective part [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]. Neglecting the topography source term for the sake of simplicity and temporarily postponing the analysis of the computation of the dispersive term, we observe that the discrete solution may be updated at each time step (or sub-step, see §3.4.1) providing the computation of the following residual :

- T F(W h , b h ) • ∇φ j dx + F ∈FT F F T F φ j ds,
for each mesh element T and basis functions φ j . To compute the surface integrals, we have to evaluate products of the form ω q (∂ m φ j ) q (F lm ) q , at each single quadrature point (in physical coordinates) indexed by q, where ω q is the corresponding quadrature weight (involving a factor coming from the determinant of the transformation's jacobian), ∂ k φ j is the partial derivative of φ j in the k th direction, (F lm ) 1≤l≤d+1, 1≤m≤d refer to the components of F(W h , b h ) regarded as a matrix and the notation ( • ) q is used for the evaluation of the quantities of interest at the quadrature node indexed by q.

As the jacobians and the derivatives values at quadrature nodes are pre-computed, it leads to 2d(d + 1)N k + C(F) floating point operations (flop in the following) at each quadrature node, where C(F) refers to the evaluation of the components (F lm ) g , and therefore [2d(d + 1)N k + C(F)]|T h |N q flop for the evaluation on Ω, with N q the number of quadrature nodes per element.

Assuming that the jumps of basis functions on each face are also pre-computed, similar considerations for the line integrals lead to [2(d+ 1)(k + 1)+ C( F)]|F h |N g flop, with N g the number of Gauss quadrature nodes per face and C( F) the number of flop needed to evaluate the numerical flux at a given node. Note that using the nodal approach and the basis functions of Lagrange, there are only k + 1 basis functions with a nonzero support on a given face. We also have to account for the (precomputed) inverse of the mass matrix products, leading to some additional [3

(N k ) 2 + 3N k (k + 1)]|T h | flop.
Remark 19. To precisely evaluate the total number of flop of the advective parts, one would have to meticulously account for the state reconstructions of §3.2.2, and the limiting procedures of §3.2.3 and §3.4.3, which affect the values of C(F) and C( F).

Remark 20. The line integrals of (25) can be exactly computed at motionless steady states with N g = k + 1 point Gauss quadrature nodes, thanks to the pre-balanced splitting [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF], allowing to reduce the number of numerical fluxes evaluation per face. Indeed, with the resulting formulation of the advective flux [START_REF] Choi | An efficient curvilinear non-hydrostatic model for simulating surface water waves[END_REF], the integrands on faces are of order 2k at equilibrium, while they are of order 3k within the formulation (16)). In the same way, the quadrature rules used for the surface integrals have to be exact only for bivariate polynomials of order 2k -1. Higher order quadrature rules would be needed to exactly compute these integrals and ensure exact well-balancing for the formulation [START_REF] Dong | A reconstructed central discontinuous galerkin-finite element method for the fully nonlinear weakly dispersive Green-Naghdi model[END_REF].

Concerning now the dispersive part, the most computationally demanding step concerns the resolution of the 4 linear systems of size N dof = N k |T h | associated with ( 37)- [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF]. Recalling the time-independent nature of the second order operator 1+αT[h b ] , the assembling of the corresponding matrix detailed in §3.3 is performed in a pre-processing step, allowing to exclude the corresponding number of flop from the computational cost per iteration. Additionally, we also pre-compute the LU factorization of this (unsymmetric) matrix, saving the associated Θ(N 3 dof ) flop. Finally, we are left with the resolution of the resulting triangular systems through forward and backward substitutions, with an overall Θ(N 2 dof ) asymptotic cost, and the computation of the corresponding second members [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF], which is straightforwardly performed by collocation at the N dof interpolation nodes.

Numerical validations

In this section, we validate the previous discrete formulation through several benchmarks. Unless stated otherwise, we use periodic boundary conditions in each direction, we set α = 1.159, ε 0 = 0.1 and the time step restriction is computed according to [START_REF] Zhang | Error estimates to smooth solutions of runge-kutta discontinuous galerkin methods for scalar conservation laws[END_REF]. Accordingly with the non-linear stability result of the previous section, we do not suppress the dispersive effects in the vicinity of dry areas. Some accuracy analysis are performed in the first two cases using the broken L 2 norm defined as follows for any arbitrary scalar valued piecewise polynomial function w h defined on T h :

w h L 2 (T h ) =   T ∈T h w h|T 2 L 2 (T )   1 2 .

Preservation of motionless steady state

This preliminary test case is devoted to check the ability of the formulation to preserve motionless steady states and accuracy validation. The computational domain is the [-1,1] × [-1,1] square, and we use an unstructured mesh of 8466 elements. The bottom elevation involves a bump and a hollow having same dimensions, respectively located at

x 1 = (x 1 , y 1 ) = (- 1 3 , - 1 3
) and

x 2 = (x 2 , y 2 ) = ( 1 3 , 1 3
), leading to the following analytic profile :

b(r 1 , r 2 ) = 1 + d e -(r1/L) 2 -d e -(r2/L) 2 , (61)
where r 1,2 are respectively the distances from x 1 and x 2 and we set d = 0.45 and L = 0.15. The reference water depth is h 0 = 1.5 m, leading to the configuration depicted on Fig. 3. Our numerical investigations confirm that this initial condition is preserved up to the machine accuracy for any value of polynomial order k. For instance, the L 2 numerical errors obtained at t = 50 s using a k = 2 approximation are respectively 3.0e-16, 7.4e-16 and 7.9e-16 for η, hu and hv.

Keeping the same computational domain and topography profile, we also perform a convergence study, using a reference solution obtained with k = 4 and a regular triangulation with space steps ∆x = ∆y = 2 -9 m. The initial free surface is set to:

ζ(t = 0, x) = a e -( x /L) 2 , a = 0.075 h 0 . (62)
The computations are performed on a sequence of regular triangular meshes with increasing refinement ranging from 2 -4 m to 2 -9 m and polynomial expansions of degrees ranging from k = 1 to k = 4, while keeping the time step constant and small enough to ensure that the leading error orders are provided by the spatial discretization. The L 2 -errors computed at t max = 0.02 s lead to the convergence curves shown on Fig. 4 for the free surface elevation and Fig. 5 for the discharge. The corresponding convergence rates obtained by linear regression are reported on each curve. The initial lack of convergence for coarse meshes might be explained, at least partially, by the aliasing errors developed in the collocation approximation of the highly nonlinear terms and third order derivatives on the topography in the dispersive part. We note that while the L 2 -errors are typically larger on the discharge than on the free surface, we asymptotically reach some convergence rates between O(h k+ 1 2 ) and O(h k+1 ) for both variables. We observe that for this problem, the resulting convergence rates are greater than the O(h k+ 1

2 ) optimal estimate one would be expected, see [START_REF] Zhang | Error estimates to smooth solutions of runge-kutta discontinuous galerkin methods for scalar conservation laws[END_REF].

k = 1 k = 2 k = 3 k = 4
k = 1 k = 2 k = 3 k = 4

Solitary wave propagation

We consider now the time evolution of a solitary wave profile defined as follows:

   h(x, t) = h 0 + εh 0 sech 2 (κ(x -ct)) , u(x, t) = c 1 - h 0 h(x, t) , (63)
with κ = 3ε 4h2 0 (1 + ε)

, and c = gh 0 (1 + ε). Note that if such profiles are exact solutions of the original Green-Naghdi equations (2), or equivalently [START_REF] Antunes | Boussinesq and Serre type models with improved linear dispersion characteristics: Applications[END_REF] with α = 1, these are only solutions of the (CG α=1 ) model up to O(µ 2 ) terms. However, for small enough values of ε, such profiles are expected to propagate over flat bottoms without noticeable deformations. We use a rectangular computational domain of 200 m length and 25 m width, and a relatively coarse unstructured mesh with faces' length ranging from 0.7 m to 1.5 m. The reference water depth is set to h 0 = 1m, the relative amplitude is set to ε = 0.2 and the initial free surface, shown on Fig. 6, is centred at x 0 = 50 m. Cross sections along the x-direction centerline, obtained with k = 3 at several times during the propagation, are shown on Fig. 9. Note that the time step computed from ( 57) is around ∆t = 0.01 s, so that approximatively 3500 iterations are necessary to reach the t max = 35s of simulation.

To further investigate the convergence properties of our approach, we now consider a sequence of regular meshes with mesh size ranging from 2 -4 m to 2 -9 m and polynomial expansions of degrees ranging from k = 1 to k = 4. The corresponding L 2 -errors computed at t max = 0.2 s are used to plot the convergence curves on Fig. 7 for the free surface elevation and Fig. 8 for the discharge. The corresponding convergence rates obtained by linear regression are also reported and we observe an irregular behavior with rates varying between O(h k ) and O(h k+1 ), generally close to O(h To complete the picture, we also investigate the conservation of mass (denoted by M (t) in the following) and energy E(t) over time by computing the time evolution of their relative error, see Fig. 10. As expected, the mass is preserved up to the machine accuracy, independently from the polynomial order k.

Unlike the classical GN equations, the (CG α=1 ) formulation is not expected to preserve the total energy of the system, but only up to O(µ 2 ) terms. These discrepancies are expected to remain relatively small in the shallow water regime µ ≪ 1. This is indeed also numerically confirmed on Fig. 10, where we can also appreciate the significant benefits of the space order accuracy increasing. Note that similar results are observed with other mesh sizes.

Run-up of a solitary wave

We investigate now the ability of the scheme in handling dry areas with a test based on the experiments of Synolakis [START_REF] Synolakis | The runup of solitary waves[END_REF]. We study the propagation, shoaling, breaking and run-up of a solitary wave over a topography with constant slope s = 1/19.85. The reference water depth is set to h 0 = 1 m, and a solitary wave profile is considered as initial condition [START_REF] Berkhoff | Verification of numerical wave propagation models for simple harmonic linear water waves[END_REF], with a relative amplitude ε = 0.28. The simulation involves a 50m × 5m basin, regularly meshed with a space step ∆x = 0.25 m and k = 2. We show on Fig. 11 some cross sections of the solution, taken at various times during the propagation and compared with the experimental data. The wave breaking is identified approximately at t * = 17 s, with the normalized time t * = t(g/h 0 ) 1/2 , and occurs between gauges #2 and #3, which is in agreement with the experiment. The whole breaking process is well reproduced, as well as the subsequent run-up phenomena.

Propagation of highly dispersive waves

The dispersive properties of the numerical model are now assessed through the study of the propagation of periodic waves over a submerged bar, following the experiments of Dingemans, see [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measure-ments[END_REF]. The com- For both test #A and #B, we set k = 2 and use a regular triangulation obtained from rectangular elements ∆x = ∆y = 0.125 m. The characteristics of the flow are quite complex here, notably due to the high non-linearities induced by the topography. The propagating waves first shoal and steepen over the submerged bar, generating higher-harmonics. These harmonics are progressively released on the downward slope, until encountering deeper waters again. In the first test, the initial amplitude is not large enough to trig the breaking of the waves and the GN equations are resolved in the whole domain. Fig. 12 indicates the location of the five wave gauges used for the first test to study the time-evolution of the free surface deformations. The results are plotted on Fig. 13. We observe a very good agreement between analytical and experimental data for the first wave gages. Of course the quality of the results at the two last gages is not as good as those obtained with optimized GN models such as [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Galan | Fully nonlinear model for water wave propagation from deep to shallow waters[END_REF], higher-level GN model [START_REF] Zhao | Application of higher-level GN theory to some wave transformation problems[END_REF] or the higher-order Boussinesq models with 2 layers of [START_REF] Lynett | Nearshore wave modeling with high-order Boussinesq-type equations[END_REF]. In particular, the extension of the present DG approach to the 3-parameters model of [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] is left for future works.

The test #B corresponds to the experiments of [START_REF] Beji | Numerical simulation of nonlinear wave propagation over a bar[END_REF]. The first reference gage is located at x = 6 m, followed by a series of seven gages regularly spaced between x = 11 m and x = 17 m. Wave breaking is observed at the level of the flat part of the bump during the propagation and to stabilize the computation, the dispersive terms are turned-off in the vicinity of troubled cells, using the strategy described in §3.4.3.

Again, a comparison between the numerical results and the experimental data is performed and shown on Fig. 14.

Periodic waves propagation over an elliptic shoal

We study now the propagation of a train of monochromatic waves evolving over a varying topography, following the experiment of Berkhoff et al [START_REF] Berkhoff | Verification of numerical wave propagation models for simple harmonic linear water waves[END_REF].

z s (x, y) =    -0.3 + 0.5 1 - x r 3.75 2 - y r 5 2 if x r 4 2 + y r 3 2 ≤ 1 0 elsewhere (64)
The reference water depth is set to h 0 = 0.45 m, and the corresponding computational domain has dimensions [-10, 12] × [-10, 10] (in m), with an extension of 5 m at inlet and outlet boundaries for the generation of incident waves and their absorption. The periodic wave train has an amplitude of a = 0.0232 m and a period of T = 1 s. Solid wall boundary conditions are used at y = 10 m and y = -10 m. We use an unstructured mesh of 25390 elements which is refined in the region where the bottom variations are expected to have the greatest impact on the wave transformations (the smallest and largest mesh face's lengths are respectively 0.2 m and 0.5 m). The computations are performed with polynomial approximations of order k = 2.

The data issued from the experiment provides the normalized time-averaged of the free surface measured along several cross sections. In our numerical experiment we focus on the following ones:

allowing a good coverage of the computational domain. Time series of the free surface elevation are hence recorded along these sections, from t = 30 s to t = 50 s and post-processed by mean of the zero up-crossing method to isolate single waves and compute the mean wave elevation. The results, classically normalized by the incoming wave amplitude a, are shown in Fig. 16 and compared with the data. The obtained results are in good agreement with the data, and qualitatively comparable, for instance, with those obtained in [START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] with other numerical methods. Note however that in [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF], the uniform cartesian mesh has 86000 rectangular elements (∆x ≈ 0.065 m and ∆y = 0.1 m) and that a 5 th order WENO scheme is used. In [START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF], 1024 × 1024 rectangular elements are needed to achieve the same qualitative agreement due to the use of a 2 nd order scheme whereas in [START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF], the mesh size in the y direction varies from h y = 0.05 m to 0.1 m and h x = 2h y .

Solitary wave propagation over a 3d reef

The following test case is based on the experiments of [START_REF] Swigler | Laboratory study investigating the 3d turbulence and kinematic properties associated with a breaking solitary wave[END_REF] and allows to study some complex wave's interactions such as shoaling, refraction, reflection, diffraction, breaking and moving shoreline in a fully two-dimensional context. The experimental domain is a 48.8 m large and 26.5 m wide basin. The topography is a triangular shaped shelf with an island feature located at the offshore point of the shelf. The island is a cone of 6 m diameter and 0.45 m height is also placed on the apex, between x = 14 m and x = 20 m (see Fig. 18). During the experiments, free surface information was recorded via 9 wave gauges, and the velocity information was recorded with 3 Acoustic Doppler Velocimeters (ADVs). The complete set up can be found in [START_REF] Swigler | Laboratory study investigating the 3d turbulence and kinematic properties associated with a breaking solitary wave[END_REF] but we also specify the gauge locations in the legend of Fig. 19, with exact values given on Fig. 22. Such a geometry motivates the use of an unstructured mesh of 45000 elements refined in the vicinity of the cone, leading to a smallest face's length of 0.12 m and a largest of 0.51 m, see Fig. 17. We set the polynomial order to k = 2. The reference water depth is h = 0.78 m and we study the propagation and transformations of a solitary wave of relative amplitude ε = 0.5. In addition to the complex processes already stated above, the transformations are quite nonlinear, making the following test particularly interesting. We show on Fig. 20-21 some snapshots of the free surface evolution before and after crossing the shelf's apex. We observe that the wave breaks at the apex slightly before t = 5 s, wrapping the cone around t = 6.5 s. We also observe refracted waves from the reef slopes and diffracted waves around the cone which converge at the rear side at approximatively t = 8.5 s. The snapshots at t = 11.5 s shows the run-up on the beach. Note that according with the positivity preservation property shown in the previous section, the simulation remains stable even with the use of a third order scheme. We also provide a comparison between the numerical results and the data taken from the experiments for both free surface evolutions at the wave gages locations on Fig. 22 and the velocity at the ADVs locations, see Fig. 23. We highlight that even with this relatively moderate level of refinement, the characteristics of the flow are well resolved. Considering [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF] (which is the only work, to our knowledge, to provide results on an unstructured mesh) we observe similar qualitative results with less elements. The results observable in [START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF][START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF] for instance seem to exhibit a slightly better agreement with the data, but using a meshes with more than 100,000 elements, highlighting the need for locally adapted meshes.

Conclusion

In this work, we introduce for the first time a fully discontinuous Galerkin formulation designed to approximate the solutions of a multi-dimensional Green-Naghdi model on arbitrary unstructured simplicial meshes. The underlying elliptic-hyperbolic decoupling approach, already investigated in [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF] for the d = 1 case, allows to regard the dispersive correction simply as an algebraic source term in the NSW equations. Such dispersive correction is obtained as the solution of scalar elliptic second order sub-problems, which are approximated using a mixed formulation and LDG fluxes. Additionally, the preservation of motionless steady states and of the positivity of the water height, under a suitable time step restriction, are ensured for the whole formulation and for any order of approximation. These properties are assessed through several benchmarks, involving nonlinear wave transformations and severe occurrence of dry areas, and some additional convergence studies are performed. The numerical method exhibits some very nice convergence properties and the numerical validations show that a good description of the wave's transformations can be achieved with a moderate number of mesh elements, thank's to the high order accuracy of the method and the ability to refine the meshes in some areas of interest.

Although the obtained numerical results clearly show very promising abilities for the study of nearshore flows, there are still several issues that need to be addressed in future works. In particular, we use a low-brow approach to stabilize the computations in the vicinity of broken waves. If such a method provides results which are in good agreement with experimental data for the cases under study in this work, we are still working to extend the method of [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] to the d = 2 case in a robust way.

Another important issue may be related to the computational cost of our approach. If we believe in the potential benefits of the use of such a non-conforming high order approach, especially concerning the robust approximations of solutions involving strong gradients occurring in the vicinity of breaking waves and possibly strong topography variations, it clearly leads to larger algebraic systems than those obtained with classical continuous approximations. Future works will therefore be devoted to go further in the development and optimization of our approach, and in particular to the reduction of the associated computational cost.

Figure 2 :

 2 Figure 2: Nodes locations for the Zhang and Shu quadrature -P 2 and P 3 cases.

Figure 3 :

 3 Figure 3: Test 1 -Motionless steady states preservation: topography and free surface

Figure 4 : 4 Figure 5 :

 445 Figure 4: Test 1 -L 2 -error for the free surface elevation vs. ∆x for k = 1, 2, 3 and k = 4 at t max = 0.02 s.

Figure 6 :

 6 Figure 6: Test 2 -Solitary wave propagation over a flat bottom: initial free surface

Figure 7 : 4 Figure 8 :

 748 Figure 7: Test 2 -L 2 -error for the free surface elevation vs. ∆x for k = 1, 2, 3 and k = 4 at t max = 0.2 s.

Figure 9 :

 9 Figure 9: Test 2 -Solitary wave propagation: cross section of the free surface elevation at several times during the propagation.

Figure 10 :

 10 Figure 10: Test2 -Solitary wave propagation: study of the total mass (M) and total energy (E) conservation. Time evolution of the relative errors |M (t) -M (0)|/M (0) and |E(t) -E(0)|/E(0) for k = 2, 3 on a mesh of 5000 elements.

Figure 12 :

 12 Figure 12: Test 4 -Propagation of highly dispersive waves: sketch of the experiment configuration and location of the gauges used for test #A.

Figure 13 :

 13 Figure 13: Test 4 -Propagation of highly dispersive waves: test #A, free surface evolution at gauges. Numerical data are denoted in plain lines.

Figure 14 :

 14 Figure 14: Test 4 -Propagation of highly dispersive waves: test #B, comparison with experimental data at gauges 1 to 8.

Figure 15 :

 15 Figure 15: Test 5 -Periodic waves propagation over an elliptic shoal: topography and view of the free surface.

Section 2 :

 2 x = 3 m , -5 m ≤ y ≤ 5 m , Section 3 : x = 5 m , -5 m ≤ y ≤ 5 m , Section 5 : x = 9 m , -5 m ≤ y ≤ 5 m , Section 7 : y = 0 m , 0 m ≤ x ≤ 10 m ,

Figure 16 :

 16 Figure 16: Test 5 -Periodic waves propagation over an elliptic shoal: comparison with experimental data along the four sections.

Figure 17 :

 17 Figure 17: Test 6 -Solitary wave propagation over a 3d reef: overview of the mesh in the refined area.

Figure 18 :

 18 Figure 18: Test 6 -Solitary wave propagation over a 3d reef: 3d view of the topography and initial condition.

Figure 19 :

 19 Figure 19: Test 6 -Solitary wave propagation over a 3d reef: location of the ADV and wave gauges.

Figure 20 :

 20 Figure 20: Test 6 -Solitary wave propagation over a 3d reef -free surface at times t = 4.5, 6.5, 8.5 and 9.5 s.

Figure 21 :

 21 Figure 21: Test 6 -Solitary wave propagation over a 3d reef -free surface at times t = 10.5, 11.5, 13 and 14 s.

Figure 22 :

 22 Figure 22: Test 6 -Solitary wave propagation over a 3d reef -Comparison with experimental solution at gauges 1 to 9.

Figure 23 :

 23 Figure 23: Test 6 -Solitary wave propagation over a 3d reef: x and y velocity components -Comparison with experimental solution at several ADVs.

) but with a sub-optimal convergence observed for the free surface with k = 4. Such a behavior is also observed with slightly larger amplitude waves.

Acknowledgements The second author acknowledges partial support from the ANR-13-BS01-0009-01 BOND.