G L Richard

A Duran

B Fabrèges

A new model of shoaling and breaking waves. Part II. Run-up and two-dimensional waves

We derive a two-dimensional depth-averaged model for coastal waves with both dispersive and dissipative effects. A tensor quantity called enstrophy models the subdepth large-scale turbulence, including its anisotropic character, and is a source of vorticity of the average flow. The small-scale turbulence is modelled through a turbulent-viscosity hypothesis. This fully nonlinear model has equivalent dispersive properties to the Green-Naghdi equations and is treated, both for the optimization of these properties and for the numerical resolution, with the same techniques which are used for the Green-Naghdi system. The model equations are solved with a discontinuous Galerkin discretization based on a decoupling between the hyperbolic and non-hydrostatic parts of the system. The predictions of the model are compared to experimental data in a wide range of physical conditions. Simulations were run in one-dimensional and two-dimensional cases, including run-up and run-down on beaches, non-trivial topographies, wave trains over a bar or propagation around an island or a reef. A very good agreement is reached in every cases, validating the predictive empirical laws for the parameters of the model. These comparisons confirm the efficiency of the present strategy, highlighting the enstrophy as a robust and reliable tool to describe wave breaking even in a two dimensional context. Compared with existing depth-averaged models, this approach is numerically robust and adds more physical effects without significant increase in numerical complexity.

Introduction

The need of accurate models in coastal engineering motivated many works in the past decades. The difficulties met by the researchers lie in the fact that the capability of the model to capture the main physical phenomena must be accompanied by an easy and reliable numerical resolution. A successful approach must combine an accurate and physically relevant model with a robust and efficient numerical scheme, both being mutually dependent. The main physical effects to model are the dispersive effects and the dissipative effects. The dispersion brings the most acute difficulties in the numerical resolution because it typically introduces third-order derivatives in the equations. Moreover the classical dispersive equations like the Green-Naghdi equations [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF]) lack dissipative terms whereas the even more classical Saint-Venant equations, also known as the nonlinear shallow water equations (Barré de Saint-Venant 1871), are nondispersive. The dispersive effects are dominant before breaking, in the so-called shoaling zone of the coastal waves propagation, and the dissipation dominates after breaking in what is called the surf zone. One of the challenge in coastal wave modelling is to derive a model capable of describing both the dispersion and the dissipation and of predicting accurately the breaking point. The zone of the beach which is not permanently covered by water is the swash zone where take place the phenomena of run-up and run-down which are also challenging for the numerical resolution.

The first dispersive equations used for coastal waves were the weakly nonlinear equations derived by [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] for the propagation of a solitary wave in water of constant depth, extended by Peregrine (1967) for water of variable depth using the depth-averaged velocity. However the weak nonlinear assumption reduces the validity of these equations. The dispersive properties were improved by different techniques, notably by [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF] and also by [START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF] who considered the velocity at an arbitray distance from the still water level instead of the depth-averaged velocity. Models which were able to dispense with the weak nonlinearity hypothesis were developed by [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF] and were adopted for coastal waves propagation under the name of fully nonlinear Boussinesq models in spite of the fact that the scope of these models exceeds by far the original Boussinesq equations. The dispersive properties of these equations, while better than those of weakly dispersive models, are not completely satisfactory because of their weakly dispersive character. They were in turn improved by various means such as considering the velocity at an arbitray depth [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part I: highly nonlinear unsteady waves[END_REF] or by using asymptotically equivalent equations [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]. New Green-Naghdi systems, asymptotically equivalent to the standard Green-Naghdi equations, but having a mathematical structure more suited to the numerical resolution, were proposed by [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF] and, in the one-dimensional case, by do [START_REF] Do Carmo | An improved Serre model: Efficient simulation and comparative evaluation[END_REF]. Rotational effects were included by [START_REF] Zhang | Boussinesq-Green-Naghdi rotational water wave theory[END_REF]. Reviews on this subject can be found in [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF] and [START_REF] Kirby | Boussinesq models and their application to coastal processes across a wide range of scales[END_REF].

Different strategies were implemented to add dissipative effects to the dispersive models in order to describe properly the breaking waves. We focus here mostly on two-dimensional models and we refer to the Part I of this work [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF] for one-dimensional models of coastal waves and to [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF] for a complete overview. In the roller model of [START_REF] Madsen | Surf zone dynamics simulated by a Boussinesq-type model. Part I: model description and cross-shore motion of regular waves[END_REF] the surface roller of breaking waves is considered as a volume of water being carried by the wave with the wave celerity. The position of the breaking point is found by a breaking criterion involving the local slope of the surface elevation. This approach requires also the geometrical determination of the roller, the determination of the roller celerity and a complex calibration of various parameters.

In another approach the energy dissipation in the breaking waves is modelled by introducing an eddy viscosity in the equations. In the model of [START_REF] Chen | Boussinesq modeling of wave transformation, breaking and runup[END_REF] two empirical parameters are used to determine the onset and cessation of breaking and the implementation of the breaking model in two horizontal dimensions requires the determination of the wave direction in order to estimate the age of a breaking event. A Smagorinsky-type subgrid model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. The basic equations[END_REF]) is used to account for the effect of the resultant eddy viscosity on the underlying flow [START_REF] Chen | Boussinesq modeling of a rip current system[END_REF]. In some other eddy-viscosity models the viscosity is calculated from a turbulent kinetic energy for which a semi-empirical transport equation with source term is solved [START_REF] Nwogu | Numerical prediction of breaking waves and currents with a Boussinesq model[END_REF][START_REF] Zhang | Rotational surf zone modeling for O(µ 4 ) Boussinesq-Green-Naghdi systems[END_REF].

Removing the dispersive terms in the Green-Naghdi equation leads to the Saint-Venant equations which produce discontinuities that dissipate energy. This is at the basis of the hybrid or switching method [START_REF] Tonelli | Simulation of wave breaking over complex bathymetries by a Boussinesq model[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF]) which switch off the dispersive terms when some breaking criterion is satisfied and which switch on these terms when a criterion for the end of breaking is activated.

For all these approaches, a breaking criterion is needed and many criteria were used such as a relative trough Froude number [START_REF] Okamoto | The Relative Trough Froude Number for initiation of wave breaking: Theory, experiments and numerical model confirmation[END_REF], a gradient of momentum [START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF], a combination of a local evaluation of the mechanical energy dissipation, a maximal front slope and a critical Froude number [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] or a combination of the surface variation and the local slope angle [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF]. The numerical shock detector of [START_REF] Krivodonova | Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws[END_REF] was used as a breaking criterion by [START_REF] Duran | Discontinuous Galerkin discretization of a new class of Green-Naghdi equations[END_REF] for the switching strategy and a numerical resolution by a Discontinuous Galerkin method.

The switching method was numerically implemented by various schemes (hybrid finitevolume/finite difference scheme in [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF], finite volume and finite element in [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF], Discontinuous Galerkin finite element scheme in Duran &[START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF][START_REF] Sharifian | A discontinuous Galerkin approach for conservative modelling of fully nonlinear and weakly dispersive wave transformations[END_REF]. The major drawbacks of the switching approach are the mesh grid sensitivity and the nonphysical oscillatory effects due to the switching of the dispersive terms (see [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF][START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF]). These numerical oscillations may increase with the mesh refinement and can be damped by a smooth transition to the hyperbolic regime. Further the numerical wave breaking detection involves the calibration of a set of empirical parameters.

Alternative approaches include the Saint-Venant equations with non-hydrostatic pressure corrections, which avoid the high-order derivatives in the Boussinesq or Green-Naghdi equations with a calculation of the vertical distribution of the non-hydrostatic pressure (see for example [START_REF] Stelling | An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation[END_REF][START_REF] Lu | A two-dimensional depth-integrated non-hydrostatic numerical model for nearshore wave propagation[END_REF], or the semi-integrated model of [START_REF] Antuono | Beyond Boussinesq-type equations: semi-integrated models for coastal dynamics[END_REF]. A more complete overview on these approaches was made by [START_REF] Kirby | Boussinesq models and their application to coastal processes across a wide range of scales[END_REF].

Although the turbulence is maybe the most striking phenomenon in a breaking wave, it is rarely taken into account directly in the models. In the eddy-viscosity models the turbulence is modelled by a turbulent-viscosity hypothesis but it is not resolved. As highlighted by [START_REF] Nadaoka | Shallow-water turbulence modeling and horizontal largeeddy computation of river flow[END_REF], the turbulence in shallow-water flows has a double-structural and strongly nonisotropic character. The double structure lies in the coexistence of a three-dimensional (3D) turbulence with length scales less than the water depth and horizontal two-dimensional (2D) eddies (with a vertical vorticity) with much larger length scales. Moreover the shallow-water turbulence may show an inverse cascade of energy or backscatter i.e. an energy transfer from the 3D turbulence toward the 2D eddies [START_REF] Nadaoka | Shallow-water turbulence modeling and horizontal largeeddy computation of river flow[END_REF][START_REF] Hinterberger | Three-dimensional and depth-averaged large-eddy simulation of some shallow water flows[END_REF]. With an eddy-viscosity hypothesis, even if calculated with a turbulent kinetic energy and an additional transport equation, not only the anisotropic character of the turbulence cannot be described but also the backscatter cannot be captured by the model.

In the Part I of this work [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF]) a filtering approach was implemented with a cutoff frequency in the inertial subrange. It follows that only the residual small-scale turbulence is modelled through a turbulent-viscosity hypothesis while the 3D subdepth turbulence is resolved. After depth-averaging the filtered equations over the depth, this subdepth turbulence was taken into account in the model by a new quantity called enstrophy. The isotropic character of the small-scale turbulence, the equality of the dissipation of the mean residual kinetic energy and its rate of production [START_REF] Lilly | The representation of small-scale turbulence in numerical simulation experiments[END_REF] and the large validity of the energy cascade hypothesis in the inertial subrange give a much greater validity to the turbulent-viscosity hypothesis for the residual turbulence with a cutoff in the inertial subrange. The existence of an explicit quantity for the subdepth large-scale turbulence is an advantage over previous approaches to model the breaking phenomenon and to describe the breaking waves. In particular, whereas finding and implementing a suitable breaking criterion is a laborious task for the models lacking a quantity describing explicitly the turbulence, the approach of [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF] leads to an easier breaking criterion based on the value of a variable of the model and, in the favourable cases, does not need a breaking criterion at all.

In the present paper the one-dimensional model of [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF] is extended to a two-dimensional model describing three-dimensional flows. Phenomena such as run-up and run-down are studied as well. In a 2D context the anisotropic character of turbulence is an important feature to model. The optimization techniques of the dispersive properties are extended to the new model in order to simulate accurately the wave propagation in the shoaling zone. An important goal is to be able to add more physical effects to the model, namely turbulence effects, without increasing significantly the complexity of the numerical resolution. Accordingly a suitable numerical scheme is extended to this model. Further the empirical laws determining the values of the model parameters are validated in a wide range of physical situations in order to give a real predictive character to the model.

The two-dimensional equations of the model are derived and discussed in §2 and the empirical laws chosen to model the eddy viscosity and the dissipation are presented in §3. The numerical implementation is explained in §4 with a formulation of the equations more suited to the numerical resolution and improving the dispersive properties. A discussion on the breaking criterion is given in §5. The numerical results on several test cases with the comparison with experimental results of the literature are presented in §6.

2 Two-dimensional depth-averaged filtered equations

Three-dimensional filtered equations

The Navier-Stokes equations of an incompressible fluid of density ρ and kinematic viscosity ν are filtered in the same manner as for the large-eddy simulation approach with a cutoff frequency in the inertial subrange. The details of this filtering approach are given in Part I of this work [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF]. The velocity field v is decomposed as v = v + v r where v is the filtered velocity field and v r is the residual (or subgrid) velocity field. The residual stress tensor is modelled by a turbulent viscosity hypothesis and the residual kinetic energy is absorbed into a modified pressure p (see for example Pope 2000 for more explanations on this approach). The filtered continuity equation is div

v = 0 (1)
and the filtered momentum equation can be written

∂v ∂t + div (v ⊗ v) = g - 1 ρ grad p + div 2 ν T D + ν∆v (2)
where g is the acceleration of the gravity, ν T is a turbulent viscosity and D is the filtered strain rate tensor D = 1 2 grad v + (grad v) T .

(3)

The kinetic energy of the filtered velocity field e f = v • v/2 satisfies the equation

∂e f ∂t + div e f v + pv ρ -2 (ν + ν T ) v • D + e p v = -ε f -P r (4)
where e p is defined by g = -grad e p , ε f = 2νD : D is the viscous dissipation in the filtered motions (the colon denotes the double dot product) and P r = 2ν T D : D is the energy transfer from the filtered motions towards the residual motions. At high Reynolds numbers, the term ε f is negligible. Denoting by ε the dissipation of the turbulent kinetic energy, a result due to [START_REF] Lilly | The representation of small-scale turbulence in numerical simulation experiments[END_REF] allows to write the equality of the dissipation of the mean residual kinetic energy and its rate of production (see also [START_REF] Pope | Turbulent Flows[END_REF][START_REF] Higgins | Energy dissipation in large-eddy simulation: dependence on flow structure and effects of eigenvector alignments[END_REF]) and therefore P r ε (the brackets denote the Reynolds averaging).

The problem is a three-dimensional flow over a variable bottom. The components of the filtered velocity field v are u and v in the horizontal directions Ox and Oy respectively and w in the vertical direction Oz. It is convenient to define the two-dimensional filtered velocity field in the horizontal plane by u = [u, v] T . The elevation of the bottom and of the free surface with respect to a horizontal datum are denoted by b(x, y) and Z(x, y, t) respectively. The water depth is h(x, y, t) = Z(x, y, t) -b(x, y). The still water depth is h 0 (x, y) and the water elevation is η(x, y, t) = h(x, y, t) -h 0 (x, y). These notations are depicted in Figure 1. In the following, the operators gradient and divergence are related to a two-dimensional space (Oxy) unless noted otherwise. The conventions used for tensor calculus are given in Appendix A.

The no-penetration boundary condition at the bottom and the kinematic boundary condition at the free surface can be written respectively

w(b) = u(b) • grad b (5) and w(Z) = ∂Z ∂t + u(Z) • grad Z. (6)
Additionally the dynamic boundary condition at the free surface is (σ • n)(Z) = 0 where σ is the Cauchy stress tensor including the turbulent viscosity effect and where n is the unit normal vector at the free surface.

Scaling

According to the shallow water hypothesis, which is here assumed to be valid, there is a small parameter µ = h * 0 /L 1 where h * 0 is a reference value of the still-water depth and where L is a characteristic lengthscale in the horizontal plane. A classical scaling is used to write the equations in dimensionless form (see [START_REF] Antuono | Beyond Boussinesq-type equations: semi-integrated models for coastal dynamics[END_REF]. The dimensionless quantities are denoted by a tilde symbol.

x = x L ỹ = y L z = z h * 0 t = µt g h * 0 h = h h * 0 p = p ρgh * 0 ũ = u gh * 0 ṽ = v gh * 0 w = w µ gh * 0 ũ = u gh * 0 b = b h * 0 Z = Z h * 0 (7)
The viscous stress tensor is defined by τ = 2ρνD = τ ij e i ⊗ e j where the vectors e i are the unit vectors (with Einstein notation and the indexes 1, 2 and 3 for x, y and z respectively).

If the Reynolds number defined by Re = h * 0 gh * 0 /ν is high enough, the viscous terms are negligible (see for example [START_REF] Antuono | Beyond Boussinesq-type equations: semi-integrated models for coastal dynamics[END_REF]. With the hypothesis of a weaklyturbulent flow (see below), these terms are negligible if Re = O(µ -3 ) [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF]) and this will be assumed in the following.

The same scaling as in [START_REF] Antuono | Beyond Boussinesq-type equations: semi-integrated models for coastal dynamics[END_REF] is used for the eddy viscosity, i.e. and that the dimensionless momentum equation in the horizontal plane can be written

∂ ũ ∂ t + div (ũ ⊗ ũ) + ∂ w ũ ∂ z + grad p = µ 2 div a + ∂a 3 ∂ z . (12)
In this expression, the two-dimensional vector a 3 = Ãr xz e x + Ãr yz e y is introduced as well as the two-dimensional tensor a = 2ν T s where

s = 1 2 grad ũ + (grad ũ) T . ( 13
)
We have

a 3 = ν T ∂ ũ ∂ z + O(µ 2 ). ( 14
)
The dimensionless momentum equation in the vertical direction Oz is

µ 2 ∂ w ∂ t + ∂ ũ w ∂ x + ∂ṽ w ∂ ỹ + ∂ w2 ∂ z = -1 - ∂ p ∂ z + µ 2 A r xz ∂ x + µ 2 A r yz ∂ ỹ + µ 2 A r zz ∂ z . (15)
The no-penetration boundary conditions writes w(b) = ũ(b) • grad b while the kinematic boundary condition is

w(Z) = ∂ h ∂ t + ũ(Z) • grad Z. ( 16
)
The dynamic boundary condition at the free surface becomes in dimensionless form

a 3 (Z) + p(Z) νT grad Z -µ 2 a(Z) • grad Z = 0, ( 17
) p(Z) + µ 2 νT a 3 (Z) • grad Z -µ 2 A r zz (Z) = 0. ( 18
)
All capillary effects are neglected. In the following the deviatoric residual stress tensor is neglected on the bottom (A r (b) 0) where there is a free-slip condition. [START_REF] Veeramony | The flow in surf-zone waves[END_REF] similarly neglected the shear stress on the bottom.

Averaged mass and momentum equations

The equations ( 11), ( 12) and ( 15) are averaged over the depth taking into account the boundary conditions. The averaged quantity corresponding to any quantity X is defined as

X = 1 h Z(x,t) b(x) X dz. ( 19
)
The two-dimensional average velocity vector is denoted by U = u and its components in the Ox and Oy directions are denoted by U and V respectively. The filtered horizontal velocity vector is decomposed as the sum of its averaged value and of a deviation u representing the large-scale turbulence and the shearing effects i.e.

u(x, y, z, t) = U (x, y, t) + u (x, y, z, t).

(20)

In the same way as in [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF], the flow is supposed to be weakly turbulent (see [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF]. This means that u = O(µ) and that in dimensionless form ũ = Ũ +µũ . This hypothesis is not very restrictive in practice since strong hydraulic jumps (with an upstream Froude number as high as 16 and consequently a very strong turbulence) can be correctly modelled this way [START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF]).

In the following the tilde symbols for dimensionless quantities are dropped to lighten the notations. After averaging over the depth, the mass equation becomes

∂h ∂t + div (hU ) = 0 (21)
With the kinematic boundary condition at the free surface, the averaging over the depth of the three first terms of the momentum equation ( 12) yields

Z b ∂u ∂t + div (u ⊗ u) + ∂wu ∂z dz = ∂hU ∂t + div Z b u ⊗ u dz. ( 22
)
The approach of [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF] is followed here for the treatment of the integral of u ⊗ u.

The enstrophy tensor is defined as

ϕ = 1 h 3 Z b u ⊗ u dz. (23)
Since, by definition, u = 0, we can write u ⊗ u = U ⊗ U + µ 2 h 2 ϕ. This gives

Z b ∂u ∂t + div (u ⊗ u) + ∂wu ∂z dz = ∂hU ∂t + div hU ⊗ U + µ 2 h 3 ϕ . ( 24
)
The equation ( 15) is needed to calculate the pressure term. First the material derivative of h is defined as ḣ

= ∂h ∂t + U • grad h. ( 25
)
Second an expression of the vertical velocity w can be obtained from the continuity equation. We get

w = (z -b) ḣ h + U • grad b + O(µ). ( 26
)
Then the left-hand part of equation ( 15) can be written

∂w ∂t + ∂uw ∂x + ∂vw ∂y + ∂w 2 ∂z = (z -b) ḧ h + D Dt (U • grad b) + O(µ) (27)
where the material derivative has the same meaning as the dot in (25) i.e. DX/Dt = ∂X/∂t + U • grad X for any scalar quantity X. This result is used to calculate the part of the pressure term that does not depend on the eddy viscosity. Using the boundary conditions, we find

Z b p dz = h 2 2 + µ 2 h 2 ḧ 3 + µ 2 h 2 2 D Dt (U • grad b) + µ 2 F 1 (28)
where F 1 is the integral of viscous terms

F 1 = Z b dz ∂ ∂x z Z A r xz dz + Z b dz ∂ ∂y z Z A r yz dz + Z b A r zz dz. ( 29
)
As in Part I of this work [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF], the turbulent viscosity is assumed to be uniform over the depth. It follows from this assumption and from the hypothesis of a weakly turbulent flow that

Z b dz ∂ ∂x z Z A r xz dz = O(µ) and Z b dz ∂ ∂y z Z A r yz dz = O(µ). ( 30
)
These two integrals can therefore be neglected. Using again the boundary conditions, this leads to the following expression for the averaged momentum balance equation

∂hU ∂t + div hU ⊗ U + µ 2 h 3 ϕ + grad h 2 2 + µ 2 h 2 ḧ 3 + µ 2 Π = -p(b)grad b + µ 2 div Z b a dz -µ 2 grad Z b A r zz dz + O(µ 3 ) (31)
where

Π = h 2 2 D Dt (U • grad b) (32)
and p(b) = h + µ 2 h ḧ/2 + µ 2 2Π /h + O(µ 3 ). The continuity equation implies that A r zz = -A r xx -A r yy . The calculation of the last integrals leads to

∂hU ∂t + div hU ⊗ U + µ 2 h 3 ϕ + grad h 2 2 + µ 2 h 2 ḧ 3 + µ 2 Π = -h grad b + µ 2 div (hA) -µ 2 f + O(µ 3 ). ( 33
)
The expression of the quantity f is

f = h ḧ 2 + D Dt (U • grad b) grad b. ( 34
)
It encompasses terms due to the variable bottom of O(µ 2 ). The tensor A satisfies the relation

A = 2ν T S + 2ν T (div U ) I ( 35
)
where I is the identity tensor (two-dimensional) and where

S = 1 2 grad U + (grad U ) T . ( 36
)
This tensor A acts like a viscous stress tensor of a compressible fluid. The averaged mass equation ( 21) of the model is analogous to the mass conservation equation of a compressible fluid, the depth h being analogous to the density. The relation ( 35) is then analogous to the constitutive equation of a Newtonian compressible fluid i.e. τ = 2µD+η(div v)I 3 where D is the strain rate tensor, I 3 the three-dimensional identity tensor, v is the velocity field (the operator divergence is here three-dimensional), µ and η being the first and second viscosities respectively.

In the case of the model, the first viscosity is equal to the turbulent viscosity ν T and the second viscosity is equal to 2ν T . Note that S is not a deviator since tr S = div U = 0. For compressible fluids, it is more convenient to define the deviatoric tensor S 0 = S -(div U )I/2 (this gives tr S 0 = 0 since tr I = 2 in a two-dimensional space). The tensor A can then be written

A = 2ν T S 0 + ζ (div U ) I (37)
where ζ = 3ν T is the sum of the first and second viscosities and is the volume viscosity of the model (in a three-dimensional space, the volume viscosity is equal to η + 2µ/3). The first viscosity and the volume viscosity are both positive which is in accordance with the second law of thermodynamics. Note that the Stokes' hypothesis, according to which the volume viscosity is equal to zero, is not satisfied in the model. The gradient term in the left-hand side of (33) includes dispersive effects originating from a non-hydrostatic correction to the pressure. It should be noted that no smallness assumption was made on the nonlinearity which implies that the model is fully nonlinear. Its dispersive properties are identical to those of the Green-Naghdi equations. An improvement of these properties is proposed in §4.

The model is anisotropic due to the enstrophy tensor ϕ. It follows from its definition (23) that the enstrophy tensor is symmetrical. However it is not an isotropic tensor and it has three independant components denoted by ϕ 11 , ϕ 12 and ϕ 22 which are defined by ϕ = ϕ 11 e x ⊗ e x + ϕ 12 e x ⊗ e y + ϕ 12 e y ⊗ e x + ϕ 22 e y ⊗ e y . This tensor represents the large-scale turbulence and the shearing effects. The large-scale turbulence of the energy-containing range has an anisotropic character which is thus taken into account in the model through the anisotropic tensor ϕ.

This is one of the main advantages of this approach in comparison with the classical approaches modelling all turbulence with an eddy viscosity. A viscosity hypothesis is a valid model for turbulence if this one is reasonably isotropic and not too far from the equilibrium between production and dissipation. These two conditions are questionable in the case of the turbulence of a breaking wave. Consequently the classical eddy viscosity approaches for breaking waves misses the anisotropic effects of the large-scale turbulence [START_REF] Nadaoka | Shallow-water turbulence modeling and horizontal largeeddy computation of river flow[END_REF]. Furthermore there is no backscatter with an eddy-viscosity model and consequently no energy transfer from the horizontal three-dimensional eddies towards the vertical two-dimensional eddies and yet this transfer can happen (see for example [START_REF] Nadaoka | Shallow-water turbulence modeling and horizontal largeeddy computation of river flow[END_REF][START_REF] Hinterberger | Three-dimensional and depth-averaged large-eddy simulation of some shallow water flows[END_REF].

With our depth-averaged large-eddy simulation (LES) approach with a cutoff frequency in the inertial subrange, only the small-scale turbulence is modelled by a turbulent viscosity hypothesis while the anisotropic large-scale turbulence of the energy-containing range is resolved. Both the isotropic and equilibrium conditions are well satisfied for the small-scale turbulence (see for example Kolmogorov's hypotheses, Kolmogorov 1941) and the absence of backscatter from the small scales towards the large scales is a very classical view in accordance with the energy cascade of [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF] which was confirmed experimentally in the case of breaking waves by [START_REF] Hattori | Experimental study on turbulence structures under breaking waves[END_REF]. The introduction of the eddy viscosity has thus better physical justifications than for the classical eddy-viscosity models. The anisotropic character of our equations has a physical basis and on the whole the resolution of the large-scale turbulence gives our model a richer physical content.

The model features six scalar unknowns which are the water depth h, the components U and V of the average velocity field in the Ox and Oy directions respectively and the three components of the enstrophy tensor. The mass (21) and momentum (33) equations provide three scalar equations. Three more equations are thus needed to close the system. In the onedimensional case [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF] the closure was obtained with the energy equation. In this two-dimensional case, the energy equation gives only one scalar equation and it is not sufficient. In the one-dimensional (1D) case the energy equation could be replaced with the enstrophy equation which was one scalar equation since the 1D enstrophy is a scalar. In the two-dimensional (2D) case the enstrophy tensor equation is used to provide the three remaining scalar equations.

Enstrophy tensor equation

The tilde symbols are also dropped in this section. The first step to derive the enstrophy tensor equation is to form the tensor product u ⊗ (12) + (12) ⊗ u. Taking into account the boundary conditions, this leads to

∂u ⊗ u ∂t + div (u ⊗ u ⊗ u) + ∂ ∂z (wu ⊗ u) + u ⊗ grad p + grad p ⊗ u = µ 2 u ⊗ div a + µ 2 diva ⊗ u + ν T u ⊗ ∂ 2 u ∂z 2 + ν T ∂ 2 u ∂z 2 ⊗ u. ( 38
)
The second step is to average this equation over the depth, again with the boundary conditions. This procedure yields

∂ ∂t hU ⊗ U + µ 2 h 3 ϕ + div hU ⊗ U ⊗ U + µ 2 U ⊗ h 3 ϕ + µ 2 h 3 ϕ ⊗ U + µ 2 div h 3 ϕ ⊗ U + µ 2 h 3 ϕ • (grad U ) T + U ⊗ grad P + grad P ⊗ U = -p(b) (U ⊗ grad b + grad b ⊗ U ) + µ 2 div (U ⊗ hA) + µ 2 [div (U ⊗ hA)] T -µ 2 h P r + O(µ 3 ) ( 39
)
where

P = h 2 2 + µ 2 h 2 ḧ 3 + µ 2 Π . ( 40
)
The tensor P r includes all dissipative effects and its expression is

P r = grad u • a + a • (grad u) T + 2ν T ∂u ∂z ⊗ ∂u ∂z + 2 (tr a) s + O(µ 2 ). ( 41
)
Note that tr P r = 2P r since the energy equation ( 4) is half the trace of the tensor equation ( 38). This tensor corresponds to a transfer from the filtered scales towards the residual scales in the same way as P r is an energy transfer from the filtered scales towards the residual scales.

Extending the result of [START_REF] Lilly | The representation of small-scale turbulence in numerical simulation experiments[END_REF] according to which P r ε, the average tensor P r is almost equal to the dissipation tensor ε which corresponds to the dissipation of the residualstress tensor σ r = -ρ(v ⊗ v -v ⊗v). The estimation P r ε made in [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF] is extended to the corresponding tensors as

P r ε . ( 42
)
The third step is to form the tensor product U ⊗ (33) + (33) ⊗ U . This gives

∂hU ⊗ U ∂t + div hU ⊗ U ⊗ U + µ 2 U ⊗ h 3 ϕ -µ 2 grad U • h 3 ϕ + µ 2 div h 3 ϕ ⊗ U + U ⊗ grad P + grad P ⊗ U = -p(b) (U ⊗ grad b + grad b ⊗ U ) + µ 2 U ⊗ div (hA) + µ 2 div (hA) ⊗ U + O(µ 3 ). ( 43
)
This equation ( 43) is finally subtracted from equation ( 39) to yield the evolution equation of the enstrophy tensor

∂hϕ ∂t + div (hϕ ⊗ U ) -2h (div U ) ϕ + grad U • hϕ + hϕ • (grad U ) T = grad U • A h + A h • (grad U ) T - ε h + O(µ). ( 44
)
It is noteworthy that this equation includes no dispersive term and no term depending on the variable bottom. Further there is no second-order nor third-order derivatives which means that the numerical resolution of this equation is a priori not especially difficult if the dissipative term is known. This equation is similar to the Reynolds-stress equation with the following differences: First there is no term involving the third-order tensor u ⊗ u ⊗ u which is of O(µ 3 ) and thus negligible. This is the main interest of Teshukov's hypothesis of a weakly turbulent flow [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF]. Second there is no velocity-pressure-gradient tensor (or no pressure-rate-of-strain tensor) because the depth-averaging procedure is based on a decomposition of the horizontal filtered velocity field but not on a decomposition of the pressure. The pressure is explicitly and consistently expressed from the momentum equation in the Oz-direction. There is thus no analogous to the fluctuating pressure field. These two differences constitute a huge simplification. Third there is no diffusive term. Lastly the production tensor has two parts. The first part,

-grad U • hϕ -hϕ • (grad U ) T , (45)
is relative to the large-scale turbulence. It originates from the depth-averaging of the filtered velocity field and thus from the part of the turbulence which is resolved. It is anisotropic. The second part,

grad U • A h + A h • (grad U ) T , (46)
is due to the residual small-scale turbulence. Modelling the residual-stress tensor with an eddyviscosity hypothesis gives this term the structure of a viscous production, although of a compressible type since div U = 0. The term -2h (div U ) ϕ is due to the equation being written for the evolution of the tensor ϕ instead of the evolution of the tensor u ⊗ u = h 2 ϕ. The left-hand part of equation ( 44) was obtained by [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF] who included no dissipation.

The compressible viscous production is new in this kind of approach.

Energy

Taking half the trace of equation ( 39) gives an energy balance equation for the model. This energy equation can also be derived by averaging over the depth the equation ( 4). This equation writes

∂ ∂t h e + µ 2 e + div hU e + µ 2 e + P U + µ 2 U • h 3 ϕ = -hU • grad b + µ 2 div (hU • A) -µ 2 h 2 tr ε + O(µ 3 ) (47)
where

e = U • U 2 + µ 2 h 2 2 tr ϕ + h 2 + µ 2 ḣ2 6 (48) and e = ḣ 2 (U • grad b) + 1 2 (U • grad b) 2 . ( 49
)
The energy of the system is e + µ 2 e where e encompasses terms of O(µ 2 ) due to the variable bottom. The term µ 2 h 2 tr ϕ/2 is a turbulent energy for the system. It includes only the largescale turbulence. The equation of tr ϕ can be found by taking the trace of equation ( 44). Alternatively it can be derived from the energy equation ( 47) and from the mass and momentum equations ( 21) and ( 33). The equation for the trace of the enstrophy tensor is thus equivalent to the energy equation. It can be written

∂htr ϕ ∂t + div (hU tr ϕ) -2h (div U ) tr ϕ + 2hϕ : grad U = 2 h A : grad U - tr ε h + O(µ). ( 50
)
The term in A : grad U is positive and results in a creation of turbulent energy since A :

grad U = 2ν T S 0 2 + 3ν T (div U ) 2 0 where S 0 = (S 0 : S 0 ) 1/2 .

Vorticity of the mean flow

Let's denote by J the determinant of the enstrophy tensor ϕ. The Cauchy inequalities imply that J 0. If the dispersive, viscous and dissipative terms are removed, the remaining system is hyperbolic [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF]) if J = 0 [START_REF] Richard | Élaboration d'un modèle d'écoulements turbulents en faible profondeur. Application au ressaut hydraulique et aux trains de rouleaux[END_REF]). If the system is restricted to its hyperbolic part, the determinant of the enstrophy tensor satisfies the equation [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF])

DJ Dt = 2Jdiv U (51)
which implies that h 2 J is conserved along the trajectories of the mean flow and can be interpreted as an entropy of the system [START_REF] Gavrilyuk | Geometrical evolution of the Reynolds stress tensor[END_REF]. Also, in the hyperbolic system, if J > 0 at a time t = 0, then J > 0 at any time t > 0. A geometrical interpretation of the evolution equation of the tensor h 2 ϕ was given in [START_REF] Gavrilyuk | Geometrical evolution of the Reynolds stress tensor[END_REF] who showed that the eigenvectors of this tensor undergo a rotation similar to a rigid body and form a natural moving frame whose evolution is determined by the mean rate of the deformation tensor. This interpretation can be further specified. Assuming that J = 0, the tensor ϕ is invertible. The evolution equation of the inverse tensor ϕ -1 can be deduced from equation ( 44). Restricting ourselves again to the hyperbolic part of the equations, this equation can be written

Dϕ -1 Dt = -2 (div U ) ϕ -1 + ϕ -1 • S + S • ϕ -1 + ϕ -1 • Ω -Ω • ϕ -1 (52)
where the tensor Ω is the mean rotation-rate tensor which is the antisymmetric part of grad U .

For comparison, we define another tensor that satisfies a similar equation. First we consider two orthogonal infinitesimal vectors dx 1 = dx 1 e 1 and dx 2 = dx 2 e 2 attached to a point of the fluid and transported by the mean flow (dx 1 > 0 and dx 2 > 0). They satisfy the equation D(dx i )/Dt = grad U • dx i , (i ∈ {1, 2}). We define the diagonal tensor P as P = dx 1 e 1 ⊗ e 1 + dx 2 e 2 ⊗ e 2 . The vectors e 1 and e 2 are eigenvectors of P and dx 1 and dx 2 are eigenvalues. The tensor Q defined by Q = (P • P) -1 can be represented by an ellipse whose semi-major axis is dx 1 and semi-minor axis is dx 2 (assuming dx 1 > dx 2 ). The following equality pertaining to the tensor Q can be derived

DQ Dt = -2 (divU ) Q + Q • S + S • Q -Q • Ω + Ω • Q. ( 53
)
The comparison between ( 52) and ( 53) shows that these equations differ only in that the mean rotation-rate tensor Ω in ( 53) is replaced by -Ω in (52). This means that the inverse enstrophy tensor is convected by the mean flow, deformed by the mean strain rate tensor and rotated by the opposite of the mean rotation-rate tensor. Moreover the tensor ϕ -1 can be represented by an ellipse whose area can be written A = π √ J. The equation ( 51) implies that

1 A DA Dt = divU . ( 54
)
The rate of change of the area of the ellipse is the rate of change of area due to the mean flow (in a three-dimensional space div U would be a rate of change of volume). The area of the ellipse associated to Q satisfies the same equation. It follows that the dilatation, convection and deformation of this ellipse is the same as an infinitesimal ellipse attached to material particles moving with the mean flow but that this ellipse undergoes a rotation with an angular velocity which is exactly the opposite of the local angular velocity of the mean flow.

The vorticity of the mean flow could be defined as Ω = rot U but, as the mean flow is two-dimensional, there is only one non-zero component which is in the Oz direction and it is better to define the vorticity as

Ω = ε : grad U . ( 55
)
In this expression, ε is the two-dimensional pseudotensor of Levi-Civita which can be written ε = ε ij e i ⊗ e j where ε ij is the two-dimensional symbol of Levi-Civita defined by

ε ij =    +1 if (i , j) = (1 , 2) -1 if (i , j) = (2 , 1) 0 if i = j. (56)
The mean rotation-rate tensor can be written Ω = -(Ω/2)e 1 ⊗ e 2 + (Ω/2)e 2 ⊗ e 1 . Thus it seems that the difference between the equations ( 52) and ( 53) is due to the vorticity of the mean flow and that the tensors ϕ -1 and Q satisfy the same equation if the mean flow is irrotational [START_REF] Debieve | Evolution of the Reynolds stress tensor in a shock wave-turbulence interaction[END_REF]. However an evolution equation of the vorticity of the mean flow can be derived from the equation ( 33). This vorticity equation can be written

∂Ω ∂t + U • grad Ω = -Ω div U + µ 2 ε : grad 1 h div 2ν T hS 0 -h 3 ϕ D - µ 2 h 2 ε : grad 3ν T h div U - h 3 2 tr ϕ ⊗ grad h + µ 2 ε : grad ḧ ⊗ grad h 3 + b 2 + µ 2 ε : grad D Dt (U • grad b) ⊗ grad h 2 + b . (57)
The decomposition of the enstrophy tensor as the sum of its isotropic (or spherical) part and of its anisotropic (or deviatoric) part is written ϕ = (trϕ)I/2 + ϕ D where ϕ D is the deviator of the enstrophy tensor. The compressible character of the equations entails the presence of source terms in this vorticity equation such as the term -Ω div U in the right-hand side of the equation. There are also baroclinic terms which in the three-dimensional usual case would be written -(1/ρ 2 ) grad ρ × grad p and which here take the form -(1/h 2 ) ε : grad P ⊗ grad h where P is some scalar field (h is here analogous to the density ρ). In particular there is a baroclinic term due to the turbulent volume viscosity ζ = 3ν T and to the trace of the enstrophy tensor. There are also source terms akin to baroclinic terms due to the dispersion and to the variable bottom. And finally there are source terms due to the first turbulent viscosity and to the deviator of the enstrophy tensor. This equation shows the big difference between the tensors ϕ -1 and Q. The equation ( 53) is only a geometrical description of the variation of Q while it is transported by the mean flow. The tensor Q has no influence on the mean flow and is only passively transported. On the contrary, the enstrophy tensor is not only transported and modified by the mean flow, it modifies the mean flow since ϕ appears in the momentum equation ( 33). This means that the equation ( 44) is a physical equation and not a mere geometrical equation such as (53) in spite of the similarity.

Even if the mean rotation-rate tensor has no influence on the trace of the enstrophy tensor nor on the energy, because the terms involving Ω in (52) are traceless, it is inconsistent to neglect these terms in the enstrophy tensor equation, assuming that the mean flow is irrotational or weakly rotational, because the enstrophy tensor is a source in the vorticity equation, notably by its deviatoric part. It is also not possible to calculate the enstrophy tensor by replacing the equation ( 44) by an equation of the kind of ( 53) with a hypothesis of a weak vorticity because this one is not tenable given that the enstrophy can create vorticity. Moreover the equation ( 53) has a completely different meaning, is not hyperbolic and has wholly different characteristics.

The presence of the enstrophy tensor, including its deviator, in the vorticity equation ( 57) shows that transfers can happen between the three-dimensional horizontal eddies modelled by the enstrophy tensor and the two-dimensional vertical eddies represented by the vorticity of the mean flow. Our depth-averaged LES approach with a cutoff in the inertial subrange follows the energy cascade from the energy-containing range towards the dissipation range but can describe transfers from the large-scale turbulence with a scale of O(h) or smaller (in the water depth) towards the vertical eddies of the average flow with an even bigger scale.

Modelling the eddy viscosity and the dissipation tensor

To close the model, two quantities remain to be specified, namely the turbulent viscosity ν T and the averaged dissipation tensor ε . Empirical laws are proposed to determine these two quantities. A correct model of the dissipation should preserve the positivity of the enstrophy tensor and thus of its determinant J and of its trace. The determinant can be written J = [(tr ϕ) 2 -ϕ:ϕ]/2. Ignoring all terms except the dissipation in the enstrophy equation ( 44

A simple way to preserve the positivity of both tr ϕ and J is to have Dtr ϕ/Dt = -tr ϕf 1 and DJ/Dt = -Jf 2 where f 1 and f 2 are positive functions depending on the variables of the mean flow. As the equations of the model must satisfy the Galilean invariance, these functions should not depend on U . Moreover, they should be invariant in every coordinate systems. This implies that they can depend on h and on the two invariants of the tensor, tr ϕ and J. The simplest way is to choose ε = ϕf (h, tr ϕ) where the function f depends only on h and tr ϕ and is obtained by a dimensional analysis. This gives f (h, tr ϕ) = C r h 2 √ tr ϕ where C r is a dimensionless quantity. The model for the average dissipation is thus

ε = C r h 2 ϕ √ tr ϕ. ( 60
)
With this choice the one-dimensional case studied in [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF] is recovered as a particular case. Further, since the turbulent kinetic energy of the model is (h 2 tr ϕ)/2, the chosen model has the same form as the model proposed by [START_REF] Rotta | Statistische Theorie nichthomogener Turbulenz[END_REF] if the dissipation and the turbulent kinetic energy are related with a mixing length proportional to the fluid depth.

The model for the turbulent viscosity is obtained by a similar approach. The eddy viscosity must be the same in every Galilean reference frames and in every coordinate systems. It is simpler to assume that ν T depends on h and tr ϕ but not on J. A dimensional analysis yields

ν T = C p h 2 √
tr ϕ where the dimensionless quantity C p can be interpreted as the inverse of a turbulent Reynolds number R leading to

ν T = h 2 R √ tr ϕ. (61)
The one-dimensional case [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF]) is also recovered as a particular case. This expression of the eddy viscosity based on the turbulent kinetic energy and on a mixing length, which is here proportional to the water depth, is a very classical one [START_REF] Kolmogorov | The equations of turbulent motion in an incompressible fluid[END_REF].

Numerical implementation

In the absence of enstrophy, the present model reduces to the classical two-dimensional Green-Naghdi equations. The introduction of this new variable entails a modified pressure law and three supplementary transport equations that can be naturally injected in the hyperbolic part of the system. As a consequence, from a general point of view, this model can be numerically treated in a straightforward manner on the basis of any existing numerical approach dedicated to the Green-Naghdi equations. Naturally, this is even simpler in the one-dimensional case, since the enstrophy tensor equation reduces to a scalar equation. On this basis, preliminary results were obtained by [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF] using an extension of the strategy proposed in Le [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF] to mild slope topographies, allowing to establish relevant empirical laws based on the study of solitary waves. In the present work we aim to pursue the efforts towards describing realistic situations, considering a two-dimensional approach on general bathymetries.

In addition, and notably for the analysis of periodic waves, we have to think about how to improve the dispersive properties of the proposed model. Before getting through the numerical implementation, we briefly discuss how to address such an issue.

Constant diagonal formulation

In view of the targeted numerical simulations, we first propose to rewrite the model in an asymptotically equivalent form, allowing both to improve the dispersive properties and provide a gain in terms of computational cost. To achieve this outcome, we exploit the works proposed by [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], followed later by [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF], mainly for purposes of computational efficiency. The present section merely consists in verifying that the current model enters into the appropriate formalism. Consequently we shall not discuss all technical aspects, and refer to the aforementioned references for more details.

As it is the case when dealing with Boussinesq-type models, the numerical challenge stands in the treatment of the discharge equation, due to the presence of high-order derivatives (up to third order) and second-order instationnary terms appearing through the material derivative. This last point has heavy practical consequences since it results in the inversion of a time-dependent elliptic operator which couples the velocity components.

The model is constituted by equations ( 21), ( 33) and ( 44). These equations, written with the components of U and ϕ, are collected for convenience' sake in Appendix B. As mentioned before, this model reduces to the standard Green-Naghdi equations with a modified pressure law, supplemented by the transport of the enstrophy tensor. We hence inherit the same technical difficulties as discussed just before, concentrated in the discharge evolution. As a consequence, in what follows we will only focus on the momentum equation.

As a first remark, adapting the notations used in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], we can express the average velocity equation (33) as

I + µ 2 T ∂U ∂t + U • grad U + grad (h + b) + µ 2 h div h 3 ϕ -hA + µ 2 Q 0 = O(µ 3 ). ( 62
)
In this expression I is the identity operator, T is a linear operator depending on h and b and defined by

T W = - 1 3h grad h 3 divW - h 2 divW grad b + 1 2h grad h 2 grad b • W + (grad b • W ) grad b (63)
for any vector field W and

Q 0 = - 1 3h grad h 3 div (U div U ) -2 (div U ) 2 + [U • grad (U • grad b)] grad b - h 2 div (U div U ) -2 (div U ) 2 grad b + 1 2h grad h 2 U • grad (U • grad b) . ( 64
)
Introducing the notation

Q 0 = Q 0 + 1 h div h 3 ϕ -hA , (65)
the previous equation becomes

I + µ 2 T ∂U ∂t + U • grad U + grad (h + b) + µ 2 Q 0 = O(µ 3 ) , (66)
and we can reproduce the arguments employed in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] to obtain the equation

I + µ 2 T ∂U ∂t + U • grad U + grad (h + b) + µ 2 Q 1 = O(µ 3 ) (67)
where

Q 1 = Q 0 -T (U • grad U ) . (68)
Denoting by I + µ 2 T -1 the inverse operator of I + µ 2 T , this equation can be written

∂U ∂t + U • grad U + I + µ 2 T -1 grad (h + b) + µ 2 Q 1 = O(µ 3 ) (69)
As highlighted by [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF], this formulation has two main advantages, namely it does not require the computation of third-order derivatives, and the presence of the operator I + µ 2 T makes the model stable with respect to high frequency perturbations. Nevertheless, a remaining major drawback is the necessity to invert this operator, since it is time-dependent and involves a coupling between the velocity components. Starting form the dimensionless mass equation and a velocity equation under the form (67), [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF] reached the following so-called "constant diagonal " formulation of the momentum equations, which reads, in the present context:

I + µ 2 αT ∂hU ∂t + div (hU ⊗ U ) + α -1 α h grad (h + b) + 1 α h grad (h + b) + µ 2 (hQ 1 + hQ 2 + Q 3 ) = O(µ 3 ), ( 70
)
where the operator T only depends on a given water depth at rest h 0 and has a time-independent diagonal structure since

T W = - 1 3 div h 3 0 grad W h 0 = - 1 3 ∂ ∂x j h 3 0 ∂ ∂x j W i h 3 0 e i (71)
for any vector field W = W i e i . The auxiliary quantities are defined as follows:

Q 2 = h grad ε • grad (h + b) • ε • grad h - 1 2h grad h 2 grad b • grad (h + b) + h 2 grad b ∆ (h + b) -[grad b • grad (h + b)] grad b + O(µ 2 ) (72)
and

Q 3 = S I + µ 2 αT -1 [h grad (h + b)] ( 73
)
where S is an operator defined as

S{W } = 1 6 grad h 2 -h 2 0 • grad W + 1 3 h 2 -h 2 0 ∆W - 1 6 ∆ h 2 -h 2 0 W (74)
for a generic smooth enough vector field W . The role of the constant α, initially introduced by [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], is precisely to improve the dispersive properties of the model (namely, to provide a better matching with respect to the linear Stokes theory). Following the works previously mentioned, we take α = 1.159 in our numerical experiments. Note that in the above formulation (70), the only difference with the original work is the expression of Q 1 (68), which through (65) also contains the conservative and viscous terms related to the enstrophy. Going back to variables with dimensions, the equation (70) becomes:

(I + αT) ∂hU ∂t + div (hU ⊗ U ) + α -1 α gh grad (h + b) + gh α grad (h + b) + hQ 1 + ghQ 2 + gQ 3 = 0. (75)
Finally, following what has been done in [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF], we reformulate the previous equation as:

∂hU ∂t + div hU ⊗ U + h 3 ϕ + gh grad (h + b) + D = 0 (76)
where

D = - gh α grad (h + b) -div h 3 ϕ + (I + αT) -1 gh α grad (h + b) + hQ 1 + ghQ 2 + gQ 3 . (77)

Pre-balanced formulation

In this work, numerical investigations will be based on the high-order discontinuous Galerkin approach developped by [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF]. This approach directly applies to (76) with a preliminary rewriting of the hyperbolic part in the prospect of well-balancing properties. This approach extends the 1D works introduced by Liang & Marche (2009) for the shallow water (Saint-Venant) equations, in the context of unstructured high-order discretizations. The key idea is to reformulate the hydrostatic pressure term as

gh grad (h + b) = g 2 grad(Z 2 -2Zb) + gZ grad b . ( 78
)
Then we denote the vector solution W = (Z, hU T , hϕ 11 , hϕ 12 , hϕ 22 ) T , and recast (76) in the compact form:

∂W ∂t + div F(W , b) + D(W , b) = S(W , b), (79)
F(W , b) =        hU hU ⊗ U + h 3 ϕ + g 2 (Z 2 -2Zb)I hϕ 11 U hϕ 12 U hϕ 22 U        , D(W, b) =       0 D 0 0 0       . ( 80
)
The dispersive R 2 -valued term D is defined by ( 77), and the source term by

S(W, b) =       0 -gZ grad b R 11 + E 11 -C r hϕ 11 √ ϕ 11 + ϕ 22 R 12 + E 12 -C r hϕ 12 √ ϕ 11 + ϕ 22 R 22 + E 22 -C r hϕ 22 √ ϕ 11 + ϕ 22       , (81)
with

R
Note that the flux function F in (79) allows a clear decoupling between the hyperbolic and dispersive parts of the system.

Discrete formulation

We now turn to the numerical scheme, which consists in a direct extension of the works of [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF]. In consequence, we shall not go into detail here, and refer to this paper for an exhaustive description. The main principles of this numerical scheme are given in Appendix C. An example of a regular triangular mesh used in the numerical computations is given in Figure 2. Overall, the main numerical challenge stands in the computation of the components of D, which involves the resolution of a global linear system. As been said before, in our case this task is considerably alleviated since the elliptic operator appearing in ( 76) is time-independent and allows to decouple the velocity components evolution. The computation of this term follows the protocol described in [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF], based on the use of Local Discontinuous Galerkin fluxes [START_REF] Cockburn | The local discontinuous Galerkin method for timedependent convection-diffusion systems[END_REF] to build the first and second-order differential operators. We refer to the above papers for technical details. As concerns the enstrophy transport, the conservative terms are treated in the hyperbolic stage, according to (79), while the associated source terms are computed in a collocated framework with direct nodal products, in the same way as those of the momentum equations. Classically, advancing in time will be carried out by standard high-order Strong Stability Preserving Runge Kutta (SSP-RK) algorithms, following the original work.

Wave breaking

Preliminary results of [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF] showed that the activation of enstrophy at the beginning of the computations may impact the flow evolution before the breaking zone, leading to underestimations of the wave amplitude. As expected, the phenomenon is all the more pronounced as high nonlinearities are involved. As a consequence, for now, and in the philosophy of other existing strategies, this compels us to introduce a breaking criterion in order to activate the enstrophy only in eligible areas. Naturally, the different detection protocols mentioned in the introduction can be used to turn on viscosity terms and generate enstrophy only when needed. As for hybrid or ad-hoc viscosity methods, this entails the introduction of a discontinuity in the momentum equations in the neighbourhood of the breaking area. From a general viewpoint, this can be a source of numerical instabilities, especially if the discretization parameters are not chosen appropriately. In the present case, the viscous terms being expressed in terms of the enstrophy, they are subject to a regular growth localized in a very thin region surrounding the breaking point, allowing to introduce the turbulent effects in a smooth way. Further, this behaviour makes the model less sensitive with respect to the mask width in which viscous terms are activated.

Another detection strategy, proposed by [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF], consists in using the enstrophy itself to detect wave breaking. Indeed, this quantity is intrinsically linked to shear effects and can be used as a relevant tool to predict the development of turbulent structures. This leads to the introduction of a new quantity, referred to as virtual enstrophy, which quantifies the amount of enstrophy that the model is potentially able to produce. As in the original work, the virtual enstrophy is computed at each time step as the real enstrophy, but without any feedback on the evolution of the other variables. Then, breaking points can be identified as those where the virtual enstrophy is high enough, which also imposes the introduction of an appropriate threshold. Once this threshold is exceeded, viscous terms can be activated in the neighbourhood of the incriminated cells in order to light the effective enstrophy to handle wave breaking, in the same manner as for the detection protocols previously mentioned. Classically, a slope limiter is applied in these areas [START_REF] Cockburn | The local discontinuous Galerkin method for timedependent convection-diffusion systems[END_REF], notably to damp the brutal entropy variations. Explicit details on the threshold calibration and other breaking parameters will be given in the following section. Note that this strategy naturally implies a non-zero value as initial data for the enstrophy. Following the original work, this value will be taken to 10 -10 s -2 in our numerical tests.

6 Numerical results

Augmented solitary wave

As a preliminary test we validate the ability of the scheme to capture an analytical solution, in the absence of topography and without breaking effects. In a recent work, [START_REF] Richard | Modelling turbulence generation in solitary waves on shear shallow water flows[END_REF] showed the existence of solitary waves with non-trivial constant enstrophy profiles as exact solutions to the model. Considering a reference enstrophy ϕ 0 , the corresponding relative amplitude ã = a/h 0 , where the amplitude a of the wave is the maximum value of η (see Figure 1), is

ã = 1 2 φ0 -(1 + 4 φ0 ) + (1 + 4 φ0 ) 2 + 4 (F r 2 -1 -3 φ0 ) φ0 , ( 84
)
where F r stands for the Froude number and φ0 = h 0 ϕ 0 /g is the dimensionless enstrophy.

Inverting the previous equation we obtain:

F r = (1 + ã) [1 + φ0 (3 + ã)] , (85)
leading to the following generalized wave celerity:

c 0 = F r gh 0 = g (h 0 + a) + ϕ 0 (h 0 + a) (3h 0 + a). ( 86
)
The free surface elevation is given by the following formula:

η(x, t) = 2a F r 2 -1 -3 φ0 F r 2 -1 -(3 + a 2 ) φ0 + (F r 2 -1 -(3 -a 2 ) φ0 ) cosh[κ(x -c 0 t -x 0 )] , (87)
where

κ = 3 (F r 2 -1 -3 φ0 ) /F r 2 , (88)
and x 0 stands for the initial location of the solitary wave. Setting the transverse velocity V and the other components of the enstrophy tensor to zero, the exact solution is given by:

       h(x, t) = h 0 + η(x, t) U (x, t) = c 0 1 - h 0 h(x, t) ϕ 11 (x, t) = ϕ 0 (89)
The computational domain consists of a 200 m long rectangular channel, meshed with a regular triangulation of characteristic size ∆x = ∆y = 0.25 m (see Figure 2). The solitary wave is initially located at x 0 = 50 m. We set h 0 = 1 m and impose F r = 1.2, ϕ 0 = 0.2 s -2 , leading to an amplitude a = 0.128 m and the initial condition depicted in Figure 3. As highlighted by [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF], this gives a solution with a smaller amplitude than the classical solitary wave solution of the Green-Naghdi equations obtained with the same Froude number. We can observe on Figure 4 a comparison between analytical and numerical results at several propagation times, obtained with the third-order scheme, highlighting a very good resolution of the wave propagation. Note that the preservation of the initial enstrophy ϕ 0 has been numerically confirmed throughout the computation up to the machine error, independently from the mesh size or the polynomial degree in the approximation space.

Wave breaking and run-up of a solitary wave

We now turn to a classical 1D test case implying wave breaking and run-up, based on the experiment of [START_REF] Synolakis | The runup of solitary waves[END_REF]. The initial condition consists of an incident wave, obtained as an exact solution of the classical Green-Naghdi equations, propagating over a beach with

(t * = t(g/h 0 ) 1/2 ).
a constant bed slope s = 1/19.85. This benchmark is widely used to exhibit the ability to capture shoaling and breaking processes, with subsequent run-up and run-down phenomena (see [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Cienfuegos | Wave-breaking model for Boussinesqtype equations including roller effects in the mass conservation equation[END_REF][START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF][START_REF] Titov | Modeling of Breaking and Nonbreaking Long-Wave Evolution and Runup Using VTCS-2[END_REF][START_REF] Tonelli | Finite volume scheme for the solution of 2D extended Boussinesq equations in the surf zone[END_REF][START_REF] Zelt | The run-up of nonbreaking and breaking solitary waves[END_REF] for instance). The numerical set-up implies a solution centred at x 0 = 10 m in a 35 m long domain, regularly meshed with a space step ∆x = ∆y = 0.0625 m, as indicated in Figure 2. A 5 m long sponge layer has been added at the left boundary to stabilize the initiation of the wave propagation. Following the experiment, the reference depth is h * 0 = 0.3 m and we consider a relative amplitude a/h * 0 = 0.28. The topography is given by: b

(x) = 0 if x ≤ 15 m s (x -15) otherwise (90)
Following the empirical laws established by [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF], the trigger threshold ψ 0 of the virtual enstrophy and the dimensionless Reynolds number are calibrated as

ψ 0 =    g h * 0 0.1 + 0.031h * 0 a if a h * 0 > 0.05 0 otherwise Re = 0.85 + 60/s if a/h * 0 > 0.05 (91)
respectively, leading to approximative values ψ 0 ≈ 6.89 s -2 and Re ≈ 3.87 in the present case.

Still following [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF], the parameter C r may rely on the value C r = 0.48. This parameter is the same in all the forthcoming simulations. We can observe on Figure 5 a comparison between experimental and numerical results at several reference times during the propagation. In accordance with the experimental observations, the wave propagates, shoals and steepens when reaching the shore, as it can be seen at times t * = t(g/h 0 ) 1/2 = 10 and t * = 15. Then nonlinear effects induced by the bed elevation trigger the breaking of the wave, occurring between t * = 15 and t * = 20. To better understand the key role played by the enstrophy variables, a particular focus of the process is available on Figure 6. We can see that a sudden production of enstrophy occurs at t * = 18.1, allowing a precise identification of the transition to the turbulent regime. Once the wave is identified as breaking, the enstrophy production is able to correctly balance the dispersive effects and we recover the expected characterics of the wave transformation, until the end of the run-up phenomenon. Note that since the enstrophy is naturally following the wave motion, there is no need to artificially turn-off irrelevant terms or introduce additional de-breaking criteria in regions which are not concerned, as is can be the case with other existing approaches. As reported by Synolakis 1987 (see also [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF], a second breaking process happens at the end of the computation, under the form of a hydraulic jump around t * = 53. As shown in Figure 7, the process in also well captured by the proposed method, highlighting the enstrophy as a relevant diagnostic quantity with respect to the general detection of shock waves. These results were obtained with the third-order scheme; identical space and time locations have been observed for the breaking points using other space orders and/or different mesh resolutions. This low sensitivity with respect to the discretization parameters is not surprising, since the breaking criterion is somehow directly resolved as a model variable. Note finally that if a precise and predictive calibration of the breaking parameters (ψ 0 , Re) has been extracted from the works of [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF], no significant variability has been observed around the corresponding reference values, which also attests to the robustness of the proposed strategy.

Wave train breaking over a bar

We now turn to a classical test based on the laboratory experiments carried out by [START_REF] Beji | Experimental investigations of wave propagation over a bar[END_REF], devoted to the study of periodic sinusoidal waves propagating over a submerged bar. The experimental setup implies a 37.7 m long wave flume equipped with a trapezoidal bar with 1/20 front and 1/10 back slopes, separated by a 2 m plane area, leading to the configuration displayed in Figure 8. The evolution of the flow is followed at eight wave gauges disposed along the channel, for which experimental data are available. The first one is placed at x = 6 m while the others are located at the level of the bar, regularly spaced from x = 11 m to x = 17 m, as indicated in Figure 8. The water depth at rest was set to h * 0 = 0.4 m, leading to a 0.1 m water depth at the top of the bar.

Several series of experiments were run by [START_REF] Beji | Experimental investigations of wave propagation over a bar[END_REF], with varying amplitudes and frequencies. The objective pursued here is to show that the current approach can also be applied to capture a complete breaking process in the context of wave trains. This leads us to consider the experimental dataset obtained with a frequency of 0.4 Hz and a wave amplitude 0.054 m, corresponding to a strongly non-linear case. In this context, the shoaling of incoming waves is followed by a spilling-type breaking at the arrival at the flat part of the bar. After passing the bar, under the combined effect of the topography and their reintroduction in deeper waters, waves experience highly nonlinear deformations, accompanied by the development of high-order harmonics, close to the dispersive limits of the model. This induces a non-trivial coupling between turbulent effects and highly nonlinear deformations throughout the simulation. Capturing these complex dynamics is therefore a quite challenging issue, generally considered as an important step in the validation of numerical methods for breaking.

Computations have been run on the domain[-10 m, 40 m], including 10 m left and right sponge layers to generate incoming waves and to allow a proper exit of the outgoing waves. A regular mesh with 1200 elements in the x-direction was used for this test. If we refer to (91), the initial non-linearity has an order of magnitude allowing to set the activation criterion ψ 0 to zero, meaning that there is no breaking criterion and that the enstrophy is computed everywhere at the beginning of the computation. In such regimes of small relative amplitudes, the dimensionless Reynolds number has a value in the range of 5 < Re < 10 to fit with numerical experiments [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF]. Based on these observations, we directly resolve the effective (and not virtual) enstrophy equation, with Re = 7.5. As previously discussed, the dissipation parameter C r is taken to 0.48. As it can be observed on Figure 9, this choice allows to capture properly the breaking mechanism. It is noteworthy that the turbulent effects are not totally damped between two successive waves leading to an accumulation of enstrophy at the level of the bar. These numerical observations are physically relevant, since they highlight that the turbulence produced by a single wave can have not entirely disappeared at the arrival of the following wave as it has been reported in several works (see for instance [START_REF] Ting | Observation of undertow and turbulence in a laboratory surf zone[END_REF].

In any case, this accumulation phenomenon is limited in practice and does not break the stability of the method. Further, as shown in Figure 10, the method offers a good agreement with the experiments, significantly better than the ones obtained with usual methods.

Looking at Figure 10 more closely, the results at Gauge 2 show that the shoaling process is well described. In particular, we recover the correct amplitudes, while small overestimations can be observed with hybrid methods (see [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF]), presumably caused by numerical instabilities propagating in the neighbourhood of the breaking area (see also [START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF]. Wave breaking occurs at the level of Gauge 3. Our results reflect the ability in capturing accurately the process, with a precision similar to recent hybrid strategies [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF][START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF]. Gauges 4 and 5 allow to examine the continuity of the breaking process. We observe an almost perfect reproduction of the wave transformation, including the free surface inflection at the front side of the waves, which is difficult to capture by the models mentioned above. The passing at Gauge 6 marks the end of the breaking process, and is accompanied by a more pronounced manifestation of the wave decomposition into secondary waves. Up to our knowledge, only the turbulent kinetic energy model used by [START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF] is able to provide such a level of agreement at this gauge. As can be seen through Gauges 7 and 8, the model is also able to faithfully describe the end of the process, more successfully than the switching strategies [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF], or ad-hoc viscosity approaches [START_REF] Klonaris | A Boussinesq-type model including wave-breaking terms in both continuity and momentum equations[END_REF]. Note finally that the presented results appear to be in the same order of quality than those obtained with direct Computational Fluid Dynamics (CFD) simulations with ad hoc eddy viscosity [START_REF] Kamath | CFD Simulations of Wave Propagation and Shoaling over a Submerged Bar[END_REF].

Tsunami wave on a conical island

Based on the laboratory experiments of [START_REF] Liu | Runup of solitary waves on a circular island[END_REF], we now investigate the evolution of a solitary wave propagating over a conical island. This test is regularly employed to study

r being the distance in m from the center of the island (x 0 , y 0 ) = (12.96 m, 13.80 m). The initial water depth is h * 0 = 0.32 m. The flow evolution can be tracked through a serie of gauges covering the experimental domain, measuring the free surface elevation. Several data sets are available, implying different initial amplitudes of the incident wave. Here we chose a relative amplitude of a/h * 0 = 0.2, corresponding to the most important initial non-linearity. In this situation, as reported in [START_REF] Titov | Modeling of Breaking and Nonbreaking Long-Wave Evolution and Runup Using VTCS-2[END_REF], wave breaking is observed all around the island. The phenomenon is however not sufficiently pronounced to threaten the numerical stability and is generally neglected in this respect. The objective here is to show that the method is able to capture these small turbulent effects, while highlighting the possible related benefits.

For the numerical simulation we used the third-order scheme on a structured mesh of 49 566 elements, corresponding to a space step ∆x = ∆y = 0.167 m approximately. Still based on the 1D experimental laws, the virtual enstrophy threshold ψ 0 is set to 7.8 s -2 . As the dimensionless Reynolds number depends on the topography (91), a local strategy can be adopted to calibrate the viscosity parameter. However, we did not observe significant differences using a constant number all over the computational domain corresponding to an intermediate value based on the lower and higher slopes.

One can observe some 3D snapshots of the solution during the propagation in Figure 11, exhibiting a good reproduction of the flow characteristics. In particular, the passing of the emerged part of the island is well resolved, as well as the junction of the two resulting lateral waves at the rear side of the cone. Time series of the free surface are given in Figure 12 at gauges number 6, 9, 16 and 22, respectively located at (9.36 m, 13.80 m), (10.36 m, 13.80 m), (12.96 m, 11.22 m) and (15.56 m, 13.80 m). We observe a good agreement with the experimental data. In particular, as reported in [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF] and [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF] instance, slight overestimations of the leading wave amplitude at WG 9 and 22 can be observed when wave breaking is not accounted for (see also [START_REF] Fuhrman | Simulation of nonlinear wave run-up with a high-order Boussinesq model[END_REF]. Injecting enstrophy in the system seems to slightly reduce these discrepancies, especially for WG 22 and the minimum run-down amplitude of the reflected wave at WG 9. As confirmed in Figure 13, the method is able to identify wave breaking at the vicinity of the island during the passing of the wave. Note that similar results were obtained with a mesh of 150 000 elements to confirm that these amplitude corrections are not due to under-resolution.

Tsunami wave propagation over a 3D reef

This last test case is extracted from the laboratory experiments described in [START_REF] Swigler | Laboratory study investigating the three-dimensional turbulence and kinematic properties associated with a breaking solitary wave[END_REF].

The set-up implies the study of a solitary wave evolving within a realistic coastal configuration, including a three-dimensional fringing reef. This is a quite demanding test case, implying highly nonlinear transformations, wave breaking and treatment of shoreline motions in the presence of steep bottom variations. These mechanisms entail complex turbulent dynamics that make this benchmark in the line of the targeted applications of the proposed model. Recent twodimensional Boussinesq-type models used this set of data to validate their ability to describe the complexity of surf-zone mechanisms such as wave breaking and bore propagation driven by strongly varying bathymetries (see for instance [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF][START_REF] Roeber | Shock-capturing Boussinesq-type model for nearshore wave processes[END_REF][START_REF] Shi | A highorder adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF]. The experimental basin is 48.8 m long and 26.5 m wide. The bottom geometry reproduces a planar beach with a triangular flat reef surmounted by an idealized island located at the center, as can be seen in Figure 14 for instance. The offshore water depth is set to h * 0 = 0.78 m, and the beach slopes extend from x = 10 m to x = 32.5 m, so that a flat bottom is recovered shortly after the shoreline. The triangular shelf has a maximum elevation of 0.71 m and is linearly connected to the beach, with steeper slopes as we get closer to the shelf edge. The island is represented by a 6 m diameter cone of 0.45 m height centred at x = 12.6 m along the y-direction centreline. The flow motion is followed by means of a measurement device composed of nine wave gauges, where the free surface was recorded, supplemented by three Acoustic Doppler Velocimeters (ADV) to capture the velocity field. Their location is indicated in Table 1.

The relative amplitude of the incoming wave is a/h * 0 = 0.5, making this test highly nonlinear. The computational domain was extended at the inlet and outlet boundaries for the generation and the absorption of the solitary wave. The presented results were obtained with the secondorder scheme on a mesh of approximatively 200 000 elements, refined in the vicinity of the apex, leading to maximum and minimum areas max T ∈T h h T = 7.23 × 10 -2 m 2 and min T ∈T h h T = 6.57 × 10 -4 m 2 respectively. Following the empirical laws recalled in (91), the breaking threshold is fixed to ψ 0 = 1 s -2 and we set Re = 3.5, which corresponds to an intermediate value based on the maximum and minimum slopes. Naturally, other reasonable choice can be made for this parameter (with a definition depending on the local slopes for instance), but these variations did not had a significant influence on the numerical results.

The Figure 14 presents several snapshots of the numerical simulation of the propagation of a solitary wave towards the coast. The color maps the value of tr ϕ, which is the turbulent energy of the system times 2/h 2 , from deep blue (almost no turbulence) to red (highly turbulent). The red parts correspond to the foamy parts of breaking waves with strong vortices and rollers.

The first important steps of the propagation, corresponding to the passing of the central cone, can be observed on Figure 14 (a), (b), (c) and(d). As it can be seen through the color mapping, wave breaking is well captured by the model since the enstrophy is triggered at this stage. In accordance with the experiments and the works mentioned above, we observe a total submersion of the island. The shoaling and breaking processes seem to be well reproduced. In particular, just after the first breaking event at the arrival of the apex, one can see the enstrophy propagating from the center of the cone to the lateral boundaries, which clearly highlights a progressive transmission of the turbulent structures along the y-direction.

Figure 14 (e), (f) and (g) focuses on the termination of the breaking process and the wave propagation over the flat part of the beach. It can be seen that the enstrophy is progressively dissipated, while some residual turbulence can still be observed around the apex. The advancing wave front is well captured, without appearent numerical instabilities. Figure 14 (h), (i) and (j) proposes snapshots corresponding to the run-down process occuring towards the end of the simulation. Again the scheme is able to detect the formation of a hydraulic jump near the shoreline, which can be seen as the two-dimensional counterpart of what has been observed in 1D with Test 2. Comparisons with experimental data at wave gauges are displayed in Figure 15 for the free surface. They exhibit the capability of the model to predict accurately the arrival times of incident and reflected waves and to provide the correct amplitudes throughout the computational domain. Similar observations can be made with the time series of the velocity obtained by ADV proposed in Figure 16.

Conclusion

We derived a new two-dimensional depth-averaged model for coastal waves. This is an extension of the one-dimensional model of [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF]. The subdepth large-scale turbulence is resolved and taken into account by a tensor quantity called enstrophy. This tensor gives an anisotropic character to the model in accordance with the anisotropy of the three-dimensional subdepth large-scale turbulence. This tensor is also a source of vorticity of the mean flow and is responsible for transfers in both directions between the large two-dimensional horizontal eddies with a vertical vorticity and the three-dimensional subdepth turbulence. The small-scale turbulence is modelled with a turbulent viscosity hypothesis.

The three equations of the model consist of a scalar equation representing the mass conservation, a vector equation expressing the momentum balance and a tensor equation for the enstrophy tensor. The equation for the trace of this tensor is equivalent to the depth-averaged kinetic energy equation. The eddy viscosity and the turbulent dissipation are obtained by analogy with classical empirical laws. The presence of viscosity and dispersion implies the absence of discontinuities in the solution.

The model is fully nonlinear and its dispersive properties are equivalent to those of the Green-Naghdi equations. These properties were further improved by the method of [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] which can be directly extended to our model. The numerical resolution was obtained by an asymptotically equivalent formulation of the model equations according to the constant diagonal method of [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF] and by the Discontinuous Galerkin scheme of [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF]. It is noteworthy that the methods developed for the Green-Naghdi equations can be readily adapted to this model since the dispersive terms are identical. The tensor equation includes no dispersive terms nor terms due to the variable bottom. The additional equations are only transport equations with source terms without high-order derivatives, not even of second order. Consequently the numerical resolution complexity is not significantly increased by comparison with the Green-Naghdi equations. Compared to the hybrid or switching methods, there is no problem of mesh grid sensitivity nor of nonphysical oscillations.

A breaking criterion is not always needed. When it is needed, the breaking criterion of [START_REF] Kazakova | A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach[END_REF] can be taken directly as well as the values for the parameters of the model. There is no need for a breaking termination criterion. There is a low sensitivity to the precise values of these parameters and also to the discretization parameters, which gives an appreciable robustness to the model. These improvements are attributable to the presence in the model of a quantity representing explicitly the turbulence, namely the enstrophy tensor.

The model was used to simulate various cases of wave propagation, such as the run-up and run-down of a wave over a beach, the propagation of a wave train over a bar and the twodimensional waves propagation around an island or around a reef and on a sloping beach. The agreement with the experimental measures is very good in all cases. In particular the breaking phenomenon, including the hydraulic jump appearing in a run-down phase, are better described than in previous depth-averaged models. The numerical simulations show that this new model is numerically robust and that it has a predictive character and a high physical content without significant increase in numerical complexity. Further developments will include particle and sediment transport for applications in coastal erosion.

Acknowledgments. The work was supported by the Service Hydrographique et Océanographique de la Marine (SHOM).

A Conventions used in tensor calculus

For any vector (first-order tensor) V = V i e i , second-order tensors A = A ij e i ⊗ e j and A = A ij e i ⊗e j and third-order tensor B = B ijk e i ⊗e j ⊗e k (using Einstein notation), the dot product is defined as

A • A = A ik A kj e i ⊗ e j , (93)
the double dot product is defined as

A : A = A ij A ji , (94)
the divergence operator is defined as

div U = ∂U i ∂x i , div A = ∂A ij ∂x j e i , div B = ∂B ijk ∂x k e i ⊗ e j , (95)
and the gradient of a vector is defined as

grad U = ∂U i ∂x j e i ⊗ e j . ( 96
)
B Equations of the model

C Numerical scheme

We consider a computational domain Ω ⊂ R 2 , discretized with a conform triangular mesh T h . Denoting by T a generic element of the triangulation and by ∂T its boundary, we introduce the broken polynomial space

P k (T h ) := v ∈ L 2 | v |T ∈ P k (T ) , ∀ T ∈ T h , (108)
where P k (T ) denotes the space of two-variable polynomials in T of degree at most k. The area of T is denoted |T | and the notations h T , p T will respectively stand for its diameter and perimeter.

In the following, the notation N k is employed for the number of freedom degrees, equal to dim(P k (T )) = (k + 1)(k + 2)/2. The approximation space is X h = P k (T ) × P k (T ) 2 × P k (T ) 3 .

The computational time interval is denoted [0, t max ] and discretized in a sequence of intermediate times (t n ) n=0,N with a local time step ∆t n = t n+1 -t n . With notations similar to the continuous frame, we hence seek for an approximate solution W h = (Z h , hU T h , hϕ 11h , hϕ 12h , hϕ 22h ) of ( 79). The resulting semi-discrete formulation can be expressed through the local statement: find W h ∈ X h such that:

T d dt W h π h dx - T F(W h , b h ) • grad π h dx + ∂T F(W h , b h ) • n ∂T π h ds + T D(W h , b h )π h dx = T S(W h , b h )π h dx , (109)
for all π h ∈ P k (T h ) and all T ∈ T h . In the formulation above, b h stands for a polynomial expansion of the topography b on P k (T h ) and n ∂T is the unit outward normal to the boundary ∂T . Given a local polynomial expansion basis {φ i } N k i , the local restriction of the solution on a given element T can be written as:

W h|T (x, t) = N k i=1 W i (t)φ i (x) , x ∈ T , t ∈ [0, t max ] , (110)
where {W i } N k i are the local expansion coefficients, so that the formulation (109) leads to In the above expression, FT,F are interface terms approximating the projection of the fluxes F(W h , b h ).n T,F along the unit outward normal corresponding to the face F of T , defined in the spirit of finite-volume methods. In [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF], the preservation of motionless steady states Z = cte, U = 0 is guaranteed for solutions of an arbitrary order based on an adaptation of the hydrostatic reconstruction [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]). As a matter of fact, it can be shown easily that the pre-balanced formalism allows a trivial treatment of such configurations, as soon as the line and surface integrals are computed exactly and the topography admits a continuous discrete representation (which is automatically ensured with nodal expansions). As a consequence, with the same choice of nodal expansion basis {φ i } N k i for the bathymetry, there is no need to introduce modified states or additional correction term in the numerical fluxes. In light of this, classical Rusanov fluxes are used to evaluate the interface terms:

F T,F = F h (W -, W + , b -, b + , n T,F ) = 1 2 F(W -, b -) • n T,F + F(W + , b + ) • n T,F -a W + -W -, (112)
where the superscripts " -/ + " refer to the interior and exterior states respectively and a = max T ∈T h λ T , λ T referring to the maximum wave speed at the level of the element T :

λ T = max ∂T U • n ∂T ± gh + 3h 2 n ∂T • ϕ • n ∂T . (113)
Based on the works of [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF], the associated CFL condition is

λ T p T |T | ∆t n ≤ 1 2k + 1 , ∀T ∈ T h , (114)
but the presence of viscous terms in the model may constrain the time step to a parabolic stability condition. This point has to be taken into account during our numerical simulations.

An implicit treatment of these terms would allow to get rid of this restriction, but at the price of partially loose the computational efficiency and the ease of implementation of the method. Such investigations are left for future works.

Figure 1 :

 1 Figure 1: Notations used in the text.

  the deviatoric part of the residual stress tensor A r = 2ρν T D = A r ij e i ⊗ e j

  tr ϕ tr εϕ : ε ) .

Figure 2 :

 2 Figure 2: Example of regular mesh.

Figure 3 :Figure 4 :

 34 Figure 3: Test 1 : 2D view of the initial condition.

Figure 5 :

 5 Figure5: Test 2 : Comparison of the 1D free surface profiles between numerical results (solid curves) and experimental data from[START_REF] Synolakis | The runup of solitary waves[END_REF] (crosses) at several times during the propagation (t * = t(g/h 0 ) 1/2 ).

Figure 6 :

 6 Figure 6: Test 2 : Virtual enstrophy profile at several times around the breaking point. Wave breaking is identified at t * = 18.1, x = 30.3m.

Figure 7 :

 7 Figure 7: Test 2 : Virtual enstrophy profile at several times around the second breaking point (hydraulic jump). Wave breaking is identified at t * = 53, x = 33.5m.

Figure 8 :

 8 Figure 8: Test 3 : Sketch of the experimental configuration of Beji & Battjes (1993) and location of the wave gauges.

Figure 9 :

 9 Figure 9: Test 3 : Free surface (top) and enstrophy (bottom) profiles at three consecutive times during the breaking process.

Figure 10 :

 10 Figure 10: Test 3 : Time series of the free surface elevation at wave gauges. Comparison between the numerical results (curve) and the experimental data of Beji & Battjes (1993) (circles).

Figure 11 :

 11 Figure 11: Test 4 : 3D view of the free surface at t = 3 s; 6 s; 9 s; 10 s; 11.5 s and 12.5 s (from left to right and top to bottom).

Figure 12 :Figure 13 :

 1213 Figure 12: Test 4 : Time evolution of the free surface at reference gauges.

Figure 14 :Figure 15 :

 1415 Figure 14: Test 5 : Snapshots of the free surface during the run-up process. Color mapping with respect to the quantity tr ϕ = ϕ 11 + ϕ 22 in s -2 .

Figure 16 :

 16 Figure 16: Test 5 : Time evolution of the velocity components at ADVs.

  h , b h ) • grad φ j dx + F ⊂∂T F F T,F φ j ds + T D(W h , b h )φ j dx = T S(W h , b h )φ j dx , 1 ≤ j ≤ N k . (111)

Table 1 :

 1 Test 5 : Wave gauges and ADV's locations.

	for

  The complete system of equations can be written

																		∂h ∂t	+	∂hU ∂x	+	∂hV ∂y	= 0	(97)
	∂hU ∂t	+	∂ ∂x	hU 2 + h 3 ϕ 11 + Π + Π +	∂ ∂y	hU V + h 3 ϕ 12
							= -gh	∂b ∂x	+	∂ ∂x		4ν T h	∂U ∂x	+ 2ν T h	∂V ∂y	+	∂ ∂y	ν T h	∂U ∂y	+ ν T h	∂V ∂x	-f x (98)
	∂hV ∂t	+	∂ ∂x	hU V + h 3 ϕ 12 +	∂ ∂y	hV 2 + h 3 ϕ 22 + Π + Π
							= -gh		∂b ∂y	+	∂ ∂x		ν T h	∂V ∂x	+ ν T h	∂U ∂y	+	∂ ∂y	2ν T h	∂U ∂x	+ 4ν T h	∂V ∂y	-f y (99)
	∂hϕ 11 ∂t	+	∂hU ϕ 11 ∂x	+	∂hV ϕ 11 ∂y		= 2hϕ 11	∂V ∂y	-2hϕ 12	∂U ∂y
						+	ν T h	8			∂U ∂x			2	+ 4	∂U ∂x	∂V ∂y	+ 2	∂U ∂y	∂V ∂x	+ 2	∂U ∂y	2	-C r hϕ 11	√	ϕ 11 + ϕ 22 (100)
	∂hϕ 12 ∂t	+	∂hU ϕ 12 ∂x	+	∂hV ϕ 12 ∂y		= hϕ 12	∂U ∂x	+	∂V ∂y	-hϕ 11	∂V ∂x	-hϕ 22	∂U ∂y	+
						ν T h	3	∂U ∂x	∂U ∂y	+	∂V ∂x	∂V ∂y	+ 5	∂U ∂x	∂V ∂x	+	∂U ∂y	∂V ∂y	-C r hϕ 12	√	ϕ 11 + ϕ 22 (101)
	∂hϕ 22 ∂t	+	∂hU ϕ 22 ∂x	+	∂hV ϕ 22 ∂y		= 2hϕ 22	∂U ∂x	-2hϕ 12	∂V ∂x
						+	ν T h	8			∂V ∂y			2	+ 4	∂U ∂x	∂V ∂y	+ 2	∂U ∂y	∂V ∂x	+ 2	∂V ∂x	2	-C r hϕ 22	√	ϕ 11 + ϕ 22 (102)
	where																	Π =	gh 2 2	+	ḧ 3 h 2	(103)
																	Π =	h 2 2	D Dt	U	∂b ∂x	+ V	∂b ∂y	(104)
													f x =	h 2 ḧ	∂b ∂x	+ h	∂b ∂x	D Dt	U	∂b ∂x	+ V	∂b ∂y	(105)
													f y =	h 2 ḧ	∂b ∂y	+ h	∂b ∂y	D Dt	U	∂b ∂x	+ V	∂b ∂y	(106)
	and																	ν T =	h 2 R	√	ϕ 11 + ϕ 22 .	(107)