Correcting surface wave bias in structure function estimates of turbulent kinetic energy dissipation rate

AbstractThe combination of acoustic Doppler current profilers and the structure function methodology provides an attractive approach to making extended time series measurements of oceanic turbulence (the rate of turbulent kinetic energy dissipation e) from moorings. However, this study shows that for deployments in the upper part of the water column, estimates of e will be biased by the vertical gradient in wave orbital velocities. To remove this bias, a modified structure function methodology is developed that exploits the differing length scale dependencies of the contributions to the structure function resulting from turbulent and wave orbital motions. The success of the modified method is demonstrated through a comparison of e estimates based on data from instruments at three depths over a 3-month period under a wide range of conditions, with appropriate scalings for wind stress and convective forcing.


Introduction
Exchanges of heat, freshwater, and trace gases between the ocean and the atmosphere are critical in regulating the climate and depend directly on the properties of the ocean surface boundary layer (OSBL) (e.g., D'Asaro 2014 ;Franks 2014;Large et al. 1994). The structure of the OSBL depends on turbulent processes that cannot be directly simulated in geographical scale numerical models and therefore have to be parameterized (Burchard et al. 2008;Belcher et al. 2012;Calvert and Siddorn 2013).
Turbulence in the OSBL is widely recognized as being produced by wind-driven surface shear stress, destabilizing surface buoyancy fluxes and (in shelf seas) tidal current shear at the bottom boundary (e.g., Brainerd and Gregg 1993;Simpson 1981). Other surface-driven processes include breaking waves (e.g., Agrawal et al. 1992;Terray et al. 1996), Langmuir circulation (e.g., Thorpe 2004), submesoscale eddies (e.g., Taylor 2016), and swell waves (e.g., Wu et al. 2015). Developing effective parameterizations for such diverse processes requires robust measurements under a wide range of environmental conditions, presenting significant observational challenges.
The structure function method is an established technique for calculating the turbulent kinetic energy (TKE) dissipation rate « from velocity profiles such as those obtained with an acoustic Doppler current profiler (ADCP) (e.g., Wiles et al. 2006;Mohrholz et al. 2008;Lucas et al. 2014;Simpson et al. 2015;McMillan and Hay 2017). The method relates « to the variance of the along-beam turbulent velocity difference evaluated over a range of separation distances. Instrument choice and configuration impose constraints on the data collected, but once configured the ADCP can be deployed to make unattended long-term observations, unlike standard microstructure techniques.
Surface waves induce orbital motions within the water column, the speed of which reduce with depth. The velocity associated with the orbital motions may be observed by the ADCP, potentially affecting the structure function and introducing bias in the « estimates. To date, the structure function technique has typically been applied to observations from sites with small amplitude surface waves or at depths unlikely to be affected by significant wave orbital velocities (Wiles et al. 2006;Lucas et al. 2014;Simpson et al. 2015;McMillan and Hay 2017).
An exception is the application of the technique by Thomson (2012) to obtain « estimates within the crests of breaking waves by mounting the ADCP onto a surfacefollowing Lagrangian float and by necessity limiting the range of separation distances over which the structure function was evaluated. Similarly, in order to measure vertical profiles of « in the near-surface under breaking waves, Sutherland and Melville (2015) adapted the technique by restricting both the range of separation distances and the time-averaging period over which the statistical properties of the structure function were evaluated. Restricting the range of separation distances minimizes the difference in the orbital velocity seen by different ADCP bins, while adopting a time-averaging period similar to or less than that of the waves will result in the wave orbital velocity being treated as a background mean flow.
Working in a shallow-water, wave-dominated environment, Whipple and Luettich (2009) assume that the velocity variance at each depth (calculated over a sampling period much longer than the wave period) is dominated by the wave orbital velocity at that depth. They fit a theoretical vertical profile based on linear wave theory to the observations in order to characterize the effective wave contribution to the structure function over a specified depth range. This is then used to remove the influence of waves and to isolate the much smaller turbulent signal. While this approach explicitly recognizes the contribution of the vertical gradient of the wave orbital velocity to the structure function, it is applicable only in situations where the wave influence dominates the structure function and does not lend itself to more general application.
The aims of this paper are, first, to demonstrate that « estimates made using the standard structure function method with ADCP data are inherently susceptible to bias in the presence of surface waves due to a contribution to the structure function from the vertical gradient in the speed of the associated wave orbital motion; and second, to present a modification to the standard method that addresses such bias. Section 2 briefly covers the underlying theory, demonstrates the standard method's bias using the wave orbital motions under synthetic monochromatic waves, and describes the proposed modified method based on the application of linear wave equations. Section 3 describes a set of long-term field observations from a shelf sea site that was used to test the standard and modified methods. Section 4 uses established similarity scaling approaches to compare the results under differing surface forcing conditions, and section 5 is a discussion of the results.

a. Structure function
The theoretical basis of the structure function technique and its derivation from the Kolmogorov similarity hypotheses is described in detail elsewhere (Sreenivasan 1991;Frisch 1995;Antonia et al. 1997;Pope 2000;Lucas et al. 2014;McMillan and Hay 2017). In summary, the technique assumes that for isotropic turbulence in high Reynolds number flows, an inertial subrange of length scales exists over which there is a conservative cascade of energy from larger to smaller motions. The statistical properties of the longitudinal turbulent velocity fluctuationdu 0 (x, r) [ u 0 (x 1 r) 2 u 0 (x), where u 0 (x) is the along-axis turbulent velocity at location x-vary as a function of the separation distance r.
Invoking Taylor's frozen field hypothesis to allow sampling of the statistical properties of the flow over time, the mean du 0 is related to « for r values within the inertial subrange as where the angle brackets indicate time averaging over a statistically valid sampling period and n is the order of the structure function (Kolmogorov 1991a,b;Pope 2000). The second-order structure function D LL (x, r) is then defined as and for values of r within the inertial subrange D LL (x, r) is related to «(x) as where C 2 is a universal constant of proportionality, frequently taken to be 2.1 based on atmospheric studies (Wiles et al. 2006;Lucas et al. 2014;Simpson et al. 2015), while McMillan and Hay (2017) use 2.0 based on both theoretical considerations and the comparison of « estimates made using the structure function and spectral integral methods. From (3), the second-order structure function exhibits a length scale dependence on r 2/3 , so a least squares linear regression of D LL (x, r) against r 2/3 , at fixed x, gives where A 0 is a measure of the Doppler and instrument noise, and A 1 is the gradient of the linear regression over the range of r evaluated. From (3), A 1 5 C 2 « 2/3 , which then gives an estimate of « at x for the sampling period as When applied to ADCP data, a sampling period of several minutes is typically used, during which multiple individual velocity profiles are collected at a frequency of 1 2 2 Hz. The along-beam velocity data are processed for each beam separately, with u 0 calculated for each bin by deducting its mean over the sampling period in order to remove the mean flow and hence any background shear. The structure function D LL is then calculated from the velocity differences at r based on multiples of the alongbeam bin size. The minimum separation is taken as two bins as a result of the lack of independence in the velocities measured in adjacent bins (Teledyne RD Instruments 2014). The squares of the velocity differences are then averaged over the sampling period as in (2). Using a central difference scheme (e.g., Wiles et al. 2006), D LL is evaluated for each bin for separation distances centered on the bin, with the r values that can be resolved dependent on the bin's position within the range of bins for which the turbulent velocity is available.
A maximum separation distance r max is specified for the regression of D LL against r 2/3 . This should be chosen to include as much of the inertial subrange as possible, although in practice the configuration of the ADCP may restrict the range over which turbulent velocities are resolved. When this is not a constraint, r max must not exceed the upper length limit of the inertial subrange, beyond which D LL is expected to tend toward a constant. The selection of r max therefore depends on both instrument constraints and the turbulent properties of the observed flow.

b. Wave orbital motion
A basic representation of deep-water surface gravity waves is to treat them as sinusoidal, with amplitude A, wavelength l, and period T, giving a radian frequency of v 5 2p/T, wavenumber k 5 2p/l, and phase speed c defined as c 2 5 v 2 /k 2 5 g/k, with g being the acceleration due to gravity.
The simplest model for the motion in the water column below such waves (e.g., Phillips 1977; Simpson and Sharples 2012) is of nonrotational circular motion with a speed at depth z (zero at the surface, positive up) of y max 5 vAe kz .
Over a vertical distance dz around depth z 0 , the difference in the speed of the orbital motion is dy max (z 0 ) 5 vA[e k(z 0 1dz/ 2 ) 2 e k(z 0 2dz/ 2 ) ] ' ky max (z 0 )dz, subject only to the adoption of the small angle approximation that sinh(kdz/2) ' kdz/2, which is valid to within 2% for dz , l/10. Hence, at all depths the vertical difference in the orbital speed varies linearly with the vertical separation distance. As illustrated in Fig. 1, this vertical variation in the speed of the orbital motion will result in a contribution to the structure function even in the absence of turbulence. Under a monochromatic wave, the along-beam velocity measured in the ADCP bins will vary sinusoidally in phase in all bins but with an amplitude that depends on the depth of the bin. Since the sampling period used to determine the structure function is normally much longer than the surface wave period (several minutes vs typically less than 15 s), the mean of the along-beam component of the wave orbital motion measured by any bin is approximately zero and will not contribute to the mean velocity deducted to calculate the fluctuating u 0 . Consequently, u 0 retains the along-beam component of the time-varying wave orbital motion. Any differences in u 0 between bins will be treated as a turbulent velocity variation when calculating D LL , potentially resulting in a bias in the calculated « estimates.
To quantify the potential bias, « values were calculated using wave orbital velocities calculated from linear wave theory for a range of monochromatic waves with amplitudes and periods representative of an exposed shelf sea environment. These synthetic wave orbital velocities were calculated for the bin locations of virtual ADCPs at depths of 20, 35, and 50 m with an upward-looking orientation, sampling via a beam with a 208 beam angle (inclination from the vertical) with 30 bins at a 0.1-m vertical bin spacing and bin 1 centered at 0.97 m from the transducer. The measurement frequency was 1 Hz with a sampling period of 300 s, resulting in 300 velocity profiles.
Assuming waves propagate in the x direction and the ADCP beam in the y 5 0 plane, the horizontal (u) and vertical (w) velocities vary as u 5 vAe kz sin(kx 2 vt) w 5 2vAe kz cos(kx 2 vt) where t denotes time.
The along-beam velocities in each bin were calculated by applying a rotation matrix based on the virtual ADCP beam geometry (Teledyne RD Instruments 2010). The structure function D LL was calculated using a central difference scheme and « estimates were determined for each bin from the regression of D LL against r 2/3 with r max equal to 2.0 m. Beam average « values were calculated as the geometric mean of the individual values for all bins for which the structure function was resolved for all r # r max . Figure 2 shows the « estimates for each of the three instruments for surface waves with amplitudes up to 2 m and periods between 7 and 13 s. The bias in « is more than 1 3 10 25 W kg 21 for an ADCP at a depth of 20 m under waves with an amplitude of 1.8 m and a period of 8 s. Even for an instrument at 50-m depth, swell waves with a period of 11-12 s and an amplitude of 1.6 m could introduce a bias of O(10 27 W kg 21 ), two orders of magnitude above the expected noise floor (Lucas et al. 2014).
The bias in « depends on the difference in the speed of the wave orbital motion over distance r max , which depends on both the amplitude and the attenuation rate of the speed of the orbital motion. Since the attenuation rate depends on wavenumber, the period of the waves contributing most to any bias will typically increase with ADCP depth.
For a spectrum of waves, linear wave theory would suggest that the along-beam velocities observed by the ADCP will be the sum of the wave orbital velocities resulting from the various component waves. While the velocity contribution from each component wave will depend on its surface properties and attenuation rate, each will exhibit the linear variation with vertical separation in (7). The composite wave orbital velocity can therefore also be expected to demonstrate a linear length scale dependency.
Though the leading order water motions associated with the surface waves are periodic and do not affect the time-averaged current profile, surface waves also produce a second-order depth-varying Lagrangian transport in their direction of propagation, the Stokes drift (e.g., Phillips 1977;Ardhuin et al. 2009). Within the structure function calculation, any nonperiodic velocity observed by an ADCP bin is considered as part of the mean flow and removed when the turbulent velocity is calculated. Asymmetric periodic flows, such as the difference between the upper and lower portions of a wave orbital motion that leads to Stokes drift, may result in a nonzero contribution to the mean flow and a contribution to the structure function based on the depth-dependent variation in the periodic motion. The Stokes drift speed decays exponentially with depth at twice the rate of the wave orbital motion (Phillips 1977). It is therefore also expected to exhibit a linear length scale dependence over a limited vertical separation distance.
Exploiting the differing length scale dependencies of the turbulent and wave-related components of the observed velocity offers the possibility of separating these two components of the structure function.
c. Modified methodology to reject impact of wave orbital motion From (1), the nth order structure function varies as r n/3 ; hence, D LL will vary linearly against r 2/3 . By contrast, from (7), the difference in the maximum wave orbital velocity magnitude dy max varies linearly with r; FIG. 1. Schematic of wave orbital motion contribution to D LL . Monochromatic deep-water surface waves of period T p drive irrotational circular motions with speed at z (zero at the surface, positive up) given by y max (z) 5 Ave kz . In the absence of any other motion, the ADCP measures only the along-beam component of the wave orbital motion, such that u(z, t) 5 y max (z) sin(vt), with the velocities being in phase between bins but varying in magnitude with bin depth. The turbulent velocity, u 0 5 u 2 hui, retains the wave orbital motion, since the bin mean over a sampling period hui TTp ' 0. The second-order structure function is the mean of the turbulent velocity variance, h(du 0 ) 2 i, for a range of separation distances; see (2). In the presence of an along-beam gradient in wave orbital motion speed, h(du 0 ) 2 i . 0 for all separation distances, resulting in an unavoidable nonturbulent contribution to D LL .
FIG. 2. Standard second-order structure function method bias in « as a result of wave orbital motion for synthetic deep-water monochromatic waves observed by virtual ADCP at depths of (a) 20, (b) 35, and (c) 50 m. Structure function D LL based on a central difference scheme with regression based on r max ; 2:0 m. Beam average « based on the geometric mean of bins for which D LL is resolved for all r # r max . ADCPs are assumed to have a sampling rate of 1 Hz, a sampling period of 5 min, a vertical bin size of 0.1 m with the first bin centered at 0.97 m from the transducer, and to be upward looking with a 208 beam angle.
hence, from (2) the contribution to D LL varies as r 2 . In the regression of D LL against r 2/3 , the contribution to the structure function from the vertical variation in wave orbital velocity will therefore increase as (r 2/3 ) 3 . The differing rates at which the contribution of the turbulent and wave orbital motion components of the structure function vary with separation distance provides the basis for the modified method. Instead of the standard least squares linear regression of D LL against r 2/3 as in (4), a least squares fit is done to determine the coefficients for the linear model: The modified method essentially assumes that the wave orbital motion and turbulence do not interact and that the associated velocities are simply additive. The contribution to D LL as a result of the vertical gradient in the speed of the wave orbital motion (contained in the A 3 coefficient) can therefore be extracted without affecting the turbulent contribution. Hence, the A 0 coefficient continues to describe the instrument and Doppler noise, and the A 1 coefficient continues to describe the turbulence, with « still calculated using (5). The effectiveness of the modified method was tested by applying it to the synthesized wave orbital velocity data described in section 2b. Figure 3 shows the regression of D LL against r 2/3 for both the standard and modified methods for the instrument at depth 35 m with a surface wave of amplitude 1 m and a period of 10 s. The standard method results in a calculated « of 1:4 3 10 27 W kg 21 and a physically meaningless negative A 0 value of 22:6 3 10 25 m 2 s 22 . By contrast, the A 0 and A 1 coefficients for the modified method correctly reflect the fact that there was no turbulent motion or system noise in the synthesized velocity data.

d. Similarity scaling
To compare the results of the standard and modified methods at different depths and under widely varying environmental conditions, two distinct surface forced regimes with established similarity scalings are considered. The relevant scaling factors are applied to « estimates calculated using both the standard and modified methods to illustrate the conformance of the results from the two methods to the standard scalings. 1) Wind stress forcing. Following Monin-Obukhov similarity theory, a local balance is assumed between « and TKE production based on a constant stress ''law of the wall'' relationship (Anis and Moum 1995;Lombardo and Gregg 1989;Brainerd and Gregg 1993;Lozovatsky et al. 2005;Tedford et al. 2014;Bogucki et al. 2015;D'Asaro 2014). This results in a scaling factor « s given by where u * is the friction velocity in the water, calculated as u * 5 (t s /r 0 ) 1/2 for surface wind stress t s and water density r 0 ; k is the von Kármán constant (0.41); and z is depth (zero at the surface, positive up). Within the mixed layer, but below the region of direct impact from breaking waves (Agrawal et al. 1992;Anis and Moum 1995), « estimates would be expected to scale as «/« s ' 1, with reported values typically in the range 1-2 based on limited duration observations (Lombardo and Gregg 1989;Lozovatsky et al. 2005;Shay and Gregg 1986;Thorpe 2005). 2) Convective forcing. By convention a positive surface buoyancy flux, B 0 . 0, indicates a loss of heat from the ocean surface to the atmosphere, increasing the ocean surface density and creating unstable conditions, leading to convection and an increase in «. Within the mixed layer, but below the Monin-Obukhov length (the depth at which wind stress forcing and convective forcing match), « is expected to be constant, reducing only at the base of the mixed layer when it encounters stratification and contributes to mixing by entrainment (Shay and Gregg 1986;Lombardo and Gregg 1989). Hence, under low wind conditions, « estimates would be expected to scale as «/B 0 ' 1, with reported values based on limited duration observations typically being in the range of 0.5-0.8 under conditions of both sustained and diurnal convection, with some indication of a time dependence as convection becomes established (Anis and Moum 1992;Brainerd and Gregg 1993;Lombardo and Gregg 1989;Shay andGregg 1984, 1986;Thorpe 2005).
Combined scalings incorporating both wind stress and convective forcing have been developed as linear combinations of the scalings for the individual forcing regimes (e.g., Lombardo and Gregg 1989;Tedford et al. 2014). However, the variation in the reported weighting coefficients suggests that the combined scaling may be less robust than the scaling for the individual regimes.
The objective of the current study is not to revisit these scalings but to use them as the basis for comparing the susceptibility of the standard and modified structure function methods to wave-induced bias. The scalings were therefore applied separately to « estimates based on field observations made under the relevant forcing conditions, and the results were compared to a default depth-constant unity reference value.

a. Dataset
The present analysis is based on observations made during the period January-March 2015 from a site in the Celtic Sea. The site has a water depth of ;150 m; is more than 200 km from any coast, removing it from the direct coastal influences; and is over 125 km from the shelf edge, minimizing the impact of any shelfbreak processes. The wave climate included both locally generated waves and remotely generated swell, unaffected by significant shoaling or coastal reflections. Three Teledyne RD Instruments 600-kHz Workhorse ADCPs were deployed on a buoyancy-tensioned mooring attached to a seabed anchor weight. The instruments were all configured in pulse-to-pulse coherent mode (mode 5) (Teledyne RD Instruments 2014) with a sampling frequency of 1 Hz and one ping per ensemble (no ensemble averaging), with a vertical bin size of 0.1 m and bin 1 centered 0.97 m from the transducer. The instruments operated for a 5-min sampling period, followed by a 15-min rest interval, resulting in three sampling periods per hour, each comprising 300 velocity profiles for each of the four beams. The uppermost instrument had a 208 beam angle and was deployed upward looking, the middle instrument also had a beam angle of 208 but was deployed downward looking, while the lowest instrument had a beam angle of 308 and was upward looking.
The mooring rotated with the tide, the depth-averaged current having spring tide maxima of ;0:5 m s 21 with a pronounced spring-neap cycle. The instruments' measurement volumes were centered at mean depths of ;24:0, 42.5, and 52.5 m. Reliable velocity measurements were typically returned for bins 1-30 for the 208 beam angle instruments and bins 1-28 for the 308 beam angle instrument, equating to bin centers at along-beam distances of ;1 to ; 4:2 m from the transducer.
Three additional moorings provided supplementary information used in this analysis. All moorings were located within 1 km of each other throughout the observation period. One of the moorings provided full water column temperature, salinity, and density (Wihsgott et al. 2016). Another was a Met Office Ocean Data Acquisition System (ODAS) buoy, which provided meteorological and wave data, including hourly measurements of average wind speeds and direction plus maximum gust speeds at 3 m above the sea surface based on sampling over a 10-min period; air and sea surface temperature; atmospheric pressure and relative humidity; plus significant wave height and average wave period based on 17.5 min of observations. The third was a U.K. Centre for Environment, Fisheries and Aquaculture Science (Cefas) SmartBuoy, which provided half-hourly sea surface temperature and salinity, plus photosynthetically active radiation (used as a proxy for solar irradiance).

b. Data analysis
Surface stress and buoyancy flux were calculated using the TOGA COARE 3 bulk flux algorithm, taking account of the heights of the instruments on the ODAS buoy (Fairall et al. 2003).
The ADCP beam coordinate u 0 were calculated independently for each bin in each beam by deducting the mean for that bin over the sampling period. Outlier values were identified by comparison with the rms value of all turbulent velocities for all bins and beams in the current sampling period and rejected. Outliers were almost exclusively in the furthest bin for which the velocity was resolved. The second-order structure function, D LL was calculated using a central difference scheme over all resolvable separation distances, r 5 r j Dr, where r j is the separation in number of bins and Dr is the along-beam bin size determined by the vertical bin size and the beam angle. For even number bin separations, r j 5 2, 4, 6, . . . around bin i: where u 0 (x i ) is the turbulent velocity in the bin centered at distance x i from the transducer. For odd number bin separations, r j 5 3, 5, 7, . . ., the average of the two possible combinations was used, so where floor( ) (ceil( )) means round down (up) to the integer.
The D LL values for all bins were used in least squares fit regressions against r 2/3 , to give a beam aggregate « value for the sampling period for both the standard (4) and modified (9) methods. The regressions were repeated for a range of r max values between 0.8 and 3.0 m (the maximum possible values given the instrument configurations). Basic result screening rejected regressions if the coefficients did not produce a strictly increasing result for r . 0. Equation (5) was used to calculate « with C 2 as 2.0. The geometric mean of the individual beam values provided a single representative « data point per sampling period for each instrument, method, and r max value over the 3 months of observations, resulting in approximately 6500 data points for each combination of instrument, method. and r max .
The adjusted coefficient of determination-R 2 adj 5 1 2 (1 2 R 2 )[m 2 1/m 2 (p 1 1)], where R 2 is the unadjusted coefficient of determination; m is the sample size, and p is the number of independent variables in the regression-was calculated for each regression. Using R 2 adj rather than R 2 allows the quality of the fit from both the standard and modified methods to be compared directly, taking into account the additional term in the modified method.

Results
The 3 months of observations included in this analysis cover a wide range of winter conditions. Throughout the period, the water column was negligibly stratified. The B 0 was characterized by a destabilizing heat flux to the atmosphere (B 0 positive) approximately 70% of the time, when the mean flux was 6 3 10 28 W kg 21 and the maximum was 1:9 3 10 27 W kg 21 . Solar irradiance resulted in intermittent diurnal stabilizing (B 0 negative) buoyancy fluxes, centered around midday and increasing in duration and maximum intensity over the period of the observations. It is anticipated that this warming may have resulted in short periods of diurnal surface stratification under low wind stress conditions; therefore, observations under these conditions were excluded from the analysis.
Wind speeds (at 3 m) had a range from 1 to 19 m s 21 with an rms of 9.2 m s 21 and maximum gusts of 28 m s 21 . Significant wave height varied between 1.2 and 14.1 m with an rms value of 5.3 m, while the average wave period varied between 4.4 and 14.4 s, with an rms of 8.0 s. The resulting surface wind stress t s varied between 2 3 10 24 and 1.2 Pa, with an rms of 0.27 Pa.
The « estimates were sorted according to the forcing conditions at the time of the observation, without any reference to adjustment time scales, resulting in the following datasets: The number of observations varied slightly between instruments and between methods, with the modified method having the same or fewer « estimates for each instrument. Observations made under conditions when t s # 0:05 Pa and B 0 # 0 (i.e., low wind and surface heating, respectively) composed 4.7% of observations and were excluded from the current analysis. The t s threshold was chosen based on the overall distributions of t s and B 0 , without any structured attempt at optimization.

a. Observation of wave orbital motion
Periodic variations were clearly apparent in much of the along-beam velocity data from each of the ADCP and were coherent across all bins in a beam. Fourier analysis typically showed a peak at or around the average surface wave period. To test whether the observations demonstrated the vertical gradient expected of wave orbital motion, the ADCP data were transformed from beam to Earth coordinates and the rms of Earth coordinate vertical velocity w rms and the difference dw rms over a vertical separation distance dz of 2.0 m was calculated for each instrument and for each 5-min sampling period. The theoretical variation in the wave orbital speed dy max was calculated over dz at each instrument's observation depth using (7), assuming monochromatic waves of amplitude equal to half of the concurrent significant wave height and with the observed average period. Figure 4 plots dw rms versus dy max together with the linear regression for each instrument. Despite the simplistic assumption of monochromatic waves in the calculation of dy max , all three instruments demonstrate a linear relationship with nearly identical coefficients over the full range of conditions. The robust correlation between dw rms and dy max , which are derived from independent datasets, indicates that wave orbital motions are producing a vertical gradient in the velocity profiles measured by the ADCP in a manner consistent with the simple theoretical model assumed.
b. Comparison of the standard and modified methods Figure 5 summarizes the results for the standard and modified methods for all three instruments and under both surface wind stress and convective forcing. All regressions are based on r max ; 2:0 m, the exact value depending on the separation distances evaluated given the ADCP geometry. The results for the two forcing processes are considered separately: 1) Wind stress forcing. The median wind stress scaled « estimates for each instrument and for both the standard and modified methods are shown in Fig. 5a, and the data are summarized in Table 1. For the standard method, the median scaled « estimates vary from 9.15 for the uppermost instrument to 1.78 for the lowest instrument, with a clear depth dependence. Over 45% of standard method « estimates at 24 m have a bias of an order of magnitude or greater compared with the default unity scaling, with .97% of observations exhibiting a bias of 2 or more. The bias decreases with depth, although over 45% of the observations at 52.5 m remain subject to a bias of 2 or more. In contrast, for the modified method, the median scaled « estimates vary between 1.11 and 0.69 for the three instruments, with no apparent depth dependence, suggesting no significant departure from the law of the wall unity scaling. 2) Convective forcing. The median surface buoyancyflux-scaled « estimates for each instrument and for both the standard and modified methods are shown in Fig. 5b, and the data are summarized in Table 2. The standard method median bias is higher for all instruments than the equivalent bias for the surfaceshear-stress-scaled observations, varying from 21.15 for the uppermost instrument to 2.21 for the lowest instrument and again demonstrating a clear depth dependence. In contrast, for the modified method, the median scaled « estimates vary between 1.36 and 0.79 for the three instruments and again exhibit no apparent depth dependency, suggesting no significant departure from the unity scaling with B 0 .

c. Method sensitivity to selection of r max
In principle, it is desirable to evaluate the structure function regression over as much of the inertial subrange as possible in order to better determine «, subject to the constraint on r max being less than the upper limit of the inertial subrange.
The sensitivity of the standard and modified methods to the choice of r max is illustrated in Fig. 6 for both wind stress and convective forcing with r max as close as possible to 1, 2, and 3 m. All of these r max values are expected to be within the inertial subrange given the water column density structure and turbulence levels. For r max ; 1 m, the regression of D LL against r 2/3 uses data for just eight separation distances (from two bins to nine bins). The number of separation distances increases approximately linearly with r max , subject to the dependence of the along-beam bin center spacing on the beam angle. For r max ; 2 m (3 m), the regression uses data for 18/16 (27/25) separation distances for the 208/308 instrument beam angles.
For the standard method, reducing r max reduces the bias but does not eliminate it. Even with r max reduced to 1 m, the median bias for observations at 24 m remains 4.2 for wind stress forcing and 8.2 for convective forcing. However, reducing r max to 1 m does reduce the median bias to less than two for the observations at 42.5 and 52.5 m for both forcing regimes.  The impact of reducing r max on the quality of the fit for the regression of D LL against r 2/3 and therefore on the confidence in the calculated « estimate is shown in Table 3 for wind stress forcing and in Table 4 for convective forcing. Reducing r max from ;2 to ;1 m dramatically reduces the mean R 2 adj values. For the modified method, varying r max has only a minimal impact on the median scaled « estimates for all three depths and for both forcing regimes. The difference in the median scaled « values is negligible for r max ; 1 m and 2 m, with the values for r max ; 3 m being fractionally lower. The R 2 adj values for the modified method consistently indicate a better fit than the standard method, although the difference is negligible for r max ; 1 m-only becoming significant with increasing r max .

d. Wave information from the modified method
The additional regression coefficient produced by the modified method (A 3 ) is expected to be dependent on the vertical difference in the speed of the wave orbital motion over the distance r max at the observation depth of the ADCP. Figure 7 plots the A 3 coefficient for each regression for each instrument against the square of the difference in the theoretical wave orbital speed based on the concurrent surface wave observations (dy max ), as described in section 4a, as well as the linear regressions for each instrument.
The scatter in Fig. 7 is considered to result from the assumption of monochromatic waves, with the average period of the surface waves not being fully representative of the spectrum of waves contributing to the vertical gradient in the wave orbital speed at the ADCP depths. However, despite this simplification, the clear linear relationship between the A 3 coefficient and (dy max ) 2 suggests that the modified method is extracting the contribution to the structure function as a result of the vertical variation in the wave orbital velocity speed as expected.  A specific dy max cannot be attributed to a unique surface wave condition, even under the assumption of monochromatic waves, since waves with different amplitudes and wavelengths could produce the same vertical velocity difference. In principle, it may be possible to use the variation of A 3 with depth to determine an ''effective'' surface monochromatic wave, but this is beyond the scope of the current study.

Discussion
While three decades of ocean turbulence measurements using ship-based microstructure profilers have provided strong quantitative links between the dissipation of turbulence kinetic energy and its forcing, the full geographic and temporal variability of turbulence, and hence mixing, remains a first-order problem in oceanographic research (Ivey et al. 2008;Moum and Rippeth 2009;Mead Silvester et al. 2014). Part of the solution to this problem has been the development of new techniques for measuring longer time series of turbulence parameters. Among the more successful techniques has been the application of moored offthe-shelf ADCPs, initially through the development of the variance method (Stacey et al. 1999;Lu and Lueck 1999;Rippeth et al. 2002), but more recently through a structure function approach (Wiles et al. 2006;Lucas et al. 2014).
In particular, the structure function technique is an attractive option, as the turbulence estimates are not sensitive to instrument motion and therefore can be made midwater column from moored platforms (Lucas et al. 2014), avoiding the specific processing to remove platform motion required for spectral techniques (Bluteau et al. 2016). Furthermore, the development of pulse-to-pulse coherent operating modes has enabled reliable estimates of « down to a noise floor estimated as ;3 3 10 210 W kg 21 (Lucas et al. 2014). However, the averaging period implicit in the structure function technique is long relative to the period of surface waves, potentially leading to a bias in « estimates as a result of the variation of the speed of the wave orbital motion with depth.
Here we have demonstrated the degree to which « is biased by the presence of surface waves using synthetic wave data. We have then developed a modified secondorder structure function method that exploits the differing length scale dependencies of the contributions resulting from turbulent and wave orbital motions in order to remove the surface wave influence. The standard and modified methods were then tested using data collected over a 3-month winter period by three ADCPs operating in pulseto-pulse coherent mode and mounted on a mooring at different depths. The observational period provided a wide range of wind, wave, and surface buoyancy flux conditions.
Estimates of « made using both the standard and modified structure function methods were then scaled using established scaling for either wind stress or convective forcing. The results using the standard method show a significant departure from the expected value under both forcing conditions. The bias is greatest for the uppermost instrument and declines significantly with depth. This accords with the hypothesis that the bias results from the vertical gradient in the speed of the wave orbital motions, which decay exponentially with depth. The median bias for convective forcing scaled « estimates were higher than those scaled for wind stress forcing at all depths, indicating that the bias from surface waves is more significant under relatively lower turbulence conditions. In contrast, the scaled « estimates obtained using the modified method collapse to approximately unity for the observations under both wind stress and convective forcing, indicating that the « profiles are in approximate accordance with the nominal scaling.
Analysis of the length scale dependence of the speed of wave orbital motions for intermediate depth waves (see the appendix) suggests that the modified method should also be effective in removing bias in « estimates from observations affected by surface waves in shallower water, providing the orbital motions match standard wave theory. However, pending evaluation against actual observations, care is needed in applying the modified method in shallow-water conditions.
These results lead to the following conclusions: d There is significant potential for bias in second-order structure function estimates of « as a result of the depth variation in the surface wave orbital velocities. d A modified method that exploits the differing length scale dependencies of the contributions to the structure function from turbulent and wave orbital motions is effective in removing the surface wave bias in the « estimates made under both wind stress and convective forced conditions. and Defra as a part of the Shelf Sea Biogeochemistry strategic research program, through Grants NE/K001760/1 and NE/K001701/1 (WP1 Candyfloss). Brian Scannell is supported by NERC Award 1500369 through the Envision doctoral training program. In addition J. Polton was supported by NERC standard grant Pycnocline Mixing in Shelf Seas (NE/L003325/1). We thank the crew of the RRS Discovery and the National Marine Facilities for their assistance in collecting the datasets. Thanks also to Jon Turton at the Met Office for supplying the ODAS buoy data, to Juliane Wihsgott for the CTD data from the moored temperature chain, and to Tom Hull for the Cefas SmartBuoy data. We also thank two anonymous reviewers for their comments, which helped to improve the manuscript.

Application with Generalized Wave Equations
The generalized equations for the motion under surface waves describe elliptical orbits with an eccentricity that depends on the wave's wavelength, the water depth, and the depth of the observation point. The horizontal and vertical velocity components under an infinitesimal monochromatic sinusoidal wave traveling in the x direction are given by where g is acceleration as a result of gravity; k is wavenumber given by k 5 2p/l and l is wavelength; v is radian frequency given by v 5 2p/T and T is wave period; z is depth, with z 5 0 at the sea surface and positive upward; h is water depth so that z 5 2h at the seabed; and t is time (Phillips 1977). A vertically oriented ADCP with a beam in the y 5 0 plane will see an along-beam velocity b 0 in the bin centered at x 5 x 0 and z 5 z 0 with contributions from both components depending on the beam angle u, which is given by The velocity difference db 0 over a vertical range dz around depth z 0 will therefore be FIG. 7. Modified method A 3 regression coefficient vs dy max for the three instruments with observations centered at depths 24.0 m (red), 42.5 m (orange), and 52.5 m (purple). Differences calculated over range dz 5 2:0 m; dw rms from Earth coordinate transformed velocities with rms over 300 profiles per 5-min sampling period; dy max based on monochromatic waves of amplitude half the observed significant wave height and with the observed average period.
OCTOBER 2017 S C A N N E L L E T A L .
where x z 0 2 dz 2 is the x coordinate of the observation bin centered at z5z 0 2 dz 2 . For u values of 208 or 308 and dz appropriate for r max values used with the structure function regression, the horizontal bin displacement x z 0 1 dz 2 2 x z 0 2 dz 2 will be l, so kx z 0 1 dz 2 ' kx z 0 2 dz 2 ' kx 0 and the orbital velocity observed in all bins is in phase. Equation (A3) then simplifies as Applying the double-angle hyperbolic identities and recognizing that cosh (sinh) is an even (odd) function, Grouping all the terms independent of dz into a function F, For kdz 1, the approximation sinh(x) ' x can be applied, giving For deep-water waves sinh[k(z 0 1 h)] ' cosh[k(z 0 1 h)] 'cosh(kh), so (A6) and (A2) become identical and (A7) becomes db 0 ' kb 0 dz, recovering (7).
More generally, (A8) suggests that while F may vary with z, db 0 will vary linearly with dz irrespective of the water depth, providing the wave orbital motion is described by (A1), subject only to the constraint of dz being small relative to l. This suggests that the modified method has the potential to be effective at removing bias resulting from wave orbital motion from « estimates over a wider range of water depths.

Testing the modified method for non-deep-water waves
It is reasonable to anticipate that there will be limits on the effectiveness of the modified method as the water depth reduces. To test this, synthetic velocity data were generated for waves with a range of wavelengths and amplitudes in different water depths in the same manner as described in section 2b but using the general wave orbital motion equations [(A1)] rather than deep-water equations [(8)].
Along-beam velocity data were calculated for a single upward-looking ADCP at a depth of 20 m, with 30 bins, with the first bin centered at 0.97 m from the transducer and with 0.1-m vertical bin center spacing. Velocities were calculated at 1-s intervals for a 5-min observation period. Surface wave wavelengths varied between 50 and 300 m and amplitudes varied between 0 and 2 m. The radian frequency was calculated from the dispersion relation c 2 5 (g/k)tanh(kh), where c is the wave phase speed. Wave orbital velocities resolved at 1-s intervals for 300 s. A background « level is imposed, varying with the surface wave amplitude from 1 3 10 210 W kg 21 for amplitude 0-m waves to 1 3 10 29 W kg 21 for amplitude 2-m waves, such that in the absence of any wave-related bias, contours 29.1, 29.2 . . .29.9 would be equally spaced horizontal lines.
OCTOBER 2017 S C A N N E L L E T A L .
The along-beam velocity data were processed to calculate the second-order structure function for separation distances up to the specified r max using a central difference scheme. A background « level was then added to the structure function so that the effectiveness of the modified method in recovering turbulence levels in the presence of wave orbital motions could be assessed. The imposed background « level varied logarithmically with wave amplitude from 1 3 10 210 to 1 3 10 29 W Kg 21 . The standard and modified methods were then used to calculate « estimates for each bin based on r max values between 1.0 and 3.0 m. An average « estimate was calculated as the geometric mean of the individual values for all bins for which the structure function was resolved for all r # r max . Figure Fig. A1 equates to a wider wave period range of 5.7-13.9 s. The figure shows that for the standard method, the bias introduced by the vertical gradient in the wave orbital speed overwhelms the imposed background «, with the level of bias for a given wavelength and amplitude increasing slightly in shallower water depths.
The results from the modified method demonstrate that the method is generally effective in recovering the imposed background « levels, with the effectiveness increasing with increasing wavelength. Reducing the water depth has only a minimal impact, with a slight improvement in effectiveness as the depth is reduced.
For the shortest wavelengths and the largest wave amplitudes, the modified method exhibits a negative bias, resulting in calculated « estimates lower than the imposed background values. This is due to the structure function regression against r 2/3 failing to separate the linear term used to calculate « from the (r 2/3 ) 3 term associated with the wave orbital motion. Increasing the imposed background level or increasing the depth of the observations reduces the effect, while increasing r max increases the effect. This effectively introduces an observation-depth-dependent limit on the method sensitivity in the presence of highfrequency waves.
The results from the tests with synthetic data demonstrate that providing the wave-induced orbital motion conforms to the standard equations, reducing the overall water depth does not significantly compromise the effectiveness of the modified method in removing bias in « estimates due to the presence of surface waves.