Fabien Marche
email: fabien.marche@umontpellier.fr

Combined Hybridizable Discontinuous Galerkin (HDG) and Runge-Kutta Discontinuous Galerkin (RK-DG) formulations for Green-Naghdi equations on unstructured meshes

Keywords: Green-Naghdi equations, high-order schemes, Hybridized Discontinuous Galerkin, free surface flows, shallow water Equations, nonlinear dispersive equations

In this paper, we introduce some new high-order discrete formulations on general unstructured meshes, especially designed for the study of irrotational free surface flows based on partial differential equations belonging to the family of fully nonlinear and weakly dispersive shallow water equations. Working with a recent family of optimized asymptotically equivalent equations, we benefit from the simplified analytical structure of the linear dispersive operators to conveniently reformulate the models as the classical nonlinear shallow water equations supplemented with several algebraic source terms, which globally account for the non-hydrostatic effects through the introduction of auxiliary coupling variables. High-order discrete approximations of the main flow variables are obtained with a RK-DG method, while the trace of the auxiliary variables are approximated on the mesh skeleton through the resolution of second-order linear elliptic sub-problems with high-order HDG formulations. The combined use of hybrid unknowns and local post-processing significantly helps to reduce the number of globally coupled unknowns in comparison with previous approaches. The proposed formulation is then extended to a more complex family of three parameters enhanced Green-Naghdi equations. The resulting numerical models are validated through several benchmarks involving nonlinear waves transformations and propagation over varying topographies, showing good convergence properties and very good agreements with several sets of experimental data.

Introduction

The propagation and transformations of incompressible, homogeneous and inviscid nonlinear surface waves in nearshore areas are ideally governed by the free surface Euler equations. In particular, the mathematical and numerical modeling of such processes require an accurate description of both dispersive and strongly nonlinear effects. In its full generality, this problem remains mathematically and numerically challenging. To acquire a better understanding of major physical processes associated with the nonlinear and non-hydrostatic propagation over uneven bottoms, several improvements have been obtained recently 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 such models can be found in [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF]. The shallow-water regime is of particular interest in nearshore oceanography: the mean water depth H 0 is assumed to be small compared to the typical wave length λ:

(shallow water regime) µ := H 2 0 λ 2 1.
In such shallow-water regime, the weakly nonlinear regime has been widely studied for several decades: the typical amplitude a of the waves are small, in the sense that (weakly nonlinear regime) ε := a h 0 = O(µ).

In such shallow water and weakly nonlinear regime, depth averaged approximations of order O(µ 2 ) of the free surface Euler equations are furnished by the Boussinesq systems. Many types of (possibly enhanced) Boussinesq-like models, which are all O(µ 2 ) accurate approximations of the free surface Euler equations, have been derived, theoretically and numerically studied, see for instance and among others [START_REF] Beji | A formal derivation and numerical modelling of the improved Boussinesq equations for varying depth[END_REF][START_REF] Boussinesq | Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics: Part 2. a slowly-varying bathymetry[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF][START_REF] Peregrine | Long waves on a beach[END_REF]. All these approximations generally differ through their linear and nonlinear dispersion properties. However, the small nonlinearity assumption is too restrictive for many applications in coastal oceanography and such Boussinesqtype models are often used outside of their range of validity, especially when nonlinear effects become overriding, like in the shoaling area and in the vicinity of breaking. In such situations, fully nonlinear models should be used instead, see for instance [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]. In this paper, we therefore focus on the less restrictive fully-nonlinear and shallow water regime: (fully nonlinear regime) ε := a h 0 = O(1).

In this regime, the corresponding equations have been first derived in the horizontal surface dimension d = 1 in [START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF][START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF], and by Su and Gardner [START_REF] Su | Korteweg-de-Vries equation and generalizations. III. Derivation of the Kortewegde-Vries equation and Burgers equation[END_REF] and Green and Naghdi [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF] in the d = 2 case. These have been recently mathematically justified in [START_REF] Alvarez-Samaniego | Large time existence for 3d water-waves and asymptotics[END_REF]. Concerning the development of numerical approximations, the Green-Naghdi (GN) equations have only recently really received attention and various numerical methods have been introduced, mostly in the surface horizontal d = 1 case, like Finite-Differences (FD) approaches [START_REF] Antunes | Boussinesq and Serre type models with improved linear dispersion characteristics: Applications[END_REF], Finite-Volumes (FV) [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][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], WENO [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Dehghan | The solution of nonlinear Green-Naghdi equation arising in water sciences via a meshless method which combines moving kriging interpolation shape functions with the weighted essentially non-oscillatory method[END_REF], pseudo-spectral (PS) [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for fully-nonlinear 1d Serre equations[END_REF], (continuous) Finite-Elements (FEM) [START_REF] Filippini | A flexible genuinely nonlinear approach for wave propagation, breaking and runup[END_REF][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], FV and FEM methods on hyperbolic relaxed approximating models [START_REF] Favrie | A rapid numerical method for solving serre-green-naghdi equations describing long free surface gravity waves[END_REF][START_REF] Guermond | Robust explicit relaxation technique for solving the Green-Naghdi equations[END_REF] and discontinuous Galerkin approaches (possibly mixed with FEM) in [START_REF] Di Pietro | Weighted interior penalty discretization of fully nonlinear and weakly dispersive free surface shallow water flows[END_REF][START_REF] Dong | A reconstructed central discontinuous galerkin-finite element method for the fully nonlinear weakly dispersive Green-Naghdi model[END_REF][START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF][START_REF] Li | High order well-balanced CDG-FE methods for shallow water waves by a Green-Naghdi model[END_REF][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] Sharifian | A discontinuous Galerkin approach for conservative modelling of fully nonlinear and weakly dispersive wave transformations[END_REF]. In the far less studied d = 2 case, several methods have been developed on cartesian meshes: FD approaches [START_REF] Antunes Do Carmo | Numerical solution of the generalized Serre equations with the Mac-Cormack finite-difference scheme[END_REF][START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF][START_REF] Zhang | Validation of Boussinesq-Green-Naghdi modeling for surf zone hydrodynamics[END_REF], FV methods [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF], hybrid FV-FD methods [START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi 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] and WENO-FD [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF], a PS 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] and more recently a Hybridizable-DG method [START_REF] Samii | An explicit hybridized discontinuous Galerkin method for Serre-Green-Naghdi wave model[END_REF] and a Central DG-FE method in [START_REF] Li | A CDG-FE method for the two-dimensional Green-Naghdi model with the enhanced dispersive property[END_REF]. Numerical approximations of GN equations on general unstructured meshes are considered in [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on unstructured simplicial meshes[END_REF].

We have also recently introduced 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 enhanced GN models in order to optimize the dispersive properties of the classical GN equations, while allowing for easier and faster numerical approximations. High-order fully discontinuous Galerkin formulations based on the Local Discontinuous Galerkin (LDG) approach have been subsequently introduced for some of these models in [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF] for the d = 1 case and [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on unstructured simplicial meshes[END_REF] for the d = 2 case. In these works, the GN equations are written as coupled nonlinear (pseudo) hyperbolic-elliptic problems, relying on the (non-dispersive) NonLinear Shallow Water (NSW) equations [START_REF] De Saint-Venant | Théorie du mouvement non-permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF] supplemented by additional algebraic source terms, which fully accounts for the O(µ 2 ) nonlinear dispersive correction. These source terms are themselves computed as the solutions of a cascade of 2 auxiliary linear second order elliptic problems. The computational cost of the solutions of these elliptic problems is significantly alleviated through the use of a new reformulation of the regularizing second order operator that allows to decouple the time evolution of the two directional components and to assemble and factorize the associated matrix in a pre-processing step, while preserving the O(µ 2 ) asymptotic consistency with the Euler equations. We observe from these works that high-order discontinuous formulations are particularly well-suited for the approximation of the solutions of nonlinear and weakly dispersive equations. Indeed, it is well known that DG methods exhibit several well-known appealing features, such as local conservation, stability, the straightforward ability to handle arbitrary high-order polynomial approximations, a great flexibility regarding the underlying geometrical discretization or the compact stencils and minimal inter-element communications allowing highly parallelizable implementations. Beyond these general features, DG formulations provide a general and unified discrete framework allowing to accurately approximate both the hyperbolic and elliptic parts of the GN equations. Targeting nearshore oceanography large-scale applications, such a framework may allow to easily handle adaptive algorithms, since refining or coarsening a grid can be achieved without enforcing the continuity property commonly associated with the conforming elements. This appealing feature may be additionally exploited together with the use of higher order approximations in regular flow areas, e.g. far from wave breaking and submersion areas. Moreover, the strong robustness of DG methods in the vicinity of sharp gradients also appears to be well suited for the modelization of wave steepening and breaking, see for instance [START_REF] Richard | A new model of shoaling and breaking waves. Part 2. run-up and two-dimensional waves[END_REF]. However, when compared with the FEM method, the DG formulations introduced in [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on unstructured simplicial meshes[END_REF] require solutions of unsymmetric systems of linear equations with more unknowns for the same grids and are therefore more expensive in terms of computational cost and storage requirements. In particular, even if the use of the optimized GN formulations of [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] greatly alleviate the computational cost usually associated with the classical GN equations, there is still a significant computational overhead when compared with the solution of the classical hyperbolic NSW equations. This overhead is mostly due to the resolution of the global elliptic problems associated with the dispersive corrections. As an answer to this weakness, and in order to lower as much as possible this computational overhead, we focus in this work on the development of a combined Hybridized Discontinuous Galerkin (HDG) and RK-DG strategy. HDG methods, relying on the hybridization through the faces of the elements in combination with a Schur complement strategy, were introduced in [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed and continuous Galerkin methods for second order elliptic problems[END_REF] in the framework of diffusion problems. These methods were shown to share appealing features with the Raviart-Thomas (RT) and Brezzi-Douglas-Marini (BDM) mixed Finite-Elements methods, allowing equally efficient implementations and achieving optimal order of convergence for all the computed fields and possibly super-convergence of some of them, while retaining the inherent advantages of DG methods. Moreover, HDG methods are known to reduce the number of globally coupled degrees of freedom, when compared to other DG methods, relying on the introduction of intermediate face-based unknowns (the hybrid unknowns) in addition to the cell-based ones and the elimination of the cell-based unknowns by a static condensation process.

During the last decades, HDG methods have been successfully adapted to a wide range of applications, see for instance [START_REF] Cockburn | Static condensation, hybridization, and the devising of the HDG methods[END_REF] for a review, and we develop, in the following, some new discrete formulations especially designed for the optimized GN equations introduced in [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF]. In these formulations, we combine:

1 a RK-DG formulation to approximate the main flow variables through the resolution of the nonlinear (pseudo)-hyperbolic part of the systems, 2 symmetric HDG formulations to approximate the trace of the auxiliary/coupling variables issued from the elliptic sub-problems written in mixed form.

Both elliptic and (pseudo)-hyperbolic problems are coupled through the local reconstruction of corresponding in-cells auxiliary variables, which are then used to compute the non-hydrostatic/dispersive source terms through integration on mesh elements. We emphasize that the approximations of auxiliary variables traces on the mesh edges are obtained through the resolution of global sparse linear systems and that, for high-order polynomial approximations (i.e. k ≥ 2), such an hybrid approach leads to a significant decreasing of the size of these systems when compared to classical DG methods, like L-DG or IP-DG, as the resulting set of algebraic equations has globally coupled degrees of freedom only on the skeleton of the computational mesh. The in-cell reconstruction of these auxiliary variables may also explicitly benefit from a local element-by-element post-processing of the potential variables in the mixed formulations to lower the order of polynomial approximation of the dispersive correction, while preserving the observed order of convergence of the main flow variables. This conveniently allows to further decrease the size of the corresponding global algebraic systems without loosing the expected accuracy. Additionally, in comparison with [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on unstructured simplicial meshes[END_REF], the considered elliptic sub-problems rely on a slightly modified formulation of the dispersive operator that leads to symmetric formulations for both local and global problems in the HDG formulations. Such an hybrid formulation is also extended to a more complex family of three parameters optimized GN equations, allowing to consider a wider window of application with respect to frequency dispersion properties. Note that we purposely avoid to consider in this work the additional issues of vorticity and flows with wave-current interactions, as well as wave breaking and occurrence of dry areas. These particular topics will be addressed in subsequent works.

The remainder of this work is organized as follows: we briefly recall the considered mathematical models together with the associated initial-boundary problems and the corresponding notations in the next section. Section §3 is devoted to the introduction of the discrete settings and discrete formulations. The proposed approach is then validated in the last section, with convergence studies and comparisons with data taken from several experiments, with a particular attention paid to the study of waves refraction, diffraction and focusing over submerged shoals.

The physical models

v = (v 1 , v 2 ) ∈ R d as v(t, x) = 1 H ζ -H0+b
U hor (t, x, z)dz, and we denote by q = Hv the corresponding horizontal momentum.

H 0 ⇣(t, x) ⌘(t, x) b(x) z x H(t, x) = ⇣(t, x) + H 0 b(x) 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " w d O O n j c X k v j W S s e 1 5 L f R + a N r P D c = " > A A A H Z 3 i c j V T b b t N A E J 3 S Q E q 4 F i S E x E u g r W Q L x 7 J z h Y e g q r z 0 s U i U V k q i 4 s s m t e r E l r 0 u T a 3 8 I H / A H 9 C / Y G b s p L k V c B T v 7 O w 5 M 2 d m J 7 F D 3 4 u l Y f z a u L d Z u P + g u P W w 9 O j x k 6 f P n m + / + B Y H S e S I Y y f w g + j U t m L h e y N x L D 3 p i 9 M w E t b Q 9 s W J f f G Z z k 8 u R R R 7 w e i r H I e i N 7 Q G I 6 / v O Z Z E 1 9 n 2 5 s + u L Q b e K A 3 j 0 H N k E o m J U m l q h q o 0 t b p a 2 u u G 8 S D y 3 A 6 9 3 E D G 7 Y Z G J i d u 2 3 4 i e O t b t v A z X y T c W x f C Q 9 k r Y R A n i W U w 7 J D M H 5 4 r z 9 t m K L W + 5 / u x H P u i H Q e + 5 / I + i + J 7 g 3 M 5 i K x x L 8 3 Y 0 a U g X a a q V K q a q T d U x d Q M v a Y q V V x a q l L P d q h c / 6 j O K J i v L 0 M x a p s 9 h Y s y 6 M X 1 l S a E I j 0 s K p P h W v E 5 6 p 9 T q V M B z N B b j S w + r u u o g Z T C Z S j r q K k L w m s k u V L X q r y a u G 9 R C b W s h J p u U g l 8 S N j q L E N a + Y Q X U t c b u Q A 2 K X Y S J r K T + L 3 M Q x 1 J d w / P U m O y O 1 n k m q Q n
z 5 j T 3 K i X e Z p I 6 l 4 L a S l S K 3 e l u J J 2 P 7 2 a q M s x a t P s N c r U W F J n 3 q o j k 8 J O 9 X E i 1 t c i g d 3 1 q W 6 x e a J 0 1

1 b u x D T 1 O h e T 7 l 4 v H D T 0 D 1 o F p 4 A T z c h L p e R T k 3 V + j m x O J R 4 u 6 W u v 6 8 9 7 7 n R l R e R 8 V 5 r 5 q N X V N f 6 m X q M Z z A / W g W e / Q m P h s p s 4 J B 9 R p k H p x M i d + 9 2 W z p 7 v G L r B T 3 n V M H N j B / L n K N j e O I A u u B C A A w k M Q c A I J N o + W B D j p w M m G B C i r w c p + i K 0 P D 4 X M I E S 7 C E 7 Q Z x A j I X + C 3 w P c N f J v S P c U 9 S Y + Q 7 m 8 f E b I b e M 3 B D t B H d l z l f m 8 4 R j k 7 d 0 R + y U Y 5 K 6 M a 5 2 H m u I X g n n 6 P 0 X b 4 r 8 X x 7 V I q E P H 7 g G D 2 s K 2 U P V O X m U h P t C y s t z V U m M E K K P b B f P I 7 Q d Z k 4 7 X W Z O z L V T d y 0 + / 8 1 I 8 t L e y b E J 3 P y 1 O l J w A d d 8 L 3 d h Q l z 7 + K Z b k I j Q V j

One parameter optimized equations

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

∂ t ζ + ∇ • (Hv) = 0, I + T [H, b] (∂ t (Hv) + ∇ • (Hv ⊗ v)) + gH∇ζ + HQ 1 [H, b](ζ, v) = 0, (3)
where the operators T [H, b] and Q 1 [H, b] are defined as follows for all regular enough R d -valued function w:

T [H, b]w = R 1 [H, b](∇ • w) + R 2 [H, b](∇b • w), (4)
Q 1 [H, b](w) = -2R 1 [H, b](∂ 1 w • ∂ 2 w ⊥ + (∇ • w) 2 ) + R 2 [H, b](w • (w • ∇)∇b), (5)
(here ∂ 1 and ∂ 2 denote space derivatives along the two horizontal directions and for a given R 2 -valued function w = (w 1 , w 2 ) , we have w ⊥ = (-w 2 , w 1 ) ) with, for all smooth enough scalar-valued function w,

R 1 [H, b]w = - 1 3H ∇(H 3 w) - H 2 w∇b, (6)
R 2 [H, b]w = 1 2H ∇(H 2 w) + w∇b. (7)
We focus in this work on the asymptotically equivalent and computationally efficient formulation of the GN equations, introduced in [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF], called GN-LM equations as a shortcut in the following. Defining the water depth at rest, which does not depend on time:

H b = H 0 -b = H -ζ,
the GN-LM equations of [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] read as follows:

∂ t ζ + ∇ • (Hv) = 0, 1 + T[H b ] (∂ t (Hv) + ∇ • (Hv ⊗ v)) + gH∇ζ + Q[H, b](ζ, v) = 0, (8)
where the linear operator T[H b ] is defined as follows for all smooth enough R-valued function w:

T[H b ]w = -C[H b ] w H b , (9)
with

C[H b ]w = ∇ • (δ[H b ]∇w), and δ[H b ] = 1 3 (H b ) 3 , (10)
and

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

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

Q 3 [H, H b ]w = 1 6 ∇w∇(H 2 -(H b ) 2 ) + H 2 -(H b ) 2 3 ∆w - 1 6 ∆(H 2 -(H b ) 2 )w. ( 13
)
We also recall that the linear dispersion properties of ( 8) can be improved in several ways by adding some terms of order O(µ 2 ) to the momentum equation, which consequently do not affect the accuracy of the model. In [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF], an asymptotically equivalent enhanced family of models parametrized by α > 0, and referred to as α-GN-LM equations in the following, is given by

   ∂ t ζ + ∇ • (Hv) = 0, 1 + αT[H b ] ∂ t (Hv) + ∇ • (Hv ⊗ v) + α -1 α gH∇ζ + 1 α gH∇ζ + Q α [H, b](ζ, v) = 0, (14)
with

Q α [H, b](ζ, v) = H(Q 1 [H, b](v) + gQ 2 [H, b](ζ)) + Q 3 [H, H b ] 1 + αT[H b ] -1 (gH∇ζ) . ( 15
)
Remark 1. 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]

-α∇ • (δ[H b ]∇ 1 H b •) 0 0 1 -α∇ • (δ[H b ]∇ 1 H b •) . ( 16
)
From a numerical viewpoint, this structure allows to compute each component of the discharge q 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 2. Another difference with (3) is the presence of the modified quadratic term

Q α [H, b],
where no computation of third order derivative is needed. The price to pay is the inversion of extra linear systems, through the computation of

Q 3 [H, H b ] 1 + αT[H b ] -1 (gH∇ζ)
. However, this extra computational cost is largely off-set by the gain obtained by using the time independent scalar operator T[H b ], as shown in [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF].

Assuming that the water depth at rest H b is bounded away from 0, we see that for any sufficiently smooth scalar-valued function v:

(1+αT[H b ])v = β[H b ] v H b -∇ • δ α [H b ]∇ v H b , (17)
with δ α [H b ] := αδ[H b ], β[H b ] := H b . ( 18
)
Defining the total free surface elevation η = H + b, denoting W = (η, q) ∈ R 1+d the corresponding collection of primal variables, introducing two novel auxiliary R d -valued variables d and m, and provided that W, d, m are regular enough, the enhanced model ( 14) can be conveniently rewritten as

∂ t W + ∇ • F(W, b) + D α (W, d) = B(W, ∇b), (19a)
β[H b ]d -∇ • δ α [H b ]∇d = 1 α gH∇η + Qα [H, b](W, m), (19b)
β[H b ]m -∇ • δ α [H b ]∇m = gH∇η, (19c)
with

F(W, b) =   q q ⊗ q η -b + 1 2 g(η 2 -2ηb)I   , (20)
B(W, ∇b) = (0, -gη∇b) , D α (W, d) = 0, H b d - 1 α gH∇η , Qα [H, b](W, m) = H(Q 1 [H, b](v) + gQ 2 [H, b](ζ)) + Q 3 [H, b](H b m), (21)
where (19b) and (19c) should be respectively intended as the definition of the auxiliary variables d and m and I refers to the second order identity tensor. This last formulation clearly highlights that the dispersive correction operator D(W, d) only acts as a source term in the momentum conservation equation, and is obtained as the solution of a cascade of auxiliary second-order elliptic sub-problems. We recall that the particular form of the second component of the nonlinear flux F(W, b) comes from the prebalanced reformulation of the model, explicited in [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, using the following splitting of the hydrostatic pressure term:

gh∇ζ = 1 2 g∇(η 2 -2ηb) + gη∇b. ( 22
)
The use of the pre-balanced formulation allows to straightforwardly construct a discrete formulation that exactly preserves motionless steady states, as detailed in §3. [START_REF] Beji | Numerical simulation of nonlinear wave propagation over a bar[END_REF]

∇η • n = 0, on ∂Ω, (23a)
q = 0, on ∂Ω, (23b)
(δ α [H b ]∇d)n = 0, on ∂Ω, (23c)
(δ α [H b ]∇m)n = 0, on ∂Ω. (23d)
Basically, a vertical wall at the boundary may be modeled only through the assumption that there is no flux at the corresponding boundary, and thus only enforcing q = 0 at the corresponding boundary. It is however generally observed and admitted that the reflection at a vertical wall is equivalent to the head-on collision of two counter-propagating solitary waves of the same shape. In practice, during the interaction of a solitary wave with a wall, the normal component of the derivative ∇η on the boundary is indeed negligible, and the additional homogeneous conditions (23c) and (23d) may be introduced relying on similar symmetry argument. Although not rigorously justified, such a choice for the reflective boundary conditions is numerically validated in §4.

Remark 4. Anticipating on the construction of the combined HDG-RKDG formulations, we observe that an equivalent mixed formulation for the elliptic sub-problems can be straightforwardly obtained as follows:

∂ t W + ∇ • F(W, b) + D α (W, d) = B(W, ∇b), (24a)
δ α [H b ] -1 S d + ∇d = 0, (24b) β[H b ]d -∇ • S d = 1 α gH∇η + Qα [H, b](W, m), (24c)
δ α [H b ] -1 S m + ∇m = 0, (24d)
β[H b ]m -∇ • S m = gH∇η, (24e)
in which S d and S m are second-order tensor fluxes respectively associated with the potentials d, m. Note that the boundary conditions (23c)-(23d) may be reformulated as

S d n = 0, on ∂Ω, (25a)
S m n = 0, on ∂Ω. ( 25b
)
Remark 5. The original GN equations (3) may be supplemented by an energy conservation property, which may be formulated as follows, for a flat bottom:

d dt E(t) = 0, E(t) := 1 2 R gζ 2 + H|v| 2 + 1 3 H|H∇ • v| 2 dx. ( 26
)
We emphasize that model [START_REF] Boussinesq | Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond[END_REF], as well as optimized model [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], only preserve this energy up to O(µ 2 ) accuracy.

Remark 6. As pointed out in [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF], a simpler model that avoids the resolution of an extra linear system could be obtained by replacing

Q 3 [H, H b ] 1 + αT[H b ] -1 (gH∇ζ) in (15) by Q 3 [H, H b ] (gH∇ζ)
, keeping the same O(µ 2 ) accuracy. This would however leads to (linear) instabilities. This is due to the fact that Q 3 (ζ) contains third order derivatives in ζ that can create high frequencies 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 details on this issue.

Three parameters optimized equations

A three parameters optimized GN model has also been introduced in [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF]. It is indeed possible to improve the dispersive properties of the previous model without modifying its asymptotic accuracy. This can be achieved by introducing a set of three parameters α, θ, γ and considering the following change of variables for the velocity:

∀θ ≥ 0, v := 1 - θ H C[H] v θ ,
leading to the following set of equations:

     1 + γT[H b ] ∂ t η + ∇ • (Hv θ ) -θ∇ • C[H]v θ = 0, 1 + α(1 + θ)T[H b ] ∂ t (Hv θ ) + ∇ • (Hv θ ⊗ v θ ) + α -1 α gH∇η + 1 α gH∇η + Q α,θ [H, b](η, v θ ) = 0, (27) with
Q α,θ [H, b](η, v) = H(Q 1 [H, b](v) + gQ 2 [H, b](η)) + (1 + θ)Q 3 [H, H b ] 1 + α(1 + θ)T[H b ] -1 (gH∇η) + θQ 4 [H](v), (28)
and

Q 4 [H]v := -∇ • (C[H]v) v + 2 3 H 2 ∇ • (Hv)∆v + H∇(H∇ • (Hv)) • ∇v + 2 3 k=1,2 H 3 (∇v k • ∇)∂ k v + Hv k ∇ (H∂ k H) • ∇ v. ( 29
)
This model is referred to as (α, θ, γ)-GN-LM equations in the following. Taking θ = γ = 0, [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF] coincides 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]. It is possible to find optimized values of (α, θ, γ) to improve the dispersion properties of the equations, see the next section. In many practical configurations, when higher harmonics are released and larger values of relative depth |k|H 0 are reached, working with the (α, θ, γ)-GN-LM equations may help to improve the quality of the flow description with a better description of higher order harmonics interactions. Again, an equivalent mixed formulation can be straightforwardly obtained:

∂ t W θ + ∇ • F(W θ , b) + D α,θ,γ (W θ , χ, d) = B(W θ , b), (30a)
δ γ [H b ] -1 g χ + ∇χ = 0, (30b)
β[H b ]χ -∇ • g χ = θ∇ • C[H]v θ , (30c)
δ α,θ [H b ] -1 S d + ∇d = 0, (30d)
β[H b ]d -∇ • S d = 1 α gH∇η + Qα,θ [H, b](W, m), (30e)
δ α,θ [H b ] -1 S m + ∇m = 0, (30f) β[H b ]m -∇ • S m = gH∇η, (30g)
in which g χ is a first-order tensor flux associated with the scalar potential χ and

D α,θ,γ (W θ , χ, d) = H b χ, H b d - 1 α gH∇η .
Considering the case of an open bounded connected domain Ω ⊂ R 2 , with boundary ∂Ω and solid-wall boundary conditions, these equations may be supplemented with the following set of boundary conditions:

∇η • n = 0, on ∂Ω, (31a)
q = 0, on ∂Ω, (31b)
g χ • n = 0, on ∂Ω, (31c)
S d n = 0, on ∂Ω, (31d)
S m n = 0, on ∂Ω. (31e)

Dispersion properties

As already highlighted in our previous studies, the presence of the operator (1 + αT[H b ]) -1 for ( 14) and 27) make these models very robust with respect to high frequency perturbations. Considering, for the sake of simplicity, the case of flat bottom, we investigate the linear behavior of small perturbation ( ζ, V ) to the motionless constant state solution (ζ = 0, v = 0). The linear equations that govern these perturbations are

(1 + α(1 + θ)T[H b ]) -1 , (1 + γT[H b ]) -1 for (
     1 - γ 3 H 2 0 ∆ ∂ t ζ + H 0 ∇ • v - θ 3 H 3 0 ∇ • (∇ • ∇ v) = 0, 1 - α(1 + θ) 3 H 2 0 ∆ ∂ t v + α -1 α g∇ ζ + 1 α g∇ ζ = 0, (32)
Looking for plane wave solutions of the form (ζ 0 , v 0 )e i(k•x-ωt) to this linearized system, one finds the following dispersion relation

ω α,θ,γ (|k|) 2 = gH 0 |k| 2 (1 + (α-1)(θ+1) 3 (|k|H 0 ) 2 )(1 + θ+γ 3 (|k|H 0 ) 2 ) (1 + γ 3 (|k|H 0 ) 2 ))(1 + α(θ+1) 3 (|k|H 0 ) 2 ) . ( 33
)
We observe of course that setting θ = γ = 0 in (33) leads to the dispersion relation associated to the linearized version of ( 14)

ω α (|k|) 2 = gH 0 |k| 2 1 + α-1 3 (|k|H 0 ) 2 1 + α 3 (|k|H 0 ) 2 , (34)
while additionally setting α = 1 allows to recover the relation associated to the classical GN equations ( 3). These relations have to be compared with the one coming from the linear Stokes theory, (see for instance [START_REF] Lannes | The water waves problem: mathematical analysis and asymptotics[END_REF] for a full derivation from the Zakharov-Craig-Sulem formulation of the free surface water waves equations):

ω S (|k|) 2 = g|k| tanh(H 0 |k|). (35)
For the applications considered here, we are interested in obtaining a model with the widest possible range of validity and the discussion concerning the choices of (α, θ, γ) that improve the dispersive properties of the model follows the usual procedure, see for instance in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[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] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF]. For the one parameter α-GN-LM equations ( 14), we take α = 1.159 throughout this article.

Concerning model [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF], the two additional parameters offer the opportunity to optimize further the dispersive properties. Optimizing the errors in term of phase and group velocities, we show on Fig. 2 the ratio of GN phase velocities on Stokes phase velocity C α,θ,γ p /C S p and the ratio of GN group velocities on Stokes group velocity C α,θ,γ g /C S g with respect to the values of |k| for three different sets of parameters: (α, θ, γ) = (1, 0, 0) (the original GN equations), (α, θ, γ) = (1.159, 0, 0) (the α-GN-LM equations) and (α, θ, γ) = (1.018, 0.191, 0.097). We see that this last set of parameters offers very good approximation properties up to kH 0 = 10 for the phase velocity and kH 0 = 6 for the group velocity with an error smaller than 2%. sever Remark 7. Several techniques are available to improve the linear frequency dispersion of the models. We choose here to use a classical trick consisting in trading a time derivative of v for a space derivative of ζ in the higher order terms, together with changing the choice of the velocity unknown. Note that we did not find any modified model with two parameters that allows to improve the dispersive properties of the one-parameter model.

3 The discrete settings diameter of the element T . For all T ∈ T h , we denote by n T the unit outward normal on ∂T , and by x T its barycenter. Mesh faces are collected in the set F h , which is partitioned as

C GN g /C S g GN -GN-LM ( , , )-GN-LM optimized in dispersion Stokes
F h = F 0 h ∪ F ∂ h
, where F 0 h collects the internal faces (also called interfaces, i.e. there exist T 1 , T 2 ∈ T h such that F ⊂ ∂T 1 ∩ ∂T 2 ) and F ∂ h the boundary faces (i.e. there exists T ∈ T h such that F ⊂ ∂T ∩∂Ω). The length of a face F ∈ F h is denoted by |F |, the barycenter of a face F ∈ F h is denoted by x F , and we let h F denote the minimum length of the mesh elements to which F belongs. 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 . For any internal face F ∈ F 0 h , we choose an arbitrarily oriented but fixed unit normal n F , and we set n F := n T F for all boundary face F ⊂ ∂T ∩ ∂Ω. Additionally, T h is partitioned as

T h = T 0 h ∪ T ∂ h , with T 0 h = {T ∈ T h , F T ∩ F ∂ h = ∅} and T ∂ h = {T ∈ T h , F T ∩ F ∂ h = ∅}.

Functional and interpolation setting

Given an integer polynomial degree k ≥ 0, we consider the following discontinuous finite elements spaces:

P k h = P k (T h ) := v ∈ L 2 (T h ) | v |T ∈ P k (T ) ∀ T ∈ T h , (36)
where P k (T ) denotes the space of bivariate polynomials in T of total degree at most k, and

P k h := v ∈ (L 2 (T h )) d | v |T ∈ (P k (T )) d ∀ T ∈ T h , (37)
P k h := V ∈ (L 2 (T h )) d×d | V |T ∈ (P k (T )) d×d ∀ T ∈ T h , (38)
In addition, we introduce finite element spaces on the (interior) mesh skeleton:

M k h := µ ∈ L 2 (F 0 h ) | µ |F ∈ P k (F ) ∀ F ∈ F 0 h , (39)
which consists of functions which are continuous inside the interfaces, but discontinuous at their borders, and

M k h := µ ∈ (L 2 (F 0 h )) d | µ |F ∈ (P k (F )) d ∀ F ∈ F 0 h . (40)
We define the following inner products

v, w Ω := Ω v w, v, w T := T v w, v, w ∂T := ∂T v w ∀T ∈ T h , v, w F := F v w ∀F ∈ F h ,
for smooth enough scalar-valued functions v, w respectively defined on Ω, T , ∂T and F . Similar inner products are also respectively defined for smooth enough vector-valued functions v, w:

v, w Ω := Ω v • w, v, w T := T v • w, v, w ∂T := ∂T v • w ∀T ∈ T h , v, w F := F v • w ∀F ∈ F h ,
and for second order tensor functions V, W:

V, W Ω := Ω V :W, V, W F := F V :W ∀F ∈ F h , V, W T := T V :W, V, W ∂T := ∂T V :W ∀T ∈ T h .
We also define the mesh elements inner products as

v, w T h := T ∈T h v, w T , v, w T h := T ∈T h v, w T , V, W T h := T ∈T h V, W T , for v, w ∈ L 2 (T h ), v, w ∈ (L 2 (T h )) d , V, W ∈ (L 2 (T h )) d×d
, and the mesh elements boundaries inner products as µ, ν ∂T h :=

T ∈T h µ, ν ∂T , µ, ν ∂T h := T ∈T h µ, ν ∂T , for µ, ν ∈ L 2 (∂T h ), µ, ν ∈ (L 2 (∂T h )) d .
The following shortcuts are also defined:

µ, ν • n ∂T h := T ∈T h µ, ν • n T ∂T , Vn, ν ∂T h := T ∈T h Vn T , ν ∂T , ν, Vn ∂T 0 h := T ∈T 0 h ν, Vn T ∂T + T ∈T ∂ h F ∈F T ∩F 0 h ν, Vn T F F , ν, Vn ∂T ∂ h := T ∈T ∂ h F ∈F T ∩F ∂ h ν, Vn T F F , for µ ∈ L 2 (∂T h ), ν ∈ (L 2 (∂T h )) d , V ∈ (L 2 (∂T h ))
d×d (these two expressions respectively refers to inner products on interfaces and boundary mesh edges). For all T ∈ T h , we denote p k T the L 2 -orthogonal projector onto P k (T ) and p k T h the L 2 -orthogonal projector onto P h . Applying p k T or p k T h to R d -valued functions, we respectively define L 2 -orthogonal projectors onto (P k (T )) d and P h . Similarly, we denote I k T the element nodal interpolation into P k (T ). The corresponding nodal distributions in elements and edges are approximate optimal nodes introduced in [START_REF] Chen | Approximate optimal points for polynomial interpolation of real functions in an interval and in a triangle[END_REF], which have better approximation properties than equidistant distributions. The global I k T h interpolation into P h is obtained by gathering the local interpolating polynomials defined on each elements. Applying I k T or I k T h to R d -valued functions, we respectively define interpolations onto (P k (T )) d and P k h . For a given final computational time t max > 0, 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+1t n =: ∆t n . More details on the computation of the time step ∆t n and on the time marching algorithms are given in §3.4.3. For any sufficiently regular scalar-valued function of time w and R d -valued function w, we let w n := w(t n ) and w n := w(t n ).

Discrete gradient, divergence and Laplace operators

Following [START_REF] Dryja | On discontinuous Galerkin methods for elliptic problems with discontinuous coefficients[END_REF], we define the jump and average operators such that, for a sufficiently smooth function ϕ and an interior vertex F ∈ F i h such that F ⊂ ∂T 1 ∩ ∂T 2 for distinct mesh elements T 1 and T 2 ,

ϕ := ϕ |T1 -ϕ |T2 , { {ϕ} } F := 1 2 (ϕ |T1 + ϕ |T2 ), (41)
In what follows, and when no confusion can arise, we omit the subscript F from both v F and { {v} } F .

To discretize the linear and nonlinear operators that appear in our models, we need discrete counterparts of the gradient, divergence and of the Laplace operators applied to discontinuous polynomial functions. For any v h ∈ P k (T h ), we define the following global lifting of the jumps of v h , see for instance [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]):

R k h ( v h ) := F ∈F h r k F ( v h ),
where, for all F ∈ F h , the local lifting operator r k F ( v h ) ∈ P k h is defined as the unique solution of the following problem:

r k F ( v h ), ψ h Ω = v h , { {ψ h } }n F F ∀ψ h ∈ P h .
Following [25, Section 2.3], we define the discrete gradient operator G k h :

P k h → P k h such that, for all v h ∈ P k h , G k h (v h ) := ∇ h v h -R k h ( v h ). ( 42
)
This gradient has better asymptotic consistency properties than the discontinuous (element-by-element) gradient ∇ h , as it accounts for the jumps of its argument through the second contribution; see [

: P k h → P k h such that, for all v h ∈ P k h , G k,⊥ h (v h ) := G k h (v h ) ⊥ .
For the sake of simplicity, for v h ∈ P k h , we will denote by ∇ 2 h v h the second order tensor obtained from the double application of the discrete gradient G k h . For any integer k ≥ 1, we also define the discrete divergence operator D k h :

P k h → P k h such that, for all v h ∈ P k h with cartesian component (v h,i ) 1≤i≤d , D k h (v h ) := d i=1 ∇ h v h,i -R k h ( v h,i ) • e i ,
where e i denotes the i th vector of the cartesian basis on R d . The discrete divergence operator enjoys consistency property similar to the ones satisfied by the previous discrete gradient; see [START_REF] Di Pietro | Discrete functional analysis tools for discontinuous Galerkin methods with application to the incompressible Navier-Stokes equations[END_REF]Proposition 6.28].

We also introduce the discrete Laplace operator L k h :

P k h → P k h such that, for all v h ∈ P k h , L k h (v h ) = D k h (G k h (v h )).
Note that, by construction, we also have

L k h (v h ) = T r(∇ 2 h (v h )).
v ∈ H 1 0 (Ω) ∩ H k+1 (Ω), it holds inf v h ∈P k (T h ) ∇v -G k h (v h ) h k , inf v h ∈P k h ∇ • v -D k h (v h ) h k , inf v h ∈P k (T h ) ∆v -L k h (v h ) h k-1 ,
where a b means a ≤ Cb with real number C > 0 independent of the meshsize h.

The discrete formulations

A combined semi-discrete in space HDG-RKDG formulation

We consider in the following (k, r) ∈ N 2 , and we assume that k ≥ 2 and r ∈ {k -1, k}. The topography parameterization b and the associated water depth at rest H b are approximated through high-order polynomial interpolation respectively as

b h = I k T h (b) and H b h = I k T h (H b ), allowing to compute the required approximations G k h (b h ), ∇ 2 h b h and ∇ 3 h b h (
third order tensor associated with third order discrete derivatives of b h ), as well as

G k h (H b h ) and L k h (H b h
). The semi-discrete in space combined HDG-DG approximation of [START_REF] Cockburn | Conditions for superconvergence of HDG methods for second-order elliptic problems[END_REF] written in the mixed formulation [START_REF] Di Pietro | Discrete functional analysis tools for discontinuous Galerkin methods with application to the incompressible Navier-Stokes equations[END_REF] reads as follows: find

W h = (η h , q h ) ∈ P k h × P k h , (S d h , d h , d h ) ∈ P r h × P r h × M r h , (S m h , m h , m h ) ∈ P r h × P r h × M r h such that: ∂ t W h , ϕ h T h + A α (W h ), ϕ h T h = 0, ∀ϕ h ∈ P k h , (43a)
δ α [H b h ] -1 S d h , Φ h T h -d h , ∇ • Φ h T h + d h , Φ h n ∂T 0 h + d h , Φ h n ∂T ∂ h = 0, ∀Φ h ∈ P r h , (43b)
β[H b h ]d h , ψ h T h + ∇ • S d h , ψ h T h + ( S d h -S d h )n, ψ h ∂T h = Q α,h , ψ h T h , ∀ψ h ∈ P r h , (43c)
S d h n, ν h ∂T 0 h = 0, ∀ν h ∈ M r h , (43d)
S d h n T := S d h n T + S ∂T (d h -d h ) on F T ∩ F 0 h , 0 on F T ∩ F ∂ h , ∀T ∈ T h , (43e)
δ α [H b h ] -1 S m h , Ψ h T h -m h , ∇ • Ψ h T h + m h , Ψ h n ∂T 0 h + m h , Ψ h n ∂T ∂ h = 0, ∀Ψ h ∈ P r h , (43f)
β[H b h ]m h , φ h T h + ∇ • S m h , φ h T h + ( S m h -S m h )n, φ h ∂T h = gH h G k h (η h ), φ h T h , ∀φ h ∈ P r h , (43g)
S m h n , µ h ∂T 0 h = 0, ∀µ h ∈ M r h , (43h)
S m h n T := S m h n T + S ∂T (m h -m h ) on F T ∩ F 0 h , 0 on F T ∩ F ∂ h , ∀T ∈ T h . (43i)
where:

(i) (43a) is a DG semi-discrete formulation associated with (19a) and the nonlinear operator A α acting on any admissible discontinuous polynomial approximation

W h ∈ P k h × P k h is defined by A α (W h ), ϕ h T h := -F(W h , b h ), ∇ϕ h T h + Fn, ϕ h ∂T h + D α,h , ϕ h T h -B(W h , b h ), ϕ h T h , ∀ϕ h ∈ P k h , (44)
where Fn ∈ (P k (∂T h )) d+1 is such that for all T ∈ T h , Fn |∂T is a numerical approximation of the trace of the normal component of the interface flux F(W h , b h )n T on ∂T , whose precise expression will be given in Section 3.4.2 below, and D h is an in-cells discontinuous polynomial approximation of the dispersive correction D(W, d), defined by:

D α,h , ϕ h T h := 0, H b h d h - 1 α gH h G k h (η h ), ϕ h T h , ∀ϕ h ∈ P k h , (45)
the definition of d h is given below, (ii) (43b)→(43e) and (43f)→(43i) are HDG discrete formulations respectively associated with the elliptic subproblems written in mixed form (24b)→(30e) and (24d)→(24e). These HDG problems allow to compute polynomial approximations of the traces of the potentials d and m on the mesh interior edges. More precisely, (43e) and (43i) should be intended as the definitions of the trace of the normal component of the numerical fluxes S d h n T and S m h n T on ∂T , for all T ∈ T h , while (43d) and (43h) are transmission equations enforcing the fact that the traces of the normal numerical fluxes S d h and S m h are single valued on internal faces. These two sets of transmission equations allow to set-up the global problems respectively leading to the computation of ( d h , m h ) ∈ (M r h ) 2 . In a second time, (43b)-( 43c) and (43f)-(43g) are local element-by-element problems defined with (S d h , d h ) and (S m h , m h ) as unknowns, once the approximations of the traces on the mesh edges d h and m h are known. Note that this formulation only accounts for the case of homogeneous Neumann boundary conditions [START_REF] Devilla | Tutorial on hybridizable discontinuous Galerkin (HDG) for second-order elliptic problems[END_REF] corresponding to global solid-wall boundary conditions on all boundary edges. Although not used in this paper, the case of (one-directional) periodic boundary conditions may also be straightforwardly implemented by enforcing accordingly the continuity of the potential and fluxes through the connected (periodic) boundaries.

(iii) the nonlinear discrete operator Q α,h occurring in (43c) aims at approximating the nonlinear term (21) as follows:

Q α,h := 1 α gH h G k h (η h ) + H h Q 1,h [H h , b h ](p k T h ( q h H h )) + gH h Q 2,h [H h , b h ](η h ) + Q 3,h [H h , H b h ](m h ),
with for all u h ∈ P k h :

H h Q h 1 [H h , b h ](u h ) := 2H 2 h G k h (H h ) (G k h (u h ) e 1 ) • (G k h (u ⊥ h ) e 2 ) + D k h (u h ) 2 + 2 3 H 3 h G k h (G k h (u h ) e 1 )(G k h (u ⊥ h ) e 2 ) + G k h (G k h (u ⊥ h ) e 2 )(G k h (u h ) e 1 ) + 2D k h (u h )G k h (D k h (u h )) + H 2 h (G k h (u h ) e 1 ) • (G k h (u ⊥ h ) e 2 ) + D k h (u h ) 2 G k h (b h ) + 1 2 H 2 h 2G k h (u h )∇ 2 h b h u h + (∇ 3 h b h u h )u h + H h u h • (u h • G k h )G k h (b h ) G k h (b h ) + H h u h • (u h • G k h )G k h (b h ) G k h (H h ), (46)
for all ζ h ∈ P k h :

H h Q 2,h [H h , b h ](ζ h ) := -H 2 h G k h G k,⊥ h (ζ h ) G k,⊥ h (H h ) -H h G k h (H h ) G k h (b h ) • G k h (ζ h ) - 1 2 H 2 h ∇ 2 h ζ h G k h (b h ) + ∇ 2 h b h G k h (ζ h ) + H h H h 2 L k h (ζ h ) -G k h (b h ) • G k h (ζ h ) G k h (b h ), (47)
and:

Q 3,h [H h , H b h ](m h ) := 2 6 m h G k h (H b h ) + H b h δ α [H b ] -1 S m h H h G k h (H h ) -H b h G k h (H b h ) + 1 3 H 2 h -(H b h ) 2 L k h (H b h )m h + H b h D k h (δ α [H b ] -1 S m h ) + 2δ α [H b ] -1 S m h G k h (H b h ) - 2 6 G k h (H h ) • G k h (H h ) + L k h (H h ) -G k h (H b h ) • G k h (H b h ) -L k h (H b h ) H b h m h . (48)
, d h ) ∈ P k-1 h × P k-1 h and (S m h , m h ) ∈ P k-1 h × P k-1 h
relying on the optimal convergence of both the primal and the dual variables, see for instance [START_REF] Cockburn | Conditions for superconvergence of HDG methods for second-order elliptic problems[END_REF] and [START_REF] Stenberg | Postprocessing schemes for some mixed finite elements[END_REF]. These are respectively defined as the unique elements of P k h satisfying respectively the local problems

∇d h , ∇ψ h T = -δ α [H b h ] -1 S d h , ∇ψ h T , ∀ψ h ∈ P k h , ∀T ∈ T h , (49a)
d h , 1 T = d h , 1 T , ∀T ∈ T h , (49b)
and

∇m h , ∇ψ h T = -δ α [H b h ]S m h , ∇ψ h T , ∀ψ h ∈ P k h , ∀T ∈ T h , (50a)
m h , 1 T = m h , 1 T , ∀T ∈ T h . (50b)
(v) S ∂T ∈ (P 0 (∂T )) d×d is a second order tensor consisting of stabilization parameters for the HDG formulation.

In the following, we simply set S ∂T = I, ensuring both stability and optimal convergence rates for all the computed fields, see for instance [START_REF] Cockburn | A projection-based error analysis of HDG methods[END_REF].

Remark 8. From a practical viewpoint, the two sets of transmission equations through internal faces (43d) and (43h) lead to the same square and symmetric system of linear equations associated with the linear operator [START_REF] Cockburn | A projection-based error analysis of HDG methods[END_REF]. As this operator does not depend on time, the corresponding sparse matrix is assembled once and for all in a preprocessing step.

Remark 9. We emphasize that ( d h , m h ) ∈ M r h × M r h are globally defined discontinuous approximations of the traces of the auxiliary variables d and m over the (interior) mesh edges, that (d h , m h ) ∈ P r h × P r h are reconstructed discontinuous approximations over the mesh elements, while (d h , m h ) are reconstructed (postprocessed if r = k -1) discontinuous approximations over the mesh elements belonging to P k h × P k h . When r = k -1, we seek for traces d h , m h and auxiliary variables d h , m h in a smaller dimension polynomial space of degree k -1. Thanks to the local post-processing, we ultimately recover approximations (m h , d h ) ∈ (P k h ) 2 , which are respectively used to compute the right hand side of the local problem (43c) and the dispersive correction D α,h occurring in the nonlinear equations (43a), see [START_REF] Lannes | The water waves problem: mathematical analysis and asymptotics[END_REF].

Remark 10. The HDG formulations (43b)→(43e) and (43f)→(43i) rely on the fact that we only consider homogeneous Neumann boundary conditions on d and m, see [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]. Hence, we borrow the alternative formulation with Neumann local problems of [START_REF] Devilla | Tutorial on hybridizable discontinuous Galerkin (HDG) for second-order elliptic problems[END_REF], which consists in prescribing the Neumann boundary conditions directly in the local problem and considering the modified definition of the numerical traces of the fluxes (43e) and (43i) instead of the usual trace definition, see [START_REF] Nguyen | An implicit high-order hybridizable discontinuous Galerkin method for linear convection-diffusion equations[END_REF]. This formulation leads to slightly smaller global problems, as we only seek for unknowns on the interior edges, and may provide a better accuracy of the solution near Neumann boundaries, as pointed out in [START_REF] Devilla | Tutorial on hybridizable discontinuous Galerkin (HDG) for second-order elliptic problems[END_REF].

Remark 11. At first sight, formulation (43a)→(43i) may appear as unattractive, involving more unknowns and approximations than the initial problem (19a)→(19c). However, since d h and m h are only defined on the mesh edges, this approach allows to significantly reduce the total number of globally coupled unknowns needed to compute the dispersive correction D α,h , while preserving the overall accuracy of the approximation and the (expected) optimal convergence rates of the solutions, see §4 for numerical validations. Indeed, the hybridization of the method and the static condensation strategy allow for a computationally efficient implementation of the method, writing the resulting problem as two (diagonal bloc) algebraic systems of linear equations where only the traces d h and m h appear as globally coupled unknowns, eliminating the variables (S d h , d h ) and (S m h , m h ) through the computation of the Schur-complement matrix system for the trace unknowns. This leads to dim(M r h ) = (r + 1) card(F 0 h ) globally coupled unknowns, while the use of (non-hybridized) SWIP-DG or L-DG methods leads to dim(P k (T h )) = (k+1)(k+2) 2 card(T h ) coupled unknowns. More details about the sparsity of the associated stiffness matrices may be found in [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed and continuous Galerkin methods for second order elliptic problems[END_REF], in which bounds on the number of non-zero coefficients are provided and compared with those of IP-DG and L-DG methods. After solving these two global problems for d h and m h , (S d h , d h ) and (S m h , m h ) can respectively be inexpensively computed in an element-by-element fashion with the local problems (43c)-(43d) and (43g)-(43h). When r = k -1, the number of globally coupled faces unknowns is further reduced (with k card(F 0 h ) coupled unknowns), while the post-processing steps allow to recover polynomial approximations of degree k in mesh elements is also inexpensively computed in an element-by-element fashion.

Interface nonlinear fluxes

The high-order reconstructed numerical flux detailed in [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF] is a good default choice to approximate the interface fluxes Fn ∈ (P k (∂T h )) d+1 , leading to a well-balanced scheme that preserves motionless steady states. Considering T ∈ T h such and F ∈ F T (we only focus on the case F ∈ F i h and do not detail the case F ∈ F b h ) we denote W -= (η -, q-) and W + = (η + , q+ ) 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 * -η -) (51)
and

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

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

ϕ h ∂T h = T ∈T h F ∈F T Fn T F , ϕ h F , (55)
and

Fn T F = F h ( W-, W+ , b, b, n T F ) + F T F , (56)
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 -λ T h (W + -W -)), (57)
with λ T h = 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 ) . (58)
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 . (59)
Note that the modified interface flux (56) only induces perturbations of order k + 1 when compared to the regular Lax-Friedrichs flux straightforwardly applied to the (not reconstructed) traces.

Time discretization

Supplementing the α-GN-LM equations with an initial data W(0, •) = W 0 , and introducing the corresponding discrete initial data W 0 h = p k T h (W 0 ), the time stepping is carried out using explicit SSP-RK schemes [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, advancing W h from time level n to n + 1 with the third-order SSP-RK scheme involves two intermediate stages denoted W n,i h , i = 1, 2 and is computed as follows :

       W n,1 h = W n h -∆t n A α (W n h ) , W n,2 h = 1 4 (3W n h + W n,1 h ) -1 4 ∆t n A α (W n,1 h ), W n+1 h = 1 3 (W n h + 2w n,2 h ) -2 3 ∆t n A α (W n,2 h ), (60)
where the corresponding time step ∆t n is computed adaptively using the following CFL condition:

∆t n ≤ 1 2k + 1 min T ∈T h ( h T λ ∂T ), (61)
with λ ∂T defined in (58).

Well-balancing for motionless steady states

The preservation of motionless steady states can be deduced from the corresponding property obtained for the nonlinear (non-dispersive) shallow water equations, see for instance [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF], provided that the problems (43b)→(43e) and (43f)→(43i) are well-defined.

Proposition 1. The formulation [START_REF] Guennebaud | Eigen v[END_REF] together with the interface fluxes discretization (56) and a first order Euler time-marching algorithm preserves the motionless steady states, provided that the integrals of (43a) 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 , (62)
with η e constant, Proof. Assuming that the equilibriumW h = (η e , 0) holds, we have to show that

A α (W h ), ϕ h T h = 0, ∀ϕ h ∈ P k h .
We first observe that the local solvers (43b)-(43c)-( 43e) and (43f)-(43g)-(43i) are well defined. Indeed, looking at (43f)-(43g), we have for homogeneous conditions (i.e. m h = 0 and gH h G k h (η h ) = 0):

δ α [H b h ] -1 S m h , Ψ h T h -m h , ∇ • Ψ h T h + m h , Ψ h n ∂T ∂ h = 0, ∀Ψ h ∈ P r h , β[H b h ]m h , φ h T h + ∇ • S m h , φ h T h + ( S m h -S m h )n, φ h ∂T h = 0, ∀φ h ∈ P r h , S m h n T := S m h n T + S ∂T m h on F T ∩ F 0 h , 0 on F T ∩ F ∂ h , ∀T ∈ T h , which gives locally, for (Ψ h , φ h ) = (S m h , m h ): δ α [H b h ] -1 S m h , S m h T -m h , ∇ • S m h T + m h , S m h n T ∂T \∩Ω = 0, ∀T ∈ T h , β[H b h ]m h , m h T + ∇ • S m h , m h T + S ∂T m h , m h ∂T \∂Ω -S m h n T , m h ∂T ∩∂Ω = 0, ∀T ∈ T h ,
where the restriction of global functions to the considered element T are not explicitly written, for the sake of simplicity. Summing both equations straightforwardly yields to S m h = 0 and m h = 0. Then, from this definition of the local solver (43f)-(43g) and the trace (43i), the existence and uniqueness of the global solution mh

∈ M r h satisfying S m h n , µ h ∂T 0 h = 0, ∀µ h ∈ M r h ,
can be obtained following, for instance, the unified formalism of [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed and continuous Galerkin methods for second order elliptic problems[END_REF] for hybridized methods and is the trivial one.

Observing that, from definitions ( 42)-( 46)-( 47)-( 48), we have Q α,h = 0 at steady states, similar considerations lead to dh = 0 and hence

D α,h , ϕ h T h = 0 ∀ϕ h ∈ P k h ,
no matter what type of in-cell reconstructions d h and m h are used. From now, prooving that

-F(W h , b h ), ∇ϕ h T h + Fn, ϕ h ∂T h -B(W h , b h ), ϕ h T h = 0, ∀ϕ h ∈ P k h ,
at equilibrium follows for instance the lines of [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on unstructured simplicial meshes[END_REF].

Remark 12. Although not detailed in this work, a positivity preservation property for the water height can also be obtained following for instance the approach of [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF], already developed in [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on unstructured simplicial meshes[END_REF].

Remark 13. The study of stability and dispersive properties of the whole discrete formulations for our models is a difficult issue, even in the simplified case of linearized equations. This is mainly due to the occurrence of the term

Q 3 [H, H b ] 1 + αT[H b ] -1 (gH∇ζ) .
The existence and uniqueness result for the discrete elliptic problems relies on the discrete stability of the chosen HDG discretization, following [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed and continuous Galerkin methods for second order elliptic problems[END_REF]. Considering the classical (non-dispersive) linearized shallow water equations, several results closely related to the discrete stability and dispersive properties of DG methods may be found in [START_REF] Le Roux | Stability/dispersion analysis of the discontinuous galerkin linearized shallowwater system[END_REF].

Boundary conditions

Boundary conditions are imposed weakly. For equations (43a), the solid-wall conditions (23a) and (23b) are enforced through suitable values of W + at virtual exterior quadrature nodes on mesh boundaries F ∂ h , allowing to compute the corresponding interface fluxes in A α (W n h ). As shown in formulation (43a)→(43i), the corresponding homogeneous Neumann boundary conditions on the elliptic equations associated with the dispersive corrections are directly enforced in the local HDG problems (43b)-(43c) and (43f)-(43g), following the approach described in [START_REF] Devilla | Tutorial on hybridizable discontinuous Galerkin (HDG) for second-order elliptic problems[END_REF]. Although not used in this work, periodic boundary conditions on opposite domain boundaries, for instance for rectangular domains, can also be obtained by enforcing the continuity of the corresponding variables and fluxes. These simple boundary conditions possibly have to be supplemented 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 use a relaxation techniques and we enforce periodic waves combined with generation/absorption by mean of a generation/relaxation areas, following the ideas of [START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] and using the relaxation functions described in [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF]. The computational domain is therefore locally extended to include a generating layer at the inlet boundaries and a sponge layer at the outlet boundaries. Note that the generating layer may also be supplemented with a sponge layer in order to avoid possibly reflected waves. The relaxed solution along the domain is defined as follows:

W rel h = F a W h + (1 -F a )F g R(t)W enf , (65)
where F a , F g stand for the absorption and generation profiles and R governs the time evolution of the generation process. Above, W enf defines the targeted (enforced) wave profile. Concerning the relaxation functions, we follow [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF] and consider the following smooth profiles:

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 set to 3. The reader is referred to the above references for more details. Note that new generating boundary conditions for GN equations based on dispersive boundary layers, inspired from [START_REF] Lannes | Generating boundary conditions for a Boussinesq system[END_REF], are currently under study.

A combined HDG-RKDG formulation for the (α, θ, γ)-GN-LM equations

The combined semi-discrete in space HDG-DG problem ( 43) can be extended to approximate the solutions of the (α, θ, γ)-GN-LM equations [START_REF] Dubiner | Spectral methods on triangles and other domains[END_REF]. In the following, we only give details concerning the additional terms and corresponding modifications and do not recall parts of the formulation similar to those of [START_REF] Guennebaud | Eigen v[END_REF]. The main differences between ( 24) and ( 30 ]. These operators are built relying on the same HDG method as in [START_REF] Guennebaud | Eigen v[END_REF], considering now two different sets of continuity equations through mesh edges (see the vectorial systems (66d)-(66h) and the scalar system (66l) besides).

Still assuming k ≥ 2 and r ∈ {k-1, k}, this extended discrete formulation reads as follows: find

W θ h = (η h , q θ h ) ∈ P k h × P k h , (g χ h , χ h , χ h ) ∈ P r h × P r h × M r h , (S d h , d h , d h ) ∈ P r h × P r h × M r h , (S m h , m h , m h ) ∈ P r h × P r h × M r h such that: ∂ t W θ h , ϕ h T h + A α,θ,γ (W θ h ), ϕ h T h = 0, ∀ϕ h ∈ P k h , (66a)
δ α,θ [H b h ] -1 S d h , Φ h T h -d h , ∇ • Φ h T h + d h , Φ h n ∂T 0 h + d h , Φ h n ∂T ∂ h = 0, ∀Φ h ∈ P r h , (66b) β[H b h ]d h , ψ h T h + ∇ • S d h , ψ h T h + ( S d h -S d h )n, ψ h ∂T h = Q α,θ,h , ψ h T h , ∀ψ h ∈ P r h , (66c)
S d h n, ν h ∂T 0 h = 0, ∀ν h ∈ M r h , (66d)
S d h n T := S d h n T + S ∂T (d h -d h ) on F T ∩ F 0 h , 0 on F T ∩ F ∂ h , ∀T ∈ T h , (66e)
δ α,θ [H b h ] -1 S m h , Ψ h T h -m h , ∇ • Ψ h T h + m h , Ψ h n ∂T 0 h + m h , Ψ h n ∂T ∂ h = 0, ∀Ψ h ∈ P r h , (66f)
β[H b h ]m h , φ h T h + ∇ • S m h , φ h T h + ( S m h -S m h )n, φ h ∂T h = gH h G k h (η h ), φ h T h , ∀φ h ∈ P r h , (66g)
S m h n , µ h ∂T 0 h = 0, ∀µ h ∈ M r h , (66h)
S m h n T := S m h n T + S ∂T (m h -m h ) on F T ∩ F 0 h , 0 on F T ∩ F ∂ h , ∀T ∈ T h , (66i)
δ γ [H b h ] -1 g χ h , π h T h -χ h , ∇ • π h T h + χ h , π h • n ∂T 0 h + χ h , π h • n ∂T ∂ h = 0, ∀π h ∈ P r h , (66j)
β[H b h ]χ h , ϕ h T h + ∇ • g χ h , ϕ h T h + ( g χ h -g χ h ) • n, ϕ h ∂T h = θM h (v θ h ), ϕ h T h , ∀ϕ h ∈ P r h , (66k)
g χ h • n, κ h ∂T 0 h = 0, ∀κ h ∈ M r h , (66l)
g χ h • n T := g χ h • n T + S ∂T (χ h -χ h ) on F T ∩ F 0 h , 0 on F T ∩ F ∂ h , ∀T ∈ T h , (66m) where:
(i) (43a) is a DG discrete formulation associated with (19a) and the discrete nonlinear operator A α,θ,γ acting on any admissible discontinuous polynomial approximation W θ h ∈ P k h × P k h is defined by

A α,θ,γ (W θ h ), ϕ h T h := -F(W θ h , b h ), ∇ϕ h T h + Fn, ϕ h ∂T h + D α,θ,γ,h , ϕ h T h -B(W θ h , b h ), ϕ h T h , ∀ϕ h ∈ P k h , (67)
where D α,θ,γ,h ∈ P k h × P k h is a numerical approximation of the dispersive correction D α,θ,γ (W θ , χ, d), defined by:

D α,θ,γ,h , ϕ h T h := H b h χ h , ϕ h T h , H b h d h - 1 α gH h G k h (η h ), ϕ h T h , ∀ϕ h ∈ P k h , (68)
with the definition of χ h given below, (ii) (66j)→(66m) is an additional HDG formulation associated with (30b)-(30c). χ h ∈ P k h is defined as follows:

-if r = k, we set χ h = χ h , -if r = k -1, χ h is
defined as a locally post-processed super-convergent discontinuous polynomial approximation of the auxiliary variable χ obtained from (g χ h , χ h )

∈ P k-1 h × P k-1 h
. More precisely, χ h is defined as the unique element of P k h satisfying

∇χ h , ∇ψ h T = -δ γ [H b h ] -1 g χ h , ∇ψ h T , ∀ψ h ∈ P k h , ∀T ∈ T h , (69a)
χ h , 1 T = χ h , 1 T . ∀T ∈ T h , (69b)
(iii) the additional third order source term M h in the mass evolution equations is defined as follows:

M h (v θ h ) := 1 3 6H h G k h (H h ) G k h (v θ h ) G k h (H h ) + 3H 2 h G k h (H h ) • L k h (v θ h ) + H 3 h D k h (L k h (v θ h )) + 3H 2 h G k h (v θ h ) : ∇ 2 h H h + G k h (H h ) • G k h (D k h (v θ h )) , (70)
where : stands for the inner product with two indices, (iv) the nonlinear discrete operator Q α,θ,γ,h occurring in (66c) aims at approximating [START_REF] Dong | A reconstructed central discontinuous galerkin-finite element method for the fully nonlinear weakly dispersive Green-Naghdi model[END_REF] as follows:

Q α,θ,γ,h := 1 α gH h G k h (η h ) + H h Q 1,h [H h , b h ](p k T h ( q θ h H h )) + gH h Q 2,h [H h , b h ](η h ) + (1 + θ)Q 3,h [H h , H b h ](m h ) + θQ 4,h [H h ](v θ h ), (71)
where Q 4,h is a discrete operator built upon definition (29) as follows:

Q 4,h [H h ](v θ h ) := -M θ,h (v θ h )v θ + 2 3 H 2 h D k h (q θ h )L k h (v θ h ) + H h G k h (v θ h ) G k h (p k T h (H h G k h (q θ h )) + 2 3 H 3 h G k h (G k h (v θ h )) : G k h (v θ h ) + q θ h,1 G k h (v θ h ) G k h • p k T h (H h G k h (H h ) • e 1 ) + q θ h,2 G k h (v θ h ) G k h • p k T h (H h G k h (H h ) • e 2 ) . (72)
We have the following well-balancing result:

Proposition 2. The formulation (66) together with the interface fluxes discretization ( 56) and a first order Euler time-marching algorithm preserves the motionless steady states, providing that the integrals of (66a) are exactly computed for the motionless steady states.

Proof. The proof follows the line of Proposition 1, accounting for the additional elliptic problem (66j)→(66m).

Numerical results

In this section, we assess the previous formulations with several benchmarks. Unless stated otherwise, we use the α-GN-LM equations with α = 1.159. We consider solid-wall boundary conditions on ∂Ω and the time step restriction is computed according to [START_REF] Nguyen | An implicit high-order hybridizable discontinuous Galerkin method for linear convection-diffusion equations[END_REF]. Some numerical accuracy and convergence analysis are performed in §4.3 using the L 2 norm defined, for any arbitrary scalar valued piecewise polynomial function w h ∈ P k h , as follows:

w h 2 T h = w h , w h T h .
When generating/absorbing layers are needed, we add a generating/absorbing layer of 3λ in at the inlet boundary and an absorbing layer of 2λ in at the outlet boundary, where λ in is the incident wave-length, see §3.4.5 for details.

For several test cases, we conveniently summarize the corresponding computational set up in Tables. We recall that k refers to the polynomial order of the in-cell DG approximations of the main flow variables (η h , q h ), r refers to the polynomial order of the HDG approximations of the auxiliary variables ( dh , mh ) (and possibly χh is the (α, θ, γ)-GN-LM equations are used) on the mesh skeleton, |T h | refers to the number of mesh elements, |F 0 h | refers to the number of interior mesh edges, h min and h max refers respectively to the minimum and maximum edges lengths of the mesh. For test cases 4.3, 4.4, 4.5 and 4.6, which rely on extruded one-dimensional channels with one dimensional flow features, we apply a rotation of angle π/4 with respect to the horizontal direction on the corresponding meshes in order to activate all the components of the flow variables during the computations.

Remarks on implementation

The combined HDG-RKDG formulations [START_REF] Guennebaud | Eigen v[END_REF] and [START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF], as well as the d = 2 extension of the SWIPDG-RKDG formulation of [START_REF] Di Pietro | Weighted interior penalty discretization of fully nonlinear and weakly dispersive free surface shallow water flows[END_REF] (which is provided in Appendix for the sake of completeness), have been implemented into our C++ WaveBox framework. The linear algebra is handled with the Eigen library [START_REF] Guennebaud | Eigen v[END_REF]. The linear systems arising from the discrete formulations are solved using the supernodal Cholesky CHOLMOD library [START_REF] Chen | Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate[END_REF]. The proposed formulations may be implemented on general polytopal unstructured meshes and do not rely on any particular properties of the polynomial basis functions. Yet, the numerical results shown in this section are obtained using structured (for the convergence studies) and unstructured simplicial meshes and the Dubiner basis functions [START_REF] Dubiner | Spectral methods on triangles and other domains[END_REF] with an implementation relying on the nodal polymorphic approach of [START_REF] Gassner | Polymorphic nodal elements and their application in discontinuous galerkin methods[END_REF] for the sake of efficiency. Every piece of stationary information is computed once, in a set-up pre-processing step. In particular, the sparse matrices issued from the global problems associated with the elemental transmission conditions (43d) and (43h), the dense matrices associated with the local solvers (43b)-( 43c) and (43f)-(43g), as well as the local post-processing ( 49)-( 50) operators are assembled and (Choleski) factorized once and for all in a preliminary set-up step. Yet, at this level of implementation, the emphasize was put on developing a unified framework with a generic interface for wave flow problems and no particular attempt has been made to improve the code performance. In particular, we do not investigate yet the solver efficiency or the parallelization possibilities of the HDG formulations for the computation of local problems. In contrast with [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on unstructured simplicial meshes[END_REF] in which a quadrature free implementation was used for the computations of the discrete version of the nonlinear operator [START_REF] De Saint-Venant | Théorie du mouvement non-permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF], we carefully compute every piece of polynomial integration with an exact quadrature rule and high-order rules are used for the integration of non-polynomial terms to lower the associated aliasing errors. For the tests shown in the following, no additional stabilization mechanism is needed. In particular, we do not need to introduce any filtering to eliminate short wave/high-frequency components generated by the nonlinear interactions.

Preservation of motionless steady state

This preliminary test case is devoted to briefly check the ability of the formulation to preserve motionless steady states. 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 = (- 1 3 , - 1 3
) and

x 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 , (73)
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. 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 = 100 s using a k = 3 approximation are respectively 2.40e-15, and 7.2e-14 for η and q.

Solitary wave propagation over a flat bottom

We consider the time evolution of a solitary wave profile over a flat bottom b = 0, defined as follows:

η(t, x) = H 0 + εH 0 sech 2 (κ(x • e 1 -ct)) , q(t, x) = c (η(t, x) -H 0 ) , (74)
with κ = 3ε 4H 2 0 (1 + ε)
and c = gH 0 (1 + ε). Note that ( 74) is an exact solution of the original GN equations

(3), but only a solution of the GN-LM models up to O(µ 2 ) terms. However, for small enough values of ε, such profiles are expected to propagate over flat bottoms without noticeable deformations, possibly with a small dispersive tails propagating backward. The computational domain is a channel [0, 200] × [0, 5] with a reference water depth H 0 = 1 m. A solitary wave of relative amplitude ε = 0.1 is initially centered at x 0 = 80 m. To study the h-convergence properties of the formulation (43), we consider a sequence of regular triangular meshes, with a level of mesh refinement defined by the number of equi-spaced segments along each side of the domain, with mesh size h ranging from h 0 to h 0 /2 p , with h 0 = 1.25 m and p = 3, together with polynomial expansions of degrees ranging from k = 1 to k = 4. The time step is chosen small enough to ensure that the dominant component of the error is related to the spatial discretization. As an exact solution is not available, we use a reference solution which computed on a refined mesh of size h = h 0 /2 4 and using k = r = 6. We compute the L 2 -errors at t max = 0.1 s, the h-convergence curves and the corresponding convergence rates obtained by linear regression for both the free surface elevation η and the discharge q. Two possible formulations are investigated: Table 1: Solitary wave propagation over a flat bottom: number of coupled freedom degrees with respect to the polynomial order of approximation and the chosen formulation for the elliptic problems Figure 3: Solitary wave propagation over a flat bottom: initial free surface k ≥ 1 and r = k: we seek the trace of the auxiliary variables in the same space as the one used for the main flow variables (η, q). The corresponding results are reported on Fig. 4 , k ≥ 2 and r = k -1: we seek the trace of the auxiliary variables in a smaller space than the one used for the main flow variables, leading to smaller global linear systems. The expected order of convergence are then recovered relying on local post-processing, see Fig. 5.

We assume k ≥ 2 in §3.4 to avoid the lowest order approximation r = k -1 = 0 for the trace space, as no super-convergent post-processing is available. However, even if not used in practice as this formulation produces too much numerical dissipation and dispersion errors for the study of dispersive problems, the choice r = k = 1 is of course allowed and leads to the expected convergence rate, as reported on Fig. 4. As expected, we observe overall rates ranging between O(h k+ 1 2 ) and O(h k+1 ) for both variables and both formulations, although the L 2 -errors on η are generally one order of magnitude smaller than those on q. For the finer mesh of the sequence (|T h | = 12800 and |F 0 h | = 18784), we provide in Table 1 a comparison of the number of freedom degrees associated with the elliptic sub-problems for both SWIP-DG with order k, HDG formulation of order r = k and HDG formulation of order r = k -1. In particular, we point out that the same level of accuracy may be reached with either choices r = k and r = k -1, but with significantly smaller global systems than with the SWIP-DG formulation.

h k = 1 k = 2 k = 3 k = 4
Figure 4: Solitary wave propagation over a flat bottom -L 2 -error for the free surface η (left) and the discharge q (right) vs. mesh size h for k = 1, 2, 3, 4 and r = k at t = 0.1 s.

Head-on collision of solitary waves

We consider now the head-on collision of two identical solitary waves propagating in opposite directions [START_REF] Craig | Solitary water wave interactions[END_REF].

The collision of the two waves implies rapidly modifies the nonlinear dispersion characteristics of the flow and the discrete formulation has to ensure the equilibrium between amplitude and frequency dispersion to propagate the wave profile at constant shape and speed. The computational domain is a channel of 300 m length and 2 m width. We initialize the computation with two solitary waves (74) of relative amplitude ε = 0.2 initially located at x = -50 m and x = 50 m and with opposite velocities. We set the order of polynomial approximations to (k = 2, r = 1) and use a regular mesh of triangular elements, with edges lengths ranging from h min = 0.2 m to h max = 0.25 m. We show in Figure 6 some snapshots of the free surface at various times during the propagation with equations [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]. The obtained numerical results are in agreement with [START_REF] Mitsotakis | On the Galerkin/Finite-Element Method for the Serre equations[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] and those exhibited in [START_REF] Craig | Solitary water wave interactions[END_REF] using Euler equations. As in the d = 1 case, we observe some minor differences with the results obtained with the original GN equations. The maximum amplitude during the collision is qualitatively similar and there is no phase shift. We notice, however, some small variations in the amplitude of the oscillations in the dispersive tails before and after the collision, mainly due to the fact that ( 74) is only an approximated solution of (14).

Shoaling and reflection of a solitary wave over a sloping beach

Next, we investigate the reflection of a solitary wave on a sloping beach. The aim of this test is to study the shoaling and full reflection of a non-breaking solitary wave propagating above a regular sloping beach, before reaching a vertical solid wall, following the set-up of [START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF]. We consider the computational domain Ω = [-55, 20] × [0, 5] (in m) and a topography defined with the following slope:

s(x) = 0 if x ≤ 0, 1/50 if x ≥ 0. ( 75
)
The water depth at rest is H 0 = 0.7 m, the solitary wave is initially centred at x = -30 m, we use a mesh with a characteristic size h = 0.2 and set (k = 2, r = 1). We show a sketch of a cut of the topography and initial free surface on Figure 7. We compute time series of the free surface at a location near the vertical wall (x = 17.75 m) and compare numerical results with experimental data taken from [START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF]. Two runs are performed with ε = 0.1 ( A) and ε = 0.171 ( B) and the results are shown on Figure 8. We clearly identify the two expected peaks corresponding respectively to the incident and reflected waves, with a very good matching between simulations and experimental data for both the wave's celerity and amplitude. In comparison with the numerical results obtained with the original GN equations, we observe a slightly more oscillating dispersive tail for the second (more nonlinear) case, in agreement with the fact that the initial solitary wave is not an exact solution of the modified α-GN-LM equations. We now study the Dingemans case [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF] which highlights the ability of the models to compute the propagation and the interaction of higher order harmonics. Using the set-up introduced in [4], we compute the propagation of regular periodic waves over a submerged bar, see Figure 9. For this test, we generate periodic waves 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.

When the incident wave encounters the upward part of the bar, it shoals and steepens, which generates higherharmonics as the nonlinearity increases. These higher-harmonics are then freely released on the downward slope and propagates after the bar. For the numerical computation, we consider a channel of 40 m long and 1 m wide. The corresponding computational set-up is provided in Table 2. With k = 2, the use of hybrid HDG-DG formulation allows to halve the number of globally coupled freedom degrees needed for the computation of the dispersive correction with the SWIP-DG method of [START_REF] Di Pietro | Weighted interior penalty discretization of fully nonlinear and weakly dispersive free surface shallow water flows[END_REF] or the L-DG method of [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on unstructured simplicial meshes[END_REF].

k r |T h | |F 0 h | h min h max 2 1 5400 8409 0.11 0.17 Comparisons are performed between the two sets of parameters corresponding to the α-GN-LM model and the (α, θ, γ)-GN-LM model with 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 Figure 10. As already pointed out in our previous study [START_REF] Di Pietro | Weighted interior penalty discretization of fully nonlinear and weakly dispersive free surface shallow water flows[END_REF], we observe that the α-GN-LM is not able to provide an accurate free surface evolution at the last wave gauge 11. This last gauge is the most discriminating one as the higher-harmonics are completely released and interact together. We observe the improvements obtained using the (α, θ, γ)-GN-LM model at the last two wave gauges. It is worth highlighting that for this case, the use of the (α, θ, γ)-LM-GN equations allows to obtain the same quality of approximation than the one obtained for instance with the two-layers model of [START_REF] Lynett | A two-layer approach to wave modeling[END_REF], with a computational overhead of less than 10%. with Γ(y) = y(6.096y). We focus here on the computation of case B: T = 2.0 s and a = 0.0075 m where T refers to the wave train period and a the amplitude, leading to µ 1 2 = 0.117 and = 0.0165. The numerical domain is defined as [-8, 36] × [0, 6.096] (meters). The computational set-up used for this test is reported in Table 3. We compute the waves propagation up to t max = 100 s and we extract time series of wave elevation at several locations along the centerline during the last 25 s. The time series of the surface elevation is analyzed in the frequency domain to obtain the first, second and third harmonic amplitudes and we consider their spatial evolution. These numerical results are compared to the experimental data on The hybrid formulation allows to obtain very good agreements with the measurements, using a very coarse mesh. The energy transfert towards the second and third harmonics through the focusing is accurately reproduced. As a comparison, the results shown in [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] and obtained with our previous FV-WENO formulation on cartesian grid relies on a mesh with equidistant spacing δ x = δ y = 0.12, and the results obtained in [START_REF] Liu | An unstructured FEM model based on Boussinesq equations and its application to the calculation of multidirectional wave run-up in a cylinder group[END_REF] and [START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] with FE formulations respectively rely on 55018 (structured) and 30705 (unstructured and adapted) triangular elements. Of course, this low number of mesh elements is partially balanced by the number of freedom degrees in each elements but we also emphasize that the chosen formulation leads to global systems coupling only 22875 freedom degrees. This number has to be compared with the 46458 global freedom degrees needed by the SWIP-DG formulation or L-DG method of [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on unstructured simplicial meshes[END_REF] to reach the same level of polynomial approximation of the dispersive correction. It is always tricky to fairly compare the efficiency of different numerical methods, as many optimization may be achieved or not relying on the chosen implementation. However, relying on the same background implementation and the same computational set-up, denoting by t IPDG and t HDG the total CPU times (elapsed times on an otherwise-idle system) respectively needed to reach the final time of the computations respectively for the SWIPDG-RKDG formulation of order k and the proposed HDG-RKDG formulation, we measure that (t IPDGt HDG )/t IPDG ≈ 0.23. This is a mean value obtained from 3 different runs, and we recall that the global matrices are assembled in a preprocessing step, so that we consider the solve rather than the matrix setup CPU time. We now 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 refraction and 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. 13). Introducing the rotated coordinates We investigate two different numerical configurations, which are summarized in Table 4. For this particular test case, the formulation resulting with the choice (k = 2, r = 1) lead to under-resolved computation, even with the use of higher order post-processing for the coupling variables, and was therefore not retained here. Time series of the free surface elevation are obtained at several locations along the sections, between t = 30 s and t = 50 s. The signal is then analyzed with the zero up-crossing method to isolate single waves and compute mean wave heights, which are normalized by the incoming (targeted) wave heights (2a). We show on Fig. 14 a comparison between numerical results and experimental data along the sections, obtained with the set-up 1.

x r = x cos(20 • ) -y sin(20 • ) , y r = x sin(20 • ) + y cos(20 • ),
Although not shown here, very similar results are obtained with the set-up 2. We observe the expected wave focusing behind the shoal and due to the sloping beach, the waves continue shoaling after passing the shoal. We obtain very satisfying results, comparable with other studies found in the literature, see for instance [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 | Time-dependent numerical code for extended Boussinesq equations[END_REF] but with significantly coarser meshes. As in the previous test case, averaging the measured CPU times on 3 different runs, we measure that (t IPDGt HDG )/t IPDG ≈ 0.14 for set-up 1, and 0.22 for set-up 2.

Waves propagation over a circular shoal: Chawla and Kirby experiments

We now present the comparison of numerical results with a test issued from the experiments performed by Chawla and Kirby in the directional wave basin at the University of Delaware in 1996, see [START_REF] Chawla | Wave transformation over a submerged shoal[END_REF]. The physical We consider the propagation of a monochromatic train of waves with H 0 = 0.45, a wave height at the input boundary of 2.36 cm and a wave period T = 1.0 s, corresponding to the Test 4 of [START_REF] Chawla | Wave transformation over a submerged shoal[END_REF]. Accounting for boundary layers, the computational domain is 25.5 m long and 18.2 m wide. The physical experiment provides measured significant wave heights along seven sections covering most the areas of the submerged shoal and behind the shoal, see again Fig. 16. A schematic view of the domain and the gauges locations is shown on Fig. 16. The numerical set-up used for this test is summarized in Table 5. period and H 0 are close or equal to those in the previous test case, the resulting wave fields for these two experiments exhibits several differences. Indeed, the focusing of waves is much stronger in this second case and occurs on the top of the shoal, while the focus of waves in the previous test case occurs at a location further behind the shoal. This focusing is well reproduced here, with results close to those shown in [START_REF] Chen | Boussinesq modeling of wave transformation, breaking, and runup. II: 2 D[END_REF] and obtained with mesh sizes of 0.05 m and 0.1 m respectively in the x and y directions. Our model reproduces very well the transverse (asymmetric) variations of the wave field resulting from the effects of combined refraction/ diffraction and the slightly off-centered shoal location.

k r |T h | |F 0 h | h min h max 2

Waves propagation over an elliptic shoal: Vincent and Briggs experiments

We conclude this section with a computation corresponding to the M1 experiment of Vincent and Briggs [START_REF] Vincent | Refraction-diffraction of irregular waves over a mound[END_REF] and which aims at showing preliminary computational results obtained with the (α, θ, γ)-GN-LM equations. Here again, we study the transformations of a non-breaking monochromatic train of waves over a submerged shoal. The experimental basin is 25 m long and 30 m large. The elliptic shoal is centered at (x 0 , y 0 ) = (6.1, 13.72) and is 6.1 m long in the x-direction and 7.92 m wide in the y-direction. Its perimeter is defined by xx 0 3.05 with H 0 = 0.4572. Following [START_REF] Lynett | A two-layer approach to wave modeling[END_REF], we generate a wave train at the left boundary with a wave height of 0.048 m and a wave period T = 1.3 s. The computational set-up used for this test is detailed in Table 6. A snapshot of the quasi-steady state free surface taken at t = 32 s is shown on Figure 17. As expected, we observe that oblique wave interactions dominate the wave-field behind the shoal, creating an irregular sea surface. We ran this test case with the two sets of equations: α-GN-LM and (α, θ, γ)-GN-LM with the parameters given in section §2.3. The significant wave heights (in cm) are extracted from time series of the free surface taken along the section located at y = 13.72 m and ranging from x = 2.0 m to x = 17 m, hence covering the wave variations over and behind the shoal. The results obtained for both models are then compared with the measurements taken from sections 7-9 of [START_REF] Vincent | Refraction-diffraction of irregular waves over a mound[END_REF] and shown on Fig 18. As in Test §4.6, we observe that the higher-frequency waves generated due to nonlinear transfers on the shoal strongly interacts behind the shoal, and if the two models give satisfactory agreements with the data from the experiment, we observe that the (α, θ, γ)-GN-LM equations provide a slightly better matching regarding the maximum wave amplification in the vicinity of the shoal. The amplitude after the shoal is over-predicted by both models. We also note that the discrete formulation associated with the (α, θ, γ)-GN-LM equations seems to be more sensible with respect to high frequency perturbations and further studies regarding additional stabilization are required for this more elaborated 3 parameters model to be used with confidence for long time simulations.

Conclusion

In this work, we further investigate the design of high-order discontinuous formulations for fully nonlinear and weakly dispersive shallow water equations. We introduce a new combined HDG-RKDG method for the asymptotically equivalent models of [START_REF] Lannes | Nonlinear wave-current interactions in shallow water[END_REF]. The main flow variables are approximated by in-elements polynomial approximations, while the coupling variables are approximated on the mesh skeleton, contributing to alleviate the computational overhead associated with the non-hydrostatic effects. This also offers the possibility to use different kind of polynomial approximations for each problem, and we take benefit from the existence of superconvergent post-processing for the auxiliary variables to further reduce the size of the global algebraic problems associated with the discrete trace space of the dispersive correction. The computational efficiency of this new hybrid approach allows to consider its extension to the further optimized three parameters (α, θ, γ)-GN-LM model of [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF]. Our numerical investigations show that the expected convergence properties are recovered, together with interesting approximation capabilities for the study of waves transformations in d = 2 horizontal configuration and an enjoyable robustness for the α-GN-LM model, as no additional stabilization mechanism is needed. On-going and future works will investigate the enforcement of positivity preservation through sub-element resolution to accurately model run-up, wave-breaking and occurence of vorticity following the preliminary works of [START_REF] Lannes | Nonlinear wave-current interactions in shallow water[END_REF] and parallelized implementations [START_REF] Filippini | UHAINA: A parallel high performance unstructured near-shore wave model[END_REF].

if F ∈ F i h is such that F = ∂T 1 ∩ ∂T 2 , ξ F κ T if F ∈ F b h is such that F = ∂T ∩ ∂Ω, (77)
where ξ is defined in Remark 15, and with the weighted average operator defined such that, for a sufficiently smooth function ϕ and an interior face F ∈ F i h such that F ⊂ ∂T 1 ∩ ∂T 2 for distinct mesh elements T 1 and T 2 , { {ϕ} } ω,F := ω 2 ϕ |T1 + ω 1 ϕ |T2 , ω i := κ Ti κ T1 + κ T2 ∀i ∈ {1, 2}.

Note that when κ ≡ C in Ω for some real number C > 0, we have ω 1 = ω 2 = 1 2 , and also the subscript ω is omitted. We also introduce the associated discrete Laplace operator L a,k h : P k h → P k h such that, for all v h ∈ P k h , L a,k h (v h ) solves -L a,k h (v h ),

ψ h T h = a h (1; v h , ψ h ) ∀ψ h ∈ P k h .
Under a mesh quasi-uniformity assumption, it can be proved that, for any v ∈ H 1 0 (Ω) ∩ H k+1 (Ω), it holds inf

v h ∈P k (T h ) ∆v -L a,k h (v h ) h k-1 .
We assume in the following k ≥ 1. The semi-discrete in space SWIP-DG approximation of ( 19) reads as follows: find W h = (η h , q h ) ∈ P k h × P k h , (d h , p h , m h ) ∈ P k h × P k h × P k h such that:

∂ t W h , ϕ h T h + A h (W h ), ϕ h T h = 0, ∀ϕ h ∈ P k h , (79)
d h , ψ h T h + gH h G k h (η h ), ψ h T h = H b h p h , ψ h T h , ∀ψ h ∈ P k h , (80)
a h (δ α [H b h ]; m h , φ h ) + β[H b h ]m h , φ h T h = gH h G k h (η h ), φ h T h , ∀φ h ∈ P k h , (81)
a h (δ α [H b h ]; p h , π h ) + β[H b h ]p h , π h T h = Qα,h , π h T h , ∀π h ∈ P k h , (82)
where (i) the discrete nonlinear operator A h in (43a) is defined by

A h (W h ), ϕ h T h := -F(W h , b h ), ∇ϕ h T h + Fn, ϕ h ∂T h + D h , ϕ h T h -B(W h , ∇b h ), ϕ h T h , ∀ϕ h ∈ P h , (83)
where Fn ∈ (P k (∂T h )) d+1 is still defined as in [START_REF] Liu | An unstructured FEM model based on Boussinesq equations and its application to the calculation of multidirectional wave run-up in a cylinder group[END_REF], and the discrete dispersive correction D h is

D h , ϕ h T h := 0, d h , ϕ h T h , ∀ϕ h ∈ P h . (84)
(ii) The discrete nonlinear operator Qα,h in ( 82) is defined by

Qα,h := gH h G k h (η h ) + H h Q 1,h [H h , b h ](p k T h ( q h H h )) + H h Q 2,h [H h , b h ](η h ) + Q 3,h [H h , H b h ](m h ),
with Q 1,h and Q 2,h defined according to [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF] and [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF] and

Q 3,h [H h , H b h ](m h ) := 2 6 m h G k h (H b h ) + H b h G k h (m h ) H h G k h (H h ) -H b h G k h (H b h ) + 1 3 H 2 h -(H b h ) 2 L a,k h (H b h )m h + H b h L a,k h (m h ) + 2G k h (m h )G k h (H b h ) - 2 6 G k h (H h ) • G k h (H h ) + L a,k h (H h ) -G k h (H b h ) • G k h (H b h ) -L a,k h (H b h ) m h . ( 85
)
Remark 14. The definition of Q 3,h implies the computation of the second order tensor G k h (m h ), as in contrast with the previous HDG formulation, the SWIP-DG method does not provide approximations of the flux variable S m . Remark 15. The penalty coefficient ξ F used in (77) denotes a user-defined parameter sufficiently large to ensure coercivity. Following [START_REF] Epshteyn | Estimation of penalty parameters for symmetric interior penalty Galerkin methods[END_REF], this coefficient is defined as

ξ F := 6δ α [H b h ] 2 δ α [H b h ] 2 k(k + 1) max T ∈T h cot θ T , (86)
where δ α [H b h ] and δ α [H b h ] are respectively the global lower and upper bounds of δ α [H b h ] on Ω and cot θ T is the cotangent of the smallest angle in a triangle T (we choose not to distinguish interior from boundary faces for the sake of simplicity).

Remark 16. An extended SWIP-DG formulation can be accordingly designed for the (α, θ, γ)-GN-LM equations. This has been implemented and used in §4.10 to estimate the computational speed-up observed when using our new hybrid formulations.

  Let us denote by x = (x, y) ∈ R d the horizontal variables, where d = 2 denotes the horizontal dimension, 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 Figure1. Denoting by U hor the horizontal component of the velocity field in the fluid domain, we define the vertically averaged horizontal velocity

Figure 1 :

 1 Figure 1: Free surface flow description: main notations

3. 1

 1 Domain partitionningLet Ω ⊂ R d , with d = 2, denotes an open bounded connected polygonal domain with boundary ∂Ω. We consider a geometrically conforming partition T h of Ω defined as a finite collection of |T h | nonempty open disjoint triangular elements T of boundary ∂T such that Ω = T ∈T h T and we set ∂T h = {∂T, T ∈ T h }. We denote |T | the area of T and the partition is characterized by the meshsize h := max T ∈T h h T , where h T is the

Figure 2 :

 2 Figure 2: Dispersion relations of optimized linearized GN equations: normalized phase velocities (top) and group velocities with respect to wave numbers.

(

  iv) d h ∈ P k h and m h ∈ P k h are discontinuous polynomial approximations of the auxiliary variables d and m (see (19b) (19c)) computed from the trace approximations d h ∈ M r h and m h ∈ M r h and which exact definition depends on the chosen value of r ∈ {k -1, k}. More precisely: if r = k, we set d h = d h and m h = m h , if r = k-1, d h and m h are defined as locally post-processed super-convergent polynomial approximations of d and m respectively obtained from (S d h

  ) are on one hand the occurrence of a new second order and scalar elliptic problem in the first equation of (30a) (the mass conservation equations, which is therefore not an exact conservation equations but an approximate conservation equation up to O(µ 2 ) terms) and on the other hand the occurence of an additional source term, namely Q 4 , in the definition of the fully nonlinear term Qα,θ [H, b]. One also has to deal with the occurence of optimization parameters γ and θ. Considering the construction of an associated discrete formulation, this leads to the construction of two different linear discrete operators, respectively associated with 1 + γT[H b ] and 1 + α(1 + θ)T[H b

  SWIP-DG (k) HDG (r=k) HDG (r=k-1

4 Figure 5 :

 45 Figure 5: Solitary wave propagation over a flat bottom -L 2 -error for the free surface η (left) and the discharge q (right) vs. mesh size h for k = 2, 3, 4 and r = k -1 with local-post-processing of the auxiliary variables at t = 0.1 s.

Figure 6 :

 6 Figure 6: Head-on collision of solitary waves: free surface elevation at several times

Figure 7 :

 7 Figure 7: Shoaling and reflection of a solitary wave on a sloping beach: topography and initial free surface

Figure 8 :Figure 9 :

 89 Figure 8: Shoaling and reflection of a solitary wave on a sloping beach: time series of the free surface at x = 17.75 m for the α-GN-LM model

Table 2 :

 2 Propagation of periodic highly dispersive waves: the Dingemans experiments -computational set-up

Figure 10 :

 10 Figure 10: Propagation of periodic highly dispersive waves: the Dingemans experiments -Time series of the free surface at the last 4 wage gauges. Comparison between experimental data at wave gauges (o) and numerical results (-).

Figure 11 :

 11 Figure 11: Waves focusing by a topographical lens: Whalin experiments -3D view of the topography

  Fig 12.

54 Table 3 :

 543 k r |T h | |F 0 h | h min h max 3 2 5162 7625 0.21 0.Wave focusing by a topographical lens: Whalin experiments -computational set-up

Figure 12 :

 12 Figure 12: Waves focusing by a topographical lens: Whalin experiments -Comparison between experiment (case B) and numerical results for the amplitudes of the first three harmonics

  the topography is defined by b = 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 13 :

 13 Figure 13: Periodic wave propagation over an elliptic shoal: : Berkhoff experiment -the topography

40 Table 4 :

 404 k r |T h | |F 0 h | h min h max set-up 1: 2 2 45562 68025 0.14 0.23 set-up 2: 3 2 15222 22645 0.22 0.Periodic waves propagation over an elliptic shoal: Berkhoff experiments -computational configurations.

Figure 14 :

 14 Figure 14: Periodic waves propagation over an elliptic shoal: Berkhoff experiments -Comparison of the computed mean wave height with the experimental data along the height wave gauges sections

25 Table 5 :

 255 2 30082 44865 0.15 0.Periodic waves propagation over a circular shoal: Chawla and Kirby experiments -computational set-up We collect time series of wave elevation at several locations corresponding to the following sections section A : {2 m ≤ x ≤ 14 m, y = 8.98 m}, section B : {x = 11.12 m, 2 m ≤ y ≤ 16 m}, section C : {x = 9.65 m, 2 m ≤ y ≤ 16 m}, section D : {x = 7.995 m, 2 m ≤ y ≤ 16 m}, section E : {x = 6.35 m, 4 m ≤ y ≤ 14 m}, section F : {x = 5.075 m, 4 m ≤ y ≤ 14 m}, section G : {x = 3.885 m, 4 m ≤ y ≤ 14 m}, between t = 30 s and t = 50 s. Again, mean wave heights are extracted and normalized by the initial incoming wave height, leading to the comparisons shown on Fig.15. As pointed out in [12], even if the values of wave

Figure 15 :

 15 Figure 15: Waves propagation over a circular shoal: Chawla and Kirby experiments -Comparison of the computed normalized mean wave height with the experimental data along the seven transects

2 + y -y 0 3.96 2 = 1 ,

 221 and the corresponding topography is defined by b(x, y) = -H 0 + 0.7620(1 -x -

Table 6 :

 6 Periodic waves propagation over an elliptic shoal: Vincent and Briggs experiments -computational set-up

Figure 17 :

 17 Figure 17: Waves propagation over an elliptic shoal: Vincent and Briggs experiments -contour view of the free surface (up) and 3d view (down) at t = 35 s

Figure 18 :

 18 Figure 18: Waves propagation over an elliptic shoal: Vincent and Briggs experiments -significant wave heights along sections 7-9

  24, Theorem 2.2] for further insight into this point. In a similar way, a discrete gradient operator G k h : P k h → P k h may also be straightforwardly defined for all R d -valued function v h ∈ P k h with applying (42) on each scalar directional component of v h , together with a discrete orthogonal gradient G k,⊥ h

Acknowledgements This work was initiated while the author was visiting Inria BSO center, CARDAMOM team, France. The author particularly thanks Inria BSO and Dr M.Ricchiuto for his kind support and hospitality. The author also thanks Pr F.J. Sayas and Dr A. Duran for useful advices regarding implementation and gratefully acknowledges the support provided under Agence Nationale de la Recherche (ANR) projects HHOMM (ANR-15-CE40-0005), NABUCO (ANR-17-CE40-0025) and CNRS-LEFE-MANU projects UHAINA and D-Wave.

Appendix -A SWIP-DG formulation for the α-GN-LM equations In this section, for the sake of completeness, we extend to the d = 2 case the Symmetric Weighted Interior Penalty Discontinuous Galerkin (SWIP-DG) discrete formulation introduced in [START_REF] Di Pietro | Weighted interior penalty discretization of fully nonlinear and weakly dispersive free surface shallow water flows[END_REF] in the d = 1 case. Let κ ∈ L ∞ (Ω) denote a uniformly positive scalar valued function and set κ T := κ |T for all T ∈ T h . We consider the following bilinear form a h (κ; •, •) on P h × P h :

with diffusion-dependent penalty coefficient