The influence of surface gravity waves on the injection of turbulence in the upper ocean

Observations were made in the near-surface layer, at about 8 m depth in 132 m deep water off the coast of Ålesund in Norway, for a duration of 2.5 months in late 2011. The measurement period covers the passage of two low pressure systems with substantial wind and wave forcing. The time series of the dissipation rate of turbulent kinetic energy,ε, and the estimates of surface gravity waves are analysed. Dissipation rates varied by 5 orders of magnitude and reached 10 −5–10−4 W kg−1 in conditions when wind speed exceeded 15 m s −1 and the significant wave height was of the order of 10 m. The data set suggests substantial injection of turbulence from breaking surface gravity waves and Langmuir turbulence. To support and interpret the observations, numerical calculations are conducted using a second-order turbulence closure scheme based on the Mellor–Yamada level 2.5 scheme, modified to incorporate the near-surface processes such as Langmuir circulation and wave breaking. The results from a run forced by observed wind and wave fields compare favourably with the observations. Comparisons with other near-surface data sets available from the literature lend confidence on our dissipation measurements and the wave-forced simulations.


Introduction
Turbulence plays a central role in air-sea interaction processes in a wide variety of temporal and spatial scales, including the exchange of heat, energy, momentum, and trace gases.It is also one of the controlling mechanisms for biological productivity such as primary production by phytoplankton, and photochemical production such as greenhouse gases.Surface gravity waves are ubiquitous phenomena spanning scales ranging from small bubbles to long wavelength swell.In certain conditions, surface waves break intermittently and influence the upper ocean mixing up to a depth of the order of one significant wave height.Furthermore, the wave-current interaction associated with large-scale coherent Langmuir cells transports momentum and energy deeper into the ocean column.The important impacts of wave dynamics on different air-sea disciplines have motivated the development of high-quality measurement systems and more sophisticated simulation techniques.
Several problems arise when measuring turbulence near the sea surface, especially in harsh and wavy conditions.Surface gravity waves exert strong disturbances on the microstructure measurements in the upper ocean.The wave orbital velocities influence measurements in both linear ways that can be removed from the measured signals using statistical filtering, and by non-linear frequency modulation that cannot be easily removed.Furthermore, bubble clouds generated from the surface wave breaking contaminate the nearsurface turbulence measurements.The lack of stable platforms makes near-surface measurements more challenging, and demands specialized sensor technologies, careful data analysis, and complex signal-processing techniques.
In order to give a realistic representation of the upper ocean mixing variability, the wave-induced mixing should be accounted for properly in the mixed layer models.There are three general sources for wave-induced mixing: wave breaking, Langmuir circulations (LC), and non-breaking waveturbulence interaction.Craig and Banner (1994) developed a two-equation closure model to account for the transfer of energy flux from waves to turbulence.This is carried out by introducing a wave-modified surface boundary condition in the ocean mixing models using, for example Mellor-Yamada (Mellor and Yamada, 1982), or k-ε models.Sullivan et al. (2007) modelled depth injection of momentum and energy by breaking waves using a large eddy simulation.Rascle et al. (2007) and Bakhoday-Paskyabi et al. (2012), using observations and model-data comparisons, showed that the observed profiles of velocity and the dissipation rate of turbulence kinetic energy (TKE), ε, are better captured if the wave forcing is taken into account in the mixing models.Kantha and Clayson (2004) revised the turbulence closure model to account for the LC effects by adding the Stokes drift production term to the TKE equations.These large, coherent structures were found to substantially modify the upper ocean mixing, dissipation, and Ekman currents.Following Kantha and Clayson (2004), Harcourt (2013) suggested a modified stability function, used in the second moment closure, from a second moment algebraic closure of the Reynolds stresses and flux equations.
Several theoretical, experimental, and observational studies have suggested that non-breaking waves can also transfer energy to the water column and influence upper ocean mixing (Kitaigorodskii et al., 1983;Babanin and Haus, 2009;Ardhuin and Jenkins, 2006;Anis and Moum, 1995;Gemmrich and Farmer, 2004;Qiao et al., 2004).Qiao et al. (2004) and Dai et al. (2010) proposed a parametrization for non-breaking wave-induced mixing.They pointed out that the wave-turbulence interaction efficiently improves the prediction of upper ocean mixing.Bakhoday-Paskyabi and Fer (2013a) applied a non-breaking wave-induced turbulence parametrization to the General Ocean Turbulence Model (GOTM) and compared simulation results with observations of ε.They pointed out that the non-breaking wave parametrization improves the prediction of ε in the upper ocean.However, the wave dynamics and their coupling with the upper ocean processes merit further investigations.
In this paper, we investigate upper ocean turbulence both experimentally and theoretically.To measure microstructure near the sea surface, we use a moored sub-surface platform equipped with different oceanographic sensors including shear probes, an acoustic Doppler velocimeter (ADV), and a high-resolution pressure sensor.Measured ε from shear probes are compared with a wave-modified mixed layer model using the level 2.5 Mellor-Yamada scheme, which contains a prognostic equation for TKE, including production and reduction of TKE.The study period is about 2.5 months, covering two strong storms.The main focus of the present modelling investigation is to examine two particular wave processes: transient wave breaking, and Stokes drift and/or Langmuir turbulence which are hypothesized to be the dominant sources of TKE during our experiment.
The field measurements and experiment site are briefly described in Sect. 2. In Sect. 3 the model formulations and a stochastic breaker model are presented.Subsequently, in Sect. 4 observations of waves and dissipation rate of TKE measurements are analysed, and model-observation compar- isons are given, especially during the storm-induced mixing events.Finally, Sect. 5 provides a summary and discussion.

Instrument
The instrument, Moored Autonomous Turbulence System (MATS) is an ocean turbulence measurement system designed to collect microstructure time series at a fixed level.The details of the instrument are described in Fer and Bakhoday-Paskyabi (2014), and only the salient features are repeated here.MATS consists of a main body platform equipped with a modified turbulence package microRider (Rockland Scientific International), a Nortek Vector ADV, and a pair of battery packs.The platform is a low-drag buoy with a main-body diameter of 46 cm, tapered to 15 cm diameter at the nose section with the turbulence sensors.In our deployment, the buoy is the upper buoyancy element of the mooring line, and a swivel allows the instrument to align with the current, pointing the sensors onto the undisturbed, free flow (Fig. 1).
The microRider is equipped with two air-foil shear probes, two fast-response FP07 thermistors, a pressure transducer, a 2-axis vibration sensor, an inclinometer (pitch and roll angles accurate to 0.1 • ), a 6-axis motion sensor, and a compass.The vector is a 6 MHz acoustic velocimeter measuring the 3-D velocity fluctuations in water.All turbulence sensors of the microRider and the sensor head of the vector protrude horizontally from the nose of the buoy pointing to the mean flow.No probe guard is installed.The sensor head of the vector is rigidly fixed to the buoy.The tip of the turbulence sensors is about 25 cm from the nose of the buoy.The shear probes are mounted orthogonally to each other to measure the ∂w/∂x and ∂υ/∂x shear components.A right-handed Cartesian coordinate system is used throughout with x pointing forward along the major axis of the instrument, y pointing to the port side of the instrument and z upward.Sampling rate is set to 512 Hz on all turbulence channels (vibration, shear and temperature gradient), 16 Hz for vector (3-D velocity) and 64 Hz for the other channels, including the high-resolution pressure sensor, compass and the motion pack.

Site and deployment
MATS was deployed off the coast of Ålesund in Norway, on 25 October 2011, 09:50 UTC at 62 • 50.28 N, 6 • 8.67 E (Fig. 2).The water depth at the position was 132 m.The instrument was recovered on 20 January 2012.MATS was set to sample 15 min bursts every hour.In total 1842 segments, of 15 min duration each, were recovered between 25 October 2011, 10:35 UTC and 10 January 2012, 03:50 UTC.
The measurement period covers two storm periods with wind speeds in excess of about 20 and 30 m s −1 (Berit and Dagmar, respectively) measured at the nearby Vigra meteorological station on 25 November and 25 December 2011 (see Sect. 4 for a description of the environmental conditions).The site is observed to be energetic.The significant wave height, H s , typical of the region, varies between 1 and 5 m, which increases to more than 12 m during the storm periods.While the hourly maximum velocity in the water column typically varies between 0.2 and 1 m s −1 , it reaches about 1.5 m s −1 during the storms (not shown).When the tidal variability is removed, the depth-averaged currents vary between 0.1 and 0.5 m s −1 , occasionally reaching values above 0.6 m s −1 .

Data processing and quality control
A detailed description of the data processing, and a discussion regarding the quality of data obtained from MATS are given in Fer and Bakhoday-Paskyabi (2014).The shear probe data voltage output is converted to shear using the known electronic constants, the sensitivity of the shear probe and the flow past the sensors measured by the vector.A smooth, 2 s low-pass filtered 3-D velocity field is used to ensure that the highly variable flow near the surface (due to wave orbital velocities) is accounted for as the flow advects past the sensors.
The dissipation rate is calculated by assuming isotropic turbulence and integrating the measured shear wavenumber spectrum in half-overlapping 60 s segments.Only the portion of the shear spectra between 1 and 20 Hz, relatively unaffected by wave motion, is used (example spectra are presented in Sect.4.2).The integration limits are described in Fer and Bakhoday-Paskyabi (2014), and the Nasmyth spectrum is used to account for the variance outside the integration band.In total, 51 520 estimates of ε are obtained.

Surface gravity wave estimate
Surface gravity waves are measured using the highresolution pressure sensor of MATS, located at approximately 8 m from the mean sea level.Using linear smallamplitude wave theory, in an irrotational wave field, ( ũ, w) = (∂φ/∂x, ∂φ/∂z), the sea surface elevation is where φ is the velocity potential, ũ and w are horizontal and vertical components of wave orbital velocities, respectively, g is the gravitational acceleration, and t is time.This surface elevation imposes a subsurface dynamic pressure p d that decays with depth as where ρ w is the water-side density.For a linear monochromatic wave, the elevation and velocity potential take the following forms: where k = 2π/λ is the wavenumber, λ is the wavelength, a is the wave amplitude, D is the mean water depth, and ω is the angular frequency obtained from the dispersion relation By ignoring the non-linear contributions in Eq. ( 2), the linearized dynamic pressure is calculated as in which K p is the pressure transfer function deduced from the linear-wave theory.Short waves are thus attenuated with depth more strongly than longer waves.The application of the transfer function leads to large corrections at high frequencies where the data are typically noisy.To avoid substantial noise growth, all frequencies beyond the peak frequency of the spectrum are swept to find a point (the cut-off frequency) at which the pressure signal first exceeds the signal noise floor by a factor of 10.The wave energy spectrum is extended by applying a f −5 shape beyond this cut-off frequency (Sect.3.1), where f is a frequency in Hz.The uncertainty in estimation of significant wave height, H s , from subsurface pressure data is shown in Fig. (3a).The uncertainty arises as a consequence of applying the linear-wave theory (Eq.6), resulting in increased relative errors for weak winds (Appendix A).
The contaminations induced by MATS pitching are removed using the vertical component of the accelerometer signal, a z .For small roll angles where ζ is the vertical displacement of MATS, θ is the platform pitch angle, and the umlaut denotes second derivative with respect to time.Hence, in the platform coordinate system, the dynamic pressure can be corrected by using The spectral distribution of p obs d , θ, and a z are shown in Fig. 11a) for an arbitrary 15 min segment.In practice, problems often arise in a calculation of ς due to an inherent low frequency drift of accelerometer sensors.This drift should be removed by using an appropriate high-pass filter (here at a period of 25 s) before integrating the accelerometer signal.After all post-processing operations, the relationship between frequency spectra of surface wave elevation and the corrected subsurface pressure is obtained as where, ς = ς/(ρ w g), p = p obs d /(ρ w g), S pς is the crossspectra between pressure data and MATS vertical displacement, and S pp and S ς ς are the auto-spectra of the observed dynamic pressure and MATS vertical displacement, respectively.We further utilize the corrected power spectrum to approximate wave bulk parameters by the means of spectral moments: Generally accepted approximations for significant wave height and the mean wave period estimations are based on the few first moments (the zeroth and second moments) of the power spectra 3 Wave-enhanced turbulence mixing In this section, we review the equations governing the wavemodified mean flow and TKE in a horizontally homogeneous ocean.For the TKE equation, we employ an improved 2.5 closure model to predict near-surface turbulence properties.In Sect.3.5, we present various parameterizations of ε used near the sea surface boundary layer.

The wave model
The wave energy evolution can be described in the form of an energy balance equation for a smooth and continuous wave energy spectrum, E(f ), in x-direction as in in which wave energy propagates at group velocity c g and changes through the contributions of spectral source and sink terms: wind energy input source term, S in , dissipation via breaking, S ds , and non-linear wave-wave interactions, S nl .
Based on the wave energy spectrum, the variance spectrum can also be specified as the energy spectrum divided by ρ w g.The total wave momentum, M(f ), is related to the wave energy by the phase speed, c p , as E(f ) = c p M(f ).Then, the wave-induced stress and dissipation stresses are calculated, using S in and S ds , as Correspondingly, the energy flux from wind to waves, in , and from waves to ocean, ds , are defined by where k denotes the wave direction.In the operational wave models, the prognostic frequency range is restricted due to some practical issues such as the computational cost.A successful candidate of spectral characteristics beyond the prognostic frequency range is the Philips spectrum where α p = 0.0081 is the Philips constant, f p is the spectral peak frequency, U 10 is the wind speed vector at 10 m height, |.| denotes the modulus of a vector, and the angular frequency ω is introduced via Eq.( 5).Equation ( 15), together with observational supports, suggests that the wave energy spectrum can be extrapolated in the diagnostic range, f > f max , by , where f max is the maximum modelled prognostic or observed frequency (Jones and Toba, 2008).

The equations of motion
The governing equations of wave-modified momentum and passive tracers are (Appendix B): where f cor is the Coriolis parameter, u and v are horizontal components of mean flow, u s and v s are the corresponding components of Stokes drift, and u , v , and w are the horizontal and vertical fluctuating velocity components, respectively.T and S are the mean temperature and salinity, and T and S are the fluctuating temperature and salinity, respectively; R is the solar radiation profile; the co-variances in Eqs. ( 16)-( 19) are turbulent fluxes that must be specified in order to close the system.Here, we use the following closures: where K m = lqS m and K h,s = lqS h,s are the eddy viscosity and diffusivity, respectively; q is the turbulent velocity, S m and S h = S s are dimensionless stability functions which are dependent on stratification as derived by Mellor and Yamada (1982) and modified by Galperin et al. (1988), and l is the mixing length in the upper ocean that can be approximated as where κ = 0.41 is the von Kármán constant and z 0 is the water-side roughness length (Sect.3.4).The Stokes drift profile in a deep, rotating ocean is where wave energy frequency spectrum governed by Eq. ( 10), and k = k k denotes the horizontal wavenumber.The terms −f cor u s and f cor v s are called Coriolis-Stokes forcing which controls the mean current profile with the potential to increase the probability of occurrence of Langmuir circulations as a function of latitude.The source terms, F x and F y , in Eqs. ( 16) and ( 17) are the horizontal components of ensemble-averaged momentum due to wave breaking.They are directly balanced by an extra term in the surface momentum boundary condition (Eq.26) to ignore the direct transport of momentum flux from air to currents.Jenkins (1987) introduced this wave-induced momentum transfer term from waves to the ocean due to dissipation of wave energy as (Appendix B) Here S ds is the wave dissipation source term in Eq. ( 10), and N controls the vertical distribution of F ds = (F x , F y ) over the water column such that In the absence of wave field information, Eq. ( 24) can be parameterized as where z is depth from the sea surface, and G m (z) denotes the vertical structure of the wave-induced momentum.Kudryavtsev et al. (2008) used a rectangular function for each monochromatic breaker at a depth proportional to its wavelength.Sullivan et al. (2004) proposed an analytical breaker momentum impulse based on measurements made by Melville and Matusov (2002).Jenkins (1987) used the depth dependence of a wave-induced momentum redistribution term in the water column by exp (−2k|z|) where k is the modulus of a horizontal wavenumber.Rascle et al. (2013) integrated the analytical breaker impulses proposed by Sullivan et al. (2004Sullivan et al. ( , 2007) ) in time and horizontal dimension to obtain the vertical profile of monochromatic breakers.Janssen (2012) suggested a depth dependency with a maximum and a vanishing first derivative at the sea surface: where z ≈ 0.5H s prescribes the vertical gradient behaviour of the wave-induced stress.In this study, we use the shape function given in Eq. ( 25) to calculate the vertical profile of the wave-induced momentum redistribution term into the water column.

Boundary conditions
At the surface, the wave-modified boundary conditions for the velocity components and passive tracers are specified as where Q t is the surface heat flux, E v − P r is evaporation minus precipitation, C p is the water heat capacity at constant pressure, S 0 is the salinity at the sea surface, and τ in and τ ds are the wave-induced stress and dissipation from Eqs. ( 11) and ( 12).τ a = (τ x wind , τ y wind ) defines the total atmospheric stress: where U 10 is wind speed vector at 10 m height, C D is the drag coefficient, and u * a is the air-side friction velocity vector.Across the air-sea boundary, due to the continuity at the interface, the air-side stress equals the water-side stress, giving where u * w is the water-side friction velocity vector.

Energy equation
The mechanical energy equation of the fluctuating motion is given by where e = 0.5(u 2 + v 2 + w 2 ) = 0.5q 2 and e = 0.5(u 2 + v 2 + w 2 ) are the TKE and its fluctuation, respectively, ρ is the fluctuating density, P curr is the shear production by mean shear, here equal to the total mechanical production, P mech , the second term is the buoyancy production/dampening term, the third term specifies the transport of TKE, and the fourth term presents the contribution of pressure transport in the TKE budget equation.The dissipation rate of TKE, ε, is calculated from the Mellor and Yamada (1982) scheme using where B = 16.6 is a dimensionless constant, and the mixing length in the upper ocean is approximated using Eq. ( 22).

Langmuir turbulence
In the wind-driven surface mixed layer, Stokes drift can influence the mean shear by the means of Coriolis-Stokes forcing and the energy budget in terms of Stokes production.Furthermore, its interaction with the small-scale vorticity within the turbulent flow via the vortex force leads to the generation of Langmuir cells that are fully turbulent eddies.Thus, the Langmuir cells can be referred to as Langmuir turbulence (McWilliams et al., 1997) and will occur when ∂u ∂z ∂u s ∂z ; and u 2 , v 2 w 2 .
Adding wave effects to momentum Eqs. ( 16) and ( 17) results in modifications to both shear production and transport terms in the energy equation (Eq.31).For example, the production term is modified as where P Stokes is the Stokes shear production term for one of the wind-driven turbulent mechanisms that is defined by It should be noted that for the horizontally homogeneous 1-D ocean models, the inclusion of the vortex force only affects the Coriolis terms, with no additional contributions in contrast to the 3-D wave-modified ocean circulation models (McWilliams et al., 1997).Therefore, in our simulations, the Stokes drift does not modify the TKE transport term.

Wave breaking
The correlation between the sea surface pressure and the water surface slope imposes an energy flux from the wind to the wave field.For a fully developed sea, the wind energy input to waves is balanced by energy flux from the breaking of waves (S ds S in ).This energy flux may be injected into the water column.The pressure-correlation term, together with the vertical transport of TKE, was parameterized using a turbulent diffusion term by Mellor and Yamada (1982).
It is also possible to parameterize the pressure-correlation term by the dissipation of surface gravity waves.Following Janssen (1999), we assume that the pressure-velocity correlation term prescribes only the work done by the wave breaking on the sea surface.Then, using potential flow theory and kinematic wave boundary condition at the sea surface ( w ∂η/∂t at z = 0), the pressure term over the wave frequency domain is determined by where p d is the wave-induced dynamic pressure from Eq. ( 6), η is the surface elevation, and angle brackets denote averaging over the wave period.The wave growth as a result of wind forcing in the angle brackets can also be related to the wave energy spectrum: where γ is the growth rate of waves by wind.Hence, in the presence of surface wave breaking, Eq. ( 34) may be expressed as Alternatively, Kraus and Turner (1967) modelled the surface flux of TKE (term 1 in Eq. 31) induced by the wave breaking proportional to the cube of air-side friction velocity alongside a proportionality factor.This wave energy factor is a dimensionless parameter (in the range of 2-5) that may depend on various dimensionless variables such as sea state, atmospheric and oceanic stability, latitude, and aerodynamic surface roughness.By converting this energy flux in terms of the water-side friction velocity, the wave energy factor, β, ranges from 50 (young and fully developed waves) to 150.The Craig and Banner (1994) parametrization prescribes the surface flux of TKE, oc , as with β 100 and incorporated it into the energy model by modifying the surface flux boundary condition (by definition, oc is negative in Eq. 35).
To calculate the breaking wave production term in terms of the pressure-velocity term, we need to determine the vertical distribution of oc in the water column.In the case of dissipation of energy by wave breaking, the profile of energy injection may be prescribed by exp(−2k|z|), suggesting that the profile follows the same vertical pattern as the Stokes drift.However, energy injected in the water column over a certain depth may be linked to the sea surface wave height (Terray et al., 1996).Therefore, we apply the profile for the injection of energy due to wave breaking as exp(−|z|/z ), where z ∼ H s .Moreover, P wb should be distributed over the water column such that Hence, by including the contribution of the wave breaking production, the total mechanical production in the TKE budget equation is modified as P mech = P curr + P Stokes + P wb .

Energy boundary conditions
By including the wave breaking production source term, there is no mechanical energy conversion at the sea surface.Thus, the upper boundary condition for TKE becomes where K q = S q lq is the vertical diffusion coefficient for the turbulent kinetic energy, S q is an empirical turbulence coefficient, and z 0 denotes the water-side roughness length (Sect.3.4).When excluding the contribution of P wb in the TKE budget Eq. ( 18), the non-zero TKE flux at the surface should be used for the energy injection via breaking waves: where oc can be estimated using S ds from Eq. ( 35), or alternatively via the Craig and Banner (1994) parametrization using Eq. ( 36).
When the water depth (here about 130 m) is substantially greater than the Stokes depth and the Ekman scale, the bottom effects can be disregarded from the numerical study of the near-surface turbulence dynamics.Thus, we use the zeroflux boundary condition at the bottom in all of our simulations: where D is the total water depth.Now, we present more details about oc when using the Craig and Banner (1994)-type parametrization (Eq.36).Terray et al. (1996) (hereafter T96) obtained experimental evidence that β cannot be treated as a constant, but rather depends on the wave age.For wave ages, A w , greater than about 10, β = 150 was consistent with their observations.Mellor and Blumberg (2004) obtained, using the observations of T96, the following expression: (39) Kantha and Clayson (2004) parameterized this quantity as in which β = 100 corresponds to A w ≈ 26.Nevertheless, when wave information is available, instead of relating the surface TKE flux to the wave stress, specifying the flux to the wave energy dissipation using Eq. ( 35) is recommended, provided that the dissipation source term could be prescribed in a proper way.In this study, we will assign the appropriate values of β wherever it is required.

Water-side roughness length
The water-side roughness length scale, z 0 , depends on the sea surface variability and the upper ocean turbulence structure.T96 parameterized z 0 from the measured significant height of waves as z 0 = 0.85H s .Soloviev and Lukas (2003), based on TOGA COARE data, suggested that z 0 = 0.6H s gives a better fit with their observations.Lewis and Belcher (2004) parameterized z 0 by where q = 0.01-0.04denotes the current to wind relative strength.Mellor and Blumberg (2004) suggested the following Charnock-type relation: using the length scale definition of l = κ max(|z|, 0.85H s ).Rascle (2007) parameterized the wave age in Eq. ( 41) as A w ≈ 30tanh(u * r /|u * w |), where they set the reference friction velocity u * r = 0.020, resulting in the correction of a wave height prediction at mid-latitude.Mellor et al. (2008) obtained the following formulation for the wavy and rough surface: and for the smooth surface as where ν is the molecular viscosity of water.Although the appropriate choice of z 0 is still under discussion, we use the parametrization by Rascle (2007) in all simulations with wave forcing.

Terray et al. (1996)
The vertical profile of ε in the absence of wave breaking is usually scaled according to the law-of-the-wall (LOW): where henceforth u * w = |u * w | is the magnitude of the waterside friction velocity (Eq.30), and z is the distance from the sea surface.When wave breaking is important, scaling of the dissipation rate using Eq. ( 44) no longer holds near the wavy surface boundary.T96 assumed that the vertically integrated ε in the water column is approximately equal to the mean energy input flux from the wind to the waves.They suggested a scaling of ε using the estimated wind energy input, F k , significant wave height, H s , and z: where F k = − oc , D is the water depth, and z t denotes the transient depth (Jones and Monismith, 2007).

Teixeira (2012)
Teixeira (2012) used the rapid distortion theory to investigate the characteristics of turbulence in the upper ocean.He attributes the high dissipation levels near the sea surface in the presence of surface gravity waves to distortion of turbulence eddies by the mean shear of wind-induced current and by straining by the Stokes drift.The dissipation enhancement arises from the same instability mechanism that is the source of Langmuir circulation, along with vigorous amplification of the turbulent shear stress.For the inviscid and non-rotating equations of motion, linearized with respect to the turbulent quantities when the surface wind and the Stokes drifts are roughly in the same direction, he obtained the evolution of dissipation in one eddy turn-over time, T L , as where U = |u| and U s = |u s |.He also provided an expression for the normalized dissipation rate by the significant wave height and the energy flux into the wave as where ζ = k w H s denotes wave steepness; c is a constant (≈ 0.65); k w is the wavenumber; and z = |z|/H s is the nondimensional depth.A(z ) is given as follows: , where La t denotes the turbulent Langmuir number that describes the importance of the Langmuir circulation on producing turbulence: (49)

Huang and Qiao (2010)
According to Huang and Qiao (2010), for conditions when the dissipation rate is balanced by the production terms in the TKE equation, the shear terms can be parameterized by where a 1 is a non-dimensional constant associated with the surface gravity waves and prescribed by regression against the observations of ε collected by Anis and Moum (1995) as where λ is wave length and β w is a dimensionless constant between 0 and 1. Huang and Qiao (2010) validated the skill of their parametrization by favourable comparisons to the observations by Anis and Moum (1995), Wüest et al. (2000), and Osborn et al. (1992).

Environmental forcing and waves
The meteorological direct forcing recorded from Vigra station and wave field inferred from MATS-corrected pressure data are shown in Fig. 4 for the period of experiment from November to the end of December 2011.To estimate wave field, we apply linear-wave theory to transfer the MATS-corrected pressure data to the surface elevation signal (Sect.2.4).Here, we ignore the Doppler shift introduced by the mean current and did not attempt to add higher-order nonlinear terms in the linear theory.To illustrate the uncertainty in using the linear theory, we assume Pierson and Moskovitz (1964) empirical spectrum for the wave energy to estimate significant wave height.The pressure depth-attenuation relative error (Eq.A2) for estimating H s is shown in Fig. 3a, together with the MATS depth variations throughout the experiment (sorted by the measured wind speed) for wind speed ranging from 5 to 22 m s −1 .The error in wave measurements from the subsurface pressure sensor increases with decreasing wind speed.For the energetic wind and wave forcing conditions, the relative error decreases substantially, giving more accurate measurements of surface gravity waves.
A dedicated wave measurement surface buoy, a wave scan (WS) system, was deployed next to MATS.The WS failed to record 2 days into the deployment.A comparison between the estimated significant wave height using corrected MATS pressure data and a 2-day record of H s (from 25 to 27 October 2011) from the WS are shown in Fig. 3b.There is a very good agreement (correlation of r 2 = 0.82) between these two measuring systems.Note that early in the record the wind speed varied between 3 and 7 m s −1 (Fig. 4b) when the relative error is large, suggesting that we can successfully measure the sea surface wave energy spectrum and the corresponding bulk parameters from MATS with better confidence during the energetic wind and wave conditions using pressure and accelerometer sensors.
The deployment period includes two cyclone events with maximum wind speeds above 15 m s −1 , associated with significant wave heights exceeding 12 m (Fig. 4b and d).Sea surface temperature (SST) and atmospheric pressure at the sea level are extracted for the selected period and location from the European Centre for Medium Range Weather Forecast archive.Figure 4a shows the atmospheric pressure, indicating an intense low pressure system passing over the measurement site.SST exceeds the air temperature, T a , on several occasions (Fig. 4c), suggesting convection, especially during the first major wind event on days between 337 and 340.The mean difference between SST and T a during this period is about 2.9 • C. SST decreases continuously throughout the experiment, and also shows abrupt drops associated with the major wind event (Fig. 4c).The wave field is characterized by waves with mean periods, T m , ranging between 7 and 11 s (Fig. 4d), corresponding to wavelengths between 40 and 500 m and significant wave height, H s , between 1.5 and 12 m.Such conditions are favourable for injection of a considerable amount of turbulence by both wave orbital motions, shear-induced turbulence, and wave breaking-induced turbulence.The wave conditions based on the linear wave theory have a mean non-dimensional depth k p D = 0.12, where k p is the peak wavenumber (Fig. 5b).The wave age (peak wave phase speed normalized by the wind speed at 10 m height, c p /U 10 ) ranges from 0.3 to 4.3, which shows a combination of young and old seas (Fig. 5a).Donelan et al. (1993) proposed that the inverse of wave age U 10 /c p = 0.83 corresponds to the wave spectrum at full development and the smaller and larger values, respectively, are identified as swell and wind sea.During our experiment, swell occurred approximately 44 % of the time.
Figure 6 shows the time series of the wind stress components, τ wind = (τ x w , τ y w , 0) , the Stokes drift speed at the sea surface, and the turbulent Langmuir number (Eq.49).During the windy period, the Stokes drift speed reaches 0.9 m s −1 (Fig. 6b).The higher wind speed values are usually associated with younger waves (lower values of c p /U 10 ), as shown in Fig. 5a.The turbulent Langmuir number is shown in Fig. 6c.La t is used to indicate the importance of Langmuir circulation (LC) on producing turbulence in the ocean boundary layer.For La t < 0.3, effects become substantial.For 54 % of the experiment, La t was below 0.3, indicating that LC can enhance upper ocean mixing.Note that this parameter is inversely related to the Stokes drift e-folding depth (Fig. 6b and c).23).Both parameters are contoured in log 10 scale in units of m 2 Hz −1 and m s −1 , respectively.
Figure 7a shows the evolution of wave energy spectrum throughout the experiment.During the major wind events when the large mechanical wind energy is transferred to the surface waves (days between 325-340 and 355-361), the increase of the spectral levels is accompanied by an increase in wave peak frequency.As the wind stress increases (Fig. 6a), the wave height and Stokes drift increase accordingly.The Stokes drift evolution therefore shows high correlation with the surface gravity wave variations in the equilibrium range; especially during the major wind events, u s becomes large because of substantial higher frequency contributions from the surface gravity waves (Figs.7b and 6b).The vertical distribution of the Stokes drift into the deeper water is controlled by the lower frequency components of the wave energy spectrum.
The time series obtained by MATS at approximately 8 m below the surface shows energetic motions in response to the environmental forcing.Figure 8a shows the root mean square (rms) of the velocity components throughout the deployment.All three components show elevated values, especially during the two storm periods.Measured pitch and roll are interconsistent and are correlated with the wind forcing.Pitch angle increases to > 17 • , accompanied by a 2 • increase in the roll angle (Fig. 8).The large inclination of the platform relative to the mean flow, hence the large angle of attack, leads to a substantial growth of uncertainty in the measurement of turbulent fluxes (see Sect. 4.3).

Example spectra from MATS
The velocity frequency spectra from a 15 min burst recorded by the Vector ADV are shown in Fig. 9.The measured veloc-  ity spectra are compared to the spectrum of vertical wave orbital velocity, estimated from the motion-corrected pressure data using the linear wave theory.All spectra show energized motions within the frequency band 0.05 to 0.5 Hz as a result of surface waves.The noise level at frequencies > 0.5 Hz is lower for the axial component due to the configuration of the acoustic beams (Fig. 9a).The broad energetic peak induced by wave motions is located within the inertial subrange, typical for the entire data set.At higher frequencies, immediately after the wave-affected frequency band, the noise level is reached, rendering it difficult to infer turbulent fluxes from the inertial subrange.We attempted to identify the inertial subrange for each spectrum of the axial component of the velocity, u, v and w by evaluating the least squares best-fit slope in the log-space, and failed to detect slopes significantly close to the −5/3 slope of the inertial subrange.This is the case for nearly all vector bursts throughout the experiment.Thus, we do not use vector data to measure turbulent fluxes using the eddy correlation technique or inertial subrange method in this study, but rely on the shear probe data at higher frequencies (Fer and Bakhoday-Paskyabi, 2014;Bakhoday-Paskyabi and Fer, 2013b).
The integrated signal from the shear probe can be compared to the the vertical component of velocity measured by the vector.Figure 9b shows a comparison between the spectrum of integrated shear probe, w shear , and the spectrum of vertical velocity component, w for a 15 min burst.This comparison represents the high correlation of the two spectra, especially at frequencies between 0.05 and 0.5 Hz when they were energized by the wave motions.There is a time-lag between w shear and w (calculated using cross-correlation in time space) that can be described as the time for an eddy to be advected toward the shear probes from the Vector sampling volume (not shown).
The shear frequency spectra obtained from the same 15 min burst are shown in Fig. 10a together with the shear spectra of the vertical component of the wave orbital velocity w o .The most energetic motions and variances occur at a period of about 10 s, which rapidly decay with frequency.There is also a high correlation between the shear spectra of the first shear probe, SH1 (∂w/∂x) and the wave-induced shear in the frequency range between 0.05 and 1 Hz.An inspection of the data from SH2 (∂v/∂x), both in terms of spectral shape and energy, shows that the probe was damaged throughout the deployment.Thus, we discard this probe data from the analysis.Figure 10b shows the spectra from the accelerometer data from the motion sensor and the vibration probe.To meet the equivalent units with the shear spectrum, we rescale the accelerometer and vibration data and describe the variance in the low frequency range in the wave band (GAz) and in the high frequency range (> 1 Hz, VAz), respectively.Here, only the vertical components (GAz and VAz) are shown; however, the spectra from the other components are similar.There is significant acceleration in the wave frequency range, on the same order as the shear probe signal.The frequency band of 1-20 Hz in the shear spectra (typical of all bursts) is less affected by the wave motions and the high frequency vibration.Therefore, in calculation of the dissipation rates, the frequency range from 1 to 20 Hz is used.This ensures that the signal is relatively uncorrupted by the substantial variance in-  2006) algorithm (SH1-c).Also shown is the empirical Nasmyth's spectrum for ǫ = 10 −6 W kg −1 .The shaded parts in each panel mark the frequency range used in calculating the dissipation rate from the shear probes, relatively unaffected by the wave motions.

41
Figure 10.Frequency spectra for a 15 min burst recorded on 2 November 2011, 13:35 UTC.(a) Spectra inferred from shear probes (dw/dt, SH1, and dv/dt, SH2) and the equivalent wave orbital velocity spectrum, dw o /dt, inferred from irrotational wave theory using the measured surface wave spectrum.(b) Spectra from shear probe SH1 compared with the vertical acceleration spectra from the gyro (GAz, representative for motion below 1 Hz) and the vibration sensor (VAz, representative of motion above 1 Hz).(c) Spectra from shear probe SH1 after removing the coherent acceleration signal using the Goodman et al. (2006) algorithm (SH1-c).Also shown is the empirical Nasmyth's spectrum for = 10 −6 W kg −1 .The shaded parts in each panel mark the frequency range used in calculating the dissipation rate from the shear probes, relatively unaffected by the wave motions.
duced by the wave orbital velocities, low frequency platform motion, and the high frequency noise due to vibration.
The response of MATS to wave forcing is further illustrated using the spectra of uncorrected pressure, P , vertical component of accelerometer, a z , and MATS pitch, θ, shown in Fig. 11a.All three spectra show a broad peak around 0.15 Hz with high coherence at frequencies between 0.05 and up to 1 Hz, corresponding to the wave field characteristics passing over the MATS pressure sensor.The peak in P and a z spectra is Doppler shifted as a result of platform waveinduced motion.Figure 11b is an example of the spectral density estimate of the raw pressure spectra, S P , the surface elevation spectra with no tail, S η ntail , the predicted surface elevation spectra using corrected pressure data together with tail correction, S η , and the corrected pressure spectra, S P cor .This figure shows that determining f cut and an appropriate shape for the diagnostic tail beyond the cut-off frequency are crucial in measuring the wave energy spectrum.

Data selection
The shear probe measurements are inaccurate when the angle of attack (AOA) of the flow exceeds about ±20 • (Osborn et al., 1992).The ADV Vector is rigidly fixed to the platform and has one axis always parallel to the shaft of the shear probes, allowing for reliable estimates of AOA.In a mean flow weaker than the wave orbital velocities, the turbulent eddies will not be advected past the instrument and the shearprobe data can be contaminated.The ratio R of the mean flow to the wave-induced flow, and AOA are calculated to detect suspicious data segments following Bakhoday-Paskyabi and Fer (2013b); Fer and Bakhoday-Paskyabi (2014).Pitch and roll behaviour of the platform, AOA, and R are used for quality screening of ε.A data segment is ignored when any of the conditions, R < 1.2, |AOA| > 7 • , mean roll greater than 5 • or rms roll greater than 1 • are met.After data screening, only 11 % of the data are available for the analysis.Figure 12 shows different quality control (QC) criteria applied for dissipation rate measurements.

Numerical calculations
In the following, we compare the observations with the results obtained from the one-dimensional mixed layer model described in Sect.3. Furthermore, we neglect the stratifica- Table 1.Summary of different simulation runs.The run without wave forcing is NW.Wave-modified simulations include runs with (1) Craig and Banner (1994) wave breaking parametrization (CB94), (2) Coriolis-Stokes forcing, wave-induced momentum redistribution term, and wave-modified boundary condition for momentum (Eq.26) (Jenkins, 1987) (J87), (3) Stokes shear production (SSP), (4) Wave-turbulence interaction (WTI) (Eq.50), and (5) wave breaking production (WBP) effects (Eq.37).tion effects in all simulations due to the lack of hydrographic profiles.The model is set up for the location of MATS, and it is assumed that all measured data used to force the model are representative of this location.We use a temporal resolution of t = 60 s, and a vertical non-equidistance resolution with a slight zooming to the sea surface.To assess wave forcing effects on upper ocean turbulence structure, we perform different simulation runs (Table 1).The model is run for (1) an unmodified Mellor-Yamada model and with no wave breaking parametrization (NW), (2) wave-modified by including Coriolis-Stokes forcing, Stokes shear production, and activating Craig and Banner (1994) wave breaking parametrization (WW1), (3) modifications by switching on the wave breaking production source term, Coriolis-Stokes forcing, and Stokes shear production (WW2), and (4) the same configuration as those used in the WW2 run, but the SSP is replaced by the wave-turbulence interaction production term  24), (b) the Stokes shear production, Eq. ( 33), (c) the wave breaking production, Eq. ( 37), based on wave energy dissipation with δ = 0.5, and (d) waveturbulence interaction source term, Eq. ( 50), during the course of the experiment.All parameters are shown in logarithmic scale in units of W kg −1 .

Run
Eq. ( 33) (WW3).In all three wave-modified runs, we use the surface roughness length introduced by Rascle ( 2007), and we apply wave-current interaction modifications for momentum equations (Eqs.24 and 26, respectively).Additionally, we use S ds to build a surface flux boundary condition for TKE in the Craig and Banner (1994) parametrization.
More details about wave energy dissipation and wind energy input source terms can be found in Appendix C.However, we briefly look at the time-depth evolution of waveinduced source terms in Fig. 13 to provide more insight into the wave-turbulence and wave-current interactions.
Figure 13a shows the magnitude of a wave-induced momentum redistribution term, |F ds |, calculated from measured wave energy spectra using Eqs.( 24) and (C2), in which the empirical coefficient δ is assigned to 0.5 to include quadratic wavenumber contribution in the wave energy dissipation calculation.At high winds with wave age less than 0.83 and La t ≤ 0.3, more momentum is forced into the deeper water as a result of wave forcing.Figure 13b shows the Stokes shear production calculated from Eq. ( 33) with a substantial temporal correlation with the wind and wave forcing.The SSP is one of the most important manifestations of wavecurrent interactions in the TKE budget equation, and plays a key role in the enhancement of ε in the upper ocean boundary layer.The SSP depends on the correlation between turbulent fluxes and Stokes drift gradients, therefore, the alignment of Stokes drift with the vertical momentum fluxes determines how this process influences the upper ocean turbulence structure.Furthermore, in the case where currents and waves are not aligned and move in opposite directions, the SSP can enhance the wave breaking processes near the sea surface (Melville, 1996).The wave breaking production term (Fig. 13c) is more energetic near the sea surface but penetrates shallower than the depth-distribution of other waveinduced source terms.The visualization of wave-turbulence interaction in Fig. 13d shows a similar time-depth evolution to SSP with enhanced energy near the sea surface.
The sensitivity of predicting upper ocean turbulence variability regarding the choice of various wave-modified models is shown in Fig. 14.The variability in the dissipation rate is highly correlated with the atmospheric forcing throughout the period of the experiment, suggesting that the mechanical forcing dominates, except for some periods with strong convection (Fig. 14a and c).We cannot investigate the convective conditions as a source of turbulence near the sea surface with the present data set due to the lack of hydrographic profiles.Including wave forcing shows that the dissipation rate is enhanced not only near the sea surface due to breaking wave effects, but also throughout the active mixed layer (Fig. 14c, d and e).All runs with wave forcing predict the elevated ε in response to the wave and wind forcing.The only difference appears when both Stokes shear production, wave-turbulence interaction, and wave breaking production terms are included (WW2 and WW3), especially for the periods before and after the first storm.Fairly good correlation with the measured ε, (Figs.16 and 17) lends some confidence on the modelled dissipation.Hence, we use the results of WW2 as representative of the runs with wave forcing in the following.
After comparing different numerical predictions of ε using 4 wave-modified runs, we now illustrate how the turbulence field evolution is sensitive to the choice of wave breaking parameterization (Sect.3.3.2).In Fig. 15, we plot the time-depth distributions of ε for different parameterizations of P wb .Figure 15a, b, and c show the energy dissipation rate distribution using (i) Eqs. ( 37) and ( 35) with δ = 0.5, (ii) Eq. (37) using wind energy input parametrization, and (iii) using Craig and Banner (1994) parametrization introduced in Eq. ( 36) with β = 100.Unsurprisingly, all simulations give the TKE maxima during the most energetic wind conditions with reasonable depth-dependent agreement in the vertical distribution of TKE.However, there are significant discrepancies in the magnitude of variability of ε between the simulation run with parametrization (i) and those predicted by (ii) and (iii) at the wave-affected surface zone.It should be noted that for β = 100 and extending the wave energy dissipation as a linear function of wavenumber (i.e.δ = 0), runs with parameterizations (i) and (iii) give a very good agreement in prediction of ε (not shown), suggesting a sensitivity of predicted ε to the choice of these empirical coefficients.Furthermore, it can be inferred from this plot that one can apply successfully the Craig and Banner (1994) parametrization for P wb , when there is no wave information.
Figure 16b compares the time evolution of the observed and modelled TKE dissipation rate with NW and WW2 simulation results together with the T96 prediction.Wave effects are highlighted during the major wind events, especially in days 328.5 and 358.5 when the significant wave height and wind speed are elevated.Although there is a high correlation between the modelled (WW2) and the measured dissipation rate, the model cannot properly reproduce ε, especially when wind speed drops below 6 m s −1 .The discrepancy can be attributed to small-scale variability due to intermittency of turbulence, unresolved convection, and likely overestimation in ε measured by MATS.The latter is mainly due to wave-induced platform motion, incomplete separation between wave orbital shear and shear probe signal, and uncertainty in using Taylor's hypothesis.Note that during stormy periods H s is on the order of 10 m, and in 60 s bursts that we use to extract dissipation, the root mean square (rms) of measured pressure reaches 3 m.Close to the sea surface, in the vicinity of strong forcing, there will be large vertical gradients in turbulence properties, for instance ε.Thus, the model output at a constant level, say 8 m from surface, cannot be directly compared to the MATS time series.Hence, uncertainty in the measurement depth is another resource of discrepancy between the model results and the measured ε.

Comparison with other data sets
We compare modelling results and measured dissipation rates with data from the literature.Although using appropriate scaling near the sea surface is a controversial topic due to its complexity as a result of intermittent wave breaking and other underlying unknown processes, we use the T96 similarity in the surf zone in which the dissipation rate is normalized by z, H s and F k ; we retain this scaling for all other data sets.Note that this choice of scaling will increase the uncertainty of analysis when the wave breaking is not an important feature in the air-sea interface.We extract the dissipation data from Burchard (2001), Gerbi et al. (2009), Jones and Monismith (2007), and Feddersen et al. (2007)   ocean boundary layer.After rejection of a large number of observed ε in the QC stage, the remaining high-quality dissipation data can properly cover the various wind and wave conditions.We scale the numerical runs and parametrization results with significant wave heights between 3 and 7 m, the typical wave conditions when ε is available.Furthermore, we identify different individual wind events in Table 2 to provide more information on MATS ability for measuring ε near the wavy surface.The simulation results, together with parameterizations of turbulent mixing (Sect.3.5), are presented in Fig. 17 where each panel presents a data set (bullets) indicated in the label, and the MATS-derived ε at six chosen wind events.
Figure 17a shows that our averaged normalized dissipation rate (with β = 100 in F k = βu 3 * w ) is in fair agreement with those observed by Burchard (2001), but with a slightly elevated energy.However, the wave-forced model, the Huang and Qiao (2010) parametrization, and the Teixeira (2012) model similarly overestimate, consistent with the MATS data.The NW run underestimates the measured ε.The Teixeira (2012) model successfully predicts our measured dissipation, with a better fit than those estimated by WW2 simulation runs for events 2 and 3 when La t is roughly less than 0.3.Our measurements of dissipation rate are thus consistent with the upper ocean turbulence structure with values much larger than those predicted from the LOW.
The data set of Gerbi et al. ( 2009) is presented in Fig. 17b.To be consistent with Gerbi et al. ( 2009), we use β = 168.Only for events 1 and 3 is our measured dissipation in marginally good agreement with ε observed by Gerbi et al. (2009).Both the model of Teixeira (2012) and the WW2 curves fit relatively better to the MATS.In Fig. 13c (β = 54), it can be seen that our measured ε exhibits more scatter than those measured by Jones and Monismith (2007); the WW2 overestimates the measured ε, while the T96 model is in better agreement with the observations of Jones and Monismith (2007), and the model of Huang and Qiao (2010) provides better fit with majority of events than those obtained with other techniques.The possible reason for these discrepancies between model results and the data of Jones and Moni-smith (2007) may be related to the wave refraction effects (by bathymetry) in their measurements.The final observational data set is from Feddersen et al. (2007), who address the nearshore vertical structure of dissipation rate of TKE (Fig. 13d).In order to fit the model of T96 to their observations, Feddersen et al. (2007) used β = 250, which we adopt in Fig. 13d.This figure shows that all model curves are in approximately good agreement with the observations of Feddersen et al. (2007), and again WW2 together with the Teixeira (2012) scaling, agree most favorably with our measurements.

Summary and conclusion
Observations were made in the near-surface layer, at approximately 8 m depth from surface, off the coast of Ålesund in Norway, for a duration of 2.5 months in late 2011.The time series of the dissipation rate of turbulent kinetic energy, ε, and the estimates of surface gravity waves are analysed.The instrument is a moored system (MATS) designed for measuring continuous time series at a fixed depth in the upper ocean boundary layer.The measurement period covers the passage of two low pressure systems, and substantial wind and wave forcing.
Measuring ε using shear probes on MATS is difficult in the wave-affected surface layer because of wave-induced platform motions, wave orbital velocity disturbances, bodyinduced motions and body vibration effects.We constrained measurements of ε by angle of attack, the platform's roll angle, and the ratio between the mean flow velocity and the wave-induced flow to avoid violation of Taylor's frozen turbulence hypothesis.After data screening, only 11 % of the data were available for analysis.The shape of the shear spectra consists of elevated variances in the frequency band between 0.05 and 0.8 Hz that is affected directly by the surface wave motions.The shear spectra are relatively uncontaminated in the frequency range of 1-20 Hz, which is used to calculate ε.
To study turbulence variations near the sea surface, we employed a second-order turbulence closure scheme of Mellor-Yamada level 2.5 type for the near-surface mixed layer.We incorporated near-surface processes to this model by including wave breaking effects and the interaction of the mean current shear with the Stokes drift.The sea state data needed to force the model were inferred from the high-resolution pressure time series measured by MATS.Using linear wave theory, the pressure data, corrected for acceleration and attenuation at measurement level, were used to calculate sea surface elevation, wave orbital velocities, and other important wave bulk parameters, including significant wave height, and mean wave period.Numerical calculations were carried out for cases without wave forcing (NW) and three different cases of wave forcing including combinations of different terms such as Coriolis-Stokes forcing, wave breaking and Stokes shear production.The results from a forcing including wave breaking, Coriolis-Stokes forcing and Stokes shear production (WW2) compared favourably with the observations.
Additionally, we used various scaling laws proposed for the depth dependence of ε such as the classical wall layer, the scalings proposed by Terray et al. (1996) and Huang and Qiao (2010), and the model by Teixeira (2012) which includes the influence of Langmuir turbulence.The wall-layer scaling underestimated the measured ε by up to two orders of magnitude.The Teixeira (2012) method, which depends on the shear induced by Stokes drift, agreed better with the measured dissipation.In this experiment, the wave field was separated into the wind sea and swell components in several occasions, however, the high-quality dissipation measurements were scarce in very energetic wind and wave conditions.The bulk of the measurements were therefore better represented using the Huang and Qiao (2010) and Teixeira (2012) methods than the Terray et al. (1996) scaling which works well when wave breaking is important.Furthermore, while the NW run results deviated significantly from the measured ε, the WW2 simulation indicated good agreement with the measurements.
There are several important caveats in our measurements and simulations of ε near the wavy surface, including wave and turbulence separation, flow distortion and contaminations induced by platform motions, uncertainties in the estimation of wave energy spectrum, in using a non-directional wave energy spectrum, in wave forcing parameterizations, in atmospheric forcing, the intermittency of turbulence, and the lack of proper coupling strategy between wave, current, and turbulence.Nevertheless, comparisons with other data sets available from literature lend confidence to our dissipation measurements and wave-forced simulations.When scaled using the wave parameters as suggested by Terray et al. (1996), comparisons showed that the dissipation is enhanced near the sea surface in the presence of surface gravity waves, and the average dissipation inferred from MATS agreed well with other data sets.

Figure 1 .
Figure 1.Picture of MATS during the deployment together with a close up of the sensors (inset).The inset on the right shows the MATS mooring in the water column (not to scale).

Figure 2 .
Figure 2. Deployment site together with the isobaths at 100 m intervals.The position of MATS (filled rectangle) and the meteorological station Vigra (filled bullet) are marked.The inset shows the location in Norway.

Figure 3 .
Fig. 3. a) Pressure depth-attenuation relative error, RHs(z, U10), in significant wave height, Hs, calculations based on the linear small-amplitude wave theory.RH s is a function of depth, z, and wind speed at a reference height of 10-m, b) Comparison between Hs measured by the Wave Scan (WS) in the vicinity of MATS (plus markers) and estimated from the MATS corrected pressure sensor data (solid line).The WS measurements are available for only 2 days.

Figure 4 .
Figure 4. Evolution of (a) air pressure; (b) wind speed U 10 and direction (thin red line), positive clockwise from north, referring to the direction from which the wind is originating; (c) sea surface temperature, SST, and air temperature, T a , (thin red line); and (d) the significant wave height, H s , and mean wave period, T m , (thin red line) measured by motion-corrected pressure data for the duration of the experiment.
Fig. 5. a) Time series of the inverse of wave age.Horizontal red dashed line is the separation between swell and wind sea.b) Time-averaged wave energy spectrum estimated from the motion-corrected pressure data during the deployment.Dashed red line highlights the peak frequency, fp.
Fig. 6. a) Time series of wind stress components τ x w and τ y w , b) the surface Stokes drift speed, and c) the Langmuir turbulence number, Lat.The red dashed line indicates the threshold below which Langmuir turbulence becomes important.

Figure 7 .
Figure 7. (a) Wave energy spectrum calculated from MATS' motion-corrected pressure data, and (b) the depth-time evolution of the Stokes drift calculated from Eq. (23).Both parameters are contoured in log 10 scale in units of m 2 Hz −1 and m s −1 , respectively.

Fig. 8 .Figure 8 .
Fig. 8. Time series of: a) the root-mean-squared (rms) current velocities, urms, vrms, and wrms, determined from ADV with sampling frequency of 16 Hz for the period of experiment, and b) pitch, θ and roll, φ, measured by the motion sensors.

Fig. 9 .Figure 9 .
Fig. 9. Frequency spectrum of u, v, and w velocities measured by the ADV for a 15-min burst starting from 02 Nov 2011, 13:35 UTC.The spectra are band-averaged in frequency in log 10 (∆f ) = 0.01 bins.Also shown is the vertical component of wave orbital velocity spectrum, wo, and b) comparison between the integrated spectra of the shear probe 1 data, wshear (black), and the vertical velocity component, w measured by the ADV (red) for the same burst.

Fig. 10 .
Fig.10.Frequency spectra for a 15-min burst recorded on 2 November 2011, 13:35 UTC.(a) Spectra inferred from shear probes (dw/dt, SH1, and dv/dt, SH2) and the equivalent wave orbital velocity spectrum, dwo/dt, inferred from irrotational wave theory using the measured surface wave spectrum.(b) Spectra from shear probe SH1 compared with the vertical acceleration spectra from the gyro (GAz, representative for motion below 1 Hz) and the vibration sensor (VAz, representative of motion above 1 Hz).(c) Spectra from shear probe SH1 after removing the coherent acceleration signal usingGoodman et al. (2006) algorithm (SH1-c).Also shown is Fig. 11.a) Power spectrum of pressure at depth of MATS (thick line), vertical component of accelerometer (dashed line), and MATS pitch angle (thin line) calculated from an arbitrary 15-min segment.The gray region shows the wave-affected frequencies.b) Comparison among corrected and uncorrected MATS pressure spectra, SP cor and SP , respectively, together with surface wave elevation spectra with tail and no tail, Sη and Sη ntail , respectively.The vertical dashed line denotes the cut-off frequency location.

Fig. 12 .Figure 12 .
Fig. 12.Time series of a) the ratio of the mean flow to the wave-induced flow, R, b) the inferred high frequency component of angle of attack, αhf that is calculated by integrating the α in the frequency band 0.5 to 2 Hz at each 60 s long segment, and c) platform's roll angle.Gray regions (dashed lines) indicate the thresholds for good quality measurements during the experiment.

Figure 13 .
Figure13.Time-depth evolution of (a) magnitude of wave-induced momentum redistribution source term, Eq. (24), (b) the Stokes shear production, Eq. (33), (c) the wave breaking production, Eq. (37), based on wave energy dissipation with δ = 0.5, and (d) waveturbulence interaction source term, Eq. (50), during the course of the experiment.All parameters are shown in logarithmic scale in units of W kg −1 .

Figure 14 .
Figure 14.Time series of (a) air-side surface friction velocity, and dissipation rate from the runs (b) NW, (c) WW1, (d) WW2, and (e) WW3 during the experiment.Dissipation rate is shown in logarithmic scale in W kg −1 .

Figure 15 .
Figure 15.The time-depth distribution of dissipation rate from run WW2 by including different breaking wave production (Eq.35) parameterizations using (a) wave energy dissipation with δ = 0.5, (b) wind energy input source term, and (c) the Craig and Banner (1994) parametrization.Dissipation rate is shown in logarithmic scale in W kg −1 .
Fig. 16.a) Time series of wind speed at 10 m height (black line), and significant wave height (red), and b) mean dissipation rate time series at the depth of MATS, measurements (square markers), simulation run NW (blue), run WW1 (black), run WW2 (red), run WW3 (dots), and estimates from T96 (triangle markers).

Table 2 .
Identified events during the period of experiment.The start and end day of the year of each event, duration, and the averaged wind speeds are listed.