Three-dimensional modelling of wave-induced current from the surf zone to the inner shelf

. We develop and implement a new method to take into account the impact of waves into the 3-D circulation model SYMPHONIE (Marsaleix et al., 2008, 2009a) following the simpliﬁed equations of Bennis et al. (2011), which use glm2z-RANS theory (Ardhuin et al.,


Abstract.
We develop and implement a new method to take into account the impact of waves into the 3-D circulation model SYMPHONIE (Marsaleix et al., 2008(Marsaleix et al., , 2009a following the simplified equations of , which use glm2z-RANS theory (Ardhuin et al., 2008c). These adiabatic equations are completed by additional parameterizations of wave breaking, bottom friction and wave-enhanced vertical mixing, making the forcing valid from the surf zone through to the open ocean. The wave forcing is performed by wave generation and propagation models WAVEWATCH III ® (Tolman, 2008(Tolman, , 2009Ardhuin et al., 2010) and SWAN (Booij et al., 1999). The model is tested and compared with other models for a plane beach test case, previously tested by Haas and Warner (2009) and Uchiyama et al. (2010). A comparison is also made with the laboratory measurements of Haller et al. (2002) of a barred beach with channels. Results fit with previous simulations performed by other models and with available observational data.
Finally, a realistic case is simulated with energetic waves travelling over a coast of the Gulf of Lion (in the northwest of the Mediterranean Sea) for which currents are available at different depths as well as an accurate bathymetric database of the 0-10 m depth range. A grid nesting approach is used to account for the different forcings acting at different spatial scales. The simulation coupling the effects of waves and currents is successful to reproduce the powerful northward littoral drift in the 0-15 m depth zone. More precisely, two distinct cases are identified: When waves have a normal angle of incidence with the coast, they are responsible for complex circulation cells and rip currents in the surf zone, and when they travel obliquely, they generate a northward littoral drift. These features are more complicated than in the test cases, due to the complex bathymetry and the consideration of wind and non-stationary processes. Wave impacts in the inner shelf are less visible since wind and regional circulation seem to be the predominant forcings. Besides, a discrepancy between model and observations is noted at that scale, possibly linked to an underestimation of the wind stress.
This three-dimensional method allows a good representation of vertical current profiles and permits the calculation of the shear stress associated with waves and currents. Future work will focus on the combination with a sediment transport model. navigational issues, coastal management and prediction of hazardous zones for swimmers. A wide variety of modelling techniques have been applied to the surf zone, based on depth-integrated equations. These include phase resolving (e.g. Chen et al., 2003;Clark et al., 2011), group-averaged (Reniers et al., 2004), or fully phase-averaged. These models are unfortunately not well adapted for continental shelf processes, which are influenced by stratification, making it difficult to model cross-shore transport phenomena uniformly from the beach to the shelf break. Recently, developed 3-D wave-current interaction theories (e.g. Mellor, 2003;McWilliams et al., 2004;Ardhuin et al., 2008c) may provide useful approaches for this problem.
Following the pioneering modelling work by Rascle (2007); Newberger and Allen (2007); Uchiyama et al. (2010), we investigate the influence of surface waves on ocean circulation in the inner shelf and surf zone. The main influences of waves on currents occur through bottom (e.g. Komar et al., 1972) and surface stresses (e.g. Donelan et al., 1993), while turbulent kinetic energy at the surface is enhanced by wave breaking (e.g. Agrawal et al., 1992). Waves are associated with mean momentum that can be observed as a surface-intensified drift velocity (Stokes, 1847). In deep water, this drift is highly correlated with the wind speed and wave height, with a magnitude of the order 6 × 10 −4 U 2 10 where U 10 , the 10 m wind speed, is in m s −1 ). In the surf zone, the drift is not correlated with wind speed and can reach as much as 30 % of the phase speed, with a strong surface intensification (Ardhuin et al., 2008c). The actual drift of water particles is the sum of this Stokes drift and the Eulerian current, with wave influences also on the Eulerian current (e.g. Xu and Bowen, 1994). Reciprocally, currents can modify waves by, refraction, partial reflection, up to blocking (Smith, 1975;Chawla and Kirby, 2002).
A first approach to the study of wave/current interactions can be to add certain effects in particular contexts (e.g. Jordà et al., 2007). For example, Mastenbroek et al. (1993) improve their numerical simulations of storm surges by introducing a wave-dependent drag coefficient for the wind.
During the 1960s, Taylor (1962) and Whitham (1962) focused on non-linear wave properties. These works then led to the radiation stress theory, which was first introduced by Longuet- Higgins and Stewart (1962), and then by Phillips (1977). This theory takes into account the excess flow of momentum due to the presence of waves in the barotropic momentum equations for the total current, thanks to the addition of radiation stress gradients. In radiation stress theory, wave and current momenta are combined and the effect of the waves is applied to this total momentum. Although this is practical for depth-integrated flows, it becomes a problem in 3-D models, in particular because the Stokes drift is not mixed and is often the main source of vertical shear near the surface, with important consequences for surface drift (e.g. Rascle and Ardhuin, 2009). Instead, the problem can be for-mulated for the current momentum only, as shown by Garrett (1976) for depth-integrated equations, and Andrews and McIntyre (1978) in the most general form. Several theories have been developed and applied for the full momentum (involving radiation stresses) or the current momentum only (in which a vortex force appears). Although much work is still to be developed for the proper treatment of turbulence in the presence of waves, several papers have established that all published theories that use radiation stresses have some errors at the leading order, which may cause spurious circulations (Ardhuin et al., 2008b;Kumar et al., 2011;. Here we shall use a formulation of the "current momentum", which is formally defined as the quasi-Eulerian momentum (Andrews and McIntyre, 1978;Jenkins, 1989), namely the Lagrangian mean velocity minus the Stokes drift.
Here we use an approximation of the exact equations from (Andrews and McIntyre, 1978) to second order in wave slope, including a transformation to cartesian coordinates (glm2-z approximation, Ardhuin et al., 2008c), in a simplified form that neglects the vertical current shear effect on the dynamic pressure . In the limit of weak vertical current shear, these equations are formally equivalent to the Eulerian-mean equations of McWilliams et al. (2004) that are based on an analytic continuation across the air-sea interface, and which have been used by Uchiyama et al. (2010).
Previous studies dealing with wave-current interaction are often focused on the surf zone (with water depths on the order of 1 m) (e.g. Uchiyama et al., 2010;Haas and Warner, 2009;Weir et al., 2011). Few studies and measurements have been dedicated to the mid-shelf zone (with water depths of order 100 m), or at least to the inner shelf (between the surf zone and mid shelf): Lentz et al. (1999Lentz et al. ( , 2008 were one of the first to study the influence of waves on the inner shelf. The purpose of this article is to extend the study of wave-current interaction to both the inner shelf and the open ocean by implementing the new set of equations of  in the primitive equation model SYMPHONIE (Marsaleix et al., 2008). By using a nested strategy, which allows studies at all scales and by completing the model with additional parameterizations of wave breaking, bottom friction and waveenhanced vertical mixing, we ensure that the forcing is valid from the surf zone through to the open ocean. We test and compare our model to measurements made on the Têt inner shelf during a typical winter storm. This inner shelf has a bathymetry made up of complex sandbar systems, therefore before tackling this real case, the accuracy of the model is first assessed in two idealized test cases.
Wave and circulation models, modified with the formulation of , are presented in Sect. 2. Section 3 describes two academic test cases of the surf zone. The first is on an idealized plane beach submitted to obliquely incident spectral waves (a case also tested by Warner, 2009, andUchiyama et al., 2010 Haller et al., 2002). Section 4 focuses on the 21 February 2004 storm in the Têt inner shelf. The simulated currents are compared to the observed ones to assess model accuracy. Finally, Sect. 5 provides a summary and conclusion.

Coastal circulation model
We used the Boussinesq hydrostatic circulation ocean model described in Marsaleix et al. (2008Marsaleix et al. ( , 2009a. Components of currents, temperature and salinity are computed on a Cgrid using an energy-conserving finite-difference method. A generalized sigma coordinate (Ulses et al., 2008c) is used in order to refine resolution near the bottom and the surface. A complete description of the bulk formulae used to compute the air/sea fluxes is given in Estournel et al. (2009). The so called SYMPHONIE model has been extensively used in studies of the Mediterranean Sea, mostly at the scale of continental shelves (Ulses, 2005;Estournel et al., 2003Estournel et al., , 2005, generally comparing satisfactorily with available in-situ observations. Leredde and Michaud (2008) however found that the model did not perform that well for the case of an extreme meteorological event reported in the Gulf of Lion in February 2007. It was concluded that the relative failure of the simulations was likely a consequence of the lack of a proper wave/current parameterization. This study incited the present one, in which we implement this particular development in our model, following the method proposed by .

General equations
The momentum equations of the coastal circulation model are rewritten in order to take into account the wave forcing. This gives the Eqs. (18)-(21) of  which govern the evolution of the quasi-Eulerian velocities (û,v,ŵ) equal to where (u, v, w) are the mean Lagrangian velocities and (U s , V s , W s ) the Stokes drift in the horizontal (x-, y-) and vertical (z-) directions. They are valid from the bottom z = −h to the local phase-averaged free surface z =η.
∂û ∂t +û ∂û ∂x +v ∂û ∂y +ŵ The evolution of C, the concentration of a passive tracer, is then governed by and the mass conservation becomes These equations were implemented in the MARS 3-D model (Lazure and Dumas, 2008). We transform them into a discrete form by using the flux-divergence form of the advection terms which can be found in most coastal hydrodynamic models (e.g. Marsaleix et al., 2008;Blumberg and Mellor, 1987;Shchepetkin and McWilliams, 2005) (see Appendix A). In addition, the Eqs.
(2)-(3) taken from  neglected the effect of the vertical current shear.
In realistic settings, a strong vertical shear can occur, so nontrivial higher-order Bernoulli head terms must be considered. Equations then become  (3) is thus replaced by a wave-induced pressure term S J and a shearinduced pressure term S shear . The depth uniform waveinduced term S J is equal to with D = η + h the water depth, g the acceleration due to gravity, E the wave energy and k the wave number. The shear-induced term is given in Eq. (40) of Ardhuin et al. (2008c), using a wave spectrum integrated form. Here, it has been replaced by a spectrum-averaged expression around the principal frequency, i.e.
The 3-D Stokes velocities being non-divergent (Uchiyama et al., 2010), we have (Bennis et al., 2011) (Eq. 18): This leads to Assuming that U s 2 + V s 2 >> W s 2 , the vertical velocity dependent terms are consequently omitted in our calculus.
Stokes velocities are given by cosh (2kD) in shallow waters for kD < 6 and (U s , V s ) = σ k(cos θ, sin θ)E2e 2·k(z−η) in deep waters for kD > 6 with σ the relative frequency and θ the angle of wave propagation.
In realistic configurations (i.e. for random waves), we replace E by the elementary variance E(θ, σ )dθ dσ and we integrate the entire expression over the spectrum of the relative frequencies and angles of wave propagation of the wave model. The WAVEWATCH III wave model, hereinafter WW3 (version 4.04-SHOM; Tolman, 2008Tolman, , 2009Ardhuin et al., 2010), provides directly the wave-induced pressure S J and the surface Stokes velocities ((U sf (k n ), V sf (k n ) = w n k n E) discretized in the frequency spectrum, so the Stokes drift can be calculated by summing these terms over the frequency spectrum: with k n the wave numbers associated to the different frequencies of the spectrum. w n are calculated by w n = √ gk n tanh(k n D) and P (z; k n ), the vertical profiles associated with the different frequencies, are defined by P (z; k n ) = 2·cosh(2k n (z+h)) cosh(2k n D) in shallow waters for k n D < 6 P (z; k n ) = 2e 2·k n (z−η) in deep waters for k n D > 6.
Stokes drift is strongly sheared at the surface so a high resolution near the surface is required.
In these equations, the wave-induced dissipation force as defined by  is split into two forces: one associated with wave-breaking dissipation (bathymetric breaking and whitecapping), and one induced by bottom dissipation, F d = F surf d + F bot d . In the absence of a known vertical profile, these two forces find themselves in the boundary conditions respectively at the surface and at the bottom as surface and bottom stresses. One can thus impose an empirical vertical profile for the two forces Uchiyama et al., 2010). On the other hand, the vertical profile of velocity is possibly not really sensitive to such issues because of the smoothing effects induced by strong vertical mixing (Rascle et al., 2006). In our case, we have chosen to consider the force associated to the bottom dissipation as a bottom stress and to impose a vertical distribution function for the force associated with wave-breaking dissipation: where κ = 0.4 is the Von Karman constant; τ wo = (τ wo,x , τ wo,y ) is the momentum flux from wave to ocean linked to wave breaking (bathymetric breaking, or whitecapping); wave-turbulence interaction and viscous effects and H sw is the significant wave height of the wind-sea only. τ wo and H sw are provided by WW3. For monochromatic waves, we link τ wo to the wave dissipation due to wave breaking b In fact, this ratio is often used in the literature, and given as a 1 ρg factor by the Simulating WAves Nearshore wave model (SWAN-version 40.72;Booij et al., 1999).

Boundary conditions
Boundary conditions then become -At the surface: with K z the vertical eddy viscosity calculated by a turbulent closure scheme representing the energy cascade towards small scales. τ a = (τ a,x , τ a,y ) the wind stress and τ aw = (τ aw,x , τ aw,y ) the momentum flux from atmosphere to wave. In fact, waves influence the flux transfers from atmosphere to ocean. A part of the atmosphere momentum flux goes directly in the ocean via τ a . Another part τ aw goes into the wave field. Then this field is subjected to dissipation and releases τ wo . At a larger scale than the surf zone, τ aw and τ wo tend to cancel each other. Actually, only a small part of τ aw (5 %) is radiated into the wave field (Ardhuin et al., 2004). WW3 provides directly τ aw . In the surf zone, the term τ wo is predominant, and τ aw is neglected.
-On the bottom: We add the momentum lost by waves due to bottom friction τ wob in the bottom boundary condition of the momentum equation. Adding this wave dissipation rate permits the reproduction of the bottom streaming flow, that has the same direction as the waves.
τ bot = (τ bot,x , τ bot,y ) is the bottom stress linked to current. We consider two different parameterizations for this term. The first one is a quadratic bottom drag parameterization and is only linked to the current through with C d the drag coefficient and V b the near bottom current. The second parameterization is a drag law function linked to waves and currents, established by Soulsby et al. (1995): with τ c the bottom stress due to current only, equal to where z 0 is a length scale representing the roughness of the boundary (here the bottom boundary) and z 1 is the distance between the first level above the bottom boundary and the bottom boundary. τ w is the bottom stress linked to waves only, given by In the following, depending on the case, we will use Eq. (21) (Sects. 3.1 and 3.2) or Eq. (22) (Sect. 4). The momentum lost by waves due to bottom friction is given by with wd the wave bottom drag calculated using the parameterization of Reniers et al. (2004): wd = 1 2 √ π ρf w |u orb | 3 , u orb the bottom wave orbital velocity calculated by and f w the wave friction factor given by Myrhaug et al. (2001): with z 0 the bottom roughness length and a bw the half orbital excursion length given by a bw = |u orb |T 2π (with T the wave period).

-Lateral boundaries:
At the open boundaries, for realistic simulations, radiation conditions from Flather (1976) are applied. Technically, we follow the Eqs. (14) of Marsaleix et al. (2006). Thus, for the sea surface elevation external variable: whereû N is the velocity normal to the boundary, and "F" refers to the external forcing term. If waves are the only external forcing: Boundary conditions (Eq. 29) are deduced from the momentum equation (Eqs. 2 and 3) and some simplifying hypotheses (steady solution, non linear terms are neglected).

Wave-induced vertical mixing
Vertical mixing is parameterized according to the k − turbulent closure scheme. The vertical eddy viscosity K z is calculated by: K z = √ 2E k l k S z and the eddy diffusivity K h = www.ocean-sci.net/8/657/2012/ Ocean Sci., 8, 657-681, 2012 √ 2E k l k S h . The turbulent length l k is related to E k the turbulent kinetic energy (TKE) and to , the dissipation rate of TKE according to S z and S h the quasi-equilibrium stability functions of Kantha and Clayson (1994) depend on E k , and the Brunt-Vaisala frequency. The equations for E k and (Burchard and Bolding, 2001) are where ] is the production term and B = g ρ K h ∂ρ ∂z is the buoyancy term. σ k = 1.3, c 0 = 0.5544, c 1 = 1.44, c 2 = 1.92 and c 3 = 1 if B ≥ 0 and c 3 = −0.52 otherwise (Warner et al., 2005).
-Bottom boundary conditions for E k : The E k bottom boundary condition is based on the assumption of the equilibrium of the production and dissipation terms (P = ).
-Surface boundary conditions for E k : Alternatively, the boundary conditions can be specified as surface flux conditions, namely: K z ∂E k ∂z = F . Where the surface flux can be computed according to Craig and Banner (1994) with τ = τ a − τ aw + τ wo , or directly prescribed from the "wave to ocean" turbulence flux computed by a wave model when available (φ oc term in WW3). The surface flux condition is believed to produce more realistic results than the P = condition (Estournel et al., 2001).
-Bottom and surface boundary conditions for : The surface and bottom conditions are computed on the first level under the surface and above the bottom boundaries. Let z 1 denote the distance between this level and the considered boundary. Boundary conditions for are obtained from E k and Eq. (30), using the latter with some appropriate hypothesis for l B a boundary length scale value. A simple formulation (Warner et al., 2005) is eventually given by l B = κ(z 1 + z 0 ), where z 0 is the length scale representing the roughness of the boundary.
Unfortunately, the formulation l B = κ(z 1 + z 0 ) potentially leads to unrealistic high values, especially when the grid resolution is low. It must indeed be realized that coastal ocean models generally have to deal with strong variations of bathymetry. For s coordinate models, this unavoidably results in a loss of resolution in the deepest areas of the numerical domain. A more complete formulation is thus used in our case. Following Estournel and Guedalia (1987), the stratification and the shear effects are taken into account through the use of the Richardson number (Ri): Moreover, at the surface, we consider that z 0 = z surf is linked to the significant wave height (Terray et al., 1996(Terray et al., , 2000. We therefore tested values between 0.8H s and 2.4H s which can be found in the literature (Rascle et al., 2006). Ideally, we should not use the significant wave height, but the significant wave height of the wind-sea only H sw . Given that the swells have a small surface slope and consequently do not break, it is more appropriate to use the wave height of the wind sea only to calculate the roughness length instead of the significant wave height H s (Rascle et al., 2008). This value is calculated according to Rascle et al. (2008) (Eq. 6), and is now available in WW3. E k and can not be lower than the minimum values E k,min = 10 −8 m 2 s −2 and min = 10 −12 m 2 s −3 . Moreover the length scale limitation suggested by Galperin et al. (1988) is transposed to the dissipation rate of TKE, that is

Wave model
In order to take into account the effects of waves in the momentum equations, some quantities provided by wave models are required: period, significant wave height, direction, wavenumber, Stokes velocities, wavelength, τ aw the momentum flux from atmosphere to wave, and τ wo the momentum flux from wave to ocean linked to wave breaking. Some of them can be directly provided by the wave model, and others can be calculated from the available parameters, depending on the wave model chosen.
In the academic case studies, presented in Sect. 3, we use the SWAN wave model, and in the realistic simulation, presented in Sect. 4, we use the WW3 model, validated at global, regional and nearshore scales. These are third generation wave-averaged models that solve the two-dimensional wave action balance equations for wave action density as a function of (θ, σ ) for the SWAN model and a function of (θ, k) for the WW3 model. In Cartesian coordinates, this equation is written as Ocean Sci., 8, 657-681, 2012 www.ocean-sci.net/8/657/2012/ with N the wave action density N = E σ , c x , c y , c σ and c θ the propagation velocities in x-, y-, σ − and θ− space, respectively. The source/sink term, S tot on the right side, is expressed in terms of energy density and represents different physical processes available in the wave model: with S in the atmospheric source function, S nl4 the nonlinear quadruplet interactions and S ds,w the dissipation by whitecapping. Other phenomena induced by the finite depth effects like S nl3 triad nonlinear wave-wave interactions, S ds,b dissipation by bottom friction and S ds,br dissipation by depthinduced breaking are taken into account. Thus, diffraction, reflection, refraction and shoaling are included. SWAN (version 40.72) accounts for all these processes. It is generally used for wave transformation at nearshore and coastal scales (Booij et al., 1999;Dufois, 2008;Rusu and Soares, 2009;Bruneau, 2009). We will use this model for the academic cases.
WW3 has been widely used at global and regional scales and its validity is now extended to nearshore scales (version 4.04) with parameterizations of wave breaking, bottom dissipation and wave dissipation, avoiding the use of a specific nearshore wave model. One can find more information about the parameterizations proposed by this version in Ardhuin et al. (2010). This model has been validated using in-situ and remote sensing data (Ardhuin et al., 2008aDelpey et al., 2010).

A normal plane beach test case
This test case consists of obliquely incident spectral waves approaching an idealized smooth plane beach. It was initially posed by Haas and Warner (2009), hereinafter named HW09 and more recently by Uchiyama et al. (2010), hereinafter called UMS10. HW09 compared two hydrodynamic models: the quasi-3-D model SHORECIRC  and the 3-D model ROMS (Shchepetkin and McWilliams, 2005), where wave forcing followed the depth-dependent radiation stress formalism of Mellor (2003). UMS10 compared these solutions with another version of ROMS where a vortex force  approach is used. All these solutions were forced rigorously by the same wave field simulated by SWAN. We suggest comparing our model to this test case to assess its validity and performance to those of previous models.
The bathymetry is a plane beach with a constant slope of 1 : 80. It has realistic dimensions (1180 m in the cross-shore direction x and 1200 m in the alongshore direction y). The coast is oriented to the west and the offshore boundary is set at x = 0 with the maximum water depth (12 m). We use the same grid spacing of 20 m in horizontal directions as in previous simulations, and 10 vertical levels. To be consistent with UMS10 and HW09, a quadratic bottom drag law (Eq. 21) is used with a drag coefficient C d sets constant equal to 0.0015.
At the offshore boundary, a Jonswap type spectral modeled wave field is imposed with a 2 m significant wave height, 10 s peak period and an incidence angle of 10 • . We use the same wave fields as HW09 and UMS10. We also neglect the roller effects and the bottom streaming. Earth rotation is excluded with the Coriolis parameter set to 0. Lateral periodic conditions are used. As a first step, we do not take into account the influence of waves on vertical mixing and on the surface roughness length. UMS10 conducted four simulations: a 2-D barotropic case (Run a) and three 3-D cases where the vertical profiles of the vector of breaking dissipation or the vertical mixing are changed (Runs b, c and d). In Run b, the vertical penetration of momentum associated with breaking waves is concentrated near the surface, whereas in Runs c and d, penetration is quite uniform along the vertical column. We impose for this test that our vertical profile for the momentum associated with breaking waves is similar to the Run b.
UMS10 also calculated an analytical solution for the barotropic velocities and the surface elevation.

Reference simulation
Waves begin to break between 500 < x < 1000 m, (as shown by the breaking dissipation rate Fig. 16a in UMS10), and the significant wave height decreases for x = 600 m. A slight setdown before the breaking point and a set-up reaching 22 cm at the shoreline are observed. After two hours of simulation time, our simulation becomes stationary. Our surface elevation agrees with both the analytical and numerical results of UMS10 (Fig. 1a). The cross-shore barotropic velocity (Fig. 1b) is the same as that of UMS10 (and equal to the depth-averaged anti-Stokes flow because of the mass balance) and the alongshore barotropic velocity (Fig. 1c)  of K z is quite different close to the surface in the surf zone, from the profile of K z in Run b, because of our difference in the turbulence closure model as pointed out above.
We will display the depth-averaged terms of the momentum equations to describe in detail how the different forces are balanced (Fig. 3)  The alongshore momentum balance is between four forces: the breaking acceleration drives the northward alongshore velocities, the southward bottom stress, the vortex force and the Eulerian advection. The vortex force is totally compensated by the Eulerian advection. Until x reaches 690 m, the vortex force is oriented southward and thereafter, it becomes positive and is oriented northward. The advection force has the same pattern but in an opposite sign. The depthaveraged cross-shore accelerations are one order of magnitude larger than the alongshore ones. Among the most important forces is the breaking acceleration, which is larger in the surf zone, directing the surface cross-shore velocity to the shoreline. This momentum input is competed with the pressure gradient force, which is negative in the surf zone. The depth uniform term of the Bernoulli's head gradient − ∂S J ∂x is less important. When waves shoal before the surf zone, this force is negative and balances with the positive pressure gradient force, creating a set-down. In the surf zone, it turns positive, whereas the pressure gradient force is negative and more important. A set-up is generated. Finally, the shearinduced term of the Bernoulli's head gradient − ∂S shear ∂x is Ocean Sci., 8, 657-681, 2012 www.ocean-sci.net/8/657/2012/ negligible except in the offshore part of the surf zone, where it is positive. As a consequence of this term, the transition from set-down to set-up is displaced further offshore and the slope of the surface is reduced. Rascle (2007) has already noticed this in his test case. It thus appears that the vertical profile of the velocities is widely dependent on the vertical profile of the breaking acceleration, on the vertical mixing, as well as the vertical profile of the vortex force which is related to the one of the Stokes velocities.

Sensitivity tests considering the surface conditions and surface roughness
As we have seen before, vertical shear is highly related to the vertical mixing. In this section, we test different surface boundary conditions in the parameterization of the turbulence closure. In the reference simulation, eddy viscosity was parameterized according to the k − turbulent closure scheme. At the surface, we test the addition of the condition of Craig and Banner (1994) (Eq. 34) and also different surface roughnesses z surf (Eq. 36). In the reference simulation, the surface roughness was fixed to 0.015 m. We therefore test the values 0.8 H s and 2.4 H s as in Rascle et al. (2006). The parameterization of Craig and Banner (1994) adds a flux of energy that slightly increases the vertical mixing only near the surface (Fig. 4a). The vertical profile of the crossshore velocity and the depth-integrated alongshore velocity are very close to the ones of the reference simulation ( Fig. 4b  and c). When the surface roughness is increased and related to the significant wave height, the vertical shear is decreased and the velocities are more depth-uniform. The alongshore velocities are thus increased and the peak is moved offshore.
In conclusion, our results agree with the previous simulations performed by other models using different theories. The littoral drift and vertical profiles are correctly reproduced by our model. Nevertheless, the sensitivity tests and the analysis show that these profiles are highly dependent on the vertical mixing and the vertical penetration of the breaking acceleration force and the Stokes velocities. Even if the model is in agreement with the others, a comparison with in-situ data or laboratory measurements is necessary to assess whether the 3-D characteristics are accurate (Sects. 3.2 and 4).  (a), the cross-shore velocities (b) and the cross-shore profile of the longshore depth-integrated velocities (c) with or without the parameterization of Craig and Banner (1994)) and different surface roughnesses in the turbulence closure. "ref" is the result for the reference simulation.

A barred beach with rip current
The purpose of this test case is to check the ability of the model to correctly reproduce the rip current phenomena, before tackling the study of the complex sandbar systems of the Têt inner shelf. We reproduce test B experiments of Haller et al. (2002) performed in the basin of the Ocean Engineering Laboratory (University of Delaware). Previous modelers have reproduced this experiment with the SHORECIRC model (Haas et al., 2003), with MARS (Bruneau, 2009), and also with ROMS using the wave forcing radiation stress approach of Haas and Warner (2009). The size of the modeled basin is 15.8 m in the cross-shore direction x and 18.6 m in the alongshore direction y. Between 1.5 m and 3 m from the wave maker, the beach slope is steep (1 : 5) but it is mild (1 : 30) for the rest of the domain. A longshore bar system made up of three bars of 7.32 m in length and 6 cm in height, separated by rip channels of 1.82 m, is located at 6 m from the coast. The grid spacing is similar to Haas et al. (2003), and is 20 cm in the horizontal direction. Seven vertical levels are used. We set wall conditions at the borders. Rotation is excluded. On the bottom, we use the Eq. (21) for the boundary condition. The wave forcing is performed by SWAN. A monochromatic wave is imposed at the offshore edge with a 0.0724 m significant wave height and 1 s peak period, perpendicular to the direction of the beach (Fig. 5a). Waves break suddenly over the bar while being more progressive through the rip channel. There is a little shoaling before the breaking point above the bar, but it is insignificant through the rip channel. As previously noted by Haas et al. (2003) and Weir et al. (2011), this is because there is no forcing by the current on waves.
The cross-shore profiles of the surface elevation are different over the bar and through the channel (Fig. 5b) and consistent with Haas and Warner (2009). In fact, there is a set-up over the bar and another one near the shoreline (corresponding to the locations where waves break), whereas in the channel, the set-up is more progressive.
In previous simulations and in the experiment, two recirculation cells of current are generated by the wave forcing ( Fig. 6): one in the surf zone with currents oriented shoreward over the bars and offshore above the channels, and another less marked between the bars and the shoreline. This second recirculation cell is made up of the excess of water brought by waves waiting to be evacuated offshore via the channels. Moreover, previous numerical simulations (e.g. Yu and Slinn, 2003;Haas et al., 2003;Haas and Warner, 2009 With a quadratic bottom drag law with a drag coefficient equal to 0.015, we obtain two rip currents and also two little recirculation cells between the bars and the shoreline. However, these circulations become quickly stationary. Decreasing this coefficient to 0.00015 and the horizontal diffusion through the Smagorinsky coefficient (divided by 5), we observe that the two rip currents continually oscillate and meander to the left and right of the channel. A time-average of the circulation during 30 min, once the oscillations are well settled, was calculated to compare with the time-averaged measurements. We observe that the rip current in the channel at the top of Fig. 6d is on average directed toward the left. This pattern is also noticed in the simulations of Haas and Warner (2009). The second time-averaged rip current is directed toward the right, but with a smaller angle. We can notice that the direction of this rip current differs from the previous simulations (Fig. 6). Different sensitivity tests have shown that this result is highly dependent on many processes, such as the bathymetry, the wave parameters, the bottom roughness, the boundary conditions, etc. Moreover, our rip currents are more extended offshore and are narrower in the channel than in the other simulations. The effects of the currents on waves not represented here, could be responsible of these differences (Smith, 2006), as Yu and Slinn (2003); Weir et al. (2011) have shown that rip currents reduce the flux of momentum from waves to currents due to wave breaking. Finally, intensities of depth-integrated currents are comparable to the data, with a maximum value equal to 0.25 m s −1 .
The consistency of the current vertical profiles is hard to check because no 3-D-measurements were performed in the experiment. They are difficult to obtain in general because of the sporadic and changing nature of currents in these kinds of systems. We will thus compare our results to the ones of Haas and Warner (2009). They pointed out that above channels, cross-shore current are sheared, being stronger in the upper part of the water column. As the rip exits the channel, the velocity stays large near the surface and is weak close to the bottom, as Haas and Svendsen (2002) had observed. In our simulation (Fig. 5c), cross-shore velocities are stronger in the channel than above the bar. We obtain the same results in and downstream from the channel, with a larger current in the upper part of the water column. We also observe that above the bar, the current is stronger in the middle of the water column, and downstream, it is oriented onshore close to the bottom and is weak near the surface. In conclusion, this simulation gives reasonable results compared to the ones of Haas and Warner (2009  These two test cases deal with littoral scales, where the wave action is the most intense. The model is able to reproduce the wave-induced processes at this scale. In the next section, a simulation at the inner-shelf scale, is performed to test the validity of our model in a region where the circulation results from a wide range of processes.

Coastal circulation and the Têt system
The Têt River discharges into the southwestern part of the Gulf of Lion (hereinafter GoL) in the northwestern Mediterranean Sea (Fig. 7). Circulation in this micro-tidal zone in front of this river is strongly controlled by wind conditions. Estournel et al. (2003) and Ulses et al. (2008a) show that two major winds, the Tramontane (NW) and Marin (SE), induce cyclonic circulation in the GoL generating a southward current along the Têt coast. During east or south-east storms, this general counterclockwise circulation is intensified in the inner shelf but is opposed by an alongshore northward littoral drift. Evidence of this drift has been provided by Anguenot and Monaco (1967) with radioactive tracers, by Delpont and Motti (1994) with aerial and SPOT images, multi-date remote sensing by Certain (2002) and by Bourrin et al. (2008), who analyzed bathymetric and sediment data. The Têt is a small river with an average water discharge of less than 10 m 3 s −1 , with exceptional peaks two orders of magnitude higher during high precipitation events (Serrat et al., 2001). The littoral zone, where it discharges, has a complex bathymetry. A recent LiDAR (Light Detection and Ranging) survey (Fig. 8) shows a sand spit developing northwards of the breakwaters of the Canet-en-Roussillon harbor, followed by a deep pit, as observed by Bourrin et al. (2008). Thereafter, complex double crescentic sandbar systems, classified as low tide terraces (LTT) by Aleman et al. (2011), are observed. They appear chaotic and are likely disturbed by the breakwaters of the harbour and the river. The internal bars have their left side more onshore than their right side, suggesting that they have been modified by the northward littoral drift.

The storm of 21 February 2004
During the sampling period, two major storms occurred, one on 4 December 2003 and another on 21 February 2004 (Guillén et al., 2006). We focus on the second storm since more data are available for this period. The storm was characterized at SOPAT by a maximum significant wave height H s > 7 m and a peak period T > 12 s, with a westward peak direction (Guillén et al., 2006). At SODAT, significant wave height reached 6 m at 5 a.m. (Figs. 10 and 11) while the wind blew out of the south-east and reached up to 16 m s −1 (Fig. 9). The water and sediment discharge of the Têt, as es-timated by Guillén et al. (2006), were very low compared to the previous storm (only 450 t of sediment compared to 20 000 t in December). Wind and waves were thus the predominant forcing during this storm. According to Guizien (2009), the return period for both storms was 10.5 years at Sète (located 100 km to the north-east) and 5 years at Banyuls (located 20 km further south). Before and after the storm, the current was southward at SODAT (11 m) (Figs. 12, 13), with low intensity (<10 cm s −1 ). At the beginning of the day (21 February), the direction of the current at SODAT turned toward the north, and the current increased in intensity to reach approximately 90 cm s −1 throughout the water column at 4 a.m. of the same day. Its intensity remained high, but then began to decrease after 4 h, remaining at moderate intensity (around 20 cm s −1 ) for 30 h while the direction turned southward. At POEM (28 m) and SOPAT (31 m) (Figs. 12, 13), the current was generally oriented southward. During the storm, when wind strengthened, it increased reaching about 50 cm s −1 at the surface, and 40 cm s −1 near the bottom. At these two offshore stations, the current remained abnormally strong (>15 cm s −1 ) for more than 50 h.

Implementation and results
We aim to accurately reproduce phenomena induced by waves and current, covering scales from the whole western Mediterranean Sea to the Têt nearshore zone. A first attempt consisted of using four nested grids for the hydrodynamic circulation model (with grid resolutions between 2.5 km and 15 m). Using this set up, spurious flows were observed near the shoreline at the northern boundary of the finest grid. This was due to the representation of the littoral drift that strongly depends on the resolution of the model. As pointed out by Davies and Jones (1996), one solution is to use an unstructured grid or a grid with a variable resolution, that covers the entire Têt inner-shelf, with a fine resolution at the Têt mouth which is gradually reduced to a coarser resolution in offshore zones. Using such grids ensures a smooth transition between offshore and nearshore zones. We choose this second approach here.

Wave model implementation and results
We use three nested grids to model the sea state, two structured grids that cover the whole western Mediterranean Sea and the Gulf of Lion, respectively, and an unstructured grid, which runs from the inner-shelf with a resolution of 550 m at the offshore boundaries to the surf zone of the Têt (Fig. 7 and Table 1). The size of the cells is 22 m near the Têt mouth. The grid is made of 64 000 nodes and 127 500 elements.
Simulations are run with WW3 for a period of two months, from 4 February 2004 to 26 March 2004 (the period for which observations are available). We use the TEST405 parameterizations as described in Ardhuin et al. (2010), which are more adapted for the younger seas that occur in the Mediterranean Sea. The wind velocities are provided by the Aladin model (a regional weather forecasting model focused on France with a resolution of 10 km) from Météo-France every 3 h, except for WW3-MEDOC where Aladin is supplemented by Arpege (a global atmospheric model from Météo-France with a grid resolution of 15 km over France). Output wave spectra are discretized over 36 directions with 10 • of resolution and 30 frequencies f n spaced with the relation f n+1 = 1.1f n from 0.05 Hz to 0.8 Hz. Bathymetry in the Têt surf zone is complex and the length of sand bars ranges between 200-300 m. To correctly reproduce the wave breaking, and consequently the wave-induced current, it is necessary to simulate the waves with a resolution coherent with the size of the bars. A resolution of 22 m is used near the Têt in this study.
We compare the wave model results to the significant wave heights and wave periods recorded by the two wave gauges Statistical results show a good agreement between the two datasets and the simulation (Table 2). For example, a correlation of 93 % is found for the significant wave height at POEM and 96 % at SODAT. During the storm period (Figs. 10 and 11), the three parameters fit well. We note, however, that significant wave heights are underestimated by the model (Fig. 11), with a bias of 20 cm at the storm apex and especially during the afternoon, with a bias of 1.5 m. We suspect this discrepancy to be linked to a wrong estimation of the wind. A comparison between wind intensity measured at the Toreilles meteorological station and the one simulated at SO-DAT ( Fig. 9) shows that the Aladin model seems to be in reasonable agreement with the data. However, this result is not convincing as the wind over the sea is expected to be stronger than the wind on land. A sensitivity test with another atmospheric model has been performed, showing that depending on the models, wave heights can be under or overestimated. Satellite wind data have been examined to find evidences of this underestimation but the absence of valid data near the coast did not allow to draw a conclusion. Finally, observations and simulations both indicate that significant wave height decays between the two sites, suggesting that wave dissipation occurs in the inner shelf zone.

Current model implementation
As for the wave model, three nested grids for the circulation model are deployed, with the focus towards the Têt nearshore (Fig. 7). All details concerning the different grids are presented in Table 3. Grid TET is a stretched curvilinear horizontal grid with a variable horizontal resolution (Madec, 2008), from 8 m × 8 m at the nearest grid point from the Têt mouth to 180 m × 180 m at the external border. Bathymetries from Berné et al. (2002) and from the LiDAR for the nearshore are used. The last one has a resolution of 5 m that is consistent with a grid resolution of 8 m. As explained above, a high resolution near the river mouth is necessary in order to reproduce all current patterns generated by the crescentic sandbars that impacted the SODAT instrument. Daily river discharges were provided by Banque Hydro and Compagnie Nationale du Rhone (http://www.hydro.eaufrance.fr/). The meteorological forcings (surface pressure, air temperature, relative humidity, wind velocity and radiative fluxes) are taken from the Aladin model every 3 h. The regional circulation model (grid MEDOC) is initialized and forced every day by the large-scale Ocean General Circulation Model (OGCM, Tonani et al., 2008). The wave forcing is not been taken into account in the circulation model at the regional scale (MEDOC) but at all other scales, every 3 h for the GoL, and every 1 h for TET. The roughness length is set to 1 cm throughout the domain.
Ocean Sci., 8, 657-681, 2012 www.ocean-sci.net/8/657/2012/  12. Comparison of the current intensity near the bottom (right) and close to the surface (left) at the three instruments, between the measured current (black) and the simulated current with (with WEC -red) and without the wave forcing (without WEC -blue).

Importance of the wave forcing
Firstly, a simulation without wave forcing (Figs. 12, 13) is performed. All other forcing terms are present, including the wind and the larger scale circulation. Simulated currents are very small, and neither littoral drift nor rip currents are observed. At SODAT (11 m), current intensity is 0.17 m s −1 close to the surface and near the bottom. It is directed southward throughout the water column. This value is not consistent with the measured values. At POEM (28 m) and SOPAT (31 m), current intensities increase similarly to the measure-ments in the first hours of the storm but decrease too soon after the apex.

Current in the surf zone
At the beginning of the storm, waves propagate from the east, with the irregularities of the bathymetry creating alongshore variations in breaking wave heights (Fig. 14), which in turn are responsible for the complex recirculation cells (Fig. 15, top) and oscillating meanders in the surf zone (Bowen, 1969). These types of meanders are often observed (e.g. Reniers et al., 2001)  made of bars and channels, the current is dominated by a rip-current flow and not a longshore drift. This is what we observe here (sections 1 and 2 in Fig. 15). A bar, at a depth of 2.5 m in section 1, is able to break waves and generates a strong feeder current that circulates through to the beach and exits offshore near section 2, where the breaking bar is too close to the shore and waves are already broken. The rip current dynamics are more complex than in the test case, and largely influenced by the Canet harbour tip. Vertical sections show that over the breaking bar (Fig. 16a), crossshore velocities are stronger close to the surface (>0.8 m s −1 ) and directed onshore almost everywhere. Above the channel (Fig. 16b), current is oriented seaward everywhere, and is stronger in the middle of the water column.
On 21 February around 2 a.m., the incident direction turns and breaking waves arrive at the coast obliquely from the east-south-east. They create a longshore northward drift almost everywhere between 5 m and 12 m of water depth (Fig. 15, bottom). Some recirculation cells are still observed in the surf zone and pertubate the path of the northward drift. Locally, the Canet harbour in the south of the domain disturbs the longshore drift by shifting it offshore. These results are consistent with the development of sand spits growing northwards at river mouths, sand bars and harbours, as observed by Delpont and Motti (1994); Bourrin et al. (2008) and in the LiDAR bathymetry. Downstream from the harbour, since the beginning of the storm, a cyclonic eddy and a return current along the northern breakwater of the harbour are generated, as discussed by Trampenau et al. (2004  These phenomena explain the strong erosion observed here and in general along the harbour side in the lee of the waves (Trampenau et al., 2004). The direction of the drift is not as simple as in the first test case, and vertical profiles of current are difficult to analyse because recirculation, rip currents and littoral drift are present and interact. In the drift, however, current is quite uniform with depth (Fig. 13, top and Fig. 17). We compare the simulated currents with the measured currents at SODAT (Figs. 12, 13). Regarding the rotation of the measured current during the rising stage of the storm (Fig. 13), it seems that the instrument is first in a zone perturbated by a rip current and then in the longshore drift. This is confirmed by the simulation. In fact, in the first hours of the storm (between 11 p.m. on 20 February and 4 a.m. on 21 February), the intensity of the rip current reaches 55 cm s −1 close to the surface (at 1 a.m. on 21 February) (Figs. 13 and 15) and with a direction globally toward the east. The measurement shows an eastward current with weaker intensities (<20 cm s −1 ). The accuracy of the rip current simulation is highly dependent on the bathymetry. The modelled bathymetry, built from the LiDAR survey con-Cross-shore velocity u (m/s) Cross-shore distance (km) Depth (m) Fig. 17. Vertical section (corresponding to section 3 in Fig. 15) of the cross-shore velocity at the storm apex. ducted in 2008, four years after this studied storm, could be the cause of the discrepancy between simulated and observed current. At 6 a.m. i.e. the storm peak, the simulated current turns towards the east-north-east. Intensity of the current is well reproduced (it reaches 90 cm s −1 at the surface and 45 cm s −1 at the bottom) with a delay of 2 h, while the simulated direction is rather eastward than northward at least near the surface. Locally, a recirculation cell strongly influences the circulation, but outside this local pattern, the northward drift dominates the circulation (Fig. 15). In addition, near the bottom, the simulated current is underestimated (45 cm s −1 against 85 cm s −1 in the reality) (Fig. 12). The misrepresentation of the bottom roughness could contribute to increase the error on the roughness of the model, which largely influences the bottom current. A sensitivity test was performed and proved that the drift intensity in the entire water column was increased when the roughness was decreased.

Current on the inner shelf
On the whole inner shelf (depth >25 m), simulated currents are southward during the entire period and are intensified during the storm.
At POEM (28 m) and SOPAT (31 m), the simulated currents with and without the wave forcing have quite the same intensities (middle of Figs. 12, 13). Theses results show that waves have little effect on current at this scale. In the first hours of the storm, simulated current fits the data, but underestimates them thereafter. The discrepancy between model results and observation may be explained by an underestimation of the wind speed, as it has been suggested in Sect. 4.2.1. A test where we increase the wind speed by a factor of 1.2 in the circulation model shows that in the surf zone results are unchanged, but on the inner shelf, current intensities reach the observed values at the surface, and are increased in the entire water www.ocean-sci.net/8/657/2012/ Ocean Sci., 8, 657-681, 2012 column. This sensitivity test reveals that circulation on the inner shelf is highly dependent on the atmospheric forcing and the regional circulation, whereas in the surf zone, processes linked to waves are the most important. Moreover, it also shows that either the atmospheric model underestimates wind speed over the sea during storms, or the calculation of the surface wind stress is not adapted when the sea roughness is increased by waves. This discrepancy in current between model and observation during a storm at coastal scales is the focus of the study of Michaud et al. (2012).

Conclusions
We have developed and implemented a new method to take into account the impact of waves on the 3-D circulation. This method can be used from the nearshore to the global scale. It is first tested on two classical academic cases. Results fit with previous simulations performed by other models and with available observational data. A realistic case was then simulated of energetic waves arriving at a coast of the northwest Mediterranean for which currents were available at different depths as well as an accurate bathymetric database of the 0-10 m depth range. A grid nesting approach was used to account for the different forcings acting at different spatial scales. The simulation coupling the effects of waves and currents is successful to reproduce the powerful northward littoral drift in the 0-15 m depth zone, while without waves, the current is slow in the opposite direction. More precisely, two distinct cases were identified: when waves have a normal angle of incidence with the coast, they are responsible for complex circulation cells and rip currents in the surf zone, and when they travel obliquely, they generate a northward littoral drift. These features are more complicated than in the test cases, due to the complex bathymetry and the consideration of wind and non-stationary processes. Wave impacts in the inner shelf are less visible since wind and regional circulation seem to be the predominant forcings. In addition, the discrepancy between model and observations is noted at that scale, possibly linked to the underestimation of the wind stress. A perspective of this study could be to fully couple wave and circulation models. This means to also take into account the effects of current on waves in the wave model. This future work will allow to estimate the potential effect of current on wave properties through blocking or refraction itself impacting the circulation through modification of the water level.
Lastly, during storm events, a classical sediment transport approach without wave forcing (e.g. Ulses et al., 2008b) does not permit the reproduction of either the northward transport in the surf zone or the transport of large amount of fine particles discharged most of the time during events combining high waves and floods. Moreover, the bottom shear stress would be strongly underestimated and then, also the possibility of resuspension for coarse sediment. In the specific case of the region studied here, we expect to extend the study of Ulses et al. (2008b) on the impact of storms on the sediment transport at regional scale to the nearshore zones. We will be then able to study the fate of sediments ranging from the river and the beach to the open ocean and so complete the study undertaken by Palanques et al. (2011).

Appendix A
Equations of We come back to the Eqs. (6) and (7).