{'fraction_non_alphanumeric': 0.036469165659008464, 'fraction_numerical': 0.026553808948004836, 'mean_word_length': 4.550469798657718, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Surface Wave Effects in the NEMO Ocean Model: Forced and Coupled Experiments

The NEMO general circulation ocean model is extended to incorporate three physical processes related to ocean surface waves, namely the surface stress (modified by growth and dissipation of the oceanic wave field), the turbulent kinetic energy flux from breaking waves, and the Stokes-Coriolis force. Experiments are done with NEMO in ocean-only (forced) mode and coupled to the ECMWF atmospheric and wave models. Ocean-only integrations are forced with fields from the ERA-Interim reanalysis. All three effects are noticeable in the extra-tropics, but the sea-state dependent turbulent kinetic energy flux yields by far the largest difference. This is partly because the control run has too vigorous deep mixing due to an empirical mixing term in NEMO. We investigate the relation between this ad hoc mixing and Langmuir turbulence and find that it is much more effective than the Langmuir parameterization used in NEMO. The biases in sea surface temperature as well as subsurface temperature are reduced, and the total ocean heat content exhibits a trend closer to that observed in a recent ocean reanalysis (ORAS4) when wave effects are included. Seasonal integrations of the coupled atmosphere-wave-ocean model consisting of NEMO, the wave model ECWAM and the atmospheric model of ECMWF similarly show that the sea surface temperature biases are greatly reduced when the mixing is controlled by the sea state and properly weighted by the thickness of the uppermost level of the ocean model. These wave-related physical processes were recently implemented in the operational coupled ensemble forecast system of ECMWF.


Introduction
Surface waves affect the ocean surface boundary layer (OSBL) through a number of processes, but perhaps most visibly through breaking waves which can be seen as whitecaps on the ocean surface (Monahan, 1971;Wu, 1979). These breaking waves enhance the turbulence in the upper part of the ocean significantly (Craig and Banner, 1994;Craig, 1996). Waves absorb energy and momentum from the wind field when they grow and in turn release it when they break Rascle et al., 2006;Ardhuin and Jenkins, 2006;Janssen, 2012). This lowers or raises the stress on the water side (i.e., the stress below the oceanic wave field) relative to the air-side stress, depending on whether the sea state is growing or decaying. Only when the wave field is in equilibrium with the energy injected by the wind will the stress on the two sides of the surface be equal.
Through the interaction with the Coriolis effect, the Stokes drift velocity associated with the wave field adds an additional term to the momentum equation. The effect was first presented by Hasselmann (1970) and has since been investigated for idealized cases by Weber (1983), Jenkins (1987), McWilliams and Restrepo (1999) and McWilliams and Sullivan (2000) among others. The force is variously known as the Stokes-Coriolis force or the Hasselmann force depending on whether it is considered to be purely an effect of the average Coriolis force acting on a particle with a Lagrangian velocity as given by the mean currents and the waves or as a tilting of the planetary vorticity (Polton et al., 2005;Broström et al., 2014). The force does not directly modify the total mass transport but it will alter the distribution of momentum over the depth of the Ekman layer (McWilliams and Restrepo, 1999;Polton, 2009).
The impact of the oceanic wave field on upper-ocean mixing and mean properties has been studied in a number of single-column mixed-layer model experiments (Craig and Banner, 1994;McWilliams and Restrepo, 1999;McWilliams and Sullivan, 2000;Burchard , 2001;Kantha and Clayson, 2004;Mellor and Blumberg, 2004;Rascle et al., 2006;Ardhuin and Jenkins, 2006;Huang et al., 2011;Janssen, 2012). Several studies have employed large eddy simulations (LES) to investigate the impact of Langmuir turbulence in the upper ocean (Skyllingstad and Denbo, 1995;McWilliams et al., 1997;Teixeira and Belcher, 2002;Polton and Belcher, 2007;Grant and Belcher, 2009), and in some cases even direct numerical simulations (DNS) have been employed (Sullivan et al., 2004). Most of these studies find that waves do indeed seem to have a rather profound impact on the upper part of the ocean, but there is still considerable disagreement about which processes are more important. So far there have been few studies of the wave impact on threedimensional ocean circulation models or fully coupled models of the ocean, the atmosphere and the oceanic wave field although the potential impact of waves on the climate system is recognized Cavaleri et al., 2012;Fan and Griffies, 2014). Fan et al. (2009) demonstrated the importance of correctly modelling momentum and energy fluxes from the wave field to the ocean under hurricane conditions. Fan and Griffies (2014) found that the introduction of Langmuir turbulence following the parameterizations by McWilliams and Sullivan (2000) and Smyth et al. (2002) as well as the parameterization of mixing by non-breaking waves suggested by Qiao et al. (2004) significantly changed the upper-ocean temperature in long-term coupled climate integrations. This latter mixing process appears similar to the mixing due to the high Reynolds numbers of the orbital motion of non-breaking waves explored by Babanin (2006) and Babanin and Haus (2009). Using a climate model of intermediate complexity,  explored three wave-related mixing processes, namely injection of turbulent kinetic energy from breaking waves, Langmuir circulation and the aforementioned mixing by non-breaking waves. Like Fan and Griffies (2014) they found that all three processes contributed to the mixed layer depth and the temperature of the mixed layer. Similarly, Huang et al. (2011) coupled WAVEWATCH III (Tolman et al., 2002) to a version of the Princeton Ocean Model (Blumberg and Mellor, 1987) and demonstrated an improved summertime temperature profile using the non-breaking parameterization of Qiao et al. (2004). They found very little direct impact by the breaking waves on the temperature.
These wave-driven processes influence the vertical structure of the temperature and current fields in the mixed layer in general, and in the upper few meters in particular. This has implications for coupled models as these processes will affect the feedback between the ocean and the atmosphere (Janssen et al., 2013). However, on shorter time scales and at higher spatial resolution it is also clear that these processes will influence the drift of objects and pollutants on the sea surface or partially or wholly submerged. This has practical importance for oil spill modelling (Hackett et al., 2006), search and rescue (Breivik and Allen, 2008;Davidson et al., 2009;Breivik et al., 2013) and dispersion of biological material (Röhrs et al., 2014). The NEMO ocean model (Madec and the NEMO team, 2012) has been coupled to the atmospheric model with wave forcing from the wave model as part of the ensemble suite of the Integrated Forecast System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF) since November 2013 (IFS Cycle 40R1). Here we describe the implementation of the three wave effects mentioned above in forced (ocean-only) integrations of NEMO using forcing from the ERA-Interim reanalysis (Dee et al., 2011) as well as their implementation in a fully coupled atmosphere-wave-ocean seasonal forecast system. The paper is organized as follows. Sec 2 describes the processes that have been implemented and lays out their actual implementation in NEMO. Sec 3 describes the results of long oceanonly integrations and compares with control runs, observations and the ORAS4 ocean reanalysis . Sec 4 describes the coupled atmosphere-wave-ocean coupling used for seasonal integrations and compares the results of a control run where no direct coupling exists between the wave model and the ocean model to a run where NEMO is forced with stresses, turbulent fluxes and Stokes drift from the wave model ECWAM (ECMWF , 2013). Sec 5 discusses the results and the deficiencies in the existing model setup. Sec 6 concludes and makes suggestions for further work on the investigation of wave effects in ocean-only as well as coupled atmosphere-wave-ocean models.

Wave effects in the Ocean Surface Boundary Layer
Introducing wave forcing in an Eulerian ocean model entails communicating the relevant forcing fields from a wave model. We start with a brief presentation of spectral wave models and how the two-dimensional wave spectrum relates to the fluxes and fields that have a bearing on the ocean surface boundary layer.

Fluxes and fields estimated from a spectral wave model
Third generation spectral wave models (Hasselmann et al., 1988;Tolman, 1991;Komen et al., 1994;Ris et al., 1999;Janssen, 2004;Tolman et al., 2002;Holthuijsen, 2007;Cavaleri et al., 2007) solve the action balance equation (see Eq (1.185) by Komen et al. (1994) and Eq (2.71) by Janssen (2004)) for the wave action density N (a function of the Cartesian co-ordinate x, frequency f and direction θ) as follows, Here, the right-hand source terms refer to wind input (in), nonlinear transfer (nl), and dissipation due to wave breaking (ds), respectively. The dissipation term may include shallow-water effects and bottom friction. The gradient in frequency represents shoaling and refraction, and ω = σ + k · u is the absolute frequency as seen by an observer standing still, whereas σ is the intrinsic frequency as seen by an observer moving with the current. We also note that where c g is the wave group veloctiy vector. The wave action density is related to the wave variance density through N = F/σ. In deep water with no current refraction Eq (1) reduces to the energy balance equation, written here in flux form. Note that the source terms in Eq (3) are related to those in the action balance equation (1) as S = σS ′ .

The air-side stress modified by surface waves
The presence of an undulating surface affects the roughness felt by the airflow. The atmospheric momentum flux to the oceanic wave field is denoted τ in . It is convenient to define an air-side friction velocity in relation to the total air-side stress, τ a , as Here, ρ a is the surface air density. Charnock (1955) was the first to relate the roughness of the sea surface to the friction velocity, where α CH is known as the Charnock constant. Janssen (1989Janssen ( , 1991 assumed that α CH is not constant but varies with the sea state, whereα CH = 0.006 (Bidlot, 2012) and the wave-induced stress, τ in , is related to the wind input to the wave field as Here ρ w is the water density. A diagnostic spectral tail proportional to ω −5 has been applied above a cutoff frequency ω w (Komen et al., 1994;ECMWF , 2013). The wave-modified drag coefficient is then where κ = 0.4 is von Kármán's constant. Note that the drag coefficient as defined here is related to the 10-m neutral wind speed, U 10N . This drag coefficient is computed by ECWAM. The Charnock parameter (6) is the main coupling mechanism between the atmosphere and the wave field in IFS, in place since 1998 (Janssen, 2004).

The water-side stress modified by surface waves
As the wind increases, the wave field responds by first growing and storing more momentum. In this phase there is a net influx of momentum to the wave field. Then, as the waves mature and the breaking intensifies, the momentum flux from the wave field to the ocean starts to close on the flux from the atmosphere to the waves. This is the equilibrium state where dissipation matches wind input, also referred to as fully developed windsea (since the waves cannot become higher), see e.g. Komen et al. (1994); World Meteorological Organization (1998); Holthuijsen (2007). Finally, as the wind dies down there will be a net outflux of momentum from the wave field, almost all of which will go to the ocean. If wind input and dissipation in the wave field were in equilibrium, the air-side stress would be equal to the total water-side stress. By water-side stress is meant the stress as seen by the Eulerian ocean, i.e., the momentum flux from the waves. However, most of the time waves are not in equilibrium (Janssen, 2012;Janssen et al., 2013), giving differences in air-side and water-side stress of the order of 5-10%, with occasional departures much larger in cases where the wind suddenly slackens. Likewise, in cases with sudden onset of strong winds the input from the wind field will be much larger than the dissipation to the ocean, lowering the waterside stress to values well below 70% of its normal ratio to the air-side stress. The water-side stress thus equals the total atmospheric stress minus the momentum flux absorbed by the wave field (positive) minus the momentum injected from breaking waves to the ocean (negative), Here, the dissipation source term is assumed to include all relevant dissipative processes. This can be written (ECMWF , 2013) The stress from waves is archived as a normalized quantity and is applied as a factor to the air-side stress in our implementation in NEMO.

Mixing parameterizations
The TKE equation with Reynolds averages can be written Here, e ≡ q 2 /2 = u ′ i u ′ i /2 is the TKE per unit mass (with q the turbulent velocity) and ǫ the dissipation rate (see e.g. Stull (1988), p 152). NEMO has the option of modelling the evolution of TKE with local closure (a prognostic equation in e only, see Stull 1988 pp 203-208 andPope 2000 pp 369-373). Assuming that the advective terms are small in comparison and making the gradient transport approximation where turbulent coefficients are proportional to the gradients in the mean quantities, we arrive at This is the standard one-equation formulation for NEMO (see the reference manual for NEMO v3.4, Madec and the NEMO team 2012, pp 176-177). Here l ǫ is the mixing length. The buoyancy term is assumed proportional to the local Brunt-Väisälä frequency, and the shear production is related to the shear of the mean flow, Finally, the mixing length is given by a relation by Blanke and Delecluse (1993), see also Gaspar et al. (1990) and Madec and the NEMO team (2012), pp 177-179. Two non-standard mixing processes present in NEMO's TKE scheme warrant our attention. The first is an artificial boost to the TKE known as the ETAU parameterization which is pegged to the surface TKE with an exponential vertical decay (see Madec and the NEMO team (2012), The depth scale h τ can vary with longitude from 0.5 m at the Equator to 30 m poleward of 44 • or be fixed at 10 m. The coefficient, here 0.05, can also be varied. The second mixing process of interest to us is a parameterization of Langmuir turbulence according to Axell (2002) which has been implemented in NEMO. The vertical velocity w LC of the Langmuir cells is assumed to peak at H LC /2, half the maximum depth to which Langmuir cells penetrate, Here v s is the surface Stokes drift speed and c LC = 0.15 is a coefficient. Axell (2002), by making an analogy with the characteristic convective velocity scale (D'Alessio et al., 1998) further assumed the Langmuir production term in the TKE equation (10) could be written This production term will attain a maximum value in the interior of the mixed layer. Craig and Banner (1994), CB94 hereafter, demonstrated that as waves break they will considerably modify the vertical dissipation profile from the traditional law-of-the-wall where dissipation ∝ z −1 (Stull , 1988). With wave breaking CB94 found dissipation ∝ z −3.4 . Terray et al. (1996) and Drennan et al. (1996) later demonstrated that the observed dissipation rates under breaking waves are indeed much higher than anticipated by the law of the wall. CB94's model has since been extended to a two-equation turbulence model by Burchard (2001), who demonstrated that the injection of turbulent kinetic energy from breaking waves was sufficient to successfully model the evolution of the mixed layer representative of North Sea conditions. CB94 suggested that the flux of turbulence kinetic energy (TKE) should be related to the water friction velocity w * as CB94 assumed that α CB was a constant ∼ 100, but noted that its range would probably be between 50 and 150, depending on the sea state (see also Mellor and Blumberg 2004). The TKE flux from breaking waves is related (see e.g. Janssen et al. 2004;Rascle et al. 2006;Janssen 2012;Janssen et al. 2013) to the dissipation source function of a spectral wave model as For consistency we have written the energy flux from the waves (thus always negative), m ≈ − ρ a /ρ w α CB , normalized by the air friction velocity u * . In NEMO the energy flux from breaking waves is introduced as a Dirichlet boundary condition on TKE, following Mellor and Blumberg (2004). It is assumed that in the wave-affected layer the mixing length can be set to a constant l w = κz w where the surface roughness length relates to the significant wave height H s as z w = 0.5H s , and that in this near-surface region diffusion balances dissipation. In this case the TKE equation takes a simple exponential solution [see Eq (10) by Mellor and Blumberg 2004], Here the inverse length scale is with S q = 0.2 and B = 16.6 given by Mellor and Yamada (1982). This is how the flux from breaking waves is implemented in NEMO v3.4. This allows the following simple boundary condition However, the inverse depth scale (20) is sea state dependent, and for a wave height of, say, 2.5 m, which is close to the global mean, λ −1 ≈ 0.5 m. Thus, e(z) varies rapidly with depth, and we have modified the boundary condition (21) by weighting the surface value by the thickness of the topmost level to attain an average value more representative of the turbulence near the surface of the model, Here L = ∆z 1 /2 is the depth of the T -point of the first level. Integrating Eq (19) is straightforward, and the average TKE boundary condition (22) becomes It is clear that as the vertical resolution increases the difference between e 1 and e 0 becomes smaller, and in the limit the two coincide. The weighting (23) is thus less important with higher vertical resolution. It is worth noting that the exponential profile (19) assumed by Mellor and Blumberg (2004) is only valid very near the surface, and in fact CB94 had already found the solution to the more general case where the mixing length is allowed to vary with depth. We have not implemented this operationally, but preliminary tests suggest that the effect is to roughly double the depth over which the TKE from breaking waves is distributed. The derivation is presented in the appendix.

The Stokes-Coriolis forcing
Waves set up a Lagrangian displacement v s in the down-wave direction known as the Stokes drift velocity (Stokes, 1847). Although it decays rapidly with depth, it can be substantial near the surface (|v s | ∼ 0.7 m s −1 ). In combination with the Earth's rotation it adds an additional veering to the upper-ocean currents known as the Stokes-Coriolis force (Hasselmann, 1970), Here f is the Coriolis frequency,ẑ is the upward unit vector, p is the pressure and τ is the stress. The full two-dimensional spectrum is in principle required to compute the Stokes drift velocity profile Janssen, 2012), This is computationally demanding and full two-dimensional wave spectra from a numerical wave model (see e.g. ECMWF 2013) may not always available. It is therefore customary to introduce a simplified, monochromatic Stokes drift profile (see e.g., Carniel et al. Here the subscript "e" distinguishes the approximate profile from the full Stokes drift velocity profile (25). The surface Stokes drift velocity vector v 0 is computed by ECWAM and is available both in ERA-Interim (Dee et al., 2011) and from the operational ECMWF forecasts (ECMWF , 2013).
To compute the profile (26) we must find the inverse depth scale k e . This is related to the transport T s through the exponential integral E 1 (Abramowitz and Stegun, 1972) and can be solved analytically (Breivik et al., 2014) to yield Rearranging we get the following expression for the inverse depth scale, Here The n-th order spectral moment is defined as The mean frequency is defined as f = m 1 /m 0 (World Meteorological Organization, 1998; Holthuijsen, 2007) and the significant wave height H s = 4 √ m 0 . We can derive the first moment from the integrated parameters of a wave model or from wave observations and find an estimate for the Stokes transport, Herek s = (sin θ s , cos θ s ) is the unit vector in the direction θ s of the Stokes transport. We approximate the Stokes transport direction by the surface Stokes drift (see Breivik et al. (2014)).

Ocean-only Forced Model Experiments
The NEMO model is run on a tripolar ORCA 1 • grid configuration with 42 vertical levels. The uppermost level is 10 m thick. The model is coupled to LIM2, a two-level thermodynamicdynamic sea ice model (Fichefet and Maqueda, 1997;Bouillon et al., 2009) and is relaxed weakly towards a climatology in temperature (3-yr e-folding time). No sea surface temperature (SST) relaxation is performed. The ORCA grid is such that the resolution is increased towards the Equator (roughly 1/3 • ) to better resolve tropical waves (see Madec and the NEMO team (2012) for details on the ORCA grid). The atmospheric and wave forcing fields have been computed from the ERA-Interim reanalysis (Simmons et al., 2007;Dee et al., 2011). ERA-Interim is a continuously updated atmospheric and wave field reanalysis starting in 1979. The resolution of the wave model is 1.0 • on the Equator but the resolution is kept approximately constant globally through the use of a quasi-regular latitude-longitude grid where grid points are progressively removed toward the poles (Janssen, 2004). Similarly, the atmospheric model fields are archived on a reduced Gaussian grid of approximately 0.75 • resolution at the Equator. Some care has to be taken when interpolating between these grids, in particular where wave parameters are interpolated from the ECWAM grid to the ORCA grid. NEMO requires fluxes to be defined in all ocean points. However, there are discrepancies between the ice coverage and the land-sea mask of the wave grid and the ocean grid. This is solved by reverting to the ECMWF drag law (Janssen, 2008;Edson et al., 2013) where ECWAM has ice or land, C D (z = 10 m) = (a + bU p 1 10 ) /U p 2 10 .
The coefficients are a = 1.03 × 10 −3 , b = 0.04 × 10 −3 , p 1 = 1.48 and p 2 = 0.21. Here U 10 is the 10-m wind speed from ERA-Interim. Where ECWAM and NEMO agree on open water, the stress is computed from the drag coefficient of ECWAM, Here, U 10N is the neutral 10-m wind speed, available on the ECWAM grid. The conversion to water-side stress is implemented as whereτ is the ratio of water-side to air-side stress (see also Eq 34). It is this parameter which is archived by ERA-Interim. A standard integration period covering the ERA-Interim period from 1979 up until the end of 2009 has been used in the following. A summary of the settings for the model runs can be found in Tables 1 and 2. Four experiments with the new wave-related effects are presented, all compared against a control experiment, CTRL, where standard settings are used for NEMO. The CTRL experiment includes a parameterization of the TKE flux from breaking waves (CB94) without explicit sea state information (see Eq (17) and also Madec and the NEMO team (2012), Sec 10.1) but has no averaging over the topmost model level (22). The stress in CTRL is computed using the ECMWF drag law (32) with air-side stress (33). In all runs the NEMOVAR observation operator y = H(x) relating model state x to observation space (y) is applied. This allows a comparison of model integrations (in our case without assimilation) against the large number of quality-controlled temperature and temperature-salinity profiles compiled in the EN3 data set (Ingleby and Huddleston, 2007). Seasonal averages (December to January, DJF, and June to August, JJA) over the period 1989-2008 are used to compare SST fields in the following.
The first wave experiment, TAUOC, uses the water-side stress described in Sec 2.3 together with the ECWAM drag coefficient (see Sec 2.2). The effect is confined mostly to areas with rapidly developing weather systems in the extra-tropics (Fig 1), where the sea state will be quite far from equilibrium. There is a slight weakening of the wind stress along the west coast of South America and along the coast of south-west Africa, leading to decreased upwelling. This is mainly a consequence of differences between the ECWAM drag coefficient and the drag law (32) caused by limited fetch near the coast. The overall effect is a slight reduction of the bias in the tropics compared to EN3 near-surface temperature measurements (0.1 K, not shown). In the tropics, differences are most likely due to differences in the drag over swell compared to the drag law (32) used for the CTRL experiment.
The second experiment, WTKE, introduces the TKE flux (18) from ECWAM (described in Sec 2.4). The differences found in the extra-tropics amount to more than 2 K (Fig 2). The large difference suggests that the standard settings of NEMO overdoes the mixing due to waves, especially with coarse vertical resolution (cf Eq 22).
The STCOR experiment introduces Stokes-Coriolis forcing from ECWAM as described in Sec 2.5. The largest impact (Fig 3) is found in areas with extra-tropical cyclones in combination with strong temperature gradients, such as across the Gulf Stream and the Kuro-Shio.
Another experiment called LOW where a Law-of-the-wall boundary condition is applied, i.e. no TKE flux from breaking waves, was also performed. The motivation is to establish a lower bound on the mixing. The difference between the WTKE and LOW runs is much smaller (not shown) than the difference between CTRL and WTKE.
We now combine the three wave experiments TAUOC, WTKE and STCOR in one experiment referred to as the WAVE run. Fig 4 shows the standard deviation (upper curves) and the biases of the CTRL run (blue) (blue) (blue) (blue) (blue) (blue) (blue) (blue) and a wave run including all three wave effects (green). The most striking feature is the large-amplitude annual cycle observed for the bias of the CTRL run. This amplitude is much weaker for the run with wave effects, and the seasonality is not at all so clear. Due to the huge differences to the mixing in the runs with and without wave effects, the heat uptake of the ocean also differs significantly. Fig 5 shows the global total heat content anomaly relative to 1979 (the start of the period) for the experiment with wave effects (green) v the CTRL experiment (blue). We also show the ocean reanalysis ORAS4 (red) which is based on an earlier version of NEMO (v3.0), see Balmaseda et al. (2013), with observations assimilated using the NEMOVAR 3D-VAR assimilation system (Mogensen et al., 2012a). The reanalysis covers the period 1957 to present, with forcing provided by ERA-40 (Uppala et al., 2005) and later ERA-Interim from 1979 and onwards. The similarity with ORAS4 is clear, and the trends almost identical. It is also clear that the impact of the wave mixing establishes itself within the first two years of the integration. The pronounced annual cycle (dashed lines) is due to the difference in oceanic volume in the southern and northern hemispheres.
To further assess the impact on the surface temperature we now compare with OI v2, an SST analysis by Reynolds et al. (2002). The comparison with the CTRL run in Fig 6a reveals large biases, especially in the summer hemisphere. These biases are reduced when wave effects are included (Fig 6b), especially in the northern extra-tropics. This is mainly due to a more correct level of mixing. The long-term SST fields from ORAS4 are very similar to OI v2 since a strong relaxation to the OI v2 gridded SST was applied as a flux correction between December 1981 and December 2009  and are not shown here.
Finally, to test the relative impact of the ETAU TKE boost (14) and the Langmuir turbulence parameterization (16), we switched off the Langmuir circulation (NOLC) and ETAU (NOETAU) in two additional experiments. As is evident from Fig 7, the Langmuir turbulence has only modest impact on the temperature in the mixed layer in the extra-tropics (50 m depth in the northern extra-tropics is shown in Fig 7). Switching off ETAU has a much larger impact on the mixing, both in terms of bias and random error.

The Coupled Atmosphere-Wave-Ocean Model
Coupling between the wave model and the ocean model was first implemented operationally in the ensemble suite of the IFS in Cycle 40R1 in November 2013, but was limited to the mixing (WTKE) and Stokes-Coriolis forcing (STCOR). The seasonal integrations described here include in addition the modified stress (TAUOC) and the ice model LIM2 but are run at lower spatial resolution than the operational ensemble forecasts. Fields are exchanged every three hours (the coupling timestep) between IFS/ECWAM and NEMO. The atmospheric model is run with a spectral truncation of T255 (corresponding to roughly 78 km) with 91 vertical levels to 1 Pa. ECWAM is run at 1.5 • resolution while NEMO is run on the ORCA 1 • grid described in Sec 3. Fig 8 shows a flow chart outlining the sequence of execution of the various components of the IFS-ECWAM-NEMO single executable over two coupling timesteps, which can be summarized as follows.
1. The atmospheric component (IFS) is integrated (internal time step 2,700 s) one coupling time step (10,800 s), yielding wind fields for ECWAM as well as radiation and evaporation minus precipitation fields for NEMO 2. ECWAM (internal time step 900 s) is tightly coupled to IFS and returns sea surface roughness to the atmospheric boundary layer at every atmospheric time step. After one NEMO coupling time step, ECWAM also gives Stokes drift velocity, turbulent kinetic energy and waterside stress to NEMO 3. NEMO is integrated (internal time step 3,600 s) one coupling time step, yielding SST and surface currents to IFS 4. All model components have now been integrated one coupling time step and the sequence begins anew Surface currents and SST are communicated back to the atmospheric model and will affect the stress and the temperature of the boundary layer. This in turn affects the oceanic wave field, but there is presently no direct feedback from NEMO to ECWAM. The coupling between the different components is described in more detail by Mogensen et al. (2012b). Here we compare a setup for seasonal integrations to seven months, with three ensemble members starting from 1 May and 1 November for the period 1993-2012. The CTRL experiment is run with a standard CB94 type wave mixing which corresponds to α CB = 100 (Eq 21), similar to the ocean-only CTRL experiment presented in Sec 3. The wind stress is computed using the ECMWF drag law (32) and 10-m wind vectors from IFS. The wave experiment includes the three processes TAUOC, WTKE and STCOR described in Sec 2 and Sec 3. Fig 9 reveals large differences in the bias in the northern extra-tropics relative to ERA-Interim for the boreal summer (JJA) at a lead time of one to three months from 1 May. Panel (a) shows the bias of the CTRL run, with large cold biases in the northern extra-tropics. These biases are broadly similar to what was found from the ocean-only (forced) runs presented in Sec 3 (cf Fig 2). Panel (b) reveals the biases in the run with wave effects to be generally much smaller, although there is a certain deterioration of the upwelling area along the coast of Baja California. Note that the cold bias in the eastern equatorial Pacific (the "cold tongue") is reduced slightly. The results are similar for the southern hemisphere summer (DJF), although the bias reduction is not as strong as for the northern hemisphere. The seasonal variation of the bias found for the forced (ocean only) runs in Fig 4 is also present in the coupled integrations, see Fig 10. Again, the bias is greatly reduced in the wave run.

Discussion
We have introduced three wave effects in NEMO, namely the sea-state dependent water-side stress, the energy flux from breaking waves and the Stokes-Coriolis force. Using ocean-only integrations and experiments with a coupled system consisting of the atmospheric model IFS, the wave model ECWAM and NEMO, we demonstrated that the impact of the wave effects is particularly noticeable in the extra-tropics. Of the three processes, the modification of the mixing (WTKE) has the largest impact (Fig 2), but as we discuss below, this is also related to the additional mixing found in NEMO. The impact of the modified stress (TAUOC) and Stokes-Coriolis (STCOR) is also significant (on the order of 0.5 K locally, see Figs 1 and 3). It is also important to note that compared to the law-of-the-wall experiment (LOW), the WTKE differs by only 0.5 K. This again suggests that the CTRL experiment has too vigorous mixing. In ocean-only integrations we see a reduction of the temperature bias in the mixed layer, particularly in the extra-tropical summer (Fig 4). This manifests itself in a more realistic oceanic heat uptake (Fig 5). The coupled seasonal integrations show a similar reduction in bias (compared to ERA-Interim, Fig 9) as the ocean-only wave run (see Fig 10). The mixing is strongly influenced by the ETAU parameterization (14). It is clear that the TKE scheme (11) has too shallow mixing without this parameterization, and it is also clear that the present parameterization of Langmuir turbulence (16) does very little due to its vertical structure which has the distinguishing feature that it will put most if not all of the enhanced turbulence deep into the mixed layer and nothing near the surface. It is thus unable to transport down enough heat to make a substantial difference. This explains why its impact on the temperature in the OSBL (Fig 7) is so much smaller than that from the ETAU term. Its vertical profile (15) is very different from Langmuir parameterizations involving the shear of the Stokes drift velocity profile (McWilliams et al., 1997;Polton and Belcher, 2007;Grant and Belcher, 2009), Here, u ′ H is the horizontal velocity fluctuations and w ′ the vertical. Due to the strong shear of the Stokes drift profile this would add a larger contribution near the surface. The ETAU profile (14) is quite similar to (35) and it appears to act as a parameterization for Langmuir turbulence with its characteristic exponential decay with depth (35), or similarly mixing by non-breaking waves (Qiao et al., 2004;Babanin, 2006;Huang et al., 2011). It facilitates deeper penetration of mixing from surface processes than what is normally assumed from breaking waves (Grant and Belcher, 2009;Belcher et al., 2012).

Concluding remarks and further work
The ocean-only integrations and coupled seasonal integrations all suggest that the right level of mixing is very important for reducing the temperature bias in the upper part of the ocean and also for the oceanic heat uptake. An important result is that introducing wave-enhanced mixing must be done in such a way that the thickness of the uppermost layer is accounted for. This is done with the present implementation by weighting with the thickness of the layer and is essential with model configurations with a thick uppermost level, e.g. ORCA1L42 as discussed here. This would not be necessary if a flux boundary condition had been used for the TKE from breaking waves, but this is not the case in NEMO, which uses the surface boundary condition first proposed by Mellor and Blumberg (2004). The impact of mixing on SST is clearly shown in Fig 6 where runs with and without wave effects are compared with the OI v2 SST analysis. That the seasonal cycle of the CTRL run is distorted by too vigorous mixing is clear, and Fig 4 shows that the annual cycle in biases extends well below the surface. The conclusions from the ocean-only experiments that temperature biases are reduced by introduction of wave-induced mixing are borne out by the seasonal coupled integrations which essentially show the same bias reduction (Fig 9). It is important to note here, though, that the additional ad hoc deep mixing in NEMO interacts with the surface processes and that without this additional mixing the model fails to mix deeply enough. We speculate that this mechanism is really masking Langmuir turbulence or mixing from non-breaking waves. More work is clearly needed with ocean circulation models and coupled models to fully answer the question of which mixing processes are dominant in the OSBL, but it is clear that getting the mixing right is a balancing act between the right deep mixing and the right mixing near the surface, and these processes are probably all wave-related.
These results are relevant for assessing the impact surface waves have on climate projections Fan and Griffies, 2014), and a natural next step would be to investigate the impact of waves on long, decadal to century-wide integrations (see also the Co-ordinated Ocean-Wave Climate Projections (COWCLIP) initiative, Hemer et al. (2012)). One candidate for forcing ocean-only integrations would be the recently completed ERA 20th century reanalysis (ERA-20C, see Poli et al. 2013, de Boisséson et al. 2014, and Dee et al. 2014. This opens up the possibility of running century-long NEMO integrations with wave effects from a state-of-the-art version of ECWAM (Bidlot, 2012) since all relevant parameters have been archived in the new reanalysis. For coupled climate projections, the nearest candidate would be EC-Earth (Hazeleger et al., 2010(Hazeleger et al., , 2012 which operates a modified version of an earlier cycle of IFS. Such experiments would help determining the importance of waves in the climate system, rather than just the impact of climate change on the wave climate.

A The dissipation profile
The exponential profile (19) for the balance of assumed by Mellor and Blumberg (2004) is only valid very near the surface where the mixing length can be assumed constant = κz w . CB94 presented the solution to the more general case where the mixing length is allowed to vary with depth. This equation has a power-law solution [cf CB94, Eq (23)], This leads to a slightly more complicated expression for the vertical average (22), than Eq (23). The consequence is that the mixing penetrates about twice as deep as in the case where the exponential approximation assumed by Mellor and Blumberg (2004).   Figure 1: Long-term SST differences between a run with water-side stress modulated by the ECWAM wave model (TAUOC experiment, see Table 2) and the CTRL run.
Figure 2: Long-term SST differences between a run with TKE flux from ECWAM (WTKE experiment, see Table 2) and the CTRL run. Note that the color scale is ±2 K. Figure 3: Long-term SST differences between a run with Stokes-Coriolis forcing (STCOR experiment, see Table 2) and the CTRL run.  (Ingleby and Huddleston, 2007) in the northern extra-tropics in the CTRL run (blue) and a run with all three wave effects switched on (green), see Table 2 for more information about the experiments. The upper curves show the standard deviation while the lower curves represent the bias. A 90-day running mean is employed.  Table 2). Full lines represent a 12-month moving average. The ORAS4 reanalysis is shown in red and the CTRL run in blue. Panel b: Total ocean heat content anomaly.   Figure 8: A flow chart of the coupling between the components of the single executable IFS-ECWAM-NEMO model setup. Two coupling time steps are shown to illustrate the sequence. This is how the operational ensemble prediction system has been run since Cycle 40R1 (November 2013). The seasonal integration experiments described here contain in addition the LIM2 ice model and are run at lower atmospheric resolution (T255). First the atmospheric component (IFS) is integrated (internal time step 2,700 s) one coupling time step (10,800 s), yielding wind fields for ECWAM as well as radiation and evaporation minus precipitation fields for NEMO. ECWAM (internal time step 900 s) is tightly coupled to IFS and gives sea surface roughness for the atmospheric boundary layer at every atmospheric time step. After one NEMO coupling time step, ECWAM also gives Stokes drift velocity, turbulent kinetic energy and waterside stress to NEMO. NEMO is integrated (internal time step 3,600 s) one coupling time step, yielding SST and surface currents to IFS. All model components have now been integrated one coupling time step (10,800 s), and the sequence begins anew with the next coupling time step.
(a) (b) Figure 9: Panel a: SST bias relative to ERA-Interim for the coupled seasonal CTRL run (JJA). Panel b: Same for the wave run. Figure 10: A comparison of surface EN3 temperature observations (Ingleby and Huddleston, 2007) in the northern extra-tropics for coupled integrations of IFS-ECWAM-NEMO as a function of lead time. The start date is 1 May. A three-member ensemble is run to a lead time of seven months for the years 1993-2008. The dashed red (lower) curve shows the bias of the run with wave effects, and the green dashed curve shows the bias of the CTRL run. The upper full lines show the root-mean-square (RMS) error.