Novel approaches to estimating the turbulent kinetic energy dissipation rate from low- and moderate-resolution velocity fluctuation time series

In this paper we propose two approaches to estimating the turbulent kinetic energy (TKE) dissipation rate, based on the zero-crossing method by Sreenivasan et al. (1983). The original formulation requires a fine resolution of the measured signal, down to the smallest dissipative scales. However, due to finite sampling frequency, as well as measurement errors, velocity time series obtained from airborne experiments are characterized by the presence of effective spectral cutoffs. In contrast to the original formulation the new approaches are suitable for use with signals originating from airborne experiments. The suitability of the new approaches is tested using measurement data obtained during the Physics of Stratocumulus Top (POST) airborne research campaign as well as synthetic turbulence data. They appear useful and complementary to existing methods. We show the number-of-crossings-based approaches respond differently to errors due to finite sampling and finite averaging than the classical power spectral method. Hence, their application for the case of short signals and small sampling frequencies is particularly interesting, as it can increase the robustness of turbulent kinetic energy dissipation rate retrieval.


Introduction
Despite the fact that turbulence is one of the key physical mechanisms responsible for many atmospheric phenomena, information on the turbulent kinetic energy (TKE) dissipation rate based on in situ airborne measurements is scarce. Research aircraft are often not equipped to measure wind fluctuations with spatial resolution better than a few tens of meters (Wendisch and Brenguier, 2013). Due to various problems related to, for example, inhomogeneity of turbulence along the aircraft track and/or artifacts related to inevitable aerodynamic problems (Khelif et al., 1999;Kalgorios and Wang, 2002;Mallaun et al., 2015), estimates of at such low resolutions using power spectral density (PSD) or structure functions are complex and far from being standardized (e.g., compare procedures in Strauss et al., 2015;Jen-La Plante et al., 2016). The following question arises: can we do any better or can we at least introduce alternative methods to increase the robustness of retrievals?
In the literature, there exist several different methods to estimate using the measured velocity signal as a starting point. One of them is the zero-or threshold-crossing method (Sreenivasan et al., 1983) which, instead of calculating the energy spectrum or velocity structure functions, requires counting of the signal zero-or threshold-crossing events (see Fig. 1a). Their mean number per unit length is related to the turbulent kinetic energy dissipation rate. The zero-crossing method is based on a direct relation between and the root mean square of the velocity derivative (Pope, 2000); hence, the measured signal should be resolved down to the smallest scales. However, this is not achievable in the case of flight measurements with moderate time resolutions. Using Taylor's hypothesis, the measured time series can be converted into a spatial signal and the sampling frequency will correspond to scales which are typically 2-3 orders of magnitude larger than the Kolmogorov scales. As a result, the number of zero crossings per unit length for such signal is Published by Copernicus Publications on behalf of the European Geosciences Union.
much smaller than the one corresponding to a high-resolution velocity signal where turbulence intensity is the same.
Interestingly, Kopeć et al. (2016) have shown that the dissipation rates estimated from such N L using very lowresolution signals, although underestimated, were proportional to calculated using structure functions scaling in the inertial range. In the follow up analyses we found that this is also the case for moderate-resolution airborne data from different sources. This led us to a question of whether it would be possible to modify the zero-crossing method such that it could also be applied to moderate-or low-resolution measurements whilst mitigating the observed underestimation at the same time. In this work we propose two possible modifications of the zero-crossing method. The first one is based on a successive filtering of a velocity signal and inertial range arguments. In the second approach we use an analytical model for the unresolved part of the spectrum and calculate a correcting factor to N L , such that the standard relation between and N L can be used.
The new approaches are tested on velocity signals obtained during the Physics of Stratocumulus Top (POST) research campaign, which was designed to investigate the marine stratocumulus clouds and the details of vertical structure of the stratocumulus-topped boundary layer (STBL; Gerber et al., 2013;Malinowski et al., 2013). The observed winds were measured using the CIRPAS Twin Otter research aircraft with sampling frequency f s = 40 Hz, which corresponds to a resolution of 2.75 m for the speed of the aircraft 55 m s −1 . Additional tests of the method with synthetic velocity signals as suggested by Frehlich et al. (2001) are also performed.
The present paper is structured as follows. In Sect. 2 we review existing methods to estimate the dissipation rate of the turbulent kinetic energy. Next, in Sect. 3, we propose the two modifications of the zero-crossing method. They are applied to a single signal from flight 13 and synthetic turbulence data and discussed in detail in Sect. 4. Next, in Sect. 5, we apply the procedures to several data sets from flights 10 and 13 to show that the results of new approaches compare favorably with those obtained from standard power spectrum and structure function methods. This is followed by "Conclusions", where the advantages of the new proposals and perspectives for further study are discussed.
2 Previous methods for the retrieval of the energy dissipation rate from measured velocity time series The need to estimate the turbulent kinetic energy dissipation rate as well as the variety of available data resulted in the formulation of a number of estimation methods. Two of the most commonly used approaches are the power spectral density and the structure function approach. Both are based on the inertial range arguments, which follow from Kolmogorov's second similarity hypothesis (Kolmogorov, 1941); hence, they are also called "indirect methods" (Albertson et al., 1997). With the assumption of local isotropy the one-dimensional longitudinal and transverse wavenumber spectra in the inertial range are given by (Monin and Yaglom, 1975;Pope, 2000) Here k 1 is the longitudinal component of the wavenumber vector k = (k 1 , k 2 , k 3 ), C 1 ≈ 0.49 and C 1 ≈ 0.65 if k 1 units are rad m −1 (cf. Pope, 2000, Eqs. 6.242, 6.243). E 11 is related to the energy spectrum function E(k): where k = |k|. As discussed in Pope (2000) experimental data confirm Eq.
(1) within 20 % of the predicted values of C 1 and C 1 over two decades of wavenumbers. Within the validity of the local isotropy assumption of Kolmogorov (1941), the energy spectrum function can be approximated by the formula (Pope, 2000) where C ≈ 1.5 as supported by experimental data, and f L and f η are nondimensional functions, which specify the shape of energy spectrum in, respectively, the energycontaining and the dissipation range. L = k 3/2 / denotes the length scale of large eddies and η = (ν 3 / ) 1/4 is the Kolmogorov length scale connected with the dissipative scales (Pope, 2000), where ν is the kinematic viscosity. The function f L tends toward unity for large kL, whereas f η tends toward unity for small kη, such that in the inertial range the formula E(k) = C 2/3 k −5/3 is recovered. Within the validity of Taylor's hypothesis, Eq. (1) can be converted to the frequency spectra, where k 1 = (2πf )/U and U is the magnitude of the vector difference between the aircraft velocity and the wind velocity, i.e., the true air speed. The vector difference is averaged along the displacement which defines k 1 . The frequency f is measured in 1 s −1 , U in m s −1 and k 1 in rad m −1 . In order to estimate the dissipation rate from the atmospheric turbulence measurements, several assumptions should be taken. Most importantly, one assumes that the turbulence is homogeneous and isotropic and that the inertial range scaling Eq. (1) holds. Then, the frequency spectrum of the longitudinal velocity component in the inertial range is (e.g., Oncley et al., 1996;Siebert et al., 2006) The value of a constant C 1 ≈ 0.49 used in this work is related to the one-sided spectra. Hence, by E 11 , E 22 or S(f ) we denote the one-sided spectra, which integrated over argument from 0 to ∞ yield the variance of the signal. With Eq. (4), the turbulent kinetic energy dissipation rate can be estimated from the power spectral density of the measured signal. Alternatively, one can consider the nth order longitudinal structure functions D n = (u L (x +r, t)−u L (x, t)) n ; here u L is the longitudinal component of velocity and r is a displacement along the direction defined by u L . In the inertial subrange, the second-and third-order structure functions are related to the dissipation rate by the following formulas (Pope, 2000): Experimental results of Saddoughi and Veeravalli (1994) indicate that C 2 ≈ 2, with an accuracy of ±15 %. Another method, also based on Eq.
(3), is the velocity variance method (Fairall et al., 1980;Bouniol et al., 2004;O'Connor et al., 2010). Let us consider a homogeneous velocity field converted to time series u(t) with the use of Taylor's hypothesis. The mean-square value of this signal u 2 (t) = u 2 is equal to the integral from 0 to ∞ of the onesided power spectral density S(f ) over the frequency space.
The signal u(t) is next filtered with a band-pass filter with cutoff numbers [f low , f up ] in the frequency space. Assuming that the filter is perfect, i.e., that it is a rectangle in the frequency space, after the filtering, a signal u f (t) with the variance is obtained. The above formula represents the portion of kinetic energy of u(t) contained in the frequencies between f low and f up . Fairall et al. (1980), Bouniol et al. (2004) andO'Connor et al. (2010) substitute Eq. (3) for S(f ) into Eq. (6) and integrate to obtain the following formula for the dissipation rate: Yet another method, also used in the atmospheric turbulence analysis (Sreenivasan et al., 1983;Katul, 2009, 2010;Wilson, 1995;Yee et al., 1995), is based on the number of zero or level crossings of the measured velocity signal. It dates back to the early work of Rice (1945), who considered a stochastic processes q and its derivative with respect to time ∂q/∂t. He then assumed that these two processes have Gaussian statistics and are independent. The formulation of this method results from investigating how frequently the signal crosses the level zero q(t) = 0 (see Fig. 1a). Working under those assumptions Rice (1945) showed that the number of upcrossings of the zero level per unit time is As (∂q/∂t) 2 is proportional to the dissipation rate of the kinetic energy, the zero-crossing method can be used to estimate this quantity. As it was argued by Sreenivasan et al. (1983), Eq. (8) also holds with less restricted assumptions, with only q having Gaussian statistics, and, moreover, even for strongly non-Gaussian velocity signals, the number of zero crossings was close to the theoretical value from Eq. (8).
For a spatially varying signal, Eq. (8) can be expressed as follows, using the characteristic wavenumber k c and the onesided wavenumber spectra (He and Yuan, 2001): The characteristic wavelength is equal to λ c = 2π/k c . Hence, the mean number of crossings (up-and downcrossings) per unit length N L with, on average, two crossings per λ c is We will now introduce the two-point longitudinal correlation of velocity R ij (r 1 e 1 ) = u i (x, t)u j (x + r 1 e 1 , t) , where e 1 is the standard basis vector, and assume that the flow is statistically stationary and homogeneous and that statistics do not depend either on time t or point x.
Using the inverse Fourier transform, the 11 component of the two-point correlation tensor R 11 and its derivatives can be written in terms of E 11 as follows (Pope, 2000): where R 11 denotes the second-order derivative of R 11 . With those relationships we can rewrite Eq. (9) in the following manner: We further define the Taylor longitudinal microscale λ f with the use of R 11 (0) and R 11 (0): Hence, Eq. (10) implies that the number of crossings per unit length is related to the longitudinal Taylor microscale λ f through Relations (11-14) are valid for any statistically homogeneous vector fields, regardless of whether or not they are isotropic (Monin and Yaglom, 1975), provided that k c is the characteristic wavenumber along the longitudinal direction. However, homogeneity alone is not a sufficient assumption to estimate the TKE dissipation rate of a 3-D turbulent field from velocity signals measured along the 1-D aircraft flight path (Chamecki and Dias, 2004). We further use the local isotropy assumption to write a relation between dissipation and the Taylor microscales (Pope, 2000): where λ g = λ f / √ 2 is the Taylor transverse microscale. Hence, finally, substituting Eq. (14) into Eq. (15) we obtain (Poggi and Katul, 2010 For the transverse velocity time series, Eq. (16) has a factor of 7.5 instead of 15.
3 New proposals for the estimation of the dissipation rate from a velocity signal with a truncated high-frequency part of the energy spectrum Based on Eqs. (9) and (10) it is clear that the number of zero crossings is related to the 11 component of the dissipation tensor D 11 (k) = 2νk 2 E 11 (k): (Pope, 2000); here β = 2.1 and η = 2 mm. It is seen that the large wavenumber (small scale) part of the spectrum has the most significant impact on the resulting value of N L . At the same time the data available from the POST measurements can only account for a small part of the total dissipation spectrum, shown qualitatively as a shaded region in Fig. 1b. The lower bound of this region follows from a finite size of the averaging window, while the upper bound is related to the finite Nyquist frequency which equals 20 Hz for the POST measurements.
If one was to use the zero-crossing method (Eq. 16) in order to estimate , it is clear that the measured number of signal zero crossings would lead to a significant underestimation of the spectrum integral as compared to the full spectrum measurements down to the very small scales. We would like to propose the reformulation of the original zero-crossing method in order to estimate the dissipation rate from the number of signal zero crossings based on a restricted range of k values available from the airborne measurements. Two proposals for such procedures are given later in the article.

Method based on successive filtering of a signal
Let us consider a signal u 1 (t) resolved in a certain range of frequencies f 0 < f < f 1 . Converting the wavenumber spectrum to the frequency spectrum, we obtain from Eq. (17) the following relation for the number of signal crossings per unit time Similarly as in the velocity variance method described in Sect. 2, let us now filter the signal using a band-pass filter characterized by a different cutoff frequency f 2 < f 1 . In such a case we obtain a different signal u 2 (t) with a reduced number of zero crossings N 2 < N 1 : If we subtract Eq. (19) from Eq. (18) we obtain In the inertial range, S(f ) is described by Eq. (4); hence, if both f 1 and f 2 belong to the inertial range, If we proceed further and filter the signal using a series of cutoff frequencies f i < f 2 , we can estimate from Eq. (21) using a linear-least-squares fitting method.
In the above derivation we assumed that the applied filter is rectangular in the frequency space. The issue of frequency response characteristics of a filter will be discussed further in Sect. 4.1.

Method based on recovering the missing part of the spectrum
In this method we attempt to account for the impact of the missing part of the dissipation spectrum by introducing a correcting factor to the number of zero crossings per unit length N L . The number of crossings per unit length is calculated from the measured signal where the fine-scale fluctuations having the highest wavenumber k cut will be denoted by N cut and the variance of this signal will be denoted by u cut 2 . From Eq. (17) it follows that N cut is related to N L by the formula We then assume a certain form of the energy spectrum, Eq.
(3). For simplicity we take f L = 1, i.e., we neglect the contribution of the largest scales to the value of the dissipation rate based on zero crossings and we consider two different forms of f η , as proposed in Pope (2000). The first is the simple exponential form with β = 2.1. The second is the more complex formula where β = 5.2 and c η = 0.4. With this, the energy spectrum reads where C = 1.5. The integral from 0 to ∞ of the dissipation spectrum 2νk 2 E(k) should be equal to , which results in β = 2.1 in Eq. (23) and provides a relation between β and c η in Eq. (24). The latter case, due to the additional degree of freedom in f η fits the experimental data better in the dissipative range (Pope, 2000). The corresponding one-dimensional spectrum E 11 can be calculated from Eq. (2): Next we change the variables in the integral Eq. (26) to ξ = βkη, introduce Eq. (26) into Eq. (22) and once again change the variables to ξ 1 = βk 1 η. As a result we obtain here C F is the correcting factor. The value of can be calculated numerically using an iterative procedure.
As a starting point for this procedure, a first guess for the Kolmogorov length η = (ν 3 / ) 1/4 should be given. With this, we calculate the correcting factor C F from Eq. (27) taking either the form of Eq. (23) or (24) for f η . Next, from Eq. (16) the value of dissipation can be estimated as We start the next iteration by calculating again the Kolmogorov length η = (ν 3 / ) 1/4 , the corrected value of C F from Eq. (27) and the new value of from Eq. (28). After several iterations the procedure converges to the final values of the dissipation rate and Kolmogorov's length η with an error defined by a prescribed norm η = |η n+1 −η n | ≤ d η . The successive steps are summarized in the form of Algorithm 1.
Algorithm 1 Procedure of iterative determination based on missing spectrum part recovery.
It should be noted that in this approach we do not have the empirical inertial range constant C, and we calculate the dissipation rate directly from the formula with viscosity, Eq. (28), as in the original zero-crossing method (see Eq. 16 and Poggi and Katul, 2010). In order to present the more detailed properties of the procedure, we used the velocity signal from one of the horizontal flight segments that took place within the turbulent atmospheric boundary layer. This segment was a part of flight 13 of the POST airborne research campaign Malinowski et al., 2013). The data were provided in the east, north, up (ENU) coordinate system. For further study we have calculated time series of the longitudinal velocity component along the track. The signals' sampling frequency was f s = 40 Hz and the duration was T = 438.75 s. The magnitude of the vector difference between the aircraft velocity and the wind velocity U , averaged over the track vector, was about 55 m s −1 and the standard deviation u = 0.28 m s −1 .
We have estimated the dissipation rate based on the number of zero crossings, according to the methods outlined in Sect. 3.1. The average dissipation rate calculated from the frequency spectrum and the structure function for the whole flight fragment Eqs. (4) and (5) was close to equal, with PSD = 2.48 × 10 −4 and SF = 2.52 × 10 −4 m 2 s −3 , respectively. These values were obtained from the linear-leastsquares fit procedure in the range f = 0.3-5 Hz for the frequency spectrum and r = 11-183 m for the structure function (see Fig. 2).
Before applying the threshold-crossing procedures the signal had to be filtered in order to eliminate errors due to largescale tendencies as well as small-scale measurements noise. For this purpose we used the sixth-order low-pass Butterworth filter (Butterworth, 1930) implemented in Matlab ® . Figure 3 presents the velocity signal over t = 50 s before fil- tering (top graph) and the same signal after filtering with f cut = 5 and f cut = 1 Hz.
The original formula of Rice (1945) was derived for the case when both the signal and its derivative have Gaussian probability density functions (PDFs) and are statistically independent. In general, such assumptions do not hold for turbulent signals. However, as discussed by Sreenivasan et al. (1983), experimental observations and further theoretical studies indicate that the formula of Rice (1945) has a more general applicability than for what it was mathematically proven for originally and is satisfied with a fair accuracy even for the case of strongly non-Gaussian signals. Figure 4a presents PDFs of the normalized original signal and the filtered signals compared with the normalized Gaussian distribution. As it is seen, filtering does not lead to significant changes in the investigated PDFs.
It is worth noting that the spectra (f 2 S(f ), Fig. 4b) display a peak at f = 10 Hz. This phenomenon has been indicated in the previous analyses of POST (Jen-La Plante et al., 2016) and appears due to measurement errors. However, as the highest cutoff frequencies used in the present study are 5 Hz, it should not affect our results.
In order to use the method based on successive signal filtering we filtered the signal with different values of f cut in the range f cut = 0.1−19 Hz. For each f cut = f i we calculated the number of zero crossings N i based on the filtered signal. The zero-crossing event was detected when the product of two consecutive values of velocity fluctuation is less than zero: v(t)v(t + t) < 0, here t = 1/f s = 0.025 s. In order to estimate the value of dissipation rate we used Eq.   filtering of a signal. This value is comparable with the estimations performed using classic methods based on the power spectra and structure functions.

Simulation analysis and error estimates
Even if the local isotropy assumption of Kolmogorov (1941) is satisfied with a good accuracy, the TKE dissipation rate estimates are subject to errors that can result from a finite sampling frequency of a signal, a finite time window, sensor bias and noise. The last of those three causes was investigated in Sreenivasan et al. (1983), where it was shown that both the variance of the noise n 2 and the variance of its derivative ṅ 2 influence the measured number of crossings. A possible remedy was proposed by Poggi and Katul (2010), who suggested to use the threshold crossings, i.e., counting the number of times a signal crosses a given threshold T = 0, instead of the zero crossings in the case of signals with low signal-to-noise ratios. As for the signal considered in the previous section the signal-to-noise ratio becomes significant at higher frequencies (above 5 Hz), see Fig. 2, which are removed by the low-pass filter used in the proposed numberof-crossings method. We also applied the method of Poggi and Katul (2010); however, as it did not lead to any systematic change in our estimates, we further present results for the zero crossings only.
In order to quantify the errors resulting from the finite sampling frequency and finite time window and to test the performance of the proposed method, we performed the simulation analysis (Frehlich et al., 2001;Sharman et al., 2014). Results are compared with the standard spectral retrieval estimates without any additional corrections. However, we note in passing that spectral methods can be improved to account for the bias errors. The example is the maximum likelihood approach (Sharman et al., 2014) where, instead of the von Kármán model, a model power spectral density S model f is used, which takes into account the procedure for generating the empirical spectrum from discrete time series of finite length. Analogous approaches could also be formulated for the methods based on the number of crossings, which is a perspective for a further study.
To test the performance of the new proposals we generated a number of artificial velocity signals with frequency spectra and two point correlation functions prescribed by the von Kármán (1948) model. The equations resulting from the application of this model to the one-sided spectra considered in this paper are written below.
R 11 (r 1 e 1 ) ≈ 0.592548 u 2 r L 0 where K 1/3 is the modified Bessel function of order 1/3. Coefficients of the Fourier series expansion of velocity signal were calculated as where i = √ −1, a and b are random numbers from the standard Gaussian distribution with zero mean and unit variance, and W j = S(f j ) f , j = 1, . . ., N . Alternatively, the coefficients W j can be calculated as the discrete Fourier transform of R 11 , as described in Frehlich et al. (2001). The artificial velocity signal is finally constructed as the discrete inverse Fourier transform of w j (see Frehlich et al., 2001).
We used artificial signals with U = 55 m s −1 and the standard deviation u = 0.28 m s −1 . Those characteristics correspond to the ones of the signal considered in the previous Sect. 4.1. We set L 0 = 83.9 m in Eq. (29) to also obtain a comparable dissipation rate estimate = 2.5 × 10 −4 m 2 s −3 . Our first aim was to test how a finite sampling rate influences the number of crossings. For this purpose in each run we created an artificial signal of length N = 2 17 points and with the sampling frequency 200 Hz (5 times larger than the sampling of the signal considered in Sect. 4.1), which resulted in signal duration t ≈ 650 s. We treated this velocity series as a reference. Next, we took every fifth sample of this signal to create a 40 Hz velocity time series. We then calculated the number of crossings, as described in Sect. 4.1, and the power spectral density. We repeated the procedure 500 times and calculated the average of the obtained profiles (see Fig. 6). Due to the finite sampling frequency we observe the effect of aliasing -spectral densities for f higher than the Nyquist frequency are added to the spectral densities at f < 20 Hz. Distortions are visible for higher frequencies both in the power spectrum, Fig. 6a, and the N 2 i u 2 i profiles, Fig. 6b. We estimated the TKE dissipation rate from the averaged profiles using the method described in Sect. 4.1,Eq. (21), keeping the lower bound of the fitting range f 1 = 0.3 Hz constant and changing the upper bound f 2 from 1 to 19 Hz. Results are presented in Fig. 7 and compared with the corresponding PSD values, using the von Kármán model as the reference model spectrum.
We observe an increase in PSD estimates with increasing f 2 and a moderate increase in NCF over the input = 2.5 × 10 −4 m 2 s −3 , which shows that the number of crossings respond to finite sampling effects differently than the power spectrum. We note here that NCF values calculated from the averaged profiles of the 200 Hz reference signal (black line in Fig. 7) seem to be slightly overpredicted in comparison to the input , especially for smaller f 2 . A possible reason is the response of the filter used in the number-of-crossings method. We attempted to estimate this error using relation (18) between the number of crossings and the dissipation spectrum.
We first integrated f 2 S f of an unfiltered signal from f 0 to f i . Next, we calculated a spectrum S filtered f of a band-passfiltered signal taking f 1 and f i as the lower and upper bounds of a filter. We integrated f 2 S filtered f over the whole available range of f . The difference between the two integrals should represent a correction due to filter response. The values estimated from the corrected N 2 i u i 2 are presented in Fig. 7 as a black dot-dashed line. As it is seen that the estimations for the lower f 2 are improved as f 2 increases, NCF seems to be underpredicted. A possible reason for this might be that the filter influences the number-of-crossings statistics somewhat differently than the spectrum alone. Next, we tested the influence of the finite temporal window on the calculated statistics. We generated 1000 artificial signals, each time changing slightly the u value in Eq. (29), which led to a change in input (see Sharman et al., 2014); the value of L 0 remained unchanged. For each signal we es-  timated PSD from the standard power spectral density using the Welch overlapped segment averaging estimator implemented in Matlab ® with a 2 8 window and NCF from the number of crossings, Eq. (21). We performed these tests for the 40 Hz signals and the fitting range 0.3-5 Hz.
We first decreased the time window, taking each time only 1/8 of the created artificial signal for the analysis, which, in terms of L 0 from Eq. (29), resulted in the signal length L ≈ 50L 0 . Results of PSD and NCF estimates as functions of the corresponding input from the theoretical profile Eq. (29) are presented in Fig. 8 (upper plots). It can be seen that the bias error is larger for PSD ; however, the scatter of NCF is larger. The linear fits and the correlation coefficients are PSD = 0.9104 − 2.32 × 10 −5 , r = 0.9898, NCF = 0.9878 + 6.80 × 10 −5 , r = 0.9343.
Hence, for the signal length comparable to the lengths from the POST campaign we can expect a small underprediction of PSD estimates due to bias error and some overprediction due to aliasing (see Fig. 7). Both result in a small overprediction of PSD (Fig. 8, left column, lower plot). As far as NCF is concerned, the simulation analysis shows that it is less sensitive to the bias error (Fig. 8, right column); however, it has a larger scatter than PSD , at least for the generated artificial velocity fields. Results for the 40 Hz signal are slightly overpredicted (Fig. 8, right column, lower plot) due to aliasing and the fact that the number-of-crossings method gives somewhat larger estimates in this fitting range (see Fig. 7). Finally, we would like to address the issue of larger scatter observed for NCF . It follows from our study that the scat- ter in NCF depends on the value of filter cutoffs in the fitting range. In the final test, we set u = 0.28 m s −1 and input epsilon = 2.5 × 10 −4 m 2 s −3 and repeated the simulation 500 times for consecutively, 1/512, 1/265, 1/128, 1/65 of the original signal of 2 17 points. This corresponds to approximately 1L 0 , 2L 0 , 4L 0 , 6L 0 for L 0 in Eq. (29). We bandpass filtered the signal and consider small fitting range of 16-18 Hz. We normalized the results obtained by the input and calculated their standard deviations. Results presented in Fig. 9 show that at least for this case the standard deviations of + NCF is comparable with the standard deviation of +

Method based on missing spectrum recovery
The measurement signal used in Sect. 4.1 was also analyzed using the second method proposed in Sect. 3.2, Eqs. (27) and (28). We consider both formulas for the function f η , Eqs. (23) and (24). The advantage of the simpler, exponential formula (23) is that the one-dimensional spectrum function E 11 (k 1 ), Eq. (26), can be written in terms of the incomplete function as follows where (a, x) = ∞ x e −t t a−1 dt.
The correcting factor (27) If Eq. (24) is used as a model for f η , both integrals in Eq. (27) must be calculated numerically. On the other hand, as discussed in Pope (2000), Eq. (24) provides a better fit of experimental data in the dissipative range. With such preparation we applied the iterative procedure, as described in Sect. 3.2. In the POST experiment the effective cutoff frequency was estimated at f cut = 5 Hz, which corresponds to k cut = (2πf )/U = 0.57 m −1 . Using the sixth-order Butterworth filter, this resulted in u 2 N 2 cut = 0.0000719 · 1 s −2 for this signal. Accordingly we used Algorithm 1 with ν = 1.5 × 10 −5 m 2 s −1 and d η = 10 −6 m. We approximated the integrals in Eq. (35) using the trapezoid rule. The results of successive approximations of C F and converge fast to a fixed value, independently of the initial guess of = 0 (Fig. 10a). The increment dk 1 in Eq. (35) was approximated by k 1 = 5 · 10 −6 m −1 . For such choice we obtained NCR = 2.61 × 10 −4 m 2 s −3 , where the subscript NCR denotes variables calculated with the number-of-crossings method based on the recovered part of the spectrum. We used this as a reference value. In order to estimate the numerical accuracy of the proposed algorithm, we calculated the error = | − NCR | for different values of k 1 (see Fig. 10b). We obtain ∼ k 1.3 1 . Next we considered Eq. (24) as a model for f η and calculated the double integral in Eq. (27) using the trapezoid rule. We obtained the corresponding value NCR = 2.58 × 10 −4 m 2 s −3 , which is very close to the estimate from the simple exponential form Eq. (23) and (35).
It is worth noting that the proposed method is accounting for a dominant (and not directly measured) part of the spectrum based on the theoretical knowledge about its shape. This knowledge is simply reduced to the form of the correcting factor C F , Eq. (27), which contains the integral of k 2 1 E 11 (k 1 ). Fig. 11 illustrates the relation between the measured and the estimated part of the spectrum for the analyzed case with both forms of the function f η , Eqs. (23) and (24). The spectral cutoff of the data considered here (5 Hz) is in the inertial range, where k 2 1 E 11 (k 1 ) with both forms of f η functions is almost indistinguishable (see Fig. 11). At the same time integrals of the remaining (recovered) parts of k 2 1 E 11 (k 1 ) are almost equal, as independently of the choice of f η both dissipative spectra 2νk 2 E(k) must integrate to . As a result, for the given spectral cutoff, NCR estimates with the simple exponential Eq. (23) and (24) forms of f η are very close. This might change for larger cutoff frequencies. We expect that, in the case where the cutoff frequency is placed in a region influenced by the form of the f η function, the spectrum with Eq. (24) will provide better estimates of the TKE dissipation rate.
2.52 × 10 −4 and NCF = 2.54 × 10 −4 m 2 s −3 . The relative differences between those estimations are less than 5 %. We finally checked the estimates of the second method using synthetic signals as described in Sect. 4.2. For the cutoff f cut = 5 Hz, we generated 500 artificial signals of length L ≈ 400L 0 and with input = 2.5 × 10 −4 m 2 s −3 . We obtained the mean NCR = 2.55 × 10 −4 m 2 s −3 and a standard deviation equal to 9 % of the input value.

Broader overview of the methods' performance
Following the findings presented in the previous section both proposed methods were tested on much larger collection of data. For this purpose we used velocity signals also obtained during the POST research campaign. We have chosen horizontal segments at various levels within the boundary layer from flights TO10 and TO13. These flights were investigated in detail by Malinowski et al. (2013), due to the fact that they represent two thermodynamically and microphysically different types of stratocumulus-topped boundary layer. The dissipation rates of turbulent kinetic energy estimated from the standard structure function method SF and dissipation rates estimated from the modified zero-crossing methods NCF and NCR introduced in Sect. 3.1 and 3.2, respectively, are compared with the results obtained from the spectral method PSD in Fig. 12. The use of the simple exponential form of f η , Eq. (23) or (24) did not lead to any visible change in the results in Fig. 12. For flight 10 we obtained the following linear fits and the correlation coefficients r: The methods based on the signal zero crossings give comparable results to those resulting from standard methods, in spite of the fact that the second method is based on different physical arguments (assumes form of the whole spectrum, including the dissipative range of frequencies).
We believe that the there is a fair consistency in those results because one should take into account that the standard frequency spectra and structure function methods calculate approximate values of . Moreover, we have indicated in Sect. 2 that the constants C 1 and C 2 in Eqs. (4) and (5) are estimated with an accuracy of ±15 %.
In the present work we proposed two novel modifications of the zero-crossing method, such that it can be applied to moderate-resolution measurements. The turbulent kinetic energy dissipation rates obtained using the proposed methods show fair agreement with results of the standard power spectrum and structure function approaches.
We note that the standard structure function and power spectra methods are often used simultaneously, for better estimates (Chamecki and Dias, 2004), in spite of the same underlying physical arguments (second similarity hypothesis of Kolmogorov, 1941). Here, the proposed approach offers yet another option. Additionally, the second method with the spectrum recovery is based on different physical arguments, as it additionally makes use of Kolmogorov's first similarity hypothesis and a model for the dissipation range of the spectrum. Nevertheless, it can be used for signals with spectral cutoffs; hence, it offers an alternative to the spectral retrieval methods.
From the perspective of practical applications we can think of several possible advantages of the zero-crossing methods. First, the number of signal zero crossings can be calculated without difficulty and the proposed procedures are easy to implement. Other advantages follow from the results of the simulation analysis performed in Sect. 4.2. For the created artificial velocity signals, the NCF estimates responded differently to errors due to finite sampling or finite time windows than PSD . These differences in errors of the numberof-crossings and the power spectral method can make the former an additional tool to improve estimates from the atmospheric measurements. Here, a further, detailed study of bias assessment and removal is needed.
Moreover, we argue that the number-of-crossings method applied to the fully resolved signals has become a fairly standard tool for estimates, which are also used in the atmospheric measurements (see e.g., Poggi and Katul, 2010). Therein, the discussed advantages of the method are that no measurements of the signal gradients (to calculate the Taylor microscale) are required, no assumptions about scaling laws in structure functions (and power spectra) are needed and no simplifications in the TKE budget are adopted (for which is computed as a residual). The methods proposed in the current paper generalize the number-of-crossings method and makes it applicable also for signals with spectral cutoff. In the second approach a certain form of the energy spectrum must be assumed in order to calculate the correcting factor C F . Nevertheless, the proposed method can be interesting in particular for data with cutoffs reaching the dissipation range but still with part of this range missing (or contaminated with noise). In such cases, using only the inertial range estimates may lead to a significant loss of information, as the data from the dissipation range are not taken into account. Finally, we can deal with a situation when the recorded amplitude of certain frequencies is deteriorated due to measurement errors; nev-ertheless, the counted number of signal zero crossings could remain unaffected. In such cases the zero-crossing method could be advantageous over the power spectrum and structure function methods.
There are several perspectives for further work. First, the proposed methods could be tested for a wider range of signals (e.g., from Eulerian measurements within the boundary layer adopting the Taylor hypothesis), characterized by different resolutions and obtained under varying atmospheric conditions, to assess the scope of their applicability. Second, as far as the model spectrum is concerned, comparison with fully resolved experimental signals or direct numerical simulations data would be valuable to test different forms of the model spectra from Pope (2000) or Bershadskii (2016).