Patrick Marsaleix
email: patrickmarsaleix@aero.obs-mip.fr

Héloïse Michaud

Claude Estournel

3D phase-resolved wave modelling with a non-hydrostatic ocean circulation model

Keywords: phase-resolved wave modelling, wave breaking, rip currents, non-hydrostatic, ocean circulation model

A phase-resolved wave model is derived from an ocean circulation model for the purpose of studying wave-current effects in nearshore zones. One challenge is to adapt the circulation model to the specificities of wave physics. This mainly concerns the consideration of nonhydrostatic effects and the parametrization of wave breaking. The non-hydrostatic pressure is calculated using the artificial compressibility method (ACM). The ACM-induced errors on wave dispersion properties are examined in detail in the context of the linear theory using idealized test cases. The possible compromise between the precision achieved on nonhydrostatic physics and the adjustable CPU cost of the ACM method is looked at in detail.

The modification of the wave characteristics by the bathymetric slope and the breaking of waves are then examined from a linear slope beach laboratory experiment. Finally the model is evaluated on the issue of rip currents and their feedback on the wave field using a laboratory experiment of a beach with a bar intersected by channels.

comparable to

), insofar as it is recommended to use the ACM method with $ such that # !

Introduction

The interactions between waves and current play a major role in coastal modelling studies and are important to correctly reproduce the nearshore circulation. Coastal engineering applications need realistic descriptions of waves and currents, their transformations towards the nearshore and behaviours with artificial coastal structures for the design of harbours or dykes to prevent overtopping and to reduce submersion or erosion risks and this is only possible with a circulation model that reproduces explicitly the waves. Taking into account the waves in the circulation models allows on the one hand the representation of the residual flows generated by the waves such as the Stokes drift [START_REF] Mellor | The Three-Dimensional Current and Surface Wave Equations[END_REF] or the rip currents [START_REF] Castelle | Rip current types, circulation and hazard[END_REF] and on the other hand, to take into account turbulence generation by waves in the turbulent closure scheme, which improves the oceanic surface and bottom boundary layers, the vertical current shear, and the vertical temperature and salinity profiles in the circulation model [START_REF] Uchiyama | Wave-current interaction in an oceanic circulation model with a vortex-force formalism: Application to the surf zone[END_REF]. In shallow waters, when the waves effects near the bottom are no longer negligible, the bottom stress is significantly increased, allowing the resuspension of sand and sediments [START_REF] Soulsby | Bed shear stresses due to combined waves and currents[END_REF], the Stokes drift, the undertow or rip currents then playing an important role in their transport in the water column.

The wave-current interaction is also decisive for the mean surface level set-up during storms, the modification of the surface stress related to the roughness created by the waves, the modelling of the marine submersion of the coastal continental surfaces [START_REF] Bertin | Importance of wave age and resonance in storm surges, The case Xynthia, Bay of Biscay[END_REF], the modification of the coastline, the evolution of the shape of the seabed [START_REF] Bouharguane | Low complexity shape optimization & a posteriori high fidelity validation[END_REF] .

At the spatial scales of the coastal circulation models (see review by [START_REF] Klingbeil | The numerics of hydrostatic structured-grid coastal ocean models: state of the art and future perspectives[END_REF], such studies are mainly done using a hydrostatic circulation model coupled with a spectral phase-averaged wave model [START_REF] Mcwilliams | An asymptotic theory for the interaction of waves and currents in coastal waters[END_REF], Ardhuin et al 2008[START_REF] Bennis | On the coupling of wave and three-dimensional circulation models: Choice of theoretical framework, practical implementation and adiabatic tests[END_REF][START_REF] Benetazzo | Wave-current interaction: Effect on the wave field in a semi-enclosed basin[END_REF][START_REF] Kumar | Midshelf to Surfzone Coupled ROMS-SWAN Model Data Comparison of Waves, Currents, and Temperature: Diagnosis of Subtidal Forcings and Response[END_REF][START_REF] Uchiyama | Three-dimensional transient rip currents: Bathymetric excitation of low-frequency intrinsic variability[END_REF]. The phaseaveraged wave model feeds the circulation model with Stokes drift, residual wave pressure, momentum and turbulence production terms related to wave erosion. In return, the circulation model gives the wave model the variations of sea surface level and surface current, the latter being able to significantly modify the wave field [START_REF] +appendices Uchiyama | Wave current interaction in nearshore shear instability analyzed with a vortex force formalism[END_REF]) as in the known example of rips current. The choice of a phase-averaged wave model is largely guided by computational time considerations as it provides mean properties of the sea state.

The spectral model can indeed be used with a horizontal resolution much lower than the wave wavelengths and therefore cover large areas with a reasonable number of grid points.

An alternative is the phase-resolved model, or deterministic model, that explicitly describes the evolution of the free surface. The SWASH models [START_REF] Zijlema | SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters[END_REF], REF/DIF [START_REF] Kirby | A parabolic equation for the combined refractiondiffraction of Stokes waves by mildly-varying topography[END_REF], FUNWAVE [START_REF] Kirby | FUNWAVE 1.0; Fully Nonlinear Boussinesq Wave Model Documentation and User's Manual[END_REF], BOSZ [START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF], NHWAVE [START_REF] Ma | Shock-capturing non-hydrostatic model for fully dispersive surface wave processes[END_REF] are widespread examples of this type of model, very often mentioned in the wave modelling literature. The processes resulting from the effects of bathymetry such as refraction, diffraction and reflections [START_REF] Magne | Topographical Scattering of Waves: Spectral Approach[END_REF], interaction with other waves or with the background circulation, infragravity waves (Bertin et al, 2018), the asymmetrical wave transformation at the approach of the breaking phase, are, in principle, better represented by this type of model [START_REF] Rusu | Modelling Waves in Open Coastal Areas and Harbors with Phase-Resolving and Phase-Averaged Models[END_REF]. The phaseresolved model is also well suited for sand mobilization, which is directly dependent on orbital velocity near the bottom [START_REF] Bouharguane | Low complexity shape optimization & a posteriori high fidelity validation[END_REF]. The phase-resolved model has a continuous approach of frequencies and directions, while the phase-averaged model solves its equations on a discrete spectrum of frequencies and directions which results in a relative inaccuracy related to frequency and direction resolution. For example, in [START_REF] Michaud | Three-dimensional modelling of wave-induced current from the surf zone to the inner shelf[END_REF] the Wavewatch III r model [START_REF] Tolman | User manual and system documentation of WAVEWATCH III (R) version 5.16[END_REF] is used with 36 discrete values of the direction leading to a truncation of the direction of the waves of the order of 10 °. However, the phase-resolved models require finer temporal and spatial discretisation, thus larger computational resources. The resolution of the horizontal grid imposes a limit on the permissible wavelengths, truncating the shortest periods of the wave spectrum.

Deterministic nearshore models should have accurate dispersion and non-linear properties and thus should be governed by the Navier-Stokes equations [START_REF] Raoult | Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments[END_REF], using either an Eulerian approach (e.g., OpenFOAM®, Zhou et al., 2017) or a Lagrangian approach (e.g. the meshless model SPH model, [START_REF] Oger | Aspects théoriques de la méthode SPH et applications à l'hydrodynamique à surface libre[END_REF]. These models are appropriate for local studies, although eventually limited by spurious diffusion or prohibitive computational time [START_REF] Benoit | Analysis of the linear version of a highly dispersive potential water wave model using a spectral approach in the vertical[END_REF]. Simplifications can be made to overcome these problems: for instance, REF/DIF based on the mild slope equation [START_REF] Berkhoff | Computation of combined refraction-diffraction[END_REF] neglects some nonlinear effects; SWASH can be used with depth-averaged equations with appropriate assumptions to introduce the non-hydrostatic pressure [START_REF] Zijlema | SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters[END_REF]; Boussinesq, [START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF] or [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF]) models take partially into account non linear or dispersive effects (e.g., Bonneton et al.,2011a, 2011b[START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF].

Within the limits of the above simplifications, phase-resolved models should theoretically be capable of simulating the low-frequency 3D processes usually modelled by circulation models, since their equations are basically the same, except for the hydrostatic assumption generally made in circulation models [START_REF] Klingbeil | The numerics of hydrostatic structured-grid coastal ocean models: state of the art and future perspectives[END_REF]. As a result, the phase-resolved wave model is able to generate residual currents such as rip currents and therefore does not require coupling with a circulation model to simulate the coupling between waves and low frequency fields. However, it requires a significant horizontal and vertical resolution (greater than the wavelength of the waves) which limits the size of the modelled domain. To date and to our knowledge, most of the applications of the phase-resolved models seem to concern domains whose horizontal dimensions are of the order of a few kilometres in each horizontal direction [START_REF] Yoon | Non-hydrostatic Modelling of Waves Transformation and Rip Current circulation: A case Study for Haeundae Beach, Korea[END_REF]. The constant rise of computers capability suggests that phase-resolved models could now be used on larger domains covering a significant part of the continental shelf. Meanwhile, ocean circulation models currently used at coastal scales could be adapted to deterministic wave modelling with, (among other objectives), the goal of linking waves and continental shelf circulations like wind-induced eddies [START_REF] Petrenko | Barotropic eastward currents in the western Gulf of Lion, north-western Mediterranean Sea, during stratified conditions[END_REF], or to better understand momentum and turbulence transfer in the surface layer [START_REF] Deigaard | Wind generation of waves: Energy and momentum transfer -An overview with physical discussion[END_REF]. This assumes, however, that their equations are suitable for phase-resolved wave modelling. One of the essential points is to be able to use these models in non-hydrostatic mode. A number of articles deal with this question, for the MITgcm model [START_REF] Marshall | Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modelling[END_REF], POM [START_REF] Kanarska | A non-hydrostatic numerical model for calculating freesurface stratified flows[END_REF], ROMS [START_REF] Kanarska | Algorithm for non-hydrostatic dynamics in the Regional Oceanic Modelling System[END_REF], GETM [START_REF] Klingbeil | Implementation of a direct nonhydrostatic pressure gradient discretisation into a layered ocean model[END_REF], SYMPHONIE (Auclair et al, 2011). However, at the resolutions where these models are most often used, the consideration of non-hydrostatic effects may not be an essential issue [START_REF] Mckiver | Investigation of model capability in capturing vertical hydrodynamic coastal processes: a case study in the north Adriatic Sea[END_REF]. The applications envisaged so far concern more the propagation of internal gravity waves [START_REF] Bordois | Nonlinear processes generated by supercritical tidal flow in shallow straits[END_REF] or the convective processes [START_REF] Paluszkiewicz | Deep convective plumes in the ocean[END_REF] than the high frequency surface waves. One of the obstacles to the deterministic modelling of waves over large domains is the cost of the non-hydrostatic solver [START_REF] Klingbeil | The numerics of hydrostatic structured-grid coastal ocean models: state of the art and future perspectives[END_REF]. As previously mentioned, in the field of phase-resolved wave modelling, the common alternative to this difficulty is to use simplified models adapted to the physical specificities of waves (Bonneton et al.,2011a(Bonneton et al., , 2011b)). In the field of ocean circulation modelling, the most widespread approach to the non-hydrostatic problem, on the contrary, is to take into account the whole complexity of the problem. It leads concretely to the resolution of a Poisson equation for the non-hydrostatic pressure. This is a global approach insofar as this leads to a system of linear equations making all the points of the numerical domain interdependent. Its resolution is in principle costly and currently constitutes a challenge in the perspective of large numerical domains of several thousand points in each direction (Roulet et al., 2017). The sigma coordinate widely used in coastal modelling also increases the cost of the resolution for two main reasons. The first reason is that the time variation of the position of the vertical levels (related to the movement of the free surface) forces to update regularly the coefficients of the main matrix of the system and its possible preconditioning. The possible representation of the wetting-drying zones adds a similar difficulty. The second reason is that the sigma coordinate transformation increases the number of nonzero coefficients of the principal matrix. Roulet et al. 2017 propose to reduce this number by reformulating the momentum equations. The same authors also evoke the perspective of the multigrid approach to significantly reduce the cost of calculation. An application in the vertical direction is proposed by [START_REF] Shi | Pressure Decimation and Interpolation (PDI) method for a baroclinic non-hydrostatic model[END_REF].

The alternative to the global approach is the local approach based on the direct resolution of non-hydrostatic pressure. Various methods have been proposed so far (Johns 1991[START_REF] Klingbeil | Implementation of a direct nonhydrostatic pressure gradient discretisation into a layered ocean model[END_REF][START_REF] Bordois | Nonlinear processes generated by supercritical tidal flow in shallow straits[END_REF][START_REF] Lee | Application of the artificial compressibility method for turbulent open channel flows[END_REF]. Among them, Lee et al, 2006 (hereafter L06) ACM (artificial compressibility method) method has a reasonable numerical cost and its implementation in a model is easy. A method derived from L06 is implemented in the SYMPHONIE circulation model [START_REF] Marsaleix | Energy conservation issues in sigma-coordinate free-surface ocean models[END_REF] and tested in the aforementioned context, namely the high frequency surface waves and the associated residual circulations. The purpose of the present paper is to assess this model on the principal issues of wave modelling. This concerns on the one hand the dispersion properties, the vertical shear of the current, the asymmetrical wave transformation at the approach of the breaking conditions, the decay of the waves in the breaking zone. On the other hand, this concerns the residual waves-related circulation such as the stokes transport, the alongshore coastal drift, the mean sea level, rip currents and their feedback on the waves.

The article is organized as follows: chapter 2 deals with the implementation of the ACM method in the ocean circulation model, chapter 3 gives an analytical solution verifying the ACM equations for low amplitude waves, chapter 4 describes academic and laboratory test cases and finally chapter 5 provides a summary and the conclusions.

Model description

Equations

For the sake of clarity the problem is first presented in a simplified way, in a vertical 2D plane and presenting only the terms that dominate wave physics:

= - - (1)
= -

where q is the non-hydrostatic pressure times the water density, u and w the horizontal and vertical current, the sea surface height. As in L06 the hypothesis of incompressibility is retained to deduce the surface elevation:

= - (3)
Equation ( 2) is in practice not calculated by the model. Instead, the combination of ( 1) and

(2) that leads to the Poisson equation is used:

+ + = - - - (4)
The left-hand side of the Poisson equation is ideally zero in the case of the incompressible hypothesis. Without questioning this hypothesis of incompressibility, the particularity of the L06 method is however to allow the left-hand side of (4) to be small and not to equal zero.

Before detailing the approach of L06, some comments can be made on the possibility of not cancelling strictly the left-hand side of (4). It can be noted, for example, that non-hydrostatic incompressible models can use iterative solvers whose convergence is de facto imperfect, precluding a strict cancellation of the left-hand side of (4). This approximation is currently accepted because the physical coherence is preserved if the convergence of the solver is sufficiently precise. We can also recall that a model such as the one used in this study does not formally represent (4) but its discrete approximation in the finite difference sense. Not strictly cancelling the left-hand side of (4) induces an error whose order of magnitude can be compared to the discretization errors of the other terms. For example, given the numerical scheme used here, the error made on the horizontal laplacian of q at the right-hand side of (4) is of the order of (APPENDIX A):

≈ ( 5
)
where k is the wavenumber. This discretization error makes it possible to consider a possible error on the left-hand side, thus relativizing the question to the transformation of the left-hand side of (4) into a term of artificial compressibility as described by L06 (see equation 2 in L06). In practice, the ACM method of L06 leads to the reformulation of (4) as follows:

! = + + (6)
Assuming a sinusoidal form for q , the artificial compressibility term at the left-hand side of ( 6) is of the same order as

" ! = # !
(c is the celerity of the waves), which, under certain conditions, can become of an order of magnitude comparable to (5) (ie

# ! << 1.
This point is closely examined later in the article.

In fact, ideally, in compressible theory, the $ constant should be the speed of the acoustic waves, actually much higher than c. However, as pointed out by [START_REF] Mahadevan | A nonhydrostatic mesoscale ocean model parts 1, 2[END_REF], this would have the consequence of requiring an extremely small time step, highly increasing the computational cost . This type of approach therefore generally uses a much lower $ value in order to increase the time step and thus makes the cost of the calculation affordable. It is therefore noted that the true compressibility of the ocean is not respected by this approach which relies more on numerical considerations than on a real physical justification (resulting in the expression "artificial compressibility" used by L06). On the other hand, apart from the equation ( 6), the other equations of the model strictly respect the incompressibility framework. Although greatly diminished by the computational time constraint just mentioned, the value of $ must nevertheless remain sufficiently large to allow rapid adjustment of equation ( 6) to variations in the hydrostatic forcing terms. This property is essential for the accuracy of the non-hydrostatic behavior of the model. In practice, $ must remain larger than the phase speed of gravity waves [START_REF] Bordois | Nonlinear processes generated by supercritical tidal flow in shallow straits[END_REF].

The resolution of the Poisson iterative equation ( 6) can be done with the same time step as the other equations. In practice its determination is done from the stability criterion corresponding to the forward-backward scheme, but considering $ rather than the speed of the gravity waves as it is larger. The time step is therefore smaller than the one of a hydrostatic model in similar conditions, but the calculation nevertheless remains affordable (as we will see in the following) because equation ( 6) is computed only once per time step of the full model.

Implementation in the SYMPHONIE model

The model used is the ocean circulation model SYMPHONIE [START_REF] Marsaleix | Energy conservation issues in sigma-coordinate free-surface ocean models[END_REF]).

As we focus on surface waves, the effect of temperature and salinity is neglected. For the sake of clarity the equations are given in a Cartesian coordinate system but the model actually uses a system of generalized sigma coordinates. The momentum equations are

given by:

= - - - + '( - -) + * + + , (7)
= - - - -' - -) + * + + , (8)
! = - * - * + + + (9)
Where * and * represent hydrostatic terms (in practice all the terms of the righthand side of ( 7) and (8) except ) and ) ). At this stage, r at the right-hand side of ( 7)

and ( 8) is equal to one but we will see later that a value other than 1 can be used to improve the accuracy of the dispersion relation and currents of short waves. The incompressibility hypothesis gives the vertical velocity and elevation of the surface:

+ + = 0 (10) = - - ( (11)
The vertical diffusion coefficient is calculated by the k-epsilon scheme [START_REF] Rodi | Examples of calculation methods for flow and mixing in stratified fluid[END_REF] and (, , , ) are breaking terms parameterized as a horizontal diffusion:

, = 0 1 + 20 1 + 3 (12) , = 0 1 + 20 1 + 3 (13)
The eddy viscosity of breaking 0 1 is calculated from a breaking criterion which will be discussed in a section to follow. Bottom and surface boundary conditions are as follows: at the surface q=0, ( , )

= 0 and at the bottom = 0 , 7 8 * + ( , )

= (, 1 , , 1 ) where 7 8 is the sea water density and (, 1 , , 1 ) are the components of the bottom stress. The bottom stress is parameterized as in [START_REF] Blumberg | A description of a three-dimensional coastal circulation model[END_REF], from a quadratic relationship of the bottom current and a drag coefficient derived from a roughness length.

The equations are calculated with a forward-backward time scheme detailed in APPENDIX B. Considering that short surface waves are intrinsically three-dimensional, the barotropic-baroclinic time splitting commonly used in circulation models for calculating separately barotropic and baroclinic velocities [START_REF] Blumberg | A description of a three-dimensional coastal circulation model[END_REF] can not applied here. In practice, all the equations are calculated with the same time step.

The time step

The time step is determined from the stability criterion specific to the forwardbackward scheme. In the case of the hydrostatic model, the criterion of stability with respect to surface waves would be:

9 +: = √ <= >?@ ( 14
)
where dx is the horizontal resolution, ℎ +: is the maximum depth of the domain and therefore < ℎ +: is the maximum theoretical phase speed expected for surface waves. In the non-hydrostatic model, the calculation of the equation for q imposes a more drastic additional criterion since it depends on $ that is greater than < ℎ +: . It is not trivial to find a simple expression similar to (14) but, unsurprisingly, it appears in simulations performed by the authors that the maximum value of the time step allowed by the non-hydrostatic model was close to √ B , ie (14) with < ℎ +: replaced by $. In practice, in all the simulations presented, the time step of the non hydrostatic model is given by a fraction (70%) of this criterion to include other stability criteria than gravity waves such as those resulting from non-linearities, namely:

9 = 0.7 √ B (15)
The determination of $ is therefore related to the determination of the time step. The larger the first, the smaller is the second. The quantity N which represents the ratio of the theoretical time step of the hydrostatic case to the time step of the non-hydrostatic case, is introduced:

E = >?@ (16)
which, using ( 14) and ( 15), is equivalent to

E = 8.F ! <= >?@ (17)
N is a key parameter of the problem. The larger it is, the smaller is the model time step and the better is the accuracy expected on the non-hydrostatic term. Since 9 +: , is the time step that would normally have the hydrostatic model in similar conditions, the ratio N is a good indicator of the relative additional cost associated with the non hydrostatic effect. The evaluation of this additional cost must also take into account that an additional equation (for q) has been added to the equations of the hydrostatic model. The latter essentially counts 4 three-dimensional equations, two for the horizontal components of the current, one for the turbulent kinetic energy and one for its dissipation rate. The relative numerical excess cost of the equation for q is therefore of the order of E G = 1.25 E .

Analytical solutions

Analytical solutions verifying model equations in the context of linear theory are now presented. We consider a 2D vertical plan, a linearized version of (7),

= - -) (18)
the continuity equation,

+ = 0 (19)
and the non-hydrostatic pressure equation ( 6). The elevation of the surface is of the form = 8 J K(L " ) . The non-hydrostatic pressure satisfying equation ( 6) is of the form:

= -( M ⁄ ) + 2O 8 J K(L " ) cosh (M ( + ℎ)) (20)
with

M = 1 - " L ! (21)
Where α appears in the left hand side of ( 6) and in the definition of N calculated by (17).

Using (20) and the surface condition q=0 leads to:

= = U (-1 + cosh (M ( + ℎ))/cosh (M ℎ)) (22)
Using ( 22) and (18) gives the horizontal current:

= L " W1 + X U YZ[\ (UL( G )) YZ[\ (UL ) -1 ] (23)
Using ( 23) and the continuity equation ( 19) to derive the vertical velocity, then applying the surface condition ^= , leads to the dispersion relation:

_ = 1 - X U ℎ + X U ` 9abℎ(M ℎ) (24)
Note that if r=1 and if N (therefore α) tends to infinity (leading to the cancellation of the lefthand side of ( 6) and M = 1) the dispersion relation ( 24) becomes equivalent to the usual dispersion relation [START_REF] Kinsman | Wind waves: Their Generation and Propagation on the Ocean Surface[END_REF]:

_ = 9abℎ( ℎ) (25)
while, using ( 23) and ( 25), we find the usual solution for the horizontal current [START_REF] Kinsman | Wind waves: Their Generation and Propagation on the Ocean Surface[END_REF]:

= _ YZ[\ (L( G )) [cd\ (L ) (26)
This result shows that the ACM method becomes equivalent to the classical nonhydrostatic incompressible method if N is sufficiently large. In the following, we distinguish the dispersion relations of the incompressible non-hydrostatic linear theory (ie ( 25)), of the non-hydrostatic ACM theory (ie ( 24)) and of the hydrostatic theory (ie ( 27)): 24) and ( 25) are actually equivalent if

_/ = < ℎ (27)
) = M e fgd\(L ) L fgd\(UL ) UL (28)
Horizontal currents ( 23) and ( 26) are identical on the surface without condition on r and µ , while near the bottom they are equivalent if Expressions ( 28) and ( 29) suggest the possibility of optimizing r according to the local values of the depth and the dominant period of waves. This point will be addressed in future applications. At this early stage of our study, r will be considered constant. Figure 2 shows the absolute value of the relative error on the wavelength,|i jkl -i| i ⁄ , as a function of the period, i jkl being deduced from (24) and L , the "true" wavelength, from (25).

) = M YZ[\ (L ) -1 YZ[\ (UL ) -1 (29)
Relative errors are more important for small periods and for large bathymetry values. Errors of the order of several % are expectable with r = 1. On the other hand, it is possible to reduce the error level below 1%, for all the values of period and bathymetry considered in Figure 2, with r = 0.9845.

The fact that ( 28) and ( 29) are not identical shows that it is not possible to strictly cancel the error on the current and on the dispersion relation with the same value of r.

However as shown in Figure 3 for h=50m and a period of 10s, ( 28) and ( 29) give very similar values of r so that it is possible to find a value of r that reduces both the current and dispersion relation errors (this will be discussed in the next section). Finally we note that a compromise on these two possible values of r is all the easier to find if N is large because ( 28) and ( 29) converge both to r=1. In conclusion, the analytical solutions corresponding to the linearized equations of the ACM method show a satisfactory representation of the non-hydrostatic effect in the range of periods concerned by the future applications of the model to the near shore zone. The accuracy however depends on the setting of parameters such as N and r. Since the accuracy of a numerical model also depends on many other factors related to the discretization of equations, such as spatio-temporal resolution and relevance of numerical schemes, the above theoretical analysis is now completed with numerical experiments.

Test cases

Wave propagation in the context of linear theory

In this section, we check the ability of the model to correctly represent the physical characteristics of the simulated wave thanks to a series of test cases consisting in exploring the dispersion properties of the waves. We consider here low amplitude waves (1 mm) in order to be within the framework of the linear approximation for which simple analytical solutions exist [START_REF] Kinsman | Wind waves: Their Generation and Propagation on the Ocean Surface[END_REF]. The objective is to define the values of N which lead to a reasonable compromise between the computation cost (consideration encouraging to lower N) and the precision of the dispersion properties (consideration encouraging to increase N).

We use a domain with a flat bottom at 50 m depth, a periodic condition in the direction normal to wave propagation. The horizontal resolution is 5m and 20 vertical levels are used.

The waves are generated at the entrance of the domain, while the fields are at rest at the initial state. The opposite boundary is closed but the dimension in the direction of propagation is large enough, so that the waves do not have time to reflect on it before the end of the simulation (10 minutes). This first numerical experiment is carried out with r = 1.

The open boundary conditions are derived from [START_REF] Marsaleix | Considerations on Open Boundary Conditions for Regional and Coastal Ocean Models[END_REF]. As recommended by [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF], an outgoing field is defined as the difference between the total field and the incoming field. We thus define anomalies of currents, sea surface height and non-hydrostatic pressure, m = -, m = -, m =where , , are prescribed by the analytical solution of the short wave linear theory. We then assume that the incoming short waves do not reflect at the coast because of dissipation and that waves propagating offshore are at lower frequencies (they may consist of infra-gravity waves if the conditions of their generations are met [START_REF] Bertin | Infragravity waves: From driving mechanisms to impacts[END_REF]). The radiative conditions, m = ± m /o and p ′ p9 ⁄ ± op ′ pr ⁄ = 0, are therefore applied using the long wave phase speed o = < ℎ.

The simulation is analyzed in order to extract the modelled wavelength and to normally expected. To our knowledge, near shore modelling studies mostly focus on areas where bathymetry does not exceed a few tens of meters. We thus assume that 50m could be a plausible upper limit for realistic models of the near shore zone (at least for the time being). In the following, the case h = 50m is the subject of greater attention.

The model errors are generally similar to those predicted by the ACM linear theory (Figure 4a and Figure 4c). A notable difference is that errors tend to zero when N increases in the case of the theory but not in the case of the model. This is explained by additional sources of error, of numerical origin, in the case of the model. These numerical errors are more important for small periods. For example, for a period of 6s, the residual error for large values of N is around 3% for h = 50m and around 2% for h = 10m (Figure 4a and Figure 4b).

The agreement between the model and the ACM linear theory is less good for h = 10m, especially for a period of 6s (Figure 4b and Figure 4d). The relative predominance of numerical errors is here related to the decrease of the wavelength caused by the decrease of h. For a period of 6s and h=10m, the wavelength is indeed around 48m, i.e. less than 10 grid points while the wavelength is around 56m i.e. more than 11 points when h = 50m.

Actually, when h decreases, causing the wavelength to decrease, the resolution of the horizontal grid eventually becomes insufficient to represent the shortest waves, unless wave breaking occurs before this critical situation.

Model errors are also sensitive to vertical resolution. The accuracy of the vertical profile of q and current depends on the number of vertical levels. The shortest wavelengths are particularly sensitive to this error because the vertical profile of q is more non-linear and is established on a smaller depth. Figure 5 shows that for the 6s period, the relative error is much smaller with 80 levels than with 10 or 20. Using more vertical levels reduces the residual error (but does not completely eliminate it). The error levels obtained with 40 levels or 80 levels are however equivalent. For larger periods (results not shown), the sensitivity of the error to the vertical resolution is smaller, plausibly because the vertical profile of the fields, whose characteristic scale of decay is related to the wavelength, is less sheared. The "minimum" error level, around 1.5% the time-stepping method itself and no longer depends on the non aspect of the error can be examined in the hydrostatic model by re assumptions we made in section

(3) the surface elevation is given by: ror level, around 1.5% in Figure 5 (case 80 levels, N=20

stepping method itself and no longer depends on the non-hydrostatic aspects. This aspect of the error can be examined in the hydrostatic model by re-invoking the simple assumptions we made in section 2.1. Combining the equation ( 1) (now without term

(3) the surface elevation is given by: 23 (horizontal axis). h=50m. (31) so that the modelled phase speed approaches ou1 -≈ o 1 -. The relative error on the phase speed is therefore of the order of = ( v/w) where n is the number of grid points per wavelength. For the shortest wavelengths considered (for which n is around 10), the relative error on c is therefore of the order of 1.5%. Given that the period is imposed by the forcing, this level of error is therefore the one expected for the modelled wavelength and, in fact, it corresponds fairly well to the residual error for N=20 and 80 levels (see Fig. 5).

In summary, the accuracy of the dispersion properties increases significantly with N until ACM-related errors become small compared to the other error sources mainly the vertical resolution and the order of precision of the numerical schemes used to calculate spatial and temporal derivatives. The determination of N would lead to choosing N around 14 in the case of the experiment considered. Another way to choose N is to define an "acceptable" precision level on the dispersion properties, depending on the period of the simulated waves. For example, a relative error of 3% is obtained for N>7 for a period of 8s, N> 5 for a period of 10s, etc…

Improved accuracy (apart from increasing N)

Distribution of vertical levels

The previous paragraph shows that the accuracy decreases appreciably with the wave period. At the shortest periods, accuracy is more sensitive to vertical resolution.

Increasing the number of vertical levels is an obvious but expensive way to reduce this second source of error. Modifying the vertical distribution of levels without changing their number is an alternative that allows a significant improvement in accuracy at small periods.

In the case of the 6s period, Figure 5 shows that 10 irregularly distributed levels (see Table 1) provide a much higher accuracy than 10 equidistant levels, and makes it comparable to the precision obtained with 20 equidistant levels. The vertical distribution proposed here is such that the accuracy is also improved or maintained for the other periods considered It should be noted that, in a realistic domain, the generalized sigma coordinate makes it possible to adapt the vertical distribution to the bathymetry. For example, in areas where bathymetry is low relative to wave wavelength, the vertical profile of current and pressure is less sheared and a more regular distribution of the levels is preferable to the irregular distribution of Table 1. Given the positive influence of increasing vertical resolution on accuracy, the increase in vertical resolution that results from the decrease in bathymetry in the sigma formulation leads to a better representation of dispersion properties when approaching the coast.

(

Influence of horizontal resolution

Increasing the horizontal resolution (without changing the other parameters, especially N) leads logically to decrease the level of error, especially at short periods.

Significantly increasing the horizontal resolution has a real but limited impact on the accuracy because of other sources of error such as the lack of vertical resolution. The period 6s is more critical than the other periods tested because, as noted above, the number of grid points per wavelength is particularly low. Doubling the resolution allows about 22 grid points per wavelength (for h=50m). Fig. 4 shows that a relative error <3% is obtained with N> 10 for a horizontal resolution of 2.5m, whereas this same level of error requires N> 14 for a horizontal resolution of 5m.

Tuning the non hydrostatic pressure gradient

A modification of the pressure gradient scheme is now proposed to use modera values while maintaining an acceptable level of accuracy on dispersion properties. We consider values of N around 5. The proposed modification is suggested by the observation that the error on the wavelength is to the fact that the modelled phase speed is systematically lower than the theoretical phase speed. In other words, the non This conjecture can be interpreted as an overestimation the pressure gradient in agreement with the conclusions drawn from the analytical solution of the linearized ACM equations ( moderating the non hydrostatic press slightly less than 1 (see the non ( 8)). Tuning the non hydrostatic pressure gradient A modification of the pressure gradient scheme is now proposed to use modera values while maintaining an acceptable level of accuracy on dispersion properties. We around 5. The proposed modification is suggested by the observation that the error on the wavelength is systematically negative (Figs 4 and5). phase speed is systematically lower than the theoretical phase speed. In other words, the non-hydrostatic effect excessively decreases the phase speed.

This conjecture can be interpreted as an overestimation of the non-hydrostatic component of in agreement with the conclusions drawn from the analytical solution of linearized ACM equations (section 3, Figure 1). The proposed improvement consists in moderating the non hydrostatic pressure gradient by multiplying the latter by a constant 1%. Each curve corresponds to a period: 14s (magenta), 12s (cyan), 10s (blue), 8s (green), 6s (red). The dashed vertical line which all relative errors are between +/-1% r Fig. 6a, obtained with N=5 (dt = 0.032s), shows that the precision is improved when r decreases to a certain level below which the sign of the error is reversed (becoming positive) and the error increases again. As expected from the ACM linear theory, the value of r that exactly cancels the error is not the same depending on the period. It is around r=0.982 for the 6s period and smaller for the other periods. The general behaviour is nevertheless quite similar for all the periods and one can identify a unique value of r which makes it possible to reduce the error below 1% for all the periods (Figure 6a, r=0.98), which corresponds to a very significant reduction of the error, especially for small periods (error of 8% and 5% for periods of 8s and 6s respectively, with r=1). The general behaviour of the model is in good agreement with the analytical solution (Figure 6b). However, since the model has other sources of errors than those of the theoretical model, Figures 6a and6b are not strictly the same and the best compromise for the value of r would be slightly larger according to (24), around 0.984 (Figure 6b). Although the theoretical model could certainly be used for a default setting of r, Figure 6a shows that errors can be further reduced with an empirical setting. This numerical estimate of the best compromise on r is therefore used in the following.

As pointed out previously, a realistic numerical domain is characterized by variable bathymetry. As 50m should be an upper limit of a near shore zone model in our future applications, we now verify that the value r=0.98, which reduces the error below 1% for all periods when h=50m, remains relevant for lower values of the water depth. These tests are done keeping the same time step (dt=0.032s) and N=5. Figure 7, obtained for h=20m, shows a precision greater than the simulation with h=50m, under the combined effect of the increase of the vertical resolution (the same number of vertical levels distributed over a shallower water column) and also the decrease in the non-hydrostatic effect caused by the decrease in the ratio between the horizontal scale and the vertical scale. If Fig. 7 suggests that the best value of r would be around 0.978 (the minimum relative error is indeed close to this value for most of the periods considered), the proposed value of 0.98 for h=50m remains relevant for shallower depths because it reduces the relative error associated with the different periods between -0.25% and 0%. Note finally that the choice of the 50m value as the ma is somehow arbitrary. There are certainly other ways of defining the shore zone, such as the depth at which the waves are modified by the seabed. In this case, a maximum value based on the relative d relevant for shallower depths because it reduces the relative error associated with the 0.25% and 0%. but h=20m. The 2 horizontal black lines indicate a relative error of + or 0.25%. The dotted vertical line corresponds to r = 0.98 for which all relative errors are between Note finally that the choice of the 50m value as the maximum value of the bathymetry arbitrary. There are certainly other ways of defining the offshore limit of the near zone, such as the depth at which the waves are modified by the seabed. In this case, a maximum value based on the relative depth, kh, may be preferable.

relevant for shallower depths because it reduces the relative error associated with the

The 2 horizontal black lines indicate a relative error of + or -0.25%. The dotted vertical line corresponds to r = 0.98 for which all relative errors are betweenximum value of the bathymetry offshore limit of the near zone, such as the depth at which the waves are modified by the seabed. In this case, The ability to strongly reduce, or even cancel, the error via the that this parameter can be used to correct errors beyond the single error caused by the ACM method. The errors associated with the vertical dis can therefore be compensated by this method. We thus examine the possible undesirable consequences of this compensating effect on other impor especially since, according to the theoreti current and the relation of dispersion is not obtained strictly with the same value of ( 28) and ( 29)). Recall also that t interest, since the intensity of the current near the bottom is crucial for the determination of the bottom friction, the resuspension of sediments near th (among others). Figure 8a shows that for the period 6s, the current more accurate than with vertical shear), the value r=0.986 profile predicted by the linear made for the period 10s (Figure 8 reduce the error on the dispersion relation and the current is which shows that the values of r The ability to strongly reduce, or even cancel, the error via the r that this parameter can be used to correct errors beyond the single error caused by the ACM method. The errors associated with the vertical discretization or the time can therefore be compensated by this method. We thus examine the possible undesirable consequences of this compensating effect on other important properties of wave physics especially since, according to the theoretical model, the cancellation of the errors on the current and the relation of dispersion is not obtained strictly with the same value of Recall also that the shape of the vertical profile of the current is of particular ce the intensity of the current near the bottom is crucial for the determination of the bottom friction, the resuspension of sediments near the bottom, the transport of sand shows that for the period 6s, r=0.98 leads to a vertical the current more accurate than with r=1 (the latter nevertheless correctly reproduces the r=0.986 being the one that would lead to a profile closest to the profile predicted by the linear non-hydrostatic incompressible theory. A similar conclusion is (red line), r = 0.98 (black dotted line), theoretical profile from the incompressible linear theory (continuous black line). For the sake of clarity the profiles corresponding to r = 0.986 (period 6s) and to to r = 0.976 (periode 10s) are not shown because they merge with the theoretical profiles.

Test case with bi-periodic domain

We are now replacing open and closed boundary conditions at both ends of the domain in the propagation direction by periodic conditions. The interest of this bi-periodic domain is that it avoids reflections at the boundaries that could degrade the results. The simulations can therefore be long, which is of interest for testing the stability of calculations in the long term. The depth is 50m. The domain has 510 points in the direction of propagation corresponding to a whole number of wavelengths, here 17, for a wavelength of 150m, a theoretical period close to 10s and a phase speed close to 15m/s. The fields are initialized from the linear theory. The amplitude of the surface elevation is small (1 mm) to remain within the framework of the linear theory. Nonlinear dissipation forces such as friction (using a quadratic parameterization) or breaking are therefore negligible. In this particular case, the Laplacian temporal filter derived from [START_REF] Marsaleix | Alternatives to the Robert-Asselin filter[END_REF], used in the equations for u, v, q to suppress the numerical mode, is the main source of dissipation. The following test quantifies its adverse effect on the long-term physical mode. Figure 9 shows the amplitude of the surface elevation at different times. Given the value of the phase speed, a distance of the order of 10km (cross shore size of realistic domains that we plan to use with this model) is covered in about ten minutes. Figure 9 shows that the attenuation is negligible over such a short time duration. An attenuation of 10% would require at least 6 hours, during which the waves would travel a distance of about 300 km. In our future realistic applications, the open boundary condition will not be located so far from the coast. The cross-shore domain dimension should be at least an order of magnitude lower, so this numerical attenuation should be less than 1%, induced by bottom friction, turbulence or wave breaking. The Hamilton and Ebersole laboratory experiment [START_REF] Hamilton | Establishing uniform longshore currents in a large-scale sediment transport facility[END_REF]Ebersole 2001 hereafter HE01, Svendsen et al 2003) focuses on the propagation of waves over a uniform linear slope beach, their breaking, and the alongshore current establishment. This case test is of particular interest to us to assess the performances of the model on the important question of the parametrization of wave breaking. Incoming waves are monochromatic, with a period of 2.5s (TEST 6N in HE01). The amplitude of the excited waves at the open boundary is adjusted to obtain a wave height corresponding to the conditions of TEST 6N, namely 0.19m at the most offshore measurement points (see the observation at 3 meters from the incoming boundary on Fig. 10). The direction of propagation is 10° from the direction normal to the beach. A description of the basin and its bathymetry is given in [START_REF] Svendsen | Modelling waves and currents at the LSTF and other laboratory facilities[END_REF] (see also Fig. 10). The particularity of the basin lies in its pumping and recirculation system made of 20 pumps that allow having periodic conditions on the lateral boundaries. In our simulation, incoming waves are specified from the linear theory at the offshore boundary. Periodic boundary conditions are applied at the lateral boundaries of the model. The dimension of the domain in the direction parallel to the beach is consistent with the condition of periodicity and the wave propagation direction of 10° in the deep part of the basin. The horizontal resolution is 0.1m and the grid has 10 vertical levels. The dimension in the cross-shore direction is 18m and the deepest bathymetry, offshore, is 0.66m. The maximum wavelength, in the deep part of the basin, is of the order of 6m. A steady state is reached after about 150s. The results presented are obtained at the end of a simulation of 500s. The roughness of the bottom corresponding to this experiment is a priori not mentioned by its authors. A roughness length of 3 × 10 z { was empirically fixed in our model in order to obtain a satisfactory agreement with the vertical profile of the measured alongshore current.

The basin dimensions are small compared to a real case. If the ratio of the horizontal scale on vertical scale is used to deduce the characteristics that this experiment would have with dimensions 50 times larger, close to the dimensions of the test cases presented in the previous section, we obtain a horizontal resolution of 5m, a maximum depth of 30m, a maximum wavelength of the order of 300m and a period of 125s. We are thus in conditions for which the problems of accuracy examined in the preceding sections are minor and make it possible to use the ACM method with a small value of N and a coefficient of reduction of the gradient of q close to one. In practice we use N=4, r=1 (as mentioned previously, we actually chose this test case for the wave breaking issue rather than the non-hydrostatic specificities discussed in the previous test case).

Wave breaking is parameterized through the horizontal diffusion of the momentum (13). The determination of the mixing coefficient, 0 1 , is an important aspect of this parameterization. Before going into the details of its calculation, we examine the order of magnitude necessary for the decay of the waves inside the breaking zone. To do this, we consider an energy balance between the gradient of the potential energy flux ) and dissipation due to horizontal diffusion -0 1

(=
(for the sake of clarity we limit ourselves to unidirectional propagation along the Ox axis). We hypothesize that within the wave breaking zone the wave height is proportional to the water depth [START_REF] Miche | Mouvements ondulatoires de la mer en profondeur croissante ou décroissante. forme limite de la houle lors de son déferlement. application aux digues maritimes. Troisième partie. Forme et propriétés des houles limites lors du déferlement[END_REF]:

| = } ( 32
)
where H is the wave height. We consider here that the parameter } is a constant of the order of 0.8 (it can nevertheless be noted that the expression of } an be more sophisticated, as in [START_REF] Guérin | Impacts of wave-induced circulation in the surf zone on wave setup[END_REF]. For convenience (and considering that shallow depths, where breaking mostly occurs, allow this approximation) the hypothesis of hydrostatic equilibrium is made in the following reasoning (although the simulations presented are still non-hydrostatic). We thus consider that the current and the surface level are connected to each other according to = = # η where o = < ℎ. It is finally assumed that the fields have a sinusoidal behaviour and that the wave amplitude 8 satisfies (32) (ie ~ = }). Averaging over a period of waves, the gradient of the potential energy flux is (APPENDIX D):

(= ) = e/ } e • ℎ / (33)
and dissipation due to horizontal mixing is

€ 1 = -0 1 • L=' # (34)
The equilibrium of ( 33) and ( 34) leads to an estimation of the average of the mixing coefficient over a period of waves:

< 0 1 >= e • e/ ℎ / " v (35)
The latter increases with the slope of the bathymetry and vanishes when the thickness of the water column is zero. The expression ( 35) is only an intermediate step towards a more complete formulation taking into account a triggering criterion based on the analysis of the instantaneous modelled fields. But before addressing this aspect, a preliminary experiment is carried out to test the relevance of the order of magnitude given by (35). Practically ( 35) is used as the value of the horizontal diffusion coefficient in (13) inside the wave breaking zone. The latter is arbitrarily defined from a threshold on bathymetry (h<0.28m) deduced from observations of wave heights (HE01). In the pre-breaking zone (h> 0.28m) 0 1 = 0. Figure 10 shows the simulated wave height and mean elevation, and the corresponding observations of HE01. As the laboratory experiment, the simulation lasts 500s. The surface height is averaged over the last 50 seconds of the simulation. The wave height is calculated from the difference between the maximum and minimum values of the surface level over the last 50 seconds of simulation. The agreement is satisfactory. The wave breaking trigger point is unsurprisingly well positioned since it is imposed by the criterion h<0.28m. The wave height is slightly overestimated in the breaking zone but the slope of decrease is parallel to that of the theoretical decay ~ = } = 0.8 (dotted line in Fig. 10). It should be noted that several simulations were performed by increasing or decreasing slightly the value of < 0 1 > and that Figure 10 corresponds to a 5% increase in the value predicted by (35).

Breaking criterion

Practically ( 35) is not really appropriate because it requires to arbitrarily specify the location of the onset of breaking. In realistic case, the breaking point of the breaking wave varies over time, in particular according to the characteristics of the waves and more generally according to the surrounding hydrodynamic conditions. An extension of (35) taking into account the instantaneous fields is therefore proposed:

0 1 =< 0 1 > ' 1 (36)
Where 0 1 is the instantaneous coefficient, < 0 1 > is the average coefficient (35) and ' 1 is a function that triggers the wave breaking according to the instantaneous fields. The average of ' 1 over a period of wave is expected to be close to 1 so that the viscosity is on average close to < 0 1 >. One consequence of this property is that 0 1 can be momentarily very high if ' 1 is negligible over a large fraction of the wave period (so that the average of 0 1 stays close to < 0 1 >). Many authors use the ratio of the wave height to the water depth (the } constant (32)) as breaking criterion (Rusu and Soares, 2012). Others use criteria based on the horizontal gradient of the modelled fields [START_REF] Lynett | Nearshore wave modelling with high-order Boussinesq-type equations[END_REF]. For example, [START_REF] Zijlema | SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters[END_REF] propose to formulate eddy viscosity from Prandtl's length and a turbulence equilibrium hypothesis which in practice links eddy viscosity to the horizontal gradient of velocities (see their equation 12). Roeber and Cheung (2012) also use the velocity gradient (see their equation 70) but in an approach that is more akin to a criterion on the Froude number, with < ℎ as an approximation of the phase speed. These authors emphasize that this criterion has the particularity to be more influential in areas where the slope of the waves is important in contrast to the † ℎ ⁄ criterion which favours the wave crests. The criterion on the gradient is therefore in principle influencing a larger proportion of the breaking zone. From the point of view of (36), an approach of this type should lead ' 1 to act on a larger time interval and thus limit the maximum values of 0 1 . On the other hand, a criterion on the wave crest should lead ' 1 to be small over a large fraction of the wave period, with the aforementioned consequence of momentarily strong 0 1 values. This consideration may be of numerical importance as there is a stability criterion which limits the eddy viscosity value for a given time step [START_REF] Zijlema | SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters[END_REF]. In the particular case of this experiment, we find that a wave breaking criterion based on † ℎ ⁄ can cause the eddy viscosity to momentarily exceed the maximum value allowed by the CFL criterion (ie 0 1 > 0 1 +: = ). One possibility of overcoming this problem is either to reduce the time step of the model, or to use a time splitting technique with a smaller time step for the calculation of the wave breaking without impacting the time step of the other equations of the model. From a numerical point of view, the approach proposed by Roeber and Cheung ( 2012) is interesting because their parameterization of the wave breaking should act on a longer duration, limiting the maximum values of the breaking viscosity. We choose this approach, mainly for the aforementioned numerical issues, but we propose to add a modification to that method in order to extend the influence of the breaking criterion to an even greater proportion of the breaking zone. This is arguably questionable from the point of view of the realism of the process [START_REF] Kennedy | Boussinesq modelling of wave transformation, breaking and runup[END_REF] but it is an interesting option considering the numerical aspects mentioned above. More specifically, our objective is to avoid reducing the time step because of a too restrictive stability condition related to the wave breaking scheme.

The triggering criterion for wave breaking is inspired from [START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF] but is modified to be active both on the wave crest and between crests. This criterion considers both velocities and their horizontal gradient. The velocity gradient is scaled by the wavenumber so that the two quantities have the same order of magnitude. A wave breaking criterion is then given by the ratio of the sum of these fields to the hydrostatic phase speed (for simplicity the following formula is given in the case of unidirectional propagation):

' 1 = o ‡ ˆ G ‰ Š ‹OE ‹@ = • w/ (37)
Where Ž is the surface current and k is the wavenumber fixed a priori according to the peak period prescribed at the open boundary. The constant o is determined empirically:

ideally, its role is to allow the average value of ' 1 to be of the order of unity over a period of waves (so that the eddy viscosity is on average close to ( 35)). It may be noted that in the idealized framework of the linear theory (where the current would be of the form 8 cos (_9 + r)) the term in the numerator of (37) would simply be 8 , independent of time, and, in a breaking situation characterized by a current of an order of magnitude comparable to that of the phase speed, the constant o would simply be of the order of ( 35). The exponent, n, is used to adjust the selectivity of the criterion. It is in particular chosen so that ' 1 tends more rapidly towards 0 upstream of the breaking zone. The simulations presented in the following use n=2. Taking into account the breaking in the turbulent closure scheme

The turbulence scheme used here is the k-epsilon scheme [START_REF] Rodi | Examples of calculation methods for flow and mixing in stratified fluid[END_REF]. A waverelated turbulence generation term may be taken into account when the circulation model is coupled to a phase-averaged wave model [START_REF] Michaud | Three-dimensional modelling of wave-induced current from the surf zone to the inner shelf[END_REF]. In the present case, this term is not taken into account because the explicit modelling of the waves leads naturally to a significant production of turbulent kinetic energy (TKE) through the vertical shear of the current. The relevance of adding an additional production term to represent the assumed increase in turbulence in the breaking zone is an open question [START_REF] Feddersen | The Effect of Wave Breaking on Surf-Zone Turbulence and Alongshore Currents: A Modelling Study[END_REF] and is therefore examined here. The HE01 experience provides the opportunity to examine this point by providing measurements of the residual onshore current at different depths and at different distances from the edge of the basin. Here, an additional production term, O 1 , representing a fraction of the kinetic energy lost by the current via the horizontal mixing associated with the breaking (ie (36), (37)) , is introduced into the TKE equation of the k-eps scheme:

O 1 = o 0 1 • + + + • (38)
Where o is a constant between 0 and 1. Equation ( 38) implies that a fraction of the breaking energy feeds the vertical mixing. The simulated alongshore residual current is obtained by averaging the alongshore current during the last 50 seconds of the simulation. Figure 12 shows the comparison with the observed current from HE01. Note that the vertical axis has been normalized so that the figure is easily comparable to that of [START_REF] Svendsen | Modelling waves and currents at the LSTF and other laboratory facilities[END_REF] (see their figure 11). Two simulations are made, one with o = 0 the other with o = 0.01. The two simulations generally give rather satisfactory results with the exception of the current at x=7.9m. Looking more closely at Fig. 12, we notice that the simulation taking into account O 1 (red lines) gives a better agreement with the observations. The simulation without (o = 0, green lines) shows indeed a current with a greater vertical gradient than that of the observation, in particular at the points located in the shallowest depths. Increasing the vertical mixing coefficient associated with zone. We therefore retain at this stage that this approach can be used to improve turbu in the wave breaking zone. However, in the state of our knowledge, the setting of (red lines) gives a better agreement with the observations. The simulation without , green lines) shows indeed a current with a greater vertical gradient than that of the t the points located in the shallowest depths. Increasing the vertical mixing coefficient associated with O 1 reduces the current shear in the breaking zone. We therefore retain at this stage that this approach can be used to improve turbu in the wave breaking zone. However, in the state of our knowledge, the setting of arbitrary and would need to be estimated on other experiments, including in situ data such as current profiles or more direct observations of turbulence [START_REF] Costa | Comparison of in situ microstructure measurements to different turbulence closure schemes in a 3-D numerical ocean circulation model[END_REF].

profile of the measured mean alongshore current (black squares), simulated , green lines), simulated with the termO 1 , o G G

. Horizontal axis: m/s.

The case of a barred beach with rip currents 42 (red lines) gives a better agreement with the observations. The simulation without O 1

, green lines) shows indeed a current with a greater vertical gradient than that of the t the points located in the shallowest depths. Increasing the reduces the current shear in the breaking zone. We therefore retain at this stage that this approach can be used to improve turbulence in the wave breaking zone. However, in the state of our knowledge, the setting of o is arbitrary and would need to be estimated on other experiments, including in situ data such lence [START_REF] Costa | Comparison of in situ microstructure measurements to different turbulence closure schemes in a 3-D numerical ocean circulation model[END_REF].

alongshore current (black squares), simulated = 0.01 (red lines).

The interaction of waves and rip currents is an example of the complexity of possible interactions between current and waves [START_REF] Weir | A vortex force analysis of the interaction of rip currents and surface gravity waves[END_REF]. The model is now applied to a laboratory experiment describing the propagation of monochromatic waves over a beach with a sandy bar and channels inducing rip currents (Haller et al, 2002, Haas andSvendsen, 2002). This experiment has already been simulated with hydrostatic circulation models coupled with a phase-averaged wave model [START_REF] Haas | Comparing a quasi-3D to a full 3D nearshore circulation model: SHORECIRC and ROMS[END_REF][START_REF] Haas | Quasi three dimensional modelling of rip current systems[END_REF][START_REF] Michaud | Three-dimensional modelling of wave-induced current from the surf zone to the inner shelf[END_REF][START_REF] Bruneau | Modélisation morphodynamique des plages sableuses[END_REF]. Like these previous authors, we are interested here in TEST B of [START_REF] Haller | Experimental study of nearshore dynamics on a barred beach with rip channels[END_REF]. The basin and its bathymetry are described in detail in the aforementioned articles. The dimensions of the basin are approximately 17m (cross-shore direction) by 18m (alongshore direction). The cross-shore slope is 1/30. A bar, located a few meters from the coast, is cut by 2 channels. Ideally, the two channels should have an identical shape and the bathymetry of the entire basin should respect an alongshore symmetry: the central bar is twice as long as the two half-bars touching the side walls and outside the bar the alongshore gradient is supposed to be nil. However, as pointed out by [START_REF] Haas | Quasi three dimensional modelling of rip current systems[END_REF], an important feature of this bathymetry is that it actually has irregularities compared to its idealized design, resulting in dissymmetry of the flow and increasing instability of the currents. The two rip currents that form at the exit of the channels are for example clearly different in their form and their intensity, and their respective temporal variations seem rather decoupled. A bathymetry interpolated from a survey of the wave basin is used here to better respect the real conditions of the experiment. The horizontal resolution is 0.06m and the vertical grid has 10 sigma levels. The time step is 1.1 × 10 seconds. Monochromatic waves of period 1s propagate in the direction perpendicular to the coast (Fig. 13). The amplitude of the excited waves at the open boundary is adjusted to obtain a wave height corresponding to the conditions of TEST B, namely about 4.1 cm at the most offshore measurement points (see Fig. 14 around x=8m). Their wavelength at the entrance of the domain (where the depth is about 0.38m) is about 1.4m (Fig. 13) in good agreement with the linear non hydrostatic wave theory. It may be noted that the hydrostatic approximation under the same conditions would lead to a wavelength of about 1.9m, ie an overestimation of about 35%. Figure 15 continued: Figure 7 in [START_REF] Haas | Comparing a quasi-3D to a full 3D nearshore circulation model: SHORECIRC and ROMS[END_REF]. Depth-integrated mean currents. The rip current slows down the waves, causing the waves to refract toward the centre of the channel (Figure 13), in agreement with the visual observation of the laboratory experiment reported by [START_REF] Haas | Quasi three dimensional modelling of rip current systems[END_REF], hereafter H03, (see their section 6). However, the rip current does not lead to a wave breaking more important than the one above the bar. The wave height actually decreases more slowly in the channel than above the bar, so that in the last meters before the coast the waves are higher in the channel axis than behind the bar (Figure 14). with the observation (Figure 14). As [START_REF] Haas | Laboratory measurements of the vertical structure of rip currents[END_REF] point out, this difference is associated with a longshore pressure gradient that drives the flow to the Laboratory experiments show that rip currents are unstable and are characterized by large temporal variations. The authors of these experiments emphasize the impossibility of reproducing twice the same series of measurements, even if the experimental conditions An exact agreement with the model is therefore not sought.

Nevertheless, we expect from the model that it gives the correct orders of magnitude of the and their variability. Figure 16 shows the temporal evolution of the current he vertical and over a wave period in the two channels. Given the above limitations, the agreement with the observations presented in H03 (see Figure 11.a in H03) The occurrence of the rip current is of the order of 200 seconds without being really regular.

The variability of the rip current is of the order of a few centimetres per second, with a maximum standard deviation around 6cm/s in the centre of the channel at y=13.5m and x=11.5m (Fig. 17). observation (Haas and Svendsen 2002, their figures 17, 18, 19). At the exit of the channel (offshore direction), the current decreases significantly.

larger and close to the bottom at x=

Conclusions

The vertical structure of the time-averaged cross-shore current in the rip channel (Fig. 18 the phase-averaged modelling of Haas and Warner intensity current around 0.2m/s between x=11m and x= point of maximum intensity (x=11.5m), the current is fairly homogeneous on the vertical except near the bottom where the current decreases strongly, in good agreement observation (Haas and Svendsen 2002, their figures 17, 18, 19). At the exit of the channel , the current decreases significantly. The vertical gradient also becomes The discrete version of the left-hand side of ( 30) is given by the left-side of:

( G )G ( ) ( ) = ( ) (9) + ( ) (9) + ⋯ C1

The right-hand side of C1 is given by a Taylor development, exponents in brackets indicating the order of the derivatives, ( ) (9) being the leading term of the errors made by the numerical scheme. Assuming a sinusoidal shape ( ( ) (9) ≈ _ (9), _ the wave frequency) leads to ( ) (9) ≈ _ (9). Similar reasoning can be done on the right-hand side of (30). The leading error term of the corresponding numerical scheme is , o , equivalent (if a sinusoidal shape is assumed, k the wavenumber) to o η and -o . Then reformulate (30) using the first two terms of the Taylor developments leads to (31).

APPENDIX D

The

Figure 1 .

 1 Figure 1. Dispersion relations at different bottom depth. Blue line: hydrostatic l

Figure 1

 1 Figure1shows these different dispersion (this choice will be discussed in the n

Figure 2 .

 2 Figure 2. Relative error on the wavelength computed by the ACM dispersion relation (24) at different

  Figure 3: r computed with (28) (Black line) and (29) (red line) as a function of N for h=50m and a

  compare it with that predicted by the theoretical relation dispersion of the incompressible linear theory (25). The modelled wavelength is diagnosed from the distance separating two successive identical values of the surface level in the direction of propagation. The grid points chosen for the calculation of this diagnosis are located beyond a distance of 625m from the open boundary and, moreover, the diagnosis is only calculated once the transient regime is reached (i.e. once the wave front has passed these points). A time average of the values obtained is finally calculated. A set of simulations is done to estimate the precision on (25) as a function of N (ranging from 2 to 20) and of the wave period. The smallest period considered here is 6 seconds, associated with a wavelength of 56m for h=50m, represented by 11 grid points given the horizontal resolution of 5m. This minimum value is considered as a reasonable threshold with respect to the capacity of the C-grid to represent the small wavelengths. The maximum period considered is 14 s which is associated with a wavelength of approximately 250 m. Periods greater than 14s are not considered because their wavelengths, much larger than the depths envisaged, render them slightly non-hydrostatic and therefore not interesting from the point of view of the authors.

Figure 4 :

 4 Figure 4: Relative error on the wavelength

Figure 5 :

 5 Figure 5: relative error on the wavelength

:

  relative error on the wavelength (vertical axis) as a function of N (horizontal axis) dashed black line indicates a relative error of -1.5%. The period is 6s. Each curve corresponds to a different number of vertical levels. From top to bottom: 80 regular levels regular levels (dashed bold black line), 20 regular levels (blue), 10 levels tightened regular levels (red).

  1.5%. The period is 6s. Each f vertical levels. From top to bottom: 80 regular levels 20 regular levels (blue), 10 levels tightened N=20), is attributable to hydrostatic aspects. This invoking the simple without term q) and = o (30) Given the numerical schemes used by our model, the error is respectively of the order of ≈ _ η and o ≈ o η for the left-and right-hand sides of (30). These errors are estimated from the first terms of a Taylor expansion (APPENDIX C). The compensation of temporal and spatial discretization errors is a characteristic of this type of discrete equation (Lemarié et al, 2015, page 127 paragraph 2). Here we see that the two above errors are compensated when the time step is at its theoretical maximum, i.e. 9 = r/o (with o = _/ ). When the time step is significantly smaller (consequence of the increase of N) the time-dependent error becomes negligible and the error is reduced to the error of the right-hand side, o , also equivalent to -o , and then the discretization of the (30) equation tends to behave like (APPENDIX C):

Figure 6 :

 6 Figure 6: Relative error on the wavelength as a function of r

(

  see the non-hydrostatic pressure gradient terms in elative error on the wavelength as a function of r. Model (a) and ACM linear the h=50m. N = 5. Horizontal black lines indicate a relative error of + or -1%. Each curve corresponds to a period: 14s (magenta), 12s (cyan), 10s (blue), 8s (green), 6s (red). The dashed vertical line (a) and r=0.9845 (b) for which all relative errors are between +/ 27 A modification of the pressure gradient scheme is now proposed to use moderate N values while maintaining an acceptable level of accuracy on dispersion properties. We around 5. The proposed modification is suggested by the observation ). This is equivalent phase speed is systematically lower than the theoretical phase hydrostatic effect excessively decreases the phase speed. hydrostatic component of in agreement with the conclusions drawn from the analytical solution of . The proposed improvement consists in ure gradient by multiplying the latter by a constant equations (7) and . Model (a) and ACM linear theory (b).

Figure 7 :

 7 Figure 7: as Figure 6a but h=20m.

Figure 8 :

 8 Figure 8: Standard deviation of the horizontal c

Figure 8b )

 8b Figure 8b) for r=0.976. The fact that a single value of

Figure 9 :

 9 Figure 9: Surface elevation (m)

Figure 10 .

 10 Figure 10. Top: Bottom depth (m). Middle: mean surface levels (m). Solid line: simulation. Squares:

Figure 11 .

 11 Figure 11. As Figure 10 but eddy viscosity calculated with (36) and (37). Top: mean surface level

Figure 11

 11 Figure11corresponds to an eddy viscosity calculated with (36) and (37) instead of (35) with

  Figure 12. Vertical profile of the measured mean without the term O 1 (o = 0, green lines), simulated with the term Vertical axis: normalized depth

  ame conditions would lead to a wavelength of about 1.9m, ie an 44 ame conditions would lead to a wavelength of about 1.9m, ie an The wave height and mean surface height calculated by the model are compared to s 6 and 16 of Haller et al, 2002). The wave n increases sharply when the waves ). The breaking conditions are quickly reached and once the ures 13 and 14). shore distance. Right: Average surface shore distance. In red: in the center of the basin (y = 9 m, with of the model is in good agreement with the description of the process given in Haller et al. (2002). In the channel ures 13 and 14, because shows the horizontal distribution of the averaged current, that is very similar to the circulation calculated with the phasepresented by Haas et al 2003 (their figure 5) and Haas and Warner

Figure 15 .

 15 Figure 15. Depth and time averaged current (m/s) computed by SYMPHONIE. Horizontal axis:

Figure 15

 15 Figure 15 continued: Figure 5 in[START_REF] Haas | Quasi three dimensional modelling of rip current systems[END_REF]. Time-averaged depth-integrated current from(a)

  The mean surface level increases at the coast, inducing a change in the average position of the water line allowed by the wetting-drying scheme of the model. The mean surface level at the coast is markedly different in the channel axis and in front of the bar, in agreement with the observation ( difference is associated with a longshore pressure gradient that drives the flow to the channels. Laboratory experiments show th by large temporal variations. The authors of these experiments emphasize the impossibility of reproducing twice the same series of measurements, even if the experimental conditions are a priori unchanged. An exact agreement with the model is therefore not sought.Nevertheless, we expect from the model that it gives the correct orders of magnitude of the rip currents and their variability. averaged over the vertical and over a wave period in the two channels. Given the above limitations, the agreement with the observations presented in H03 (see Figure11.a in H03) is good. The rip current at y=4.6m than the rip current at y=13.6m around 20 cm/s. The amplitude and time scale of its variability are also close to the behaviour described in H03.

Figure 16 .

 16 Figure 16. Cross-shore depth-averaged

  y=4.6m is a little overestimated but, as in H03, it is mostly low y=13.6m. The latter has the good order of magnitude, with maximum around 20 cm/s. The amplitude and time scale of its variability are also close to the averaged phase-averaged current (cm / s) in the two channels at x=10.8 m as a function of time (s). Periods shorter than 1s are filtered. Negative values in the direction of the wave generator. y = 13.6m (black line), y = 4.6m (red line).50). As[START_REF] Haas | Laboratory measurements of the vertical structure of rip currents[END_REF] point out, this difference is associated with a longshore pressure gradient that drives the flow to the at rip currents are unstable and are characterized by large temporal variations. The authors of these experiments emphasize the impossibility of reproducing twice the same series of measurements, even if the experimental conditions An exact agreement with the model is therefore not sought.Nevertheless, we expect from the model that it gives the correct orders of magnitude of the shows the temporal evolution of the current he vertical and over a wave period in the two channels. Given the above limitations, the agreement with the observations presented in H03 (see Figure11.a in H03) is a little overestimated but, as in H03, it is mostly lower . The latter has the good order of magnitude, with maximum around 20 cm/s. The amplitude and time scale of its variability are also close to the current (cm / s) in the two channels at x=10.8 m as a function of time (s). Periods shorter than 1s are filtered. Negative values in the

Figure 17 .

 17 Figure 17. Standard deviation of the . Standard deviation of the cross-shore depth-averaged current (m/s)

Fig. 18 )

 18 Fig.18), in good agreement with the observations (Haas and Svendsen 2002, their figures 13

Figure 18 .

 18 Figure 18. time-averaged cross

  er and close to the bottom at x=9m the current is even reversed , in good agreement with the observations[START_REF] Haas | Laboratory measurements of the vertical structure of rip currents[END_REF] figures 13 averaged cross-shore current (m/s) at y=13.6m. Negative values are in the offshore 53 ent in the rip channel (Fig.18) of Haas and Warner (2009, their Figure intensity current around 0.2m/s between x=11m and x=12m. At the 11.5m), the current is fairly homogeneous on the vertical except near the bottom where the current decreases strongly, in good agreement with the observation (Haas and Svendsen 2002, their figures 17, 18, 19). At the exit of the channel The vertical gradient also becomes (positive values in , in good agreement with the observations (Haas and Svendsen 2002, their figures 13 Negative values are in the offshore the non-hydrostatic behaviour depends on the fact that this adjustment is faster than the variation of the other fields. A meticulous adjustment of the time step of the model is in practice necessary, taking into account various parameters, including the wave propagation speed and a possible tuning consisting in reducing the non-hydrostatic pressure gradient. So far the reduction factor is a constant finding of a compromise in a set of possible values taking into account a range of periods and depths. Future development will adapt this tuning to the wave peak period and bathymetry using the analytical reduction factor solution that this study established from the linearized ACM equations.The model is tested on idealized simulations and laboratory experiments. The ability of the non-hydrostatic method to accurately represent the dispersion properties of surface waves is first demonstrated in the context of the linear theory.[START_REF] Hamilton | Establishing uniform longshore currents in a large-scale sediment transport facility[END_REF] laboratory experiment is used to evaluate the model on the issue of wave height amplification associated with beach slope as well as wave breaking. The parameterization of the wave-breaking is based on momentum horizontal diffusion, trigger conditions being related to a Froude number. The tuning of the horizontal viscosity is consistent with considerations on the mechanical energy balance. The reasoning can be extended to the possible production of turbulent kinetic energy in the turbulent closure scheme, leading in practice to the increase in vertical turbulent viscosity. If the vertical profile of the residual alongshore current can apparently be slightly improved, on the other hand, the fraction of the kinetic energy lost by the breaking wave, which can in this way be injected into the vertical turbulence, remains uncertain for the authors and requires, obviously, to be better understood. Along these lines, a part of the wave spectrum will still be missing in our future realistic applications because of the limited grid resolution (the smallest periods being concerned) with, as a result, a possible underestimation of the bottom stress in shallow areas that should be taken into account.Finally, the[START_REF] Haller | Experimental study of nearshore dynamics on a barred beach with rip channels[END_REF] laboratory experiment, with its bar and its two channels, makes it possible to evaluate the model on the rip current problem. The model correctly reproduces rips currents and their high temporal variability. The feedback of the rip

  potential energy flux is " = with = /o. A sinusoidal form sinis assumed leading to " = # 8 -˜b Ψ . A wave period averaging gives " = # 8 . Within the wave breaking zone we assume ~ = }. Using o = < ℎ leads to "

Table 1 :

 1 results not shown). Thickness of the ten layers (1 corresponds to the bottom) as a fraction of the water depth.

	levels	Relative thickness	⁄	ℎ
	10	0.03852
	9	0.04623
	8	0.05547
	7	0.06657
	6	0.07988

Acknowledgements

The authors acknowledge the support of the SIROCCO team http://sirocco.obs-mip.fr/, the NUMEROFIX project (NUmerical ModElling of Real Ocean surFace mIXing) funded by LEFE/GMMC, the HPC resources from CALMIP (Grant 2018-P1325) and from GENCI and CINES (Grand Equipement National de Calcul Intensif, project A0040110088).

Short black line with label: scale of the wavelength at the entrance of the basin. Horizontal axis: cross-shore distance (m) from the wave generator. Vertical axis: alongshore distance (m) Consistently with [START_REF] Haas | Laboratory measurements of the vertical structure of rip currents[END_REF] and [START_REF] Haas | Quasi three dimensional modelling of rip current systems[END_REF], the simulation lasts 30 minutes to highlight the observed temporal variations of rip currents. The wave height is obtained from the difference between the minimum and maximum of the surface level obtained over time intervals of 1 second. An average wave height is then calculated by averaging these values obtained every second over the simulation except the first minute to avoid the spin-up phase. Table 2 summarizes the characteristics of the numerical grid and the calculation time. The latter can be roughly extrapolated to more realistic applications considering that the scale ratio between this simulation and that of a realistic case (with for example dx=3m) is about 50, in other words an area of about 1km x 1km and a simulation duration of about 1 day. The A second breaking zone is highlighted ne agreement with the description of the process given in [START_REF] Haller | Experimental study of nearshore dynamics on a barred beach with rip channels[END_REF]. In the channel axis, the wave height increases before the channel, as shown in of the blocking by the opposite rip current depth-averaged current, that is very similar to the circulation calculated with the phase averaged modelling presented by Haas et al 2003 (their figure 5) and Haas and Warner 2009 (their figure 7).

The wave height and mean surface height calculated by the model are compared to the observations of the TEST B case (Figures 6 and 16 of Haller et al, 2002). The wave height first increases slowly towards the coast, then increases sharply when the waves Figure 14). The breaking conditions are quickly reached and once the bar crossed, the wave height is much lower and relatively constant (Figures 13 and 14 A second breaking zone is highlighted near the coast. This behaviour of the model is in good agreement with the description of the process given in [START_REF] Haller | Experimental study of nearshore dynamics on a barred beach with rip channels[END_REF]. In the channel axis, the wave height increases before the channel, as shown in Figures 13 and 14 With the exception of the post-breaking zone that can be processed in some models by a hydrostatic approach, the phase-resolved wave modelling requires a non-hydrostatic model. As the representation of the non-hydrostatic problem in all its complexity seems to be inevitably very expensive to date [START_REF] Klingbeil | The numerics of hydrostatic structured-grid coastal ocean models: state of the art and future perspectives[END_REF], a scheme allowing a reasonable compromise between, on the one hand the extra cost of the computation, and on the other hand an acceptable loss of accuracy of non-hydrostatic properties of waves, has been introduced into the circulation model. It is shown that relative errors of less than a few % (on the dispersion relation and the vertical current profile) can be obtained with non-hydrostatic simulations costing 5 to 6 times more than a hydrostatic model calculated under the same conditions. The non-hydrostatic additional cost remains obviously high but low enough to consider simulations of realist events such as storms on very large areas and it is anyway worth noting that wave-current effect modelling based on a hydrostatic circulation model normally couple the latter with a phase-averaged wave model whose cost is possibly not negligible.

The non-hydrostatic method implemented in the circulation model is derived from the method proposed by [START_REF] Lee | Application of the artificial compressibility method for turbulent open channel flows[END_REF]. In practice, the classical Poisson problem is transformed into a pressure propagation equation. The resolution is iterative, local, and does not need to solve a system of linear equations. An essential aspect of this method is to properly control the various parameters related to the iterative aspect of the problem. The equation of pressure propagation depends on a constant (the latter being comparable to a propagation velocity) which determines the speed of the adjustment of the non-hydrostatic pressure to a change of the hydrostatic terms in the momentum equations. The precision of currents on the waves, in terms of wave amplification and refraction, is well represented. The onset of wave breaking is also well simulated, either above the bar or in the channel axis.

This paper does not consider some 3D aspects as the effect of tracers on the pressure gradient and turbulence, which will be taken into account in future applications of this model. This also includes sediments that may be of great importance for modelling of rip currents as shown by [START_REF] Ma | A wave-resolving model for nearshore suspended sediment transport[END_REF]. These authors show that a wave-resolving non hydrostatic model is an efficient approach for the problem of resuspension of sediments by waves, and the possible feedback introduced by the modification of the equation of state by sediments, which can have a significant influence on the pressure gradient and the level of turbulence. In practice our wave model will be coupled with the MUSTANG sediment model (Le [START_REF] Hir | Dynamics of sand and mud mixtures: A multiprocess-based modelling strategy[END_REF] The next step is now to implement our model on realistic beaches, where we will confront the results with in situ data. Infragravity waves and the risk of submersion are possible applications. Some developments are still needed to provide a realistic wave spectrum to the boundary conditions of the phase-resolved wave model. The phaseaveraged model, which can be applied over large domains at a reasonable cost, appears to be a good candidate for providing the boundary conditions to the phase-resolved model. The phase-averaged wave model can also provide an estimate of the missing part of the wave spectrum in the deterministic model, helping the former to improve bottom turbulence in shallow areas. Finally, wave growth by local wind action will eventually be considered in the perspective of large areas [START_REF] Deigaard | Wind generation of waves: Energy and momentum transfer -An overview with physical discussion[END_REF].

APPENDIX A

The discrete version of at the left-hand side of ( 4) is given by the left-side of:

The right-hand side of A1 is given by a Taylor development, being the leading term of the errors made by the numerical scheme. Assuming a sinusoidal shape ( ≈ , k the wave number) leads to (5).

APPENDIX B: time stepping scheme

The equations are calculated with a forward-backward time stepping scheme. The alternation of the fields along the time axis is presented in Figure 19. The height of the surface at time t, , is calculated first, knowing and velocities at time t-½. Equation ( 9) gives knowing , , and hydrostatic tendencies at time t-1. Note that the vertical laplacian of q is computed with an implicit centric scheme of the type '' + ' so as not to limit the time step when the vertical mesh becomes very small (which becomes the case in sigma coordinates when the bathymetry tends to zero). The vertical velocity at time t-½, used for the advection terms of the equations for u and v, is computed from u and v at time t-½. Finally, the velocities at time t + ½ are calculated from the velocities, the sea surface height and the non hydrostatic pressure of the previous times. Vertical turbulent diffusion terms are calculated with an implicit scheme.

APPENDIX C

Consider a constant diffusion coefficient < 0 1 >. The dissipation term due to horizontal diffusion is € 1 = -< 0 1 > . Assuming a constant wavenumber leads to = 8 o›-Ψ

and € 1 = -< 0 1 > 8 o›-Ψ. A wave period averaging gives

. Using ~ = } leads to: