Calculation of turbulent dissipation rate with acoustic Doppler velocimeter

Turbulent dissipation rates are calculated from an acoustic Doppler velocimeter by fitting the measured wavenumber spectrum to a universal turbulence spectrum. A combination of the Butterworth filter and empirical model decomposition is employed to reduce Doppler noise and high‐frequency fluctuations. Different from the classical inertial subrange dissipation method fitting to the Kolmogorov − 5/3 slope, we propose the method here which can make longer available spectrum bands. We analyzed 408 bursts with turbulent dissipation rates ranging from 10−8 to 10−5 W kg−1 as measured in the coastal ocean of the South China Sea and found for all of these bursts, features of the clean spectrum can be resolved to the dissipation range of turbulence, in which about 15% bursts can be resolved to the Kolmogorov wavenumber, and 51% to 1/2 Kolmogorov wavenumber. Comparisons of this method with the inertial subrange method indicated that its estimated turbulent dissipation rates were somewhat smaller than those from the inertial subrange method in most of the bursts. Their ratios had a mean value of 0.77, which is typical for oceanic turbulence measurements.

Oceanic turbulence, and associated vertical mixing, plays an important role in the modulation of oceanic circulation and climate. Proper measurement and calculation of the turbulent dissipation rate are essential for the parameterization of vertical mixing, which is crucially important for the development of oceanic circulation models. The acoustic Doppler velocimeter (ADV) is extensively used to oceanic turbulence measurements. It can resolve accurately turbulence intensity of the vertical component of measured velocities (Voulgaris and Trowbridge 1998). Observations by Lien and D'Asaro (2006) showed that the velocity spectrum of the ADV can resolve turbulent dissipation rates as low as 10 28 W kg 21 .
Typically, the turbulent dissipation rate is estimated using the inertial subrange dissipation method for the ADV (Liu and Wei 2007;Lozovatsky et al. 2008;Bluteau et al. 2011). For homogeneous isotropic turbulence with high-Reynolds numbers, the slope of the velocity spectrum remains almost constant in an intermediate range, named as the inertial subrange, where motions are determined by inertial effects and viscous effects are negligible. Experimental results by Saddoughi and Veeravalli (1994) showed that the values of this slope differ by no more than 10% throughout the inertial subrange. Subsequently, the turbulent dissipation rate E can be obtained by fitting the inertial subrange of the measured spectrum A k r ð Þ to the theoretical slope, i.e., the Kolmogorov 25/3 curve: in which k r is the radian wavenumber and a is an empirical constant, which is usually set to 0.71 for the vertical velocity (Sreenivasan 1995;Lien and D'Asaro 2006). The ADV measures a series of velocities at a single point. The wavenumber spectrum A k r ð Þ is usually transformed to the frequency spectrum A f ð Þ using Taylor's frozen turbulence hypothesis (Taylor 1935) to give: where f is the frequency and U is the mean velocity. Thus, E is estimated by rewriting Eq. 2 (Liu and Wei 2007) as Here, the angle bracket represents averaging in the inertial subrange.
As its name suggests, the inertial subrange dissipation method is only valid in the inertial subrange. At wavenumbers beyond this subrange, the shape of the spectrum deviates greatly from the 2 5/3 curve in log-log space (Pope 2000). In the coastal ocean, the effects of surface waves always superimpose on the inertial subrange, whose extent is therefore significantly reduced. The turbulent anisotropy induced by the buoyancy and velocity shear can also reduce the extent of the inertial subrange from the low wavenumber end of spectrum (Gargett et al. 1984;Saddoughi and Veeravalli 1994). In some cases, its effect is beyond the range contaminated by surface waves (Doron et al. 2001). Thus, the inertial subrange may be quite short and identifying it becomes challenging. Under extreme conditions, the inertial subrange vanishes (Gargett et al. 1984;Luznik et al. 2007) and then the inertial subrange method is unusable for estimating the dissipation rate (Bluteau et al. 2011).
Compared with the inertial subrange, the dissipation range of the turbulent spectrum is less susceptible to the turbulent anisotropy, as well as contaminations from wave motions (Gargett et al. 1984;Saddoughi and Veeravalli 1994;Doron et al. 2001). The Panchev-Kesich (P-K) spectrum (Panchev and Kesich 1969) is believed to be a theoretical form of the universal turbulence spectrum and represents an extension to the Kolmogorov's turbulence theory (Wolk et al. 2002). In this paper, we attempt to calculate the turbulent dissipation rate by fitting the measured shear spectrum of the ADV to the P-K spectrum.

De-spiking
It is well known that turbulence measurements from the ADV suffer from random spikes, Doppler noise, installation vibrations, and other unknown disturbances (Nikora and Goring 1998). Therefore, adequate postprocessing, including de-spiking and removal of noise, is necessary in order to obtain clean spectra for the ADV datasets.
The spikes in the ADV datasets are mainly induced by Doppler signal aliasing (Voulgaris and Trowbridge 1998). The true 3D phase space method obtained by Mori et al. (2007) is used to remove spikes in time-series of the ADV velocity data in this paper. This method is developed from the phase-space thresholding method of Goring and Nikora (2002). It can significantly remove the spikes in comparison with the original method.

Removal of noise
Different from the spikes, Doppler noise is inherent to the measuring principle of the ADV (Voulgaris and Trowbridge 1998). The empirical model decomposition (EMD) is an adaptive and temporally local analytic method that is especially effective for nonstationary and nonlinear data. It can decompose a complicated dataset into a series of intrinsic mode functions (IMFs) and a residual trend (Huang et al. 1998). This method is a general data-driven method without any a priori basis functions, and the signal is decomposed based on the local characteristics of the data. It has been successfully applied to the eddy-covariance analysis of air-sea interface turbulent fluxes (Wang et al. 2013), ADV data processing (Qiao et al. 2016), and geophysical studies (Huang and Wu 2008) Qiao et al. (2016) showed that the EMD can effectively remove the noise in the ADV data, in which the ADV vertical velocity was decomposed into 12 IMFs, and components 1 and 2 were treated as noise because they have characteristics of white noise. In our ADV datasets, the first few components of IMFs are also associated with Doppler noise and high-frequency contamination. Therefore, it is necessary to remove these components to obtain a clean spectrum from the original ADV dataset.
For the raw shear spectrum of the ADV, the highfrequency noises are usually too strong, so that a filter must be used before conducting the EMD. Otherwise, these noises can spread to the whole band of spectrum when applying the EMD. A low-pass Butterworth filter is applied in our analysis, which is usually used in data processing of turbulence (Roget et al. 2006). Its cutoff wavenumber is usually set to the wavenumber at which the original spectrum begins to deviate significantly from the theoretical P-K curves. In this study, it is simply set to the wavenumber at which the predicted P-K spectrum from the original spectrum deviates from the 99% confidence interval of the original spectrum.

Calculation of the turbulent dissipation rate
The 1D P-K spectrum proposed by Roget et al. (2006) is used, whose analytical approximation equation is where k n 5 k/k c , k is the cyclic wavenumber in cpm, is the Kolmogorov wavenumber (Wolk et al. 2002), and m is the kinematic viscosity dependent on the temperature using a polynomial approximation: m5 1:79274720:05126103 T10:0005918645T 2 À Á 310 26 (6) where, T is the temperature of sea water in 8C (Carniel et al. 2012). Thus, it is about 8.4 3 10 27 m 2 /s for the seawater of 278C. Compared with the Nasmyth universal spectrum (Nasmyth 1970), the P-K spectrum contains more power at lower wavenumbers and rolls off slightly faster at high wavenumbers (Roget et al. 2006). The turbulent dissipation rate E is estimated by fitting the measured ADV shear spectrum to the theoretical P-K spectrum. The ADV measures instantaneous 3D velocities. For every burst, the double rotation method is applied to rotate velocities into an along-stream component u, a cross-stream component v, and a vertical component w (Voulgaris and Trowbridge 1998;Wilczak et al. 2001). Due to the geometry structure of ADV sensor, the noise of its horizontal velocity components is much larger than that of vertical component (Lohrmann et al. 1994;Nikora and Goring 1998).Therefore, the w is used to calculate the power spectrum and shear spectrum. Both of them are calculated using the Welch's method, in which the 180 s dataset is divided into eight segments with 50% overlap, tapered with a Hamming window. A fast Fourier transform of 1024 points is used in the power spectral density estimate.
Under the Taylor's frozen turbulence hypothesis, the velocity shear can be calculated from @w @x 5 @w U@t (7) where, @t is the sampling interval and U is the mean advection velocity across the sampling volume of ADV. The measured frequency spectrum w f ð Þ can be transformed into wavenumber spectrum w k ð Þ using (Wolk et al. 2002) The upper and lower bounds of the wavenumber resolved by the ADV shear spectrum are defined as the wavenumbers at which the predicted P-K spectrum deviated from the 95% confidence interval of the clean spectrum. In this prescribed wavenumber range, the least square method is used for fitting the clean spectrum to the theoretical P-K spectrum shown by Eq. 4 and then estimating the turbulent dissipation rate.

Assessment
The datasets were collected near an Observation Platform of Marine Meteorology tower in the northern coast of the South China Sea (Qiao et al. 2016). This platform is located 7 km to the nearest shore with a water depth of 16 m. A Nortek ADV with an acoustic frequency of 6 MHz is used to measure instantaneous 3D velocities. This ADV is fixed to a bottom-mounted frame. Its probe pointes upward and is 0.73 m above the seafloor. The sampling volume comes from a cylinder with a diameter of 15 mm, a length of 19 mm, and a distance from the probe of 0.15 m. This somewhat large volume can provide substantial scattering material for the acoustic systems. Its sampling frequency is set to a rate of 64 Hz. The sampling time is set to 180 s for every 1800s. Thus, there are 1,1520 samples for every burst.
The 482 bursts were collected from 24 October to 3 November 2015. During this period, wave measurements from a Sentinel V ADCP were available. Of these bursts, 74 are discarded because their spectrum does not have a distinct 2 5/3 slope in the theoretical inertial subrange. These poor bursts are usually associated with weak flows or strong vibrations of the ADV. The remaining 408 bursts have mean velocities ranging from 0.04 to 0.23 m s 21 (Fig. 1A) and sea water temperatures from 27.0 to 27.68C.
During these 408 bursts, the surface waves have significant wave heights ranging from 0.45 to 1.49 m (Fig. 1B), and wave periods from 4.1 to 7.9 s (Fig. 1C). For flows affected by surface waves, the instantaneous velocity can be decomposed into three components as where, u w is the wave orbital velocity and u 0 is the turbulence velocity. Lumley and Terray (1983) proposed that the advection velocity U should be replaced by u w in the frequency-wavenumber transformation when the former is smaller than the latter. The EMD can be used to separate approximately the wave motions from the turbulence (Qiao et al. 2016). As a whole, U is larger than or close to u w at the sampling depth during these bursts (Fig. 1A). Therefore, U is used to estimate frequency spectra in our analysis, i.e., Eq. 2. A typical time series of vertical velocity from the ADV is shown in Fig. 2A. During this burst, the water depth is 15.4 m, and the mean velocity and sea water temperature across the sampling volume are 0.10 m s 21 and 27.48C, respectively. At the surface, surface waves are weak with significant wave height of 0.59 m and wave period of 6.7 s. Similar to previous studies (Goring and Nikora 2002;Cea et al. 2007;Jesson et al. 2013), there are obvious spikes in the time series for the raw data. For the de-spiked vertical velocity, its mean w is 1.5 3 10 23 m s 21 . This nonzero vertical velocity may be associated with changes of sea level induced by diurnal cycles of the tide. Thus, a double rotation is used to make w50 before calculating the velocity shear and spectrum (Fig. 2B).
The velocity spectrum calculated from the double-rotated and de-spiked vertical velocity is shown in Fig. 3A. It has an obvious peak at low frequencies due to the effect of surface waves, and a noise floor at high frequencies. As previous studies pointed out, identifying an inertial subrange band is challenging for the ADV spectrum (Bluteau et al. 2011). It seems that the spectrum was in good agreement with the expected 2 5/3 slope at frequencies from 0.3 to 3 Hz. However, studies have suggested the theoretical upper bound of inertial subrange is at 0.1 k c (Pope 2000). This upper bound is only about 15 cpm or 1.5 Hz in this burst. Thus, the estimated dissipation rate is about 7.4 3 10 27 W kg 21 by Eq. 3.
The time series of velocity shear for this burst is shown in Fig. 2C, and its wavenumber spectrum is shown in Fig. 3B, which presents a rapid increase at high wavenumbers. The values at high wavenumbers are larger than that at low wavenumbers by up to two orders of magnitude. The high-frequency noises are significantly reduced after applying a first-order lowpass Butterworth filter with a cutoff frequency of 2.4 Hz or a wavenumber of 24 cpm ( Fig. 2D; pink line in Fig. 3B).
The filtered velocity shear can be decomposed into different IMFs by the EMD (Fig. 4). When the first component of IMFs (mode 1) is removed from this filtered dataset, the noises are further reduced (green line in Fig. 2D). The original shear spectrum follows that predicted by the P-K spectrum at wavenumbers from 3 to 20 cpm and deviates slightly from it from 20 to 30 cpm (blue line in Fig. 3B). However, the clean spectrum is in good agreement with the P-K spectrum till the wavenumber up to 88 cpm (green line in Fig. 3B). The Kolmogorov wavenumber k c is about 148 cpm in this burst. This implies that the clean spectrum can resolve up to 0.6 k c . Fitting the clean spectrum to the P-K spectrum from 3 to 88 cpm, we can get the dissipation rate E PK 5 4.3 3 10 27 W kg 21 , which is somewhat smaller than that from the inertial subrange method.
Besides of Eq. 8, the shear spectrum w k ð Þ can be obtained directly from the velocity spectrum S(f) as In our analysis, the shear spectrum defined above is very similar as that by Eq. 8 at low frequencies (Fig. 5). At high frequencies, however, the former seems to be slightly larger than the latter. Although the EMD can simply reduce the noise  , the blue curve is the original shear spectrum, the pink curve is the spectrum filtered by a low-pass Butterworth filter, and the green curve is the clean spectrum removed furtherly mode 1 of IMFs. The grey dashed curves are the theoretical P-K spectrum, and the thick curves indicate its corresponding dissipation rates (in W kg 21 ). The thick red curve is the predicted P-K spectrum from the clean spectrum. The k low and k up are lower and upper bounds of the wavenumber range for calculating the turbulent dissipation rate by fitting the clean spectrum to the P-K spectrum, respectively. The k cut is the cutoff wavenumber of the Butterworth filter.
from original ADV velocity, it fails to make longer available spectrum bands in the shear spectrum. There has still considerable noise at high-frequency region when the first mode of IMFs is removed, while a part of turbulence energy is removed when the second mode is further removed. The velocity and shear spectra of the other two typical bursts are shown in Figs. 6 and 7. Different from that in Fig. 3, the clean spectrum in Fig. 6 can be obtained only by filtering the original spectrum using a low-pass Butterworth filter. In Fig. 7, however, it is obtained by furtherly removing the first two components of IMFs (modes 1 and 2). This implies that the combination of the Butterworth filter and the EMD is not fixed for different bursts. Figure 8 shows the relationship of the combinations and the turbulent dissipation rates for 408 effective bursts in our measurements. The dissipation rates of these bursts span about three orders of magnitude with a range from 1 3 10 28 to 1 3 10 25 W kg 21 . It can be seen that at relatively strong turbulence, only a low-pass Butterworth filter is needed to reconstruct the clean spectrum; at low-turbulence energy, however, both the filter and the EMD are needed.

Discussion
Doppler noise has properties similar to that of white noise and introduces significant error in the measured velocity and shear spectra. Durgesh et al. (2014) proposed two methods, including the noise auto-correlation approach and the proper orthogonal decomposition approach, to attenuate Doppler noise and improve the velocity spectrum. Hejazi et al. (2016) used a linear correlation algorithm to lower the noise level and remove the spikes. There have also noisereduction methods requiring two independent and simultaneous measurements of the same vertical velocity component (Hurther and Lemmin 2001;Doroudian et al. 2010), which apply to four beam Doppler systems. Although these methods can decrease noise contaminations of spectra and then obtain a wider inertial subrange, there still has considerable energy of noise at high frequencies of corrected spectra (Khorsandi et al. 2012;Durgesh et al. 2014).
The developed method can produce a longer available spectrum band for calculations of turbulent dissipation rate. Figure 9 shows the cumulative percent of k up /k c for these 408 bursts. In our measurements, the advection velocities are relatively weak with values from 0.04 to 0.23 m s 21 (Fig. 1A),   spectrum of the shear. In (A), the blue curve is the original velocity spectrum, the pink curve is the spectrum removed the first mode of IMFs, the green curve is that removed its first two IMFs, and the red line is the theoretical Kolmogorov 25/3 curve. In (B), the blue, pink, and green curves corresponding to those in (A) are the shear spectra calculated by Eq. (11), respectively; the red curve is the original shear spectrum obtained from Eq. (8). which are helpful in producing a relatively highwavenumber band via Eq. 8. For these bursts, they all can resolve to 0.3 k c . The peaks in the P-K spectrum represent transitions from the inertial subrange to the dissipation range, which are at about 0.125 k c (Gregg 1999). This implies that the clean spectra of all the bursts cover part of the inertial subrange of turbulence but extend well into the dissipation range with suitable removal of noise. The P-K spectrum accounts for 96% of the total spectral variance at wavenumbers below 0.5 k c (Moum et al. 1995). The clean spectra of 51% bursts can resolve to 0.5 k c in our measurements.
Moreover, there are about 15% bursts whose clean spectra resolve to k c (see an example in Fig. 10). A direct comparison of the turbulent dissipation rates estimated by fitting to the P-K spectrum (E PK ) and those from the inertial subrange method (E in ) is shown in Fig. 11. The E PK is slightly smaller than E in in most of the bursts. Their ratios range from 0.5 to 1.4 with a mean value of 0.77 for all the bursts. This result is similar as that by Doron et al. (2001). Using a submersible particle image velocimetry (PIV) system, they compared different methods of estimating the turbulent energy dissipation and found that the estimate from fitting the universal spectrum was smaller than that by the inertial subrange method if the extended and interpolated data was used, and then was closer to that by the direct method. They gave two causes for this discrepancy: (1) the former has a longer available wavenumber band; (2)  dissipation rates. In the vertical coordinate, "1" denotes that only a Butterworth filter is needed for removal of noise, "2" denotes a combination of the Butterworth filter and removing mode 1 of IMFs, and "3" denotes a combination of the Butterworth filter and removing modes 1-2. bursts. The k up is the maximum wavenumber resolved by the clean spectrum, and the k c is the Kolmogorov wavenumber. contamination by wave motions is relatively weak in the dissipation range in contrast to the inertial subrange. Both of them can reduce uncertainties of turbulent estimate. It should be noted that accurate measurements of the turbulent dissipation rates are still challenging. Even for airfoil type shear sensors, the ideal systematic bias of the E calculation is estimated to be within a factor of two due to cumulative uncertainties (Moum et al. 1995). In fact, there may exist larger uncertainties in routine measurements (Stips and Prandke 2000). Thus, the discrepancy between E PK and E in is within the range of typical values for oceanic turbulence measurements.

Conclusions
The ADV is an alternative instrument for oceanic turbulence measurements. In this paper, we assess the capability of estimating turbulent dissipation rates from the ADV by fitting to the universal turbulence spectrum. In this method, the removal of noise of the original shear data is the most critical part. A combination of the Butterworth filter and the EMD is suggested to reduce Doppler noise and highfrequency fluctuations. Compared with the classical inertial subrange dissipation method fitting to the Kolmogorov 2 5/3 slope, this method usually produces a longer available spectrum band, which contributes to enhancing the confidence of calculation of turbulent dissipation rates.
The 408 efficient bursts collected in a coastal area are analyzed using this method. These bursts have mean velocities ranging from 0.04 to 0.23 m s 21 , and turbulent dissipation rates from 10 28 to 10 25 W kg 21 . The results show that the method of removal of noise depends on the magnitude of the turbulent dissipation rate of the burst. For relatively strong turbulence, a low-pass Butterworth filter is enough to reconstruct the clean spectrum. For low turbulence, however, both the filter and the EMD are needed.
Under suitable removal of noise, the clean spectra can resolve to the dissipation range of oceanic turbulence for all of these bursts, in which about 51% bursts can resolve to 0.5 k c , and 15% bursts to k c in our measurements. Comparisons of this method with the inertial subrange method indicate that the turbulent dissipation rates estimated by the former are somewhat smaller than those from the latter in most of the bursts. Their ratios range from 0.5 to 1.4 with a mean value of 0.77. These discrepancies are typical for oceanic turbulence measurements.  Fig. 11. Scatter comparison of the turbulent dissipation rates estimated by fitting to the P-K spectrum (E PK ) and those from the inertial subrange method (E in ). The blue, pink and green dots are corresponding to the method of removal of noise 1, 2, and 3 in Fig. 8, respectively.