Development of a two-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 upper-oceanic flow dynamics with surface waves. The focus is on the impact of such coupling on upper-ocean properties (temperature and currents) and mixed layer depth (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 modifications of the tracer–continuity equation and turbulent closure scheme (a one-equation turbulent kinetic energy – TKE – closure here). To assess the new developments, we perform a set of sensitivity experiments with a global oceanic configuration at $\mathrm{1}/\mathrm{4}{}^{\circ }$ resolution coupled with a wave model configured at $\mathrm{1}/\mathrm{2}{}^{\circ }$ resolution. Numerical simulations show a global increase in wind stress due to the interaction with waves (via the Charnock coefficient), particularly at high latitudes, resulting in increased surface currents. The modifications brought to the TKE closure scheme and the inclusion of a parameterization for Langmuir turbulence lead to a significant increase in the mixing, thus helping to deepen the MLD. This deepening is mainly located in the Southern Hemisphere and results in reduced sea surface currents and temperatures.

1 Introduction

An accurate representation of ocean surface waves has long been recognized as essential for a wide range of applications from marine meteorology to ocean and coastal engineering. Waves also play an important role in the short-term forecasting of extratropical and tropical cyclones by regulating sea surface roughness (Janssen, 2008; Chen and Curnic, 2015; Hwang, 2015). More recently, the impact of waves on oceanic circulation at the global scale has triggered interest from the research and operational community (e.g., Hasselmann, 1991; Rascle and Ardhuin, 2009; D'Asaro et al., 2014; Fan and Griffies, 2014; Li et al., 2016; Law Chune and Aouf, 2018). In particular, surface waves are important for an accurate representation of air–sea interactions, and their effect on fluxes of mass, momentum, and energy through the wavy boundary layer must be taken into account in ocean–atmosphere coupled models. For example, the momentum flux through the air–sea interface has traditionally been parameterized using near-surface winds (typically at 10 m) and the atmospheric surface layer stability (Fairall et al., 2003; Large and Yeager, 2009; Brodeau et al., 2016). The physics of the coupling depend on the kinematics and dynamics of the wave field. This includes a wide range of processes from wind–wave growth, nonlinear wave–wave interaction, and 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, 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. Indeed, their large eddy simulations (LESs) 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 cell 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 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 cell parameterizations leads to excessive mixing, particularly in winter.

Most previous studies of the impact of ocean–wave interactions at the global scale have used 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 system with 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; WAVEWATCH III® Development Group, 2016) as the surface wave component. NEMO and WW3 are coupled using the OASIS Model Coupling Toolkit (OASIS-MCT; Valcke, 2013; Craig et al., 2017), which is widely used in the climate and operational community. The various steps for our implementation are the following: (i) the 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; and (iv) a test of the implementation based on a realistic global configuration at $\mathrm{1}/\mathrm{4}{}^{\circ }$ for the ocean and $\mathrm{1}/\mathrm{2}{}^{\circ }$ for the waves.

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 Sect. 2. This includes the addition of the Stokes–Coriolis force, the vortex force, and the wave-induced pressure gradient. In Sect. 3 our modeling system coupling the NEMO oceanic model and the WW3 wave model via the OASIS-MCT coupler is described in detail. Numerical simulations are presented in Sect. 4 using a global configuration at $\mathrm{1}/\mathrm{4}{}^{\circ }$ for the oceanic model and $\mathrm{1}/\mathrm{2}{}^{\circ }$ for the wave model. Using sensitivity runs, we assess those global configurations with particular emphasis on the impact of wave–ocean interactions on mixed layer depth, sea surface temperature and currents, turbulent kinetic energy (TKE) injection, and kinetic energy. Finally, in Sect. 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

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 the first one here, which will be used for the numerical simulations in Sect. 4. With ${\mathit{u}}_{\mathrm{h}}=(u,v)$ the horizontal velocity vector, ω the dia-surface velocity component, θ the potential temperature, ρ the density, ζ the relative vorticity, p_h the hydrostatic pressure, and p_s the surface pressure, the Reynolds-averaged equations (with 〈⋅〉 the averaging operator omitted here for simplicity) are as follows.

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\partial }_{t}u=& +(f+\mathit{\zeta })v-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\partial }_x\Vert {\mathit{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}-{\displaystyle \frac{\mathit{\omega }}{e_{\mathrm{3}}}}{\partial }_{k}u\\ & -{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left({\partial }_x(p_{\mathrm{s}}+p_{\mathrm{h}})-\left({\partial }_{x}z\right){\displaystyle \frac{\left({\partial }_k{p}_{\mathrm{h}}\right)}{e_{\mathrm{3}}}}\right)\\ & -{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k〈u^{\prime }{\mathit{\omega }}^{\prime }〉+F^u\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\partial }_{t}v=& -(f+\mathit{\zeta })u-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\partial }_y\Vert {\mathit{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}-{\displaystyle \frac{\mathit{\omega }}{e_{\mathrm{3}}}}{\partial }_{k}v\\ & -{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left({\partial }_y(p_{\mathrm{s}}+p_{\mathrm{h}})-\left({\partial }_{y}z\right){\displaystyle \frac{\left({\partial }_k{p}_{\mathrm{h}}\right)}{e_{\mathrm{3}}}}\right)\\ & -{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k〈v^{\prime }{\mathit{\omega }}^{\prime }〉+F^v\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}{\partial }_t\left(e_{\mathrm{3}}\mathit{\theta }\right)=& -{\partial }_x\left(e_{\mathrm{3}}\mathit{\theta }u\right)-{\partial }_y\left(e_{\mathrm{3}}\mathit{\theta }v\right)-{\partial }_k\left(\mathit{\theta }\mathit{\omega }\right)\\ & -{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k〈{\mathit{\theta }}^{\prime }{\mathit{\omega }}^{\prime }〉+F^{\mathit{\theta }}\end{array}\text{(4)}& {\displaystyle }{\partial }_t{e}_{\mathrm{3}}=-{\partial }_x\left(e_{\mathrm{3}}u\right)-{\partial }_y\left(e_{\mathrm{3}}v\right)-{\partial }_k\mathit{\omega }\text{(5)}& {\displaystyle }{\partial }_{k}p=-\mathit{\rho }g{e}_{\mathrm{3}}\end{array}$$

Here, k is a nondimensional vertical coordinate, the lateral derivatives ∂_x and ∂_y have to be considered along the model coordinate, and e₃ is the vertical scale factor given by $e_{\mathrm{3}}={\partial }_{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,

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& w(z=\mathit{\eta })={\partial }_t\mathit{\eta }+{\left.u\right|}_{z=\mathit{\eta }}{\partial }_x\mathit{\eta }+{\left.v\right|}_{z=\mathit{\eta }}{\partial }_y\mathit{\eta }+(E-P),\\ & w(z=-H)=-{\left.u\right|}_{z=-H}{\partial }_{x}H-{\left.v\right|}_{z=-H}{\partial }_{y}H,\end{array}\end{array}$$

with η the height of the sea surface and (E−P) the mass flux across the sea surface due to precipitation and evaporation, a momentum surface boundary condition for the Reynolds stress vertical terms,

$$-{\left.〈u^{\prime }{\mathit{\omega }}^{\prime }〉\right|}_{z=\mathit{\eta }}={\displaystyle \frac{{\mathit{\tau }}_u^{\mathrm{oce}}}{{\mathit{\rho }}_{\mathrm{0}}}},{\left.-〈v^{\prime }{\mathit{\omega }}^{\prime }〉\right|}_{z=\mathit{\eta }}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{v}}^{\mathrm{oce}}}{{\mathit{\rho }}_{\mathrm{0}}}},$$

with ${\mathit{\tau }}^{\mathrm{oce}}=({\mathit{\tau }}_u^{\mathrm{oce}},{\mathit{\tau }}_{\mathrm{v}}^{\mathrm{oce}})$ a wind stress vector that represents the part of the stress that drives the ocean, and a dynamic 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 $\mathit{\omega }(z=\mathit{\eta })=(E-P)$ and $\mathit{\omega }(z=-H)=\mathrm{0}$. We do not explicitly include the boundary conditions for the tracer equations here since they are unchanged from classical primitive equation models in the presence of wave motions. As mentioned earlier, in Eqs. (1) to (5) prognostic variables have to be interpreted in an Eulerian mean sense even if the averaging operator is not explicitly included.

2.1 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 (Ardhuin et al., 2008) 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), and Moghimi et al. (2013). The three-component Stokes drift vector is ${\mathit{u}}^{\mathrm{s}}=({\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}},{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}},{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}})$ and is non-divergent at lowest order (Ardhuin et al., 2008, 2017b). The coupled wave–current equations for the Eulerian mean velocity and tracers in a vector-invariant form (the equivalent flux form is given in Appendix A) are as follows.

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}{\partial }_{t}u=& +(f+\mathit{\zeta })(v+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}})-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\partial }_x\Vert {\mathit{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}-{\displaystyle \frac{(\mathit{\omega }+{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}})}{e_{\mathrm{3}}}}{\partial }_{k}u\\ & -{\displaystyle \frac{\left({\partial }_{x}z\right)}{e_{\mathrm{3}}}}\left({\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}{\partial }_{k}u+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}{\partial }_{k}v\right)-{\displaystyle \frac{{\partial }_x(p_{\mathrm{s}}+{\stackrel{\mathrm{̃}}{p}}^J)}{{\mathit{\rho }}_{\mathrm{0}}}}\\ & -{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left({\partial }_x(p_{\mathrm{h}}+{\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}})-\left({\partial }_{x}z\right){\displaystyle \frac{{\partial }_k(p_{\mathrm{h}}+{\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}})}{e_{\mathrm{3}}}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k〈u^{\prime }{\mathit{\omega }}^{\prime }〉+F^u+{\stackrel{\mathrm{̃}}{F}}^u\end{array}\text{(8)}& {\displaystyle }\begin{array}{rl}{\partial }_{t}v=& -(f+\mathit{\zeta })(u+{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}})-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\partial }_y\Vert {\mathit{u}}_{\mathrm{h}}{\Vert }^{\mathrm{2}}-{\displaystyle \frac{(\mathit{\omega }+{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}})}{e_{\mathrm{3}}}}{\partial }_{k}v\\ & -{\displaystyle \frac{\left({\partial }_{y}z\right)}{e_{\mathrm{3}}}}\left({\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}{\partial }_{k}u+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}{\partial }_{k}v\right)-{\displaystyle \frac{{\partial }_y(p_{\mathrm{s}}+{\stackrel{\mathrm{̃}}{p}}^J)}{{\mathit{\rho }}_{\mathrm{0}}}}\\ & -{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left({\partial }_y(p_{\mathrm{h}}+{\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}})-\left({\partial }_{y}z\right){\displaystyle \frac{{\partial }_k(p_{\mathrm{h}}+{\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}})}{e_{\mathrm{3}}}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k〈v^{\prime }{\mathit{\omega }}^{\prime }〉+F^v+{\stackrel{\mathrm{̃}}{F}}^v\end{array}\text{(9)}& {\displaystyle }\begin{array}{rl}{\partial }_t\left(e_{\mathrm{3}}\mathit{\theta }\right)=& -{\partial }_x\left(e_{\mathrm{3}}\mathit{\theta }\right(u+{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}})-{\partial }_y(e_{\mathrm{3}}\mathit{\theta }(v+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}))\\ & -{\partial }_k\left(\mathit{\theta }\right(\mathit{\omega }+{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}}\left)\right)-{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k〈{\mathit{\theta }}^{\prime }{\mathit{\omega }}^{\prime }〉+F^{\mathit{\theta }}\end{array}\text{(10)}& {\displaystyle }{\partial }_t{e}_{\mathrm{3}}=-{\partial }_x\left(e_{\mathrm{3}}\right(u+{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}\left)\right)-{\partial }_y\left(e_{\mathrm{3}}\right(v+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}\left)\right)-{\partial }_k(\mathit{\omega }+{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}})\text{(11)}& {\displaystyle }{\partial }_k\left(p_{\mathrm{h}}+{\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}}\right)=-\mathit{\rho }g{e}_{\mathrm{3}}+{\mathit{\rho }}_{\mathrm{0}}\left({\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}{\partial }_{k}u+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}{\partial }_{k}v\right)\end{array}$$

Here, wave-induced terms are represented with tildes. The ${\stackrel{\mathrm{̃}}{F}}^u$ and ${\stackrel{\mathrm{̃}}{F}}^v$ terms represent the sink–source of wave momentum 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 𝓦_St−Cor, the vortex force 𝓦_VF, and a wave-induced pressure 𝓦_Prs:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\mathcal{W}}_{\mathrm{St}-\mathrm{Cor}}=\left(\begin{array}{c}f{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}\\ -f{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}\\ \mathrm{0}\end{array}\right),\\ & {\mathcal{W}}_{\mathrm{VF}}=\left(\begin{array}{c}\mathit{\zeta }{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}-\frac{{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}}}{e_{\mathrm{3}}}{\partial }_{k}u-\frac{\left({\partial }_{x}z\right)}{e_{\mathrm{3}}}\left({\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}{\partial }_{k}u+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}{\partial }_{k}v\right)\\ -\mathit{\zeta }{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}-\frac{{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}}}{e_{\mathrm{3}}}{\partial }_{k}v-\frac{\left({\partial }_{y}z\right)}{e_{\mathrm{3}}}\left({\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}{\partial }_{k}u+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}{\partial }_{k}v\right)\\ \frac{{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}}{e_{\mathrm{3}}}{\partial }_{k}u+\frac{{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}}{e_{\mathrm{3}}}{\partial }_{k}v\end{array}\right),\\ & {\mathcal{W}}_{\mathrm{Prs}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left(\begin{array}{c}{\partial }_x\left({\stackrel{\mathrm{̃}}{p}}^J+{\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}}\right)-\left({\partial }_{x}z\right)\frac{{\partial }_k\left({\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}}\right)}{e_{\mathrm{3}}}\\ {\partial }_y\left({\stackrel{\mathrm{̃}}{p}}^J+{\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}}\right)-\left({\partial }_{y}z\right)\frac{{\partial }_k\left({\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}}\right)}{e_{\mathrm{3}}}\\ \frac{\mathrm{1}}{e_{\mathrm{3}}}{\partial }_k{\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}}\end{array}\right),\end{array}\end{array}$$

where the terms involving horizontal derivatives of ω have been neglected in 𝓦_VF. In 𝓦_Prs, the ${\stackrel{\mathrm{̃}}{p}}^J$ term corresponds to a depth-uniform wave-induced kinematic pressure term¹, while ${\stackrel{\mathrm{̃}}{p}}^{\mathrm{Shear}}$ is a shear-induced three-dimensional pressure term² associated with the vertical component of the vortex force. The vortex force contribution 𝓦_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 pressure term ${\stackrel{\mathrm{̃}}{p}}^{\mathrm{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 generally weak vertical shears in the wave mixed layer. The effect of that term was also found to be much weaker than ${\stackrel{\mathrm{̃}}{p}}^J$ in shallow coastal environments, except in the surf zone. This assumption has the advantage of leaving the hydrostatic relation (Eq. 11) unchanged. Our implementation of wave-induced terms in NEMO is in line with Bennis et al. (2011) and corresponds to the simplified form of Eq. (12):

$$\begin{array}{rl}& {\mathcal{W}}_{\mathrm{St}-\mathrm{Cor}}=\left(\begin{array}{c}f{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}\\ -f{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}\\ \mathrm{0}\end{array}\right),{\mathcal{W}}_{\mathrm{VF}}=\left(\begin{array}{c}\mathit{\zeta }{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}-\frac{{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}}}{e_{\mathrm{3}}}{\partial }_{k}u\\ -\mathit{\zeta }{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}-\frac{{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}}}{e_{\mathrm{3}}}{\partial }_{k}v\\ \mathrm{0}\end{array}\right),\\ & {\mathcal{W}}_{\mathrm{Prs}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left(\begin{array}{c}{\partial }_x{\stackrel{\mathrm{̃}}{p}}^J\\ {\partial }_y{\stackrel{\mathrm{̃}}{p}}^J\\ \mathrm{0}\end{array}\right).\end{array}$$

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 tracer 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 been modified accordingly since the surface kinematic boundary condition (Eq. 6) in terms of vertical velocities w and associated ${\stackrel{\mathrm{̃}}{w}}^{\mathrm{s}}$ now reads

$$\begin{array}{rl}w+{\stackrel{\mathrm{̃}}{w}}^{\mathrm{s}}=& {\partial }_t\mathit{\eta }+({\left.u\right|}_{z=\mathit{\eta }}+{\left.{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}\right|}_{z=\mathit{\eta }}){\partial }_x\mathit{\eta }\\ & +({\left.v\right|}_{z=\mathit{\eta }}+{\left.{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}\right|}_{z=\mathit{\eta }}){\partial }_y\mathit{\eta }+(E-P)\end{array}$$

to express the fact that there is a source of mass at the surface that compensates for the convergence of the Stokes drift; hence, the barotropic mode is

$$\begin{array}{}\text{(13)}& \left\{\begin{array}{l}{\partial }_t\mathit{\eta }=-{\partial }_x\left((H+\mathit{\eta })(\stackrel{\mathrm{‾}}{u}+{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{u}}}^{\mathrm{s}})\right)\\ -{\partial }_y\left((H+\mathit{\eta })(\stackrel{\mathrm{‾}}{v}+{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{v}}}^{\mathrm{s}})\right)+P-E,\\ {\partial }_t\stackrel{\mathrm{‾}}{u}=+f\stackrel{\mathrm{‾}}{v}-g{\partial }_x\mathit{\eta }-\frac{C_{\mathrm{b},x}}{(H+\mathit{\eta })}\stackrel{\mathrm{‾}}{u}+\stackrel{\mathrm{‾}}{G^x}+\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{̃}}{G}}^x}\\ {\partial }_t\stackrel{\mathrm{‾}}{v}=-f\stackrel{\mathrm{‾}}{u}-g{\partial }_y\mathit{\eta }-\frac{C_{\mathrm{b},y}}{(H+\mathit{\eta })}\stackrel{\mathrm{‾}}{v}+\stackrel{\mathrm{‾}}{G^y}+\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{̃}}{G}}^y},\end{array}\right.\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathit{\varphi }}=\frac{\mathrm{1}}{H+\mathit{\eta }}{\int }_{-H}^{\mathit{\eta }}\mathit{\varphi }\mathrm{d}z$, ${\mathbf{C}}_{\mathrm{b}}=(C_{\mathrm{b},x},C_{\mathrm{b},y})$ represents the bottom drag coefficients, and $\mathbf{G}=(\stackrel{\mathrm{‾}}{G^x},\stackrel{\mathrm{‾}}{G^x})$ is the usual NEMO forcing term containing coupling terms from the baroclinic mode and slowly varying barotropic terms (including nonlinear advective terms) held constant during the barotropic integration to gain efficiency. In Eq. (13), $\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{̃}}{G}}^x}$ and $\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{̃}}{G}}^y}$ contain the additional wave-induced barotropic forcing terms corresponding to the vertical integral of the 𝓦_St−Cor and 𝓦_VF terms, which are also held constant during the barotropic integration. A thorough analysis of 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.

2.2 Computation and discretization of Stokes drift velocity profile

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 ${\mathit{u}}_{\mathrm{h}}^{\mathrm{s}}\left(\mathit{\eta }\right)$ and the norm of the Stokes volume transport $\Vert {\mathbf{T}}^{\mathrm{s}}\Vert$. In Breivik et al. (2014, 2016), Stokes drift velocity profiles are derived under the deepwater approximation in the general form ${\mathit{u}}_{\mathrm{h}}^{\mathrm{s}}\left(z\right)={\mathit{u}}_{\mathrm{h}}^{\mathrm{s}}\left(\mathit{\eta }\right)\mathcal{S}(z,k_e)$, with k_e a depth-independent spatial wavenumber chosen such that the norm of the depth-integrated Stokes transport (assuming an ocean of infinite depth) is equal to $\Vert {\mathbf{T}}^{\mathrm{s}}\Vert$. The functions 𝒮_B14(z,k_e) from Breivik et al. (2014) and 𝒮_B16(z,k_e) from Breivik et al. (2016) for $z\in [-H,\mathit{\eta }]$ are given by

$$\begin{array}{rl}& {\mathcal{S}}_{\mathrm{B}\mathrm{14}}(z,k_e)=\left({\displaystyle \frac{e^{\mathrm{2}{k}_e(z-\mathit{\eta })}}{\mathrm{1}-\mathrm{8}{k}_e(z-\mathit{\eta })}}\right),\\ & {\mathcal{S}}_{\mathrm{B}\mathrm{16}}(z,k_e)=e^{\mathrm{2}{k}_e(z-\mathit{\eta })}-\sqrt{\mathrm{2}{k}_e\mathit{\pi }(\mathit{\eta }-z)}\mathrm{erfc}\left(\sqrt{\mathrm{2}{k}_e(\mathit{\eta }-z)}\right).\end{array}$$

with erfc the complementary error function. It can be easily shown that for an ocean of infinite depth, the vertical integrals of those functions are respectively equal to $\frac{\mathrm{1}}{\mathrm{6}{k}_e}$ for 𝒮_B16 and $\frac{\mathrm{1.34089}}{\mathrm{8}{k}_e}\approx \frac{\mathrm{1}}{\mathrm{5.97}{k}_e}$ for 𝒮_B14. Standard computations of Stokes drift in numerical models are done in a finite-difference sense; however, due to the fast decay of ${\mathit{u}}_{\mathrm{h}}^{\mathrm{s}}\left(z\right)$ with depth, a finite-volume approach seems more adequate in this case.

$$\begin{array}{}\text{(14)}& \begin{array}{rl}({\mathit{u}}_{\mathrm{h}}^{\mathrm{s}}{)}_k& ={\displaystyle \frac{{\mathit{u}}_{\mathrm{h}}^{\mathrm{s}}\left(\mathit{\eta }\right)}{(e_{\mathrm{3}}{)}_k}}\underset{z_{k-\mathrm{1}/\mathrm{2}}}{\overset{z_{k+\mathrm{1}/\mathrm{2}}}{\int }}\mathcal{S}(z,k_e)\mathrm{dz}\\ & ={\displaystyle \frac{{\mathit{u}}_{\mathrm{h}}^{\mathrm{s}}\left(\mathit{\eta }\right)}{(e_{\mathrm{3}}{)}_k}}\left[\mathcal{I}(z_{k+\mathrm{1}/\mathrm{2}},k_e)-\mathcal{I}(z_{k-\mathrm{1}/\mathrm{2}},k_e)\right]\end{array}\end{array}$$

Such a finite-volume interpretation of the Stokes drift velocity can also be found in Li et al. (2017) and Wu et al. (2019). The 𝒮_B16 function is more adapted for this kind of approach since the primitive function only requires special functions available in the Fortran standard.

$${\mathcal{I}}_{\mathrm{B}\mathrm{16}}(z,k_e)={\displaystyle \frac{\mathrm{1}}{\mathrm{6}{k}_e}}\left[e^{\mathrm{2}{k}_e(z-\mathit{\eta })}+\mathrm{4}{k}_e(z-\mathit{\eta }){\mathcal{S}}_{\mathrm{B}\mathrm{16}}(z,k_e)\right]$$

Since NEMO is discretized on an Arakawa C grid, the components of the Stokes drift velocity must be evaluated at cell interfaces, and a simple average weighted by layer thicknesses is used:

$$\begin{array}{rl}& {\stackrel{\mathrm{̃}}{u}}_{i+\mathrm{1}/\mathrm{2},j,k}^{\mathrm{s}}={\displaystyle \frac{(e_{\mathrm{3}}{)}_{i,j,k}{u}_{i,j,k}^{\mathrm{s}}+(e_{\mathrm{3}}{)}_{i+\mathrm{1},j,k}{u}_{i+\mathrm{1},j,k}^{\mathrm{s}}}{\mathrm{2}(e_{\mathrm{3}}{)}_{i+\mathrm{1}/\mathrm{2},j,k}}},\\ & {\stackrel{\mathrm{̃}}{v}}_{i,j+\mathrm{1}/\mathrm{2},k}^{\mathrm{s}}={\displaystyle \frac{(e_{\mathrm{3}}{)}_{i,j,k}{v}_{i,j,k}^{\mathrm{s}}+(e_{\mathrm{3}}{)}_{i,j+\mathrm{1},k}{v}_{i+\mathrm{1},j,k}^{\mathrm{s}}}{\mathrm{2}(e_{\mathrm{3}}{)}_{i,j+\mathrm{1}/\mathrm{2},k}}}.\end{array}$$

Note that no explicit computation of the vertical component of the Stokes drift is necessary since in Eqs. (7)–(11) ${\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}}$ only appears summed with ω such that the relevant variable is $\mathit{\omega }+{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}}$ as a whole. This quantity is diagnosed from the continuity equation (Eq. 10, where the temporal evolution of vertical scale factors ∂_te₃ is given by the free-surface evolution when a quasi-Eulerian vertical coordinate is used; e.g., z^⋆ or terrain-following coordinates).

Figure 1 (a) Reconstructed zonal component of a Stokes drift profile for $\Vert {\mathbf{T}}^{\mathrm{s}}\Vert =\mathrm{0.4}$ m² s⁻¹, $u^{\mathrm{s}}(z=\mathit{\eta })=\mathrm{0.1}$ m s⁻¹, and $v^{\mathrm{s}}(z=\mathit{\eta })=\mathrm{0}$ m s⁻¹ for a 1 m resolution vertical grid using the Breivik et al. (2014) function (black dots), the Breivik et al. (2016) function (grey dots), and the finite-volume Breivik et al. (2016) function (black vertical lines). (b) Their continuous counterparts.

As illustrated in Fig. 1, for the typical vertical resolution used in most global models the properties of the discretized Stokes profiles can be very different from their continuous counterparts. Indeed, the 𝒮_B16(z,k_e) function has been considered superior to 𝒮_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 1 m of depth turn out to be larger for 𝒮_B14(z,k_e) compared to 𝒮_B16(z,k_e). In this case, the fast variations of 𝒮_B16(z,k_e) near the surface cannot be represented by the computational vertical grid. A vertical resolution finer than the one currently used in most global ocean models near the surface would be required to properly represent the Stokes drift shear.

2.3 Subgrid-scale physics

2.3.1 Turbulent kinetic energy prognostic equation and boundary conditions

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 ${\mathcal{W}}_{\mathrm{VF}}^z={\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}{\partial }_{z}u+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}{\partial }_{z}v$ in the hydrostatic relation (Eq. 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 as ${\mathit{\varphi }}^{\prime }=〈\mathit{\varphi }〉-\mathit{\varphi },(\mathit{\varphi }=p,\mathit{\rho },u,v)$, associated with the ${\mathcal{W}}_{\mathrm{VF}}^z$ term are

$$({\mathcal{W}}_{\mathrm{VF}}^z{)}^{\prime }={\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}{\partial }_z{u}^{\prime }+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}{\partial }_z{v}^{\prime }.$$

After multiplication by w^′ and averaging, we obtain

$$〈w^{\prime }({\mathcal{W}}_{\mathrm{VF}}^z{)}^{\prime }〉={\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}{\partial }_z〈u^{\prime }{w}^{\prime }〉+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}{\partial }_z〈v^{\prime }{w}^{\prime }〉-{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}〈u^{\prime }{\partial }_z{w}^{\prime }〉-{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}〈v^{\prime }{\partial }_z{w}^{\prime }〉,$$

where the last two terms on the right-hand side cancel, with similar terms appearing when forming the equations for $〈u^{\prime }{\partial }_t{u}^{\prime }〉$ and $〈v^{\prime }{\partial }_t{v}^{\prime }〉$ (see Eqs. 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^{\mathrm{s}}{\partial }_z〈u^{\prime }{w}^{\prime }〉$ and $v^{\mathrm{s}}{\partial }_z〈v^{\prime }{w}^{\prime }〉$, which can be further rewritten as

$$\begin{array}{rl}& {\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}{\partial }_z〈u^{\prime }{w}^{\prime }〉=-〈u^{\prime }{w}^{\prime }〉{\partial }_z{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}+{\partial }_z\left({\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}〈u^{\prime }{w}^{\prime }〉\right),\\ & {\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}{\partial }_z〈v^{\prime }{w}^{\prime }〉=-〈v^{\prime }{w}^{\prime }〉{\partial }_z{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}+{\partial }_z\left({\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}〈v^{\prime }{w}^{\prime }〉\right).\end{array}$$

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 that is usually parameterized as $-K_e{\partial }_{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

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\partial }_{t}e=& {\displaystyle \frac{A^{\mathrm{vm}}}{e_{\mathrm{3}}^{\mathrm{2}}}}\left[({\partial }_{k}u{)}^{\mathrm{2}}+({\partial }_{k}v{)}^{\mathrm{2}}+\left({\partial }_{k}u\right)\left({\partial }_k{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}\right)+\left({\partial }_{k}v\right)\left({\partial }_k{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}\right)\right]\\ & -A^{\mathrm{vt}}{N}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k\left[{\displaystyle \frac{A^{\mathrm{vm}}}{e_{\mathrm{3}}}}{\partial }_{k}e\right]-c_{\mathit{ϵ}}{\displaystyle \frac{e^{\mathrm{3}/\mathrm{2}}}{l_{\mathit{\epsilon }}^{\mathrm{2}}}}\end{array},\end{array}$$

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_{\mathit{\epsilon }}\approx \mathrm{1}/\sqrt{\mathrm{2}}$). Once the value of e is know, eddy diffusivity and viscosity are given by

$$A^{\mathrm{vm}}=C_{\mathrm{m}}{l}_{\mathrm{m}}\sqrt{e},A^{\mathrm{vt}}=A^{\mathrm{vm}}/\mathrm{Prt},$$

with Prt the Prandtl number (see Sect. 10.1.3 in Madec, 2012, for the detailed computation of Prt), l_m a mixing length scale, and C_m a constant.

In addition to the modification of the shear production term in the TKE equation, the wave will affect the surface boundary condition 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,

$$\begin{array}{}\text{(16)}& {\left({\displaystyle \frac{A^{\mathrm{vm}}}{e_{\mathrm{3}}}}{\partial }_{k}e\right)}_{z=z_{\mathrm{1}}}=-{\mathit{\rho }}_{\mathrm{0}}g\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{S}_{\mathrm{oce}}\mathrm{d}\mathit{\omega }\mathrm{d}\mathit{\theta }={\mathrm{\Phi }}_{\mathrm{oce}},\end{array}$$

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, 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 topmost grid cell (i.e., at z=z₁). This amounts to interpreting the half-grid cell at the top as a constant flux layer, which is consistent with 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, respectively estimating 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_{\mathrm{up}}\left(z\right)=l_{\mathrm{dwn}}\left(z\right)=\sqrt{\mathrm{2}e\left(z\right)/N^{\mathrm{2}}\left(z\right)}$. 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 water column such as the thermocline. This limitation amounts to controlling the vertical gradients of l_up(z) and l_dwn(z) such that they are not larger that the variations of depth (Madec, 2012)

$${\partial }_k\left|l_{\cdot }\right|\le e_{\mathrm{3}},l_{\cdot }=l_{\mathrm{up}},l_{\mathrm{dwn}}$$

Then, the dissipative and mixing length scales are given by $l_{\mathrm{m}}=\sqrt{l_{\mathrm{up}}{l}_{\mathrm{dwn}}}$ and $l_{\mathit{\epsilon }}=min\left(l_{\mathrm{up}},l_{\mathrm{dwn}}\right)$. Following Redelsperger et al. (2001) (their Sect. 4.2.3), a boundary condition consistent with the Monin–Obukhov similarity theory for the length scale l_dwn (while l_up necessitates only a bottom boundary condition) is

$$l_{\mathrm{dwn}}(z=\mathit{\eta })=\mathit{\kappa }{\displaystyle \frac{(C_{\mathrm{m}}{c}_{\mathit{\epsilon }}{)}^{\mathrm{1}/\mathrm{4}}}{C_{\mathrm{m}}}}{z}_{\mathrm{0}},$$

with κ the von Karman constant and C_m and c_ε the constant parameters in the TKE closure. The surface roughness length z₀ can be directly estimated from the significant wave height provided by the wave model as z₀=1.6H_s (Rascle et al., 2008, their Eq. 5), which provides a proxy for the scale of the breaking waves. Note that in our study, no explicit parameterization of the 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 the effects of wave-related mixing processes such as Langmuir turbulence. For this reason, it is not used in the present simulations.

2.3.2 Langmuir turbulence parameterization

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 (15). P_LC is defined as

$$P_{\mathrm{LC}}={\displaystyle \frac{w_{\mathrm{LC}}^{\mathrm{3}}}{d_{\mathrm{LC}}}},$$

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

$$\begin{array}{rl}& w_{\mathrm{LC}}=\left\{\begin{array}{ll}{c}_{\mathrm{LC}}\Vert {\widehat{\mathit{u}}}_{\mathrm{LC}}^{\mathrm{s}}\Vert \sin \left(-\frac{\mathit{\pi }z}{d_{\mathrm{LC}}}\right),& \text{if}-z\le d_{\mathrm{LC}}\\ \mathrm{0},& \text{otherwise}\end{array}\right.,\\ & -\underset{-d_{\mathrm{LC}}}{\overset{\mathit{\eta }}{\int }}{N}^{\mathrm{2}}\left(z\right)z\mathrm{d}z={\displaystyle \frac{\Vert {\widehat{\mathit{u}}}_{\mathrm{LC}}^{\mathrm{s}}{\Vert }^{\mathrm{2}}}{\mathrm{2}}},\end{array}$$

where $\Vert {\widehat{\mathit{u}}}_{\mathrm{LC}}^{\mathrm{s}}\Vert$ is the portion of the surface Stokes drift contributing to Langmuir cell intensity and c_LC a constant parameter. In the absence of information about the wave field it is generally assumed that $\Vert {\widehat{\mathit{u}}}_{\mathrm{LC}}^{\mathrm{s}}\Vert \propto \sqrt{\Vert \mathit{\tau }\Vert }$. 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 $\Vert {\widehat{\mathit{u}}}_{\mathrm{LC}}^{\mathrm{s}}\Vert$ via

$$\begin{array}{}\text{(17)}& \Vert {\widehat{\mathit{u}}}_{\mathrm{LC}}^{\mathrm{s}}\Vert =max\left\{{\mathit{u}}^{\mathrm{s}}\left(\mathit{\eta }\right)\cdot {\mathit{e}}_{\mathit{\tau }},\mathrm{0}\right\},\end{array}$$

with e_τ the unit vector in the wind stress direction. The difference between the surface Stokes drift $\Vert {\mathit{u}}^{\mathrm{s}}\left(\mathit{\eta }\right)\Vert$ and $\Vert {\widehat{\mathit{u}}}_{\mathrm{LC}}^{\mathrm{s}}\Vert$ given by Eq. (17) is shown in Fig. 2 and compared to the usual parameterization of $\Vert {\mathit{u}}^{\mathrm{s}}\left(\mathit{\eta }\right)\Vert$ as $\mathrm{0.377}\sqrt{\Vert {\mathit{\tau }}^{\mathrm{oce}}\Vert /{\mathit{\rho }}_{\mathrm{0}}}$ in the uncoupled case (see Madec, 2012). The modulation of $\Vert {\widehat{\mathit{u}}}_{\mathrm{LC}}^{\mathrm{s}}\Vert$ depending on the wind stress orientation significantly reduces the input of the surface Stokes drift contributing to Langmuir cell intensity, especially in the Southern Ocean, while other regions are less affected. Finally, a value for the parameter c_LC must be chosen. Based on single-column experiments detailed in Appendix B, we find that parameter values in the range 0.15–0.3 provide satisfactory results compared to the LESs of Noh et al. (2016) and will be considered for the numerical experiments discussed later in Sect. 4.2.

Figure 2 Annual average of the surface Stokes drift module $\Vert {\mathit{u}}^{\mathrm{s}}\left(\mathit{\eta }\right)\Vert$ (m s⁻¹) (a), the portion of the Stokes drift aligned with the wind, as given in Eq. (17) (b), and the surface Stokes drift as parameterized by $\mathrm{0.377}\sqrt{\Vert {\mathit{\tau }}^{\mathrm{oce}}\Vert /{\mathit{\rho }}_{\mathrm{0}}}$ in the uncoupled case (c).

While the Axell (2002) parameterization was already implemented in NEMO, there are three major novelties in our implementation: (i) the online coupled strategy allows us to use the surface Stokes drift directly delivered by the wave model instead of the original value empirically estimated from the wind speed (e.g., 1.6 % of the 10 m wind). (ii) We only considered the component of the Stokes drift aligned with the wind, and (iii) based on a series of single-column simulations (see Appendix B) the coefficient c_LC evaluated as 0.15 by Axell (2002) is set to a 0.3 value. Those changes, together with the new surface boundary condition for the TKE equation, lead to a deeper penetration of the TKE inside the mixed layer and, as shown in Sect. 4.2.3, greatly improved the MLD distribution.

3 Modeling system and coupling strategy

Our coupled model is based on the NEMO oceanic model, the WW3 wave model, and the OASIS library for data exchange and synchronization between the two components.

3.1 Numerical models and coupling infrastructure

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 Eq. 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 wave propagation, 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 Eqs. (7) and (13). Moreover, the modifications to the standard NEMO one-equation TKE closure scheme are given in Sect. 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 only the details sufficient to understand the coupling of waves with the oceanic model here, and an exhaustive description of WW3 is given by WAVEWATCH III® Development Group (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 (the subscript w is used here to avoid confusion with previously introduced notations):

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\partial }_t{N}_{\mathrm{w}}+{\partial }_{\mathit{\varphi }}\left(\dot{\mathit{\varphi }}{N}_{\mathrm{w}}\right)+{\partial }_{\mathit{\lambda }}\left(\dot{\mathit{\lambda }}{N}_{\mathrm{w}}\right)+{\partial }_{k_{\mathrm{w}}}\left({\dot{k}}_{\mathrm{w}}{N}_{\mathrm{w}}\right)\\ & +{\partial }_{{\mathit{\theta }}_{\mathrm{w}}}\left({\dot{\mathit{\theta }}}_{\mathrm{w}}{N}_{\mathrm{w}}\right)={\displaystyle \frac{S}{\mathit{\sigma }}},\end{array}\end{array}$$

where λ is longitude, ϕ is latitude, and S is the net spectral source term that includes the sum of the rate of change of the surface 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 wave–ocean coupling. Indeed, as shown earlier in Eq. (16), the S_oce term is used to compute the TKE flux transmitted to the ocean, and the S_in term enters the computation of the wave-supported stress. They are computed here following Ardhuin et al. (2010b). In Eq. (18), the dot variables correspond to a propagation speed given by the following.

$$\begin{array}{}\text{(19)}& {\displaystyle }\dot{\mathit{\varphi }}=\left(c_{\mathrm{g}}\cos {\mathit{\theta }}_{\mathrm{w}}+{\left.v\right|}_{z=\mathit{\eta }}\right)R^{-\mathrm{1}}\text{(20)}& {\displaystyle }\dot{\mathit{\lambda }}=\left(c_{\mathrm{g}}\sin {\mathit{\theta }}_{\mathrm{w}}+{\left.u\right|}_{z=\mathit{\eta }}\right){\left(R\cos \mathit{\varphi }\right)}^{-\mathrm{1}}\text{(21)}& {\displaystyle }\begin{array}{rl}{\dot{\mathit{\theta }}}_{\mathrm{w}}=& c_{\mathrm{g}}\sin {\mathit{\theta }}_{\mathrm{w}}\tan \mathit{\varphi }{R}^{-\mathrm{1}}+\sin {\mathit{\theta }}_{\mathrm{w}}{\displaystyle \frac{\partial {\mathit{\omega }}_{\mathrm{w}}}{\partial \mathit{\varphi }}}\\ & -{\displaystyle \frac{\cos {\mathit{\theta }}_{\mathrm{w}}}{\cos \mathit{\varphi }}}{\displaystyle \frac{\partial {\mathit{\omega }}_{\mathrm{w}}}{\partial \mathit{\lambda }}}{\left(k_{\mathrm{w}}R\right)}^{-\mathrm{1}}\end{array}\text{(22)}& {\displaystyle }{\dot{k}}_{\mathrm{w}}=-{\displaystyle \frac{\partial \mathit{\sigma }}{\partial H}}{\displaystyle \frac{\mathit{k}}{k_{\mathrm{w}}}}\mathbf{\cdot }\mathbf{\nabla }D-\mathit{k}\mathbf{\cdot }\mathbf{\nabla }{\mathit{u}}_{\mathrm{h}}(z=\mathit{\eta })\end{array}$$

Here, R is the Earth's radius, ${\mathit{u}}_{\mathrm{h}}(z=\mathit{\eta })=({\left.u\right|}_{z=\mathit{\eta }},{\left.v\right|}_{z=\mathit{\eta }})$ represents the surface currents provided by the ocean model, c_g is the group velocity, ω_w the absolute radian frequency, and H the mean water depth. Equation (18) is solved for each spectral component (k_w,θ_w) coupled by the advection and source terms. Equations (19)–(22) show how the oceanic currents affect the advection of the wave action density; there are also indirect effects via the source term (Ardhuin et al., 2009).

The coupler: OASIS-MCT

The practical coupling between NEMO and WW3 has been implemented using the OASIS-MCT (Valcke, 2013; Craig et al., 2017) software primarily developed for use in multicomponent 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 whereby the matrix containing the mapping weights is computed offline once for all. In practice, after compiling OASIS-MCT, the resulting library is linked to the component models so that they have access 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.

3.2 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 the atmosphere is consumed by the wave field and contributes to the growing waves (the so-called wave-supported stress); 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 m wind ${\mathit{u}}_{\mathrm{10}}^{\mathrm{atm}}$ is sent to the wave model, which internally computes the dimensionless Charnock parameter α_ch characterizing the sea surface roughness (Charnock, 1955; Janssen, 2009). This information is used by the wave model to compute its own atmospheric wind stress ${\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{atm}}$ assuming neutral stratification, i.e., ${\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{atm}}={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{DN}}\left({\mathit{\alpha }}_{\mathrm{ch}}\right)\Vert {\mathit{u}}_{\mathrm{10}}^{\mathrm{atm}}\Vert {\mathit{u}}_{\mathrm{10}}^{\mathrm{atm}}$ 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 ${\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{oce}}$. Using the latest available values of α_ch, ${\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{atm}}$, ${\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{oce}}$, and ${\mathit{u}}_{\mathrm{10}}^{\mathrm{atm}}$, the oceanic model computes an atmospheric wind stress τ^atm using its own bulk formulation, and the local value of the momentum flux going into the water column is diagnosed as

$$\begin{array}{}\text{(23)}& {\mathit{\tau }}^{\mathrm{oce}}={\mathit{\tau }}^{\mathrm{atm}}-\left({\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{atm}}-{\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{oce}}\right),\end{array}$$

where the τ_ww3 quantities are interpolated from the wave grid to the oceanic grid. In NEMO, the wind stress is computed using the IFS (Integrated Forecasting System: https://www.ecmwf.int/en/forecasts/documentation-and-support/changes-ecmwf-model/ifs-documentation, last access: 2 July 2020) bulk formulation such as implemented in the AeroBulk (https://github.com/brodeau/aerobulk, last access: 2 July 2020) package (Brodeau et al., 2016). In particular, the roughness length that enters the definition of the drag coefficient is a function of the Charnock parameter α_ch,

$$z_{\mathrm{0}}={\mathit{\alpha }}_{\mathrm{ch}}{\displaystyle \frac{u_{\star }^{\mathrm{2}}}{g}}+{\mathit{\alpha }}_m{\displaystyle \frac{\mathit{\nu }}{u_{\star }}},$$

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 ${\mathit{\alpha }}_{\mathrm{ch}}^{\mathrm{0}}=\mathrm{0.018}$. In our implementation, the momentum fluxes are computed using the absolute wind ${\mathit{u}}_{\mathrm{10}}^{\mathrm{atm}}$ at 10 m rather than the relative wind ${\mathit{u}}_{\mathrm{10}}^{\mathrm{atm}}-{\mathit{u}}_{\mathrm{h}}(z=\mathit{\eta })$. Indeed, several recent studies have emphasized that the use of relative winds is relevant only when a full 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 us 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 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. However, it was necessary in practice since the WW3 results were very sensitive to the bulk formulation, and at the same time it was not conceivable to use the WW3 bulk formulation to force the ocean model because the latter ignores the effect of stratification in the atmospheric surface layer. Previous implementations in NEMO (e.g., Breivik et al., 2015; Alari et al., 2016; Staneva et al., 2017; Law Chune and Aouf, 2018; Wu et al., 2019) assumed that the wave field only acts on the norm of τ^atm and not on its orientation. Instead of Eq. (23), the atmospheric wind stress was corrected as follows.

$${\mathit{\tau }}^{\mathrm{oce}}={\mathit{\tau }}^{\mathrm{atm}}\left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{oce}}}{{\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{atm}}}}\right)$$

However, this approach potentially leads to artificially large values of τ^oce when ${\mathit{\tau }}_{\mathrm{ww}\mathrm{3}}^{\mathrm{atm}}$ is small, and it does not take into account the slight change in τ^oce direction induced by the waves.

3.3 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, meaning that no 3D arrays are exchanged through the coupler. All 2D interpolations 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\in \mathbb{N},n_t\ge \mathrm{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 does not include an explicit coupling between waves and sea ice, while it is known that waves lead to ice breakup, pancake ice formation, and associated enhancement of both freezing and melting; 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).

Table 1 Variables exchanged between NEMO (O) and WW3 (W) via the OASIS-MCT coupler. The 10 m wind ${\mathit{u}}_{\mathrm{10}}^{\mathrm{atm}}$ is interpolated online by WW3 and does not go through the OASIS-MCT coupler.

4 Global $\mathrm{1}/\mathrm{4}{}^{\circ }$ coupled wave–ocean simulations

4.1 Experimental setup and experiments

4.1.1 The global coupled ORCA25 configuration

The wave hindcasts presented here are all based on the WW3 model in its version 6.02 configured with a single grid at 0.5^∘ resolution in longitude and latitude. A spectral grid with 24 directions and 31 frequencies is exponentially spaced over the interval $[f_{min},f_{max}]$ with $f_{min}=\mathrm{0.037}$ Hz and $f_{max}=\mathrm{0.7}$ Hz. A one-step monotonic third-order coupled space–time advection scheme (also called the 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 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 errors are usually compensated for by a proper adjustment of the dissipation source term (Banner and Young, 1994; Ardhuin et al., 2007). As mentioned earlier in Sect. 3.2, the model was run with 10 m winds, without any air–sea stability correction. No wave measurements were assimilated in the model, but the stand-alone wave model was developed based on spectral buoy and synthetic aperture radar (SAR) data (Ardhuin et al., 2010b) and calibrated against altimeter data by adjusting the wind–wave 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 $\mathrm{1}/\mathrm{4}{}^{\circ }$ 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. The LIM3 sea ice model is used for the sea ice dynamics and thermodynamics (Rousset et al., 2015). The vertical mixing coefficients are obtained from the one-equation TKE scheme described in Sect. 2.3, and the convective processes are mimicked using an enhanced vertical diffusion parameterization that increases vertical diffusivity to 10 m² s⁻¹ at which static instability occurs. Water density is computed from temperature and salinity through the use of a polynomial formulation of the UNESCO (1983) nonlinear equation of state (Roquet et al., 2015). The vector-invariant form of momentum advection uses Arakawa and Lamb (1981) for the vorticity and a specific formulation to control the Hollingsworth instability (Ducousso et al., 2017). Momentum lateral viscosity is biharmonic and acts along geopotential surfaces. It is set to a value of 1.5×10¹¹ m⁴ s⁻¹ at the Equator and varies proportionally to Δx³ away from the Equator. Advection of tracers is performed with a flux-corrected transport (FCT) scheme (Lévy et al., 2001), and lateral diffusion of tracers is harmonic and acts along an iso-neutral surface. It is set to a value of 300 m² s⁻¹ at the Equator, which varies proportionally to Δx. The bottom friction is nonlinear and the lateral boundary condition is free-slip. In this setup, the baroclinic time step is set to Δt_nemo=900 s and a barotropic time step 30 times smaller. Compared to the standard uncoupled ORCA025 configuration, the additional computational cost associated with WW3 and the exchanges through the coupler is about 20 %.

4.1.2 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 ERA-Interim mean states for rainfall as well as shortwave and longwave radiative fluxes are consistent with satellite observations from the Remote Sensing Systems (RSS) Passive Microwave Water Cycle (PMWC) product (Hilburn, 2009) and GEWEX SRB 3.1 data (https://gewex-srb.larc.nasa.gov/, last access: 2 July 2020). Momentum and heat turbulent surface fluxes are computed using the IFS bulk formulation from the AeroBulk package (Brodeau et al., 2016) using air temperature and humidity at 2 m, mean sea level pressure, and 10 m winds.

4.1.3 Sensitivity experiments and objectives

Sensitivity experiments have been conducted to check the proper implementation of various components of the present coupled modeling system. For the sake of clarity, our developments are split into four components: (i) the modification of the wind stress by waves through the Charnock parameter and the inclusion of wave-supported stress, (ii) the modifications of the NEMO governing equations through the Stokes–Coriolis, vortex force, and wave-induced surface pressure terms, (iii) the addition of a Langmuir turbulence parameterization, and (iv) the modifications to the TKE scheme. As summarized in Table 2, sensitivity experiments are designed in such a way to incrementally increase the level of complexity and test the effect of each component. The No_CPL experiment corresponds to the classical NEMO setup in which the 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 follows: $e(z=\mathit{\eta })=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{15.8}{\mathit{\alpha }}_{\mathrm{CB}}{)}^{\mathrm{2}/\mathrm{3}}\frac{\Vert {\mathit{\tau }}^{\mathrm{atm}}\Vert }{{\mathit{\rho }}_{\mathrm{0}}}$ with α_CB=100. Based on the results of Mellor and Blumberg (2004) we expect that in the uncoupled case the nature of the boundary condition (i.e., Dirichlet vs. Neumann) does not significantly impact numerical solutions³. The WS_CPL experiment is identical as No_CPL except that the wave coupling is introduced within the wind stress computation, as described in Sect. 3.2. The ST_CPL experiment is as WS_CPL except that all terms relative to the Stokes drift described in Sect. 2.1 are added in NEMO. TKE_CPL corresponds to ST_CPL but with the modified TKE scheme described in Sect. 2.3.1. All_CPL (1 and 2) experiments are like TKE_CPL but with a fully modified TKE scheme including the Langmuir cell parameterization described in Sect. 2.3.2. All those simulations have been performed for 2 years (2013–2014), during which 2013 is spin-up and only 2014 is analyzed. We considered 2 years sufficient 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.

Table 2 Various model configurations analyzed in Sect. 4.2.

4.2 Numerical results

4.2.1 Wave 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 the wave patterns more closely. 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 the uncoupled case (${\mathit{\alpha }}_{\mathrm{ch}}^{\mathrm{0}}=\mathrm{0.018}$) are shown. As expected, the Charnock parameter tends to be stronger in the area where the waves are higher. Generally an increase in the Charnock parameter is observed in the northern and southern basin, while there is a net decrease in α_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 ${\mathit{\alpha }}_{\mathrm{ch}}^{\mathrm{0}}$ are very latitudinal with very few longitudinal variations.

Figure 3 (a, c) Seasonal averages of significant wave height (in meters) for January–February–March (JFM, panel a) and July–August–September (JAS, panel c). (b, d) Seasonal average of the difference between the Charnock parameter as computed by the wave model and the default value ${\mathit{\alpha }}_{\mathrm{ch}}^{\mathrm{0}}=\mathrm{0.018}$ for JFM (b) and JAS (d).

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 wind stress $\Vert {\mathit{\tau }}^{\mathrm{oce}}\Vert$ between those two cases is primarily due to the use of different Charnock parameters and the inclusion of the wave-supported stress.

Figure 4 (a) Drag coefficient (C_D) as a function of the 10 m wind speed $\Vert {\mathit{u}}_{\mathrm{10}}^{\mathrm{atm}}\Vert$ and (b) wind stress norm $\Vert {\mathit{\tau }}^{\mathrm{oce}}\Vert$ as a function of $\Vert {\mathit{u}}_{\mathrm{10}}^{\mathrm{atm}}\Vert$ (black curves represent the mean value, while the vertical bars represent the standard deviation).

Figure 4a illustrates that the Charnock parameter mostly affects the drag coefficient C_D, and hence the surface wind stress, for large winds. The ocean–wave coupling does not lead to appreciable differences in the drag coefficient C_D for wind speeds lower than 8 m s⁻¹. 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. Figure 4b shows how the wind stress is modified by this increase in 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 in wind stress associated with the increased drag coefficient is always larger than the decrease associated with 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 seen in our simulations are very localized in time and space and it is thus difficult to conclude on their overall effect on upper-ocean dynamics such as Ekman pumping and the surface currents.

Figure 5 Wind stress difference $\Vert {\mathit{\tau }}^{\mathrm{oce}}\Vert -\Vert {\mathit{\tau }}^{\mathrm{atm}}\Vert$ (N m⁻²) due to the correction made for growing waves for the WS_CPL experiment as a function of the 10 m wind speed.

4.2.2 Wave impact on surface TKE injection

As described in Sect. 2.3, in the ocean–wave 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 in which a Dirichlet condition is imposed. We aim here to assess the impact on the order of magnitude of the near-surface TKE. Since the Neumann boundary condition is applied at the center of the topmost grid box (i.e., approximately at 50 cm of depth), we compare in Fig. 6 the TKE value at 1 m of depth between the coupled (All_CPL2) and the uncoupled (No_CPL) case. Positive values mean that near-surface TKE is larger in the coupled simulation.

Figure 6 Seasonal difference of 1 m depth turbulent kinetic energy (m² s⁻²) between the coupled case (All_CPL2) and the uncoupled case (No_CPL). (a) January, February, and March (JFM); (b) July, August, and September (JAS).

It shows an almost homogeneous increase in the TKE (up to more than 100 %) in the extratropical areas. While low seasonal variability in the extratropical areas is visible in Fig. 6, a spatial averaging by hemisphere (Fig. 7) highlights seasonal variability with a strong increase in both the near-surface TKE value and the TKE difference between the two experiments during winter. In Figs. 6 and 7 (and also in the remainder of the paper), the spatial averaging is made between 25 and 60^∘ S in the Southern Hemisphere and between 25 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 in the surface TKE injection associated with waves is expected to contribute to an overall increase in mixed layer depth provided that the mixing length diagnosed by the turbulent closure scheme allows for the effective propagation of this additional TKE deeper in the mixed layer.

Figure 7 Spatially averaged turbulent kinetic energy (m² s⁻²) at 1 m of depth over (a) the Southern Hemisphere and (b) the Northern Hemisphere.

4.2.3 Wave 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. Figure 8 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_{\mathrm{MLD}}^{\mathrm{noCPL}}-h_{\mathrm{MLD}}^{\mathrm{CPL}})/h_{\mathrm{MLD}}^{\mathrm{noCPL}}$ with h_MLD considered negative downward). It shows a significant deepening of the mixed layer at high latitudes in the coupled case with only a very few localized mixed layer shallowing up to 60 %, mainly in the Southern Hemisphere.

Figure 8 (a, b) Seasonally averaged MLD differences (All_CPL2-No_CPL) relative to the uncoupled simulation No_CPL. Red corresponds to a deeper MLD for All_CPL2.

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 dataset from de Boyer Montégut et al. (2004). The MLD has been computed as the depth at which the density is 3 % smaller that the density at 10 m as in de Boyer Montégut et al. (2004). Figure 9 represents the spatially averaged MLD; 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. 9a), 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. 9b) the situation is rather different.

Figure 9 Spatially averaged MLD for (a) the Northern Hemisphere and (b) the Southern Hemisphere.

From January to July, the deepening of the 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, thereby increasing the bias between the data and model. Since mesoscale activity makes 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 the year 2014 (panel a) and during summer only (panel b). In both cases 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 the 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. 10b), the probability density function (PDF) in the coupled case matches the one computed from ARGO data almost perfectly. Despite the fact that we did not activate the 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.

Figure 10 Mixed layer depth probability density function for (a) the full 2014 year and (b) summer 2014.

To better understand which components of the wave–ocean coupling are responsible for this improvement, the summer PDF in the Southern Hemisphere has been computed for each of the experiments described in Table 2. Results are shown in Fig. 11. First of all, it can be seen that all the wave–ocean interactions described in previous sections lead to an improvement in terms of mixed layer depth distribution compared to the uncoupled 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.

Figure 11 Mixed layer depth probability density function in the Southern Hemisphere during summer months. The details of each experiment can be found in Table 2.

4.2.4 Wave 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). Figure 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. 12b) the summer warm bias is reduced by half in the coupled simulation and a slight warming occurs during the winter. While the summer surface 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.

Figure 12 Time series of the spatially averaged sea surface temperature (^∘C); (a) Northern Hemisphere and (b) Southern Hemisphere.

To better characterize the wave impact on the SST, we show in Fig. 13a the difference in terms 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 cold bias (see Fig. 13b), especially in the Pacific Ocean, and the warm bias in the North Pacific is significantly reduced.

Figure 13 (a) Annual average of the differences between No_CPL and OSTIA sea surface temperatures (^∘C) for the year 2014 (positive when the model is warmer). (b) Annual average of the difference between All_CPL2 and No_CPL (positive when All_CPL2 is warmer).

As already noticed by Law Chune and Aouf (2018) the warming in the equatorial and tropical regions mainly results from a lower wind stress caused by a value of the Charnock parameter lower than the value used in the uncoupled case (see Fig. 3b, d). A consequence is a decrease in the drag coefficient, leading to smaller turbulent exchange coefficients and reducing the heat flux. As mentioned above, in extratropical regions, some warm bias tends to be partially reduced by the extra mixing induced by the waves at high latitude and/or 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 Hovmöller diagram like the ones shown in Fig. 14. For instance, it can be seen that the summer warm bias in the Northern Hemisphere (Fig. 14a) coincides well with the cooling induced by the coupling with waves (Fig. 14b). Similarly we can also observe a warming in the tropical and equatorial regions (Fig. 14b) corresponding to the cold bias seen in Fig. 14a. In the southern extratropical region, a summer cooling is observed. It is induced by the wave coupling, whereas Fig. 14a shows a slight warm bias. During winter we can observe a warming in Fig. 14b north of 60^∘ S, which again partially corresponds to a cold bias in Fig. 14a.

Figure 14 Hovmöller diagram of the longitudinally averaged sea surface temperature (^∘C) differences between (a) No_CPL and OSTIA and (b) between All_CPL2 and No_CPL.

Figure 15 Time series of the spatially averaged surface kinetic energy (m² s⁻²) for (a) the Northern Hemisphere and (b) the Southern Hemisphere.

4.2.5 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. Whatever the hemisphere there is a net decrease in surface KE (up to 20 % in the south) when a coupling with the waves is included. This decrease in surface kinetic energy reflects a decrease in surface current magnitudes. 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 current magnitudes 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 in velocities in the whole boundary layer as previously shown by Rascle et al. (2008) (orange lines in Fig. 16). The inclusion of 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.

Figure 16 Zonally averaged zonal (a) and meridional (b) currents (m s⁻¹) between 60 and 25^∘ S as a function of depth (m) for the simulations described in Table 2.

5 Conclusions

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 detail: (i) the 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) to include wave effects and a parameterization of Langmuir turbulence; (iii) development of a coupling interface based on the OASIS-MCT software for the exchange of data between the two models; and (iv) tests of our developments on a realistic global configuration at $\mathrm{1}/{\mathrm{4}}^{\circ }$ for the ocean coupled to a $\mathrm{1}/{\mathrm{2}}^{\circ }$ resolution wave model. Compared to an ocean-only simulation, the coupling with a wave model (with a resolution twice as coarse as the oceanic model) leads to an additional computational cost of about 20 %.

Following McWilliams et al. (2004) and Ardhuin et al. (2008), in the weak vertical current shear limit, the wave-induced terms implemented in NEMO include the Stokes–Coriolis force, the vortex force, Stokes advection in tracer and continuity equations, and a wave-induced surface pressure term. The prognostic equation for TKE also includes an additional forcing term associated with the Stokes drift vertical shear and various modifications of its boundary condition described in Sect. 2.3.

The development of a coupling infrastructure based on OASIS-MCT has several advantages as it allows for an efficient data exchange (including the treatment of nonconformities 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 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 the ETAU parameterization). This parameterization that allows users to empirically propagate the surface TKE at depth using a prescribed shape function is a pragmatic way to cure the shallow mixed 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), and Staneva et al. (2017) have used the ETAU parameterization in their setup. However, as suggested by Breivik et al. (2015), we can speculate that such a 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 one-equation TKE scheme described in Sect. 2.3 is more physically justifiable than the ETAU parameterization and requires much less parameter tuning.

The numerical experiments based on the ORCA25 configuration discussed in Sect. 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 middle and high latitudes, whereas accounting for the portion of the wind stress consumed by the waves has a small impact (in our experiments it leads to a maximum 2 % decrease in 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, 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. Since the comparison with ARGO data can be tricky due to the scarcity of data, we looked at the results in terms of mixed layer depth (MLD) probability density functions. This allowed us 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 the Langmuir cell parameterization of Axell (2002), which is consistent with recent results obtained by Reichl et al. (2016) and Ali et al. (2019) using a K-profile parameterization (KPP) closure scheme.

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 better agreement with the OSTIA SST analysis. More generally, although the global SST biases are not totally compensated for, they tend to be reduced when considering the effect of waves (see Sect. 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 $\mathrm{1}/{\mathrm{4}}^{\circ }$ 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 (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–current interactions in the forthcoming NEMO official release. We can speculate that the ocean–wave 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 two-way coupling compared to a one-way coupling on the oceanic and wave solution must be investigated in the future. Indeed, the improvements of the quality of surface wave simulations associated with 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–100 km) 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.

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 preferred to the vector-invariant form presented in Eqs. (7) and (8). For simplicity, we consider just the i component in horizontal curvilinear coordinates and the z coordinate in the vertical (results will be extended to the j component and to a generalized vertical coordinate). Consistently with the notations of Madec (2012), e₁ and e₂ are the horizontal scale factors. We denote ${\mathrm{A}}_{\mathrm{v}}^u$ the extra term needed to guarantee equivalence between the flux formulation and the vector-invariant form. ${\mathrm{A}}_{\mathrm{v}}^u$ is defined such that

$$\mathrm{\nabla }\cdot \left({\mathit{u}}^{\mathrm{s}}u\right)+{\mathrm{A}}_{\mathrm{v}}^u=-\mathit{\zeta }{v}^{\mathrm{s}}+{\displaystyle \frac{w^{\mathrm{s}}}{e_{\mathrm{3}}}}{\partial }_{k}u.$$

Since $\mathrm{\nabla }\cdot {\mathit{u}}^{\mathrm{s}}=\mathrm{0}$, we have $\mathrm{\nabla }\cdot \left({\mathit{u}}^{\mathrm{s}}u\right)={\mathit{u}}^{\mathrm{s}}\cdot \mathrm{\nabla }u$ and thus

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{e}_{\mathrm{1}}{e}_{\mathrm{2}}{\mathrm{A}}_{\mathrm{v}}^u=& -v^{\mathrm{s}}\left[{\partial }_i\left(e_{\mathrm{2}}v\right)-{\partial }_j\left(e_{\mathrm{1}}u\right)\right]+{\displaystyle \frac{e_{\mathrm{1}}{e}_{\mathrm{2}}}{e_{\mathrm{3}}}}{w}^{\mathrm{s}}{\partial }_{k}u\\ & -\left[e_{\mathrm{2}}{u}^{\mathrm{s}}{\partial }_{i}u+e_{\mathrm{1}}{v}^{\mathrm{s}}{\partial }_{j}u+{\displaystyle \frac{e_{\mathrm{1}}{e}_{\mathrm{2}}}{e_{\mathrm{3}}}}{w}^{\mathrm{s}}{\partial }_{k}u\right]\\ & =-v^{\mathrm{s}}\left[v{\partial }_i{e}_{\mathrm{2}}-u{\partial }_j{e}_{\mathrm{1}}+e_{\mathrm{2}}{\partial }_{i}v-e_{\mathrm{1}}{\partial }_{j}u\right]\\ & -e_{\mathrm{2}}{u}^{\mathrm{s}}{\partial }_{i}u-e_{\mathrm{1}}{v}^{\mathrm{s}}{\partial }_{j}u.\end{array}\end{array}$$

Hence,

$${\mathrm{A}}_{\mathrm{v}}^u=\underset{\text{Metric term on Stokes drift}}{\underbrace{-{\displaystyle \frac{v^{\mathrm{s}}}{e_{\mathrm{1}}{e}_{\mathrm{2}}}}\left(v{\partial }_i{e}_{\mathrm{2}}-u{\partial }_j{e}_{\mathrm{1}}\right)}}-\underset{\text{Additional term}}{\underbrace{\left({\displaystyle \frac{u^{\mathrm{s}}}{e_{\mathrm{1}}}}{\partial }_{i}u+{\displaystyle \frac{v^{\mathrm{s}}}{e_{\mathrm{1}}}}{\partial }_{i}v\right)}}.$$

The same computation for the j component leads to the following equations in generalized vertical coordinates.

$$\begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}& {\partial }_t\left(e_{\mathrm{3}}u\right)=-{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{1}}{e}_{\mathrm{2}}}}\left[{\partial }_i\left(e_{\mathrm{2}}\right(u+{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}\left)u\right)+{\partial }_j\left(e_{\mathrm{1}}\right(v+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}\left)u\right)\right]\\ & +{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k\left(\right(\mathit{\omega }+{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}}\left)u\right)+\left[f+{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{1}}{e}_{\mathrm{2}}}}\left(v{\partial }_i{e}_{\mathrm{2}}-u{\partial }_j{e}_{\mathrm{1}}\right)\right]\\ & (v+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}})+{\displaystyle \frac{{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}}{e_{\mathrm{1}}}}{\left({\partial }_{i}u\right)}_z+{\displaystyle \frac{{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}}{e_{\mathrm{1}}}}{\left({\partial }_{i}v\right)}_z-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}{e}_{\mathrm{1}}}}{\partial }_i(p_{\mathrm{s}}+{\stackrel{\mathrm{̃}}{p}}^J)\\ & -{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}{e}_{\mathrm{1}}}}{\left({\partial }_i{p}_{\mathrm{h}}\right)}_z+{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k〈u^{\prime }{\mathit{\omega }}^{\prime }〉+F^u+{\stackrel{\mathrm{̃}}{F}}^u{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_t\left(e_{\mathrm{3}}v\right)\\ & =-{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{1}}{e}_{\mathrm{2}}}}\left[{\partial }_i\left(e_{\mathrm{2}}\right(u+{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}\left)v\right)+{\partial }_j\left(e_{\mathrm{1}}\right(v+{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}\left)v\right)\right]\\ & +{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k\left(\right(\mathit{\omega }+{\stackrel{\mathrm{̃}}{\mathit{\omega }}}^{\mathrm{s}}\left)v\right)-\left[f+{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{1}}{e}_{\mathrm{2}}}}\left(v{\partial }_i{e}_{\mathrm{2}}-u{\partial }_j{e}_{\mathrm{1}}\right)\right]\\ & (u+{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}})+{\displaystyle \frac{{\stackrel{\mathrm{̃}}{u}}^{\mathrm{s}}}{e_{\mathrm{2}}}}{\left({\partial }_{j}u\right)}_z+{\displaystyle \frac{{\stackrel{\mathrm{̃}}{v}}^{\mathrm{s}}}{e_{\mathrm{2}}}}{\left({\partial }_{j}v\right)}_z-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}{e}_{\mathrm{2}}}}{\partial }_j(p_{\mathrm{s}}+{\stackrel{\mathrm{̃}}{p}}^J)\\ & -{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}{e}_{\mathrm{2}}}}{\left({\partial }_j{p}_{\mathrm{h}}\right)}_z+{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{3}}}}{\partial }_k〈v^{\prime }{\mathit{\omega }}^{\prime }〉+F^v+{\stackrel{\mathrm{̃}}{F}}^v\end{array}$$

Here, ${\left({\partial }_i•\right)}_z$ and ${\left({\partial }_j•\right)}_z$ are derivatives along the z coordinate.

Appendix B: Sensitivity to the c_LC parameter from single-column experiments

Single-column experiments based on Noh et al. (2016) have been performed to study the behavior of the NEMO vertical closure with the Langmuir cell parameterization of Axell (2002). In the Noh et al. (2016) experiments the initial condition is given by

$$u(z,t)=v(z,t)=\mathrm{0},\mathit{\theta }(z,t)=min\left\{T_{\mathrm{0}}-N_{\mathrm{0}}^{\mathrm{2}}{\displaystyle \frac{(z-\mathrm{5.0})}{\mathit{\alpha }g}},T_{\mathrm{0}}\right\},$$

with α the thermal expansion coefficient in the equation of state defined as $\mathit{\rho }=-\mathit{\alpha }{\mathit{\rho }}_{\mathrm{0}}(T-T_{\mathrm{0}})$ with ρ₀=1024 kg m⁻³. A zonal wind is imposed with $u_{\star }=\mathrm{0.02}$ m s⁻¹, and the Stokes drift is given by

$${\mathit{u}}_s=(u_s,\mathrm{0}),u_s={\left({\displaystyle \frac{\mathrm{2}\mathit{\pi }a}{\mathit{\lambda }}}\right)}^{\mathrm{2}}\sqrt{{\displaystyle \frac{g\mathit{\lambda }}{\mathrm{2}\mathit{\pi }}}}{e}^{-\mathrm{4}\mathit{\pi }z/\mathit{\lambda }}.$$

The various parameter values are

$$\begin{array}{rl}& f_{\mathrm{cor}}={\mathrm{10}}^{-\mathrm{4}}{\mathrm{s}}^{-\mathrm{1}},h_{max}=\mathrm{120}\mathrm{m},T_{\mathrm{0}}=\mathrm{16}{}^{\circ }\mathrm{C},\\ & N_{\mathrm{0}}^{\mathrm{2}}={\mathrm{10}}^{-\mathrm{5}}{\mathrm{s}}^{-\mathrm{2}},\end{array}$$

with 96 vertical levels for the discretization and 16 h simulations. We only consider the case with a=1 m and λ=40 m, which gives a turbulent Langmuir number of La_t≈0.32. Numerical results are shown in Fig. B1 (upper panels) and are consistent with the results of Noh et al. (2016) with a deepening of the oceanic mixing length of about 10 m when Langmuir turbulence is accounted for (see LES results in Fig. 3 in Noh et al., 2016). For C_LC=0.15 in the Axell (2002) parameterization, the deepening is too weak, while for C_LC=0.3 it is closer to the Noh et al. (2016) LES results. Note that for those experiments, the value of d_LC is almost identical to the mixed layer depth. Figure B1 (lower panels) illustrates the fact that for a stronger stratification (i.e., with $N_{\mathrm{0}}^{\mathrm{2}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻² instead of $N_{\mathrm{0}}^{\mathrm{2}}={\mathrm{10}}^{-\mathrm{5}}$ s⁻²) the effect of Langmuir turbulence on mixed layer depth is negligible. Indeed, in this case Langmuir cells do not provide enough mixing to erode the stratification.

Figure B1 Solution obtained for the Noh et al. (2016) single-column experiment after 16 h for different parameter values in the Axell (2002) Langmuir cell parameterization in the case $N_{\mathrm{0}}^{\mathrm{2}}={\mathrm{10}}^{-\mathrm{5}}$ s⁻² (upper panels) and $N_{\mathrm{0}}^{\mathrm{2}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻² (lower panels).