Development of a 2-way coupled ocean-wave model: assessment on a global NEMO(v3.6)-WW3(v6.02) coupled configuration

This paper describes the implementation of a coupling between a three-dimensional ocean general circulation model (NEMO) and a wave model (WW3) to represent the interactions of the upper oceanic ﬂow dynamics with surface waves. The focus is on the impact of such coupling on upper-ocean properties (temperature and currents) and mixed-layer depths (MLD) at global eddying scales. A generic coupling interface has been developed and the NEMO governing equations and boundary conditions have been adapted to include wave-induced terms following the approach of McWilliams et al. (2004) and Ardhuin et al. (2008). In particular, the contributions of Stokes-Coriolis, Vortex and surface pressure forces have been implemented on top of the necessary modiﬁcations of the tracer/continuity equation and turbulent closure scheme (a 1-equation TKE closure here). To assess the new developments, we perform a set of sensitivity experiments with a global oceanic conﬁgura-tion at 1 / 4 o resolution coupled with a wave model conﬁgured at 1 / 2 o resolution. Numerical simulations show a global increase of wind-stress

taken into account in ocean-atmosphere coupled models. For example, the momentum flux through the air-sea interface has traditionally been parameterized using the near surface winds (typically at 10 meter) and the atmospheric surface layer stability (Fairall et al., 2003;Large and Yeager, 2009;Brodeau et al., 2016). The physics of the coupling depends on the kinematics and dynamics of the wave field. This includes a wide range of processes from wind-wave growth, nonlinear wave-wave interaction, 25 wave-current interaction to wave dissipation. Such complex processes can only be adequately represented by a wave model.
Besides affecting the air-sea fluxes, waves define the mixing in the oceanic surface boundary layer (OSBL) via breaking and Langmuir turbulence. For example, Belcher et al. (2012) showed that Langmuir turbulence should be important over wide areas of the global ocean and more particularly in the Southern ocean. In this region, they show that the inclusion of the effect of surface waves on the upper-ocean mixing during summertime allows for a reduction of systematic biases in the OSBL depth. 30 Indeed their large eddy simulations (LES) suggest that under certain circumstances wave forcing can lead to large changes in the mixing profile throughout the OSBL and in the entrainment flux at the base of the OSBL. They concluded that wave forcing is always important when compared to buoyancy forcing, even in winter. Moreover, Polonichko (1997) and Van Roekel et al. (2012) emphasized the fact that the Langmuir cells intensity strongly depends on the alignment between the Stokes drift and wind direction. Langmuir turbulence is maximum when wind and waves are aligned and becomes weaker as the misalignment 35 becomes larger. Li et al. (2017) highlighted that ignoring the alignment of wind and waves (i.e. assuming that wind and waves are systematically aligned) in the Langmuir cells parameterizations leads to excessive mixing particularly in winter.
Most previous studies of the impact of ocean-wave interactions at global scale have been using an offline one-way coupling and included only parts of the wave-induced terms in the oceanic model governing equations (e.g. Breivik et al., 2015;Law Chune and Aouf, 2018). In this study, the objective is to introduce a new online two-way coupled ocean wave modeling 40 system with a great flexibility to be relevant for a large range of applications from climate modeling to regional short-term process studies. This modeling system is based on the Nucleus for European Modelling of the Ocean (NEMO, Madec, 2012) as the oceanic compartment and WAVEWATCH III® (hereinafter WW3, The WAVEWATCH III ® Development Group, 2016) as the surface wave component. NEMO and WW3 are coupled using the OASIS Model Coupling Toolkit (OASIS3-MCT, Valcke, 2012;Craig et al., 2017) which is widely used in the climate and operational community. The various steps for our 45 implementation are the following (i) inclusion of all wave-induced terms in NEMO, only neglecting the terms relevant for the surf zone which is outside the scope here (ii) modification of the NEMO subgrid scale physics (including the bulk formulation) to include wave effects and a parameterization for Langmuir turbulence (iii) development of the OASIS interface within NEMO and WW3 for the exchange of data between the models (iv) test of the implementation based on a realistic global configuration at 1/4 o for the ocean and 1/2 o for the waves. 50 To go into the details of those different steps, the paper is organized as follows. The modifications brought to the oceanic model primitive equations, their boundary conditions, and the subgrid scale physics to account for wave-ocean interactions are described in Sec. 2. This includes the addition of the Stokes-Coriolis force, the Vortex force, and the wave-induced pressure gradient. In Sec. 3 our modeling system coupling the NEMO oceanic model and the WW3 wave model via the OASIS3-MCT coupler is described in details. Numerical simulations are presented in Sec. 4 using a global configuration at 1/4 o for the oceanic 55 model and 1/2 o for the wave model. Using sensitivity runs, we assess those global configurations with particular emphasis on the impact of wave-ocean interaction on mixed-layer depth, sea-surface temperature and currents, turbulent kinetic energy (TKE) injection, and kinetic energy. Finally, in Sec. 5, we summarize our findings and provide overall comments on the impact of two-way ocean-wave coupling in global configurations at eddy-permitting resolution.
2 Inclusion of wave-induced terms in the oceanic model NEMO 60 In order to set the necessary notations, we start by introducing the classical primitive equations solved by the NEMO ocean model. Note that between the two possible options to formulate the momentum equations, namely the so-called "vector invariant" and "flux" forms, we present here the first one which will be used for the numerical simulations in Sec. 4. With u h = (u, v) the horizontal velocity vector, ω the dia-level velocity component, θ the potential temperature, ρ the density, ζ the relative vorticity, p h the hydrostatic pressure, p s the surface pressure, the Reynolds-averaged equations (with · the averaging operator, omitted here for simplicity) are Here k is a non-dimensional vertical coordinate, lateral derivatives ∂ x and ∂ y have to be considered along the model coordinate, and e 3 is the vertical scale factor given by e 3 = ∂ k z, where z is the local depth and ρ is given by an equation of state (Roquet et al., 2015). The necessary boundary conditions include a kinematic surface and bottom boundary condition which can be expressed in terms of the vertical velocity w with η the height of the sea-surface and (E − P ) the mass flux across the sea-surface due to precipitations and evaporation, a momentum surface boundary condition for the Reynolds stress vertical terms with τ oce = (τ oce u , τ oce v ) a wind-stress vector which represents the part of the stress that drives the ocean, and a dynamic 80 boundary condition on the free surface leading to the continuity of pressure across the air-sea interface. The kinematic boundary conditions (6) for w(z = η) and w(z = −H) translate into ω(z = η) = (E − P ) and ω(z = −H) = 0. We do not include explicitly here the boundary conditions for the tracer equations since they are unchanged from classical primitive equations models in the presence of wave motions. As mentioned earlier, in equations (1) to (5) prognostic variables have to be interpreted in an Eulerian-mean sense even if the averaging operator is not explicitly included.

Modification of governing equations and boundary conditions
Asymptotic expansions of the wave effects based on Eulerian velocities (McWilliams et al., 2004) or Lagrangian mean equations  lead to the same self-consistent set of equations for weak vertical current shears. These are further applied and discussed by Uchiyama et al. (2010), Bennis et al. (2011), Michaud et al. (2012, or Moghimi et al. (2013). The 3-component Stokes drift vector is u s = ( u s , v s , ω s ), and is non-divergent at lowest order (Ardhuin et al., , 2017b where wave-induced terms are represented with tildes. The F u and F v terms represent the sink/source of wave-momentum 100 due to breaking, bottom friction and wave-turbulence interaction. These terms will be neglected since they are expected to play a significant role only in the surf zone. The other extra contributions to the momentum equations include the Stokes-Coriolis force W St−Cor , the vortex force W VF , and a wave-induced pressure W Prs where the terms involving horizontal derivatives of ω have been neglected in W VF . In W Prs , the p J term corresponds to a depth uniform wave-induced kinematic pressure term 1 , while p Shear is a shear-induced three-dimensional pressure term 2 associated with the vertical component of the vortex force. The vortex force contribution W VF can be further simplified by neglecting the terms involving the vertical shear. In particular, the vertical component of the vortex force is absorbed in a 110 1 In the notations of Ardhuin et al. (2008) this term corresponds to p J = ρ 0 S J 2 In the notations of Ardhuin et al. (2008) this term corresponds to p Shear = ρ 0 S Shear pressure term p Shear (that gives the S shear term in the notations of Ardhuin et al., 2008). That particular term was neglected in Bennis et al. (2011) because of the general weak vertical shears in the wave-mixed layer. The effect of that term was also found to be much weaker than p J in shallow coastal environments, except in the surf zone. This assumption has the advantage to leave the hydrostatic relation (11) unchanged. Our implementation of wave-induced terms in NEMO is inline with Bennis et al. (2011) and corresponds to the simplified form of (12) Because of geostrophy, it is obvious that the addition of the Stokes-Coriolis force requires the effect of the Stokes drift on the mass and tracers advection to be taken into account. Regarding the joint modification of the tracers and continuity equations, it is clear that constancy preservation is maintained (i.e. a constant tracer field should remain constant during the advective transport) and that an additional wave related forcing must be added to the barotropic mode. The NEMO barotropic mode has 120 been modified accordingly since the surface kinematic boundary condition (6) in terms of vertical velocities w and associated w s now reads to express the fact that there is a source of mass at the surface that compensates the convergence of the Stokes drift, hence the barotropic mode is is the usual NEMO forcing term containing coupling terms from the baroclinic mode as well as slowly varying barotropic terms (including nonlinear advective terms) held constant during the barotropic integration to gain efficiency. In (13), the G x and G y contain the additional waveinduced barotropic forcing terms corresponding to the vertical integral of the W St−Cor and W VF terms which are also held 130 constant during the barotropic integration. A thorough analysis on the impact of the additional wave-induced terms on energy transfers within an oceanic model can be found in Suzuki and Fox-Kemper (2016). Note however that the study of Suzuki and Fox-Kemper (2016) is based on the Craik-Leibovich equations which are a special case of the more general wave-averaged primitive equations. Those sets of equations are equivalent to each other only at lowest order in vertical shear.

135
Reconstructing the full Stokes drift profile u s in the ocean circulation model would require obtaining the surface spectra of the Stokes drift from the wave model. Instead, profiles are generally reconstructed considering a few important parameters, including the Stokes drift surface value u s h (η) and the norm of the Stokes volume transport T s . In Breivik et al. (2014) and Breivik et al. (2016), Stokes drift velocity profiles are derived under the deep-water approximation in the general form u s h (z) = u s h (η)S(z, k e ) with k e a depth-independent spatial wavenumber chosen such that the norm of the depth integrated 140 Stokes transport (assuming an ocean of infinite depth) is equal to T s . The functions S B14 (z, k e ) from Breivik et al. (2014) and S B16 (z, k e ) from Breivik et al. (2016) for z ∈ [−H, η] are given by with erfc the complementary error function. It can be easily shown that for an ocean of infinite depth, the vertical integral of those functions are respectively equal to 1 6ke for S B16 and 1.34089 8ke ≈ 1 5.97ke for S B14 . Standard computations of Stokes drift 145 in numerical models are done in a finite difference sense, however due to the fast decay of u s h (z) with depth, a finite volume approach seems more adequate, in this case Such finite-volume interpretation of the Stokes drift velocity can also be found in Li et al. (2017) and Wu et al. (2019). The S B16 function is more adapted for this kind of approach since the primitive function does only require special functions available in 150 the fortran standard Since NEMO is discretized on an Arakawa C-grid, the components of the Stokes drift velocity must be evaluated at cell interfaces, a simple average weighted by layer thicknesses is used : Note that no explicit computation of the vertical component of the Stokes drift is necessary since in (7)-(11) ω s only appears summed with ω such that the relevant variable is ω + ω s as a whole. This quantity is diagnosed from the continuity equation (10) where the temporal evolution of vertical scale factors ∂ t e 3 is given by the free-surface evolution when a quasi-Eulerian vertical coordinate is used (e.g. z or terrain-following coordinates).
As illustrated in Fig. 1, for the typical vertical resolution used in most global models the properties of the discretized Stokes 160 profiles can be very different from their continuous counterparts. Indeed, the S B16 (z, k e ) function has been considered superior to S B14 (z, k e ) because the vertical shear near the surface is expected to be better reproduced. However in Fig. 1 it is shown that this is no longer the case at a discrete level since the discrete vertical gradients at one meter depth turns out to be larger for S B14 (z, k e ) compared to S B16 (z, k e ). In this case, the fast variations of S B16 (z, k e ) near the surface can not be represented by the computational vertical grid. A vertical resolution finer than the one currently used in most global ocean models near the Under the assumption of horizontal homogeneity, generally retained in general circulation models, the contribution from Stokes drift to the Turbulent Kinetic Energy (TKE) prognostic equation arises from the Vortex force vertical term W z VF = u s ∂ z u + 170 v s ∂ z v in the hydrostatic relation (11). Mimicking the way the TKE equation is usually derived (see e.g. Tennekes and Lumley, 1972) and using an averaging operator · satisfying the "Reynolds properties", we find that the turbulent fluctuations, defined after multiplication by w and averaging we obtain where the last two terms in the right-hand-side cancel with similar terms appearing when forming the equations for u ∂ t u and v ∂ t v (see eqn (A.7) and (A.8) in Skyllingstad and Denbo (1995)). The extra terms associated with the Stokes drift in the horizontally homogeneous TKE equation are thus u s ∂ z u w and v s ∂ z v w which can be further rewritten as The first term will modify the shear production term, it can also be derived by taking the Lagrangian mean of the wave-resolved TKE equation (Ardhuin and Jenkins, 2006). The second will enter the TKE transport term which is usually parameterized as −K e ∂ z e. The prognostic equation for the turbulent kinetic energy e in NEMO under the assumption that K e = A vm , with A vm the eddy viscosity, is thus with A vt the turbulent diffusivity, N the local Brunt-Väisälä Frequency, l ε a dissipative length scale, and c ε a constant parameter (generally such that c ε ≈ 1/ √ 2). Once the value of e is know, eddy diffusivity/viscosity are given by with Prt the Prandtl number (see Sec. 10.1.3 in Madec (2012) for the detailed computation of Prt), l m a mixing length scale, and C m a constant.

190
In addition to the modification of the shear production term in the TKE equation, the wave will affect the surface boundary condition both for e, l m , and l ε . The Dirichlet boundary condition traditionally used in NEMO for the TKE variable is modified into a Neumann boundary condition meaning that the injection of TKE at the surface is given by the dissipation of the wave field via the wave-ocean S oce term,

195
which is a sink term in the wave model energy balance equation usually dominated by wave breaking, converted into an ocean turbulence source term. In practice, this sum of S oce is obtained as a residual of the source term integration, hence it also includes unresolved fluxes of energy to the high frequency tail of the wave model. Due to the placement at cell interfaces of the TKE variable on the computational grid, the TKE flux is not applied at the free-surface but at the center of the top-most grid cell (i.e. at z = z 1 ). This amounts to interpret the half grid cell at the top as a constant flux layer which is consistent with 200 the surface layer Monin-Obukhov theory.
The length scales l m and l ε are computed via two intermediate length scales l up and l dwn estimating respectively the maximum upward and downward displacement of a water parcel with a given initial kinetic energy. l up and l dwn are first initialized to the length scale proposed by Deardorff (1980), l up (z) = l dwn (z) = 2e(z)/N 2 (z). The resulting length scales are then limited not only by the distance to the surface and to the bottom but also by the distance to a strongly stratified portion of the 205 water column such as the thermocline. This limitation amounts to control the vertical gradients of l up (z) and l dwn (z) such that they are not larger that the variations of depth (Madec, 2012) ∂ k |l · | ≤ e 3 , l · = l up , l dwn with κ the von Karman constant and C m , c ε the constant parameters in the TKE closure. The surface roughness length z 0 can be directly estimated from the significant wave height provided by the wave model as z 0 = 1.6H s  eqn (5)) which provides a proxy for the scale of the breaking waves. Note that in our study, no explicit parameterization of the 215 mixing induced by near-inertial waves has been added (Rodgers et al., 2014). As highlighted by Breivik et al. (2015), without activating this ad hoc parameterization in the standard NEMO TKE scheme, the model does not mix deeply enough. They also speculated that this ad hoc mixing could mask effects of wave-related mixing processes such as Langmuir turbulence. For this reason, it is thus not used in the present simulations.

Langmuir turbulence parameterization 220
Langmuir mixing is parameterized following the approach of Axell (2002). This parameterization takes the form of an additional source term P LC in the TKE equation (14). P LC is defined as where w LC represents the vertical velocity profile associated with Langmuir cells and d LC their expected depth. Following Axell (2002), w LC and d LC are given by where u s LC is the portion of the surface Stokes drift contributing to Langmuir cells intensity and c LC a constant parameter. In the absence of information about the wave field it is generally assumed that u s LC ∝ τ . As mentioned in the introduction, Polonichko (1997) and Van Roekel et al. (2012) showed that the intensity of Langmuir cells is largely influenced by the angle between the Stokes drift and the wind direction. To reflect this dependency we account for this angle in our definition of u s with e τ the unit vector in the wind-stress direction. The difference between the surface Stokes drift u s (η) and u s LC given by (16) is shown in Fig. 2 and compared to the usual parameterization of u s (η) as 0.377 τ oce /ρ 0 in the uncoupled case (see Madec, 2012). The modulation of u s LC depending on the wind-stress orientation significantly reduces the input of   4.2.3, greatly improved the MLD distribution.

Modeling system and coupling strategy
Our coupled model is based on the NEMO oceanic model, the WW3 wave model, and the OASIS library for the data exchange and synchronization between both components.
3.1 Numerical models and coupling infrastructure 250 The ocean model : NEMO NEMO is a state-of-the-art primitive-equation, split-explicit, free-surface oceanic model whose equations are formulated both in the vector invariant and flux forms (see (1) for the vector invariant form). The equations are discretized using a generalized vertical coordinate featuring, among others, the z -coordinate with partial step bathymetry and the σ-coordinate as well as a mixture of both (Madec, 2012). For efficiency and accuracy in the representation of external gravity waves propagation, 255 model equations are split between a barotropic mode and a baroclinic mode to allow the possibility to adopt specific numerical treatments in each mode. The NEMO equations are spatially discretized on an Arakawa C-grid in the horizontal and a Lorenz grid in the vertical, and the time dimension is discretized using a Leapfrog scheme with a modified Robert-Asselin filter to damp the spurious numerical mode associated with Leapfrog (Leclair and Madec, 2009). For the current study the NEMO equations have been modified to include wave effects as described in (7) and (13). Moreover the modifications to the standard 260 NEMO 1-equation TKE closure scheme are given in Sec. 2.3.

The wave model : WW3
The NEMO ocean model has been coupled to the WW3 wave model. In numerical models, waves are generally described using several phase and amplitude parameters. We provide here only the few sufficient details to understand the coupling of waves with the oceanic model, an exhaustive description of WW3 is given by The WAVEWATCH III ® Development Group 265 (2016). WW3 integrates the wave action equation (Komen et al., 1994), with the spectral density of wave action N w (k w , θ w ), discretized in wavenumber k w and wave propagation direction θ w for the spectral space (subscripts w are used here to avoid confusion with previously introduced notations).
where λ is longitude, φ is latitude, and S is the net spectral source term that includes the sum of rate of change of the surface 270 elevation variance due to interactions with the atmosphere via wind-wave generation and swell dissipation (S atm ), nonlinear wave-wave interaction (S nl ), and interaction with the upper ocean that is generally dominated by wave breaking (S oce ). Those parameterized source terms are important in waves-ocean coupling. Indeed, as shown earlier in (15) the S oce term is used to compute the TKE flux transmitted to the ocean, and the S in term enters in the computation of the wave-supported stress. They are here computed following Ardhuin et al. (2010b). In (17), the dot variables correspond to a propagation speed given by where R is earth radius, u h (z = η) = ( u| z=η , v| z=η ) are the surface currents provided by the ocean model, c g is the group 280 velocity, ω w the absolute radian frequency, and H the mean water depth. Equation (17) is solved for each spectral component (k w , θ w ) coupled by the advection and source terms. Equations (18)-(21) show how the oceanic currents affect the advection of the wave action density, there are also indirect effects via the source term ).

The coupler : OASIS3-MCT
The practical coupling between NEMO and WW3 has been implemented using the OASIS3-MCT (Valcke, 2012;Craig et al., 285 2017) software primarily developed for use in multi-component climate models. This software provides the tools to couple various models at low implementation and performance overhead. In particular, thanks to MCT (Jacob et al., 2005), it includes the parallelization of the coupling communications and runtime grid interpolations. For efficiency, interpolations are formulated in the form of a matrix-vector multiplication where the matrix containing the mapping weights is computed offline once for all. In practice, after compiling OASIS3-MCT, the resulting library is linked to the component models so that they have access 290 to the specific interpolation and data exchange subroutines. Now that we have described the different components involved in our coupled system, we go into the details of the nature of the data exchanged between both models.

Oceanic surface momentum flux computation
Surface waves affect the momentum exchange between the ocean and the atmosphere in two different ways. First the modification of surface roughness acts on the incoming atmospheric momentum flux τ atm . Second, a part of the momentum flux from 295 the atmosphere is consumed by the wave field and contributes to the growing waves (the so-called wave-supported stress) and conversely the waves release momentum to the ocean when they break and dissipate. This implies that the wind-stress transferred to the oceanic model (we call it τ oce ) is different from the atmospheric wind-stress τ atm . These two coupling processes are taken into account in our coupled framework.
The 10 meters wind u atm 10 is sent to the wave model which internally computes the dimensionless Charnock parameter α ch 300 characterizing the sea surface roughness (Charnock, 1955;Janssen, 2009). Those informations are used by the wave model to compute compute its own atmospheric wind-stress τ atm ww3 assuming neutral stratification, i.e. τ atm ww3 = ρ a C DN (α ch ) u atm 10 u atm 10 with C DN a neutral drag coefficient which is function of the Charnock parameter. Then the wave model computes the momentum flux transferred to the ocean τ oce ww3 . Using the latest available values of α ch , τ atm ww3 , τ oce ww3 , and u atm 10 , the oceanic model computes an atmospheric wind-stress τ atm using its own bulk formulation and the local value of the momentum flux going 305 into the water column is diagnosed as where the τ ww3 quantities are interpolated from the wave grid to the oceanic grid. In NEMO, the wind-stress is computed using the IFS 3 bulk formulation such as implemented in the AeroBulk 4 package (Brodeau et al., 2016). In particular the roughness length which enters in the definition of the drag coefficient is function of the Charnock parameter α ch where α m = 0.11, u is the friction velocity, and ν the air kinematic viscosity whose contribution is significant only asymptotically at very low wind speed. Note that in the uncoupled case the default value of the Charnock parameter is α 0 ch = 0.018. In our implementation, the momentum fluxes are computed using the absolute wind u atm 10 at 10m rather than the relative wind u atm 10 − u h (z = η). Indeed, several recent studies have emphasized that the use of relative winds is relevant only when a full 315 coupling with an atmospheric model is available since in a forced mode it leads to an unrealistically large loss of oceanic eddy kinetic energy (e.g. Renault et al., 2016). This is not a limitation of our approach since a simple modification of a namelist parameter allows to run with relative winds but this case is not investigated in the present study.
In our coupling strategy two different values of the atmospheric wind-stress and of the wave to ocean wind stress are computed with two different bulk formulations. This strategy is not fully satisfactory since it breaks the momentum conservation. and not on its orientation. Instead of (22), the atmospheric wind stress was corrected as follow: However, this approach potentially leads to artificially large values of τ oce when τ atm ww3 is small and it does not take into account the slight change in τ oce direction induced by the waves.

Additional details about the practical implementation
In Table 1 the different variables exchanged between the oceanic and wave models are given. All variables are 2D variables 330 meaning that no 3D arrays are exchanged through the coupler. All 2D interpolation are made through a distance weighted bilinear interpolation. The time discretization steps ∆t ww3 for WW3 and ∆t nemo for NEMO are generally different with ∆t ww3 > ∆t nemo and chosen such that ∆t ww3 = n t ∆t nemo (n t ∈ N, n t ≥ 1). In this case, coupling fields from NEMO to WW3 are averaged in time between two exchanges, while fields from WW3 to NEMO are sent every ∆t ww3 steps and therefore updated every n t time steps in NEMO. If ∆t ww3 > ∆t nemo , the coupler time-step is set to ∆t ww3 . Our current implementation 335 does not include an explicit coupling between waves and sea-ice while it is known that waves lead to ice break-up, pancake ice formation and associated enhancement of both freezing and melting and, in return, this wave dissipation in ice-covered water (e.g. Stopa et al., 2018) leads to ice drift. Such explicit coupling is currently under development within the NEMO framework (Boutin et al., 2019).  scheme (a.k.a. Ultimate Quickest scheme) is used with a specific procedure to alleviate the so-called garden sprinkler effect (Tolman et al., 2002). As suggested in Phillips (1984), the dissipation induced by wave breaking is proportional to the local saturation spectrum (see also Ardhuin et al., 2010a;Rascle and Ardhuin, 2013). The wind input growth rate at high frequency is based on the formulation of Janssen (1991) with an additional "sheltering" term to reduce the effective winds for the shorter waves (Chen and Belcher, 2000;Banner and Morison, 2010). For the computation of nonlinear wave-wave interactions, the 350 discrete interaction approximation of Hasselmann et al. (1985) is used. This last approximation is known to be inaccurate but it is thought that the associated error are usually compensated by a proper adjustment of the dissipation source term (Banner and Young, 1994;Ardhuin et al., 2007). As mentioned earlier in Sec. 3.2, the model was run with 10 meter winds, without any airsea stability correction. No wave measurements were assimilated in the model but the stand-alone wave model was developed based on spectral buoy and SAR data (Ardhuin et al., 2010b), and calibrated against altimeter data by adjusting the wind-wave 355 coupling parameter (Rascle and Ardhuin, 2013). The WW3 time-step for the global configurations is ∆t ww3 = 3600 s.
For the oceanic component, we use a global ORCA025 configuration at a 1/4 o horizontal resolution (Barnier et al., 2006).
The vertical grid is designed with 75 vertical z-levels with vertical spacing increasing with depth. Grid thickness is about 1 m near the surface and increases with depth to reach 200 m at the bottom. Partial steps are used to represent the bathymetry. coefficients are obtained from the 1-equation TKE scheme described in Sec. 2.3 and the convective processes are mimicked using an enhanced vertical diffusion parameterization which increases vertical diffusivity to 10 m 2 s −1 where static instability occurs. Water density is computed from temperature and salinity through the use of a polynomial formulation of the UNESCO (1983) non-linear equation of state (Roquet et al., 2015).The vector-invariant form momentum advection is using Arakawa and Lamb (1981) for the vorticity and a specific formulation to control the Hollingsworth instability (Ducousso et al., 2017).

365
Momentum lateral viscosity is biharmonic and acts along geopotential surfaces. It is set to a value of 1.5 × 10 11 m 4 .s −1 at the equator and vary proportionally to ∆x 3 away from the equator. Advection of tracers is performed with a flux-correctedtransport (FCT) scheme (Lévy et al., 2001), and lateral diffusion of tracers is harmonic and acts along iso-neutral surface. It is set to a value of 300 m 2 s −1 at the equator which varies proportionally to ∆x. The bottom friction is non-linear and the lateral boundary condition is free-slip. In this setup, the baroclinic time step is set to ∆t nemo = 900s, and a barotropic time step 30 370 times smaller. Compared to the standard uncoupled ORCA025 configuration, the additional computational cost associated to WW3 and to the exchanges through the coupler is about 20%.

Atmospheric forcings
The atmospheric fields used to force both ocean and wave models are based on the ECMWF (European Centre for Medium-Range Weather Forecasts) ERA-Interim reanalysis (Dee et al., 2011). Corrections have been applied to guarantee that the 375 ERA-Interim mean state for rainfalls, shortwave and longwave radiative fluxes are consistent with satellite observations from Remote Sensing Systems (RSS) Passive Microwave Water Cycle (PMWC) product (Hilburn, 2009) and GEWEX SRB 3.1 data 5 . Momentum and heat turbulent surface fluxes are computed using the IFS bulk formulation from AeroBulk package (Brodeau et al., 2016) using air temperature and humidity at 2 meters, mean sea level pressure and 10 meter winds. The No_CPL experiment corresponds to the classical NEMO setup where wave effect is parameterized through a wind-stress dependent TKE surface boundary condition as suggested by Craig and Banner (1994). In this approach, a Dirichlet surface boundary condition is used and expressed as follow: e(z = η) = 1 2 (15.8α CB ) 2/3 τ atm ρ 0 with α CB = 100. Based on the results of Mellor and Blumberg (2004)  Neumann) does not significantly impact numerical solutions 6 . The WS_CPL experiment is identical as No_CPL except that the wave coupling is introduced within the wind-stress computation, as described in Sec. 3.2. ST_CPL experiment is as WS_CPL except that all terms relative to the Stokes drift described in Sec. 2.1 are added in NEMO. TKE_CPL corresponds to ST_CPL but with the modified TKE scheme described in 2.3.1. All_CPL(1&2) experiments are like TKE_CPL but with a fully modified TKE scheme including the Langmuir cells parameterization described in 2.3.2. All those simulations have been 395 performed for 2 years (2013)(2014) where 2013 is let as spinup and only 2014 is analysed. We considered 2 years were enough to illustrate the fact that our developments were actually producing the expected results. Integrating longer in time could also lead to drifts in the stratification independently from the wave effects and could thus distort our interpretation. In any case, it must be clear that the objective here is not to go through a thorough physical analysis of coupled solutions but to check and validate our numerical developments.

Waves impact on oceanic Wind stress
The wave distribution being inhomogeneous on the globe, it is expected that with the wave-modified wind stress parameterization the stress should follow more closely the wave patterns. In Fig. 3, the seasonal average of the significant wave height and of the difference between the Charnock coefficient computed by the wave model and the default constant value used in 405 the uncoupled case (α 0 ch = 0.018) are shown. As expected, the Charnock parameter tends to be stronger in the area where the waves are the higher. Generally an increase of the Charnock parameter is observed in the northern and southern basin while there is a net decrease of α ch near the equator. There is also a strong seasonality in the northern hemisphere with a reduction in summer and a strong increase in winter. The differences between α ch and α 0 ch are very latitudinal with very few longitudinal variations.  To isolate the effect of the Charnock parameter we compare the results obtained in the No_CPL and WS_CPL experiments.
Those two experiments show relatively similar sea surface temperature patterns meaning that the modification of the windstress τ oce between those two cases is primarily due to the use of different Charnock parameters and the inclusion of the wave-supported stress. Fig. 4 (panel a) illustrates that the Charnock parameter mostly affects the drag coefficient C D , hence the surface wind-stress, for large winds. The ocean-wave coupling does not lead to appreciable differences in the drag coefficient 415 C D for wind speeds lower than 8 m s −1 . On the contrary, since large values of the Charnock parameter are observed for large wind speeds, the coupling significantly increases the drag (as well as its variance) at high winds. Fig. 4 (panel b) shows how the wind-stress is modified by this increase of the drag coefficient jointly with the wave-supported stress which tends to decrease the wind-stress magnitude (Fig. 5). At low wind speed the wind-stress magnitude is not affected by the coupling with waves while for strong winds the increase of wind-stress associated with the increased drag coefficient is always larger than 420 the decrease associated to the wave-supported stress. This latter effect reduces the wind stress by no more than 2%, for the characteristic scales of our study, this correction is thus almost negligible. The wind-stress changes due to the coupling with

Waves impact on surface TKE injection 425
As described in section 2.3, in the ocean-waves coupled case, the surface boundary condition for the TKE equation is a Neumann condition whose value is directly given by the wave model, unlike the uncoupled case where a Dirichlet condition is imposed. We aim here at assessing the impact on the order of magnitude of the near-surface TKE. Since the Neumann boundary condition is applied at the center of the top-most grid box (i.e. approximately at 50 cm depth), we compare in Fig.   6 the TKE value at 1 meter depth between the coupled (All_Cpl2) and the uncoupled (No_CPL) case. Positive values means 430 that near-surface TKE is larger in the coupled simulation. It shows an almost homogeneous increase of the TKE (up to more than 100%) in the extra-tropical areas. While low seasonal variability in the extra-tropical areas is visible in Fig. 6, a spatial averaging by hemisphere (Fig. 7) highlights seasonal variability with a strong increase in both near-surface TKE value and TKE difference between both experiments during winter. In Fig. 6, 7 (and also in the remainder of the paper), the spatial averaging is made between 25 S and 60 S in the southern hemisphere and between 25 N and 60 N in the northern hemisphere to avoid any conflicts with sea-ice and to remove the equatorial region from the comparison. The increase of the surface TKE injection associated with waves is expected to contribute to an overall increase of mixed layer depth provided that the mixing length diagnosed by the turbulent closure scheme allows to effectively propagate this additional TKE deeper in the mixed layer.

Waves impact on Mixed layer depth
In this section, we evaluate the wave effect on vertical mixing using the mixed layer depth (MLD) as a relevant metric. Fig. 8 440 represents the seasonally averaged difference in MLD between the coupled (All_CPL2) and the uncoupled (No_CPL) case relative to the No_CPL case (i.e. (h nocpl mld − h cpl mld )/h nocpl mld with h mld considered negative downward). It shows a significant deepening of the mixed layer at high latitudes in the coupled case with only very few localized mixed layer shallowing up to 60% mainly in the southern hemisphere. To assess whether the overall deepening of the mixed layer is realistic, we make a comparison with available observations. Available observations for 2014 were extracted following an updated data set from de 445 Boyer Montégut et al. (2004). The MLD depth as been computed as being the depth where the density is 3% smaller that the density at 10m as in de Boyer Montégut et al. (2004). Fig. 9 represents the spatially averaged MLD where the blue line is the spatially averaged MLD obtained from ARGO floats (available during the same period) in both hemispheres. In the northern hemisphere (Fig. 9, a), there is only a slight improvement compared to data during winter and late summer when implementing the coupling with waves. In the southern hemisphere ( Fig. 9, b) the situation is rather different. From January to July, the 450 deepening of MLD induced by the wave coupling significantly reduces the bias between the model and ARGO data. From July to December, results in the coupled case show an overestimation of MLDs which were already too deep in the uncoupled case, therefore increasing the bias between data and model. Since mesoscale activity make direct comparisons to data unreliable for such a short period of time, we compare the normalized distribution of MLD between the different simulations and available ARGO data. Results are presented in Fig. 10 for year 2014 (panel (a)) and during summer only (panel (b)). In both cases 455 the improvement in the northern hemisphere is very modest. As far as the southern hemisphere is concerned the coupling with waves leads to a significant improvement compared to MLD derived from ARGO floats despite the fact that there are still too many low MLD values in the range 50 − 100 m. In comparison with the uncoupled case there is a more realistic spreading toward deeper mixed layer depths. More particularly in summer (Fig. 10, b), the probability density function (PDF) in the coupled case matches almost perfectly the one computed from ARGO data. Despite the fact that we did not activate the 460 ad-hoc extra mixing induced by near-inertial waves (Rodgers et al., 2014) our implementation of the wave-ocean interaction leads to a significant deepening of the MLD in a realistic way. To better understand which components of the wave-ocean   coupling are responsible for this improvement, the summer PDF in the South hemisphere has been computed for each of the experiments described in Tab. 2. Results are shown in Fig. 11. First of all, it can be seen that all the wave-ocean interaction described in previous sections lead to an improvement in terms of mixed layer depth distribution compared to the uncoupled 465 case. Indeed, the modification of the wind stress by the wave field introduced in WS_CPL, increases both surface currents and near surface TKE values resulting in a slight deepening of the MLD. Adding the Stokes drift related terms in the primitive equations contributes only modestly to the deepening of the MLD while most of the improvement results from the modified TKE scheme with some slight improvement when the Langmuir parameterization is activated. It is somewhat reassuring to see that the better agreement with ARGO data is obtained when all components of the coupling are activated.

Waves impact on sea-surface temperature
Since the near-surface mixing is strengthened by the coupling we can expect an impact on sea-surface temperature (SST). Fig.   12 represents the time series of SST for each hemisphere. The Northern hemisphere is characterized by a warm bias during summer with a very slight improvement when coupling with waves. In the Southern hemisphere (Fig. 12, b) the summer warm bias is reduced by half in the coupled simulation and a slight warming occurs during the winter. While the summer surface 475 cooling might be linked to the mixed-layer deepening, the winter warming might be rather linked to advection as observed by Alari et al. (2016) for the Baltic sea. It could also result from an increased heat content during summer leading to higher SST during winter. To better characterize the wave impact on the SST, we show in Fig. 13 (panel a) the difference in term of annual mean between the No_CPL experiment and OSTIA analysis exhibiting a cold bias in the No_CPL simulation in equatorial and tropical regions and a warm bias in the northern part of the Pacific ocean. The coupling with waves tends to diminish the 480 cold bias (see Fig. 13, b) especially in the Pacific ocean and the warm bias in the north Pacific is significantly reduced. As wind stress caused by a value of the Charnock parameter lower than the value used in the uncoupled case (see Fig. 3,b,d). A consequence is a decrease of the drag coefficient leading to smaller turbulent exchange coefficients reducing the heat flux. As mentionned above, in extra-tropical regions, some warm bias tend to be partially reduced by the extra mixing induced by the 485 waves at high latitude or/and by the increased turbulent transfer coefficient. The tendency of the wave coupling to improve the near-surface temperature distribution can also be verified on a time-latitude Hovmuller diagram like the ones shown in Fig. 14. For instance, it can be seen that the summer warm bias in the northern hemisphere (Fig. 14,a) coincides well with the cooling induced by the coupling with waves (Fig. 14,b). Similarly we can also observe a warming in the tropical and equatorial regions ( Fig. 14,b) corresponding to the cold bias seen in Fig. 14 (panel a). In the southern extra-tropical region, a summer cooling is 490 observed. It is induced by the wave coupling whereas Fig. 14 (panel a) shows a slight warm bias. During winter we can observe north of 60 S a warming in Fig. 14 (panel b) which again partially corresponds to a cold bias in Fig. 14 (panel a).

Surface current and Kinetic Energy
The last aspect of our solutions we would like to evaluate is the impact of the surface waves on surface currents and kinetic energy (KE). To do so, we show in Fig. 15 time series of the spatially averaged surface kinetic energy for both hemispheres.

495
Whatever the hemisphere there is a net decrease of surface KE (up to 20% in the south) when a coupling with the waves is included. This decrease of surface kinetic energy reflects a decrease of surface currents magnitude. Indeed, as detailed in Fig.   16 which represents the vertical profile of the horizontal components of the current in the oceanic surface boundary layer, the coupling with waves decreases both the surface currents magnitude and the shear. While currents from the WS_CPL are increased due to increased wind stress, the Stokes Coriolis force when included in momentum equations leads to a decrease of 500 velocities in the whole boundary layer as previously shown by Rascle et al. (2008) (orange lines in Fig. 16). Inclusion of the vertical mixing due to waves and Langmuir circulation attenuates the currents in the surface layer, resulting in further reduced surface currents and stronger currents at the bottom of the boundary layer (purple lines in Fig. 16). This concludes our checking of the proper functioning of the coupling with waves as described in the present paper.

505
In this paper we have described the implementation of an online coupling between the oceanic model NEMO and the wave model WW3. The impact of such coupling on the model solutions has been assessed from the oceanic point of view for a global configuration. In particular, the following steps to set up the coupled model have been discussed in details (i) inclusion of all wave-induced terms in NEMO primitive equations, only neglecting the terms relevant for the surf zone which is outside the scope of the NEMO community, (ii) modification of the subgrid scale vertical physics (including the bulk formulation) 510 to include wave effects and a parameterization of Langmuir turbulence, (iii) development of a coupling interface based on the OASIS3-MCT software for the exchange of data between both models, and (iv) tests of our developments on a realistic global configuration at 1/4 • for the ocean coupled to a 1/2 • resolution wave model. Compared to an ocean-only simulation, the coupling with a wave model (with a resolution twice coarser than the oceanic model) leads to an additional computational cost of about 20%.

515
Following McWilliams et al. (2004) and Ardhuin et al. (2008), in the weak vertical current shears limit, the wave-induced terms implemented in NEMO include the Stokes-Coriolis force, the vortex force, Stokes advection in tracer and continuity equations as well as a wave-induced surface pressure term. The prognostic equation for TKE also includes an additional forcing term associated with the Stokes drift vertical shear as well as various modifications of its boundary condition described in Sec. 2.3. 520 The development of a coupling infrastructure based on OASIS3-MCT has several advantages as it allows for an efficient data exchange (including the treatment of non-conformities between the computational grids) but also for versatility in the inclusion of a wave model in existing ocean-atmosphere or ocean-only models. At a practical level, the OASIS interface we have implemented in NEMO is similar to other interfaces (e.g. toward atmospheric models) existing in the code which is important for maintenance and for further developments. It paves the way for a seamless and more systematic inclusion of the 525 coupling with waves for NEMO users.
Unlike most previous studies of wave-ocean coupling using NEMO, we have shown that satisfactory results can be obtained from the TKE vertical turbulent closure scheme without activating the ad hoc parameterization for the mixing induced by near-inertial waves, surface waves and swell (known as ETAU parameterization). This parameterization which amounts to empirically propagate the surface TKE at depth using a prescribed shape function is a pragmatic way to cure the shallow mixed 530 layer depths in the southern ocean found in simulations ignoring wave effects. Previous studies of wave-ocean coupling by Breivik et al. (2015), Alari et al. (2016) or Staneva et al. (2017) have been using the ETAU parameterization in their setup. However, as suggested by Breivik et al. (2015), we can speculate that such parameterization could mask the impact of the wave coupling even though it turned out to be necessary to obtain realistic mixed layer depths. We believe that our modification of the standard NEMO 1-equation TKE scheme described in Sec. 2.3 is more physically justifiable than the ETAU parameterization 535 and requires much less parameter tuning.
The numerical experiments based on the ORCA25 configuration discussed in Sec. 4.2 were meant to check that our developments were having the expected impact on numerical solutions. First, we confirmed that using the Charnock parameter computed in the wave model instead of a constant value globally increases the wind-stress magnitude, particularly at mid and high latitudes whereas accounting for the portion of the wind-stress consumed by the waves has a small impact (in our 540 experiments it leads to a maximum of 2% decrease of the wind-stress). Second, using the mixed layer depth as an indicator to assess the amount of vertical mixing, the modifications brought to the NEMO turbulence scheme (i.e the new boundary condition for TKE and for the mixing length, the addition of the Stokes shear in the TKE equation, and the modified Axell (2002) parameterization for Langmuir cells) lead to an important extra mixing contributing to a deepening of the surface mixed layer particularly in the southern hemisphere. When compared to ARGO data it shows a significant improvement during the summer, 545 while during the winter the extra wave-induced mixing deepens the already too deep mixed layer. Note that the Fox-Kemper et al. (2008) parameterization to account for the restratification induced by mixed layer instabilities (Boccaletti et al., 2007;Couvelard et al., 2015) during the winter was not used in our experiments. This parameterization induces even more shallow summer mixed layer depths. As far as the northern hemisphere is concerned, coupled results show an improvement when compared to ARGO for winter with a deepening of the mixed layer while in summer results are similar to the uncoupled case.

550
Since the comparison with ARGO data can be tricky due to the scarcity of the data, we looked at the results in terms of mixed layer depths (MLD) probability density functions. This allowed to highlight the significant improvement in MLD distribution when coupling with the waves. Furthermore, we noticed that all components of the ocean-wave coupling act to deepen the mixed layer and therefore have a cumulative effect. However the main contributor is the fully modified TKE scheme including Langmuir cell parameterization of Axell (2002) which is consistent with recent results obtained by Reichl et al. (2016) and Ali Since the magnitude of the vertical mixing is increased by the coupling with waves we expect an impact on sea surface temperature and currents. Indeed, the summer deepening of the mixed layer in the southern hemisphere leads to colder sea surface temperatures resulting in a better agreement with OSTIA SST analysis. More generally, although the global SST biases are not totally compensated, they tend to be reduced when considering the effect of waves (see Sec. 4.2.4). The currents in the oceanic surface boundary layer are reduced by the Stokes Coriolis force (which counteracts the Ekman current, Rascle et al., 2008). They are also affected by the increased vertical mixing which tends to reduce the surface currents (and thus the surface kinetic energy) and strengthen the currents at the base of the surface boundary layer. The reduction of surface kinetic energy due to the wave-ocean coupling in the global 1/4 • resolution configuration is of the same order of magnitude as the reduction observed when accounting for surface currents in the computation of the wind stress in a coupled ocean-atmosphere model 565 (e.g. Renault et al., 2016). A fully coupled ocean-wave-atmosphere model would thus be necessary to properly disentangle the different contributions at play impacting the oceanic surface kinetic energy. Even if additional diagnostics on various configurations at different resolutions are still needed to exhaustively evaluate the impact of each component of the ocean wave coupling, the results presented in the paper confirm the robustness of our developments and our implementation will serve as a starting point for the inclusion of wave-currents interactions in the forthcoming NEMO official release. We can 570 speculate that the ocean-waves coupled ORCA025 configuration might become a standard component of future Coupled Model Intercomparison Project (CMIP) exercises. We already mentioned as a perspective the addition of a coupling with an interactive atmospheric boundary layer either via a full atmospheric model or a simplified boundary layer model (e.g. Lemarié et al., 2020). Furthermore, the gain of an online 2-way coupling compared to a 1-way coupling on the oceanic as well as on the wave solution must be investigated in the future. Indeed, the improvements of the quality of surface waves simulations associated to 575 a coupling with large-scale oceanic currents are well documented particularly in the Agulhas current (Irvine and Tilley, 1988) and in the Gulf Stream (Mapp et al., 1985). Ardhuin et al. (2017a) have also shown a strong impact of small-scale currents (10-100km) on wave height variability at the same scales. We can therefore expect improvements for both wave and ocean forecasts when the coupling is implemented in an operational context.
Code and data availability. The changes to the NEMO code have been made on the standard NEMO code (nemo_v3_6_STABLE). The 580 code can be downloaded from the NEMO website (http://www.nemo-ocean.eu/, last access: 11 July 2019). The NEMO code modified to include wave-ocean coupling terms and the OASIS interface is available in the zenodo archive (https://doi.org/10.5281/zenodo.3331463, Couvelard (2019)). The WW3 code version 6.02 has been used without further modifications and can be downloaded from the NOAA github repository (https://github.com/NOAA-EMC/WW3, last access: 11 July 2019). Our modifications of the OASIS interface in the WW3 code have already been integrated in the official release. The OASIS3_MCT code is also freely available (https://portal.enes.org/oasis/, last access: 585 11 July 2019). The exact versions of the WW3 and OASIS3_MCT codes that were used have also been made available in the zenodo archive (https://doi.org/10.5281/zenodo.3331463, Couvelard (2019)) The initial and forcing data for both the oceanic and wave model, analysis scripts, namelists and data used to produce the figures are also available in the zenodo archive.

Appendix A: Flux-form wave-averaged momentum equations
In this appendix we describe the necessary changes when a flux formulation for advective terms in the momentum equations is 590 preferred to the vector invariant form presented in (7) and (8). For simplicity, we consider just the i-component in horizontal curvilinear coordinates and a z-coordinate in the vertical (results will be extended to the j-component and to generalized vertical coordinate). Consistently with the notations of Madec (2012), e 1 and e 2 are the horizontal scale factors. We note A u v the extra term needed to guarantee the equivalence between the flux formulation and the vector-invariant form. A u v is defined such that 595 ∇ · (u s u) + A u v = −ζv s + w s e 3 ∂ k u.
Since ∇ · u s = 0 we have ∇ · (u s u) = u s · ∇u, and thus e 1 e 2 A u v = −v s [∂ i (e 2 v) − ∂ j (e 1 u)] + e 1 e 2 e 3 w s ∂ k u − e 2 u s ∂ i u + e 1 v s ∂ j u + e 1 e 2 e 3 w s ∂ k u = −v s [v∂ i e 2 − u∂ j e 1 + e 2 ∂ i v − e 1 ∂ j u] − e 2 u s ∂ i u − e 1 v s ∂ j u hence 600 A u v = − v s e 1 e 2 (v∂ i e 2 − u∂ j e 1 )

Additional term
Same computation for the j-component leads to the following equations in generalized vertical coordinates