M Kazakova

G L Richard

A new model of shoaling and breaking waves: One-dimensional solitary wave on a mild sloping beach

We present a new approach to model coastal waves in the shoaling and surf zones. The model can be described as a depth-averaged large-eddy simulation model with a cutoff in the inertial subrange. The large-scale turbulence is explicitly resolved through an extra variable called enstrophy while the small-scale turbulence is modelled with a turbulent viscosity hypothesis. The equations are derived by averaging the mass, momentum and kinetic energy equations assuming a shallow water flow. The model is fully nonlinear and has the same dispersive properties as the Green-Naghdi equations. It is validated by numerical tests and by comparison with experimental results of the literature on the propagation of a onedimensional solitary wave over a mild sloping beach. The wave breaking is characterized by a sudden increase of the enstrophy which allows us to propose a breaking criterion based on the new concept of virtual enstrophy. The model features three empirical parameters. The first one governs the turbulent dissipation and was found to be a constant. The eddy viscosity is determined by a turbulent Reynolds number depending only on the bottom slope. The third parameter defines the breaking criterion and depends only on the wave initial nonlinearity. These dependences give a predictive character to the model which is suitable for further developments.

Introduction

Mathematical modelling of coastal wave propagation is a quite challenging issue since it is difficult to describe in the same model the dispersive effects in the shoaling zone and the energy dissipation of breakers in the surf zone. As it is impractical to solve the full Navier-Stokes system over any significant domain, approximate models must be used. Many asymptotic models were proposed in the last decades for the coastal wave propagation.

The first attempt to describe a surf zone wave propagation was made in the context of the nonlinear shallow-water equations, also called the Saint-Venant equations (Barré de Saint-Venant 1871). [START_REF] Meyter | Run-up on beaches[END_REF] reviewed analytic solutions of the shallow-water equations over a beach. [START_REF] Stocker | Water waves-The mathematical theory with applications[END_REF] presented perhaps the first numerical solution of the shallow-water equations for a sloping beach using the method of characteristics. His method produces quite accurate results for the simple cases. [START_REF] Freeman | Wave breakers on a beach and surges on a dry bed[END_REF] and later [START_REF] Iwasaki | On the shoreline and leading front conditions of tsunami waves in the light of the method of characteristics[END_REF] improved this method. However when the characteristic lines cross at the breaking point the wave evolution has to be treated in a different manner. Further for non-uniformly sloping beaches, this method is cumbersome.

In order to be able to calculate the wave breaking, [START_REF] Hibbert | Surf and run-up on a beach: a uniform bore[END_REF] proposed a groundbreaking method based on the numerical solution of the nonlinear shallow-water equations in their conservative form. This method gives physically realistic results. However it is not very robust in the run-up phase [START_REF] Synolakis | Discussion of "Wave Reflection and Run-Up on Rough Slopes[END_REF]). Furthermore, it is important to take into account a large energy dissipation following the wave breaking. Consequently subsequent version of this algorithm [START_REF] Packwood | Surf and run-up on beaches: models of viscous effects[END_REF][START_REF] Kobayashi | Wave reflection and run-up on rough slopes[END_REF] added artificial dissipative terms such as viscosity and friction for an accurate wave amplitude description.

Indeed the non-linear shallow-water equations give a rather good description of breaking waves representing them as shocks, and the run-up simulations generally give accurate results. However due to the absence of dispersive effects the wave profile is not described precisely particularly in the shoaling zone. Moreover this model predicts that all waves break which is obviously incorrect. Thereby, a more relevant modelling of wave shoaling before breaking is based on various Boussinesq-type models [START_REF] Brocchini | Nonlinear shallow water equations modeling for coastal engineering[END_REF]). Non-breaking numerical solutions for the Boussinesq equations are reviewed by [START_REF] Liu | Report on the international workshop on long-wave run-up[END_REF]. Since the 1990's considerable efforts were made to improve the dispersive properties of the Boussinesq model [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part II: a slowly varying bathymetry[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF][START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part I: highly nonlinear unsteady waves[END_REF].

Most of Boussinesq models are derived in the assumption of a weak nonlinearity (the wave amplitude is much smaller than the characteristic water depth). The natural idea is to use a fully nonlinear model, as the one derived by [START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF] and [START_REF] Su | Korteweg-de Vries equation and generalizations, III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF] in the onedimensional (1D) case, or the two-dimensional (2D) fully nonlinear model presented in [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF]. This model can be derived from the Euler equations as an asymptotic model in the shallow water regime without any hypothesis on the wave amplitude (see e.g. [START_REF] Lannes | The Water Waves Problem: mathematical analysis and asymptotics[END_REF]). This fully-nonlinear model has better dispersive properties than the Boussinesq equations, and the water wave profile is better described. Except for being formulated in terms of the velocity vector at an arbitrary level, the equations of [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part I: highly nonlinear unsteady waves[END_REF] mentioned above are basically equivalent to the 2D Green-Naghdi equations. In fact the Green-Naghdi equations are often known in the literature as the fully-nonlinear Boussinesq-type model. The dispersive properties of the Green-Naghdi equations were further improved by [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] and Chazel et al. (2011).

However the weakly nonlinear or fully nonlinear Boussinesq type equations do not include dissipation due to wave breaking, and thus become invalid in the surf zone. To extend the validity of those equations, [START_REF] Heitner | Numerical model for tsunami run-up[END_REF] introduced an artificial viscosity term into the momentum balance equation. This approach was followed by many researchers. [START_REF] Zelt | The run-up of nonbreaking and breaking solitary waves[END_REF] used an eddy viscosity formulation together with a Lagrangian Boussinesq model. The inclusion of this term in the momentum equation helps to control the energy dissipation, and it must be calibrated with experiments. Similar techniques were used in [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], [START_REF] Wei | Generation of waves in Boussinesq models using a source function method[END_REF] and [START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking and runup. I: 1[END_REF]. Another idea is to use the so-called roller models which include dissipation through an extra convective term in the momentum equation (see for example [START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF][START_REF] Svendsen | A Boussinesq breaking wave model with vorticity[END_REF][START_REF] Karambas | Breaking waves in the surf and swash zone[END_REF][START_REF] Briganti | Boussinesq modeling of breaking waves: description of turbulence[END_REF][START_REF] Musumeci | The flow in the surf zone: a fully nonlinear Boussinesq-type of approach[END_REF][START_REF] Dimas | Surface roller model for the numerical simulation of spilling wave breaking over constant slope beach[END_REF][START_REF] Cienfuegos | Wave-breaking model for Boussinesqtype equations including roller effects in the mass conservation equation[END_REF][START_REF] Viviano | A nonlinear rotational, quasi-2DH, numerical model for spilling wave propagation[END_REF]. The basis of these models is the concept of surface roller due to [START_REF] Svendsen | Wave heights and set-up in a surf zone[END_REF]. In this approach additional terms are included not only in the momentum equation, but also in the mass balance equation and the thickness of the roller region has to be estimated.

In the so-called switching or hybrid approach the dispersive terms are suppressed in the breaking regions. Contrary to the nonlinear shallow-water system, the Boussinesq equations do not admit discontinuities because they are dispersive. The removal of the dispersive terms reduces the system to the nonlinear shallow water equations. This relatively simple idea was developed in many studies (for example [START_REF] Tonelli | Simulation of wave breaking over complex bathymetries by a Boussinesq model[END_REF][START_REF] Shi | A highorder adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF][START_REF] Duran | Discontinuous Galerkin discretization of a new class of Green-Naghdi equations[END_REF], 2017[START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF][START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF].

Generally speaking, two large groups of methods can be distinguished. The first one includes additional terms such as viscous terms and the second one uses hybrid methods. In both approaches some parameters need to be calibrated. Moreover some criteria must be defined to initiate or terminate the breaking process. Generally those criteria are based on different empirical relations. For example the breaking process is activated if the wave local energy dissipation forms a peak, or if the wave slope reaches a critical value.

Detailed analysis of the different approaches cited above and an estimation of the influence of the numerical dissipation can be found in [START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF]. The authors highlight that despite the absence of an explanation for the physical meaning of the coefficients in the eddy viscous terms, this approach shows little sensitivity to the grid, while hybrid methods are not very robust and lead to numerical instabilities that depend on the mesh size.

The aim of the present study is to derive a dispersive model for breaking wave propagation which resolves the large-scale turbulence. In most models the turbulence is not resolved but modelled with a turbulent viscosity hypothesis. In some models, this turbulent viscosity is calculated with a turbulent kinetic energy for which a transport equation with source terms is solved [START_REF] Nwogu | Numerical prediction of breaking waves and currents with a Boussinesq model[END_REF][START_REF] Karambas | Breaking waves in the surf and swash zone[END_REF], Zhang et al. 2014). However the validity of a turbulent viscosity hypothesis is questionable (see the discussion in [START_REF] Pope | Turbulent Flows[END_REF]. It is often assumed a local balance between production and dissipation, especially in the case of a mixinglength approach [START_REF] Nwogu | Numerical prediction of breaking waves and currents with a Boussinesq model[END_REF], which is not assured in a breaking wave. Further the anisotropic character of shallow-water turbulence cannot be well described by such an approach [START_REF] Nadaoka | Shallow-water turbulence modeling and horizontal largeeddy computation of river flow[END_REF]. [START_REF] Kim | A depth-integrated model for weakly dispersive, turbulent, and rotational fluid flows[END_REF] distinguished between a vertical and a horizontal eddy viscosity, the latter being estimated with the Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. The basic equations[END_REF] involving the horizontal grid size. [START_REF] Nadaoka | Shallow-water turbulence modeling and horizontal largeeddy computation of river flow[END_REF] and [START_REF] Hinterberger | Three-dimensional and depth-averaged large-eddy simulation of some shallow water flows[END_REF] used a filtering approach but resolved only the two-dimensional eddies with a length scale larger than the fluid depth while the subdepth scales were not resolved and were modelled with a turbulent viscosity hypothesis. Note that these models were not especially designed for coastal waves but for other shallow-water flows. [START_REF] Madsen | Subgrid modelling in depth integrated flows[END_REF] used a filtering approach with a grid size comparable to the water depth and also an eddy viscosity hypothesis calculated with the Smagorinsky model. This model was applied to nearshore tidal currents. As pointed out by [START_REF] Nadaoka | Shallow-water turbulence modeling and horizontal largeeddy computation of river flow[END_REF], whereas the subgrid turbulence in the original approach of [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. The basic equations[END_REF] is reasonably isotropic and locally in equilibrium since the cutoff frequency is in the inertial subrange, it is not the case for the subdepth turbulence. It follows that the validity of a turbulent viscosity hypothesis and of the Smagorinsky model is again questionable for a cutoff at or near the water depth length scale. Further an eddy viscosity approach implies that the energy can be transferred only from the large scales toward the small scales i.e. there is no backscatter. However an energy transfer from the subdepth turbulence toward the two-dimensional eddies of larger lengthscale is possible [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]. To take such a transfer into account, [START_REF] Hinterberger | Three-dimensional and depth-averaged large-eddy simulation of some shallow water flows[END_REF] introduced a stochastic force in their model.

In this work we use a filtering method with a cutoff frequency in the inertial subrange as in the classical large-eddy simulation (LES) approach, before averaging the resulting filtered equations over the depth. It is thus much safer to assume that the residual (or subgrid) turbulence is isotropic and locally in equilibrium. This method resolves the large-scale subdepth turbulence and has the advantage to capture the anisotropic character of this turbulence, as well as the possible energy transfer towards the larger scales. Of course, these effects will appear in the two-dimensional (2D) model for three-dimensional flows which will be the subject of a work to come. However the other advantage of this method is that it can be assumed that the energy transfer rate from the filtered scales towards the residual scales is almost equal to the dissipation i.e. that there is a close balance between the production of the residual kinetic energy and the dissipation [START_REF] Lilly | The representation of small-scale turbulence in numerical simulation experiments[END_REF]). This result is important for the derivation of the 1D model.

In the derivation of Boussinesq-type models or of the Green-Naghdi model, there is an assumption on the velocity profile. In many cases it is assumed either that the horizontal velocity component is uniform or almost uniform over the whole water depth or that the flow is irrotational, which leads to a parabolic velocity profile. Other assumptions were formulated. For example, the roller model of [START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF] supposes that the velocity profile is flat below the roller region and that the velocity in the roller region is also constant and equal to the wave celerity which has thus to be estimated. More complex velocity profiles including rotational effects were also proposed in the literature. [START_REF] Zhang | Rotational surf zone modeling for O(µ 4 ) Boussinesq-Green-Naghdi systems[END_REF] used an asymptotic method with a polynomial expansion of the horizontal velocity to include the rotational effects. In some models, the vorticity equation in the fluid depth was also explicitly included to calculate the rotational effects on the velocity profile. The variation of the vorticity with the depth was calculated analytically [START_REF] Veeramony | The flow in surf-zone waves[END_REF][START_REF] Musumeci | The flow in the surf zone: a fully nonlinear Boussinesq-type of approach[END_REF] or numerically [START_REF] Briganti | Boussinesq modeling of breaking waves: description of turbulence[END_REF]). Most often these velocity profiles are either incorrect for breaking waves (notably the flat profile or the hypothesis of an irrotational flow) or complex and leading to rather complicated models. Further a velocity profile implies that the turbulence is not resolved and included only with a turbulent viscosity hypothesis. If the large-scale turbulence is resolved, a velocity profile is not practicable.

An alternative approach is to get rid entirely of any velocity profile by adding a new variable to the equations. The addition of an extra variable implies to add also an equation to the model. For a 1D model, this new equation is given by the averaged kinetic energy equation. [START_REF] Svendsen | A turbulent bore on a beach[END_REF] use the energy as the extra variable and solved the energy equation. However because they divided vertically the flow into a turbulent region and a non-turbulent region both having a variable thickness, a velocity profile had to be assumed. In this case, this profile was a third-order polynomial. Moreover the pressure was taken hydrostatic which means that there was no dispersive effect and that the model was restricted to turbulent bores. A 2D hyperbolic model of shear flows without any velocity profile assumption was derived in [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF] from the Euler equations also with the assumption of a hydrostatic pressure. When augmented with dissipation [START_REF] Richard | A new model of roll waves : comparison with Brock's experiments[END_REF][START_REF] Lannes | The Water Waves Problem: mathematical analysis and asymptotics[END_REF] it gives an excellent description of roll waves and classical hydraulic jumps. The closure of the model was obtained with the hypothesis of weakly-sheared flows which is not very restrictive in practice since highly turbulent hydraulic jumps can be correctly described. The addition of a non-hydrostatic correction to the pressure and thus of dispersive effects was proposed in the conservative case in [START_REF] Castro | Fully nonlinear long-waves models in presence of vorticity[END_REF] and [START_REF] Richard | Modelling turbulence generation in solitary waves on shear shallow water flows[END_REF]. [START_REF] Gavrilyuk | Spilling breakers in shallow water: applications to Favre waves and to the shoaling and breaking of solitary waves[END_REF] first included both dissipative and dispersive effects in the framework of this approach with a two-layer model, the upper turbulent layer including shearing effects and the irrotational lower layer being described by the Green-Naghdi model. In the present work this method is extended to include dissipation and dispersion in a one-layer model.

The paper is organised as follows, in §2 the filtered equations are presented. Then the depthaverage model is derived in §3 and the empirical laws for the eddy viscosity and the turbulent dissipation are discussed in §4. The last two sections deal with numerical tests. The numerical implementation is presented §5. In §6 the model is validated by comparison with experimental results.

Filtered conservation equations

We study the propagation of coastal waves from the shoaling zone to the shoreline. The flow can be highly turbulent, especially in the breaking zone and thereafter in the surf zone. A suitable model of turbulence is thus needed to capture important physical features of the breaking waves.

The classical approach to turbulence is based on the Reynolds decomposition. The velocity field is written as the sum of the mean velocity and of the turbulent fluctuation. The mean kinetic energy of the flow can be decomposed into the kinetic energy of the mean flow and the turbulent kinetic energy. Turbulent processes usually remove energy from the mean flow and transfer it to the fluctuating velocity field. This transfer is called production since most of the time it refers to a loss of the mean kinetic energy and a production of turbulent kinetic energy.

The turbulent motions range in size from large scales, which are of the order of the characteristic lengthscales of the flow, to small scales. According to the energy cascade hypothesis [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF] and to the Kolmogorov hypothesis [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluids for very large Reynolds numbers[END_REF], the production transfers energy first to the large-scale motions. This energy is then transferred to smaller and smaller scales until the Kolmogorov scales where the energy is dissipated by viscous processes. The large scales contain most of the energy and are therefore called the energy-containing range. The dissipating small scales are called the dissipation range. Between these two ranges, at sufficiently high Reynolds numbers lies the inertial subrange where the lengthscales are high enough so that the viscous effects are negligible but small enough compared to the lengthscales of the flow such that there is almost no production. In the inertial subrange, the energy is transferred by inviscid processes toward the smaller scales. Another feature of turbulent motions according to Kolmogorov hypothesis is that the large eddies are anisotropic whereas the small eddies are isotropic.

We use an approach similar to the large-eddy simulation (LES) method. The velocity field v is filtered to decompose the velocity into a filtered velocity field v and a residual velocity field

v r v = v + v r . ( 1
)
The difference between the Reynolds decomposition and this filtering operation is that the former decomposes the velocity field into a non-turbulent (mean) field and a turbulent field whereas in the latter the filtered velocity field includes the large-scale turbulence and the residual velocity field includes the small-scale turbulence. Ideally, the filtering operation is a low-pass filter that allows to resolve turbulent motions of a scale greater than some specified length chosen in the inertial subrange. It follows that the anisotropic energy-containing range is in the filtered field and that the isotropic dissipation range is in the residual field. Another difference with the Reynolds decomposition is that, in general, the filtered residual velocity is not equal to zero i.e. v r = 0, although if the filter is a projector such as the sharp spectral filter, v = v and v r = 0 (for more details on all this approach, see [START_REF] Pope | Turbulent Flows[END_REF].

The filtering operation is applied to the Navier-Stokes equations of an incompressible fluid of density ρ and kinematic viscosity ν. The filter is supposed to be homogeneous. The filtered continuity equation becomes simply div v = 0.

(2)

The filtered momentum equation is

∂v ∂t + div (v ⊗ v) = g - 1 ρ grad p + ν∆v (3)
where p is the filtered pressure and ⊗ the tensorial product. The residual stress tensor is defined by

σ r = -ρ (v ⊗ v -v ⊗ v) . (4)
The residual kinetic energy is defined from the trace of this tensor as

k r = - tr σ r 2ρ . (5)
The residual stress tensor can be decomposed into an isotropic part and a deviatoric anisotropic part A r as

σ r = - 2 3 ρk r I + A r . ( 6
)
The tensor I is the identity tensor. The residual stress tensor is then modelled by a turbulentviscosity hypothesis. This usual hypothesis implies that there is no backscatter in the model i.e.

the energy is transferred only from the large scales toward the small scales. Denoting by D the filtered strain rate tensor defined by

D = 1 2 grad v + (grad v) T , (7)
the anisotropic residual stress tensor is written

A r = 2ρν T D (8)
where ν T is a turbulent viscosity. The isotropic part of the residual stress tensor is then absorbed into a modified pressure which is denoted simply by p to lighten the notations and which is

p = p + 2 3 ρk r . ( 9
)
The filtered momentum equation is finally written

∂v ∂t + div (v ⊗ v) = g - 1 ρ grad p + div 2 ν T D + ν∆v. ( 10
)
The filtered specific kinetic energy is

e k = 1 2 v • v. ( 11
)
It follows from the definition ( 5) that

e k = e f + k r (12)
where e f is the kinetic energy of the filtered velocity field

e f = 1 2 v • v. ( 13
)
Writing g = -grad e p , the latter satisfies the equation

∂e f ∂t + div e f v + pv ρ - v • A r ρ -2νv • D + e p v = -ε f -P r (14)
where

ε f = 2νD : D (15)
and

P r = 2ν T D : D. ( 16
)
The dot and the colon mean the dot product and the double dot product respectively. For a very high Reynolds number, the terms involving the molecular viscosity can be neglected. This is the case of the term ε f which represents the viscous dissipation from the filtered velocity field. Most of the viscous dissipation processes take place in the dissipation range. The dominant dissipative term in the equation of e f is P r which represents an energy transfer from the filtered motions towards the residual motions i.e. from the large scales towards the small scales. Since most of the energy is contained in the energy-containing range, the mean kinetic energy e is almost equal to the mean kinetic energy of the filtered velocity field e f (the brackets denote the Reynolds averaging). The dissipation due to the mean flow being negligible at high Reynolds numbers, the dissipation of the mean kinetic energy is almost entirely due to the dissipation of the turbulent kinetic energy, usually denoted by ε and called simply dissipation. Consequently, the mean dissipation of e f , denoted by P r is nearly equal to the dissipation. This also implies 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], 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])

P r ε. ( 17
)
The classical closure of this approach is due to [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. The basic equations[END_REF] who modelled the eddy viscosity by analogy with the mixing-length hypothesis. We do not follow this path because the filtered equations will be averaged over the fluid depth. We will use instead a simpler closure by assuming that the turbulent viscosity is uniform over the water depth but otherwise variable in space and time. The empirical law chosen for this eddy viscosity is discussed in §4.

Averaged conservation equations

Governing equations

We study the particular case of a two-dimensional flow over a variable bottom. The notations are presented in Figure 1. The components of the velocity field are u in the horizontal direction Ox and w in the vertical direction Oz. The bottom topography is measured by its elevation b(x) over a horizontal datum. The elevation of the free surface over this horizontal datum is Z(x, t) while the water depth is h(x, t) = Z(x, t) -b(x). The still water depth is denoted by h 0 (x) and the water elevation over this level is η(x, t). In the case of a solitary wave, the amplitude a of the wave is defined as the maximum value of the elevation i.e. a = η max .

The filtered mass conservation equation is

∂u ∂x + ∂w ∂z = 0. ( 18
)
The filtered momentum balance equation writes in the Ox direction

∂u ∂t + ∂u 2 ∂x + ∂uw ∂z = - 1 ρ ∂p ∂x + 1 ρ ∂A r xx ∂x + ∂A r xz ∂z + ν ∂ 2 u ∂x 2 + ∂ 2 u ∂z 2 (19)
and in the Oz direction

∂w ∂t + ∂uw ∂x + ∂w 2 ∂z = -g - 1 ρ ∂p ∂z + 1 ρ ∂A r xz ∂x + ∂A r zz ∂z + ν ∂ 2 w ∂x 2 + ∂ 2 w ∂z 2 . ( 20
)
The balance equation for the kinetic energy of the filtered motions can be written The boundary conditions are the no-penetration condition at the bottom

∂ ∂t u 2 2 + w 2 2 + ∂ ∂x u u 2 2 + w 2 2 + gz + pu ρ - A r xx u ρ - A r xz w ρ - τ xx u ρ - τ xz w ρ + ∂ ∂z w u 2 2 + w 2 2 + gz + pw ρ - A r xz u ρ - A r zz w ρ - τ xz u ρ - τ zz w ρ = -ε f -P r . (
w(b) = u(b)b , (22)
where b = db/dx, the kinematic boundary condition at the free surface

w(Z) = ∂h ∂t + u (Z) ∂Z ∂x (23)
and the dynamic boundary condition at the free surface

(σ • n) (Z) = 0 (24)
where σ = -pI + A r + τ is the Cauchy stress tensor and where n is the unit normal vector at the free surface.

Scaling

The model is derived within the scope of the shallow-water approximation. If h * 0 is a reference value of the still-water depth and if L is a characteristic length of variation of the flow parameters in the horizontal direction, then there is a small parameter

µ = h * 0 L 1. ( 25
)
The equations are then written in dimensionless form. The dimensionless quantities are denoted by a tilde symbol. The following scaling is classical for this situation [START_REF] Antuono | Beyond Boussinesq-type equations: semi-integrated models for coastal dynamics[END_REF]):

x = x L ; z = z h * 0 ; ũ = u gh * 0 ; w = w µ gh * 0 ; t = µt g h * 0 ; h = h h * 0 ; p = p ρgh * 0 . ( 26
)
For the viscous stress tensor, the scaling is

τxx = Lτ xx ρν gh * 0 ; τzz = Lτ zz ρν gh * 0 ; τxz = τ xz ρν h * 0 g . ( 27
)
The Reynolds number is defined by Re = h * 0 gh * 0 /ν. The turbulent viscosity is supposed to be of O(µ) and it is written in dimensionless form as in [START_REF] Antuono | Beyond Boussinesq-type equations: semi-integrated models for coastal dynamics[END_REF]:

νT = ν T µh * 0 gh * 0 . ( 28
)
Consequently the components of the deviatoric part of the residual stress tensor are scaled as

Ãr xx = A r xx µ 2 ρgh * 0 ; Ãr zz = A r zz µ 2 ρgh * 0 ; Ãr xz = A r xz µρgh * 0 ; (29)
We can write

Ãr xx = 2ν T ∂ ũ ∂ x ; Ãr zz = -Ãr xx ; Ãr xz = νT ∂ ũ ∂ z + µ 2 ∂ w ∂ x ; (30)
The dimensionless mass conservation equation is simply

∂ ũ ∂ x + ∂ w ∂ z = 0. ( 31
)
The momentum equation in the Ox direction becomes in dimensionless form

∂ ũ ∂ t + ∂ ũ2 ∂ x + ∂ ũ w ∂ z = - ∂ p ∂ x + µ 2 ∂ Ãr xx ∂ x + ∂ Ãr xz ∂ z + µ Re ∂ τxx ∂ x + 1 µRe ∂ τxz ∂ z . (32)
In the Oz direction, the momentum equation can be written

µ 2 ∂ w ∂ t + ∂ ũ w ∂ x + ∂ w2 ∂ z = -1 - ∂ p ∂ z + µ 2 ∂ Ãr xz ∂ x + µ 2 ∂ Ãr zz ∂ z + µ Re ∂ τxz ∂ x + µ Re ∂ τzz ∂ z . ( 33
)
It follows from the chosen scaling that the dimensionless viscous dissipation in the filtered velocity field can be defined as εf = h * 0 ε f /(gν) and that the dimensionless energy transfer toward the residual motion is P r = P r /(µg gh * 0 ). The dimensionless energy equation can then be written

∂ ∂ t ũ2 2 + µ 2 w2 2 + ∂ ∂ x ũ ũ2 2 + µ 2 w2 2 + z + pũ -µ 2 Ãr xx ũ + Ãr xz w - µ Re (τ xx ũ + τxz w) + ∂ ∂ z w ũ2 2 + µ 2 w2 2 + z + p w -Ãr xz ũ -µ 2 Ãr zz w - 1 µRe τxz ũ - µ Re τzz w = - 1 µRe εf -P r . ( 34
)
The boundary conditions are also written in dimensionless form with b = b/h * 0 and Z = Z/h * 0 . If we assume that b varies horizontally on a characteristic length of O(L), then b = b /µ. The no-penetration condition on the bottom and the kinematic condition at the free surface write simply

w(b) = ũ(b) b ; w(Z) = ∂ h ∂ t + ũ(Z) ∂ Z ∂ x , (35)
whereas the dynamic boundary condition at the free surface (24) gives two scalar equations

Ãr xz (Z) + p(Z) -µ 2 Ãr xx (Z) ∂ Z ∂ x = 0 (36) p(Z) -µ 2 Ãr zz (Z) + µ 2 Ãr xz (Z) ∂ Z ∂ x = 0. ( 37
)
There is no surface tension and no shear stress condition imposed at the free surface.

The filtered equations are averaged over the depth. For any quantity A, its average value over the depth is denoted by A and defined as

A = 1 h Z(x,t) b(x) A dz. ( 38
)
The averaged horizontal velocity is denoted by U = u and in dimensionless form Ũ = ũ .

The filtered horizontal velocity is decomposed as

u(x, z, t) = U (x, t) + u (x, z, t). (39)
Equivalently, u = U + u + u r . The term u represents the deviation of the filtered horizontal velocity from its averaged value. It includes the anisotropic large-scale turbulence which contains most of the turbulent energy and the depth variations of the Reynolds-averaged velocity field (shearing effects). On the other hand, the residual velocity field u r includes the isotropic smallscale turbulence where most of the dissipation takes place. The model is derived in the hypothesis of a weakly turbulent (and weakly sheared) flow which means that u is of O(µ). This hypothesis is very similar in effect as the hypothesis of a weakly-sheared flow [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF]) with a value of Teshukov's exponent β equal to 1. Therefore we write

ũ(x, z, t) = Ũ (x, t) + µũ (x, z, t). (40)
In the following derivation of the model, all terms of an order up to O(µ 2 ) will be kept and the terms of O(µ 3 ) will be neglected.

The Reynolds number is supposed to be high enough so that all viscous terms can be neglected. This hypothesis is usual for this kind of problems (see for example [START_REF] Antuono | Beyond Boussinesq-type equations: semi-integrated models for coastal dynamics[END_REF]. It can be written

Re = O(µ -3 ). The terms of O(µ/Re) are then of O(µ 4 ). The terms of O(1/(µRe)) involve τxz whose dominant term is ∂ ũ/∂ z. Since ∂ Ũ /∂ z = 0, τxz is in fact of O(µ)
and the terms of O(1/(µRe)) are of O(µ 3 ) and therefore negligible. The term εf is negligible for the same reason.

The effect of the hypothesis of a weakly turbulent flow can also be applied to the scaling of D and P r . In D the dominant component is normally, in dimensionless form

Dxz = 1 2 ∂ ũ ∂ z + µ 2 ∂ w ∂ x . (41)
However, ∂ ũ/∂ z = µ∂ ũ /∂ z. This implies that P r = P r µ 3 g gh * 0 .

(42)

Consequently, the last term in (34) writes in fact -µ 2 P r . It is thus of O(µ 2 ) and must be kept.

As in [START_REF] Veeramony | The flow in surf-zone waves[END_REF], the shear stress on the bottom is neglected and this implies a free-slip condition at the sea bed. The whole deviatoric residual stress tensor is neglected on the bottom i.e. A r (b) 0.

Averaging procedure

To lighten the notations, the tilde are dropped in this section dealing entirely with dimensionless quantities. Averaging the mass conservation equation (31) over the depth, taking into account the boundary conditions, gives ∂h ∂t

+ ∂hU ∂x = 0. ( 43
)
The continuity equation ( 31) allows us to find an expression for the vertical filtered velocity

w = (z -b) ḣ h + U b + O (µ) (44)
where the dot symbol denotes the material derivative along the average flow

ḣ = ∂h ∂t + U ∂h ∂x . ( 45
)
The momentum equation in the Oz direction (33) gives first

∂p ∂z = -1 -µ 2 (z -b) ḧ h -µ 2 Db U Dt + µ 2 ∂A r xz ∂x + µ 2 ∂A r zz ∂z . ( 46
)
and then, with the dynamic boundary condition (37), an expression of the pressure

p = Z -z -µ 2 ḧ h z 2 -Z 2 2 -b (z -Z) -µ 2 Db U Dt (z -Z) + µ 2 ∂ ∂x z Z A r xz dz + µ 2 A r zz . ( 47
)
In these expressions, both the notations Da/Dt and ȧ denote the material derivative of a quantity a along the average flow (see ( 45)). With the dynamic boundary condition ( 36) and since we took A r (b) 0, the integration of the pressure terms and of the residual stress terms gives

- Z b ∂p ∂x dz + Z b µ 2 ∂A r xx ∂x dz + Z b ∂A r xz ∂z dz = - ∂ ∂x Z b p dz -p(b)b + µ 2 ∂ ∂x Z b A r xx dz. ( 48
)
The integration of the pressure leads to

Z b p dz = h 2 2 + µ 2 h 2 ḧ 3 + µ 2 h 2 2 Db U Dt + µ 2 Z b dz ∂ ∂x z Z A r xz dz + µ 2 Z b A r zz dz. ( 49
)
We can write

z Z A r xz dz = z Z ν T ∂u ∂z dz + O(µ 2 ). ( 50
)
Since the turbulent viscosity is uniform over the fluid depth,

z Z ν T ∂u ∂z dz = ν T (u -u(Z)) . ( 51
)
The hypothesis of a weakly turbulent flow implies that u -u(Z) = O(µ). The corresponding term in ( 49) is thus of O(µ 3 ) and negligible. Then

Z b A r xx dz = - Z b A r zz dz. ( 52
)
With the expression of A r xx , we get

Z b A r xx dz = 2ν T Z b ∂u ∂x dz. ( 53
)
The hypothesis of a weakly turbulent flow allows us to write u

(Z) = U + O(µ) and u(b) = U + O(µ). This leads to Z b A r xx dz = 2ν T h ∂U ∂x . ( 54
)
The calculation of the pressure at the sea bed gives

p(b)b = hb + µ 2 hb ḧ 2 + Db U Dt . ( 55
)
Finally the integration of the right-hand side of equation ( 32) gives

- Z b ∂p ∂x dz + Z b µ 2 ∂A r xx ∂x dz + Z b ∂A r xz ∂z dz = - ∂ ∂x h 2 2 + µ 2 h 2 ḧ 3 + µ 2 h 2 2 Db U Dt -4ν T h ∂U ∂x -hb -µ 2 hb ḧ 2 + Db U Dt . (56)
Taking into account the boundary conditions, the integration of the left-hand side of equation ( 32) leads to

Z b ∂u ∂t + ∂u 2 ∂x + ∂uw ∂z dz = ∂hU ∂t + ∂ ∂x h u 2 . ( 57
)
The treatment of u 2 is the same as in [START_REF] Richard | A new model of roll waves : comparison with Brock's experiments[END_REF]. First, u 2 is written U 2 +µ 2 u 2 since, by definition, u = 0. Second, a new variable, called enstrophy, is defined as

ϕ = u 2 h 2 . ( 58
)
The averaged momentum equation can finally be written

∂hU ∂t + ∂ ∂x hU 2 + µ 2 h 3 ϕ + h 2 2 + µ 2 Π + µ 2 Π = -hb -µ 2 f + O(µ 3 ) ( 59
)
where

Π = h 2 ḧ 3 -4ν T h ∂U ∂x , (60)
Π = h 2 2 Db U Dt ( 61
)
and

f = hb ḧ 2 + Db U Dt . ( 62
)
We will say that the sea bed has a mild slope if b varies horizontally on a characteristic length of O(µL). In this case, the terms in µ 2 Π and -µ 2 f become of O(µ 3 ) and are negligible. The averaged momentum equation reduces to

∂hU ∂t + ∂ ∂x hU 2 + µ 2 h 3 ϕ + h 2 2 + µ 2 h 2 ḧ 3 -4ν T h ∂U ∂x = -hb + O(µ 3 ). ( 63
)
Since the model features three variables h, U and ϕ, three equations are needed. Two equations are provided by the averaged mass conservation equation and by the averaged momentum equation and the third equation is the averaged energy equation. With the boundary conditions, we can integrate

Z b ∂ ∂t u 2 2 + ∂ ∂x u 3 2 + ∂ ∂z wu 2 2 dz = ∂ ∂t h 2 u 2 + ∂ ∂x h 2 u 3 . ( 64
)
We can expand as above u 2 = U 2 +µ 2 h 2 ϕ. The same method yields u 3 = U 3 +3µ 2 h 2 U ϕ+ O(µ 3 ). As in [START_REF] Teshukov | Gas-dynamics analogy for vortex free-boundary flows[END_REF], the hypothesis of weakly turbulent flows allows us to neglect u 3 and thus to close the problem. Another integration combined with the boundary conditions and with the expression (44) of the vertical velocity leads to

µ 2 Z b ∂ ∂t w 2 2 + ∂ ∂x uw 2 2 + ∂ ∂z w 3 2 dz = µ 2 ∂ ∂t h ḣ2 6 + hU 2 b 2 2 + hb U ḣ 2 + µ 2 ∂ ∂x U h ḣ2 6 + hU 2 b 2 2 + hb U ḣ 2 . ( 65
)
With the expression (47) of the pressure and with the boundary conditions, we get

Z b ∂ ∂x uz + pu -µ 2 A r xx u -µ 2 A r xz w + ∂ ∂z wz + pw -A r xz u -µ 2 A r zz w dz = ∂ ∂x U h 2 2 + µ 2 h 2 ḧ 3 + µ 2 h 2 2 Db U Dt -µ 2 4ν T h ∂U ∂x + hb U. ( 66
)
The averaged value of P r over the depth is estimated by the averaged value over the depth of its mean value (in the Reynolds sense) P r . Then the approximation (17) gives

P r ε . (67)
Physically, this means that the energy dissipation of the model is mostly due to the dissipation of the turbulent kinetic energy in the dissipation range i.e. in the small scales. The expression of ε is discussed in §4. The averaged energy equation can finally be written

∂ ∂t h e + µ 2 e + ∂ ∂x hU e + µ 2 e + U µ 2 h 3 ϕ + h 2 2 + µ 2 Π + µ 2 Π = -hb U -µ 2 h ε + O(µ 3 ) (68)
where Π and Π are given by ( 60) and ( 61) respectively and where

e = U 2 2 + µ 2 h 2 ϕ 2 + µ 2 ḣ2 6 + h 2 ; e = b 2 U 2 2 + ḣ 2 b U. ( 69
)
The factor µ 2 in front of h ε comes from the hypothesis of a weakly turbulent flow as explained at the end of §3.2 (see relation ( 42)). In the case of a mild slope, e and Π become negligible and the energy equation reduces to

∂he ∂t + ∂ ∂x hU e + U µ 2 h 3 ϕ + h 2 2 + µ 2 Π = -hb U -µ 2 h ε + O(µ 3 ). ( 70
)

Enstrophy equation

The model is composed of the mass equation ( 43), the momentum equation (59 or 63 for a mild slope) and the energy equation (68 or 70 for a mild slope). From these three equations, an evolution equation for the enstrophy can be derived. From equation (59), we first obtain the equation

∂ ∂t hU 2 2 + ∂ ∂x hU 3 2 + U ∂ ∂x µ 2 h 3 ϕ + h 2 2 + µ 2 Π + µ 2 Π = -hb U -µ 2 U f . (71)
Subtracting this equation from equation ( 68) leads to

∂ ∂t µ 2 h 3 ϕ 2 + µ 2 h ḣ2 6 + h 2 2 + µ 2 he + ∂ ∂x µ 2 h 3 U ϕ 2 + h 2 U 2 + µ 2 hU ḣ2 6 + µ 2 hU e + U µ 2 h 3 ϕ + h 2 2 + µ 2 Π + µ 2 Π -U ∂ ∂x µ 2 h 3 ϕ + h 2 2 + µ 2 Π + µ 2 Π = -µ 2 h ε + µ 2 U f . (72)
The mass equation implies that

ḣ = -h ∂U ∂x . ( 73
)
The following equalities are satisfied

∂he ∂t + ∂ ∂x hU e + U Π = U ∂Π ∂x + U f , (74)
∂ ∂t h ḣ2 6 + ∂ ∂x hU ḣ2 6 + U h 2 ḧ 3 = U ∂ ∂x h 2 ḧ 3 , (75)
∂ ∂t h 2 2 + ∂h 2 U ∂x = U ∂ ∂x h 2 2 (76) and ∂ ∂t h 3 ϕ 2 + ∂ ∂x 3h 3 U ϕ 2 -U ∂ ∂x h 3 ϕ = h 3 2 Dϕ Dt . (77)
The evolution equation of the enstrophy can finally be written

∂hϕ ∂t + ∂hU ϕ ∂x = 8ν T h ∂U ∂x 2 - 2 h ε . ( 78
)
The enstrophy is related to the averaged large-scale turbulent kinetic energy and to the shearing effects of the Reynolds-averaged flow. The above equation shows that it is created by the effect of the turbulent viscosity and dissipated by the dissipation of the turbulent kinetic energy. The latter takes place in the small scales after the energy had been transferred from the large scales.

The underlying conservative hyperbolic system has the mathematical structure of the Euler equations of compressible fluids. For this system, the enstrophy is analogous to the entropy. If this hyperbolic system had to be resolved, the appearance of discontinuities would impose to solve the mass, momentum and energy equations since these discontinuities would conserve the energy and would create enstrophy. Because of the diffusive and dispersive terms of our model, no discontinuity can arise. It is thus equivalent in theory to solve the mass, momentum and energy equations or to solve the mass, momentum and enstrophy equation. However, in practice, the enstrophy equation ( 78) is much simpler than the energy equation ( 68). The enstrophy equation has no dispersive term and no term depending on the topography. Numerically, this equation is much easier to handle than the energy equation. The substitution of the enstrophy equation for the energy equation in order to solve numerically the equations was already made by [START_REF] Lannes | Nonlinear wave-current interactions in shallow water[END_REF] in their dispersive model including terms up to O(µ 3 ) but no dissipation nor viscosity. Take care of the notation differences since their µ is for us µ 2 and the quantity they denote by Ẽ is hϕ.

Dissipation and eddy viscosity

The closure of the model composed of equations ( 43), (59) (or (63)) and ( 78) requires to know the dissipation ε and the eddy viscosity ν T . We can assume that these quantities can depend on h and on the enstrophy ϕ but not on the average velocity U in order to satisfy easily the principle of Galilean invariance (ϕ is a Galilean invariant quantity). Including a dependence on U would in fact requires a dependence on a velocity difference U -v * where v * would be a suitable velocity. The only obvious choice for v * would be the bottom velocity but it is difficult to explain a dependence on the bottom velocity while the turbulent and viscous processes close to the bottom are completely neglected in this model. Another choice for the turbulent viscosity would be a dependence on h and g writing ν T = C ν h √ gh as in [START_REF] Musumeci | The flow in the surf zone: a fully nonlinear Boussinesq-type of approach[END_REF], where C ν is a dimensionless constant. However such an expression would imply that the turbulent viscosity is greater where the depth is greater, thus in the shoaling zone, and smaller where the wave breaks and in the surf zone whereas the opposite variation would be expected. This would make the model highly dependent on the breaking criterion. It seems much more preferable to be inspired by the models using a dependence of the viscosity with the turbulent kinetic energy.

With a dependence on h and ϕ only, a dimensional analysis shows that

ν T = C p h 2 √ ϕ, (79)
where C p is a dimensionless quantity, and that

ε = C r 2 h 2 ϕ 3/2 (80)
where C r is another dimensionless quantity. Note that C p can be interpreted as the inverse of a Reynolds number R such that

ν T = h 2 √ ϕ R . ( 81
)
With these choices, the model has the mathematical structure of the turbulent-kinetic-energy model (TKE or km ) proposed by [START_REF] Kolmogorov | The equations of turbulent motion in an incompressible fluid[END_REF] and [START_REF] Prandtl | Über ein neues Formelsystem für die ausgebildete Turbulenz[END_REF] who suggested to base the turbulent viscosity and the dissipation on the turbulent kinetic energy and the mixing length. In our model, the enstrophy and the water depth replace the turbulent kinetic energy and the mixing length respectively (the turbulent kinetic energy is in fact homogeneous to h 2 ϕ) with the difference that the problem of completeness of the TKE model does not apply to this model since h is variable of the model and does not need to be specified. Another difference is that, although the equation of the enstrophy is analogous to the equation of the turbulent kinetic energy of the TKE model, the mass equation ( 43) shows that our model is analogous to the equations of a compressible fluid (the water depth h being the analogous of the density).

The relation (80) implies that the dimensionless form of ε is

ε = ε µ 3 g gh * 0 . ( 82
)
This is in accordance with (42). The relation ( 81) is also in agreement with the scaling (28) of the eddy viscosity. It follows that the quantities C p , R and C r are all of O(1). In this approach the scaling and the dependence of the turbulent viscosity have a clear physical meaning and interpretation. There is no need to impose a variation of the turbulent viscosity with the depth which would greatly complicate the model for a doubtful benefit. The viscous terms are needed for their diffusive role and to create enstrophy. The only important thing is to be able to predict the value of the eddy viscosity in all cases and the hypothesis of a viscosity which is constant over the depth is sufficient for this purpose as it is shown in the following sections.

Numerical resolution

The resulting system for a mild slope topography can be written in dimensional form

                     ∂h ∂t + ∂hU ∂x = 0, ∂hU ∂t + ∂ ∂x hU 2 + gh 2 2 + h 3 ϕ + h 2 ḧ 3 = ∂ ∂x 4 R h 3 √ ϕ ∂U ∂x -ghb , ∂hϕ ∂t + ∂hU ϕ ∂x = 8h √ ϕ R ∂U ∂x 2 -C r hϕ 3/2 (83)
This system is dispersive, and the discretization of the high-order derivatives requires a special approach. It should be mentioned that, from the point of view of the dispersive terms, the proposed system (83) can be treated in the same manner as the well-known Green-Naghdi system. The same techniques for the numerical resolution of the Green-Naghdi equations can be used for the evaluation of the high-order derivatives or for the preservation of the steady states.

In the last two decades the Green-Naghdi system received a lot of attention due to its improved dispersive properties and the following list of works is not intended to be exhaustive. Different approaches were applied for the one-dimensional system. A finite-difference method is proposed in Antunes Do [START_REF] Do Carmo | Numerical solution of the generalized Serre equations with the MacCormack finite difference scheme[END_REF] and in Antunes Do Carmo (2013). A finite-volume method is developed in [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: model development and analysis[END_REF][START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: boundary conditions and validation[END_REF], Le [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF] and [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF]. [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF] used a pseudo-spectral approach. Hybrid methods are implemented in [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF], [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] and [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF]. For the continuous finite element method and for the discontinuous Galerkin method we refer to [START_REF] Mitsotakis | On the Galerkin/finite-element method for the Serre equations[END_REF][START_REF] Mitsotakis | A modified Galerkin/finite element method for the numerical solution of the Serre-Green-Naghdi system[END_REF] for the former and to [START_REF] Panda | Discontinuous Galerkin methods for solving Boussinesq-Green-Naghdi equations in resolving non-linear and dispersive surface water waves[END_REF], [START_REF] Dong | A reconstructed central discontinuous Galerkin-finite element method for the fully nonlinear weakly dispersive Green-Naghdi model[END_REF] and [START_REF] Duran | Discontinuous Galerkin discretization of a new class of Green-Naghdi equations[END_REF] for the latter. Two-dimensional numerical studies including finitedifference, finite-volume and finite-element methods, as well as hybrid methods can be found in Le [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF], [START_REF] Shi | A highorder adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF], Antunes Do Carmo (2013), [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF], [START_REF] Popinet | A quadtree-adaptative multigrid solver for the Serre-Green-Naghdi equations[END_REF] and [START_REF] Duran | A discontinuous Galerkin method for a new class of Green-Naghdi equations on simplicial unstructured meshes[END_REF]. In this study we concentrate only on the one-dimensional equations to validate the proposed model. Moreover in the following numerical tests, a slowly-varying topography was assumed.

The numerical integration of the Green-Naghdi equations generally includes an elliptic operator inversion. The change of variables proposed in Le [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF] allows to split the numerical resolution into two steps. Firstly, the system is rewritten as a hyperbolic system of equations and treated with a Godunov's type method. Then an elliptic equation for the velocity is solved. We adopted this strategy for the model (83) due to its straightforward realization.

Numerical scheme

The characteristics of the hyperbolic part of this system are

λ 1 = U, λ 2,3 = U ± gh + 3h 2 ϕ. ( 84
)
Following Le [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF] we introduce the new variables

K = U + 1 3h ∂(h 2 ḣ) ∂x , ( 85
) α = - 2 3 h 3 ∂(hU ) ∂x . ( 86
)
Under this change of variables the system (83) can be rewritten

                 ∂h ∂t + ∂hU ∂x = 0 ∂hK ∂t + ∂ ∂x hU K + gh 2 2 + h 3 ϕ + α = ∂ ∂x 4 R h 3 √ ϕ ∂U ∂x -ghb ∂hϕ ∂t + ∂hU ϕ ∂x = 8h √ ϕ R ∂U ∂x 2 -C r hϕ 3/2 (87)
This system is solved for the conservative variables h, hK, hϕ by a finite-volume scheme (Godunovtype) with the HLL approximate Riemann solver for the hyperbolic part of the equations and by a finite-difference method for the source terms. The viscous and topography source terms being not stiff an explicit scheme was used. The stability of this scheme in the presence of viscosity is ensured by a Courant-Friedrichs-Levy (CFL) condition.

The topography source term in the momentum equation should be discretized in a special manner to preserve the steady state h + b = const, U = K = 0 and ϕ = 0. To do so, we followed the method of [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] to obtain a well-balanced reconstruction.

The second step of the method of Le [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF] is the resolution of the elliptic equation ( 85) to calculate U from the values of h and K given by the first step. We rewrite this equation in terms of hU and hK hK = hU -1 3

∂ ∂x h 2 ∂(hU ) ∂x + 1 6 ∂ ∂x hU ∂h 2 ∂x .
This equation is solved numerically at each time step by using a simple second-order finitedifference method. The inversion of a tridiagonal matrix is needed, which is done with the Thomas algorithm. This step is the most time-consuming of all the numerical procedure.

A scheme of the second order in space and time was obtained with the Monotone Upstream-Centred Scheme for Conservation Laws (MUSCL) method coupled with the Heun method.

Convergence of the numerical scheme

The numerical scheme is validated in the case of the classical soliton solution. The Green-Naghdi system admits an exact solution which corresponds to a wave propagating without deformation. This solution was found by [START_REF] Su | Korteweg-de Vries equation and generalizations, III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF] and coincide with the expression derived by Rayleigh (1876). This particular kind of solitary wave was first observed by [START_REF] Russell | Report on waves[END_REF]. The model ( 83) admits also such a type of solution on a flat horizontal bottom in the absence of dissipation. The exact expressions for the water depth h ex , the fluid velocity U ex and the enstrophy ϕ ex were derived in a more general case by [START_REF] Richard | Modelling turbulence generation in solitary waves on shear shallow water flows[END_REF]. In this particular case, these expressions reduce to ϕ ex = ϕ 0 = constant,

h ex h * 0 = 1 + 2ã F r 2 -1 -3 φ F r 2 -1 -(3 + ã2 ) φ + [F r 2 -1 -(3 -ã2 ) φ] cosh [κ (x -c 0 t -x 0 ) /h * 0 ] (88)
and

U ex = c 0 1 - h * 0 h ex , (89)
where x 0 is the initial abscissa of the wave maximum depth position, h * 0 is the still water depth at an infinite distance from the wave, c 0 is the celerity of the soliton, ã = δ = a/h * 0 its dimensionless amplitude (and the nonlinearity parameter), F r = |c 0 |/ gh * 0 a Froude number representing the dimensionless wave celerity, φ = ϕh * 0 /g the dimensionless enstrophy and where The celerity of the soliton can be written

κ = 3 F r 2 (F r 2 -1 -3 φ). ( 90
)
c 0 = (h * 0 + a) [g + ϕ (3h * 0 + a)]. (91)
The dimensionless amplitude of the wave ã is

ã = 1 2 φ -(1 + 4 φ) + (1 + 4 φ) 2 + 4 φ (F r 2 -1 -3 φ) . (92)
This soliton solution reduces to the expression of Rayleigh ( 1876) and [START_REF] Su | Korteweg-de Vries equation and generalizations, III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF], corresponding to the 1D Green-Naghdi equations, if ϕ 0 = 0. If ϕ 0 = 0, this solution describes a wave with a smaller amplitude for the same Froude number compared to the Rayleigh solution (see [START_REF] Richard | Modelling turbulence generation in solitary waves on shear shallow water flows[END_REF]. The depth profile of a soliton calculated for F r = 1.095, h * 0 = 1 m and ϕ 0 = 0.05 s -2 is presented in Figure 2 where it is compared to the depth profile of the Rayleigh soliton (ϕ 0 = 0) for the same reference water depth and Froude number. The main effect of the enstrophy is to decrease the soliton amplitude.

The numerical scheme is used to calculate the propagation of a soliton over a flat horizontal bottom with F r = 1.095, h * 0 = 1 m and ϕ 0 = 0.05 s -2 . The depth profiles are presented in Figure 3 at different times and compared to the exact solution. No deviation can be observed from the theoretical profile. It is thus numerically established that the proposed numerical scheme preserves the dynamics of the solitons.

This case is used to study the convergence of the algorithm. The deviation of the calculated solution h(t, x) from the exact soliton solution h ex (t, x) is estimated for an increasing level of space refinement. The discrete analogue of the norm L ∞ in time and L 2 in space of the relative error is plotted in Figure 4. This proves numerically that the scheme is of the second order.

Breaking criterion

Turbulence generation

The phenomenon of wave breaking is characterized by a rapid increase of turbulence and vorticity and by a change of the shape of the wave profile which becomes more asymetrical. During the breaking process the kinetic energy of the mean flow is first converted into turbulent kinetic energy of organized vortices before being dissipated into small-scale turbulence (see for example [START_REF] Lubin | Three-dimensional large-eddy simulation of air entrainment under plunging breaking waves[END_REF]). In the framework of our model this means that the breaking process corresponds to a rapid increase of the enstrophy which represents the turbulence of large scale. The spectral analyses resulting from the experimental measures of [START_REF] Hattori | Experimental study on turbulence structures under breaking waves[END_REF] confirm the energy cascade picture in the case of breaking waves and also the formation of vortices.

Before the beginning of the breaking process the flow can be considered to be irrotational except close to the bottom and close to the free surface. In the case of breaking waves the main source of turbulence is the injection of turbulence at the free surface [START_REF] Hattori | Experimental study on turbulence structures under breaking waves[END_REF]. The existence of a small vorticity near the free surface is due to the surface curvature (Longuet-Higgins 1992) and it is thus mainly localized at the wave crest. However [START_REF] Lin | A numerical study of breaking waves in the surf zone[END_REF] showed that, compared to the vorticity generated by the breaking wave, this small crest vorticity is not only much smaller but also it has the opposite sign. This implies that the characteristic vorticity and vortices of breaking waves are generated during the breaking of the wave. In the shoaling zone, the interior of the flow is irrotational and in this case the Corrsin-Kistler equation [START_REF] Corrsin | Free-stream boundaries turbulent flows[END_REF] shows that the Reynolds stresses have no effect on the mean velocity field. The vorticity is created on the free surface during the wave breaking and is transported into the interior region of the wave [START_REF] Lin | A numerical study of breaking waves in the surf zone[END_REF]. The surface-generated turbulence penetrates downward too during breaking and can reach the bottom wall [START_REF] Hattori | Experimental study on turbulence structures under breaking waves[END_REF]. That the origin of vorticity and turbulence generation is located at or near the free surface during breaking is further confirmed by the major role of air entrainment in turbulence generation [START_REF] Lubin | Three-dimensional large-eddy simulation of air entrainment under plunging breaking waves[END_REF].

The precise location of vorticity and turbulence generation was determined by several experimental and numerical works. In the case of a plunging breaker [START_REF] Kimmoun | A particle image velocimetry investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF] identified the generation of a first vortex at the beginning of the plunging phase when the free surface steepens strongly. They observed that this vortex is not located at the top of the crest but rather at the toe of the front of the wave. Other vortices are generated in the front part of the wave during a series of splash-up events [START_REF] Lubin | Three-dimensional large-eddy simulation of air entrainment under plunging breaking waves[END_REF][START_REF] Kimmoun | A particle image velocimetry investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. The production of turbulent kinetic energy is also located at the front face of the crest and it remains attached there during the whole wave propagation from the breaking point to the swash zone [START_REF] Kimmoun | A particle image velocimetry investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. In a spilling breaker the major source of turbulence is at the foot of the breaker when the water falling, tumbling, or sliding, down the front meets the undisturbed water [START_REF] Brocchini | The dynamics of strong turbulence at free surfaces. Part 2. Free-surface boundary conditions[END_REF]. The experiments and the numerical simulations show that the vorticity generation and the turbulence generation take place during breaking in approximately the same front region of the wave extending from the foot of the breaker to the front side of the breaking crest. This region is depicted in grey in Figure 5 in the case of a plunging breaker (inspired by the results of Lin &[START_REF] Lin | A numerical study of breaking waves in the surf zone[END_REF][START_REF] Kimmoun | A particle image velocimetry investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF][START_REF] Kimmoun | A particle image velocimetry investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. Figure 5(a) shows the steepening phase at the beginning of the plunging process and Figure 5(b) shows the wave after the second splash-up event. Vorticity and turbulent kinetic energy are also measured behind this region because both quantities are transported after they are created.

A depth-averaged model cannot describe all the details of the breaking process and especially not the splash-up sequences. The vorticity and turbulence generation during the steepening of the free surface is taken into account by an enstrophy creation caused by the production term in the enstrophy equation ( 78). When the free surface steepens, (∂U/∂x) 2 and consequently the production term increase and enstrophy is created. However the particular generation mechanism that takes place at the front side of the crest and at the toe of the front cannot be captured by a depth-averaged model. The information on the sign of ∂U/∂x is lost in the production term. Consequently enstrophy can be produced not only at the front side of the wave but also at the rear side, particularly if the wave profile is weakly asymetrical. In the case of a soliton the wave profile is completely symetrical. In the absence of a description of the generation of vorticity and turbulence on the front surface of the wave, there is a production on both sides of the crest. A sloping bottom creates some asymmetry in the wave but this does not always prevent an enstrophy production in the rear of the wave. After breaking, thanks to the dissipation, the wave becomes asymetrical enough to remove the problem but a criterion must be added to prevent a non-physical behaviour in the breaking zone.

Since the turbulent fluctuations have no effect on the mean velocity field if the flow is irrotational [START_REF] Corrsin | Free-stream boundaries turbulent flows[END_REF], there should be no enstrophy and no eddy viscosity in the shoaling zone before the beginning of the surface generation of vorticity at the wave front. This gives the grounds for the introduction of a breaking criterion. The first element of a breaking criterion is to impose that the enstrophy production can happen only at the front side of a wave. This is obtained simply by taking ν T = 0 if ∂U/∂x > 0. This condition is equivalent to ḣ < 0 since the mass equation implies that the material derivative of the water depth is ḣ = -h∂U/∂x. The front part of the wave is thus defined as ḣ > 0. There is many cases where this condition is needless and, even when it is useful, its omission causes only a transient problem that is rapidly corrected by the action of the dissipation. However, since its implementation is easy, it seems that there is no reason not to include it. It remains to be seen if this condition is sufficient or if a second element is needed for the breaking criterion.

Virtual enstrophy

The nonlinearity parameter of the initial wave is defined by δ = a * /h * 0 where a * is the initial amplitude of the wave and where h * 0 is the still water level at the initial position of the wave. When δ < 0.05 the model ( 83) can be used directly without any breaking criterion other than the condition of §5.3.1. The breaking position is predicted correctly due to the sudden increase of the enstrophy. An exemple of a solitary wave with δ = 0.048 propagating on a constant slope equal to tgβ = 1/60 is presented in Figure 6(a) where the maximal value of the enstrophy is given as a function of the position during the wave propagation. This is a simulation of an experiment of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF] whose geometrical characteristics are presented in Figure 7. In this case, the still water depth at the initial position was h * 0 = 1.2 m. The values of the parameters of the model are C r = 0.48 and R = 6 and these values can be used for all solitary waves with δ < 0.05 on this slope. The increase of the enstrophy corresponds to the breaking point since it implies also a brutal increase of the large-scale turbulent energy and a diminution of the wave amplitude. As shown in Figure 6(a), the enstrophy keeps a very small value from the beginning of the waves propagation until the breaking region where the sudden increase takes place. In the shoaling zone the values of the enstrophy range from 10 -7 s -2 to 10 -4 s -2 whereas in the surf zone, the enstrophy is of the order of 1 s -2 or 10 s -2 , thus several orders of magnitude higher. This evolution is also presented in Figure 6(b) with another scale showing that the enstrophy is much smaller than 10 -3 s -2 in the shoaling zone. The extreme steepness of the enstrophy curve at the breaking point shows that the enstrophy is a relevant quantity to characterize the wave breaking.

However, for solitary waves with a nonlinearity greater than about 0.05, the initial evolution of the enstrophy causes an important attenuation of the wave amplitude. Compared to the case with δ < 0.05, the increase of the enstrophy at the breaking point is less steep and the enstrophy is also less small in the shoaling zone. The enstrophy is still a relevant quantity for breaking but its values are not small enough before breaking to be entirely negligible and this is the cause of the wave attenuation. Therefore for such conditions the enstrophy must not be created in the shoaling zone. This implies that the turbulent viscosity should be activated only when the breaking is likely to occur and that a breaking criterion is needed for these waves.

There are several breaking criteria based on empirical relations for the horizontal velocity, the free-surface gradient or the local energy dissipation (see e.g. [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]. For other criteria we refer to [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF], [START_REF] Bacigaluppi | A 1D stabilized finite element model for nonhydrostatic wave breaking and run-up[END_REF] or [START_REF] Duran | Discontinuous Galerkin discretization of a new class of Green-Naghdi equations[END_REF]. Although most of those criteria can be adapted to our model, we introduce a new criterion specific to this model and based on the sudden increase of the enstrophy at the breaking point. The idea is to calculate the enstrophy produced by the wave while preventing any feedback on the wave's characteristics to avoid the spurious amplitude attenuation. A new quantity denoted by ψ is introduced. It follows the same equation as the enstrophy ϕ except that the turbulent viscosity is always activated in the equation of ψ if ∂U/∂x < 0. The equation of ψ is thus

∂hψ ∂t + ∂hU ψ ∂x = B 8h √ ψ R ∂U ∂x 2 -C r hψ 3/2 , (93)
where B = 0 if ∂U/∂x > 0, otherwise B = 1 (see §5.3.1). The equations (83) of the model are unchanged and they don't depend on ψ. This means that ψ has no effect on the wave propagation. This new quantity, homogeneous to the enstrophy, represents the amount of enstrophy that the wave is potentially able to create and therefore it is called virtual enstrophy.

The virtual enstrophy undergoes a sudden when the is about to break and this can be used to characterize the breaking point. The eddy viscosity in the equation of the real enstrophy is set equal to zero as long as the virtual enstrophy did not reach a threshold value ψ 0 . This means that there is no enstrophy creation before the breaking point and thus no amplitude attenuation in the shoaling zone. The cell where ψ becomes equal to ψ 0 indicates the beginning of the breaking process. The breaking point can be defined as this position or as the position where the wave amplitude begins to decrease which happens shortly after the former position. Once ψ > ψ 0 , the eddy viscosity in the real enstrophy equation takes its normal value (taking into account the condition of §5.3.1). Consequently the real enstrophy increases rapidly and the expected breaking process is obtained. Note that the initial values of the virtual enstrophy and of the real enstrophy should never be set equal to zero because in this case the turbulent viscosity would be also equal to zero. Instead of zero, a very small initial uniform value ranging from 10 -12 s -2 to 10 -7 s -2 must be given to ψ and ϕ in all cells.

A criterion to remove the eddy viscosity if the wave stops breaking is easy to define along the same guidelines but for the moment it proved to be useless in all test cases.

The precise threshold value of the virtual enstrophy is found by comparison with the experiments and is presented below. The advantage of the proposed criterion is that its implementation is straightforward. Since this criterion depends on the value of one of the variables of the model, Figure 8: Variation of the maximal value ψ max of the virtual enstrophy for a soliton calculated with R = 2, h * 0 = 2, 2 m and a nonlinearity parameter δ equal to 0.3 (thin solid curve), 0.5 (dotdashed curve), 0.7 (dashed curve), 0.78 (thick solid curve) and 0.9 (dotted curve).

it does not depend on the mesh size nor on the order of the scheme not even on the numerical scheme provided that the calculation is close to the convergence.

The evolution of the maximum value of the virtual enstrophy for a solitary wave with a nonlinearity parameter δ = 0.137 and the conditions of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF] with h * 0 = 2.2m is presented in Figure 6(a). The increase of ψ just before the breaking point is rapid but it is less steep than for δ = 0.048. With the different scale presented in Figure 6(b), it can be seen that the maximum value of ψ increases first from its initial value 10 -9 s -2 to a constant value of about 0.2 • 10 -3 s -2 corresponding to the part of the propagation over a horizontal bottom. When the wave arrives over the sloping bottom (see Figure 7), the maximum value of the virtual enstrophy is no longer constant and increases until the threshold value for breaking is reached (see Figure 6(b)). Although the maximum value of ψ in the shoaling zone is small, it is several orders of magnitude greater than for δ = 0.048.

Propagation of a soliton in a constant water depth

When used with the breaking criterion involving the virtual enstrophy, the model ( 83), including all dissipation terms, allows a soliton to propagate indefinitely over a flat horizontal bottom provided that the virtual enstrophy is smaller than the threshold value (see §5.3.2). If this condition is satisfied, the eddy viscosity remains equal to zero and the system reduces to its conservative part and thus admits a soliton solution. It is known that a soliton cannot physically exist if the nonlinearity parameter of the wave is greater than some limit value. The most commonly accepted limit value is δ = 0.78 (McCowan 1894). For a nonlinearity greater than this limit, the wave breaks and dissipation appears.

The propagation of a soliton is simulated numerically taking the expression of the soliton of the Green-Naghdi equations (Rayleigh's soliton) as the initial wave conditions and an enstrophy ϕ equal to 0 (alternatively an infinitesimal value can be used without detriment). The virtual enstrophy ψ is initially set to an infinitesimal value, for example of the order of 10 -9 s -2 . The calculated evolution of the maximum virtual enstrophy is presented in Figure 8 for R = 2, h * 0 = 2.2 m and a nonlinearity ranging from 0.3 to 0.9. The maximum virtual enstrophy increases initially quickly and then becomes forever constant. The value of this constant depends on the nonlinearity parameter of the wave and on the value of R. The maximal value of the virtual enstrophy is greater if the nonlinearity is greater (see Figure 8) or if R is smaller. Note that this constant value does not depend on the initial value of the virtual enstrophy provided that [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF].

it is small enough. There is much leeway to choose this initial value though, since according to our tests any value smaller than about 10 -6 s -2 is acceptable. Although it is not shown on the figure, the entire virtual enstrophy profile becomes also invariant during the propagation after the transient regime and not only the maximal value. This assumes that the maximum virtual enstrophy in the permanent regime is smaller than the threshold value. Otherwise the eddy viscosity is activated and the wave breaks. Consequently the threshold value must be chosen equal to the maximum virtual enstrophy value of the most nonlinear existing soliton. Taking the limit of [START_REF] Mccowan | On the highest wave of permanent type[END_REF], the threshold value corresponds to δ = 0.78. However, as the maximum virtual enstrophy value depends also on R, a more complete study of wave breaking on various slopes and of various nonlinearities must be made before being able to propose empirical laws for R and for the threshold values ψ 0 . In particular, since R controls not only the increase of the virtual enstrophy but also the dissipation in breaking waves, the determination of the law governing R requires to study also the breaking waves.

6 Application to a mild-slope topography 6.1 Comparison with the experiments of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF] We study now the propagation of a solitary wave over a mild-slope topography. The experiments of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF] are used to validate the model (83). In these experiments a solitary wave is generated with a high-resolution wavemaker and it propagates in a channel 300 m long, 5.0 m wide and 5.2 m deep. A plane beach with a slope equal to tgβ = 1/60 starts 50 m after the wavemaker. The setup is shown on Figure 7. The initial still water depth is denoted by h * 0 , the initial wave amplitude by a * and the unperturbed water depth at the breaking point by h b . Several tests were made for an initial depth h * 0 equal to 1.2 m, 2.2 m or 2.9 m and for different values of the nonlinearity parameter δ ranging from 0.019 to 0.338. The breaking position is given for each test case and the evolution in time of the water depth at different positions along the channel are available. The experimental conditions of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF] which we used for our numerical simulations are given in Table 1 with the values of the relevant parameters. These 14 cases were chosen to cover a wide range of nonlinearities and depths. Each case, except for one, is referred by the same trial number as in [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF].

The localization of the breaking points is a quite challenging issue. Moreover the definition of the breaking point itself can have several interpretations. For example, the amplitude decreases dramatically for a strongly nonlinear wave when breaking occurs. In [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF] the breaking point is defined as the location where the front of the wave becomes nearly vertical tangent and where bubbles subsequently appear. We can relate this to the beginning of a rapid turbulent energy increase i.e. to a sudden enstrophy growth in the proposed approach. This idea allows us to determine an appropriate threshold value ψ 0 for the virtual enstrophy in order to reach a good agreement with the experimental breaking points. In our numerical simulations, the breaking position can be defined either at the place where the criterion on ψ 0 is first satisfied or where the wave amplitude starts to decrease. The second position follows closely the first one but takes place nevertheless a little bit later. Given this uncertainty in the definition of the breaking point, some deviations to the experimental values is inevitable. The values of the model parameters were chosen to be able to predict the breaking position with a good accuracy in all cases. The deviation on the abscissa of the breaking point is always smaller than 1 m (compare to the wavelength which ranges from 20 m to 100 m approximately).

The breaking position depends on both the threshold value ψ 0 and on the Reynolds number R but the amplitude variation of the wave after breaking depends on R but not on ψ 0 . The value of C r has an effect only in the front region of the wave where there is a high value of the enstrophy. Afterwards only R has a significant effect on the wave profile. The values of the model parameters can thus be found by comparison of the calculated depth profiles and breaking point positions with the corresponding experimental values. In order to obtain a predictive model, the values of the parameters were sought as functions of the bathymetry and of the initial nonlinearity only. It is important to note that the results presented here were not obtained with an ad-hoc tuning of the parameters in each case but with the predictive laws which were brought out from the comparisons with the experimental measures and which are given below.

A good agreement of the numerical simulations with the experimental depth profiles and breaking positions was found with the following values of the model's parameters. Firstly C r was found to have the same value C r = 0.48 in all cases. Preliminary calculations with this model in the case of roll waves and hydraulic jumps (not yet ready for publication) seem to give also this same value. The value C r = 0.48 could be universal. Secondly as explained above, no breaking criterion is needed for a nonlinearity parameter smaller than 0.05 (apart from the condition of §5.3.1). This means that ψ 0 = 0 if δ < 0.05. In this case, the value R = 6 gives a good agreement on the waves profiles after breaking. There is not a great sensitivity to the value of R however and a value of R slightly different from 6 can be used. Thirdly for higher nonlinearities, the virtual enstrophy threshold was found to depend on the nonlinearity if δ > 0.05. In this expression the dimensionless virtual enstrophy threshold is defined as ψ0 = ψ 0 h * 0 /g. The corresponding value of R for δ < 0.05 is constant and lies between 1.6 and 2.0. The constant value R = 1.7 is used in this section.

A representative series of results is presented in Figure 9 No attempt is made here to improve the dispersive properties. The model is fully nonlinear and has the same dispersive properties as the equations of Green-Naghdi. Its dispersive properties are thus better than those of weakly-nonlinear Boussinesq models. However the dispersive properties can be improved in the same way as those of the Green-Naghdi model as in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] or [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF]. Because in this part the standard dispersive properties of the Green-Naghdi system are kept, a small deviation to the experimental profiles in the shoaling zone can sometimes be observed. For example, in Figure 9 the calculated wave amplitude is smaller than the experimental one at x = 192.36 m just before the breaking point.

Amplitude evolution and breaking depth

The model predictions on the evolution of the wave amplitude are compared to the power-law of [START_REF] Synolakis | Evolution of maximum amplitude of solitary waves on plane beaches[END_REF] a

h b = h 0 h b n ( 95
)
where a is the amplitude of the wave at a position where the still water depth is h 0 , h b is the water depth at the breaking point and n is an exponent depending on the region. [START_REF] Synolakis | Evolution of maximum amplitude of solitary waves on plane beaches[END_REF] showed that four regions exist, two in the shoaling zone before breaking and two after breaking. During the wave propagation shoreward there is first the zone of gradual shoaling where n = -1/4 (corresponding to Green's law), then the zone of rapid shoaling where n = -1 (Boussinesq's law) followed by the zone of rapid decay (n = 4) and by the zone of 95) in each zone including in the fifth zone where we find a value of the exponent n close to 1/4. The curve is a little bit less sharp than the experimental one at the breaking point. This is probably due to the dispersive properties of the model which are identical to the Green-Naghdi equations and which are not optimal as it is explained above. It results that the amplitude is slightly underestimated just before the breaking point. In all other regions, the predicted wave amplitude evolution is very similar to the experimental results of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF].

The breaking depth and thus the breaking position predicted by the model is checked by comparison with the empirical law proposed by Grilli et al. (1997) who defined the slope parameter for solitary waves S 0 = 1.521 tgβ/ √ δ. They found that the breaking types depend on the value of this parameter. According to this criterion the waves studied by [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF] are all of the plunging type since 0.025 < S 0 < 0.30 (the wave with the highest nonlinearity has δ = 0.338 and S 0 = 0.044). However [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF] observed that for S 0 < 0.05 the breaker type is not a pure plunging type but that the spilling breaker prevails. They also found that the properties of amplitude evolution are not significantly affected by the breaker types. Grilli et al. (1997) found that the breaking depth was correctly predicted by the law

h b h * 0 = 0.149 (S 0 /δ) 0.523 (96)
if S 0 < 0.30. On the whole this law was also confirmed by the experiments of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF].

The comparison of the breaking depths predicted by our model with this law are presented in Figure 12. Our numerical results are very similar to the experimental results of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF]. There is an overall agreement with the law (96). The deviations follow exactly the same trends as the experimental measures. In particular, as noted by [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF], the law (96) under-predicts slightly the breaking water depth for S 0 /δ < 2. This small deviation has the greatest value, both in the experiments and in the numerical simulations, for S 0 /δ close to 1. Conversely, for large values of S 0 /δ (between 9 and 10), the numerical results are smaller than the values predicted by the law (96). In the experiments of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF], most values of h b /h * 0 for S 0 /δ between 9 and 10 are also smaller than the predictions of this law. This comparison shows that the model gives an accurate prediction on the breaking depth and position, in a very good agreement with the experiments.

Dependence on the slope

The influence of the slope on the values of the virtual enstrophy threshold ψ 0 and on R is studied from recent experimental works by [START_REF] Fuchs | Solitary impulse wave transformation to overland flow[END_REF] and [START_REF] Hafsteinsson | Solitary wave run-up: wave breaking and bore propagation[END_REF] for various slope angles (β = 1 • , 3 • , 6 • and 11 • ) and a nonlinearity parameter ranging from 0.2 to 0.7. The value of the bottom slope angle for the numerical simulations of the experiments of [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF] is β = 0.95 • which is very close to β = 1 • . The numerical tests confirm that the proposed laws for ψ 0 and R (see §6.1) give also satisfactory predictions for the case β = 1 • of [START_REF] Hafsteinsson | Solitary wave run-up: wave breaking and bore propagation[END_REF].

However numerical tests show that these values are not valid for other slope angles. Since the increase of the virtual enstrophy is faster if R is smaller, a given breaking point can be obtained with different sets of values for ψ 0 and R. The wave breaks later if ψ 0 or R are increased. However given the steepness of the variation of the virtual enstrophy near the breaking point (see Figure 6), the threshold value ψ 0 has a much smaller effect on the breaking position than R. The value of ψ 0 must be strongly modified change significantly the breaking point. It is therefore more practical to keep the same value of ψ 0 as above and to suppose that only R depends on the slope. This choice is also simpler to implement in the numerical scheme because in this case R is locally determined by the bathymetry. Another criterion to determine R is the comparison with the depth profiles of the breaking waves and this method was used in §6.1. However for higher slopes, this method cannot be used. For β = 11 • the wave does not break at all [START_REF] Fuchs | Solitary impulse wave transformation to overland flow[END_REF]. For β = 6 • the wave breaks so close to the shore than no bore is generated and for β = 3 • the bore lasts only for a short moment [START_REF] Hafsteinsson | Solitary wave run-up: wave breaking and bore propagation[END_REF]). According to our numerical tests, supposing that R depends on the slope but not ψ 0 is not only the simplest choice but it works also better.

The assumption that R depends only on the slope and that ψ0 depends only on the nonlinearity as in §6.1 gives immediately the value of R in the case of a horizontal bottom if the value of McCowan (1894) of the maximum height of a soliton is taken (see §5.3.3). With δ = 0.78, the law (94) gives the dimensionless threshold value ψ0 = 0.14. Then the numerical calculation shows that the highest nonlinearity for a non-breaking soliton is δ = 0.78 if and only if R = 0.85, which is thus the value of R for tgβ = 0.

For each nonlinearity parameter of the waves studied by [START_REF] Hafsteinsson | Solitary wave run-up: wave breaking and bore propagation[END_REF], the value of ψ0 is found by (94) and the value of R is sought in order to obtain breaking depths and positions in agreement with the experimental measures. This gives values of R for β = 1 98), a solitary wave on an incline slope of 11 • does not break in accordance with [START_REF] Fuchs | Solitary impulse wave transformation to overland flow[END_REF].

Conclusion

We derived a new model for coastal waves capable of describing shoaling and breaking waves. This one-dimensional model resolves the large-scale turbulence while the small-scale turbulence is modelled through a turbulent-viscosity hypothesis. The three equations of the model are obtained by averaging over the depth the balance equations for mass, momentum and kinetic energy with a shallow-water hypothesis. This approach can be seen as a depth-averaged largeeddy simulation model with a cutoff frequency in the inertial subrange. No assumption is made on the order of magnitude of the nonlinearity. As a result, the model is fully nonlinear and its dispersive properties are identical to those of the Green-Naghdi equations. The model's variables are the fluid depth, the average velocity and the enstrophy, which includes the largescale turbulence effects. The viscous and dispersive terms preclude discontinuities arising and the enstrophy transport equation can be conveniently substituted for the energy equation. The breaking point is characterized by a sudden increase of the enstrophy. In some cases no breaking criterion is necessary to predict the wave breaking and the propagation in the surf zone. In other cases, a breaking criterion must be implemented to prevent a decrease of the wave amplitude before breaking. Most breaking criterions can be used although a breaking criterion based on a virtual enstrophy is proposed. The model is validated in the case of a solitary wave on a mild sloping beach by comparisons with experimental results of the literature, most notably from [START_REF] Hsiao | On the evolution run-up of breaking solitary waves on a mild sloping beach[END_REF]. Numerical simulations were made over a wide range of values of the nonlinearity parameter and of the slope. The three dimensionless parameters are a dissipation coefficient C r which has the constant value 0.48, a turbulent Reynolds number depending only on the slope and a threshold value, if the breaking criterion based on the virtual enstrophy is used, which depends only on the initial nonlinearity parameter of the wave. If the bottom topography and the initial wave conditions are known, the model can thus predict the whole solitary wave propagation.

The continuation of this work will extend the model to the two-dimensional case and will include complex bathymetries, run-up, run-down, wave trains and improved dispersive properties. Further developments are envisaged to model particles and sediment transport for applications in coastal erosion.

Figure 1 :

 1 Figure 1: Notations used in the text.

  21) In this equation, the filtered viscous stress tensor is τ = 2ρνD. The components of the involved symetrical tensors are defined as A r = A r xx e x ⊗ e x + A r xz e x ⊗ e z + A r xz e z ⊗ e x + A r zz e z ⊗ e z and τ = τ xx e x ⊗ e x + τ xz e x ⊗ e z + τ xz e z ⊗ e x + τ zz e z ⊗ e z , where e x and e z are the unit vectors in the Ox and Oz directions.

Figure 2 :

 2 Figure 2: Depth profile of a soliton for F r = 1.095 and h 0 = 1 m if ϕ 0 = 0.05 s -2 (solid line) and if ϕ = 0 (dashed line).

Figure 3 :

 3 Figure 3: Numerical solution of a soliton calculated for F r = 1.095, h 0 = 1 m and ϕ 0 = 0.05 s -2 (solid line) and comparison to the exact solution (crosses) at different times: t = 0 s, t = 10 s, t = 20 s, t = 30 s and t = 40 s.

Figure 4 :

 4 Figure 4: Convergence of the numerical scheme used for the calculation of a soliton for F r = 1.095, h 0 = 1 m and ϕ 0 = 0.05 s -2 .

Figure 5 :

 5 Figure 5: Location of turbulence generation and of vorticity generation for a plunging breaker (in grey). The direction of propagation is from left to right. (a) steepening phase at the beginning of the plunging process. (b) breaking wave after the second splash-up event (from experimental results of Kimmoun & Branger 2007).

Figure 6 :

 6 Figure 6: Evolution of the maximal value of the enstrophy for δ = 0.048 (solid curve) and of the virtual enstrophy for δ = 0.137 (dashed curve). The thin vertical lines show the breaking positions in each case. (a) propagation until the wave breaking. (b) detailed zoom on the initial values of the virtual enstrophy.

Figure 7 :

 7 Figure 7: Experimental setup. A solitary wave of initial amplitude a * in an initial still water depth h * 0 arrives over a sloping beach with an incline angle β. The depth at the breaking point is h b .

Figure 9 :

 9 Figure 9: Calculated variation in time of the wave amplitude (solid curve) and comparison with the experimental results in the case of the trial 54 of Hsiao et al. (2008). The dashed curve is the calculated enstrophy variation.

Figure 10 :

 10 Figure 10: Calculated variation in time of the wave amplitude (solid curve) and comparison with the experimental results in the case of the trial 19 of Hsiao et al. (2008). The dashed curve is the calculated enstrophy variation.

  in the case of the trial 54 of Hsiao et al. (2008) (h * 0 = 2.9 m and δ = 0.086). Three results in the shoaling zone are presented at the abscissa 176.55 m, 187.06 m and 192.36 m followed by three results in the surf zone at the abscissa 202.88 m, 205.57 m and 216.09 m. The breaking point was measured by Hsiao et al. (2008) at 193 m. The breaking criterion is first satisfied at 192.7 m and the amplitude starts to decrease at 193.5 m. In each case the variation of the wave amplitude as a function of time is shown as well as the enstrophy variation (dashed curve). A breaking wave can be easily identified by the presence of a non-zero value of the enstrophy. The results in the case of a wave with a high initial nonlinearity (Trial 19 of Hsiao et al. 2008, δ = 0.338 and h * 0 = 1.2 m) are presented in Figure 10. Only two different representative positions in the surf zone are shown. The asymmetry of the depth profile and the wave amplitude are both correctly described.

Figure 11 :

 11 Figure 11: Evolution of the amplitude for a solitary wave: numerical results (points) and law 95 (solid lines). (a) h * 0 = 1.2 m and δ = 0.338. (b) h * 0 = 2.2 m and δ = 0.152.

Figure 12 :

 12 Figure 12: Comparison of the numerical breaking depth (points) with the prediction by (96).

  • , β = 3 • and β = 6 • besides the value for β = 0.95 • of Hsiao et al. (2008) already known.

Figure 13 :

 13 Figure 13: Numerical results for β = 1 • (•), β = 3 • ( ) and β = 6 • (•) and comparison with the law (96) (solid curve).

Table 1 :

 1 Conditions of the numerical simulations. Trial numbers refer to

	Trial	3		14	15	17	19	21
	h * 0 (m)	1.2	1.2	1.2	1.2	1.2	1.2	2.2
	δ	0.048 0.112 0.227 0.248 0.2875 0.338 0.0195
	Trial	25	31	37	41	43	49	54
	h * 0 (m)	2.2	2.2	2.2	2.2	2.2	2.9	2.9
	δ	0.054 0.069 0.120 0.137 0.152 0.019 0.086

Acknowledgement

We thank Philippe Bonneton and Arnaud Duran for helpful suggestions and discussions. The work was supported by the Service Hydrographique et Océanographique de la Marine (SHOM).