Anne-Claire Bennis

Fabrice Ardhuin

Franck Dumas

On the coupling of wave and three-dimensional circulation models : Choice of theoretical framework, practical implementation and adiabatic tests

Keywords: wave-current coupling, radiation stresses, MARS3D, WAVEWATCH III

Many theoretical approaches and implementations have been proposed for the coupling of the three-dimensional ocean circulation with waves. The theoretical models are reviewed and it is shown that the formulation in terms of the quasi-Eulerian velocity circumvents the essential difficulty of alternative formulations for the Lagrangian mean velocity. Namely, models based on this Lagrangian velocity require an estimation of wave-induced motions to first order in the horizontal gradients of the wave field in order to estimate the vertical flux of wave pseudo-momentum. So far, only three-dimensional wave models have been able to provide these estimates, and all published theories based on the simpler Airy theory are not consistent at the leading order, because they ignore or incorrectly estimate the vertical momentum flux. With an adiabatic example on a sloping bottom it is shown that this inconsistency produces very large spurious velocities. These errors are independent of the slope for the inviscid case, and are still significant when a realistic vertical mixing is applied. A quick diagnostic of the potential accuracy of a theoretical model is the vertical profile of the wave-induced forcing terms: if it is not uniform over depth in adiabatic conditions then it will produce spurious artificial flow patterns in conditions with shoaling waves. Although conceptually more challenging, the quasi-Eulerian velocity theories only introduce minor modifications of the solution procedure for the standard primitive equations: a modification of the surface boundary condition for the mass conservation, the addition of the Stokes drift in the tracer advection equations, and sources of momentum and turbulent kinetic energy with associated surface and bottom fluxes. All the necessary modifications of primitive equation models are

Introduction

Three-dimensional oceanic flows can be strongly forced or modified by waves, in particular in the nearshore (e.g. [START_REF] Longuet-Higgins | Longshore currents generated by obliquely incident sea waves, 1[END_REF]Newberger and Allen, 2007a) and the coastal ocean [START_REF] Lentz | Observations and a model of undertow over the inner continental shelf[END_REF]). Yet the numerical modelling of these complex flows is only slowly coming of age, with recent works using quasi-three dimensional [START_REF] Haas | Quasi-threedimensional modeling of rip current systems[END_REF] or fully three dimensional models [START_REF] Uchiyama | Wave-current interaction in nearshore shear instability analyzed with a vortex force formalism[END_REF]. Although models capable of resolving the wave motion are becoming feasible on small scales, as shown by [START_REF] Lubin | Threedimensional large eddy simulation of air entrainment under plunging breaking waves[END_REF], the demands of coastal zone management in terms of coverage and resolution are still barely met by depth-integrated models in which the wave motion is averaged over the phase of at least the short waves (e.g. [START_REF] Reniers | Morphodynamic modeling of an embayed beach under wave group forcing[END_REF]). There has thus been a large effort to develop models of intermediate complexity, capable of resolving the vertical structure of the mean flow which may be needed to account for mixing and dispersion [START_REF] Svendsen | Nearshore mixing and dispersion[END_REF] while still keeping the hydrostatic approximation for the mean flow. Similarly, a proper representation of near-surface currents and drift requires the introduction of wave effects, in particular the Stokes drift and wave-induced mixing (e.g. [START_REF] Rascle | Drift and mixing under the ocean surface revisited. stratified conditions and model-data comparisons[END_REF]. For these applications, the wide community of users of numerical models for the ocean circulation such as POM and ROMS [START_REF] Blumberg | A description of a three-dimensional coastal ocean model[END_REF][START_REF] Shchepetkin | A method for computing horizontal pressure-gradient force in an oceanic model with nonaligned vertical coordinate[END_REF] is calling for minor modifications to make these models capable of representing waves effects.

A large body of often conflicting theoretical results have been published on the form of the wave-modified primitive equations that would be suitable for such models. We may cite, in chronological order, [START_REF] Dolata | Wave setup and wave-induced currents in coastal zones[END_REF], [START_REF] Jenkins | The use of a wave prediction model for driving a nearsurface current model[END_REF], [START_REF] Weber | Transient ocean currents induced by wind and growing waves[END_REF], [START_REF] Rivero | On the vertical distribution of u w[END_REF], [START_REF] Péchon | Numerical modelling of the three-dimensional wave-driven currents in the surf zone[END_REF], [START_REF] Groeneweg | Changes in the mean velocity profiles in the combined wave-current motion described in GLM formulation[END_REF], [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF], [START_REF] Perrie | The impact of waves on surface currents[END_REF], [START_REF] Mcwilliams | An asymptotic theory for the interaction of waves and currents in coastal waters[END_REF], [START_REF] Xia | Vertical variation in radiation stress and wave-induced current[END_REF], Newberger and Allen (2007b), Ardhuin et al. (2008b), Ardhuin et al. (2008a), Mellor (2008) ... Although each of these work is based on a particular set of hypotheses, for example some assume a horizontally uniform wave field, it is expected that all theories should agree on the most simple cases that they are supposed to cover. In practice this is not the case because some of these theoretical models are not self-consistent.

Our purpose is not to blame this or that author for inconsistencies, limitations, or mistakes. The goal of the present paper is rather to help people implement correctly the effect ot waves in a primitive equation model. Here we explain why some of the wave-averaged equations may appear different but still represent correctly the same reality, and we give some constraints that should be obeyed by wave-averaged equations so that model developers can make their own judgement of the published equations, even before implementing them. One basic test of the wave-induced forcing should be that it is uniform over the vertical in adiabatic conditions with a stationary wave forcing, for example for waves shoaling on a slope. Recent published equations from Mellor (2003, hereinafter M03) and Mellor (2011a) to Mellor (2008, hereinafter M08) and Mellor (2011b) and [START_REF] Xia | Vertical variation in radiation stress and wave-induced current[END_REF] do not verify this condition, as illustrated by Figure 1. This inappropriate forcing profile is due to a minor inconsistency in the M03 derivation or more fundamental problems in M08. This problem was first pointed out by Ardhuin et al. (2008a, hereinafter AJB08), and its practical consequences are shown for the first time. Although errors in the underlying theoretical model are likely to be dwarfed by parametrization errors in the case of strongly dissipative environments like the surf zone, they may still explain some of the differences found between various models (e.g. [START_REF] Haas | Comparing a quasi-3d to a full 3d nearshore circulation model: SHORECIRC and ROMS[END_REF].

We also take the present opportunity to present the approximated Generalized Lagrangian Mean equations (glm2 -z) by Ardhuin et al. (2008b) (hereinafter ARB08) in a more readable form, giving details on how they were implemented in the MARS3D flow model [START_REF] Lazure | An external-internal mode coupling for a 3d hydrodynamical model for applications at regional scale (MARS)[END_REF]. These steps have already been pioneered by [START_REF] Mcwilliams | An asymptotic theory for the interaction of waves and currents in coastal waters[END_REF] and [START_REF] Uchiyama | Wave-current interaction in nearshore shear instability analyzed with a vortex force formalism[END_REF] with equations that are mathematically consistent with those in Ardhuin et al. (2008b). These more theoretical presentations are often obscured by their generality and completeness. We thus here present the equations and implementation in the most simple form, warning the reader when the simplification causes a loss of generality. The present paper is thus an introduction to Ardhuin et al. (2008b) and[START_REF] Mcwilliams | An asymptotic theory for the interaction of waves and currents in coastal waters[END_REF] oriented towards practical implementations. For a physical discussion one may read [START_REF] Lane | Wave-current interaction: A comparison of radiation-stress and vortex-force representations[END_REF]. We focus here on non-dissipative conditions where exact solutions are most easily found, and we refer to [START_REF] Uchiyama | Wave-current interaction in an oceanic circulation model with a vortex-force formalism: Application to the surf zone[END_REF] for further discussion of wave breaking, mixing and bottom friction parametrizations.

!! !"#$ " "#$ ! !#$ % %#$ &'!" !( !! !"#$ " )*+,!-./

Theoretical analysis of wave-averaged equations for the Lagrangian velocity

A brief review

Theories can be categorized according to two criteria (Figure 2). First, the equations are depth-integrated or not. For the depth-integrated equations, with some approximations related to wave non-linearity, the problem is clear. Second, the momentum balance can take two forms, which are equivalent [START_REF] Longuet-Higgins | Radiation stress in water waves, a physical discussion with applications[END_REF][START_REF] Garrett | Generation of Langmuir circulations by surface waves -a feedback mechanism[END_REF][START_REF] Smith | Wave-current interactions in finite-depth[END_REF]. One form of the equations is for the total momentum in the water column M ( M in Phillips 1977, page 62), and the other is for the momentum of the mean flow only M m ( M in Phillips 1977, page 61). M, is related to the mass transport velocity, which naturally arises when working with fluxes of wellmixed solutes. For a mean current u that is uniform below the level of wave troughs, the other momentum variable is approximately M m = ρuD where D is the mean water depth and ρ is a depth-averaged water density. The wave momentum1 is simply the difference M w = M -M m , and it is a horizontal vector. Hence, M m may be more closely related to the mean current measurable by a fixed instrument. This interpretation, however, has to be considered with caution, since there may be different ways to extend the definition of u from the wave troughs to the wave crest level, which is needed for a model that includes surf zones and adjacent shallow areas.

These details cannot be ignored when considering the three-dimensional sets of equations, for the total momentum ρU, where U is the Lagrangian velocity or the mean flow momentum ρ u, where u is the quasi-Eulerian velocity, see eq. ( 2). [START_REF] Andrews | An exact theory of nonlinear waves on a Lagrangian-mean flow[END_REF] have derived exact equations for U and u that are equivalent.

Thus, at the same order of approximation, we only have one set of equation for U and another for u. Both sets of equations are mathematically equivalent. Any other equation must have some internal inconsistencies. We shall illustrate this statement for M03, because these are the most widely used equations, and we give enough details on other sets of equations, such as [START_REF] Xia | Vertical variation in radiation stress and wave-induced current[END_REF], so as to make the problem obvious.

A generic recipe for three-dimensional theories

Any theory for wave-averaged equations goes through three steps. First, the control volume in which the momentum is averaged must be defined. For depth-integrated equations, this is simply the full water column. For three-dimensional equations, there is a vertical discretization. The control volume may be moving, following all or part of the wave motion (Figure 3). The mean position x = (x, y, z) of the volume, is associated with the actual position ξ(x, y, z, t) = (ξ 1 (x, y, z, t), ξ 2 (x, y, z, t), ξ 3 (x, y, z, t)) such as

x + ξ(x, y, z, t) = (x + ξ 1 (x, y, z, t), y + ξ 2 (x, y, z, t), z + ξ 3 (x, y, z, t)) . (1)
The boundaries of the volume at time t 0 are the sides of an elementary cube centered of x and deformed by the transformation x → x + ξ. 2004) and Ardhuin et al. (2008b), the second derivative in the vertical current profile is neglected in the second, while the first is only adiabatic. Names that appear in red correspond to theories that are not fully consistent with their originating hypotheses. In the case of 3D theories for the total momentum, the problem generally comes from the vertical flux of momentum and may only arise on a sloping bottom, not explicitly considered by [START_REF] Groeneweg | Wave-current interactions in a generalized Lagrangian mean formulation[END_REF].

For the Eulerian average, the volume does not move, it keeps its original cubic shape and orientation, and ξ = (ξ 1 , ξ 2 , ξ 3 ) = 0. For the average defined by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF], the horizontal displacements are zero, ξ 1 = 0 and ξ 2 = 0, and ξ 3 follows the wave motion. This is the simplest possible form of a wave-following coordinate, and it thus has a great appeal.

Second, the momentum balance of the volume is the result of body forces and momentum fluxes through the boundary. These fluxes involves an advective part (ρu2 for the horizontal advection of horizontal momentum and ρuw for the vertical advection) and a stress part that is the product of the stress tensor with the unit vector normal to the boundary. Neglecting shear stresses for simplicity leaves only the pressure p for the stress. In general, the sides of the control volume are sloping, which is very important for the exact definition of the momentum fluxes: the momentum flux associated to pressure forces is not a diagonal tensor in general coordinates. For example, if and only if a surface is not horizontal or vertical, there is a pressure-induced vertical flux of the horizontal momentum through it.

For ξ = 0, the control volumes have tilted facets, with a flux of horizontal momentum through the sloping bottom and top. The corresponding flux of the x-component of the momentum is p∂ξ 3 /∂x. Fluxes due to sloping iso-coordinates are often forgotten (e.g. [START_REF] Xia | Vertical variation in radiation stress and wave-induced current[END_REF], Mellor 2008b 2 ) or poorly approximated (e.g. [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF].

Last, the averaged equations can be transformed to another coordinate system, such as terrain-following coordinates. Averaging can cause an implicit distortion of the coordinates, for example the mean flow may appear divergent although the original flow is not. This effect may be easily corrected (see e.g. Ardhuin et al., 2008b, for the transformation of approximated Generalized Lagrangian Mean equations to Cartesian and terrain-following coordinates).

Because of the large difference in density between air and water, the use of a standard Eulerian average, used by [START_REF] Rivero | On the vertical distribution of u w[END_REF] or Newberger and Allen (2007b), is problematic in the region between crests and troughs, where both air and water are to be found. A strict Eulerian average produces a continuously varying density ρ from about 1.29 kg.m -3 to 1026 kg.m -3 , which is not compatible with the usual primitive equation models as it would introduce a strong diffusion of properties, including heat content, between the atmosphere and ocean. Mathematical extension of the velocity field across the interface have been used by [START_REF] Mcwilliams | An asymptotic theory for the interaction of waves and currents in coastal waters[END_REF], but it provides quantities that have, a priori, no simple physical interpretation since they are not given by a known averaging operator. Yet, Ardhuin et al. (2008b) showed that the resulting velocity actually corresponds to the quasi-Eulerian velocity ( u, v, w) first introduced by [START_REF] Jenkins | The use of a wave prediction model for driving a nearsurface current model[END_REF]: this is the mean Lagrangian velocity (U, V, W ) minus the wave-induced drift (U s , V s , W s ),

( u, v, w) = (U, V, W ) -(U s , V s , W s ).
(2)

This definition requires a wave-following coordinate system. The averaging is also connected to the choice of the momentum variable, which varies between different theories, as summarized in Figure 2.

-p dξ 1 /dx

ψ=270°ψ =0 ψ=180°ψ =90°- p dξ 3 /dx u(w-dξ 3 /dt) u(u-dξ 1 /dt)
Figure 3: Illustration of the momentum fluxes into a moving control volume defined by the (possibly moving) position vector of the control volume, Ξ(x, y, z, t) = x+ξ (where x is the mean position and ξ is the volume control displacement). Here the control volume moves in a circle and is shown for 4 wave phase, ψ = 0 • , ψ = 90 • , ψ = 180 • and ψ = 270 • . u is the full flow velocity (including waves and current). The fluxes are the sum of advective fluxes through the facets and the pressure acting on the possibly sloping facets. For the advective part, the velocity is relative to the moving facet, which produces terms like -u∂ξ 3 /∂t.

Momentum equations

In the previous section, we defined the control volume and we presented the different averages existing in the literature . We write now the generic equation for the mean momentum (e.g. mean Lagrangian velocity or quasi-Eulerian velocity) obtained after averaging on the control volume. The mean can be Eulerian, GLM or following ξ 3 . Neglecting buoyancy, Earth rotation, mixing and viscous effect, the generic conservation equation for the mean momentum ρu (where u = (u, w) and (•) is the averaging operator) in the control volume and with motion restricted to a vertical plane is

∂(ρu) ∂t + ∂(ρu u) ∂x + ∂(ρu w) ∂z = F, (3)
with F a forcing term due to waves and water levels which explicitly given below.

The rate of change of the momentum (the first term on the left hand side of eq. ( 3) equals body forces (such as gravity, absent for the horizontal momentum ρu), plus the divergence of fluxes of momentum. The mean flow advective fluxes appear in the left hand side of (3), the other fluxes constitute the forcing term F . For a generic control volume, F is

F = F uu + F uw + F px + F p3 , (4)
(see [START_REF] Brekhovskikh | Mechanics of Continua and Waves Dynamics[END_REF] for more details)

with the horizontal advection

F uu = ∂(ρ u 2 ) ∂x . (5)
u is the horizontal velocity associated to the wave propagation such as

u = kac [exp(kz) cos(kx -σt)] . ( 6
)
where k is the wavenumber, a is the wave amplitude, c is the phase speed, σ is the frequency.

The vertical advection (where w and ξ3 are, respectively, the vertical velocity w = kac [exp(kz) sin(kx -σt)] and the vertical position associated to the wave propagation)

F uw = ∂ w -∂ ξ 3 /∂t ρ u ∂(z + ξ 3 ) , (7)
is important for fixed (Eulerian) control volumes (e.g. [START_REF] Rivero | On the vertical distribution of u w[END_REF] but is negligible for the M03-AJB08 volume, designed to make w -∂ ξ 3 /∂t as small as possible.

The last two terms are the pressure gradient across the sides of the volume, here assumed vertical3 ,

F px = ∂p ∂x , (8)
and the vertical divergence of the pressure-induced flux through the sloping iso-surfaces of the vertical coordinate,

F p3 = ∂S p3 ∂z = ∂ ∂z p(∂ ξ 3 /∂x) , (9)
which is zero in Eulerian averages, for which ξ 3 = 0 (where ξ 3 and p are respectively the vertical displacement and pressure induced by waves). We recall that p can be different for different control volumes, i.e. in the case of M03-AJB08 the pressure is nearly hydrostatic, which is not the case of a fixed control volume for which a -ρ w 2 correction occurs. This derivation was correctly done by M03 for his control volume. Using the notations of that paper, s = ξ 3 = ξ 3 , ξ 1 = 0, ξ 2 = 0. On the contrary, [START_REF] Xia | Vertical variation in radiation stress and wave-induced current[END_REF] completely omitted the vertical flux. In the case of Mellor (2008) the average of the vertical flux gradient F p3 was not estimated and replaced by an Eulerian average of the pressure gradient, which is not consistent with the Lagrangian average of the other terms, as described by [START_REF] Bennis | Comments on the depth-dependent current and wave interaction equations: A revision[END_REF]. We will thus not discuss further that paper.

Estimation of p(∂ ξ 3 /∂x):

The fundamental problem of the Lagrangian approach This vertical flux term redistributes wave momentum over the vertical. Thus, if that term is omitted, as in [START_REF] Xia | Vertical variation in radiation stress and wave-induced current[END_REF], the depth-integrated equations can still be correct. However, omitting this term, or giving an inconsistent approximation for it, generally produces a non-uniform profile of the net wave-induced force F px + F p3 in the absence of wave breaking (see Figure 1). In the adiabatic limit, without mixing or friction, there is no possible balance of this force with the pressure gradient induced by the sea surface slope, which is not a physical situation. Here is a mathematical description of the problem.

The equation of motion are given by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] (his equation ( 34a)). When we neglected the Coriolis force, density stratification, mixing and after applying the continuity equation, the equations correspond to equation (3) for the particular case of the control volume that follows only the vertical wave-induced motions and with motion restricted to a vertical plane. The equation is similar to the following equation and rewritten in a terrain following vertical coordinate, ς, it is

∂U ∂t + U ∂U ∂x + W D ∂U ∂ς = F D . ( 10
)
where (U, W ) is the Lagrangian-mean flow field.

The force on the right hand side is given by Equation ( 4) with F uw = 0,

F px = -gD ∂ η ∂x
(where η is the mean surface elevation), F uu and F p3 are given in Appendix D. Over a sloping bottom in finite water depth, F uu is of the order of the bottom slope ε 2 , and a consistent solution of eq. ( 3) requires an approximate of F p3 to the same order. Because this is the product of p and s, each of these two terms must be estimated to first order in ε 2 . Airy theory, which is the solution for waves over a flat bottom and is used by M03 and almost all parameterizations, is thus insufficient. By analogy with the problem of estimating the reflection of waves over a sloping bottom [START_REF] Meyer | Surface wave reflection by underwater ridges[END_REF], we conjecture that there is no general analytical expression for the vertical flux F p3 in terms of local wave, water depth and current properties, even in the limit of small bottom slopes.

So far we have no proof for this conjecture, but it is supported by several facts. First of all, all published analytical expressions are based on Airy theory and either ignore the vertical flux F p3 , this is the case of [START_REF] Walstra | Calculation of wavedriven currents in a 3D mean flow model[END_REF] and [START_REF] Xia | Vertical variation in radiation stress and wave-induced current[END_REF], or give erroneous expressions for it, as in [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF], Mellor (2008), Mellor (2011a), Mellor (2011b). Their error is easily seen by comparing, for our adiabatic test case with shoaling waves, the net force F , including the pressure gradient due to the surface slope. All of these theories give a strong vertical shear for F as in figure 1, whereas it should be uniformly zero, allowing a stationary solution. We observe that only Ardhuin et al. (2008b) proposed a correct estimation of F p3 , but they used the now classical approach of approximating the three-dimensional wave problem as a series of modes with different vertical structures and amplitudes that vary on the horizontal. For irrotational flows, the exact solution to the Laplace equation and boundary conditions is recovered in the limit of an infinite number of modes. The first mode corresponds to the Airy wave solution while the other evanescent mode are not propagating modes, but they correspond to the near-local adjustment of the wave field to the topography. Contrary to the wave action equation which is hyperbolic, the equations for the amplitudes of the evanescent modes are a set of coupled elliptic equations. As a result, these modes are not functions of the local water depth only but they are generally defined by the water depths over the entire domain. Very often the adjustment is only significant over one wavelength and for steep slopes, as in the simplest problem of wave propagation over a step (e.g. [START_REF] Rey | Propagation of surface gravity waves over a rectangular submerged bar[END_REF][START_REF] Rhee | On the transmission of water waves over a shelf[END_REF]. In that case, the adaptation of the iso-potential lines clearly extends over a finite distance from the step showing that the solution is indeed not only a function of the local (flat) bottom properties and amplitude of the freely propagating waves. The same is true, although less conspicuous, for a smooth bottom (e.g. [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small amplitude water waves over variable bathymetry regions[END_REF].

Obtaining the solution of the coupled equations is a formidable computing task (e.g. [START_REF] Magne | Evolution of surface gravity waves over a submarine canyon[END_REF], even if the series of modes is made to converge faster, as is the case with the vertical mode decomposition proposed by [START_REF] Chandrasekera | Extended linear refractiondiffraction model[END_REF]. What looks like a mathematical curiosity has dramatic consequences in terms of wave-driven flows, with very large spurious errors for poorly approximated profiles of F p3 , as we shall see in section 3. Given this impasse, the only option available today is thus to work with the quasi-Eulerian velocity.

Equations for the quasi-Eulerian velocity

In order to facilitate the implementation of quasi-Eulerian equations, we give here a short and simplified account of Ardhuin et al. (2008b). Starting from the general equations of [START_REF] Andrews | An exact theory of nonlinear waves on a Lagrangian-mean flow[END_REF], Ardhuin et al. (2008b) have given an approximation to second order in the wave non-linearity and transformed the equations with a change of the vertical coordinate, so that the Jacobian associated with the averaging procedure is equal to one, and both the resulting quasi-Eulerian flow field ( u, v, w) and Lagrangian-mean flow field (U, V, W ) are non-divergent.

The quasi-Eulerian flow field does not contain wave momentum, and this is the main difference with the M03 theory. Solving for ( u, v, w) removes the problem of the approximation of the vertical fluxes of wave momentum in S p3 (equation ( 9)) because the tricky part is a vertical flux of wave momentum: without wave momentum there is no problem anymore in the vertical fluxes of the quasi-Eulerian momentum. The influence of waves on the mean quasi-Eulerian current appear as forcing terms (A1, A2, B1, B2, C1, C2, C3 in the equations ( 11),( 12),( 15)).

In order to simplify the equations we generally give the wave forcing expressions for monochromatic waves as a function of the surface elevation variance E. In the case of quasi-linear random waves the corresponding forcing is simply the sum of the monochromatic wave forcing with E replaced by the elementary elevation variance E(f, θ)df dθ (where f and θ are, respectively, the wave frequency and direction), as detailed in Appendix C. The following GLM equations use the Cartesian z coordinate. However, the most of coastal hydrodynamical models (e.g. [START_REF] Marsaleix | Energy conservation issues in sigma-coordinate free-surface ocean models[END_REF] use equations in terrain-following coordinates (Appendix A), further transformed in flux form (Appendix B). This last form is most compact because the vertical stokes drift W s disappear from it.

Momentum, mass, and tracer conservation

For simplicity we neglect the effect of the vertical current shear and partial standing waves in the wave forcing term, so that eq. ( 42) in Ardhuin et al. (2008b

) becomes 4 ∂ u ∂t + u ∂ u ∂x + v ∂ u ∂y + w ∂ u ∂z -f v + 1 ρ ∂p H ∂x = f + ∂ v ∂x - ∂ u ∂y V s -W s ∂ u ∂z - ∂J ∂x A1 + F m,x + F d,x + F b,x A2 , (11)
and

∂ v ∂t + u ∂ v ∂x + v ∂ v ∂y + w ∂ v ∂z + f u + 1 ρ ∂p H ∂y = -f + ∂ v ∂x - ∂ u ∂y U s -W s ∂ v ∂z - ∂J ∂y B1 + F m,y + F d,y + F b,y B2 ,( 12
)
where the left hand side is the classical primitive equation model for the quasi-eulerian velocity ( u, v, w) with p H the hydrostatic pressure, ( F m,x , F m,y ) the mixing effects (that redistribute momentum). The right hand side contains the forcing terms (see A1, B1, A2, B2 terms, with only the A1 and B1 terms for adiabatic case) where ( F d,x , F d,y ) the source of quasi-Eulerian momentum that is equal to the sink of wave momentum due to breaking and wave-turbulence interaction, ( F b,x , F b,y ) the source of quasi-Eulerian momentum that is equal to the sink of wave momentum due to bottom friction, which should only be included when the wave bottom boundary layer is resolved, J the wave-induced mean pressure (eq. ( 20)), and (U s , V s , W s ) the three-dimensional Stokes drift5 . Mixing is also influenced by waves, but this aspect will not be discussed here (see [START_REF] Craig | Modeling wave-enhanced turbulence in the ocean surface layer[END_REF][START_REF] Groeneweg | Changes in the mean velocity profiles in the combined wave-current motion described in GLM formulation[END_REF][START_REF] Rascle | Drift and mixing under the ocean surface revisited. stratified conditions and model-data comparisons[END_REF][START_REF] Uchiyama | Wave-current interaction in an oceanic circulation model with a vortex-force formalism: Application to the surf zone[END_REF]. The second lines in eq. ( 11)-( 12) contain the vortex force introduced by [START_REF] Garrett | Generation of Langmuir circulations by surface waves -a feedback mechanism[END_REF] in this context, and further discussed by [START_REF] Lane | Wave-current interaction: A comparison of radiation-stress and vortex-force representations[END_REF] and [START_REF] Smith | Wave-current interactions in finite-depth[END_REF].

The mass conservation is

∂ u ∂x + ∂ v ∂y + ∂ w ∂z = 0, ( 13
)
and the evolution of a conservative passive tracer concentration C is,

∂C ∂t + ∂ ∂x   ( u + U s )C C1   + ∂ ∂y   ( v + V s )C C2   + ∂ ∂z   ( w + W s )C C3   (14) = ∂ ∂x γ x ∂C ∂x + ∂ ∂y γ y ∂C ∂y + ∂ ∂z γ z ∂C ∂z . ( 15
)
where γ = (γ x , γ y , γ z ) is the turbulent diffusivity.

All four conservation equations are valid from the bottom z = -h to the local phase-averaged free surface η.

However, with the mode splitting, there is another important modification that is made through the barotropic mode. The surface kinematic boundary condition is given by

∂ η ∂t + ( u + U s ) ∂ η ∂x + ( v + V s ) ∂ η ∂y = w + W s . ( 16
)
It thus appears that, as in [START_REF] Mcwilliams | An asymptotic theory for the interaction of waves and currents in coastal waters[END_REF] or Newberger and Allen (2007b), there is a source of mass at the surface that compensates the convergence of the Stokes drift. In surface-following coordinates there is no velocity through the surface and w + W s vanish, leaving only the convergence of the Stokes drift to force the usual mass conservation equation (see Appendix A: eq. (A.3)).

For tracers, because the equations are unchanged (only for the explicit presence of the Lagrangian mean velocity), the boundary conditions are unchanged from classical primitive equation models.

Wave-induced forcing terms

The three-component Stokes drift (U s , V s , W s ), wave-induced pressure term J, and momentum source due to wave dissipation ( F d,x , F d,y ), can all be computed from only a few local parameters. These include the wave-induced surface elevation variance E, the phase-averaged water depth D = h + η, the wavenumber vector k = k(cos θ, sin θ), the intrinsic radian frequency σ = gk tanh(kD), the water depth D. For random waves, these expressions are easily extended by summing over the spectrum and replacing E by the elementary variancce E(f, θ)dθdf (see Appendix C).

The horizontal Stokes drift vector (U s , V s ) is given by,

(U s , V s ) = σk(cos θ, sin θ)E cosh(2kz + 2kh) sinh 2 (kD) . ( 17
)
At the lowest order (e.g. Ardhuin et al., 2008a), the full Stokes drift flow is non-divergent and verify,

∂U s ∂x + ∂V s ∂y + ∂W s ∂z = 0. ( 18
)
As a result, the less well-known vertical Stokes drift component is given by the horizontal divergence of (U s , V s ),

W s (z) = -U s | z=-h ∂h ∂x -V s | z=-h ∂h ∂y - z -h ∂U s ∂x + ∂V s ∂y dz. ( 19
)
In adiabatic conditions, the only other term is the wave-induced mean pressure J,

J = g kE sinh(2kD) . ( 20
)
In the coupled system, the horizontal Stokes velocity is computed in the coupler from the frequency spectrum of the surface Stokes drift, which is provided by the wave model, so that the wave model does not need to know the depths of the flow model levels. This also allows to force the flow model with a stored wave output that is independent of the flow model vertical resolution (see http://tinyurl.com/2wr6hoa for details). The vertical component W s is obtained by solving equation ( 19).

No definite theory exists for the force induced by wave dissipation ( F d,x , F d,y ), as only the depth-integrated force is known (e.g. [START_REF] Smith | Wave-current interactions in finite-depth[END_REF]. Still it is important to follow some general principles. Namely, the energy lost by waves when breaking is the source of surface turbulence and the corresponding wave momentum is the force ( F d,x , F d,y ). One should be careful to avoid double counting by making sure that the wind to wave momentum flux is not included in the wind stress used as a surface boudary condition for the flow model.

An empirical parametrization for the vertical profile must be used. We may clearly distinguish between the force due to wave breaking and that due to bottom dissipation [START_REF] Walstra | Calculation of wavedriven currents in a 3D mean flow model[END_REF]. We know S oc the amount of energy given up by waves as they break, either in finite depth or deep water (e.g. [START_REF] Thornton | Transformation of wave height distribution[END_REF][START_REF] Ardhuin | Semi-empirical dissipation source functions for wind-wave models: part I, definition, calibration and validation[END_REF], and S bf the loss of energy due to bottom friction (e.g. [START_REF] Ardhuin | Swell transformation across the continental shelf. part I: Attenuation and directional broadening[END_REF]. With a strong vertical mixing due to breaking waves the vertical distribution of the momentum source is not very important [START_REF] Rascle | Drift and mixing under the ocean surface. a coherent one-dimensional description with application to unstratified conditions[END_REF]. Wave dissipation may be provided as surface stress, with a vertical profile given by the delta function

δ z,b η , ( F d,x , F d,y )(z) = (τ oc,x , τ oc,y )δ z,b η = k σ (cos θ, sin θ)S oc (f, θ)δ z,b η df dθ, (21)
where S oc (f, θ) is the spectral density of the waves-to-ocean energy flux, equal to the dissipation source function in the spectral wave energy balance (e.g. [START_REF] Ardhuin | Semi-empirical dissipation source functions for wind-wave models: part I, definition, calibration and validation[END_REF].

Boundary conditions at the bottom

The bottom friction is absent in the test cases presented here. However, for case with bottom friction, the following equations can be used.

Starting from the bottom, at z = -h, for a non-resolved wave bottom boundary layer, the momentum lost by waves due to bottom friction is lost in the bottom (Longuet-Higgins, 2005) and should not be added in the water column (ie. ( F b,x , F b,y ) = (0, 0)), and the horizontal velocity should be prescribed as velocity at the bottom given by the streaming solution of [START_REF] Longuet-Higgins | Mass transport under water waves[END_REF], still approximately valid for turbulent bottom boundary layers (e.g. [START_REF] Marin | Eddy viscosity and Eulerian drift over rippled beds in waves[END_REF],

( u, v)| z=-h = 1.5 (U s , V s )| z=-h (22)
and the vertical velocity is naturally

w = -u ∂h ∂x -v ∂h ∂y . ( 23
)
We note that many parameterizations have been proposed for the bottom stress K z ∂ u ∂z where K z is the (varying) eddy viscosity. Eq. ( 22) could be used for the quasi-Eulerian velocity at the bottom.

If the wave bottom boundary layer were resolved then the appropriate bottom boundary condition is: ( u, v)| z=-h = (0, 0). In this case, one should introduce near the bottom the source of momentum

( F b,x , F b,y )(z) = k σ (cos θ, sin θ)S bf (f, θ)G(z)df dθ, ( 24
)
where G(z) is a function that integrates to 1 across the wave bottom boundary layer. This may be re-written

( F b,x , F b,y )(z) = (τ wb,x , τ wb,y )G(z). ( 25
)
The wave bottom stress vector (τ wb,x , τ wb,y ) corresponds to the momentum lost by the wave field via bottom friction and can be computed by the wave model. According to [START_REF] Walstra | Calculation of wavedriven currents in a 3D mean flow model[END_REF], the vertical distribution function G(z) can be defined such as

   G(z) = 1 - D + z -η δ wbbl if -D + η + δ wbbl ≥ z ≥ -D + η, G(z) = 0 if z > -D + η + δ wbbl . ( 26
)
where δ wbbl is the wave bottom boundary layer thickness.

Another important condition is the energy lost by waves due to bottom friction, which is a source of turbulent kinetic energy in the bottom boundary layer. The total energy is the same integral as eq. ( 24), without the k σ factor, and replacing the profile G(z) by a parametrization follow that may [START_REF] Mellor | Oscillatory bottom boundary layers[END_REF].

Boundary conditions at the surface

At the surface, the stresses are imposed, giving the upper boundary condition for the turbulent momentum flux,

K z ∂ u ∂z = τ a,x -τ aw,x (27)
where τ a,x and τ aw,x are, respectively, the x-component of the wind stress and of the wave-supported stress

(τ aw,x , τ aw,y ) = k σ (cos θ, sin θ)S atm (f, θ)df dθ, (28)
where S atm (f, θ) is the spectral density of the wind to wave energy flux (e.g. [START_REF] Ardhuin | Observation and estimation of Lagrangian, Stokes and Eulerian currents induced by wind and waves at the sea surface[END_REF], approximately equal to the input source function in the spectral wave energy balance. Here again the other boundary condition for the flux of turbulent kinetic energy is given by the loss of wave energy due to breaking, and wave-turbulence interaction, and the same integral as ( 21)

without the k σ factor. This flux may also be distributed as a near-surface source.

Lateral boundary conditions

When open boundary conditions are used, one may impose a zero mass flux to facilitate the numerical convergence [START_REF] Rascle | Impact of waves on the ocean circulation (impact des vagues sur la circulation océanique)[END_REF], which takes the form, ( u, v) = (-U s , -V s ). ( 29)

where (•) denotes the depth-integrated variable.

Some details about the pratical implementation

The quasi-Eulerian equations are implemented in the flow model MARS3D [START_REF] Lazure | An external-internal mode coupling for a 3d hydrodynamical model for applications at regional scale (MARS)[END_REF]. MARS3D is coupled with the wave model WAVEWATCH III [START_REF] Tolman | Effects of observation errors in linear regression and bin-average analysis[END_REF][START_REF] Tolman | User manual and system documentation of WAVEWATCH-III TM version 3[END_REF] thanks to the coupler PALM [START_REF] Buis | PALM: A computational framework for assembling high performance computing applications[END_REF]. Parts of the additional terms in the momentum equations are not computed inside the flow model. The terms ∂J ∂x , ∂J ∂y , F d,x , F d,y are computed inside the coupler with fields from the wave model, which are the wave mean induced pressure, the significant wave height, the wave-to-ocean stress and with fields from the circulation model: the mean sea surface elevation and the sigma levels. The other terms for the momentum equations are calculated inside the MARS3D model because the knowledge of the current, the coriolis term, the thickness of wave bottom boundary layer and the roughness length is required. In this view, the following fields from waves are transmitted by the coupler to MARS3D: the three components of the Stokes drift, the norm of the orbital velocity near the bottom, the root mean square wave height, the intrinsic radian frequency, the wave number, the wave direction, the wave dissipation due to bottom friction, the wave mean induced pressure derivatives, the F d,x and F d,y terms and one term used with open boundary conditions that is equal to -J g . The fields exchanged between the models are shown on Figure 4.

We modify the depth-integrated continuity equation in order to implement the new surface boundary condition (see eq (A.3)). We use the depthintegrated horizontal components of the Stokes drift which are computed inside the MARS3D model from the Stokes drift, provided by the coupler.

The coupled model can be used in one-way or two-way coupling mode. Only the communications drawn with solid arrows (see Figure 4) are used when the one-way mode is applied. For the two-way coupling mode, all communications are activated.The coupling procedure starts with the hydrodynamical model which sends the currents and sea surface height to the wave model via the coupler. Then, the wave model computes all the wave parameters and sends them to the coupler. Using the sea surface height and the sigma levels which come from hydrodynamical model, the forcing terms are computed and sent to MARS3D. The exchange between the models are done at a coupling time step which can be larger than the time steps of the individual models. The coupled model has been designed to run in parallel across multiple processors in order to reduce the computing time. For the moment the spatial grids are the same for the both models and we have not used the internal re-gridding capabilities of PALM.

Summary of new terms introduced

The forcing of the wave field on the ocean circulation requires the knowledge of all the fields listed in Table 1.

Compared to equations for the Lagrangian mean velocity, such as those by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF], the amount of wave forcing data to be transferred is signif- 2D 25

Table 1: List of wave-forcing terms required to force an ocean circulation model solving for the quasi-Eulerian velocity. The J term is a 2D field when the effect of the vertical shear of the quasi-Eulerian current is neglected, as done here. In general J is a 3D forcing field (Ardhuin et al., 2008a). The terms τ wb,x and τ wb,y are only used when the wave bottom boundary layer is resolved.

icantly reduced, since the latter form requires the 3D fields S xx , S yy and S xy , as well as the 3D fields U s and V s to correct the velocities before applying the turbulence closure [START_REF] Walstra | Calculation of wavedriven currents in a 3D mean flow model[END_REF]. This lower complexity of the quasi-Eulerian equations for the 3D case is contrary to the 2D case, in which seven 2D fields are needed, versus 3 to 5 (if properly dealing with the bottom boundary condition) for the depth-integrated Lagrangian equations. In both cases, for a full consistency of the ocean circulation and wave model, one should also use the wind stress of the wave model, as discussed by [START_REF] Janssen | Impact of the sea state on the atmosphere and ocean[END_REF], and a proxy of the breaking wave heights, possibly the wind sea wave height [START_REF] Rascle | A global wave parameter database for geophysical applications. part 1: wave-current-turbulence interaction parameters for the open ocean based on traditional parameterizations[END_REF].

The wave model can also be used to provide energy fluxes for the surface flux of turbulent kinetic energy [START_REF] Janssen | Impact of the sea state on the atmosphere and ocean[END_REF], or the near-bottom flux of turbulent kinetic energy (TKE) due to bottom friction [START_REF] Mellor | Oscillatory bottom boundary layers[END_REF]. For the adiabatic conditions considered here all the stresses τ aw , τ wb , τ oc are zero, together with these fluxes of TKE.

4. An adiabatic test case for 3D wave-current models

Description of numerical set-up

This adiabatic test has a known numerical solution. It is adapted from Ardhuin et al. (2008b) and corresponds to steady monochromatic waves shoaling from 4 to 6 m depth on a slope without breaking nor bottom friction, and for an inviscid fluid. Here the bottom is symmetric with the bottom sloping back down to 6 m, in order to allow periodic boundary conditions if needed (see Figure 6). Both waves and bottom topography are uniform in the y-direction. The flow is confined to a channel with free-slip boundary conditions (North and South boundaries). The East and West boundaries are open (see Figure 6).

Several wave conditions are tested:

-Significant wave height: H s = 0.34 m or H s = 1.02 m, -Wave period: T = 5.24 s or T = 13 s, -Wave direction: θ = 90 • , -Wave steepness: 1 = 0.0266 or 1 = 0.0798.

Our MARSD model configuration uses 100 sigma levels regularly spaced, 5 active points in the y-direction and 78 active points in the x-direction. The time step was set to 0.05 s for the H s = 1.02 m tests (and 1 s for H s = 0.34 m). The ARB08 model has been tested with one active point in y-direction and the time step was set to 1 s for all numerical simulations. For the sake of simplicity, the wave model time step is taken equal to the flow model time step.

For these shoaling waves the group velocity varies a little (5.4%) from 4.89 m.s -1 to 4.64 m.s -1 , due to the fact that the non-dimensional depth kD is close to unity. Because the current is much less than the group speed, the waves propagate with a nearly constant energy flux, resulting in a small increase of wave amplitude, by 2.7%, in the shallower part of the domain.The Eulerian analysis of such a situation was given by [START_REF] Longuet-Higgins | On the wave-induced difference in mean sea level between the two sides of a submerged breakwater[END_REF]. With H s = 0.34 m and T = 5.24 s, the Longuet-Higgins solution gives a mean water level 0.32 mm lower in the shallow region, and both studies by [START_REF] Rivero | On the vertical distribution of u w[END_REF] and [START_REF] Lane | Wave-current interaction: A comparison of radiation-stress and vortex-force representations[END_REF], clearly show that there is no other dynamical effect: the Eulerian mean current is steady and simply compensates for the divergence of the wave-induced mass transport (see Figure 7). Because the relative variation in phase speed is more important, from 6.54 to 5.65 m.s -1 , it produces a strong divergence of the Stokes drift (see Figure 5), which accelerates in shallow water. The quasi-Eulerian velocity is irrotational, thus nearly depth-uniform, and compensates the Stokes drift divergence by a strong convergence. This situation is a stationary solution. Eulerian and Lagrangian velocities Although the inconsistency of most equations for the Lagrangian velocity have been described in Ardhuin et al. (2008b), no result have yet been published on the actual consequences of this inconsistency. We chose to illustrate these consequences with the equations by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] that have been used in other studies, unfortunately without a detailed analysis.

A first important test of the model is the solution for the mean sea surface elevation, which is given by [START_REF] Longuet-Higgins | On the wave-induced difference in mean sea level between the two sides of a submerged breakwater[END_REF]. This is correctly reproduced by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] equations and the quasi-Eulerian equations. We will now consider the currents. The reference Lagrangian velocity using for the comparison is obtained from the quasi-Eulerian analysis as U = u + U s (see Figure 7). The reference solution exhibits a vertical shear that is entirely due to the Stokes drift (see Figure 5) and the quasi-Eulerian velocity is homogeneous over the water column. The MARS3D model is now used to solve either the M03 equations or the eqs. ( 11) and ( 12). In the case of the M03 equations the wave induced forcing has a strong vertical gradient which is unrealistic (see Figure 1). As a result, a strong surface velocity develops, associated to a counter-current below of comparable magnitude (Figure 8). This circulation pattern is very different from the known correct solution (Figure 7). In spite of the small bottom slope and wave steepness, in only 15 minutes the resulting velocity reaches 17 cm s -1 , which is about 10 times the correct solution shown in Figure 7. Further, if the model is integrated for a longer time, the region of positive acceleration on the up-slope meets the region of negative acceleration on the down-slope, resulting in large vertical velocities and further strange model adjustments. This erroneous velocity is clearly driven by the erroneous wave-induced force term F , which is proportional to gDε 2 1 ε 2 , where ε 1 is the wave steepness, ε 2 is the bottom slope (Ardhuin et al., 2008a). For the bottom shape and wave period chosen here, the maximum value of F is 0.29gDε 2 1 ε 2 . Obviously, the depth dependence of F plays an important role. Since F becomes depth-uniform for kD → 0, one may expect that the problem could vanish in shallow water. Unfortunately, in practice, the velocity at which the current first stabilizes (here after 15 minutes), is independent of ε 2 , provided that the change in water depth remains the same. If the bottom topography is stretched by a factor 1/α in the x direction, the slope increases by a factor α and the change advection compensates the local increase of F . Mathematically, equation (10) follows a Froude scaling: when x is replaced by x = αx and t by t = α 2 t, the equation is unchanged if F = αF , and thus U (x , t ) = U (x, t). As a result, for any wave field approaching the shore from deep water, even on a very gently sloping continental shelf, there will be a very large spurious onshore velocity at the surface. Based on the present case, this velocity can exceed 10 times the Stokes drift. This momentum is generated where kD ∼ 1, and self-advects onshore. Obviously, some realistic mixing will reduce this effect. Using a realistic constant eddy viscosity of 2.8×10 -3 m 2 .s -1 only reduces the current by about a factor 2 to 3 (see Table 3). This factor depends on the wave amplitude since the introduction of viscosity breaks the Froude scaling. We also test the ARB08 model mainly based on equations ( 11) and ( 12). This model solve now the quasi-Eulerian velocity (Lagrangian velocity minus Stokes velocity). In order to compare with the Lagrangian reference solution (see Figure 7, first panel), we add the Stokes velocity to the model solution.

K z = 0 m 2 .s -1 . H s (m) T p (s) K z (m 2 .s -1 ) resulting U (m.s -
The ARB08 model give a quasi-Eulerian current solution u (see Figure 9), and thus differs from the M03 model which give the Lagrangian current. The dynamic effect of the waves is conveyed by dynamic forcing terms (see terms A 1 , A 2 , B 1 , B 2 in equations ( 11) and ( 12)) and an equivalent mass source at the surface (eq. ( 16)). The quasi-Eulerian current is nearly depth-uniform as expected (Figure 9). As a result, the Lagrangian current given by the ARB08 model (Figure 9, first panel) is similar to the reference current (Figure 7, first panel). The flow structure and the intensity of the flow are within a fraction of a percent. So, we can conclude that the ARB08 model correctly simulates for this case the three-dimensional flow in presence of waves.

Conclusion

It was demonstrated here that equations for the three-dimensional waveforced circulation that are formulated in terms of the Lagrangian mean velocity (total momentum) and use analytical functions of the local wave field and topography can produce very large spurious velocities. This result was anticipated by Ardhuin et al. (2008b) who showed that the vertical flux of momentum is a priori not a function of the local water depth, and can be estimated from non-local evanescent wave modes that contribute to the general solution of the Laplace equation over a sloping bottom. The magnitude of the problem is revealed by the present study. These come from poor approximation of the vertical flux of wave momentum, and they persist errors when a reasonable vertical mixing is included. These error may become negligible in the surf zone, but they may also play a big part in the differences in vertical velocity profiles reported by [START_REF] Haas | Comparing a quasi-3d to a full 3d nearshore circulation model: SHORECIRC and ROMS[END_REF], when comparing a version of ROMS solving the [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] equations with SHORECIRC (see their figure 4). We also wish to point out another common source of differences between model results. Some models, like SHORCIRC solve for the quasi-Eulerian mean velocity U -U s while ROMS solved for the Lagrangian mean velocity U . The difference between the two is the Stokes drift, which can be very large in the surf zone, up to 30% of the wave phase speed (Ardhuin et al., 2008b).

From the present model results, we conclude that today there is no acceptable short-cut to a three-dimensional equation for the Lagrangian velocity U : the only possibility would be to solve for the wave motion to first order in the bottom slope. This can be done with a model of the kind developed by [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small amplitude water waves over variable bathymetry regions[END_REF] and [START_REF] Gerosthathis | Coupledmode, phase-resolving model for the transformation of wave spectrum over steep 3d topography. a parallel-architecture implementation[END_REF], with at least 10 vertical modes. Given the large effort required for a 4 by 4 km region with only 3 modes [START_REF] Magne | Evolution of surface gravity waves over a submarine canyon[END_REF], this is hardly a practical solution. The only practical solution avalaible today is thus the use of a momentum equation for the quasi-Eulerian velocity, such as proposed by [START_REF] Mcwilliams | An asymptotic theory for the interaction of waves and currents in coastal waters[END_REF], Newberger and Allen (2007b), or Ardhuin et al. (2008a). This approach has been applied to surf zone problems by [START_REF] Rascle | Impact of waves on the ocean circulation (impact des vagues sur la circulation océanique)[END_REF] and [START_REF] Uchiyama | Wave-current interaction in an oceanic circulation model with a vortex-force formalism: Application to the surf zone[END_REF].

INSU and grant ANR-BLAN-08-0330-01, and F.A. is supported by a FP7-ERC grant 240009 "IOWAGA" and the NOPP program through ONR grant N00014-10-1-0383. where

Appendix

• ς = z -η D
is the sigma coordinate with η the mean elevation, h the bottom depth and D = η + h the mean water column depth,

• W = ∂ς ∂t + u ∂ς ∂x + v ∂ς ∂y + w D , • ∂ς ∂x = - 1 D • ∂ η ∂x - ς D • ∂D ∂x , • ∂ς ∂y = - 1 D • ∂ η ∂y - ς D • ∂D ∂y , • ∂ς ∂t = - (ς + 1) D • ∂ η ∂t ,
• ( F b,x , F b,y ) are only used when the wave bottom boundary layer is resolved.

The depth-integrated continuity equation becomes where ω is defined by equation (A.5) and ( F b,x , F b,y ) are only used when the wave bottom boundary layer is resolved.

The continuity equation becomes

∂ η ∂t + ∂ [D( u + U s )] ∂x + ∂ [D( v + V s )] ∂y + ∂ω ∂ς = 0, (B.3)
and the evolution of a conservative passive tracer concentration C is defined by equation (A.4).

We notice that only the horizontal components of stokes drift (U s , V s ) must be known with this formulation. The W s depending terms are removed.

Appendix C. Explicit form of random wave forcing terms for the quasi-Eulerian velocity

For random waves, eq. ( 17) becomes (U s , V s ) = σk(cos θ, sin θ)E(f, θ) cosh(2kz + 2kh) sinh 2 (kD) df dθ, (C.1)

where E(f, θ) is the spectral density of the surface wave elevation variance, usually known as the wave spectrum, the state variable of most numerical wave models, and the wave-induced pressure term becomes, J = g kE(f, θ) sinh 2kD df dθ. (C.2)

Figure 1 :

 1 Figure 1: Example profiles of net wave-induced forcing over a sloping bottom in the test case proposed by Ardhuin et al. (2008a) (H s = 1.02 m, T = 5.26 s, no mixing. The bottom slope is 2 = 0.0798. More details in section 4). ς is the terrain-following vertical coordinate with ς = 0 at the surface and ς = -1 at the bottom.

Figure 2 :

 2 Figure 2: Relationships between wave-averaged theories according to their choice of momentum variable and depth integration. An arrow from a to b indicates a derivation link: b can be derived from a. In the case of McWilliams et al. (2004) and Ardhuin et al.(2008b), the second derivative in the vertical current profile is neglected in the second, while the first is only adiabatic. Names that appear in red correspond to theories that are not fully consistent with their originating hypotheses. In the case of 3D theories for the total momentum, the problem generally comes from the vertical flux of momentum and may only arise on a sloping bottom, not explicitly considered by[START_REF] Groeneweg | Wave-current interactions in a generalized Lagrangian mean formulation[END_REF].

Figure 4 :

 4 Figure 4: Coupling process between the wave model (WAVEWATCH III) and the hydrodynamical model (MARS3D) via the coupler (PALM).

Figure 5 :

 5 Figure 5: Stokes velocity in x-direction for H s = 1.02 m, T = 5.26 s and K z = 0 m 2 .s -1 .

Figure 6 :

 6 Figure 6: Bathymetry and configuration

Figure 7 :

 7 Figure7: Solution given by quasi-Eulerian analysis for the inviscid sloping bottom case with H s = 1.02 m, T = 5.24 s and without mixing. Lagrangian velocity U (first panel, contours are equally spaced from -0.01 to 0.025 m s -1 ) and quasi-Eulerian velocity U -U s (second panel). The thick black line is the bottom elevation.

Figure 8 :

 8 Figure 8: M03 model: Solution given by the coupled model after 15 minutes of integration: Lagrangian velocity U (first panel) and Lagrangian velocity minus Stokes drift U -U s (second panel). The thick black line is the bottom elevation. H s = 1.02 m, T = 5.26 s, K z = 0 m 2 .s -1 .

Figure 9 :

 9 Figure 9: ABR08 model: velocity u (first panel) and Quasi-Eulerian velocity u(m/s) (second panel) which is equal to the Lagrangian velocity minus the Stokes drift. H s = 1.02 m, T = 5.24 s and K z = 0 m 2 .s -1 . The thick black line is the bottom elevation.

•

  denotes depth-integrated variable.The evolution of a conservative passive tracer concentration C is, The modified vertical velocity ω is defined byω = w + W s + D ∂ς ∂t + D( u + U s ) ∂ς ∂x + D( v + V s )

Table 2 :

 2 Overview on the numerical simulations

Table 3 :

 3 M03 model: Surface velocity at x =200 m for different model settings. The settings corresponding to the test in Ardhuin et al. (2008b) are given in the second line H s = 0.34 m, T = 5.6 s, K z =0 m 2 s -1 . The surface velocity values are written for T = 900 s where H s = 1.02 m and for T = 2700 s where H s = 0.34 m.

  A. glm2 equations in sigma coordinates : Momentum, mass, tracer conservation Let (x, y, z, t) denote the cartesian coordinate system and (x , y , ς, t ) the sigma coordinate system.

	∂ u ∂t	+ u	∂ u ∂x	+ v	∂ u ∂y		+ W	∂ u ∂ς	-f v +	1 ρ	∂p H ∂x	+	∂p H ∂ς	•	∂ς ∂x
		= f +		∂ v ∂x	+	∂ v ∂ς	•	∂ς ∂x	V s -	∂ u ∂y	+	∂ u ∂ς	•	∂ς ∂y	V s
		-	W s D	•	∂ u ∂ς	-	∂J ∂x	-	∂J ∂ς	•	∂ς ∂x	+ F d,x + F m,x + F b,x ,	(A.1)
	and
	∂ v ∂t	+ u	∂ v ∂x	+ v	∂ v ∂y	+ W	∂ v ∂ς	+ f u +	1 ρ	∂p H ∂y	+	∂p H ∂ς	•	∂ς ∂y
		= -f +		∂ v ∂x	+	∂ v ∂ς	•	∂ς ∂x	U s +	∂ u ∂y	+	∂ u ∂ς	•	∂ς ∂y	U s
		-	W s D	•	∂ v ∂ς	-	∂J ∂y		-	∂J ∂ς	•	∂ς ∂y	+ F

d,y + F m,y + F b,y . (A.2)

  • F cm is the turbulent mixing.

	∂D v ∂t	+	∂[D( u + U s ) v] ∂x	+	∂[D( v + V s ) v] ∂y	+	∂(ω v) ∂ς	+ Df u +	D ρ	•	∂p H ∂y	+	∂p H ∂ς	∂ς ∂y
		= V s D	∂ v ∂y	+	∂ v ∂ς	∂ς ∂y	-U s D f -	∂ u ∂y	-	∂ u ∂ς	∂ς ∂y	-D	∂J ∂y	-D	∂J ∂ς	∂ς ∂y
			+ F d,y + F m,y + F b,y .						(B.2)
	Appendix B. Flux formulations of the quasi-Eulerian glm2 equa-
				tions in sigma coordinates : Momentum, mass, tracer
				conservation
	∂D u ∂t	+	∂[D( u + U s ) u] ∂x	+	∂[D( v + V s ) u] ∂y	+	∂(ω u) ∂ς	-Df v +	D ρ	•	∂p H ∂x	+	∂p H ∂ς	∂ς ∂x
		= U s D	∂ u ∂x	+	∂ u ∂ς	∂ς ∂x	+ V s D f +	∂ v ∂x	+	∂ v ∂ς	∂ς ∂x	-D	∂J ∂x	-D	∂J ∂ς	∂ς ∂x
			+ F d,x + F m,x + F b,x ,						(B.1)
	and

The wave pseudo-momentum is defined as a quantity that only involves the zeromean displacement of the water particles, and may differ from other definitions that could include the mean flow response, as explained by[START_REF] Mcintyre | On the 'wave momentum' myth[END_REF]. For simplicity, we shall call 'momentum' the pseudo-momentum.

In that work, the control volume is identical to the one in M03, but an Eulerian average of the pressure is taken, making the averaged equations inconsistent. The consequences of this error are discussed in[START_REF] Bennis | Comments on the depth-dependent current and wave interaction equations: A revision[END_REF].

These sides actually have an angle in the case of the Generalized Lagrangian Mean (GLM), giving more complex equations.

Although derived via a different approach, the resulting equations are mathematically equivalent to the ones given by[START_REF] Uchiyama | Wave-current interaction in an oceanic circulation model with a vortex-force formalism: Application to the surf zone[END_REF], except for the Bernoulli head term (our J) in which we have neglected the effects of the current profile. Here they were derived from a quasi-Eulerian average of the flow, as detailed inArdhuin et al. (2008b), which has a clear physical interpretation as the Lagrangian mean minus the Stokes drift.

Although the vertical component of the Stokes drift may not be familiar to the reader, it appears, just like the horizontal components, in the general definition of the wave pseudomomentum[START_REF] Andrews | An exact theory of nonlinear waves on a Lagrangian-mean flow[END_REF] Ardhuin et al., 2008b). In particular for inviscid conditions over a sloping bottom it is physically obvious that the drift of water particles must follow the bottom and thus must have a vertical component. In practice W s can be computed from (U s , V s ) as the full Stokes drift flow is approximately non-divergent(Ardhuin et al., 2008b)

Acknowledgments

The Authors express their gratitude to P. Marsaleix for his proof-reading and advice. A-C. B. acknowledges the support of a post-doctoral grant from

Appendix D. Forcing terms for the Lagrangian mean velocity

The force that is given by the horizontal divergence of the flux of horizontal momentum (S xx ) is

Using Airy theory, S xx is given by,

and the vertical profile function F CS changes with f and is defined by

with similar definitions for F SS (respectively F CC ), replacing cosh in the numerator (respectively sinh in the denominator) by sinh (respectively cosh).

The horizontal force that is given by the vertical divergence of the flux of vertical momentum (S x3 ) is

In this case, Airy theory is insufficient for a consistent approximation. Yet [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] still used Airy theory, thus producing the erroneous expression,