P Bacigaluppi

M Ricchiuto

P Bonneton

Upwind Stabilized Finite Element Modelling of Non-hydrostatic Wave Breaking and Run-up

Keywords: Wave propagation, wave breaking, shock-capturing, stabilized finite elements, SUPG scheme, Boussinesq equations, shallow water equations, wet/dry fronts, wave breaking model

We discuss a methodology to model the propagation, wave breaking and run-up of waves in coastal zones. We represent the different coastal phenomena through the coupling of nonlinear shallow water equations with the extended Boussinesq equations of Madsen and Sørensen. Each of the involved equations has a major role in describing a particular physical behaviour of the wave: the latter equations permit to model the propagation, while the non-linear shallow water ones lead waves to locally converge into discontinuities. We start from the third-order stabilized finite element scheme for the Boussinesq equations, developed in a previous scientific work (Ricchiuto and Filippini, J.Comput.Phys. 2014) and develop a non-linear variant, and detach the dispersive from the shallow water terms. A shock-capturing technique based on local non-linear mass lumping that permits in the shallow water regions to degrade locally the scheme to a first-order one across bores (shocks) and dry fronts is proposed. As for the detection of the breaking fronts, the shallow water areas, this involves physics based breaking criteria. We present different definitions of the breaking criterion, including a local implementation of the convective criterion of (Bjørkavåg and H.

Upwind Stabilized Finite Element Modelling of Non-hydrostatic Wave Breaking and Run-up

Résumé : On décrit une approche pour la simulation de la propagation et déferlement des vagues en proche côte basée sur la couplage entre les équations de Boussinesq ameliorées de Madsen and Sørensen, pour la propagation, et les equations Shallow Water, pour le déferlement et le runup. La contruction de ce mod le hybride passe d'abord la proposition une variante non-linéaire du schéma élément finis stabilisé de (Ricchiuto and Filippini, J.Comput.Phys. 2014) capable de résoudre les chocs de maniere monotone. Cela est obtenu par un operateur locale de condensation de la matrice de masse qui réduit le schéma de ( Mots-clés : Propagation de vagues, déferlement, capture de choc, éléments finis, stabilisation, equations de Boussinesq, Shallow Water, inondation, runup, détection fronts déferlants

Introduction

The aim of the present paper is to present a model able to describe the propagation and transformation of waves in near shore areas. The focus will be on the phenomena of wave breaking, shoaling and run-up and run-down. The characteristics of the waves such as refraction, diffraction, reflection and non-linear interaction, will have to be accurately described and, in particular, a new approach to detect wave breaking, which corresponds to the steepening of the wave's crest which turns into a roller, characterised by a turbulent motion, will be presented.

To this purpose, a non-hydrostatic model has to be considered, such as Peregrine's equations [START_REF] Peregrine | Long waves on a beach[END_REF], which were the first capable to describe through two dependent variables, the surface displacement and the depth averaged horizontal velocity, the refraction and diffraction effects modifying the shape of a wave in presence of variable depths [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF][START_REF] Wei | A coastal process model based on time-domain Boussinesq equations[END_REF]. Following Peregrine's model, various authors derived similar and improved equations, such as, for example, Nwogu who derived a similar model considering the velocity at a certain depth [START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF] and Madsen and Sørensen who extended Peregrine's equations to intermediate depths through the introduction of terms of high order, resulting thus formally to zero for the accuracy of the model [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF][START_REF] Wei | A coastal process model based on time-domain Boussinesq equations[END_REF]. Both Nwogu and Madsen and Sørensen enhanced Boussinesq equations, although having different descriptions for the dispersive terms, have the dispersive relations derived by a Padé approximation of the full linear dispersion relation, giving an accuracy of O(σ 4 ) [START_REF] Wei | Time-dependent numerical code for extended Boussinesq equations[END_REF], with σ the ratio between the water depth and the wavelength. Further models taking into account, beside the dispersive nature of the wave, also the non-linear one, have been derived for example by Serre Green-Naghdi [START_REF] Kim | A derivation of the Green-Naghdi equations for irrotational flows[END_REF][START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF], which is constituted by a hyperbolic part, corresponding to the non-linear Shallow Water equations, plus dispersive terms, and in the FUNWAVE model from Wei's fully non-linear Boussinesq equations [START_REF] Kirby | FUNWAVE 1.0 Fully nonlinear Boussinesq wave model documentation and user's manual[END_REF][START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]. These models alone are not able to approximate the wave breaking phenomena, as the dispersive terms become negligible in this case, having the non-linear terms a predominance due to the dissipating energy characterising the formation of a roller.

Starting from Filippini's model [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF] given by the enhanced Boussinesq equations of Madsen and Sørensen, in the present paper it has been chosen to perform the wave breaking with a hybrid equation.

In the following, the development of a monotone Shallow Water equations management and of a criterion able to detect breaking waves will be presented. When referring to monotone Shallow Water equations, monotone systems are being meant, and are obtained, in general, through the application of flux-limiters, as in Gasdynamics, which are mainly slope-smoothers, capable of combining a low order method for treating discontinuities and high order ones, when the solution is smooth, so that, when occurring spurious oscillations, this smoothing tool provides for their elimination from the solution. Sweby and LeVeque in [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF] and [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], respectively, have resumed some of the main techniques, as the van Leer flux-limiter, used by many authors in literature as in [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF][START_REF] Shiach | A temporally second-order accurate Godunov-type scheme for solving the extended Boussinesq equations[END_REF]. Other examples are the Van Albada-Van Leer edge-based non-linear slope limiter [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF][START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF], the Minmod limiting function [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF][START_REF] Shi | A high-order adaptive timestepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF], the Superbee limiter [START_REF] Roeber | Schock-capturing Boussinesqtype model for nearshore wave processes[END_REF], also for the two-dimensional case through a two-dimensional slope limiter [START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF], all of which have been dealing with the fluxes represented through finite volume methods. A further example, with the use of a Residual Distribution (RD) scheme, is given by [START_REF] Abgrall | Toward the ultimate conservative scheme : Following the quest[END_REF][START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for steady flow problems[END_REF][START_REF] Ricchiuto | Stabilized Residual Distribution for shallow water simulations[END_REF], where a nonlinear limiter, similar to a multidimensional version of the Minmod, in order to capture shocks, has been considered. Generally, almost all, when not all, authors apply the limiting function to the fluxes, while in this paper a different approach is presented, based on the mass-lumping, with limiters applied to the mass terms, which for the 1D case gives excellent results. Moreover, the scheme has to represent shocks correctly, preventing the solution of being entropy violating, [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF][START_REF] Pelanti | A review of entropy fixes as applied to Roe's linearization[END_REF], with a good solution, for example, in the Harten's family of entropy fix [START_REF] Harten | Self-adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF], applied like proposed in [START_REF] Madrane | Entropy stability of Roe-type upwind finite volume methods on unstructured grids[END_REF].

Inria

Beside the treatment of wave breaking, the second aim is to be able to approximate the motion of run-up of waves on shores, which requires the design of the model for the treatment of dry areas and especially the presence of wet/dry fronts [START_REF] Brufau | A numerical model for the flooding and drying of irregular domains[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][START_REF] Fuhrman | Simulation of nonlinear wave run-up with a high-order Boussinesq model[END_REF][START_REF] Kurganov | A second-order well-balanced positivity preserving centralupwind scheme for the Saint-Venant system[END_REF][START_REF] Lynett | Modeling wave runup with depth-integrated equations[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF]. In particular, a good description which considers also the modifications to perform to the terms related to the bathymetry changes, in agreement with the adopted scheme was found in [START_REF] Castro | The numerical treatment of wet/dry fronts in shallow flows: Application to one-layer and two-layer systems[END_REF].

The accomplishment of treating wave breaking has been up to now a riddling quest, with many different approaches. Historically, in the 90s, two main different criteria have been developed [START_REF] Okamoto | The relative trough Froude number for initiation of wave breaking: theory experiments and numerical model confirmation[END_REF].

The first was characterised by the analysis of the information of a full phase of the wave, such as the period and wave length, which brought to consider this type as a phase-averaged, resulting though in a very expensive and not very stable approach. Some time later, a second type, a phase-resolving approach was developed [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF][START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking and run-up. I: 1D[END_REF][START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF], which resulted more stable and accurate, beside the less expensive computational costs, as it used the information obtained by the solution, such as e.g. the free surface elevation. In this second approach, two main families can be further distinguished in literature: the hybrid ( [START_REF] Bonneton | Modelling of periodic wave transformation in the inner surf zone[END_REF][START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF]) and the eddy viscosity ( [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF][START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking and run-up. I: 1D[END_REF][START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF][START_REF] Sørensen | Boussinesq-type modelling using an unstructured finite element technique[END_REF]) formulation. This work considers the phase-resolving hybrid approach. Two breaking criteria are tested. The first is the one proposed by Tonelli and Petti in [START_REF] Tonelli | Hybrid finite volume-finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Tonelli | Finite volume scheme for the solution of 2D extended Boussinesq equations in the surf zone[END_REF][START_REF] Tonelli | Simulation of wave breaking over complex bathymetries by a Boussinesq model[END_REF][START_REF] Tonelli | Shock-capturing Boussinesq model for irregular wave propagation[END_REF][START_REF] Tonelli | Numerical simulation of wave overtopping at coastal dikes and low-crested structures by means of shock capturing Boussinesq model[END_REF], based on the use of the non-linearity coefficient as a trigger for wave breaking. The second criterion is the new hybrid criterion of Kazolea et al. [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF][START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF]. In this case wave breaking is detected using the wave's slope and the wave's vertical velocity.

We present different definitions of the breaking criterion, including a local implementation of the convective criterion of Bjørkavåg and H. Kalisch, Phys.Letters A 2011), and discuss in some detail the implementation of the shock capturing technique. Moreover, a new criterion based on a local implementation of the convective criterion in [START_REF] Bjørkavåg | Wave breaking in Boussinesq models for undular bores[END_REF], is introduced. The new idea is to consider the ratio between the crest's velocity and wave's celerity as a monitor for wave breaking, leading to the introduction of a physical Froude number given by the ratio of these two quantities.

The remainder of the paper is organised as follows: Section 2 is devoted to the presentation of the model equations and its discretisation, while in section 3 the techniques used to solve the Shallow Water equations numerically are presented. Further, in chapter 4 the proposed hybrid approach, together with the adopted wave breaking criteria are discussed. An extensive numerical validation and physical investigation of the different wave breaking strategies is reported finally in section 5.

Modelling approach and equations

Consider, starting from Peregrine's equations and following the procedure discussed in [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF], the one-dimensional depth-averaged system of the Madsen and Sørensen extended Boussinesq equations that read:

         ∂ t η + ∂ x q = 0 ∂ t q + ∂ x (uq) + gH∂ x η + gHC f u -Bh 2 ∂ xxt q -βgh 3 ∂ xxx η+ - 1 3 h∂ x h∂ xt q -2βgh 2 ∂ x h∂ xx η = 0 (1a) (1b)
This system of equations is written (w.r.t figure 1) in terms of conservative variables η(x, t) and q(x, t), that are respectively the surface elevation and the discharge defined as q = Hu, which measures the flow rate through a unit vertical surface. The depth averaged velocity is u, while H(x, t) = η(x, t) + h(x) corresponds to the water depth with h(x) the depth of the bathymetry w.r.t the still water level.

Figure 1: Depth averaged equations : notation

The parameters of β and B characterize the optimization of the dispersion properties of a linearized model with respect to the wave theory of Airy studied extensively in [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF]. In particular, for (1) the values are β = 1 15 and B = β + 1 3 , allowing to match a fourth order Padé development of the dispersion parameter [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics[END_REF]. The equations 1 proposed in [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF] approximate in an accurate way the physics of the propagation of a wave in the near shore region, when waves are relatively long, with a ratio water-depth over wavelength σ 2 ≪ 1. This model is no more valid when approaching the shoreline and the wave steepens causing non-linear effects becoming dominating possibly ending up in a wave breaking (ǫ = A/d ∼ 1). Several approaches exist, see [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF] for a recent review. The extensive study of [START_REF] Bonneton | Modelling of periodic wave transformation in the inner surf zone[END_REF], indicate that the behavior of a non-hydrostatic wave in a breaking region can be modelled by the dissipation across non-linear discontinuities of hyperbolic models such as the shallow water equations. When considering a solitary wave running on an uneven bathymetry with A/d ∼ 1, what would be expected for a real wave, would be the formation of a breaking front where the wave rolls over the submerged shelf, forming a socalled roller. The roller in fact, represents the moment when the wave's crest breaks, after having steepened and grown and, in particular, it starts in the moment that a given point overwhelms a point situated lower on the surface of the moving front. When computing this situation through the approximation given by the Madsen and Sørensen equations, the water, that would physically form a crest separates from the major wave, resulting in multiple waves which are less steeper than the first wave and which become more distinct, during the continuing of the propagation, as each has its own frequency dispersion and propagates with its own velocity. Indeed, when this behavior occurs, non-linear effects should predominate on the dispersive ones. A good model to approximate this occurrence is given by the non-linear shallow water equations

∂ t η + ∂ x q = 0 ∂ t q + ∂ x (uq) + gH∂ x η + gHC f u = 0 (2a) (2b)
These equations are in fact able to dissipate energy, through the formation of a hydraulic bore (a shock), allowing to approximate what physically happens during the breaking. Due to the extreme high dissipative behavior of these last modelling equations, the non-linear shallow water model is inadequate for the description of the motion of the wave on a regular bed. The physics of the propagation, wave breaking and run-up of a wave should be consequently modelled through a coupling of both the Madsen and Sorensen's equations and the sole shallow water ones, recognizing if a wave breaking is occurring and in that case switch from the dispersive dominated equation of Boussinesq to the non-linear shallow water model.

Inria 3 Non-Linear shallow water equations scheme implementation

We start by considering the frictionless non-linear shallow water equations (2) that read

   ∂ t η + ∂ x q = 0 ∂ t q + ∂ x (uq) + ∂ x (g H 2 2 ) -gH∂ x h + gHC f u = 0 (3a) (3b)
The discretisation is accomplished through an upwind stabilized Galerkin finite element method (SUPG) for the space, and the Crank-Nicholson method for the time. These choices are based on the previous investigation in [START_REF] Ricchiuto | Upwind residual discretization of enhanced boussinesq equations for wave propagation over complex bathymetries[END_REF], where different residual based techniques have been analyzed. Using the fact that ∂ t h = 0, we can rewrite the hyperbolic non-linear system of partial differential equations in a vectorial conservative form, with the vector of unknowns W = W(x, t) = [H(x, t) q(x, t)] T = [η q] T depending on x and t, the flux vector given by F (W) = [q (uq +g H 2

2 )] T and the source term as

S = S(W, x) = [0 (-gH∂ x h)] T ∂ t W + ∂ x F (W) = -S (4)
The ( 4) with initial and boundary conditions represents a well posed problem and the equivalent quasi-linear form, which might be considered in case of smooth solutions is given by

∂ t W + ∂F (W) ∂W ∂ x W = -S (5)
Considering the Jacobian matrix A = ∂F (W) ∂W , given c 2 = gH, its eigenvalues are:

Λ = u -c 0 0 u + c (6)
and consequently the right and left eigenvectors :

R = 1 1 u -c u + c (7)
L = 1 2c u + c -1 c -u 1 (8)
The analytical form of the Jacobian matrix can be shown to be :

A = RΛL = 0 1 c 2 -u 2 2u (9)
This characterization of the PDE assumes a certain smoothness of the solution W and of its derivatives. However, it is known that non-linear conservation laws as (4) develop discontinuous solutions in finite time. In order to be able to recognize more general admissible solutions, and to reduce the regularity constraints on W, the weak form of 4 is considered. It is assumed that W is very smooth w.r.t. time, so that for example for a fixed x ∈ [0, L], we have

W(t, x) ∈ C 1 ([0, T ]).
this allows to consider a point-wise-in-time weak formulation in space which reads:

L 0 v∂ t Wdx - L 0 F (W)∂ x vdx + L 0 vSdx + [v F ] L 0 = 0 (10)
Here, as in the space-time formulation, the term F allows to introduce the boundary conditions on the left and right ends of the spatial domain. In the following periodic boundary conditions will be assumed, which means that the term [v F ] L 0 will be considered zero.

Following now the standard derivation of the C 0 Galerkin finite element method [START_REF] Quarteroni | Modellistica Numerica per Problemi Differenziali[END_REF], we introduce a discretisation of the spatial domain Ω = [0, L] consisting of N constant cells, with the space between two nodes composing a cell being ∆x. On this grid we consider the space of all the piecewise linear polynomials V h . This space is spanned by a basis of functions ϕ i (x) (i=0, ...,N), which are piecewise linear and continuous on the constant partition of the spacial domain [START_REF] Quarteroni | Modellistica Numerica per Problemi Differenziali[END_REF] and with the property of being zero on the nodes i-1 and i+1 as well as on the rest of the nodes except for the elements between i-1 and i+1 and the node i, where it assumes the value of one. In particular, we have that

W h (x, t) ≈ N +1 j=1 ϕ j (x)W j (t) (11)
Requiring [START_REF] Castro | The numerical treatment of wet/dry fronts in shallow flows: Application to one-layer and two-layer systems[END_REF] to be valid for an

arbitrary element v = v h ∈ V h , is equivalent to require that L 0 ϕ i ∂ t W h dx - L 0 ∂ x ϕ i F (W h )dx + L 0 ϕ i Sdx = 0, ∀ i (12)
The integration on the spatial domain Ω can be consequently split into a sum of integrals over each single element.The i-th equation can be recast as

i+1 j=i-1 xi+1 xi-1 ϕ i ϕ j dx dW j dt - i+1 j=i-1 xi+1 xi-1 ϕ j ∂ x ϕ i dx F (W j ) + xi+1 xi-1 ϕ i S dx = 0 (13)
In order to solve (13) the Simpson's quadrature formula [START_REF] Quarteroni | Matematica Numerica[END_REF] is applied. The resulting terms corresponding to the integral of the first term are the components of the Galerkin mass matrix, and the obtained system is an ordinary differential equation in time, where the equation for each node of the spatial domain results given by a centered difference scheme:

∆x 6 d(W i-1 ) dt + 2∆x 3 d(W i ) dt + ∆x 6 d(W i+1 ) dt + 1 2 (F i+1 -F i-1 ) + + ∆x 2 (S i+ 1 2 -S i-1 2 ) = 0 (14)
where

S i-1 2 = xi xi-1 ϕ i S dx
and similarly for S i+ 1 2 .

As widely discussed in [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF], the accuracy of the Galerkin finite element method is very good when considering the Boussinesq equations. However, the use of the Galerkin scheme in the hyperbolic limit of the shallow water equations, will lead to an unstable discretization, unless the time integration is properly chosen [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF][START_REF] Quarteroni | Modellistica Numerica per Problemi Differenziali[END_REF]. These instabilities will manifest as unbounded Inria spurious modes (oscillations).

In order to eliminate these spurious modes stabilizing terms are introduced. An example of a scheme using such a stabilization mechanism is the Streamline Upwind-Petrov Galerkin weighted residual method (SU /PG). The idea behind this method lays in the choice of the test functions, which are considered equal to the finite element basis plus a streamline upwind perturbation. This means that the element corresponding to the upstream term will weight more than the downstream one: ϕ toti = ϕ i + δϕ i . Due to the intrinsic hyperbolic nature of the problem, the streamline upwind weighting operator is computed through the Jacobian matrix A and is defined as:

δϕ i = A ∂ x ϕ i τ (15)
with τ a scaling parameter with the dimension of a time associated to the characteristic speeds, which is chosen using the classical expression [START_REF] Quarteroni | Modellistica Numerica per Problemi Differenziali[END_REF]

τ ± = 1 2 ∆x |λ ± |
Introducing the residual of equation ( 4):

R(W h ) = ∂ t W h + ∂ x F h + S h (16)
such that R(W) = 0 if W is the exact solution of the scheme can be compactly written as:

Ω ϕ i Rdx + Ω δϕ i Rdx = 0 (17)
Defining

sign(A) = R λ - |λ -| 0 0 λ + |λ + | L (18)
and noting that ∂

x ϕ i = 1/∆x if x ∈ [x i-1 , x i ], while ∂ x ϕ i = -1/∆x if x ∈ [x i-1 , x i ],
the (SU /PG) method can thus be written as the centered scheme obtained by the Galerkin finite element method plus upwinding term related to the adjacent left and right cells i -1 2 and i + 1 2 of node i :

Galerkin + i i-1 1 2 sign(A)Rdx - i+1 i 1 2 sign(A)Rdx = 0 (19)
Using a local linearisation to evaluate the Jacobian A in the stabilization terms, allows to recast the scheme as the third upwind stabilized finite element scheme [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF]:

∆x 6 dW i-1 dt + 2∆x 3
dW i dt + ∆x 6
dW i+1 dt + + 1 2 (F i+1 -F i-1 ) + ∆x 2 (S i+ 1 2 + S i-1 2 )+ + sign(A i-1 2 ) 2 ∆x dW i-1/2 dt - sign(A i+ 1 2 ) 2 ∆x dW i+1/2 dt + + sign(A i-1 2 ) 2 (F i -F i-1 ) - sign(A i+ 1 2 ) 2 (F i+1 -F i )+ + sign(A i-1 2 ) 2 i i-1 Sdx - sign(A i+ 1 2 ) 2 i+1 i Sdx = 0 (20)

RR n°8536

For the time discretisation, the second-order implicit method of Crank-Nicholson was chosen, according to the results obtained in Filippini's thesis [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF]. Differently from [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF], however, here we have used a variable time step, computed using a "user defined" Courant-Friedrichs-Lewy (CFL) number ν , and according to the formula

ν = a △ t △x (21)
with a = max(λ j ), the spectral radius of the Jacobian matrix. At each time step, the maximum of the eigenvalues defined as λ = |u| + √ gH will be computed and the ∆t computed from the CFL formula.

∆t = ν∆x max i (|u i | + √ gH i ) (22)
The time will be updated at each time step with the new computed ∆t. Denoting with n the time step, the applied time integration method is given by:

u n+1 -u n △t + f ( u n+1 + u n 2 ) = 0 (23)
This leads to a non-linear system of algebraic equations of the type R(W) = 0, which have to be solved by some iterative method. This is accomplished through the Newton's iterative method, which, introducing the solution vector W = η n+1 q n+1 T , and the initial state of the iterations W 0 = (η n q n ) T , and W k the solution computed at iteration k, can be written as :

M (W k )(W k+1 -W k ) = -R(W k |W 0 ) (24)
with M the Jacobian matrix defined as M = ∂R(W k |W0)

∂W

. From now on we will refer to -R(W k |W 0 ) as to the Right Hand Side (RHS). Resuming the steps to be done in order to obtain a solution at a generic time step are the following:

• an physical initial condition at time t = 0, W 0 = [η 0 q 0 ] T is given • a cycle is started which considers the the time iteration, and which is exited only in case the maximum time has been reached

• the initial solution is set as initial guess:

W 0 = W n • F (W 0 |W n ) is computed • the mass matrix M n = M (W 0 ) = ∂R(W0|W0) ∂W is computed
• the Newton iterative method is applied in a while cycle that will continue until either the convergence R(W k |W 0 ) < ǫ or the maximum number of iterations k max is reached, and at each iteration compute :

W k+1 = W k -M -1 n F (W k |W n )
• the new time step is computed through the maximum eigenvalue with the formula [START_REF] Harten | Self-adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF].

As this represents the most intensive computational part of the entire solution, the Jacobian has been taken as frozen through each time step, with this choice being a compromise between keeping the matrix constant through all time iterations, and the one that computes it at each step also of the convergence cycle.

Inria

Shock Capturing Improvement

In correspondence of discontinuities the method has been modified by introducing a lumping on the mass matrix. This has been achieved via a discontinuity sensor δ = δ(u) that switches on and off the terms of the mass. The mass-lumping on the system can be obtained starting from:

Ω ϕ i ∂ t W h dx - Ω F h ∂ x ϕ i dx + Ω ϕ i S h dx + sign(A i-1 2 ) 2 Φ i-1 2 - sign(A i+ 1 2 ) 2 Φ i+ 1 2 = 0 (25)
The first terms can be rewritten as

∆x dW i dt + Ω ϕ i ∂ t W h dx -∆x dW i dt = ∆x dW i dt + ∆x 6 dW i-1 dt - ∆x 3
dW i dt + ∆x 6 dW i+1 dt (26)
while the term of the Galerkin part:

Ω ϕ i ∂ x F (W h ) ≈ F (W i+1 ) -F (W i-1 ) 2 (27)
The condensed mass matrix which is represented by the terms in [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF] will have on the node of the diagonal the sum of all the elements corresponding to the nodes on the same line, and all other non-diagonal elements will have to be zero. Considering for example a single cell, in order to do mass lumping, a local smoothness monitor δ is introduced in [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF]. According to the value assumed by the parameter, in case it equals zero a scheme of the first order is adopted, else the high order scheme is assembled. Thus delta in equation [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF], with the introduction of ( 26) and ( 27) leads to :

∆x dW i dt + δ i-1/2 { ∆x 6 [ dW i-1 dt - dW i dt ] + ∆x 2 sign(A i-1 2 ) dW i-1 2 dt } + δ i+1/2 { ∆x 6 [ dW i+1 dt - dW i dt ] - ∆x 2 sign(A i+ 1 2 ) dW i+ 1 2 dt } + F (W i+1 ) -F (W i-1 ) 2 + ∆x 2 (S i-1 2 + S i+ 1 2 )+ + sign(A i-1 2 ) 2 (F i -F i-1 ) - sign(A i+ 1 2 ) 2 (F i+1 -F i )+ + sign(A i-1 2 ) 2 i i-1 S h - sign(A i+ 1 2 ) 2 i+1 i S h = 0 (28)
If no mass lumping is performed the resulting scheme is the Streamline Upwind Petrov-Galerkin scheme (SUPG), which is third order accurate, as seen in [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF]. In case when δ assumes a zero value, a first order Upwind Flux-Splitting Scheme is obtained. In the scalar case, this scheme is local extremum diminishing (LED), in particular it can be recast as

du i dt + i =j c ij (u i -u j ) = 0 (29)
where all the c ij s are positive. The property of the LED as explained in [START_REF] Deconinck | Residual Distribution schemes: foundation and analysis[END_REF] states that the local maximum of the solution has the characteristic of not increasing while the minimum will not decrease. More particularly, as we have shown in this section, when integrated with the Crank-Nicholson scheme in time, the solution obtained with this scheme enjoys a discrete maximum principle. This means that the first order scheme is able to preserve the monotonicity of the solution. For the system of the shallow water equations, this is not possible to be said, as in general the property of a system are not possible to be analyzed. In practice, however, the first order upwind scheme shows a monotonicity preserving behavior. Moreover, we will also show that it preserves the non-negativity of the depth in flooding areas. Note also that, the mass lumping on the Galerkin scheme in case of the chosen function space P 1 introduces a degree of dissipation, but in general does not imply a reduction of the asymptotic overall accuracy which remains above second order, as shown in [START_REF] Ricchiuto | Explicit Runge-Kutta Residual Distribution schemes for time dependent problems: Second order case[END_REF], where a similar mass lumping has been performed on a residual distribution discrete approximation. However, for the SUPG scheme mass lumping allows to go from a linear third order scheme (in 1D) to a linear monotonicity preserving first order one. Note that the discontinuity sensor is normally applied to the reconstruction used to evaluate the fluxes, or on the fluxes themselves. Our idea is instead to control the form of the mass matrix, while maintaining the form of the discretization of the spatial derivatives. This works very well as a shock-capturing technique in the one dimensional case. In two dimensions, a more complex construction will be necessary to obtain a non-oscillatory discretization starting from a high order Petrov-Galerkin scheme.

Implemented Limiters

To define δ, we have tested different kinds of limiters which have been proposed in literature, normally implemented as slope-limiters or flux limiters. In addition, a different typology of sensor, based on the smoothness of the dependent variable, which in our case could be either the free surface height η or its motion q (or both) has been introduced. In general, slope-limiters are designed to avoid, in case of a discontinuity, possible numerical non-physical oscillations by smearing out the solution, taking, for example, the minimum between two adjacent cells, so that in case of a peak it is possible to neglect it, maintaining a high order accuracy [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF].

It is also of great importance, not to smooth exaggeratedly the solution, as the discontinuities have still to be identifiable and sharp. Three different type of smoothness sensor functions have been implemented in the code : the Monotonized Central-difference limiter [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], the Superbee limiter [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], and a smoothness sensor we devised. In all three cases, a sensor δ i is obtained for each node. In each cell we set δ i+1/2 = min(δ i , δ i+1 ) and δ i-1/2 = min(δ i -1, δ i ).

Monotonized Central-difference Limiter

The monotonized central-difference limiter (MC) is, according to [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], highly used and compatible with many kinds of problems. It compares the central difference due to the method of Fromm with two times each slope of the adjacent cell of node i.

σ i = min( η i+1 -η i-1 2∆x , 2 η i -η i-1 ∆x , 2 η i+1 -η i ∆x ) (30)
Then a control is performed, in order not to smooth each node, but only where the |σ i | > 0.5, in case that |σ i | > 1, δi = 0, which means that the mass matrix is lumped, else, in case 0.5

< |σ i | < 1, δi = 1 -|σ i |.
Moreover, in order to not introduce brutal changes in the slope, a further smoothing is done on δ i , in order to switch gradually from the complete mass lumping to a non condensed one.

δ i = 1 4 δi-1 + 1 2 δi + 1 4 δi+1 (31)
Inria .

Figure 2: The Smoothness function of the Superbee and the MC limiter are compared and are shown how they behave in confront to the shaded area that represents the total variation diminishing [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF]. Note that the Φ corresponds to our σ, while θ is the argument that is considered in the minimum function of σ

Superbee Limiter

The Superbee limiter, introduced by Roe, is capable to maintain a second-order accuracy. Similarly to the MC one, it compares the slopes of the cells of the adjacent nodes of i, and, according to [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF] it takes the maximum of the minimum of the slope of right cell and twice that of the left one and that of twice the right one and the left one. The original Superbee limiter has been adapted to the scheme modifying slightly the computation of σ. Instead of choosing between the minimum of the right and twice the left slope and the other way round, firstly, the ratio of the left and right cell is taken twice and, then, the minimum between this and one is chosen. Note that a variable, ǫ u of the order of 10 -6 has been introduced, so that in case the slope of the right cell is constant, no operation between invalid operands is done.

σi = min(2 η i -η i-1 η i+1 -η i + ǫ u , 1) (32)
This permits to choose the σ in the maximum between zero and the σ of [START_REF] Kurganov | A second-order well-balanced positivity preserving centralupwind scheme for the Saint-Venant system[END_REF].

σ i = max(0, σi ) (33)
Finally, in order not to introduce sudden variances and in order to be able to smooth correctly, the minimum of sigma of the three nodes composing the two implied cells is chosen.

δ i = min(σ i-1 , σ i , σ i+1 ) (34)

Smoothness Sensor

Since in our applications we want to preserve as much as possible smooth extrema, while still avoiding oscillations near strong discontinuities, we have considered a different type of limiter based on smoothness analysis. The main idea is to consider two different approximations of a derivative of the solution. To fix ideas, let us consider the case of the second order derivative, although the first order derivative would be similar. We consider two approximations :

(∂ xx η) O4 = - 1 12∆x 2 η i+2 + 4 3∆x 2 η i+1 - 5 2∆x 2 η i + 4 3∆x 2 η i-1 - 1 12∆x 2 η i-2 RR n°8536
and

(∂ xx η) O2 = 1 ∆x 2 η i+1 - 2 ∆x 2 η i + 1 ∆x 2 η i-1
Clearly for a smooth variation of η we have

(∂ xx η) O4 = ∂ xx η + O(∆x 4 ) , (∂ xx η) O2 = ∂ xx η + O(∆x 2 )
On the other hand, for a discontinuity of finite amplitude ∆η around node i we have

(∂ xx η) O4 ≈ 15∆η 12∆x 2 , (∂ xx η) O2 ≈ ∆η ∆x 2
So, in the smooth case we have

|η i+2 -4η i+1 + 6η i -4η i-1 + η i-2 | 12∆x 2 = | (∂ xx η) O4 -(∂ xx η) O2 | = O(∆x 2 )
while across a discontinuity

|η i+2 -4η i+1 + 6η i -4η i-1 + η i-2 | 12∆x 2 = | (∂ xx η) O4 -(∂ xx η) O2 | ≈ ∆η 4∆x 2
Introducing a parameter ǫ = 10 -3 , to prevent obtaining invalid solutions in case the difference between the nodes is zero, the following formula is used:

σ i = min(1, ǫ + |ηi-ηi-1| ∆x + |ηi-ηi+1| ∆x ǫ + |ηi+2-4ηi+1+6ηi-4ηi-1+ηi-2| 12∆x 2 ) (35)
In the smooth case, while the numerator in the right slot is bounded, the denominator is of O(∆x 2 ), so the value of 1 will be retained. Across a discontinuity, the numerator gives ∆η/∆x, while the denominator gives approximately ∆η/4∆x 2 , so σ i = 4∆x ≪ 1. In order to maximize the smoothness of the numerical solution, and to avoid abrupt variations of the sensor we then set

δ i = 1 4 δi-1 + 1 2 δi + 1 4 δi+1 , δi =    σ i if σ i < 0.5 1 otherwise (36)

Preliminary Numerical Results

We present a preliminary verification and comparison of the analytical solution and the limiters implemented on a Riemann problem with initial solution of η l = 0[m] and η r = -0.5[m] on a constant bathymetry of 1[m] height, a spatial domain [0,100], and a time span of 30s with a CFL of 0.3. The solution has been compared at time t = 20/ g/h 0 . For the analytical solution we have referred to [START_REF] Ambrosi | Approximation of shallow water equations by roe's riemann solver[END_REF][START_REF] Delestre | Swashes: a compilation of shallow water analytic solutions for hydraulic and environmental studies[END_REF][START_REF] Stoker | Water Waves: The Mathematical Theory with Applications[END_REF], implementing it as reported in [START_REF] Delestre | Swashes: a compilation of shallow water analytic solutions for hydraulic and environmental studies[END_REF]. These results demonstrate that the mass lumping increases the quality of the solution. The figures in figure 3 show how numerical oscillations are eliminated and the solution results smoothed through the application of the correction to the scheme.

Inria (a) (b) (c) (d)

Entropy Fix

As thoroughly discussed in Leveque's book [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF] the major problem of the computed solution is that the second law of the thermodynamics, which states, that the entropy variation is major or equal to zero is violated, as in the case of an expansion, the scheme leads to a diminishing entropy, which is absolutely non physical. Moreover, the proposed scheme is in conservation form and the limit solution is the weak one, which means that the only violation possible is that of the rarefaction, where the Rankine-Hugoniot condition is satisfied, while the entropy inequality not [START_REF] Harten | Self-adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF]. In [START_REF] Pelanti | A review of entropy fixes as applied to Roe's linearization[END_REF], the entropic solution is seen as a limit solution when thermal conductivity and viscosity disappear. In our case, which is one-dimensional, the expansion might be seen as a discontinuity, and in particular, this discontinuity has the property in the transonic case of remaining captured in a cell, which means mainly, that the initial discontinuity is not able to be dissipated by the scheme and has therefore to be smoothed. In literature various methods have been studied in order to work out this problem, for example introducing numerical viscosity or doing the so-called entropy fix. In our case we will focus on the entropy fix and in particular on one of the Harten's family, namely the second Harten and Hyman's (1983) [START_REF] Harten | Self-adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF] entropy fix. Harten's entropy correction method consists in preventing the eigenvalues of the Jacobian matrix A from becoming too small or to vanish. The idea is that since the matrix A controls the dissipative effects, if the eigenvalues are too small, as in the problem considered , the discontinuity is not correctly dissipated. The entropy fix has the task of identifying the small eigenvalues and substitute them with another value. In order to find the region of interest a control over the eigenvalues of the Jacobian matrix A is performed, allowing to check if if |λ| < ǫ, with ǫ a small positive number. Following the developments of [START_REF] Harten | Self-adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF][START_REF] Kermani | Modified entropy correction formula for the Roe scheme[END_REF][START_REF] Pelanti | A review of entropy fixes as applied to Roe's linearization[END_REF] and [START_REF] Madrane | Entropy stability of Roe-type upwind finite volume methods on unstructured grids[END_REF], we have implemented the entropy correction. The criteria that states if the value is to small or not, is given by an ǫ for each absolute value of the computed eigenvalue, given by:

ǫ 1 = 0.5 max(0, max(λ 1a -λ 1b , λ 1b -λ 1a )); ǫ 2 = 0.5 max(0, max(λ 2a -λ 2b , λ 2b -λ 2a )); (37)
And we then define for j = 1, 2 [START_REF] Harten | Self-adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF][START_REF] Kermani | Modified entropy correction formula for the Roe scheme[END_REF][START_REF] Pelanti | A review of entropy fixes as applied to Roe's linearization[END_REF]:

v j =        |λ j | if | λj | ≥ ǫ j λ2 j 2ǫ j + ǫ j 2 otherwise (38)
modifying the sign(Λ):

Λ =     λ 1 v 1 0 0 λ 2 v 2     (39)
Rewriting the dissipation of the scheme in a non conservative form |sign(A)|∆W , which is not eliminated in case of a stationary expansion and permitting thus to dissipate it correctly. In particular, starting from [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF], considering the upwind part sign(A)(∂ x F h + S) neglecting the integrals for brevity, the right hand side of the computation considers substantially the flux difference between two nodes, which has the major problem that vanishes in correspondence of a steady discontinuity. So we have modified this term as:

|A| △ W -sign(A)c 2 △ h (40)
with sign(A) as defined in [START_REF] Dingemans | Water Wave Propagation Over Uneven Bottoms: Linear wave propagation[END_REF] and

|A| = R v 1 0 0 v 2 L ( 41
)

Positivity and Wet/dry Fronts

In case that the Riemann problem considers an initial solution with the position of the free surface at the same height of the beach, which means that the water depth is zero and there Inria is a wetting front which moves over a dry bathymetry, a correction has to be considered. The first problem considered is to avoid unphysical velocity values in regions close to the wet/dry interface. These values might arise from a zero by zero type division q/H. Following [START_REF] Ricchiuto | Stabilized Residual Distribution for shallow water simulations[END_REF], we introduce two small quantities, which represent limit values which have to be small compared to the mesh size :

ǫ u = ∆x 2 L and ǫ h = 10 -2 ǫ u (42)
The first quantity, is used to mass flux and velocity to zero (which, as it will be explained later, is the consistent limit value in presence of friction). So, in case H i results smaller than ǫ u , q i as well as u i are set to zero, as in [START_REF] Brufau | A numerical model for the flooding and drying of irregular domains[END_REF][START_REF] Ricchiuto | Stabilized Residual Distribution for shallow water simulations[END_REF]. The second quantity, ǫ h , is used to detect and flag nodes as being dry.

Next, we modify the scheme in order to make sure that the first order discretization is used in cells near the wet/dry front. As we will shortly show, this is enough to guarantee that H stays non-negative. Nevertheless, in practice, a control on nodes with negative water depth columns H i < 0 is performed in order to avoid non-physical behaviours. In particular, if H < 0 η is corrected with η i = -h i in order to result in an empty dry node H i = 0.

Exponential Function

In order to properly detect regions in proximity of dry areas, we use the exponential filter proposed in [START_REF] Ricchiuto | Stabilized Residual Distribution for shallow water simulations[END_REF]. The exponential considers in a stencil composed by five nodes i -2, i -1, i, i + 1, i + 2 over which we compute the minimum of the water column depth H min . We now define

y = H MAX -ǫ h max(ǫ h , (H min -ǫ h )) (43)
with H MAX computed as the max over the whole grid. We then define the filter :

f h = e -( 10 N y) 2 (44)
with N = L/h the number of cells. Each of the considered limiters is multiplied with this function, so that the first order scheme is recovered near dry areas.

Positivity preservation for the first order scheme

The first order upwind method we are considering is very similar to the Roe scheme described in [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamic -A practical introducion[END_REF] and used in [START_REF] Brufau | A numerical model for the flooding and drying of irregular domains[END_REF][START_REF] Castro | The numerical treatment of wet/dry fronts in shallow flows: Application to one-layer and two-layer systems[END_REF][START_REF] Delis | A robust high-resolution finite volume scheme for the simulation of long waves over complex domains[END_REF][START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF]. The main difference is given by the used averaging. We consider simple arithmetic averages with

u i+1/2 = q H = (q i+1 + q i )/2 (H i+1 + H i )/2
which is that used in the Q-scheme of van Leer. In Roe's approximation, used in [START_REF] Brufau | A numerical model for the flooding and drying of irregular domains[END_REF][START_REF] Castro | The numerical treatment of wet/dry fronts in shallow flows: Application to one-layer and two-layer systems[END_REF][START_REF] Delis | A robust high-resolution finite volume scheme for the simulation of long waves over complex domains[END_REF][START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF], they have

u Roe i+ 1 2 = H i+1 q i+1 + √ H i q i H i+1 + √ H i
We want to find the conditions under which our first order scheme preserves the positivity of H.

In the case of the scalar advection equation, we have seen that the scheme reduces to a classical upwind method, which is a LED one and preserves the monotonicity if integrated with a SSP method. In particular we have seen that it verifies a discrete maximum principle when integrated in time with Crank-Nicholson.

In the case of shallow water, we would like at least to make sure that given

H n i-1 ≥ 0, H n i ≥ 0 and H n i+1 ≥ 0, we will have H n+1 i ≥ 0.
To check if this is indeed the case, we consider the simpler case of explicit Euler time integration, and analyse the evolution equation of H :

H n+1 i = 1 - ∆t 2∆x (α i-1 2 i + α i+ 1 2 i ) H n i - ∆t 2∆x c i-1 2 i α i-1 H n i-1 + c i+ 1 2 i α i+1 H n i+1 + ∆t 2∆x c i-1 2 i ξ∆x u + c |u + c| - u -c |u -c| i-1 2 - ∆t 2∆x c i+ 1 2 i ξ∆x u + c |u + c| - u -c |u -c| i+ 1 2 (45)
For a flat bathymetry, that is for ξ = 0, the request that H n+1 i ≥ 0 is generally always satisfied provided that both α i-1 and α i+1 are non-negative, and under a CFL condition dictated by

α i-1 2 i +α i+ 1 2 i
. In appendix A ( A) it is shown that, under a constraint on the maximum acceleration in drying cells, with an explicit Euler integration method the scheme using the averages we have chosen permits to preserve H n+1 i ≥ 0 provided the time step is smaller than a limit temporal step ∆t positivity , which turns out to be always larger than the limit value ∆t stability = ∆x/ max i (|u i | + c i ) used to compute the time step :

∆t ≤ ∆t stability ≤ ∆t positivity

Well-balancedness and and wet/dry fronts

Having proved that the first order scheme is able to preserve the positivity, and bearing in mind the changes discussed in section in 3.3, it is necessary to verify that the terms due to the flux and source in the upwind are numerically balanced. The well-balancedness is the property to be able to consider in a correct way bathymetry variations, also in presence of dry areas [START_REF] Ricchiuto | Stabilized Residual Distribution for shallow water simulations[END_REF]. In particular, the scheme is said to have the C-property, if it preserves exactly the steady state lake at rest state corresponding to u = 0 and η = 0 [START_REF] Ricchiuto | On the C-property and generalized C-property of residual distribution for the shallow water equations[END_REF]. Even for schemes that do verify this property submerged regions, the presence of dry areas may generate spurious numerical waves [START_REF] Castro | The numerical treatment of wet/dry fronts in shallow flows: Application to one-layer and two-layer systems[END_REF][START_REF] Delis | A robust high-resolution finite volume scheme for the simulation of long waves over complex domains[END_REF][START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF][START_REF] Ricchiuto | Stabilized Residual Distribution for shallow water simulations[END_REF].

Two main elements have to be checked : the approximation of the source terms in the Galerkin component S i±1/2 in (28); the balance between bathymetry terms and gradients of the unknowns in the modified dissipation terms [START_REF] Melville | The surface velocity field in steep and breaking waves[END_REF]. The first term has been coded as:

S i+1/2 = - gH i+1/2 2 (h i+1 -h i ) = - gH i+1/2 2 ∆h i+1/2
and similarly for S i-1/2 . One easily checks that when summing these contributions to the centered flux difference part of the hydrostatic flux, the result is of course zero if the free surface is flat. Concerning the upwind stabilization part, we start by recalling the components of the extended form of (40), we evaluate the last expressions in correspondence of the lake at rest state u = 0 and η = 0. We obtain from the first as well as the second equation zero. The terms in the upwind part after the entropy fix modification are thus well-balanced in case of completely dry or wet nodes.

In presence of dry areas, this is not necessarily true anymore. Several conditions have to be maintained. The first is to ensure that ∆η is seen to be zero at the dry front [START_REF] Ricchiuto | Stabilized Residual Distribution for shallow water simulations[END_REF]. As the value of η in the first dry node only depends on the bathymetry, this is not necessarily the case, and an Inria artificial hydrostatic force might be seen by the scheme which can induce oscillations or spurious mass flows. In fact, referring to [START_REF] Castro | The numerical treatment of wet/dry fronts in shallow flows: Application to one-layer and two-layer systems[END_REF][START_REF] Delis | A robust high-resolution finite volume scheme for the simulation of long waves over complex domains[END_REF], the mass flow across a wet-dry front has to be zero and preserve the C-property of the scheme in all considered conditions which include the run-up on a non-constant slope. To satisfy these conditions, the bathymetry term referred to as ∆h has to be modified in correspondence of the presence of a wet/dry front. Following [START_REF] Castro | The numerical treatment of wet/dry fronts in shallow flows: Application to one-layer and two-layer systems[END_REF][START_REF] Delis | A robust high-resolution finite volume scheme for the simulation of long waves over complex domains[END_REF][START_REF] Ricchiuto | Stabilized Residual Distribution for shallow water simulations[END_REF], the source term due to the bathymetry changes is computed as

∆h i+1/2 =          h i+1 -h i if H i , H i+1 < ǫ h H i+1 -H i if H i+1 < ǫ h and h i+1 -h i ≤ H i H i+1 -H i if H i < ǫ h and h i+1 -h i ≥ H i+1
The last two conditions, correspond to detecting the presence of dry nodes "higher" than the last wet node (adverse slopes. In particular, in case

H i+1 < ǫ h if H i -h i ≤ -h i+1
, which means that if there is a positive slope which goes 'above' the reference of still water level and thus a possibly dry area, or in case of the node in i having the water depth below a certain level, and a descending slope with -h i >= H i+1 -h i+1 , ∆h is redefined as ∆h = H i+1 -H i . This allows to restore the condition η i = η i+1 .

Friction term discretization

When treating the water running over shores, and in proximity of dry areas, and in general, when the water depth is low, the friction term in the horizontal direction must be take into account. So, before proceeding to the final evaluation of the scheme, we discuss how this term has been included in the scheme. Introducing the Manning coefficient n:

C f = n 2 |u| H 4 3 (46)
The friction term in the shallow water equation results in:

f = gHC f u = gn 2 u|u| H 1 3 = gn 2 q|q| H 7 3 (47)
Note that this term tends to infinity as H → 0. In particular, the velocity ODE associated to this term (e.g. in the context of a splitting scheme) is

du dt = - gn 2 H 4/3
|u|u which shows that the consistent limit for H → 0 is u = 0. So, following [START_REF] Ricchiuto | Stabilized Residual Distribution for shallow water simulations[END_REF], this term will be computed only in case that the water depth column is higher than ǫ u , and consistently introduced both in the Galerkin and in the Upwind part, and, of course in the Jacobian of the Newton method described in section 4.2. Recasting the equations as

∂ t η + ∂ x q = 0 ∂ t q + f + ∂ x (uq) + g∂ x H 2 2 -gH∂ x h = 0 (48)

Numerical verification and validation for the Shallow Water equations

Wave run-up over a constant slope

To further verify and validate our implementation, we consider the tests of the run-up of a solitary wave with A=0.28m over a slope of 1:19.85 and initial depth of h 0 = 1m. The domain goes from -20 to 100 m, and the solitary wave positioned at half a wavelength from the beginning of the slope. A manning coefficient of 0.01 is used, and a CFL=0.3 and ∆x = 0.1. Thee parameters are the same used in [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF] (see also [START_REF] Delis | A robust high-resolution finite volume scheme for the simulation of long waves over complex domains[END_REF][START_REF] Synolakis | The runup of solitary waves[END_REF]).

Figure 4: Sketch of the test for the run-up of a solitary wave taken from [START_REF] Synolakis | The runup of solitary waves[END_REF] The experimental solution has been compared to the limiter given by the exponential function multiplied with each of the three considered limiters of Superbee, Smoothing Sensor and Monotonized Central-difference. X 0 = X a -1 [START_REF] Fuhrman | Simulation of nonlinear wave run-up with a high-order Boussinesq model[END_REF]. Considering the slope starting in X a = 24.85m, X b = 58m and X c = X d = 65m. The solitary wave is taken with an amplitude A = 0.28m using the generator of the solitary wave implemented in [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF] has been used. Comparing the three limiters individually with the first order scheme and with the exponential function limiter, far from the dry front this last gives back the third order linear scheme, which, as can be seen on figure 5, gives many oscillations when the initial solitary wave has converged into a shock. The first order scheme results to smooth too much the solution, while the limiters show an overall similar behavior. The crest that is identifiable on the wet/dry front in 5 is due to the presence of the friction, which gives a delay in the solution, and is normal. The smoothness sensor provides practically the same result as the high order scheme, with the exception of the shock which is captured in a perfectly monotone manner.

3 4 A h 0 acosh( √
The experimental data with the three implemented limiters each multiplied with the exponential function are compared further in 6. It emerges that in the first considered instant of the run-up the solution does not correctly reproduce the physical behavior of the wave (see 6(a)). This is principally due to the fact that the wave is in an area where the non-linear effects should not yet predominate, and the fact that the scheme considers only the shallow water equations leads to the formation of a shock which should not yet occur. This shows the need to introduce dispersive term, as, e.g., those displayed in the extended Boussinesq equations, allowing to better propagate smooth waves in near shores. The presence of the early developed shock (due to the adopted model) leads indeed to a shorter wave height when compared to the experimental data. In contrast to this initial difference, the agreement at later times is very good. This is probably a sign that, once the wave has broken, the Shallow Water equations give a good approximation of the wave dynamics. In particular the run-up and backwash phases are well reproduced, including the formation of a hydraulic jump during backwash. As discussed in [START_REF] Okamoto | The relative trough Froude number for initiation of wave breaking: theory experiments and numerical model confirmation[END_REF], wave breaking is a fundamental hydraulic phenomenon, especially in the nearshore region. Wave breaking drives the wave set-up and set-down (mean water level increase/decrease respectively), the wave run-up, the longshore currents, the rip currents, nearshore circulations,.. The wave force acting against coastal structures is also maximized at the wave breaking location, which is one of the most fundamental factors for practical applications in coastal engineering. Concerning the modelling of breaking waves, in the chapter devoted to the shallow water equations we have seen, through the example of the run-up of a solitary wave on a shore, that smooth waves converge into a shock too early when using this model. The Boussinesq equations, on the other hand, if not properly stabilized with some dissipative mechanism, cannot represent the dissipation of potential energy typical of wave breaking. This results in over-shoaling waves, becoming unstable due to the dispersive nature of the equations. Two issues have to be dealt with : how to detect breaking fronts and what is the modelling approach to be used.

As well explained in [START_REF] Okamoto | The relative trough Froude number for initiation of wave breaking: theory experiments and numerical model confirmation[END_REF], historically the first detection criteria were based on quantities computed using information over one full phase of the wave. The breaking formulations are described by water depth and wave characteristics, such as wave height, wave period, and wave length, at the wave breaking location. Wave characteristics are representative values over one wave cycle, i.e. phase-averaged. With the possibility of simulating shorter waves, thanks to the greater computational power available, new detection criteria have been proposed, based for example on local geometrical features of the free surface [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF][START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking and run-up. I: 1D[END_REF][START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF][START_REF] Sørensen | Boussinesq-type modelling using an unstructured finite element technique[END_REF]. These "phase-resolving" criteria are more stable, accurate, and, of course, computationally more convenient. They do however require that a sufficiently accurate description of the wave profile is available.

Once a breaking wave is detected, a proper mathematica has to be chosen. As already said, the Boussinesq equations give a good representation of the wave's propagation, shoaling, diffraction, etc, but break down when approaching too much the shore. In order to restore the correct physics, wave breaking has been incorporated into Boussinesq models by means of different artificial techniques, all introducing some dissipative mechanism allowing to reproduce the energy transformation process taking place during breaking. In this context, two main families of breaking models exist. Either a properly designed eddy viscosity term is added to the Boussinesq equations, or a hybrid formulation is used, in which one reverts to the shallow water equations.

In the first case, the momentum (or mass and momentum) dissipation allows to simulate the physical transformation of potential energy [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF][START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking and run-up. I: 1D[END_REF][START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF][START_REF] Sørensen | Boussinesq-type modelling using an unstructured finite element technique[END_REF], in the second case this is achieved via the energy dissipation through a non-linear shallow water shock [START_REF] Bonneton | Modelling of periodic wave transformation in the inner surf zone[END_REF][START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF].

The approach followed in this paper is the second one. The main reason for this choice is that, even theoretically, the Boussinesq equations degenerate into the non-linear shallow water equations as the dispersive terms become negligible in comparison to the non-linear ones.

Recalling the definition of the non-dimensional parameters ǫ = A/H (amplitude over depth), representing the non-linearity of the wave, and σ = H/λ (depth over wavelength), we can say that it makes perfectly sense to use the shallow water equations when ǫ = O(1) and σ → 0, as it is the case of non-linear (tall) waves close to the shore. In particular, this happens for example when hydraulic jumps occur at the transition from subcritical to supercritical flow. In order to classify these discontinuities, a Froude number, which is an analogous tool as the Mach number in Gasdynamics, is defined

F r = u √ gH (49)
When the transition is mild, which corresponds to a F r ≈ 1, on the water surface a train of smooth waves develops, and the Boussinesq model is able to give a good representation of the wave properties. Considering for example a bore between two steady states, we have a hydraulic jump or translating bore, and even if apparently different oscillation can be seen, it is actually only one jump. For such a simple configuration the mass can be written as H 1 u 1 = H 2 u 2 , having denoted with 1 the quantities at the toe of the wave's front and with 2 those of the wave's crest, and the same for the momentum, and the problem is controlled by a single number given by H 2 /H 1 , or equivalently by the Froude number [START_REF] Ricchiuto | Explicit Runge-Kutta Residual Distribution schemes for time dependent problems: Second order case[END_REF]. When the Froude number ranges between 1 and 1.3 the wave is not dissipative and no wave breaking occurs. In this case, it is possible to get a jump, but in fact it is a non dissipating jump, which means that there is no actual breaking. Above 1.3 the wave starts to break and for large Froude numbers a turbulent wave front develops. In fact, as the transition gets stronger, rollers start to develop and the jump becomes turbulent, and further, direct jumps are formed, which can be studied through the non-linear shallow water equations [START_REF] Tonelli | Simulation of wave breaking over complex bathymetries by a Boussinesq model[END_REF]. The transition between undular jumps and fully developed jumps can be compared to the change from non breaking and breaking waves and has been identified to happen in a range of the Froude number between 1.3 and 1.6 [START_REF] Okamoto | The relative trough Froude number for initiation of wave breaking: theory experiments and numerical model confirmation[END_REF][START_REF] Tonelli | Hybrid finite volume-finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Tonelli | Finite volume scheme for the solution of 2D extended Boussinesq equations in the surf zone[END_REF][START_REF] Tonelli | Simulation of wave breaking over complex bathymetries by a Boussinesq model[END_REF][START_REF] Tonelli | Shock-capturing Boussinesq model for irregular wave propagation[END_REF].

As remarked above, the information necessary to apply phase-resolving criteria requires a sufficiently accurate description of the waves. We pass from a non-dissipative, energy conserving set of equations to a hyperbolic non-linear system, which will immediately turn the smooth data into a discontinuity. This poses constraints which are different in the two limits.

Upwind Finite Element Discretization of the Hybrid Equations

Concerning the implementation of this approach, the idea is to start from the shallow water equations, and add the extra dispersive terms giving the enhanced Boussinesq system in regions in which the wave is not breaking, undergoing a hydraulic jump, or running up a shore (wet/dry fronts). We introduce a flag f M S assuming either the value 0 or 1, and allowing to switch on and off the terms related to the enhanced Boussinesq equations. This implementation is thus able to treat in a physically correct way the propagation and to change from the Boussinesq equations to the non-linear shallow water equations when non-linear terms become dominating. In this way, we locally switch to the hyperbolic non-linear model, whose shocks will provide the dissipative behavior sought. In the equations this gives:

∂ t q + ∂ x (uq) + gH∂ x η + gHC f u + f M S [-Bh 2 ∂ xxt q -βgh 3 ∂ xxx η+ - 1 3 h∂ x h∂ xt q -2βgh 2 ∂ x h∂ xx η] = 0 (50)
The discretization of these equations is done all along with the third order upwind stabilized SUPG finite element scheme, combined with the non-linear shock-capturing approach discussed in chapter 4, in the shallow water regions. The scheme resulting from all the modifications Inria introduced in this paper can be finally written as

∆x dWi dt + F (Wi+1) -F (Wi-1) 2 + ∆x 2 (S i-1 2 + S i+ 1 2 ) + ∆x fi + |A i-1 2 | 2 (Wi -Wi-1) - |A i+ 1 2 | 2 (Wi+1 -Wi) + sign(A i-1 2 ) 2 i i-1 S h - sign(A i+ 1 2 ) 2 i+1 i S h +δ i-1/2 { ∆x 6 [ dWi-1 dt - dWi dt + fi-1 -fi] + ∆x 2 sign(A i-1 2 ) dW i-1 2 dt + f i-1/2 } +δ i+1/2 { ∆x 6 [ dWi+1 dt - dWi dt + fi+1 -fi] - ∆x 2 sign(A i+ 1 2 ) dW i+ 1 2 dt + f i+1/2 } +fMS i (Ψ Ψ Ψ i-1/2 i + Ψ Ψ Ψ i+1/2 i ) = 0 (51)
having denoted by f the friction term array (0, f ) t (see section 4.6), and where the terms

Ψ Ψ Ψ i±1/2 i
represent the contributions of the left and right cell to the discretization of the dispersive terms in node i (see [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF] for further details). In particular, following [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF] these terms are composed of a Galerkin component in the second equation given by

i+1 i-1 ∂ x (Bh 2 ϕ i )∂ xt q h -ϕ i (βgh 3 ∂ x w η h - 1 3 h∂ x h∂ xt q h -2βgh 2 ∂ x h w η h ) dx
plus SUPG upwind stabilization terms reading, e.g.

in cell i + 1/2, sign(A i+1/2 ) 2      0 -Bh 2 i+1/2 dw q i+1 dt - dw q i dt - h i+1/2 ∂ x h i+1/2 3
dq i+1 dt - dq i dt -βgh 3 i+1/2 (w η i+1 -w η i ) -2βgh 2 i+1/2 ∂ x h i+1/2 w η i+1/2     
with h 2 i+1/2 cell average values computed using Gauss quadrature, and with w η and w q auxiliary variables obtained from

∆x w η i = - i+1 i-1 ∂ x η h ∂ x ϕ i dx, ∆x w q i = - i+1 i-1 q h ∂ x ϕ i dx
The interested reader can refer to [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] for more details. As a last remark, note that the shock detection is only activated in shallow water regions, so we have set

δ i = 1 if f M Si = 1.
The only remaining element to be determined is the criterion with which one model or the other is applied : a wave breaking criterion. In the following, three criteria will be presented, which will be referred to as the local, the hybrid and the physical criterion.

Wave Breaking Criteria

The most fundamental description of the wave breaking condition as a hydrodynamic phenomenon is that the wave becomes unstable when the particle velocity at the crest exceeds the celerity of the wave. This explains why wave breaking occurs but does not tell us where or when it would happen in the field. As recalled earlier, in literature two major families of wave breaking criteria can be been found: the phase-averaged and the phase-resolving. Phase-averaged wave breaking criteria result from early studies dating back even to the end of the 19th century [START_REF] Mccowan | On the highest wave of permanent type[END_REF][START_REF] Mitchell | On the highest waves in water[END_REF][START_REF] Stokes | On the theory of oscillatory waves[END_REF]. They are called phase-averaged because these formulas include wave characteristics, i.e. the wave height and length or period, and these wave characteristics are representative of one full phase of the wave. This implies some form of parametrization of the wave which is representative of one full period. The phase-averaged type can be furthermore divided into two categories. One type predicts the wave height at the breaking location, the other type describes the geometrical limit of the wave at the breaking location, as for example in the wave steepness criterion, or and the water depth limit criterion.

The second typology, the phase resolving models, use some information at certain locations of the wave for the parameters of the breaking position equation. Using detailed information within the phase of the wave, such as e.g. an angle or a vertical speed of the free surface elevation in a certain point [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF][START_REF] Okamoto | The relative trough Froude number for initiation of wave breaking: theory experiments and numerical model confirmation[END_REF][START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF], these models are expected to provide a more accurate breaking trigger, and for this reason they are the ones used in this work. Several criteria can be found in literature, and the interested reader can refer to [START_REF] Bonneton | Fully nonlinear weakly dispersive modelling of wave transformation, breaking and runup[END_REF][START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF] and references therein for an overview. The three criteria tested in this work are described in the following sections.

Local Wave Breaking Criteria

In [START_REF] Tonelli | Hybrid finite volume-finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Tonelli | Finite volume scheme for the solution of 2D extended Boussinesq equations in the surf zone[END_REF][START_REF] Tonelli | Simulation of wave breaking over complex bathymetries by a Boussinesq model[END_REF][START_REF] Tonelli | Shock-capturing Boussinesq model for irregular wave propagation[END_REF], Tonelli and Petti presented a criteria based on the parameter ǫ which considers the wave amplitude over the water depth as

ǫ = A h (52)
that grows with non-linearity. The ǫ has not been subject to any calibration [START_REF] Tonelli | Finite volume scheme for the solution of 2D extended Boussinesq equations in the surf zone[END_REF], but has been set to 0.8 according to the physical phenomenon of the limiting condition of the Froude number of 1.6 [START_REF] Okamoto | The relative trough Froude number for initiation of wave breaking: theory experiments and numerical model confirmation[END_REF][START_REF] Tonelli | Hybrid finite volume-finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Tonelli | Shock-capturing Boussinesq model for irregular wave propagation[END_REF], and considering that the approximate limit value in order to obtain a still stable solitary wave is ǫ = 0.78 [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF]. Clearly, this criterion acts locally on each node and has no memory of what has happened the instant before, nor in the neighbouring nodes. To avoid intermittency, and not to introduce instabilities in the scheme, once the ǫ assumes a value greater than 0.8, in correspondence of which he dispersive terms are switched off, in [START_REF] Tonelli | Simulation of wave breaking over complex bathymetries by a Boussinesq model[END_REF][START_REF] Tonelli | Shock-capturing Boussinesq model for irregular wave propagation[END_REF] in order to apply the Boussinesq equations again, ǫ has to be in a range between 0.55 ÷ 0.25.

Practical Implementation:

For each node in the spatial domain the formula in 52 can be seen as

ǫ = η h (53)
which implemented as

ǫ i = |η i | |h i | + ǫ h with ǫ h = 10 -2 ∆x 2 L ( 54
)
with the correction of ǫ h in order to not have incorrect solutions as the bathymetry approaches the reference level and gets zero. For each node i of the domain it is verified if ǫ i > ǫ c , with ǫ c normally put to 0.8, unless otherwise specified. The case of the water travelling in the positive direction, taken from left to right, and the negative are distinguished, as once the node with the exceeding quantity has been identified, the front has to be detected according to which direction the wave is travelling, and this is achieved simply by controlling if q i is ≥ 0 or ≤ 0. The nodes on the front of the breaking wave are flagged until the value of the nodes are below ǫ db = 0.3.

Inria

The obtained flagged nodes will correspond to a value of zero, which inserted in (50), will switch off the terms corresponding to the dispersive terms.

Finally, for this criteria, that will be referred to in this paper from now on as the local wave breaking criteria it is possible to say that the driving force is non-linearity, which, as discussed in section 5.2, is also related to the ratio between the difference of the water height at the crest and the toe of the front of the wave over the water depth at the wave's crest. Clearly, this criteria do not analyze the solution wave by wave, considering for each the elevation, but simplifies the problem to a generalized typology of wave. This however, will not necessarily work for waves travelling over arbitrary bathymetries. For example, in presence of wave set-up (or set-down), or in other words in presence of a change in the reference or average water level, this criteria will fail to recognize this change. A typical example is given by periodic wave run-up on a shore which very often presents an important set-up which will perturb the validity of this criteria. In conclusion, the local criteria is an excellent tool to test rapidly a model and give a first glance at the qualitative behavior occurring for a certain configuration, but is not at all reliable for detailed analysis.

Hybrid Wave Breaking Criteria

From the approach of the ǫ to detect when a wave becomes unstable due to the starting of a breaking, two further phase resolving type models have been identified in [START_REF] Okamoto | The relative trough Froude number for initiation of wave breaking: theory experiments and numerical model confirmation[END_REF]: the criterion of Schäffer [START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF] and of the eddy viscosity model of Kennedy [START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking and run-up. I: 1D[END_REF]. These two criteria are the building blocks the second criterion tested. In [START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF] the breaking is detected through the local slope angle of the front of the wave. In particular Schäffer introduced a formula with the tangent of the angle of the slope of the front being a function of time and two constant angles chosen to represent the angle at which the break starts, and one to determine the end. Similarly, Kennedy developed in [START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking and run-up. I: 1D[END_REF] an initial and final threshold value for the free surface elevation's vertical speed. This second criterion is also known as an eddy viscosity model, to whose category also belongs for example Roeber's criterion of detecting the wave breaking through a local momentum gradient ∂ x (uH) as detailed in [START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF]. From these two criteria, in [START_REF] Tomasicchio | The BCI criterion for the initiation of breaking process in boussinesq-type equations wave models[END_REF] a further criterion was developed: considering the non-linear transport equation for the free surface

∂ t η + c∂ x η = 0 (55)
the vertical speed of the free surface ∂ t η results equal to the negative product of the wave celerity and the free surface slope, which results to be the method of the characteristic for unsteady free surface flow problems. Moreover, in [START_REF] Tomasicchio | The BCI criterion for the initiation of breaking process in boussinesq-type equations wave models[END_REF] this has been rewritten as tan φ = ηt c obtaining through the thresholds that had been considered by Kennedy in [START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking and run-up. I: 1D[END_REF] and Schäffer in [START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF] a criterion with φ ∈ [8.5 • -33 • ]. A further approach based on [START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF] and [START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking and run-up. I: 1D[END_REF] has been analysed in [START_REF] Cienfuegos | Wave-breaking model for Boussinesq-type equations including roller effects in the mass conservation equation[END_REF], where starting from ∂ t η + c∂ x η = 0 and considering ∂ x η ≃ -tan φ the criterion of |∂ x η| ≥ tan φ was developed. Moreover, many other authors have used and considered the local steepness as a limiting tool to identify the breaking, as in [START_REF] Sørensen | Surf zone dynamics simulated by a Boussinesq type model. III. wave-induced horizontal nearshore circulations[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] for example. Kazolea in her PhD thesis [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF] considers two criteria which originate partially from Schäffer and from Kennedy. In fact, her criteria, which are those we will refer to as the hybrid wave breaking criteria is composed by two parts. The first is similar to Kennedy's and considers the vertical velocity component. If this exceeds a value which is proportional to the wave phase celerity and affected by the scale at which the wave is considered, the wave is considered as a breaking one and is given by

∂ t η > γ g|h| with 0.35 ≤ γ ≤ 0.65 (56)
where, as in [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF], γ is set to 0.6, unless otherwise specified. This criterion in inefficient in presence of steady hydraulic jumps for which ∂ t η ≈ 0. In order to be able to detect possible breaking fronts also for these occurrence, a further criterion is introduced, referring to the critical slope angle of Schäffer, and is that similar to that developed in [START_REF] Tomasicchio | The BCI criterion for the initiation of breaking process in boussinesq-type equations wave models[END_REF] with

|∂ x η| ≥ tan φ with φ ∈ [14 • , 33 • ] (57)
This criteria depends strongly on the breaking type, if spilling or plunging. Battjes in [START_REF] Battjes | Surf similarity[END_REF] describes the plunging as preceding the spilling, in accordance to the increasing wave steepness and decreasing slope angle and the parameter to which he refers to has a similar representation as [START_REF] Shi | A high-order adaptive timestepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF]. Following Kazolea in [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF], in the following paper, if not specified otherwise, the value has been set to φ = 30, similar to [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF]. Finally, having defined the two criteria which, in case of breaking one or both will be activated, the area where the non-linear shallow water equations have to be applied has to be identified. Following [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF] and [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF] and using the notation of figure 7, the magnitude of the roller, where a breaking has been detected and through which length the dispersive terms have been switched off, and in particular has been found experimentally to be with l r = 2.9(H 2 -H 1 ) and numerical tests showed that l N SW = 2.5l r .

Figure 7: Sketch of the lengths that determine the roller magnitude.

Lastly, in order to be able to identify when the wave stops to break, the Froude number

F r = 2 H2 H1 + 1 2 -1 8 (58)
with H 1 the water depth of the toe of the wave front and H 2 the crest's, is computed and if F r < F r c , the dispersive terms are switched on. The critical Froude is set to F r c = 1.3 (see [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF]).

Practical Implementation:

The implementation of this criteria has been subdivided in three phases, as it has been seen, that including phase 2 the flagged area was more reliable, because it eliminated the possibility Inria of having flagged a partial zone instead of the complete roller, due to the presence of a spurious oscillation or a saddle point:

PHASE 1: Initial Flagging Equation (56) turns into |η n+1 i -η i | ∆t ≥ γ g|h i | (59) while (57) becomes max(d x η) ≥ tan φ ( 60
)
with max(d x η)given by:

max η n+1 i-1 -η n+1 i+1 2∆x , η n+1 i -η n+1 i+1 ∆x , η n+1 i-1 -η n+1 i ∆x
If at least one of these criteria is verified, the node is flagged.

PHASE 2: Flagged region enhancement

In each package of possible breaking waves, identified through the previous flagging performed in phase 1, the maximum and minimum values of the free surface elevation is found. Moreover, taking the point in space corresponding to the maximum and minimum, the middle point is considered and there half of the length of lNSW = 2.5 2.9(H 2 -H 1 ) is computed. If the length obtained is larger than one stencil of our scheme, which is 4∆x corresponding to 5 points, the nodes within a distance 0.5 lNSW from the middle point are flagged.

PHASE 3: Computation of breaking length

As in phase 2, the positions of the maximum and minimum of the free surface, x maxj and x minj , in each wave package are found and the length l N SW is re-evaluated together with the Froude number [START_REF] Shiach | A temporally second-order accurate Godunov-type scheme for solving the extended Boussinesq equations[END_REF]. To make sure that breaking is applied on the right side of the wave, we compute the celerity

c b = q(x maxj ) -q(x minj ) ǫ u + η(x maxj ) -η(x minj ) (61)
with j denoting the wave package and with ǫ u = ∆x 2 L . Finally if, the de-breaking Froude number (58) is larger than 1.3, if l N SW is larger than 4∆x (one stencil), and if c b (x maxj -x minj ) < 0 (breaking on the front side of the wave), then the nodes within a distance 0.5 l N SW from the middle point are flagged and treated as shallow water nodes (f M Si = 0).

Physical Wave Breaking Criteria

As mentioned many times, a key role in wave breaking is played by the Froude number. Equivalent expressions of this quantity are used in the criteria already seen. In this work we also propose to test a physical definition of this quantity, based directly on its definition [START_REF] Ricchiuto | Explicit Runge-Kutta Residual Distribution schemes for time dependent problems: Second order case[END_REF] and referring to a local implementation of the convective criterion in [START_REF] Bjørkavåg | Wave breaking in boussinesq models for undular bores[END_REF]. As said more than once, wave breaking occurs when the free surface velocity exceeds the wave celerity. So we propose to study the quantity

F r = u s c b ( 62
)
with u s the free surface velocity at the crest of the wave, and c b the wave celerity. If F r > 1 then the wave should break. From now on this is referred to as the physical wave breaking criteria.

Using the fact that ∂ x u ≈ 0 at the wave's crest, and for constant bathymetry-slopes, as is roughly the case in the surf region, the expression for u s can be shown to be [START_REF] Bjørkavåg | Wave breaking in boussinesq models for undular bores[END_REF][START_REF] Cienfuegos | Numerical Modelling of Two Dimensional Water Wave Propagation Processes and Topographically Induced Breaking[END_REF]

u s = ū - 2 3 H 2 ∂ 2 ū ∂x 2 (63)
Concerning c b , an expression can be obtained by combining the expressions

∂ t η ≈ -c b ∂ x η and ∂ t η = -∂ x q, giving c b ≈ ∂ x q/∂ x η.
In practice, the celerity is computed as

c b = ∂ x q ∂ x η (64)
Practical Implementation: As in the hybrid criteria, we have three phases consisting in the initial detection of the breaking fronts, in the refinement of this initial guess, and in the final application of the breaking criterion. Phases 1 and 2 are performed exactly as in the hybrid criteria. In addition, even if not necessary, we compute u s over all the spatial domain (see 4.6). This is then used as a starting point and smoothed with one iteration of a 5 points moving average (smoothing over one Boussinesq stencil). Clearly in the future a more local procedure should be devised, allowing to compute u s only in the points necessary for the smoothing. For each identified front, the celerity ( 64) is approximated as

c b = ∆q ∆η (65)
where the ∆s are computed between wave crest and wave toe. Foe each front, taking the value of the smoothed u s in the point of the local maximum of η, we compute

F r = |u smax | |c b | + ǫ u ( 66
)
If Fr>1, and if the l N SW (computed as in the hybrid breaking criteria) is larger than 4∆x, and if u s (x maxj -x minj ) < 0 (breaking on the front side of the wave), than the nodes within a distance 0.5 l N SW from the middle point are flagged and treated as shallow water nodes (f M Si = 0).

Froude number approximation

Concerning the u s the approximation of (63) two different approximations have been compared: the finite differences and the finite element method. In the finite difference approach (63) reads

u s (x i ) ≈ ūi - 2 3 H 2 i u i+1 -2u i + u i-1 ∆x 2 (67)
while in the second case we have

us(xi) ≈ ūi + 2 3∆x i+1 i-1 ∂x(ϕiH 2 )∂x ū ≈ ūi + 2 3∆x 2 i i-1 ϕiH∂xh∂x ū + i i-1 H 2 ∂xϕi∂x ū +2 i+1 i ϕiH∂xh∂x ū + i+1 i H 2 ∂xϕi∂x ū ( 68
)
which is computed with quadrature formulas leading to an exact evaluation.

In order to see the differences and to choose a correct approximation let's consider a simple Inria example with a solitary wave propagating over a constant bathymetry and compare the analytical solution of the solitary wave, taken from [START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] with that of the finite differences and the finite elements with and without smoothing for both u s , and for c b evaluated as in [START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF]. Using the differential equations defining the the solitary wave of [START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF], we can derive an exact analytical expression for the free surface velocity

us a = q h 0 + q c a -2 H 2 3 c a h 0 + q c a - q (h 0 + q c a ) 2 d q c a ; ( 69
)
where c a is the analytical celerity [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF]

c a = gh 0 A 2 (A + 3h 0 ) 6h 2 0 (A -h 0 log( h0 + A h0 )) (70)
and where d q is the second derivative of q which is given by [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF]

d q = c a q - q 2 h 0 + q c a -g q 2 2c 2 a -gh 0 q c a Bh 2 0 c a -βg h 3 0 c a
The differences between the numerical approximations of the quantities used in the breaking criteria, and their analytical counterparts are shown in figure 9. In correspondence of the initial time iterations, during which the initial exact profile stabilizes onto a numerical one, the differences appear to be very high, and the approximations of u s oscillate quite a bit due to the presence of the second order derivative. After this initial transient, the computed c b stabilizes onto a constant value with a 0.36% error w.r.t the exact one, while the time history of the u s and Froude number are quite close tothe one obtained by applying the analytical formula to the solution. Comparing the Froude number given by the ratio of u s and c b for the finite element with and without smoothing, it is possible to see, that, the difference between the finite element and the analytical solution for the Froude number, have in general a difference of 0.53%, which is in fact quite satisfactory. This basically means that the physical Froude number to be considered in our computation should be rather 0.9947 than 1, but in the overall, the assumption is very accurate also for the unity Froude number. Note also that the finite element approximation fives much lower oscillations compared to the finite difference one. However, some smoothing is still necessary to minimize the intermittency of the breaking. The choice in smoothing the computed u s with the finite element approach is explained in figure 8. We present the profile of u s in correspondence of a very steep breaking front interacting with another wave. The figure shows the influence of performing more than one iterations of smoothing over 5 nodes. This results of course in a smoother curve, but that the peak are smeared too much, which would result eventually in an under-predicted Froude. The simple singular moving average over five nodes is the best compromise between smoothing and maintaining a reasonable height of the peak.

Wei's solitary waves shoaling on slopes test

In order to study the physics of the wave breaking, we consider the benchmarks proposed in Wei et al. [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF], involving solitary waves shoaling over different constant slopes. These tests are used to investigate the physical behavior obtained with the different breaking models, from the point of view of the shape of the wave and considering the physics of breaking. For this we will explicitly use the definition of breaking as the moment when the velocity of the peak of the wave is greater than that from the celerity on the front of the wave, and verify how close to this limit each criterion gets, as done in [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF].

The set-up of the first test is the following. A solitary wave of amplitude A = 0.2m is considered over a spatial domain of [-12, 40]m, with a bathymetry which starts at h 0 = 1m depth and a constant slope of 1:35 which begins in x = 0m, ends in x = 35m and continues with a zero slope until the domain's end. The solitary wave starts at the toe of the slope in x = 0m. The total considered time is 30s. We run the simulation with CFL= 0.3 on a mesh with ∆x = 0.05. Absorbing boundary conditions are set on the left and right ends of the domain using sponge layers 3m wide (see [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] for details).

The wave's evolution is described as a function of the dimensionless time :

t ′ = t g h 0
In particular, the shape of the wave has been studied at four different instances t ′ 1 = 16.24, t ′ 2 = 20.64, t ′ 3 = 24.03 and t ′ 4 = 25.94 in accordance to [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]. Note that t ′ 4 is the instant at which the start of breaking is observed in the reference.

As a first test, we compare the wave shapes obtained using two different limiters : the Superbee, and the Smoothness sensor. As the results of figure 10 show, the Superbee cuts the smooth peaks of the wave giving unphysical plateaux, while this is not the case with the Smoothness sensor which provides smooth peaks. This test has led us to choose, for all the following studies, the Smoothness sensor. The wave shapes obtained with the Smoothness sensor for different breaking criteria are reported in figure 11, where, in the top picture, we report the reference results, taken from [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF] and obtained with a fully non-linear non-hydrostatic model with dispersion and shoaling properties much closer to the theoretical ones than those of the Madsen and Sørensen equations used here.

From the pictures, we clearly see that the local breaking criterion results in waves which are too short and flat, which is a sign of early breaking. The hybrid criterion gives shapes much closer to the reference ones, even though still providing a wave profile at time t ′ 4 too flat, which is again the sign of an early breaking. When compared to the reference, the physical criterion gives the best results. To better understand these wave profiles, we have looked at the times at which wave breaking starts for all the models. Breaking with the local criterion occurs at t ′ = 23.02, with the hybrid at t ′ = 25.48, while with the physical criterion at t ′ = 28.04. In order to see, which possible Froude number for the physical criterion corresponds to the breaking time reported in [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF], we repeated the test Inria several times, resulting in the value Fr c =0.75 as it gives the wave breaking at t ′ = 25.91 which is practically the same as in the reference, and, as seen from figure 11, provides wave profiles identical to those obtained with Fr c =1. Clearly that the local criterion breaks too early and provides too much dissipation and too low wave heights, while the hybrid criterion results to be better, being also not too far away from the reference. The physical criterion with Fr c =1 gives an over-estimation of the breaking time w.r.t. the criterion used in [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]. The wave shapes at the first breaking time are compared for completeness in figure 12. The asymmetric shapes observed are those typical of shoaling waves. The result of the local criterion (top-right) is the one with the lowest height, while the other criteria show much higher peaks and a much larger wave skewness, closer to what is presented in [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]. The physical criterion with Fr c =1 is the one providing the highest and most asymmetric wave. For completeness, as plot overlaying the different wave profiles is reported on figure 13. The observations made are of course confirmed by this direct comparison. In order to better understand our results, following [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF] we have performed a comparison between temporal evolution of the wave celerity and of the particle velocity of the wave's crest. For the analysis, the horizontal particle velocity has been evaluated with the formula (68) with a 5 points smoothing (moving averages). The result has been made non-dimensional as u ′ s = u s / √ gh 0 . Next, we have computed the wave celerity c w . Under the assumption that η(x-c w (t)t), and q(x-c w (t)t). Note that this is not true, as in our study of propagation we have a wave which deforms, so c w is not only a function of time, but also of space. We will nevertheless assume that between two different instances the wave can be considered self-similar, or that at least this hypothesis can be used for the wave front. Proceeding as done in section 5, we end up with the approximate formula for c w = ∆q/∆η, which we have simply evaluated in this post-processing phase as the ratio q max /η max . We have introduced the non-dimensional celerity c ′ w = c w / √ gh 0 . The results obtained are summarized in figure 14. We recall that breaking is supposed to start in the moment when the free surface velocity is greater than the celerity. We see from the figures that the trends obtained are very similar to those reported in the reference [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF], perhaps with a slight difference in the values of the wave celerity computed in [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF] ĉ′ w = dxc dt / √ gh 0 , with x c the position of the wave crest. These results give a deeper view of the physics represented by each breaking model. The figures show that all the models give an increasing particle velocity at the crest, and a decreasing wave celerity at the front, which is what one should expect in shoaling conditions. The red circles in our results represent the start up of breaking (represented as a circle in the reference as well). The results of the local breaking criterion clearly show that breaking occurs too early when the particle velocity of the crest is about 50% of the wave celerity. Once more, the hybrid criterion seems to be better, predicting breaking when the particle velocity is slightly below 75% of the wave celerity. The physical criterion is of course the only one able to provide breaking exactly when the crest velocity equals the wave celerity, thus giving a physically correct description of wave breaking. For completeness, we report the plot relative to the physical criterion with Fr c =0.75 which is very similar to the one of the hybrid criterion, with breaking at a moment when the particle velocity is slightly above 75% of the wave celerity. This analysis shows that the physical criterion is the one based on the most sound physics. The hybrid criterion, however, is not so far away.

From the point of view of the physics, we can also remark the presence of a sudden increase in particle velocity u s , visible in all the plots, after breaking starts. In our case this is certainly a consequence of the wave profile sharpening, and eventually convergent into a shock, and of the asymptotic formula used to compute u s based on the curvature of the depth averaged velocity profile (see [START_REF] Bjørkavåg | Wave breaking in Boussinesq models for undular bores[END_REF][START_REF] Cienfuegos | Numerical Modelling of Two Dimensional Water Wave Propagation Processes and Topographically Induced Breaking[END_REF], and cf. equation ( 68)). However, experimental investigations of breaking waves have shown very similar velocity excursions in correspondence of breaking waves, as shown by the result reported in figure 15 taken from [START_REF] Melville | The surface velocity field in steep and breaking waves[END_REF] Figure 15: Comparison between the horizontal velocity (continuous line) and the free surface elevation (dashed line) of a broken and unbroken wave [START_REF] Melville | The surface velocity field in steep and breaking waves[END_REF] .

We consider now the same problem setup, but with a solitary of amplitude A = 0.3m shoaling on a slope 1:15. Following [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF] we start by looking at the dimensionless time instants t ′ 1 = 3.23, t ′ 2 = 6, t ′ 3 = 8.4 and t ′ 4 = 11.32. Differently from the previous case, these instants do not include a breaking point for the reference [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]. In our case all the criteria provide breaking before t ′ 4 . In particular, for the local criteria the breaking takes place at t ′ = 7.67, the hybrid predicts the breaking at t ′ = 9.32, the physical criterion predict breaking at t ′ = 10.99, and the physical criterion with Fr c =0.75 gives t ′ = 10.12 as the beginning of wave breaking. The comparison of the resulting wave profiles is reported in figure 16 together with the result taken by [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF]. Clearly, for this case both the local and the hybrid criteria break very early, giving at time t ′ 4 very low wave heights with a relatively sharp wave front already almost converged into a shock. Note that, for both profiles, the breaking region obtained does not contain the smooth maximum behind the sharp front of the wave. This gives to profiles different from the sawtooth wave one might expect, and gives the profiles shown in figure 16. The physical criterion, on the other hand, breaks much later, when the solitary has shoaled more and has already given a peaky and very asymmetric profile, closer to what is observed for breaking waves. As before, we compare the time evolution of the crest free surface velocity with the wave celerity. The results are summarized in figure 18. The results show that, as in the previous case, the local criterion predicts breaking at a stage where the free surface velocity is less than 50% (about 45%) of the wave celerity, while the hybrid criterion breaks at u s ≈ 0.65 c w , much earlier the in the previous configuration. As before, the physical criterion gives the a description which is physically the most correct. Note that, while both results obtained with the physical criterion show the same sudden increase in u s observed in the previous configuration, this is not seen for the local and hybrid criteria. Our explanation of this behavior is that, while the physical criterion gives wave profiles with very sharp fronts, the local and hybrid criteria develop profiles retaining a smooth extremum behind the bore for long time after the breaking. So when evaluating u s we miss the sharp front of the wave. In conclusion, it seems that the local breaking criterion consistently underpredicts the wave heights due to a very early breaking. The physical criterion seems definitely the best one, due to the good physical shape of the wave, and to a correct description of the physics of the breaking phenomenon. The hybrid criterion is not so far away, even though it still predicts breaking too early. The physical criterion with Fr c =0.75, chosen by matching the the breaking time of [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF] in the first test, is still better than the hybrid one.

Run-up of a Solitary Wave

We repeat the test of a solitary wave run-up on a sloping beach already seen in section 4.7. In particular, we consider again a solitary wave with amplitude A = 0.28m on a sloping beach with 1:19.85 to show the improvement obtained with the hybrid approach. The free surface values of η obtained have been considered at different instants, representing the wave shortly before the beginning of the shoaling, the shoaling and run-up and finally the back-wash, which are exactly the same times considered for the non-linear shallow water model. The results are reported on figure 19. The non-dimensional time t ′ is defined as t ′ = t g/h 0 .

During the propagation of the wave, slightly before the beginning of the shoaling, at t ′ 15 the shape of the wave corresponds perfectly in height and shape to the one given by the experimental data. The three criteria produce exactly the same solution, as none of them has predicted breaking yet. In comparison to the shape that had been obtained in 6 , the wave is higher and does not present a shock on the front. A slight phase advance seems to be present, which we assume begin related to the celerity of the solitary of [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] which is.

c = gh0 A h0 2 1 + A 3h 0 A h 0 -ln(1 + A h 0 ) > gh0 (71)
The difference between these two celerities is of a 10% for the considered test. The next considered instant is t ′ = 25 and displays the beginning of the run-up. Again the shape of the solutions are perfectly in agreement with the data. A very small phase advance is perhaps still visible. At t ′ = 30 the computed wave profile is till very close to the experimental data. In particular, it is possible to observe that during the run-up the water tongue has the same thickness observed in the experiments. At t ′ = 45 the solution still agrees quite well with the data, however a hydraulic jump is already visible in the numerical results, which has not developed yet in the experiments. At t ′ = 55 the experiments show an hydraulic jump during back-wash which is well reproduced by the simulations. The local criterion gives some weak oscillations, which might be due to a mis-capturing of the breaking across the jump. Not also that in this case the water tongue results thinner, despite the inclusion of the friction term, whose negligence had been alleged in [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II boundary conditions and validation[END_REF] to be the possible cause for this small disagreement with the data. Finally, at t ′ = 80 the numerical wave also has a proper behaviour and, in contrast to the other considered cases, seems even slightly delayed, which might be due to the friction terms that slows the solution close by the boundary of the bathymetry.

We find the overall behaviour of the solution in very good agreement with the data. The improvement w.r.t. the results shown in section 4.7 is impressive. The propagation of the wave is now physically correct, and we still capture nicely the run-up and back-wash of the wave. No great differences are observed between the three implemented breaking criteria, beside some occasional slight discrepancy, as for example the weak oscillations obtained with the local criterion in correspondence of the hydraulic jump. On the whole, they are almost equivalent for this test.

Periodic Wave over a Submerged Bar

We consider now the first of two benchmarks involving periodic waves, instead of a solitary wave. The breaking of such waves over a submerged bar is analyzed, referring to the set-up of Beji and Battjes laboratory experiment as reported in [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF]. This test, according to [START_REF] Shiach | A temporally second-order accurate Godunov-type scheme for solving the extended Boussinesq equations[END_REF], has been first done in Dingemans [START_REF] Dingemans | Water Wave Propagation Over Uneven Bottoms: Linear wave propagation[END_REF] to verify Delft Hydraulics numerical model HISWA and was then repeated by Beji and Battjes. The test consists in a periodic wave with amplitude and period A = 0.027m and T = 2.525s. The wave is generated in X 0 = 22m over a 0.4m initial high water level h 0 with a slope of 1:20 that begins at X a = 28m, continues until X b = 34m where a plateau of 0.1m depth starts for 2m in length and finally a downward slope of 1:10 until X d = 39m. For the simulations, the domain ranges from [0, 60]m with a ∆x = 0.04, a 3m long sponge layer and CF L = 0.3. The parameters of the local criteria where set, inspired by [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF], at ǫ = 0.6, while for the physical criteria the pre-flagging has been performed with γ = 0.3, Φ c = 15 and F r = 0.75, and the hybrid one with γ = 0.3, Φ c = 30 and the de-breaking at F r = 1.3. These parameters have been chosen following [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF]. In order to be able to compare the numerical data with the available experimental data, four gauges have been placed at the toe of the slope, at the end of the upward slope, in the middle of the plateau of the bar and at the end of it, which means in G 1 = 28m, G 2 = 34m, G 3 = 35m and G 4 = 36m. To generate periodic waves, we have used the internal wave generator proposed by [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF] to which we refer for details (see also [START_REF] Walkley | A numerical Model for Extended Boussinesq Shallow-Water Equations[END_REF]). As done in [23,[START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF], for the comparison with the experimental data the first four waves have been considered. This is done in order to avoid the effects of the reflection due to the boundary condition on the right. Indeed, besides reducing the wave heights, these reflections cause a phase error, which is especially visible in G 4 = 36m. For the comparison, the numerical data have been phase-calibrated only for the first gauge, and the same calibration has been applied to the other three. The temporal evolution of the wave height in the four gauges is reported in figure 23. The resulting comparison shows that, already for the second gauge, the shape of the wave obtained with the local and physical criteria result slightly higher than that of the experimental data and the hybrid criteria, even though the phase seems to be quite in agreement for each of the considered models. In G 3 in 23(c) the dashed line is higher than the other criteria and has an undulating behaviour right after the major peak. A similar behavior is observed in the last gauge, for which the local criterion over-predicts the main peak's height and provides a much more oscillatory signal. As for the physical and the hybrid criteria, a much better agreement with the experimental data is observed. Both criteria allow a better prediction of the waves' height reduction due to breaking. This agreement shows up also for the last gauge in correspondence of the beginning of the downward slope, which can be seen as the critical point of comparison, showing that a criterion has allowed to "dissipate enough of the waves energy". The hybrid criterion appears to be the best in agreement with the experimental data. The physical one provides a weaker and more intermittent breaking. To compare the breaking behavior during the time necessary for one crest to travel the length of the plateau, we report on figure 21 the wave profiles at the first breaking instance, at a intermediate time (same for all criteria), and at the last seen breaking instance for the considered wave. The hybrid criterion is the first to show wave breaking. Setting t = 0 at the first breaking instance for this criterion, the local criterion breaks relatively late at t = 0.88s while the physical criterion breaks at t = 0.47s. The central pictures, referring to t = 0.93s, shows that the local criterion has stopped breaking, while for both the hybrid and physical criteria we see the expected wave attenuation provided by wave breaking. The local criterion "flashes" again at t = 2.08s probably due to the decrease in bathymetry, but this sudden re-breaking has virtually no effect n the wave shape. On the contrary, the hybrid and the physical criteria keep on breaking until t = 1.96s and t = 1.79s, respectively. This gives the expected lowering in wave height and strength, as clearly shown fro the right-most top and bottom pictures in the figure .  Clearly, even with Fr c =0.75 the physical criterion gives a weaker breaking, which is shown by higher and more oscillatory wave profiles in both figures 23 and 21. For completeness, we report the evolution observed for Fr c =1 in the physical criterion, case in which no breaking is observed at all. The results are compared to those obtained with Fr c =0.75 in figure 22. In absence of breaking we obtain higher and more oscillating waves. In conclusion to this test, the hybrid criteria seems the best which approaches the physical behaviour of the wave's shape. However note that to obtain this result it is critical to lower the value of the constant γ to 0.3. Similar results can be obtained with the physical criterion if the value of Fr c is chosen sufficiently low.

Run-up of a Periodic Wave

Inspired by the test known as the spilling breaking type test of Hansen and Svendsen as depicted in [START_REF] Hansen | Regular waves in shoaling water: experimental data[END_REF][START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF], the shoaling and breaking over a shore of a set of regular waves given by a periodic source of a so-called spilling and spilling-plunging wave typology are considered. The spatial domain for the test spans from 0 to 65m, with 5m-long sponge layers. From an initial constant water depth of h 0 = 0.36m a 1:32.26 slope begins in X a = 34.775m. Periodic waves are generated with the internal wave generator of [START_REF] Filippini | Residual based discretization of Madsen and Sørensen system of Boussinesq equations for non-hydrostatic wave propagation[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF], which has been positioned at x = 20m, at a distance of 14.775m from the beginning of the slope, as in [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF]. As in the reference, the simulation is run with grid size of ∆x = 0.025, and with the CFL set to 0.35. Concerning the breaking criteria, we have set γ = 0.6 and Φ c = 30 • for both the physical and hybrid criteria, and Fr c =1 and Fr=1.3 respectively, while ǫ = 0.8 for the local one. Concerning the wave characteristics, we simulate the case 051041, with period and amplitude of T = 2s and A = 0.018m , corresponding to an Ursell number U = 4.8077, and the case 031041 with period and amplitude of T = 3.333333s and A = 0.0215m , corresponding to an Ursell number U = 17.5588. Figure 24 reports the wave profiles at the first and last breaking instants detected by each criterion for a single wave. Note that the shallow wave region is the one between the two vertical lines. The local criterion, at least with the implementation described in section 4.4.1, always breaks if h < 0 since in this case |η| = |H -h| > |h|. For all three criteria, breaking starts approximately at the same time, and in particular it is possible to see how the local criteria works in comparison to the other two, as the flagged length is much smaller. The first wave to stop to break is the hybrid one, while the last the local which, as already remarked, breaks in all the region in which h < 0. On figure 25 we compare the spatial distribution of wave height and mean water level (setup) with the experimental data of [START_REF] Hansen | Regular waves in shoaling water: experimental data[END_REF]. The results are qualitatively close to those reported in [START_REF] Kazolea | Mathematical and computational modeling for the generation and propagation of waves in marine and coastal environments[END_REF]. In terms of wave height, the results are not far from the data. However, we can see from the change in slope in the computed results that wave breaking is predicted too early by both the local and hybrid criterion. The physical criterion seems instead to predict a late breaking. Even so, the wave heights are under-predicted by all criteria, including the physical one. This, according to [START_REF] Veeramony | The flow in surf-zone waves[END_REF], might be due to the fact that the potential part of the model is based on the lowest order of the weakly non linear theory of Madsen and Sørensen found in [START_REF] Sørensen | Surf zone dynamics simulated by a Boussinesq type model. III. wave-induced horizontal nearshore circulations[END_REF]. In terms of set-up, we can clearly see the large over-prediction provided by the local criterion, probably due to the fact that for h < 0 breaking is activated independently on the wave profile. The results of the other two criteria are very close to the data, especially those of the physical criterion. For completeness Inria we also report the results obtained with the physical criterion for Fr c =0.75, which are somewhat similar to those obtained with the hybrid criterion. In particular, as also visible from the plots on figure 26 Fr c =0.75 yields an earlier breaking, as expected. The results of the test number 031041 are summarized in figures 28, 29, and 27, in terms respectively of wave profiles (with the shallow water region between the vertical lines), wave heights and mean water level distribution, and comparison between wave profiles obtained with the physical criterion for Fr c =0.75 and Fr c =1. The results are qualitatively similar to those of the previous configuration. The local criterion clearly provides an unphysical breaking in all the region where h < 0, while both the hybrid an physical criterion give more physically correct results. This is not only visible in the wave profiles of figure 28, but also in the set-up plots on figure 29 showing a large overestimation of the setup for the local criterion. The hybrid and physical criterion give instead, a very good estimation of the wave set-up. In terms of wave heights, the results of this case are somewhat better than those of the previous one. Al the criteria give now an early prediction of the wave breaking, but the underestimation of the wave amplitudes is much less pronounced than in the previous case. The results of the physical criterion are for this case clearly the ones closest to the experiments. Also for this case, we report for completeness the results of the physical criterion with Fr c =0.75. As before, we can see that the wave height and wave set-up plots (bottom row on figure 29 are somewhat similar to those of the hybrid criterion, with perhaps an even earlier prediction of breaking. This is also clearly visible from the plots of figure 27. In conclusion, for this test the physical criterion gives the best results, both in terms of wave height and set-up.

Inria 6 Conclusions

A non-hydrostatic hybrid model for the propagation, breaking and run-up of waves was presented. Starting from the adopted modelling approach, the limits given by the choice of using either the extended Boussinesq model equations or the shallow water ones were outlined, highlighting the dispersive nature of the enhanced Boussinesq equations, which is not capable of dispersing energy in concomitance of strongly non-linear phenomena, such as the propagation over an irregular bed topography, and the non-linear character of the Shallow Water equations that dissipates energy regardless of the bathymetry variations. A hybrid approach was proposed to combine the best of both models. Considering the model of the frictionless Shallow Water equations, section 3 presented first the adopted discretisation in space with the upwind stabilized finite element method and the Crank-Nicholson method for the time, showing the structure of the adopted scheme. Shock capturing techniques, in order to avoid numerical oscillations, were included for the advective equation in presence of a discontinuity and were introduced in the non-linear Shallow Water scheme through mass-lumping. The change from the first-order scheme back to the third-order is regulated by a suitably designed limiter defined on each node of the domain. In particular, among the three implemented limiters, the Monotonized Central-difference, the Superbee and the self devised smoothing sensor, which have all in common the evaluation of the slope of the free surface elevation, the smoothing sensor was found to be the best compromise between smoothing and the capability to maintain an accurate representation of the motion of the waves, observation driven by considering a Riemann problem and comparing the analytical solution with the numerical one. Further, investigating a second typology of Riemann problem, the scheme was shown to be entropy violating, and therefore an entropy fix of the Harten and Hyman's family was introduced to correct the value of the eigenvalues of the Jacobian matrix in case of very small values. Moreover, the dissipation of the scheme was rewritten in a non conservative form, in order to prevent, in case of stationary expansions, the neglecting of the associated terms and thus allowing a correct dissipation of the discontinuity. Moreover, being the scheme not able to treat wet/dry fronts, an exponential function has been developed and considered in addition to the implemented flux-limiters. The consideration of possible dry fronts required an inspection, developed in Appendix A, of the conditions guaranteeing the preservation of the positivity of the solution. This property has resulted to be always verified with the chosen time step depending on the CFL number . The scheme was also checked for the well-balancedness of the flux and source terms, and the resulting unbalance led to a re-implementation of the terms related to the bed-topographies changes, in particular enabling the scheme of correct treating of the transition of water over dry bathymetries. Finally, after having introduced in the scheme the friction terms, the numerical results of a Riemann problem, as well as the comparison with experimental data of the run-up of two different solitary waves on a shore were presented. In particular, from a considered test with the implemented scheme outlined the bad approximation for the shoaling, motivating the introduction of the hybrid system of equations in section 4. After an introductory part of the existing different wave breaking approaches the coupling of the non-linear Shallow Water equations with the dispersive terms, which together form the extended Boussinesq equations of Madsen and Sørensen were presented. The description of three considered breaking criteria followed. From the study of two existing wave breaking criteria, the third one was proposed by means of the physics, which considers the Froude number. Through different shoaling tests in section 5, the physics and behaviour of the three implemented criteria in the hybrid scheme has been studied, focusing also on the analysis of the velocity of the free surface on the crest and the celerity, which are the two quantities characterising the physical wave breaking criterion, showing the accuracy of all three criteria with the available experimental data, and in particular highlighting the good physical approach of the physical wave breaking criterion. Moreover, these observations have been confirmed through an extensively testing on different benchmarks with excellent agreement with the experiments. This new hybrid equation approach with the use of an upwind stabilized finite element method and wave breaking criteria have been shown in the overall to represent an accurate and physically reliable model. Finally, the positive results of the final tests confirmed the solidity of the physical wave breaking criterion, to be able to detect very good breaking fronts, which is seen as a promising preamble for future works.

Future works

In future works the obtained hybrid equation model will be extended from one to two spatial dimensions. The representation of the shock capturing scheme will be a demanding quest, as in contrast to the monodimensional case, the two dimensional problem will require beside a mass lumping also an additional term which might have the appearance of a dissipative term. Moreover, the way to implement the breaking criteria will have to be studied. The first two criteria are not expected to rise extra difficulties, as they have already been implemented in the two dimensional case, although in a context of a different scheme model. The physical criterion is foreseen to be the most riddling one, as for example the choice of node for the computation of the celerity is not intuitive, since in the two-dimensional case often it is not possible to speak in a proper way of a real moving front. Beside the coupling of the non-linear Shallow Water equations and the Madsen and Sørensen extended Boussinesq ones, further models are going to be implemented, such as the Serre Green-Naghdi one, which, beside the good dispersive terms, includes terms that are able to represent very well the non-linearity, although resulting slightly more complicated than our considered scheme. Also the time integration will be improved reaching at least the third order. Similarly, higher order finite elements could be used for the spatial numerical approximation, as for example with Lagrange P 2 and P 3 , hermite finite elements or other high order polynomials. Finally, the scheme could apply adaptive un-structured meshes, in order to represent with a better resolution the motion of the wave front and of the water line during the wave's propagation, breaking and run-up over shores.

Inria

A Positivity of First-order Scheme

In the following, the conditions mentioned in 3 for which if H n i-1 ≥ 0, H n i ≥ 0 and H n i+1 ≥ 0 we obtain H n+1 i ≥ 0, which means that the scheme of the first order preserves the positivity of the depth, are here presented.

∆x(H n+1 i -H n i ) + T i-1/2 + T i+1/2 = 0 (72)
where:

T i-1 2 = ∆t 2 {(q i -q i-1 )+ + 1 2c (|u + c| -|u -c|)(q i -q i-1 )+ + 1 2c (H i -H i-1 )(|u -c|(u + c) -|u + c|(u -c))+ - c 2 (h i -h i-1 )( u + c |u + c| - u -c |u -c| )} (73)
and

T i+ 1 2 = ∆t 2 {(q i+1 -q i )+ - 1 2c (|u + c| -|u -c|)(q i+1 -q i )+ + 1 2c (H i+1 -H i )(|u -c|(u + c) -|u + c|(u -c))+ + c 2 (h i+1 -h i )( u + c |u + c| - u -c |u -c| )} (74)
with all quantities without sub-script evaluated using local averages. The proof will consist in two parts: in the first one it will be assumed to have a constant bathymetry, so that the last term in ( 73) and ( 74) disappears, and once the conditions for this case have been found, the conditions for a changing bathymetry are considered. Moreover, the cases of wetting and drying cells will separately considered and for each case the conditions in case of a sub-critical or supercritical flux will be analysed.

Hypothesis: constant Bathymetry

We start by rewriting [START_REF] Veeramony | The flow in surf-zone waves[END_REF] as:

H n+1 i = 1- ∆t 2∆x (c i-1 2 i α i-1 2 i + c i+ 1 2 i α i+ 1 2 i ) H n i + ∆t 2∆x c i-1 2 i α i-1 2 i-1 H n i-1 + c i+ 1 2 i α i+ 1 2 i+1 H n i+1 (75)
having set

T i±1/2 = ∆t 2 c i±1/2 i α i±1/2 i H n i -α i±1/2 i±1 H n i±1
Provided that α i-1 and or α i+1 are non-negative, we can always satisfy the condition that H n+1 i ≥ 0, provided that the time step is small enough. In case either α i-1 , or α i+1 is negative, this is not true anymore, and we may have H n+1 i < 0 whatever the time step. In the following these conditions are studied.

One dry node

Considering now the case for one dry nodes H n i+1 , q n i+1 = 0 a, we will look for the conditions that verify

H n+1 i ≥ 0. 2T i+ 1 2 ∆t = {-q i + 1 2 q i (|1 + F | -|1 -F |)+ - c i+ 1 2 2 H i (|1 -F |(1 + F ) -|1 + F i+ 1 2 |(1 -F ))} =c i+ 1 2 α i i+ 1 2 H i (81)
with, similarly to (77) and (78), F = u i /c and F = u/c. In this case, we get again that F = F , so that:

α i+1/2 i = -F + F 2 (|1 + F | -|1 -F |)+ - 1 2 (|1 -F |(1 + F ) + |1 + F |(1 -F )) (82)
Which, in our hypothesis that the node belonging to i + 1 is dry, has only the terms relative to node i. Concerning the term T 1 , this time two terms, i and i -1, are both non zero, and from (72):

2T i-1 2 ∆t = c i-1 2 (α i-1 2 i H i -α i-1 2 i-1 H i-1 ) (83)
with (using the same notation as before):

α i-1 2 i-1 = F + F 2 (|1 + F | -|1 -F |)+ + 1 2 (|1 -F |(1 + F ) + |1 + F |(1 -F )) (84)
RR n°8536

for ξ ≤ 0 (back wash case) there is no problem since this contribution is always increasing. If ξ > 0 (uphill supercritical drying) the scheme becomes (cf. previous analysis)

H n+1 i = ∆t ∆x c i-1/2 i F H n i-1 -∆t|ξ|c i-1/2 = ∆t ∆x c i-1/2 F H n i-1 -|ξ|∆x
which leads again to a restriction on the spatial resolution allowing to capture this (unphysical) uphill supercritical drying with positive depths :

∆x ≤ F H i-1 |ξ|

Figure 3 :

 3 Figure 3: Comparison between different choices of limiters: (a) analytical solution, first-order and third-order scheme; (b) Analytical solution, Smoothness sensor, Superbee and Monotonized Central-difference limiter; (c), (d) zooms of (b).

Figure 5 :Figure 6 :

 56 Figure 5: Comparison for the runup of the solitary wave with A=0.28m between the first order scheme and the exponential function with the limiters: Superbee (top), Smoothness Sensor (center) and Monotonized Central Differences (bottom) at two different instants t g h0 = 15 and t g h0 = 25.

Inria 4 Hybrid

 4 Equations for Wave Breaking Treatment 4.1 Hybrid Shallow Water -Boussinesq coupling for wave breaking

Figure 8 :

 8 Figure 8: Comparison of different approximations for the u s computation in a breaking area.

Figure 9 :Inria 5

 95 Figure 9: Comparison between different approaches with the analytical solution for (a) the velocity of the free surface corresponding to the wave's crest (b) celerity of the wave (c) Froude number.

Figure 10 :

 10 Figure 10: Comparison of wave shapes at four different times computed with the Superbee limiter (top) and the Smoothness sensor (bottom), using the physical criterion

Figure 11 :

 11 Figure 11: Parametrical study for the physical breaking criterion: different Froude numbers for A=0.2m and 1:35 slope -(a) Reference results, (b) Local criterion, (c) Hybrid criterion, (d) Physical Criterion, (e) Physical criterion with Fr c =0.75.

Figure 12 :

 12 Figure 12: Wave breaking instants for A=0.2m and slope 1:35 -(a) Local (left) and hybrid (right) wave breaking criteria; (b) Physical wave breaking criteria with respectively Fr=0.75 (left), Fr=1 (right).

Figure 13 :

 13 Figure 13: Direct comparison of the wave profiles. (a) Wei's results for A=0.2 and slope 1:35, (b) local, hybrid and physical criteria superposed, (c) local, hybrid and physical with Fr c =0.75 superposed.

Inria

Figure 14 :

 14 Figure 14: Parametric study of the breaking criteria: time evolution of wave celerity and crest free surface particle velocity for A=0.2m and 1:35 slope-(a) Reference results, (b) local criterion. (c) hybrid criterion, (d) physical criterion, (e) physical criterion with Fr c =0.75.

Inria

Figure 16 :

 16 Figure 16: Direct comparison of the wave profiles. (a) Wei's results for A=0.3 and slope 1:15, (b) local, hybrid and physical criteria superposed, (c) local, hybrid and physical with Fr c =0.75 superposed.

Figure 17 :

 17 Figure 17: Wave breaking instants for A=0.3m and slope 1:15 -Local (left), hybrid (middle) and physical with Fr=1 (right).

Figure 18 :

 18 Figure 18: Parametrical study for the physical breaking criterion: time trend of the celerity and velocity profiles in the wave's crest for A=0.3m and 1:15 slope-(a) Reference results (e) Fr=0.75, (d) Fr=1.

Figure 19 :

 19 Figure 19: Run-up of a solitary wave with A=0.28m. Results obtained with local, hybrid, and physical breaking criteria. (a) t ′ = 15 (left) and t ′ = 25 (right), (b) t ′ = 30 (left) and t ′ = 45 (right), (c) t ′ = 55 (left) and t ′ = 80 (right).

Inria

Figure 20 :

 20 Figure 20: Sketch of the submerged bar test case.

Figure 21 :

 21 Figure 21: Comparison of the instances when the breaking first occurs on a wave (left) at a midtime (center, same for all criteria) and the last breaking point (right), using in (a) the hybrid, (b) the local, and (c) the physical criterion.

Figure 22 :

 22 Figure 22: Comparison of the instances when the breaking first occurs on a wave (left) at a midtime (center) and the last breaking point (right) for the physical criteria with (a) Fr=0.75 and (b) Fr=1.

Inria

Figure 23 :

 23 Figure 23: Time evolution of the free surface in the four gauges. (a) Gauge G 1 (beginning of the upward slope), (b) Gauge G 2 (beginning of the plateau), (b) Gauge G 3 (midpoint of the plateau), (d) Gauge G 4 (beginning of the downward slope).

Figure 24 :

 24 Figure 24: Hansen and Svendsen test 051041. Snapshots of the first and last instant of breaking of the (a) local, (b) hybrid and (c) physical (Fr=1) breaking criteria.

Figure 25 :

 25 Figure 25: Hansen and Svendsen test 051041. Wave height (left) and mean water level (set-up) (right) for (a), (b) hybrid, (c) physical criteria, and for the physical criterion with Fr c =0.75 (d).

Figure 26 :

 26 Figure 26: Hansen and Svendsen test 051041. Snapshots of the first and last instant of breaking for the physical breaking with (a) Fr=0.75 and (b) Fr=1.

Inria

Figure 29 :

 29 Figure 29: Hansen and Svendsen test 031041. Wave height (left) and mean water level (setup) (right) for (a), (b) hybrid, (c) physical criteria, and for the physical criterion with Fr c =0.75 (d).

•

  |1 + F | = 1 + F • |1 -F | = -1 + FAnd (80) results 2F > 0 which is always fulfilled.Drying of the cell: F < 0 Sub-critical Case: F > -1 (1+F > 0) This implies that• |1 + F | = 1 + F • |1 -F | = 1 -FAnd (80) results 1 + F > 0 which is always fulfilled due to the sub-critical condition.Super-critical Case: F < -1 (1+F < 0) This implies that• |1 + F | = -1 -F • |1 -F | = 1 -FAnd (80) results zero which means that this is always fulfilled.

  Ricchiuto and Filippini, J.Comput.Phys. 2014) au schéma de Roe classique. Le couplage entre le modèle Boussinesq et Shallow Water est en suite étudié. En particulier, on considère différent critére physiques de détection de fronts déferlants. En particulier, on présente une implementation numérique locale du critère convectif de (Bjørkavåg and H. Kalisch, Phys.Letters A 2011), qui est comparée au critères proposés dans (Kazolea et. al, J.Comput.Phys. 2014) et (Tonelli and Petti, J.Hydr.Res. 2011). Le modèle obtenu est validé sur des nombreux benchmarks avec données expérimentales.

RR n°8536

Bacigaluppi & Ricchiuto & Bonneton

Inria

351, Cours de la Libération Bâtiment A 29

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

One dry cell

Considering now the case for two dry nodes H n i+1 , q n i+1 = 0 and H n i , q n i = 0, we will look for the conditions that verify H n+1 i ≥ 0. First of all T 2 = 0 due to the dry cell, now studying the terms i -1 of (73):

defining the Foude number as:

We recall that the local average of the velocity is obtained as

which in this case leads to u = u i-1 resulting in F = F and moreover

We can thus rewriting (78) as :

In order to verify the positivity of H n+1 i , we will have to find the values which guarantee that α i-1 ≥ 0.

which results from (79).

Wetting of the cell: F > 0 Sub-critical Case: F < 1 (1-F > 0) This implies that

And (80) results 1 + F > 0 which is always fulfilled.

Super-critical Case: F > 1 (1-F < 0) This implies that Inria so that the values that fulfil the request have to give α i-

Wetting of the cell:

Drying of the cell: F < 0

Sub-critical Case: F > -1 (1+F > 0) Doing the same considerations as for the two dry nodes, it is obtained that:

which is always positive if if the flow is sub critical in both nodes, i.e. 1 + F ≥ 0, and is satisfied in presence of a smoooth trans-critical acceleration, provided that

which will always satisfied if the mesh size is fine enough.

Super-critical Case: F < -1 (1+F < 0)

i-1 = 0, which means that it is always fine.

Note: The analysis of the positivity of α

and is omitted for brevity.

Concerning the coefficients of node i, this is the most general case in which we have the contribution from both the terms T i+1/2 and T i-1/2 :

Inria with F = u i /c. These terms will lead to a time step limitation associated to

This condition is analyzed in the following. Note that the time step used in the computations verifies ∆t ∆x max i

Wetting of the cell: F > 0 Sub-critical Case: F < 1 (1-F > 0)

Under our assumptions we have F < 1, u < c, that lead to F u < u → c -F u > c -u > 0, leading to:

This two considerations bring to the conclusion that α i < 2u i + 2c i , which leads to the constraint

and is in practice always verified. Super-critical Case: F > 1 (1-F < 0)

which is thus also always verified.

Drying of the cell: F < 0 Sub-critical Case: F > -1 (1+F > 0) In this case we can again write

The analysis is identical to the sub-critical wetting case and is omitted for brevity.

Super-critical Case

and thus always verified.

Hypothesis: ∆h = 0

We now the typical case of a bathymetry with a constant slope ξ so that h i = h i-1 + ξ∆x. Only the case of a limit wet dry cell is considered The first order scheme becomes:

are the term of the equation that have to be studied in order to obtain the restrictions that have to be verified to preserve the positivity of the depth.

Sub-critical Case:

Wetting of the cell with ξ > 0 and two dry nodes (H i = H i+1 = 0): (87) becomes, under the consideration that

which can be considered fulfilled for any slope, negative and positive. Omitting for brevity the other cases for drying and for only one dry node, it is possible to say that in the sub-critical case the non-negativity of the term related to the bathymetry is always fulfilled.

Super-critical Case:

Wetting of the cell with two dry nodes: Considering similarly to the previous case (87):

for ξ ≥ 0 (downhill flooding case) there is no problem since this contribution is always increasing. If ξ < 0 (run up) the scheme becomes (cf. previous analysis)

which leads to a restriction on the spatial resolution allowing to capture supercritical run up with positive depths :

Drying of the cell with two dry nodes:

In case of drying cells, which means u i < 0, the term becomes