## Abstract

In this work, direct numerical simulation (DNS) of the stratocumulus cloud-top mixing layer is used to test various approaches to estimate the turbulence kinetic energy (TKE) dissipation rate *ε* from one-dimensional (1D) intersections that resemble experimental series. Results of these estimates are compared with “true” (DNS) values of *ε* in buoyant and inhomogeneous atmospheric flows. We focus on recently proposed methods of the TKE dissipation-rate retrievals based on zero crossings and recovering the missing part of the spectrum. These methods are tested on fully resolved turbulence fields and compared to standard retrievals from power spectra and structure functions. Anisotropy of turbulence due to buoyancy is shown to influence retrievals based on the vertical velocity component. TKE dissipation-rate estimates from the number of crossings correspond well to spectral estimates. The method based on the recovery of the missing part of the spectrum works best for Pope’s model of the dissipation spectrum and is sensitive to external intermittency. This allows for characterization of external intermittency by the Taylor-to-Liepmann scale ratio. Further improvements of this method are possible when the variance of the velocity derivative is used instead of the number of zero crossings per unit length. In conclusion, the new methods of TKE dissipation-rate retrieval from 1D series provide a valuable complement to standard approaches.

## 1. Introduction

Turbulence contributes to many atmospheric phenomena, including atmospheric convection and clouds. An important quantity that characterizes the smallest scales of such flows is the mean turbulence kinetic energy (TKE) dissipation rate *ε*. In formulating subgrid models for large-eddy simulation (Moeng and Sullivan 1994; Patton et al. 1998) or Lagrangian trajectory analysis of passive scalars (Poggi and Katul 2006), a robust estimation of the TKE dissipation-rate profile is needed.

Several methods have been proposed to calculate *ε* from one-dimensional (1D) velocity time series by making use of the local isotropy assumption. Indirect methods are based on the inertial-range arguments that follow from Kolmogorov’s hypotheses (Kolmogorov 1941; Albertson et al. 1997). Such methods are commonly used in the analysis of low- and moderate-resolution velocity time series of in situ airborne measurements (Sharman et al. 2014; Kopeć et al. 2016a). In the case of fully resolved velocity signals, the direct methods, based on measuring the mean variance of velocity fluctuation gradients, can be applied. Alternatively, Sreenivasan et al. (1983) proposed the zero-crossing approach, which requires counting the number of times per unit length the velocity signal crosses the zero threshold, denoted by , as shown in Fig. 1. The so-called Liepmann scale, defined as , is assumed to be equal to the transverse Taylor’s microscale , which is used to calculate *ε*. Since in signals with spectral cutoff are much smaller than in fully resolved signals, Wacławczyk et al. (2017) proposed two possible modifications to the zero-crossing method in order to estimate *ε* from moderate-resolution measurement data.

The first method is based on a successive filtering of the velocity signal, assuming that turbulence is homogeneous and isotropic and that the inertial scaling of −5/3 holds. In the second approach, an analytical model for the unresolved section of the spectrum is used to calculate a correcting factor to , so the actual relation between *ε* and can be used. Wacławczyk et al. (2017) validated these approaches on the data obtained during the Physics of Stratocumulus Top (POST) research campaign designed to investigate the marine stratocumulus clouds as well as the details of the vertical structure of the stratocumulus-topped boundary layer (Gerber et al. 2013; Malinowski et al. 2013).

The abovementioned methods for *ε* retrieval are based on the local isotropy assumption. However, this assumption might not always be fulfilled in real atmospheric conditions (Chamecki and Dias 2004; Jen-La Plante et al. 2016). Buoyancy is one reason for anisotropy in atmospheric flows. The energy spectra of buoyancy-driven turbulence has been studied by several authors (Bolgiano 1959; Lumley 1964; Lindborg 2006; Waite 2011; Kumar et al. 2014). First Bolgiano (1959) and Obukhov (1959) proposed the energy spectrum should scale as in stably stratified flows [referred to as Bolgiano–Obukhov (BO) scaling], where *k* is the wavenumber. Such scaling was later assumed to also hold in thermally driven flows; however, a fine-resolution simulation performed by Kumar et al. (2014) revealed turbulent convection exhibit the Kolmogorov spectrum. This was also confirmed by a direct numerical simulation (DNS) study of Rayleigh–Bénard convection (Verma et al. 2017) at a Prantl number . Here, *ν* and *κ* are, respectively, the kinematic and thermal diffusivities of a fluid. Kimura and Herring (2012) investigated homogeneous incompressible turbulence, subjected to a range of degrees of stratification, using the pseudospectral DNS method. The authors argued that due to the anisotropy of the flow a single mean dissipation rate cannot provide a universal Kolmogorov constant.

Physically complex atmospheric turbulence is not only inhomogeneous or buoyancy driven, but also includes the coexistence of laminar and turbulent regions called external intermittency (Lighthill 1956; Kurowski et al. 2009). The volume fraction occupied by a turbulent flow is called the intermittency factor *γ*. The motivation of this work is to investigate how the presence of anisotropy due to buoyancy and external intermittency affects the various retrieval techniques of *ε* in the atmospheric configurations, including the novel ones based on the number of crossings. Moreover, as the data used for the retrieval techniques are not idealized as DNS output, analysis of low-pass filtered velocity time series is undertaken, as measured by an artificial aircraft flying through the cloud, to assess performance of the methods. All the *ε* estimates are compared with actual *ε* values from DNS of the mixing layer at the stratocumulus cloud top. In spite of the inhomogeneity and physical complexity of the flow, the calculated *ε* profiles generally agree with DNS values within a certain degree of accuracy. The observed deviations follow from the physical complexity of the flow and low Reynolds number (Re) of the DNS as compared to real atmospheric conditions. The latter issue makes the spectral retrieval methods difficult due to the relatively short inertial range. Further, an additional source of errors includes the deviations of the Taylor-to-Liepmann scale ratio from unity, as the assumption lies behind the number of crossing method (Sreenivasan et al. 1983). In this work, we show that in case of externally intermittent flows.

Due to the abovementioned difficulties, the present work focuses on the second method proposed in Wacławczyk et al. (2017), based on an analytical model to resolving the missing part of the spectrum. We propose its alternative form, replacing the Liepmann scale with the Taylor microscale. Results obtained with this new approach compare favorably with the DNS over a wide range of cutoff wavenumbers.

This paper is structured as follows: In section 2 we describe the current state of knowledge and propose modification of the iterative method. The setup of the case study used in our analysis is explained in section 3. In sections 4 and 5 results of the TKE dissipation-rate retrieval are presented. Section 6 provides conclusions of the analysis results.

## 2. TKE dissipation-rate estimates from 1D signals

### a. Direct and indirect methods

The TKE dissipation rate is defined as (e.g., Pope 2000)

where is the fluctuating strain rate tensor, denotes the *i*th component of fluctuating velocity, and is the ensemble average operator. The exact definition cannot be used to estimate *ε* in case only 1D intersections of turbulent velocity field are available from experiments. Additionally, the resolution of the measured signals can be deteriorated due to finite sampling frequency of a sensor, as well as measurement errors.

The methods used to retrieve the TKE dissipation rate from 1D signals can be divided into two categories: direct and indirect. In direct methods the gradients of velocity are measured. Indirect methods relate the small-scale phenomenon of dissipation with inertial-range scales, as predicted by Kolmogorov’s second hypothesis (Kolmogorov 1941). Additionally, all methods are based on the local isotropy assumption (Kolmogorov 1941).

The two most common indirect approaches use an inertial-range scaling form of the power spectra and structure functions. In the homogeneous and isotropic turbulence, the following energy spectrum is assumed (Pope 2000):

where the constant is derived from experimental data and and are nondimensional functions. The term should be equal to the TKE dissipation rate *ε* if the second similarity hypothesis is satisfied.

Functions and specify the shape of the energy spectrum in the energy-containing range and the dissipation range respectively; *L* is the length scale of large eddies and is the Kolmogorov length scale, which is connected with the dissipative scales. The function tends to unity for small while tends to unity for large , such that in the inertial range, the formula is recovered. Pope (2000) proposed the following forms of functions and :

where is a positive constant, , and

where and . If , Eq. (4) reduces to the exponential spectrum:

where . An alternative model spectrum for is the Pao spectrum defined as

where . One-dimensional longitudinal and transverse energy spectra and , respectively, are related to the energy spectrum function (Pope 2000):

In the inertial range, the spectra follow Kolmogorov’s −5/3 law:

where and . Equations (8) make it possible to estimate the TKE dissipation rate from the inertial-range profile of the one-dimensional energy spectra.

Alternatively, the profiles of the second- and third-order longitudinal structure functions can be used to calculate *ε*. The *n*th-order structure function reads . Here is the longitudinal component of the velocity fluctuation vector. In the inertial subrange, the second- and third-order structure functions are related to the dissipation rate by the formulas (e.g., Pope 2000)

where and should approximate *ε*.

In the case of a turbulent signal resolving the smallest scales, the “direct” relation between the TKE dissipation rate and the longitudinal, or transverse, Taylor microscale can be used:

The longitudinal Taylor microscale equals

and the transverse microscale is . In case of isotropy, coincides with *ε*.

Other direct methods for calculating TKE dissipation rate, based on number of zero crossings, have been proposed by Sreenivasan et al. (1983) and are used by many authors (see, e.g., Poggi and Katul 2009, 2010; Wilson 1995; Yee et al. 1995). The zero-crossings method was first introduced by Rice (1945). It assumes that the stochastic process *q* and its derivative have Gaussian statistics and are statistically independent. Following this, the square of the number of zero crossings per unit time is

Sreenivasan et al. (1983) considered the Liepmann scale defined as

where is a number of zero crossings of per unit length. Using Eq. (12) Sreenivasan et al. (1983) assumed that ; hence

For this, it was argued that Eq. (12) also holds for strongly non-Gaussian velocity signals (or for non-Gaussian derivative of the time series). This implies that strong departures from Gaussianity do not necessarily yield values appreciably different from unity for the ratio of .

Based on this result and Eq. (10), Sreenivasan et al. (1983) proposed a formula for calculating *ε*, applicable to fully resolved signals (measured down to the smallest dissipative eddies), which reads

or for the number of crossings calculated from the time series of transverse velocity fluctuation component :

The number of crossings is related to the energy spectra by the formula

In case a signal is low-pass filtered the number of zero crossings per unit length depends on cutoff wavenumber. Hence, Wacławczyk et al. (2017) proposed a possible modification for zero-crossing method to retrieve *ε* from the restricted range of *k* values. The motivation was to increase robustness of *ε* retrieval using different statistics. Two procedures formulated in Wacławczyk et al. (2017) are discussed in more detail in section 2b, below.

### b. Methods based on number of crossings

If we assume that the applied filter is rectangular in the wavenumber space, then from Eqs. (17) and (8) the TKE dissipation rate can be estimated from

where is the variance and is the number of crossings per unit length of a signal filtered with a cutoff wavenumber which is inside the inertial range. Filtering the signal with a series of cutoff wavenumbers , can be estimated from Eq. (18) using a linear least squares fitting method, and used as a proxy for the TKE dissipation rate *ε*.

We note in passing that the scaling of with was also investigated by Mazellier and Vassilicos (2008) to estimate the dissipation-rate constant in Taylor’s formula , where is the longitudinal integral length scale of turbulence.

The second method is based on recovering the missing part of the spectrum in the inertial and dissipative range, by introducing a correcting factor to the number of crossings per unit length. As such, this method can be treated as a smooth blending between indirect and direct methods as it recovers the former as the filter cutoff moves into the inertial range and the latter as the filter cutoff moves into the dissipative range.

The number of crossings per unit length *N*_{cut} is calculated from the low-pass filtered signal, where the finescale fluctuations have the highest wavenumber *k*_{cut}. Assuming again that the filter is rectangular in the wavenumber space, Eq. (17) written for the filtered signal is

where is the correcting factor. Assuming the energy spectrum can be described by Eq. (2) with and using relation (7) between and we obtain

where , , and is the Kolmogorov length. Introducing Eq. (21) into Eq. (20) and changing the variables to and , the correcting factor is obtained as

To calculate from Eq. (22) a value of *η* should first be specified; hence, an iterative procedure was proposed in Wacławczyk et al. (2017). It starts with an initial guess of the TKE dissipation rate . With this, the corresponding value of the Kolmogorov length is calculated and introduced into Eq. (22) for . The TKE dissipation rate after the first iteration is found from Eq. (23). The procedure can be repeated; that is, the next approximation of can be calculated and substituted into Eq. (22). After several iterations the procedure converges to the final value of that should approximate the TKE dissipation rate *ε* with an error defined by a prescribed form , where is a given error value.

In this method, the cutoff *k*_{cut} may be placed in the inertial or dissipative range. In the latter case, the spectral retrieval methods may lead to loss of certain information as they are based on the inertial-range scaling only. In Wacławczyk et al. (2017), performance of the new methods was tested on measurement data obtained during the POST airborne research campaign (Gerber et al. 2013; Malinowski et al. 2013) with the cutoff placed well in the inertial range. It was shown that estimates obtained with the new methods were comparable with results of standard retrieval techniques; however, differing responses to errors due to finite sampling and finite averaging windows were observed. Hence, the new methods can complement the standard techniques to increase robustness of *ε* retrieval.

### c. Alternative formulation of the iterative method

Estimates of from Eqs. (22) and (23) may be deteriorated if the ratio of Taylor’s microscale to the number of crossings microscale deviate from unity [see Eq. (14)]. For this reason, we propose a different formulation of this method.

For low-pass-filtered signals (and filters rectangular in the wavenumber space), Eq. (24) becomes

Hence, is related to by the formula

If we introduce Eq. (21) for into Eq. (26), we obtain the same correcting factor as in Eq. (22). Based on Eqs. (12), (23), and (26), the value of dissipation rate is

An iterative procedure, similar to the one described in section 2b will be used to calculate . With a first guess of , the correcting factor will be calculated from Eq. (22) and introduced into Eq. (27) to calculate the new value of . The procedure can be continued until the condition is satisfied.

In this work, we will investigate and compare the performance of both approaches from Wacławczyk et al. (2017) and the new Eq. (27) with different model assumptions for as written in Eqs. (4) and (6) with DNS. As the DNS data contain complete information about turbulence, it will be possible to assess how the *ε* estimates change with changing cutoff wavenumber.

## 3. Stratocumulus cloud-top mixing-layer simulation for DYCOMS II RF01 case

As a test case, we consider a cloud-top mixing layer. This system mimics the cloud-top region of stratocumulus clouds and proves convenient in studying some aspects associated with submeter scales, like evaporative cooling, as simulations of the complete boundary layer cannot reach these small grid spacings (Mellado et al. 2010; Mellado 2017; Mellado et al. 2018). The system consists of two horizontal layers of moist air: an upper region, which is warm and unsaturated and represents the free troposphere, and a lower region, which is cool and saturated and represents the cloud. In-cloud turbulence and the vertical wind shear across the cloud top creates the cloud-top mixing layer that is illustrated in Fig. 2. In-cloud turbulence is driven by the longwave radiative cooling of the cloud top and by the evaporative cooling caused by the mixing of cloudy and tropospheric air. Radiative cooling is characterized by the net upward radiative flux and the radiative extinction length , over which that cooling concentrates. We consider the first research flight of the DYCOMS II field campaign as reference, and we use the measurement-based estimates and (Stevens et al. 2003). The radiative properties imply a reference buoyancy flux , where *g* is the gravitational acceleration, and a reference velocity scale . The Reynolds number in the simulation is , which is about 300 times smaller than that in the atmosphere. The velocity variation across the cloud-top region is 3 m s^{−2} (Faloona et al. 2005).

The horizontal size of the computational domain is . The domain is discretized with 5120 × 5120 × 2048 points in the streamwise, spanwise, and vertical directions, which assures fine resolution of the flow down to the smallest dissipative eddies with a characteristic size . The system is statistically homogeneous over the horizontal planes; the data along these planes are used to construct the different statistics, which depend on the vertical coordinate *z* and the time *t*. Further details of the simulations can be found in Schulz and Mellado (2018). We investigated the three velocity components in four horizontal planes at heights , where corresponds to the height of minimum buoyancy flux and corresponds to the height of maximum buoyancy flux.

Figure 3 includes vertical profiles of the mean velocity; in the streamwise (*x*) direction; the root-mean-square (rms) of the three velocity components , , and in the streamwise, spanwise, and vertical directions, respectively; and the budget of the turbulence kinetic energy. The angle brackets indicate horizontal average. The viscous dissipation rate of the TKE is defined by Eq. (1), the buoyancy flux is , and the shear production term is . It is worth noting that profiles of rms velocity component fluctuations in Fig. 3 indicating anisotropy of turbulence are in agreement with the measurement data reported in Jen-La Plante et al. (2016). Moreover, profiles of budget terms of turbulence kinetic energy are consistent with observations reported in Brost et al. (1982) and results of high-resolution large-eddy simulations (LESs) (Kopeć et al. 2016b; Heinze et al. 2015).

## 4. TKE dissipation-rate estimates from inertial-range scaling

### a. DNS signals

The methods related to the local isotropy assumption and inertial-range scaling are commonly used to analyze 1D signals from airborne measurements. At the same time, turbulent flows in clouds or atmospheric boundary layers are in fact inhomogeneous and buoyant. The purpose of this analysis is to check how predictions of these methods, when applied to DNS data, compare with the true value of calculated from Eq. (1). All analyses were done with MATLAB software.

We first investigated 1D spectra of three velocity components *u*, *υ*, and *w* (see Figs. 4–6, respectively). To calculate the compensated spectra, we multiplied each , , and by a corresponding at and by or . Spectra calculated in the *x* direction were additionally averaged in the spanwise (*y*) direction. Similarly, spectra calculated in the *y* direction were averaged in the streamwise (*x*) direction. The horizontal lines in Figs. 4–6 are equal to the constant coefficients from Eq. (8), for the longitudinal and for the transverse spectra.

We observe similar profiles of corresponding compensated spectra at planes , , . Results calculated at , placed in the upper part of stratocumulus cloud, are clearly different. This region of the flow is affected by the presence of shear and stable stratification (see Fig. 3) as well as external intermittency. The contribution of large-scale instabilities induced by the shear is visible in the profile of at plane as a maximum at small *k*. Moreover, differences are observed between different types of spectra. The longitudinal ones, and (Figs. 4a and 5b), seem to closely follow Kolmogorov’s K41 theory in a certain range of wavenumbers. This is in spite of the relatively low Re of the considered flow, where a clear separation between the dissipative and energy-containing scales may not be attained. At the same time, the transverse spectra and (Figs. 4b and 5a) seem to scale with or where *a* and *b* are somewhat smaller than 5/3. Interestingly, inertial range with the scaling close to and can be distinguished for the transverse spectra of the third velocity component in Fig. 6; however, the constant is larger than the value 0.65 at planes , , and . Kaiser and Fedorovich (1998) argued that the excess of over value typical for the isotropic turbulence could be caused by the presence of dominating buoyant forcing, which favors vertical motions. In such a case, pressure fluctuations were insufficient to isotropize turbulence. Spectral anisotropy of velocity component structure in the layer indicates stable stratification above the cloud top, and is in agreement with the experimental data and LES reported in Pedersen et al. (2018).

In the following, we investigate to what extent deviations from the K41 theory observed in Figs. 4–6 affect estimations of *ε*. To estimate from power spectra [Eqs. (8)] we fit a line with −5/3 slope on a logarithmic plot. The log of the intercept is equivalent to or . Figure 7 provides the scaling of with filter cutoff and [see Eq. (18)]. For each , a corresponding of a filtered signal was calculated. The linear fit slope is equivalent to , out of which the dissipation rate was calculated. We used the sixth-order Butterworth filter to calculate successive , ensuring there was no difference in results between the fifth- and the sixth-order filters, which was also reported in Mazellier and Vassilicos (2008). The frequency response characteristic of the filter was investigated in Wacławczyk et al. (2017) on artificial velocity time series. In this case, the filter led to small (~4%) overprediction of . In this work we neglected the effect of the filter on *ε* estimates.

Next, we estimate and from the second- and third-order structure functions as written in Eq. (9). We obtained inertial-range values by fitting linearly a slope of 2/3 to the second-order structure function as shown in Fig. 8. An analogous procedure is performed to calculate .

In isotropic turbulence, all estimates of *ε* should, theoretically, be equal. It would hence seem appropriate to use the same fitting ranges for , , and , appropriately converted from *k* space to *r* space (i.e., with ). In practice, the fitting was difficult due to the short inertial ranges—an attribute of low-Re flows. Although inertial ranges of atmospheric high-Re flows cover a few decades of *k* numbers, finite averaging windows and finite resolutions of the measurements cause analogous problems—results are in fact dependent on the chosen fitting ranges. As pointed out in Hou et al. (1998) and confirmed by other authors (Chamecki and Dias 2004), the effect of finite power-law range in the spectral space results in a much shorter power-law range in the physical space. Also, in our case, the fitting ranges optimal for were different than those for power spectra. To show differences in the estimates, we compare results of different methods and different fitting ranges in detail. We consider planes and , as they are regions of maximum and minimum buoyancy flux, respectively.

Table 1 presents results for the plane which is placed in the turbulent region inside the cloud. The true value of *ε* at this plane equals . The first fitting ranges presented in Table 1 seemed to be optimal for the investigated spectra; the second seemed to be optimal for the second-order structure functions . We present dissipation-rate estimates from the power spectrum , number of crossings , and second- and third-order structure functions and . We observed a certain discrepancy of results, also between and that are standard methods of estimating TKE dissipation rate. Moreover, are overpredicted and , , and underpredicted when compared to .

The TKE dissipation-rate estimates from 1D signals, , , and at planes , are compared in Fig. 9a. The fitting ranges were chosen to match the inertial range of structure functions. It is seen that are in most cases smaller than . The linear fit for the presented data is

An underprediction of versus was also observed by Jen-La Plante et al. (2016) in the POST measurements of stratocumulus clouds in the well-mixed cloud-top layer (therein called CTL). Results from POST for this part of the cloud are presented in Fig. 9b and the linear fit is

Having exact DNS data at hand we can observe that the true dissipation rate can differ both from and (see Fig. 9a). Here, the linear fit line is flat:

Table 2 shows the corresponding fits for the horizontal profile placed in the upper part of the stratocumulus cloud. Here, the discrepancies between and estimates from 1D intersections of the velocity field are larger. Results differ also between the horizontal velocity components (*u* in *x*, *u* in *y*, *υ* in *x*, and *υ* in *y*) which makes the local isotropy assumption questionable.

As it is seen in Tables 1 and 2, the estimates of *ε* from the vertical velocity component *w* differ from those based on horizontal components. They are overpredicted in comparison to at plane that is placed inside the CTL, where buoyancy is a source of turbulence generation (see Table 1) and greatly underpredicted at plane (see Table 2) where the negative buoyancy damps the vertical velocity fluctuations. In the analysis of POST data by Jen-La Plante et al. (2016), the layer of the cloud was referred to as moist and sheared cloud-top mixing sublayer (CTMSL), and likewise, *ε* estimates based on *w* were underpredicted in this region as compared to estimates from *u* (see Fig. 7 therein). In the subsequent sections of this work, we will estimate *ε* based on the four signals with horizontal velocity components, in order to obtain results closer to .

We compare *ε* estimates using the two different fitting ranges from Tables 1 and 2 (averaged over the four signals *u* in *x*, *u* in *y*, *υ* in *x*, and *υ* in *y*) in Figs. 10a and 10b. The structure function’s fitting ranges give better results. We can observe that calculated using Eq. (18) agrees closely with at . Our results of from the third-order function are underestimates [also reported by Chamecki and Dias (2004)], which is contrary to the idea that these estimates are preferable due to their analytically derived constant (4/5). The maximum value of the TKE dissipation rate was found in the cloud-top mixing sublayer , which agrees with the experiment reported in Jen-La Plante et al. (2016). However, it is seen in Fig. 10a that all estimates are underpredicted in comparison to in this nonisotropic, shear-influenced part of the cloud.

### b. Moderate- and low-resolution signals

Signals available from in situ airborne measurements are far from the idealized fully resolved DNS data. Finite sampling frequency of a sensor and measurement errors induce effective spectral cutoff of velocity time series. To investigate the influence of the finite sampling on the TKE dissipation-rate estimates, we perform the following tests of DNS data. We consider a virtual aircraft that measures velocity signal with effective cutoff wavenumbers, , placed, approximately, within or close to the inertial range. For each , if the aircraft flies in the streamwise *x* direction we create a new path every 100 grid points in the *y* direction, such that 52 signals are collected. Similarly, if the aircraft flies in the *y* direction, we create 52 paths in the *x* direction. In the first test, we average the obtained power spectra and profiles and calculate TKE dissipation rates and . Finite sampling frequency causes aliasing; that is, spectral densities for *k* higher than the wavenumber are added to the spectral densities at . This causes a bias error of the TKE dissipation-rate estimates. As seen in Fig. 11 the bias of is smaller than the bias of ; in particular, results of for are still close to . This result is in line with previous error analyses on artificially generated velocity time series by Wacławczyk et al. (2017), where a smaller bias error but somewhat larger scatter of was observed in comparison to . In the present analysis, we do not introduce any additional corrections to the power spectra or to the profiles, to reduce the bias (such methods can be formulated for high-Re flows; see Sharman et al. 2014).

Next, in order to test scatter of the results we estimate and from each velocity signal separately. For a given we used the same fitting range for all signals. Increasing we extended the lower bound of the fitting range; that is, for *k*_{cut} = 0.62, 1.25, 2.5, and 5 m^{−1}, the fitting ranges were, respectively, *k* = [0.3, 0.6], [0.3, 0.8], [0.3, 1.0], and [0.3, 1.2] m^{−1}. Figure 12 presents values of versus calculated at the plane for *k*_{cut} = 5 and 0.62 m^{−1}. In Fig. 12a, we observe somewhat larger scatter of ; however, as decreases, scatter of both and becomes comparable (see Fig. 12b).

An additional method makes it possible to decrease statistical uncertainties of the TKE dissipation-rate estimates. Figure 13 presents standard errors of and separately and the error of the mean calculated from a twice larger sample that contains both and . The estimates were done for the 2 × 52 1D intersections (in the *x* and *y* direction) of the velocity field measured by the virtual aircraft along the plane . We assume that the samples are uncorrelated; hence, the standard error equals , where std is the standard deviation and *N* is the size of the sample. As seen in Fig. 13, , calculated from the twice larger sample containing both and , is smaller than either or individually.

The obtained results confirm the method based on the number of crossings responds differently to errors due to finite sampling than the spectral retrieval technique. It can also complement standard approaches to reduce the standard error of the mean TKE dissipation rate.

## 5. TKE dissipation-rate estimation with the direct and iterative methods

### a. Direct methods

Results discussed in section 4 reveal TKE dissipation-rate recovery based on inertial-range arguments is difficult in the considered flow case. The first source of error relates to the relatively low Re of DNS simulations. The available DNS data allow the estimation of *ε* from the direct methods. We calculated from the Sreenivasan–Rice Eqs. (15) and (16) (Sreenivasan et al. 1983) and compared it with from Eq. (10). Results for planes and are given in Tables 1 and 2. At plane , the discrepancy between and is caused by the values that deviate from unity. Possible reasons for this could be the low Re number of the considered flow, and strong non-Gaussianity of the pdfs of velocity derivatives. The estimates of from *u* and *υ* velocity components given in Table 1 are close to , while estimates from the vertical component are overpredicted. This shows that, even in the core region of the model cloud, the local isotropy assumption is not satisfied.

At the plane we observe considerable discrepancies between and and large values of (see Table 2). Large were also reported by Kailasnath and Sreenivasan (1993) in the upper part of the boundary layer, affected by the external intermittency. This leads us to the idea that is an indicator of external intermittency, as we will describe in more detail in section 5c. Values of and averaged over the four horizontal signals *u* in *x*, *u* in *y*, *υ* in *x*, and *υ* in *y* are additionally compared with in Fig. 14. As before, the largest difference is observed at plane where the given 1D intersections of the velocity field are clearly not sufficient to estimate the TKE dissipation rate.

### b. Two formulations of the iterative method

In this section, we consider the second, iterative approach from Wacławczyk et al. (2017), described in section 2b, where the cutoff can be moved toward the dissipative part of the spectrum. We test different models for the function in Eq. (22), be it the Pope [Eq. (4)], Pao [Eq. (6)], or exponential model [Eq. (5)]. Moreover, we discuss results of both formulations of the iterative method, the one based on number of crossings [Eq. (23)], as proposed originally in Wacławczyk et al. (2017), and the new, alternative form based on the variance of velocity derivative [see Eq. (27)].

Figure 15 presents model spectra of *u* in *x* for the horizontal plane . As it is seen, the Pope formulation [Eq. (4)] provides a much better fit with the DNS spectra than the Pao [Eq. (6)] or exponential models [Eq. (5)].

Next, we investigate a DNS signal that is first low-pass filtered with the use of a sixth-order Butterworth filter with a given . According to the procedure described in Wacławczyk et al. (2017), in order to estimate *ε* in the iterative method, a first guess for the Kolmogorov length is made. We take , however, independently of the initial guess the procedure always converges to the same value of or . We calculate the correcting factor from Eq. (22), and next, the value of dissipation rate is estimated with Eq. (23) or (27). We approximate the integrals in Eq. (22) with the trapezoidal rule. We repeat the procedure, as described in section 2b, until the condition with is satisfied. Convergence is reached in all simulations before the tenth iteration.

Figure 16 shows the difference between calculated from Eq. (23) and estimated using Eq. (27). The latter compares more favorably with the DNS results over a wide range of numbers. Again, the best agreement is observed for the Pope model spectrum [Eq. (4)], which is the most favorable choice for the iterative method.

Figure 14 shows and calculated according to the Pope model for , which is within the dissipative range, as a function of vertical coordinate . The difference in results between both formulations can be explained by the fact that the Rice formula [Eq. (12)] is only approximately satisfied for the considered signals. Let us define the length scale of the filtered signal, analogous to the Taylor microscale [Eq. (11)] :

Analogously, will denote the Liepmann scale calculated for the filtered signal. We note here that in case of the airborne measurements of high-Re turbulence with cutoff wavenumbers placed in the inertial range, investigated in Wacławczyk et al. (2017), the condition was satisfied with a good accuracy. As far as the present DNS data are concerned, changes with (results not shown here). It is closer to 1 if is placed in the inertial range but increases with increasing toward values presented in Tables 1 and 2. This fact is the source of existing discrepancies between and , and as seen in Fig. 16b, compares markedly better with .

### c. ratio as the intermittency measure

The motivation of the present subsection is to understand the reason for the strong deviations of the ratio from unity. The data seem to suggest that this deviation is caused by the external or global intermittency connected with the existence of laminar spots within the turbulent flow.

In the literature, several different methods were proposed to differentiate between rotational (turbulent) and irrotational (nonturbulent) parts of a measured velocity signal (Zhang et al. 1996). Each requires definition of an indicator function *q*, a criterion function , and a threshold level . The flow is assumed turbulent when . If instantaneous values of vorticity are known from measurements or DNS data, the enstrophy can be used as the criterion function *q*. However, *γ* estimation based on vorticity may also be subject to error due to the presence of mean gradients or nonturbulent wavelike motions that spuriously increase above the threshold (Ansorge and Mellado 2016).

As a first approximation, we assume that in the externally intermittent flow, the statistics will change to and . Moreover, the laminar part of the signal does not significantly contribute to the number of crossings; hence, we will detect crossings per unit length in the intermittent signal. With this, the Taylor microscale, the Liepmann scale, and their ratio will change to

where the subscripts *I* are related to the statistics in the intermittent flow. If , then in the intermittent flow, . We note, however, that in our case is around 0.8 even in the core region of the flow. We will compare predictions of Eq. (32), that is, the ratio

(the subscript *T* is used to denote mean value in the turbulent, core region of the flow), with *γ*. For this purpose, we first subtracted the mean from the instantaneous vorticity. The enstrophy based on the fluctuating vorticity was our criterion function. Regions where was smaller than a certain threshold value were identified as “laminar spots.” We calculated *γ*, as the mean volume fraction of turbulent flow at a given vertical height, by averaging in the streamwise direction and, additionally, in the spanwise direction over four planes *x*/*L*_{0} = 0, 13.5, 27, and 40.5. In Fig. 17 this profile is compared with the calculated ratio from Eq. (33). We utilized .

We observe favorable agreement between both curves, at least for larger *γ* values. Discrepancies for small *γ* are due to numerical errors, as both and are small in these regions. The results suggest that the Liepmann-to-Taylor scale ratio, calculated from 1D intersections of the velocity field is a good indicator of external intermittency. In other words, discrepancies between and reported in section 5b reveal the presence of external intermittency in a given flow region.

## 6. Conclusions

In this work, we focus on scaling of the energy spectra of turbulent flows in stratocumulus clouds and investigate different methods of TKE dissipation-rate retrieval from 1D intersections of the flow domain. We investigate data from numerical experiments in the stratocumulus cloud-top mixing-layer simulations. In such experiments, high Re observed in nature could not be reached; however, we argue model assumptions can still be tested, enabling conclusions applicable to “real world” flows to be drawn. Finite sampling frequency of a sensor and measurement errors deteriorate results of airborne experiments. Comparison with high-resolution numerical simulations might help to estimate the role of resulting effective cutoff frequencies and aliasing.

The investigated flow case appeared largely influenced by buoyancy effects that cause deviations from the Kolmogorov scaling. This, in turn, results in errors of the TKE dissipation-rate retrieval based on local isotropy assumption. We found the longitudinal spectra of horizontal velocity components and are comparable to Kolmogorov scaling over a certain range of wavenumbers, unlike the transverse spectra. The 1D spectra of the vertical component show −5/3 scaling range; however, the constant is larger than the isotropic value. As a result, TKE dissipation-rate estimates from *u*, *υ*, and *w* velocity components differ, which withstands the local isotropy assumption. We also show that estimates in the upper section of the cloud are subject to large errors, as the buoyancy flux is minimum and stable stratification strongly hinders vertical motions.

In this work, we investigated different methods of TKE dissipation-rate retrieval, including the two approaches based on the number of crossings per length proposed in Wacławczyk et al. (2017). The first method used the inertial-range arguments and provided scaling of in this range. From results presented in section 4, we can conclude the performance of this approach is comparable with standard spectral retrieval methods. Moreover, we investigated velocity signals with effective spectral cutoffs measured by a “virtual aircraft” flying through the stratocumulus cloud. We showed that estimated from the number of crossings had smaller bias error than calculated from energy spectra. On the other hand, standard deviations of results were larger than that of for two higher cutoffs. Still, an additional method of the TKE dissipation-rate retrieval makes it possible to reduce the standard error of a mean estimated from a finite-size sample.

The second method proposed in Wacławczyk et al. (2017) was based on the recovery of the missing part of the spectrum, that is, the part with *k* higher than the cutoff wavenumber . It is based on a model for the inertial and dissipative parts of the spectrum. Hence, it could be used for signals with placed in the dissipative range. We showed that Pope’s model for the dissipative part of the spectrum provides the best fit to the DNS data. As moves toward the high-wavenumber part of the spectrum, estimated deteriorated. As identified, the discrepancies follow from the deviations of the Taylor-to-Liepmann scale from unity. We also showed that this ratio could be used as a certain indicator of the external intermittency.

We proposed an alternative formulation of the second method, where the variance of velocity derivative is used instead of the number of crossings per length. The remaining procedure is consistent; that is, the correction factor for the missing part of the spectrum and *ε* are calculated iteratively. Results compare very favorably with the DNS data. This also suggests that the dissipative part of the spectrum has a universal form with a prescribed dependence on *ε*.

This study revealed that novel methods for TKE dissipation-rate retrieval can complement standard approaches. A perspective for a further study is to test their performance on a larger set of experimental data.

## Acknowledgments

This work received funding from the EU Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie Actions, Grant Agreement 675675. MW and SPM acknowledge matching funds from the Polish Ministry of Science and Higher Education 341832/PnH/2016. The authors acknowledge the proofreading of the manuscript by Kristin Goździkowska.

## REFERENCES

*Turbulent Flows.*Cambridge University Press, 771 pp., https://doi.org/10.1017/CBO9780511840531.

*Concentration Fluctuations and Averaging Time in Vapor Clouds.*American Institute of Chemical Engineers, 181 pp., https://doi.org/10.1002/9780470937976.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).