A wavelet-entropy based segmentation of turbulence measurements from a moored shear probe near the wavy sea surface

In this study, we explore the applicability of a wavelet-entropy based segmentation technique in reduction of motion-induced contaminations in time-domain from subsurface turbulence measurements made by a moving shear probe. After the quality screening of data, the Shannon entropy procedure is combined with a time-dependent adaptive wavelet thresholding method to split each 60-s long shear segment into a number of motion-reduced subblocks. The wavelet-entropy strategy leads to preventing the false detection effect caused by applying either wavelet de-noising or Shannon entropy alone for conditions where the turbulence (strongly) overlap with scales induced by waves or platform motions. The longest stationary subblock, with a size greater than 16-s, is then used to extract the Turbulent Kinetic Energy (TKE) dissipation rate, ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document}. Efficiency of the proposed method is verified by comparing with ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document} measurements made by a nearby free-falling microstructure profiler. While the quality of observations is constrained by a number of factors such as sensors’ angle of attack, and the wave kinematical and dynamical effects, results demonstrate significant improvements, by approximately a factor of 5–10, compared with ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document} measurements from each 60-s segment using the Goodman et al. [13] method. Furthermore, the magnitudes of the motion-corrected ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document} using the proposed method is largely consistent with the scaling suggested by Terray et al. [30].


Introduction
The study of oceanic turbulence near the sea surface is of great importance in order to understand the exchange processes of heat, momentum, energy, mass, and gas transfer rate taking place across the air-sea interface. Stewart and Grant [29] were the first who successfully measured the turbulent fluctuations in a tidal channel based on the inertial and inertial-convective subranges. Several towed, moored, and profiling instruments thereafter developed to efficiently estimate important turbulent quantities such as dissipation rates of Turbulent Kinetic Energy (TKE), , temperature variance, and eddy diffusivity throughout the water column [8,14,18,25,28,32]. Despite extensive advances in measuring technology over the last 40 years, accurate and reliable estimates of turbulent quantities, e.g.
, near the wavy sea surface is still a challenging problem due to the diverse sources of disturbance from the wave orbital velocities (in particular for the weak mean current), the wave-induced platform instabilities (large Angle-of-Attack, AOA), sensor limitations to operate appropriately in such energetic environment, and the lack of adequate statistically reliable signal processing tools to separate wave-related contaminations from the turbulent fluctuating motions [10]. On the other hand, besides large aliasing errors pertinent to the wave-induced sensor's instability in measurements from the moored floating systems (i.e. deterministic wave effect) and problems related to the mooring-line response [26], waves contribute in the production of boundary layer turbulence by breaking, with effects confined to the top few meters decaying very rapidly away from the surface, and Stokes drift of (irrotational) waves which tilts vertical vorticity into the downwave direction, forming counter-rotating vortices called Langmuir circulations [7]. Alternatively, waves can directly inject energy to the upper ocean boundary layer when they are not truly irrotational so that their orbital motions can induce turbulence (i.e. stochastic wave effect).
To measure directly TKE dissipation, shear probes can be used as a more precise complement (due to their high sampling frequency with low noise level) of the (indirect) acoustic-based measuring techniques such as Acoustic Doppler Velocimeter (ADV) and Acoustic Doppler Current Profiler (ADCP). The accuracy of turbulence measurements from shear signals depends, however, on the portion of resolved shear features and scales within the inertial subrange, in particular at low wavenumbers. Prior to obtain values, the quality of the collected shear time series needs to be exactly verified in terms of stationarity of the signal, small AOA of sensor's tip relative to the mean flow and the applicability of Taylor's frozen turbulence hypothesis. In particular for the unattended moored platforms, these constraints may significantly confine plausible quality to a small subset of environmental conditions [10,22]. For the (statistically) reliable estimates of , the qualified data segments should also be long enough to contain a sufficient range of turbulent scales. The corresponding shear spectra after conversion from frequency to wavenumber domain are then corrected for both the spatial averaging due to the physical size of the airfoil probe [20] and for the vibration of sensor using the method outlined in Goodman et al. [13]. However, the shape of corrected shear spectra still exhibit direct and indirect (wave-current-turbulence interaction) contributions from the surface gravity waves which may greatly influence the quality of TKE dissipation calculations in the near surface marine boundary layer.
In this study, we use shear datasets collected from a Moored Autonomous Turbulence System (MATS) deployed in the close vicinity of air-sea interface, fixed at an approximately 6.5 m below the sea surface (Sect. 2). Near the wavy interface, the recorded shear data reveal strong anisotropy (i.e. departure from isotropy in both large and small scales of motions) and statistically non-stationary wave-and motion-induced (organised) structures. These features may not be appropriately detected/removed by applying some known procedures to separate the turbulent fluctuations and sensor motions (vibrations), for example the multivariate technique of Goodman et al. [13] which reduces the effect of platform vibration. In such cases, thresholding of the time sequence is a reliable technique to separate the background turbulence information from the wave-and the motion-induced disturbances. However, the application of the global thresholding strategy is limited for our shear datasets since the probability of finding long-duration segments ( ∼ 16-s) using the global thresholding scheme significantly decreases due to timevarying and non-uniform nature of burst-like (deterministic wave/motion related) disturbances. Therefore, we propose a segmentation algorithm in time-domain based on an adaptive wavelet thresholding method coupled with a Shannon energy technique (i.e. AWST) to split each shear time series into a finite number of (quasi) stationary and motion-decontaminated subblocks (Sect. 3). The longest stationary subblock with a size greater than 16-s is then used to estimate . The results from MATS's shear probe (only for the shear probe 1, hereafter SH1) are further compared with measurements from the nearby free-falling microstructure profiler at depths close to the depth of MATS, predictions via the "law of the wall" scaling, and scaling suggested by Terray et al. [30], Sect. 4.

Data collection and forcing conditions
Observations of upper ocean microstructure were carried out during a cruise of the Research Vessel Håkon Mosby, between 28 and 30 November 2012. The measurement site was approximately 30 km southwest of Bergen, Norway in ∼ 20 m water depth (Fig. 1a). A loosely tethered free-falling profiler MSS-90L (ISW Wassermesstechnik, Germany, MSS hereafter) and the MATS were used to sample ocean turbulence using two shear probes and two fast-response temperature sensors. MATS was set to acquire data at approximately 6.5 m below the sea surface for this deployment at sampling frequency of 512 Hz for fast channels (i.e. 3D accelerometers, two airfoil shear channels, and two fast-response temperature sensors) with 15-min bursts followed by 1-min of file book-keeping (Fig. 1b). Each 15-min burst was then segmented into half-overlapping 60-s sub time sequences, which are used for our analysis to extract . A Nortek ADV with sampling frequency of 16 Hz mounted horizontally from the nose of MATS's buoy provided 3D velocity measurements used to convert the shear frequency spectrum to the wavenumber spectrum using the Taylor's frozen turbulence hypothesis. Wave bulk parameters were measured by a bottom-mounted 1 MHz Nortek AWAC which recorded, at 1 Hz for 1024-s at each hour, both currents and the distance from the sea surface via Acoustic Surface Tracking (AST).
The dominant wind speeds observed from ship's meteorological mast were approximately between 2 and 10 m s −1 with almost steady direction from southwest (Fig. 2a). Few hours after beginning of experiment, the wind increased to > 8 m s −1 until almost 17 hours where  Time series of a wind speed (black line) and wind direction (dashed red line) measured from the ships's meteorological mast at 15 m height and scaled to 10 m height using the COARE 3.0 algorithm at neutral stability condition [9]; b wave bulk parameters measured from the AWAC deployed in the close vicinity of MATS's location at a water depth of approximately 20 m; and c wave age, A w = c p ∕u * a , where c p is the peak phase speed calculated from the dispersion relation for the shallow water and u * a is the friction velocity at air the wind speed dropped to < 6 m s −1 accompanied by an increase in the wave peak period. Significant wave height, H s , varied between 1.3 and 1.8 m, and the mean wave period of approximately 8-s was on average the dominant wave period during the experiment. The wave age at early in the deployment shows developed seas ( A w ≥ 80 ) and decreases to about 40 when wind speed drops at day 1.2 ( Fig. 2c).

Methodology
Estimation of turbulence using time series of shear data acquired from a moored subsurface buoy near the wavy surface is not a trivial task because of the presence of strong non-turbulent structures resulting from platform's vibration, variable resonance response of the subsurface buoy to the environmental forcing, mooring-line tension, and wave-induced contaminations. While using appropriate swivels will reduce both the tension from mooring-line and the resonance response of the buoy, the (deterministic) motion-induced vibrations may be effectively reduced by using the multi-variate spectral coherence technique of Goodman et al. [13], hereinafter G06. Platform motion can further change the flow velocities near the measuring probes, i.e. time-varying flow distortion [24]. Surface gravity waves contribute in production of both turbulent (stochastic) and non-turbulent (deterministic) scales near the sea surface, and are also able to induce motions to the measuring sensors. Furthermore, the buoy motion induces scales interacting strongly with wave orbital velocities and may cause further challenges in the separation of the two. Although the velocities induced by the movement of sensors can be properly removed/reduced using the data from the inertial motion package affixed to the MATS [11], the application of the method to our data cannot be addressed with confidence. This is because the power spectra of ADV data contaminated significantly by both the wave-induced disturbances at low-frequency components of motions and noise at the high-frequency portion of vast majority of spectra. Therefore, we could not detect credible near-inertial subrange above the noise level for the large number of ADV measurements in order to estimate the dissipation rate of TKE, Fig. 3.
In this study, the shear signal segmentation aims at segmenting the signal into a series of motion-reduced subblocks as a first step in estimation of TKE dissipation rates from (non-stationary) shear data. In order to identify the non-turbulent wave-, vibration-, or motion-induced structures from the background signal, a four-stage procedure is designed as: (i) the choice of a near optimal wavelet family and quality check of shear data; (ii) an adaptive wavelet thresholding module to detect wavelike motions of time-varying amplitudes and periods from the quality controlled shear signal; (iii) a Shannon entropy to split the non-turbulent burst-like features from the turbulent structures and a non-parametric stationarity check of the segmented blocks which their lengths are greater than 16-s; and (iv) estimation of dissipation rates of TKE from the largest stationary segment of shear probe time series, see Fig. 4. A search algorithm for the selection of the near optimum wavelet family is given in Appendix 1, and Appendix 2 briefly explains the non-parametric stationarity test used in this study.

Shear spectra and wave advection
The voltage output of shear channels from the MicroRider (MR) are converted to shear signals in physical units using some known electronic constants and the speed of flow passed the sensors from the Vector data, read more details in Wolk et al. [32] and Bakhoday-Paskyabi and Fer [4]. According to the latter in order to calculate the rate of TKE dissipation, the ADV time series are used to convert the shear spectra from the frequency domain (here w (f ) and v (f ) ) to the wavenumber space through Taylor's hypothesis of frozen turbulence: Typical frequency spectra of u, v, and w components of velocity measured from ADV; u , v , and w respectively. Also shown contains wave-affected subrange and −5∕3 spectral slope for identifying inertial subrange. The vertical dashed lines show the lower f l p = 0.5f p , and the upper f u p = f p + 0.5 bounds for the waveaffected subrange (coloured area), here f p = 1∕T p where T p denotes the wave peak period where i = 2 and 3 stand for w and v components of shear signals, U d is the mean flow magnitude (in rms sense) past the sensor averaged over the 60-s duration of the segment estimated from ADV data, f is the frequency, and k = f ∕U d denotes the wavenumber. In an oscillatory boundary layer, it is, however, not clear how to directly use Eq. (1) due to the advection of turbulence eddies by a steady background mean drift and an unsteady waveinduced motion, as well as the difficulty in separation of wave motions from turbulence scales over a wide range of frequencies/wavenumbers. According to Lumley and Terray [19], in the presence of waves, majority of the bulk of turbulence energy is shifted to frequencies higher than the wave dominant frequency, f p , compared with the same turbulence advected only by the background steady current. The broadening of peaks in the turbulence spectrum is controlled primarily by the harmonics of the wave dominant peak frequency so that the oscillatory characteristics appeared in the turbulence spectrum show dependency to the values of where ̃ denotes the magnitude of the wave-induced standard deviation. The smaller values of R, the more significant effects are expected from the wave advection on the shape of spectrum close to the f p , and at high frequencies much above f p and roll-off frequency of energy-containing eddies when advected by the ambient current alone. Therefore, the rates of dissipation extracted from the −5∕3 fits to the velocity spectra within the inertial subrange are overestimated if the wave effects are not taken into account. The elevated energy levels appeared in the shape of velocity spectra for small values of R (i.e. ≪ 1 ) can be described as a correction factor I for the −5∕3 region of the velocity spectra not affected by waves as where k = (k 1 , k 2 , k 3 ) is the wavenumber vector with magnitude of k, lm indicates the Kronecker delta function, and ū i and ̃2 i are the mean current velocities and the waveinduced variances, respectively, for i = 1, 2, 3 by assuming ū 3 = 0 . If we assume spectral shape at all frequencies can be represented by f −5∕3 , a simplified expression for I by setting ̃3 to zero is obtained as [31]: where is a non-dimensional frequency, and is the angle between wave field and current. In this paper, we reduce the effects of wave advection to satisfy Taylor's hypothesis of frozen turbulence by limiting the quality data for conditions where R < 1.15 , and only use expressions (3) and (4) to illustrate the effects of wave advection on the shape of shear spectrum when converted to the velocity spectrum (Appendix 4, and inset in Fig. 10b). Correction of shear data to account for the wave advection and its effect on the estimation of are addressed in more details in another investigation.

Quality control criteria
Time series of shear probes collected from the MATS during this experiment show substantial anomalies in response to the background current, surface waves, bodyinduced motions, and flow distortion, Fig. 5a, b. These disturbances are often time-varying, non-Gaussian, and correlated with suitable information (i.e. turbulent motions) in the signal. Some of them such as the sensor's vibration can be then removed/reduced using information from MATS's acceleration datasets (Fig. 5c, d) recorded by the (4) vibration sensors using the Goodman et al. [13] method at an appropriate frequency band. Additionally, shear probes to measure rely on the potential flow assumption in order to estimate the hydrostatic turbulence-induced lift force. To not violate this principle assumption, we extend the quality criteria by calculating the instantaneous AOA=tan −1 w∕u h from the MATS's ADV velocity measurements, here w and u h are the vertical velocity and the segment-averaged along-axis velocity, respectively. Note that the observed values for the AOA are not representative for the true values of the AOA since the velocity measurements are by themselves subject to the motion-induced biases. In summary, we label each segment as good data satisfying any of conditions including AOA HF ≤ 7 o , the de-trended roll angle (in rms-sense) less than 1 o , R > 1.15 , and mean current speed greater than 0.05 m s −1 (see also [2,10]).

Adaptive wavelet-based de-noising
Wavelet Transform (WT) as a mathematical microscope is a powerful alternative to the Fourier method to analyse the non-stationary signals. In contrast to sines and cosines as basis functions for the classical Fourier series, wavelets are based on translations and dilations of a single fixed function , so-called the mother wavelet, as where * represents the complex conjugate, s(t) is a squareintegrable function at time t, a and b are dilation and translation parameters, respectively [21]. Equation (5) suggests that wavelets are local in both time (via translations) and frequency (via dialations). For most practical purposes, it is sufficient to use Discrete form of WT (i.e. DWT) which for a certain choice of (here Daubechies wavelets) based on the multi-resolution analysis decomposes the signal into approximation and detail subsets at multiscale levels: where s is a time series, here shear signal, N denotes the number of sampling points, a J (k) are the approximation coefficients at decomposition level J, and d j (k) are the detail coefficients at jth level, Fig. 6a.
In the development of wavelet application for the shear data collected underneath the wavy air-sea interface, we gain the fact that the explosive peaks and noisy structures in the measured (shear) signal die out swiftly by increasing wavelet scales, while the information related to the organised wave-and motion-induced structures may have significant components over multiple scales [15,16], see Fig. 7. This is the primary motivation to develop an adaptive wavelet-based thresholding algorithm in time-domain to separate the turbulence from the motion-induced structures.
The shear spectra of the SH1 corresponding to the low and high wind conditions show low frequency peaks below approximately 2 cycle per meter (cpm) correspond to the gravity wave motions (e.g. Fig. 10). In the absence of superposition of turbulence and different modes of Fig. 6 a A discrete wavelet decomposition tree with depth of 2 including two leaves at each level, i.e. detail coefficients, d 1 , d 2 , and approximation coefficients, a 1 , a 2 . h and g are low-pass and high-pass wavelet filters, respectively. The procedure is performed by applying a downsampling of 2; and b the flowchart of the adaptive wavelet thresholding algorithm at ith iteration. Here, T j,i is the optimum threshold at wavelet level j and iteration i, and m j, are non-characteristic and characteristic coefficients (in the wavelet space), respectively. f i and m i are the reconstructed time series (using IWT) of the characteristic and non-characteristic coefficients in the physical space, respectively wavelike motions above 2 cpm, we are led to an unfair assumption that there might be a spectral gap (or a phase relationship) between waves and turbulence, in particular for the wind-generated waves, for wavenumbers ≫ 2 cpm. While the spectral gap between waves and turbulence is not practically possible and there is clearly association of waves with turbulence of different scales, we conceptually assume that the observed shear signal s(t) can be decomposed as where f(t) represents the non-turbulent wave/motion induced component, and m(t) represents a mixture of uncorrelated noise and turbulent fluctuating motions. . It is clear that motion/wave-like features have significant contributions over multiple scales. For the sake of simplicity, we plot the coefficients at few levels For the adaptive sub-level thresholding method in this paper, we utilise the 16th order Daubechies wavelet with J = 8 levels of decomposition (see Appendix 1). The iterative process is outlined in Algorthim 1 in which the shear probe data-set is first subject to quality check and the main loop is then started by decomposing the signal using DWT into wavelet (approximation and detail) coefficients, see Eq. 6. By assuming f = 0, the wavelet coefficients at scale j and iteration i are passed through the following threshold criterion: to construct f j i (t) (characteristic coefficients with values larger than the threshold) and m j i (t) (non-characteristic coefficients with values smaller than the threshold) at each wavelet level j ≤ J . Here, and are the mean value and the standard deviation of the coefficients. is an adjustment factor which for our data varies between 2 and 4 (here β is set to 4 based on trial and error) to keep T large enough at different scales, see Fig. 6.
After performing the thresholding operation over all j, the Inverse Wavelet Transform (IWT) is used to reconstruct the characteristic coefficients, m i , and the where * = 0.05 is a recommended value. After termination of the iterative process, the reconstructed characteristic time-sequence, f, is first normalised as s norm [n] = f [n]∕ max(f ) and is then passed to the segmentation process.

A time-dependent Shannon Entropy Segmentation (SES)
To characterise the randomness and chaotic nature of the filtered signal, we use the average Shannon energy by the use of sliding windows associated with a random variable as follows where N denotes the size of sliding window and s j [n] is the normalised de-noised (characteristic) signal, i.e. s norm , at jth window, see Fig. 8. We use windows with length of 4 × t with 50% overlap, where t denotes the sampling interval. The energy representation in Eq. (10), which extracts the envelope of s(t), attenuates the low-intensity components of the signal and enhances the medium energy components, see Fig. 8. The non-turbulent (wave-and vibration-induced) patches are then detected by all the mean-removed, Ẽ s = (E s − E s ) , and normalised Shannon energy, Ẽ s ←Ẽ s ∕ max(Ẽ s ) , time sequences using the following thresholding role: Here overbar denotes a temporal mean and is a thresholding parameter which can be considered as a constant (zero/near-zero) value (we use = 0 ) or determined adaptively. A segment with maximum length (greater than 16-s long) satisfying the stationary criterion (Appendix 2) is then selected to estimate the dissipation rates of TKE by utilising the Goodman et al. [13] method to remove the effects of residual vibrations. Overall, the method can successfully separate the wavelike motions by thresholding the values of Ẽ s when turbulence scales interacting weakly with the wavelike structures/objects (i.e. quasi-deterministic wave effect). The selected segments may, however, exhibit contributions from the irrotational wave motions and complex wave-turbulence interaction (i.e. stochastic wave effect).
In Eq. (10), one can utilise time series of the original (normalised) signal without applying the WT (i.e. SES method). Applying both methods (i.e. SES and AWT) to each 60-s shear data leads to results with good quality for conditions where the selected segments are long enough (e.g. ≥ 16-s). For our datasets, the SES method produces largely shorter duration segments compared with those detected from the AWST scheme and reveals more sensitivity to the strength and nonlinearity of the motioninduced structures.

Dissipation rate
Under the local isotropy assumption, the stationary segment of shear signal after applying the AWST algorithm (i.e. s(t) = w∕ x for SH1) is used to calculate the dissipation rate of TKE by where x denotes the horizontal axis, w is vertical velocity component, w is the observed shear spectrum in the wavenumber domain computed using Welch's averaged periodogram [4]. is the kinematic viscosity as a function of the water temperature, k is the wavenumber, k min and k max = (2 ) −1 ( ∕ 3 ) 1∕4 are the bounds for the integration, and the overbar in Eq. (12) represents a time average operator. The lower integration limit is set to 7 cpm and the shear spectrum is corrected for the shear probe's spatial response limitation with a cut-off wavenumber of k c = 48 cpm (by the transfer function H(k) = 1∕(1 + (k∕k c ) 2 ) for Rockland's scientific shear probes). An initial estimate of is derived from the observed shear spectrum to estimate the upper limit of the integration by resolving the viscous subrange and minimising the contribution from the noise-dominated portion of the spectrum [10]. Due to the energy elevation near the surface, most of the variance in the viscous subrange is moved beyond the k = 150 cpm (towards higher wavenumbers) which is the maximum wavenumber for reliable measurement of using the integration technique. Therefore, the upper limit of integration is iteratively determined and unresolved wavenumbers are computed using the variance of the empirical Nasmyth spectrum [23]. For an in-depth explanation of extraction from Eq. (12), we refer reader to Fer and Bakhoday-Paskyabi [10] and Bluteau et al. [6].

Data and quality screening
After applying data screening criteria listed in Sect. 3.2, the total number of SH1 segments decreases from 2324 to 669 at which about 67% of the removed segments are because of R criterion alone. Figure 9a shows the high frequency component of AOA, , with the grey region indicating good quality data so that AOA HF ≤ 7 o , where w * and u * h are the band-passed vertical and horizontal velocity components for the frequencies between 0.5-2 Hz. Additionally, we discard all segments with substantial (rms) roll angle greater than 1 o (Fig. 9b) and small values of the ratio between the mean flow and the wave-induced flow, i.e. R < 1.15 in which values of ̃ are extracted from the integration of the ADV's velocity spectra within the wave-affected frequency band (Fig. 9c). The frequency of the wave peak, above 0.08 Hz, measured from the AWAC data is used to identify the left, f l p , and the right, f u p , bounds of the wave frequency range. Here, the lower frequency is specified as 50% of the peak wave frequency and the upper frequency bound is considered as f u p = f p + 0.5 Hz. Due to the uncertainty in calculation of R, resulting from the complicated interactions between the non-turbulent waves and ADV motion within the wave-affected frequency band (and much higher and lower than f p ), we recalculate R using a spectral model by fitting each observed spectrum to the Kaimal et al. [17] empirical model (Fig. 16). The method has been proven to be less influenced by the contaminations induced from the wave motions within the inertial subrange [12], see also Appendix 3. Figure 9c suggests that the vast majority of energy elevation within the wave-affected frequency band is associated with the wave orbital motions (red dotted line), and the integration-based method (blue line) is then a reliable approach to estimate ̃.

Shear spectra and the AWST method
The wavenumber shape of the ensemble averaged 60-s shear spectrum satisfying the quality criteria, w , is shown in Fig. 10a along with the spectra from the not calibrated vibration sensors. Accelerometer spectra show high levels of energy elevation at high and low wavenumbers corresponding to the sensor's vibration noise and contaminations induced by the surface waves, respectively. The coherent variances with accelerometer data have been further removed from the shear spectrum using the G06 method. The averaged spectrum can be then divided into four wavenumber (frequency) bands including: (1) wave-dominated/energy-containing subrange (red area); (2) transition subrange associated with the nonlinear interactions between waves and turbulence (with a noncascading behaviour); (3) a narrowband inertial subrange over a sufficient number of wavenumbers (green area) with a very little contribution from the dynamical waveinduced shear, but may still suffer from the kinematical wave effects; and (4) a subrange associated with the dissipative scales and sensor's noise. Shown in this figure also consists of a wavenumber spectrum of the vertical wave orbital accelerametion, w∕ t , at the depth of MATS (calculated by the linear wave theory from the high resolution measurements of pressure, see [4]). Strong variances in the wave band rapidly decay with wavenumber within the transition subrange. As energy levels in the wave-affected frequency band increase, the transition subrange is further expanded towards the higher wavenumbers, and strongly influences the energy broadening and distribution within the inertial subrange in almost kinematical sense. This wave effect may constrain the low wavenumber cutoff to a large value when is used to calculate from Eq. (12). On the basis of Fig. 10a, we illustrate the one-dimensional energy spectrum, w (k) = k −2 w (k) , to further investigate the effects of wave-induced motions on the measured shear spectrum. In the presence of waves, the energycontaining rolloff, expected to be shifted towards the lower wavenumbers when R decreases, is swamped by the wave-induced variances with no apparent flat plateau at very low wavenumbers. The vertical/horizontal excursion of wave orbits are much larger than the sizes of energycontaining eddies and turbulence energy in this subrange (and the transient subrange) is significantly redistributed and modulated by the wave-induced scales (wave bias). This makes the separation between turbulence and waves very difficult and results in an elevation of energy levels at wavenumbers much farther than wave peak wavenumber due to the advection of turbulence by waves (i.e. green area wavenumbers ≾ 18 cpm with a slope slightly deviated Fig. 10 a The ensemble averaged shear spectrum over all qualified 60-s segments (thick red curve) and averaged acceleration spectra from vibration sensors, VA y and VA z respectively. The grey curve indicates the Nasmyth's turbulence spectrum. The wavenumber band of 7-40 cpm is passed into a routine in order to iteratively adjust the integration band; and b the averaged power spectrum converted from the averaged shear spectrum, i.e. w (k) = k −2 w (k) when H s = 1.4 m, T p = 11.8 s, R = 1.75 , and U d = 0.32 m s −1 . The inset demonstrates the diagram of the wave advection modifier I as a function of R for all 60-s segments from the −5∕3 power law, blue line). The extent at which wave advection controls the level of energy elevation may be explained by the correction factor I(R, ) defined in Eq. (4). For very small values of R (Appendix 4), applying I causes efficient corrections in the shape of energy spectrum and in the value of extracted . A scatterplot of I against R for all data shows that 0.8 < R < 4 (inset in Fig. 10b) with < 90 o , for the vast majority of qualified segments. Therefore, in our analysis, the uncertainty due to the wave advection for qualified segments that pass the criterion R > 1.15 produces little error in estimation of . It is pointed out that a precise interpretation of the measured shear/energy spectra underneath gravity waves requires further in-depth analysis for taking into account different aspects of advection of turbulent eddies by directional (rotational/irrotational) waves.
Prior to the application of the vibration removal technique, i.e. G06, each 60-s long segment is encountered with the AWST algorithm based on the Daubechies wavelet of order 16 with decomposition level of J = 8 . The performance of this de-noising algorithm depends directly on the choice of the decomposition level and type of selected wavelet family (Appendix 1). While small J leads to an incomplete de-noising, the large value for J results in signal distortion. Figure 11a, black curve, depicts the filtered time series after applying the adaptive wavelet algorithm for a segment of data taken when H s = 1.4 m. The Shannon envelop of the characteristic coefficients in the physical space reflects then the amplitude deflections and durations of the non-turbulent disturbances for the 60-s segment, i.e. Fig. 11a (red). Moreover, we have applied a global thresholding strategy based on the Shannon entropy segmentation (SES) of the original signal, i.e. Eq. (10), in the absence of applying the adaptive wavelet thresholding procedure, Fig. 11b (red). Considering the strength of the non-turbulent events and the value of the global threshold (here = 0 ) to constrain time series of entropy metric Ẽ s , both the AWST and the SES methods can identify the burstlike disturbances with sufficient quality for each 60-s time series. However, using the SES-based global thresholding strategy generally decreases the probability of detecting long-duration segments for our shear datasets due to time-varying (in amplitude/period) and the non-uniform nature of disturbances, and typical disability of the nonadaptive thresholding schemes in capturing these sharpedged features (Fig. 11c). After detection and elimination of the motion-related data points, all selected segments run through the stationary check procedure using the Reverse Arrangement (RA) test (Appendix 2). The largest stationary and coherent time-block with a length greater than 16-s, Fig. 11c (red), is then selected to estimate . Among wide range of wave-related (deterministic/stochastic) motions, the AWST method can capture those scales rising from weak wave-turbulence interactions and from deterministic wave-induced platform motions. Figure 12 illustrates example shear time series and their corresponding spectra for two different bursts from the shear probe 1 (SH1) satisfying the quality screening Fig. 11 a Results of applying the adaptive wavelet thresholding method on a 60-s segment of SH1 for the sea-state condition when H s = 1.3 m. Black curve indicates the application of the DWT (i.e. Db16, see Table 2) for J = 8 and red curve represents the Shannon envelop detection method for the wavelet coefficients; b similar to (a) except using SES method instead of the AWST scheme to calculate time series of entropy metric Ẽ s ; c the longest detected segments from both the AWST method (red area) and the SES method (blue area).The black lines in this figure show the original 60-s long time-sequence criteria, i.e. AOA HF , roll angle, flow distortion effect, turbulence advection by waves, and the segment length limit [2]. In Fig. 12a, a 22-s time-sequence from the corresponding 60-s segment is selected by the AWST segmentation algorithm as the longest less motion-contaminated portion of the time series. Due to the high level of turbulence under wavy conditions, the minimum length of data-segment (i.e 16-s) is long enough for statically significant calculation of from Eq. (12) using appropriate theoretical integration limits. Figure 12b consists of the spectra from the (non-calibrated) vibration sensors with arbitrary spectral levels in order to only identify/ remove the shear variances that are coherent with vibrations (purple and green curves) using the Goodman et al. [13] approach. Both components of acceleration show elevated levels of vibration at high ( > 40 cpm) and low ( < 3 cpm) wavenumbers. In Fig. 12b, we also compare the wavenumber domain spectra of the entire 60-s time series cleaned with the G06 method (red) and the 22-s segment (yellow) selected by the AWST method. Note that the reduced variances resulting from the application of G06 method are partly due to the elevated levels of vibration at low wavenumbers. This low wavenumbers variance-removal using the G06 method does not cause any significant underestimation in the estimated values of dissipation rates compared with the true values since the wave-induced variances are the most dominant features in the wave-affected subrange. The AWST-based Fig. 12 a Results of applying the adaptive wavelet thresholding method for a 60-s segment from the SH1 when H s = 1.2 m. Blue curve indicates the original signal, and the red curve represents the longest detected subblock with the size of 17-s using the AWST algorithm; b the corresponding vibration corrected wavenumber spectrum using Goodman et al. [13] method (red line) for 60-s segment and the one based on the AWST without applying the Goodman et al. [13] method (yellow) for the 17-s time-sequence; c a 60-s segment from another burst of the SH1 (blue line) together with the AWST-selected 31-s block (red line); and d the corresponding wavenumber spectra for the 60-s segment shown in Fig. 12c. The grey curves illustrate the Nasmyth spectral shapes. For the AWST method, we apply the Daubechies wavelet with the order of 16 and J = 8 (Appendix 1). The red areas in this figure demonstrate the wave-affected band. Figure 12b, d consist of the spectra from the (non-calibrated) vibration sensors with arbitrary spectral levels. The arrows show the input integration bounds to be used in the iterative algorithm to correct for the unresolved portion of shear spectra spectrum removes elevated variances over wavenumbers before the spectrum starts to roll off, here < 30 cpm. Furthermore, both shear spectra roll off with the same steep compared with the Nasmyth's empirical spectrum (grey line), identifying the inertial subrange, and the viscous dip wavenumber is also identifiable from the spectra of the AWST-based time segments (the black arrow). Figure 5c and d show the same analysis carried out for another 15-min burst at which the wavenumber spectrum of a 31-s segment from the AWST method conforms better to the empirical spectrum (grey line) compared with the vibration decontaminated spectrum of the 60-s time series.

Turbulence dissipation rate
Measurements of turbulence profiles from the MSS profiler were carried out within 3 periods in November 2012: period 1 (days 1.33-1.82 since November 28); period 2 (days 2.25-2.26 since November 28); and period 3 (days 2.29-2.33 since November 28). Total number of 216 casts has been conducted down to the depth of 50 m, approximately 200-300 m away from the MATS location. The buoyancy of the profiler was carefully balanced in the water column to keep an approximately constant sinking velocity ranging from 0.1 to 0.2 m s −1 . Figure 13a, c shows the time-averaged wave energy frequency spectra at each period of MSS deployment. The significant wave height does not vary very much and the values for the wave peak period, ranging between 10-11 s, suggest a swelldominated sea-state during the profiling experiments. Due to strong turbulence generation near the drifting ship and wave-induced motions near the sea surface, the data below 9 m from the surface are considered to be reliable for calculations (Fig. 13d, e, and f ). Furthermore, we use time series of averaged over a vertical segment for depths between 9 and 11 m below the sea surface (green areas), to reduce the effects of wave-induced forcing, and other disturbances from surface in the extraction of dissipation rates for comparison with results from the MATS's The values for are calculated over segments of 1024 samples with a 50% overlap in order to preserve a good tradeoff between turbulence homogeneity and its statistical significance. The green areas demonstrate the vertical layers, at depths between 9-11 m below the sea surface, from which the extracted values of are compared with those from the MATS' SH1 measurements. However, this vertical distance between two measuring systems increases the uncertainty in the statistically significant comparisons of two measurements due to the depth-dependent variability of surface gravity waves (depth-attenuation of wave field), see also Fig. 3 in Bakhoday-Paskyabi and Fer [3]. The processing of the shear probe data from the MSS is similar to the procedure discussed for MR (Sect. 3.5) using the isotropic Eq. (12), but for the wavenumber spectra of the vertical shear using u∕ z and v∕ z . In Fig. 13 d-f, we determine profiles of measured by the MSS shear spectra from 50% overlapping 8-s segments, using an FFT with length of 1024 sampling points (i.e. 2-s). Moreover, the turbulence in shallow water during our profiling experiments did not reach anymore to the values of below the detection limit of MSS's shear probes, i.e in the order of ≈ 10 −10 m 2 s −3 .
The evolution of dissipation rates collected by MSS in the depth interval 50 to 0 m during the first profiling period (i.e. the first 190 casts) is qualitatively shown in Fig. 14a along with an example of a representative MSS cast (profile), dashed vertical line. We set manually the depth segment in this example in order to reduce the spatial uncertainty resulting from the vertical distance between two measuring systems (Fig. 14c). To compare the quality of the AWST method in reducing the  Fig. 14a. This plot contains also the detected segment (red line) by the AWST method; c vertical profile of MSS's SH1 (black line) and the manually selected depth bin (red line); d the spectrum of the representative MSS cast at depths between 6.7-11.3 m, and the cleaned MAST's SH1 spectrum using G06 technique. The dashed grey line indicates the Nasmyth spectrum for the value of estimated from the MSS's SH1 time series; and e the shear spectrum of the representative MSS cast, and the MAST's SH1 spectrum by applying the AWST method with no vibration removal motion-induced contaminations, a 60-s time series of the MATS's SH1 for almost the same time as the representative MSS cast is shown in Fig. 14b. The influences of motion and vibration are apparent in the form of several non-turbulent structures, and the AWST method shows a successful detection of a less motion-affected segment, red line. The shear spectra from the MATS's 60-s time series and the MSS are shown in Fig. 14d along with the Nasmyth spectrum for quality check of the observed spectra. Both spectra for wavenumbers between 7-30 cpm may be considered to be of acceptable quality for the isotropic turbulence, i.e. inertial subrange, but the discrepancies between two may be explained by the spatial separation between two systems, the intermittency of turbulence, and contaminations induced by waves and platform motions (for MATS measurements). When applying the AWT method alone (not AWT+SES), the discrepancy between two microstructure measurements substantially decreases (Fig. 14e) which again confirm the ability of the method in reducing motion-related (deterministic) disturbances at wavenumbers before the spectral roll-off at higher wavenumbers. Figure 15 shows the hourly averaged time series of u * w and H s (Fig. 15a), and dissipation rates (Fig. 15b) over the whole period of experiment for the qualified and stationary segments. Here, the number of data points at each hourly bin are averaged and the accuracy of the expected value of with 95% confidence limits is estimated by the maximum likelihood estimator from a log-normal distribution [1,4]. Also shown in Fig. 15b are the dissipation rates estimated from the available free-falling MSS profiler data averaged at the depths very close to the depth of MATS, i.e. MSS , as ground-truth estimates (circle markers). The agreement between MSS and the MATS' dissipation rates for 16-s segments, AWST MATS , suggests that the segmentation technique can efficiently reduce (remove) the motion-related contaminations in time-domain from the moored shear observations near the sea surface. It is also able to improve the accuracy of estimations by a factor of 5-10 compared with those calculated from the application of Goodman et al. [13] method to each 60-s time sequence, MATS . Comparison of MATS's dissipation measurements with those predicted by a wall-layer estimate at the depth of MATS: Fig. 15 Time series of hourly averaged: a water-side friction velocity (black line) and significant waveheight (red markers); and b dissipation rates of TKE from the cleaned 60-s long segments using Goodman et al. [13] method (black markers); from the longest subblock at each 60-s long segment extracted from the application of the AWST method along with the Goodman et al. [13] scheme to reduce the effects of residual sensor's vibration (blue markers); and dissipation rates measured from the MSS averaged at depths between 9-11 m below the sea surface (open circles). Grey bars denote the 95% confidence limits. Thick red line shows the prediction of from the Terray et al. [30] scaling and thin dashed line indicates estimates of from the LOW. All qualified segments are subject to high-passed filter with a cut-off frequency of 0.3 Hz before calculating shows TKE enhancement due to the proximity of the MATS to the sea surface. Here = 0.41 is the von Kármán constant, u * w denotes the water-side friction velocity (Fig. 15a), and z is the depth (of MATS) relative to the mean sea level. For the TKE dissipation generated by the wave breaking, Terray et al. [30] (hereafter T96 scaling) suggested the following scaling at depths which decays as z −2 : where F k = F u 3 * w is the wind energy input to the wave field, and F is a function of wave age. While observations contain waves with variable ages mostly greater than 25 (Fig. 2-c), we simplify the dependence of F to the wave age by using a constant value of 156, red curve in Fig. 15 (note that the value of F is 78 based on the best fit with our observed ), see also [3]. Comparisons show a broad agreement with those estimated from the T96 scaling ( Fig. 15) except for the period when the wind speed drops from 12 m s −1 to approximately 6 m s −1 and the platform experiences significant wave-induced motions (i.e AOA HF > 7 o and the rms-value of roll ≥ 1 o , Fig. 2c, d). Investigation of responsible mechanisms for this energy enhancement, at days between 1.9-2.2, is not the scope of present study and will be addressed elsewhere. Moreover, reprocessing the entire qualified data, using the application of the AWT method alone leads to approximately 2 times larger values of (not shown).
While conducting statistical tests between two sample distributions of depends on the vertical/horizontal separation between two sample spaces, the depth-attenuation of surface waves, and turbulence intermittency, we employ a two-sample t-test to check the null hypothesis that two sets of samples have identical means (i.e. both sample spaces originate from the same population) for a MSS-depth segment of 9-11 m. By assuming both sample sets have Gaussian distribution, an appropriate test statistic of the difference between the two means is given by where A and B are the means of the MSS and the MATS extracted , respectively, A and B denote the variances of the MSS and the MATS estimated , respectively, and n A and n B are sample sizes. For the significance level , we reject the null hypothesis if the measured t is greater than the critical t-value. Table 1 comprises the results of two-tailed t-test with n A − n B − 2 degrees of freedom at a 5% significance level over the depth range of 9-11 m. The calculated t for both sample spaces of are smaller than the critical value, t * , suggesting that the null hypothesis cannot be rejected without another cause (i.e. at a 5% significance level, the two means are identical in statistical sense). The resulting p-values suggest somewhat identical chance of observing extreme t values in the measured shear time series. Note that the values given in Table 1 are sensitive to the satisfaction of underlying t-test assumptions and the appropriate choice of the MSS depth-range. Hence, the null hypothesis may be rejected if we extend the thickness of the MSS depth-range to include deeper depths.

Conclusions
As a result of kinematical (wave advection) and dynamical (large orbital velocities) interactions between waves and turbulence in the upper ocean boundary layer, the shape of shear spectra comprises a narrowband inertial subrange (typically in the 7-30 cpm wavenumber range) together with appearance of significant wave-induced variances over a broad range of frequencies/wavenumbers. Such permanent wave-and motion-induced contaminations have found to be responsible for the overestimation in the values of dissipation rates of TKE from measurements carried out by the Moored Autonomous Turbulence System (MATS).
In this study, we examined the application of an adaptive amplitude thresholding method based on wavelet transformation and Shannon entropy envelope extraction to detect and reduce non-stationary and non-turbulent structures from the MATS's shear probe measurements near the sea surface. The algorithm segmented each 60-s time-sequence (affected by the wave-and motioninduced disturbances) into a finite number of subblocks from which the longest stationary block, with a size larger than 16-s, was used along with the Goodman et al. [13] method to provide a statistically significant and reliable estimate for the TKE dissipation, . Comparisons with the measurements from a nearby deployed free-falling microstructure profiler demonstrated that the proposed method could improve the accuracy of measurements from a subsurface moored instrument by a factor of 5-10 when comparing with those calculated from cleaned shear spectra using only Goodman et al. [13] method. We attributed the discrepancies to the intermittency of turbulence, different sampling methodology, instability of measuring sensors, (stochastic/deterministic) wave-induced contaminations, and variability in the depth-attenuation of wave field. Comparisons of the observed against those from the classical wall-layer and the scaling suggested by Terray et al. [30] confirmed again skill of the AWST algorithm in reduction of the (deterministic) wave-related disturbances from turbulence measurements near the wavy interface at frequencies/wavenumbers before the spectra start to roll off. We further conducted an analysis by applying the Shannon entropy segmentation (i.e. SES) without using the adaptive wavelet thresholding to identify the nonturbulent features from the turbulent motions. While both approaches (i.e. AWT and SES) operate almost equally-well for long segments, using each of them alone to perform segmentation decreases the probability of detection of long enough subblocks with sufficient quality (i.e. false detection) in our datasets. It should be noted that calculations of dissipation rates using 16-s long segments may not adequately resolve the turbulent intermittency and important turbulent scales depending on environmental conditions, e.g. small away from the wavy interface. Furthermore, the effectiveness of the proposed method depends on the choice of the optimum wavelet family, its maximum level of decomposition, and other model tuneable parameters (i.e. , * , and ). We found that the method is able to return reliable estimates of for conditions where the accelerometer data coherent with the shear signals are not available or problematic. The method with appropriate values of its parameters is also applicable for any datasets (e.g.. ADV) containing strong contributions from the non-turbulent structures interacting weakly with turbulence. constructed by the following Shannon entropy measure for each selected wavelet family from our library: The value of S from Eq. (16) reflects the sparsity of wavelet coefficients for a particular wavelet function at which the lower value corresponds to the higher degree of energy concentration over frequencies associated with the non-turbulent structures. We have applied the Shannon entropy measure to all 60-s shear time series. Table 1 contains the average values of S over all segments and indicates that the Db32 mother wavelet is the (near) optimal choice based on our library and particular value of J.
Since both Db32 and Db16 perform equally-well in our configuration, we use Db16 wavelet family in our analysis in order to increase the robustness of method.

Appendix 2: Stationarity: reverse arrangements test
While a large number of methods for stationarity test have been suggested in the literature (e.g. Runs test, unit root tests, the autocorrelation test, and wavelet-based test), we use a non-parametric stationarity test in weak sense socalled reverse arrangements test [5] in this study. Let s 1 , s 2 , ⋯ s N be an arbitrary (discrete) shear time series, where N denotes the number of sampling points. A Reverse Arrangement (RA), h ij for j > i , occurs when a subsequent observation, s j , in shear vector is smaller than the current shear observation, s i . The number of RAs for the ith entry is determined by The total number RAs for all observations is then calculated from A = ∑ N k=1 A i . By comparing the value of A with that of estimated from a realisation of a stationary random process, the stationarity of the vector is checked according to the RAs test at significant level (here = 0.05 ). For example, the number of RAs for a stationary time series with length N follow approximately normal distribution with mean A = N(N − 1)∕4 and variance

Appendix 3: Wave-induced variance: spectral model
To assess the efficiency of the integration method for estimating the wave-induced variances, we alternatively use a spectral model which is less affected by the complications of kinematic interaction of non-turbulent waves, platform motion, and turbulent fluctuations. As shown in Fig. 16, the (demeaned) vertical velocity spectrum of the ADV can be divided into: (i) a low frequency area (i.e. f ≪ f p ); (ii) a frequency band contaminated by waves and platform motions at a frequency band which is identified by ̃I = [0.5f p , f p + 0.5] , where f p in Hz denotes the frequency of the wave peak estimated from the AWAC data; (iii) inertial subrange; and (iv) noise affected subrange. For three velocity components, we fit the empirical Kaimal spectrum to the one-sided observed spectra, where denotes either vertical or horizontal turbulent velocities [12]: using a two-parameter ( ̂2-k o ) least-squared minimisation, where k indicates the wavenumber and k o is the wavenumber associated with the dominant energy containing scales, and ̂2 is the fitted turbulent velocity variance [17]. Fitting procedure is performed for frequencies between Vertical velocity frequency spectrum, w (f ) , of a 60-s segment measured from the MATS-mounted ADV (black) and the spectrum extracted from fitting empirical Kaimal spectrum, ̂w , to the observed spectrum at frequencies much lower and higher than the wave-affected band. The vertical dashed lines show the lower f l p = 0.5f p , and the upper f u p = f p + 0.5 bounds for the wave-affected subrange (coloured area). The wave peak frequency is estimated from the AWAC data 0.05 and 3 Hz by omitting the spectral energy within the wave-affected frequency band Ĩ (Fig. 16, red line). Here, the energy at frequencies much higher than the wave band is assumed to not be affected by the unsteady wave advection, and we select few sampling points for curve-fitting at frequencies much below the wave band to guarantee minimum contribution of the kinematical wave effects. The procedure is significantly sensitive to the higher noise floor known for the horizontal components of velocity. By assuming no correlation between the non-turbulent wave scales and turbulence motions, the wave-induced variance is calculated as where 2 is the contaminated observed variance. The total s t a n d a rd d e v i a t i o n i s t h e n c a l c u l a te d by ̃= √̃2 u +̃2 v +̃2 w .

Appendix 4: Wave advection
In the absence of dynamic coupling between waves and turbulence, Lumley and Terray [19] showed that the random waves are able to modify the slope of inertial subrange in frequencies too far from the wave peak frequency. Trowbridge and Elgar [31] measured by a single point measuring system in a dominantly oscillatory flow. They applied the modifier I(R, ) (i.e. Eq. 4) to correct for the elevation of energy levels in the inertial subrange due to the wave advection in order to estimate , see variation of I for different values of R and in Figure 17a.
Assuming wave and turbulence are independent processes and using Taylor's frozen turbulence hypothesis, for locally homogeneous and isotropic turbulence, the timedependent covariance tensor under the influence of random waves (propagating along the x-axis) is obtained as follows where i = √ −1 is the imaginary unit, dk = dk 1 dk 2 dk 3 , and c mn is a wave-related characteristic function that to the first order can be determined from the cross-spectrum of the wave orbital excursion in the x m and x n directions based on the linear wave theory [19]. The spectral tensor of the isotropic turbulence is defined by lj (k) = E(k) 4 k 2 ( lj k 2 − k l k j ), w h e r e t h e m o d e l e n e r g y s p e c t r u m E(k) = 1.5 2∕3 k −5∕3 f ( k)f ( k) is a function of and the large scale turbulence length scale , = 3∕4 −1∕4 is the Kolmogorov length scale, and f and f are universal shape functions [27]. According to the derivation of Lumley and Terray [19], we use = 1 m, = 1.4 × 10 −6 m 2 s −1 , and = 10 −6 m 2 s −3 to generate energy power spectrum by incorporating the effects of wave advection as In calculation of synthetic power spectrum, here only for the vertical velocity w (f ) = S 33 (f ) , from Eq. (22), the wave fields are generated from the Joint North Sea Wave Project (JONSWAP) fully-developed spectrum for the wave peak period of T p = 8 s and significant wave heights of H s ∼ 0 m and H s = 2 m, respectively. Figure 17b shows spectra (22) S lj (f ) = ∫ lj (t)e −2 if dt. Fig. 17 a The curve of the modifier I(R, ) for different values of R and calculated from Eq. (4); and b effects of wave advection on the shape of vertical velocity spectra for steady (blue line) and unsteady (black line) currents. Applying the correction factor I removes/reduces the advection effect for frequencies much higher than f p (red line). Wave energy spectrum is calculated by the JON-SWAP spectrum [3] and the wave-related characteristic functions are calculated using the linear wave theory and orbital velocity cross-spectra of vertical velocity in the absence of wave dynamical/ kinematical effects (reference spectrum, blue line), and for R = 0.3 (black line). The shape of the power spectrum at frequencies much higher than f p clearly shows departure from the reference spectrum with an elevation in the level of energy. Furthermore, less variances become appeared at frequencies below f p compared with the reference spectrum as a result of wave advection. Using the modifier I(R, ) shows successful removal of energy elevation induced by the wave advection. On the basis of these results, it is concluded that the wave advection for the small values of R varies the properties of turbulence and the shape of shear spectrum at frequencies much higher than the wave band and about f p (Fig. 17b). Therefore, an appropriate frequency/wavenumber correction is required to systematically adjust the spectral levels at frequencies affected kinematically by the wave orbital motions.