Assessment of the Precision of Spectral Model Turbulence Analysis Techniques Using Direct Numerical Simulation Data

The spectral model turbulence analysis technique is widely used to derive kinetic energy dissipation rates of turbulent structures (ɛ) from different in situ measurements in the Earth's atmosphere. The essence of this method is to fit a model spectrum to measured spectra of velocity or scalar quantity fluctuations and thereby to derive ɛ only from wavenumber dependence of turbulence spectra. Owing to the simplicity of spectral model of Heisenberg (1948), https://doi.org/10.1007/bf01668899 its application dominates in the literature. Making use of direct numerical simulations which are able to resolve turbulence spectra down to the smallest scales in dissipation range, we advance the spectral model technique by quantifying uncertainties for two spectral models, the Heisenberg (1948), https://doi.org/10.1007/bf01668899 and the Tatarskii (1971) model, depending on (a) resolution of measurements, (b) stage of turbulence evolution, (c) model used. We show that the model of Tatarskii (1971) can yield more accurate results and reveals higher sensitivity to the lowest ɛ‐values. This study shows that the spectral model technique can reliably derive ɛ if measured spectra only resolve half‐decade of power change within the viscous (viscous‐convective) subrange. In summary, we give some practical recommendations on how to derive the most precise and detailed turbulence dissipation field from in situ measurements depending on their quality. We also supply program code of the spectral models used in this study in Python, IDL, and Matlab.


Plain Language Summary
Turbulence plays a central role in most geophysical fluids, but our understanding of it remains limited. Atmospheric turbulence plays roles as diverse as dispersion of pollutants in the boundary layer to strong influences on the thermal and wind fields on global scales from the surface to above 100 km. It also is key to transports in, and the large-scale circulation and structure of, Earth's oceans. Measurements quantifying turbulence intensities and their environments are key to understanding its many effects but remain challenging. In situ measurements of various quantities enable estimates of turbulence intensities but must be calibrated to be of optimal benefit. This study employs a direct numerical simulation of Kelvin-Helmholtz instabilities that quantifies the associated turbulence dynamics exactly over the range of scales simulated to evaluate theoretical spectral forms enabling the best estimates of the known turbulence intensities.
In the last decades, direct numeric simulations (DNS) were successfully used to characterize the structure, dynamics, and anisotropy of turbulence (e.g., Fritts et al., 2003Fritts et al., , 2006. Early DNS studies only captured limited inertial range turbulence dynamics, nevertheless enabled an assessment of the vorticity dynamics driving the turbulence cascade Arendt, 1998;Arendt et al., 1997;Fritts et al., 1998). The next-generation DNS already allowed for simulations of turbulence field down to fine scales within dissipation range with sufficient details (e.g., Fritts et al., 2009aFritts et al., , 2009b. The highly resolved velocity field produced by such DNS allows for a precise and detailed derivation of kinetic energy dissipation rate with a spatial resolution close to those achieved for the velocity field. Since also scalar fields of potential temperature fluctuations are calculated in these DNS, it is possible to relate them to the dissipation of the kinetic energy. Results of high-resolution DNS of gravity wave (GW) instabilities and the produced volumetric data are shown and discussed in detail in for example, Fritts et al. (2009aFritts et al. ( , 2009b. More recent DNS by  and  studied multiscale dynamics accompanying GW instability arising as a result of GW-fine structure (GW-FS) interactions. These simulations enlightened differences in morphologies of dissipation fields at different stages of evolution accompanying different types of interactions. Such simulations reproduced the fine structure of the velocity and dissipation fields and their evolution in time and were successfully used to explain observations of mesosphere/lower thermosphere (MLT) dynamics (Fritts et al., 2017).
Our goal in this work is to apply the spectral model analysis technique to the fluctuation fields derived in the DNS and thereby to derive the energy dissipation rates exactly as it is done for in situ measurements. The derived dissipation fields are to be compared with those ones calculated in DNS. This will yield an assessment of the biases introduced by the spectral model analysis technique.
It is worth noting that the direct measurement of turbulence energy dissipation rates is rather challenging, especially in the natural environment (i.e., atmosphere and ocean). This makes the DNS a unique and valuable tool for the validation of such data analysis techniques and quantification of their precision.
This article does not aim at discussing the merits of theories and underlying assumptions, but to assess the precision and compare uncertainties of the spectral model technique when applying particular spectral models to the analysis of in situ measurements. For detailed discussion and comparison of those assumptions and gained results the reader is referred to for example, Hinze (1975), Reid (1960), and Tatarskii (1971) and to a number of more focused works that address specific topics of analysis techniques or review articles cited in our manuscript.
The article is structured as follows. In the next section, the spectral model analysis technique is briefly described and the main equations are summarized. The DNS data itself and how the analysis is applied to these data are described in Sections 3 and 4, respectively. The results of this analysis are described in Section 5 and critically discussed in Section 6. In Section 7, we summarize the main results. Lübken (1992) developed a practical algorithm to derive turbulence kinetic energy dissipation rate, ɛ, from a measured universal equilibrium range spectrum. The universal equilibrium range of turbulent spectrum includes inertial subrange, where energy transfer occurs from large to small-scales (from low to high wavenumbers) and the inertial forces dominate the motion, and all scales smaller than that (e.g., Hinze, 1975 Due to its advantages, this method has been widely recognized as the so far most accurate approach for dissipation rate measurements in the middle atmosphere (e.g., Hocking, 1999) and has as such been utilized in a large number of sounding rocket experiments (>40 to date, e.g., Strelnikov et al., 2013). Previously in situ turbulence measurements in MLT made use of absolute power values within inertial subrange which resulted in much-enlarged error margins (e.g., Blix et al., 1990;Lübken, 1993;Thrane et al., 1985).

Spectral Model Technique
Here, we shortly summarize the theoretical basis for the spectral model technique. This technique utilizes a single expression spectral model which must simultaneously describe both inertial (inertial-convective) and viscous (viscous-diffusive) subranges for velocity (scalar) fluctuations fields. That is why this method is called the spectral model technique.
Such simple spectral models which can provide suitable estimates for the 1D velocity spectrum E(k) or scalar spectrum E ϑ (k) as a function of energy dissipation rate at a range of wavenumbers evolved from a series of works (e.g., Driscoll & Kennedy, 1981, 1983Heisenberg, 1948;Lübken, 1992;Lübken et al., 1993;Tatarskii, 1971, and references therein). These models for example, of Driscoll and Kennedy (1985), Heisenberg (1948), and Tatarskii (1971) are based on an assumed form for the spectral energy transfer rate (see e.g., Hinze, 1975, for details) and showed a good agreement with universal equilibrium range spectral data measured in the Earth atmosphere (e.g., Lehmacher et al., 2018;Lehmacher & Lübken, 1995;Lübken, 1992Lübken, , 1997Lübken et al., 1993Lübken et al., , 2002Rapp et al., 2004;Strelnikov et al., 2003Strelnikov et al., , 2013. In general case spectra of the scalar field at high wavenumbers beyond the inertial subrange additionally, depend on scalar properties described by the dimensionless numbers Sc or Pr. The Schmidt number Sc = ν/D θ and the molecular Prandtl number Pr = ν/α, describe the ratio of the kinematic viscosity ν to diffusion of atmospheric constituents (D θ diffusion coefficient of tracer θ) or temperature (α is thermal diffusivity), respectively. Batchelor (1959) derived asymptotic expressions for scalar spectra for cases of very high and very low Sc (Pr). These results can be further used to derive an Sc-(Pr-) dependent spectral model (e.g., Driscoll & Kennedy, 1985;Hill, 1978). We do not consider the Sc (Pr) dependencies in this work but only treat cases where the Sc (Pr) value is close to unity, which covers a large enough range of scalar fields and available measurements. Also, in what follows we only deal with a scalar spectrum and the velocity spectrum can be treated in a similar way.
Several works suggested an interpolation formula that describes both inertial-convective and viscous-diffusive subranges (e.g., Driscoll & Kennedy, 1985;Grant et al., 1962a;Heisenberg, 1948;Novikov, 1961;Smith & Reynolds, 1991;Tatarskii, 1971). The spectral model technique aimed at the derivation of the kinetic energy dissipation rate ɛ from a measured spectrum E ϑ . Lübken's idea was to only use the scale (wavenumber) dependence of the spectrum E ϑ (k) and not its absolute level. By fitting a model spectrum to the measured one the scale (wavenumber) of the transition between the inertial-convective and viscous-diffusive subranges, l 0 = 2π/k 0 (inner scale), can be derived quite precisely. Energy dissipation rate is then directly derived from the inner scale l 0 . The advantage of this approach is that normalization of the spectrum does not affect the ɛ-derivation results. In other words, there is no need for precise measurements of absolute values of fluctuations, but only relative ones. By applying some algebra Lübken adapted the original interpolation formulas to the form applicable to measurements. Thus, the adapted Heisenberg (1948) spectrum reads : where k 0 = 2π/l 0 is the wavenumber for inner scale l 0 , a 2 = 1.74 and f a = 2 are constants discussed in Lübken (1992) and in Section 6, and Γ is gamma function.
Similarly, the model of Novikov (1961), also described in the book of Tatarskii (1971) and, after Lübken (1992Lübken ( , 1997 often referred to as "Tatarskii model" is described by the equation (Lübken, 1992): where = ∕ ( 0.033 ⋅ 2 ) 3 is normalized kinetic energy dissipation rate, y = k/k 0 is a dimensionless wavenumber, k 0 is the wavenumber for inner scale l 0 , = The key feature of the adapted models is that they explicitly include l 0 (ɛ) dependence in the form: where η is Kolmogorov scale and the dimensionless constant C is model dependent.
There are different approaches how to deriving the constant C. Thus for example, Obukhov (1949) defined the inner scale l 0 as an intersection of asymptotic extensions of the structure functions (which can be related to the spectrum) in inertial and viscous subranges. Gurvich et al. (1967) suggested deriving this constant empirically based on measured spectra. Lübken utilized relation between second derivative of structure function at zero, H ϑ (0) and 3D spectrum Φ ϑ (e.g., Gurvich et al., 1976;Hinze, 1975;Tatarskii, 1971): The 3D spectrum and its 1D intersection with all the assumptions mentioned above are related via (e.g., Gurvich et al., 1976;Hinze, 1975;Tatarskii, 1971): Combining Equations 4 and 5 Lübken (1992) and Lübken et al. (1993) derived: where superscript H and T denotes Heisenberg or Tatarskii model, respectively. Lübken (1992Lübken ( , 1997 and Lübken et al. (1993) applied the spectral model technique using models of Heisenberg (1948) and Tatarskii (1971), that is, Equations 1 and 2, to relative fluctuations of neutral air density measured in mesosphere. Based on a limited set of data Lübken (1997) and Lübken et al. (1993) showed that application of these models reveals values of the derived energy dissipation rates which are close to each other. Since then mostly the model of Heisenberg (1948) has been applied by scientific community for derivation of turbulence energy dissipation rate, ɛ, based on the Lübken's spectral model technique (e.g., Blix et al., 2003;Chandra et al., 2012;Croskey et al., 2004;Das et al., 2009;Kelley et al., 2003;Lehmacher et al., 2006Lehmacher et al., , 2018Triplett et al., 2018). The main reason for that was the relative simplicity of implementation and numerical stability of the Heisenberg (1948) model. Staszak et al. (2021) and Strelnikov et al. (2017Strelnikov et al. ( , 2019) applied Lübken's technique utilizing both Heisenberg (1948) and Tatarskii (1971) models and showed that the results can reveal considerable discrepancies as far as absolute ɛ-values are concerned, however yielding very similar relative vertical structure and variability.

DNS Data
In this work, we make use of the DNS by  and  where they studied spanwise-and domain-averaged turbulence evolutions and statistics which yields knowledge on the evolution of turbulent patches as a whole, as well as their morphological and dynamical properties. In particular,  studied influences of FS orientation and character on GWs, instability, and turbulence evolutions arising in these flows. Figure 1 shows an example of 2D slices taken from 3D fields obtained by . The dimensions of the shown surfaces are normalized to the vertical size of the simulation domain. For a typical GW-breakdown scenario in the mesosphere, this vertical size will be 3-15 km. The shown 2D-fields are tilted at an angle of ∼5° for consistency and comparability with figures in Fritts et al. ( , 2009aFritts et al. ( , 2009b, , and As in their previous studies,  and  solve the nonlinear Navier-Stokes equations subject to the Boussinesq approximation in a Cartesian domain aligned along with the phase of the primary GW. The equations were non-dimensionalized with respect to the GW vertical wavelength λ z and the buoyancy period, T b = 2π/N. In those DNS, the following parameters were used: a kinematic viscosity ν = 1 m 2 s −1 and a Prandtl number Pr = 1; a sufficiently high value of Reynolds number Re = 2 ∕ = 2 × 10 5 appropriate for a GW in the mesosphere having λ z ∼3-15 km. The mesh size of the DNS at time = 11.5 and 20.0 T b is (NX, NY, NZ)=(4,320, 1,080,2,160) and (2,400,600,1,200), respectively, where X, Y, and Z denotes the streamwise, spanwise, and vertical direction, respectively. The dimension of the computational domain (xL, yL, zL) = (1.98, 0.5, 0.994987)λ z , where λ z is the vertical wavelength of the primary GW. The rate of dissipation of turbulent kinetic energy ɛ is calculated as follows: where ν is kinematic viscosity and s ij is the rates of strain: This results that the dimensions of the ɛ DNS -fields are 2/3 of the dimensions for the velocity (or potential temperature) fields.
An example of distributions of the three parameters obtained in the DNS in vertical-streamwise surfaces, that is, the data to be analyzed in this work is shown in Figure 1. These data were taken at DNS time of t = 11.5 T b when the structures were in their well-developed mature state. In this work, we analyze snapshots of the DNS data taken at different times which includes different stages of turbulence evolution. In the next sections, we will demonstrate the results of fluctuations data analysis using two DNS times t = 11.5 T b and t = 20.0 T b . This will mainly show two largely different stages of fully developed and strongly decayed turbulence from the domain-average point of view. However, the same data also include, as their internal parts, portions of newly created, developed, and decayed structures in smaller regions of the simulation domain. We will address this in detail in Section 6.
In situ measurements (either from rockets, aircraft, or balloons) do only measure a single profile across the 2D-field shown in Figure 1a or 1b. Such a profile is a subject for further analysis using the described in Section 2 spectral model technique.

Analysis Approach
For an incompressible flow (i.e., for motions significantly slower than the speed of sound) under Boussinesq approximation, relative density fluctuations (originally studied by Lübken, 1992) reveal the same structuring as relative fluctuations of potential temperature (e.g., Nappo, 2002): where θ′ and are fluctuations and mean of the potential temperature; ρ′ and are fluctuations and mean values of air density.

of 23
This implies that by analyzing the potential temperature fluctuations derived in this DNS we can directly draw conclusions on the spectral model technique originally introduced by Lübken (1992). By taking a profile from the simulated fluctuations of potential temperature ( Figure 1b) and applying Lübken's spectral model analysis technique (Section 2) one can derive a profile of the turbulence kinetic energy dissipation rate, ɛ. The latter, in turn, can be compared with the profile directly calculated in DNS (Figures 1b and 1c).
As mentioned in Section 3, the original DNS data are dimensionless. To make it representative of MLT dynamics one has to scale the computational domain by a vertical wavelength of GW, λ z . The kinetic energy dissipation rate can be scaled to the real physical units by the factor = 2 ∕ 3 Fritts et al., 2017). For the data demonstrated in this work, we used λ z = 10 km and T b = 5 min, which are quite typical for the MLT region. We also applied different scaling factors using λ z = 5-15 km and T b = 3-8 min and did not note differences in the resulting dissipation fields.
Thus, our analysis approach is as follows. A profile of potential temperature fluctuations taken from the DNS data represents the density fluctuations measured in situ by, for example, a rocket-borne instrument. This profile is to be analyzed by the spectral model technique, yielding a profile of the turbulence kinetic energy dissipation rates, ɛ. We will apply two spectral models, the Heisenberg (1948) and the Tatarskii (1971) model, thereby deriving profiles of ɛ H and ɛ T , respectively. The derived profiles will be compared with the profile of the energy dissipation rate calculated in the DNS, ɛ DNS .
As noted by Fritts et al. (2017), their DNS studies show that a single (or even several sporadic) ɛ-profile(s) cannot adequately characterize the turbulence field in terms of its mean or highest values. Therefore, it makes more sense to obtain some statistics by analyzing vertical-streamwise cross-sections, similar to those shown in Figure 1b, by subsequently deriving ɛ-profiles and, thereby constructing ɛ H -and ɛ T -surfaces for comparison with the ɛ DNS -surface ( Figure 1c). This will also yield a statistical basis for the assessment of biases introduced by the fluctuation data analysis technique.
The exact analysis technique is described in detail by, for example, Strelnikov et al. (2003) or Strelnikov et al. (2013). It is based on theory and models developed by Lübken (1992) and Lübken et al. (1993) and summarized in Section 2, but utilizes wavelet spectral analysis technique instead of the Fourier transform originally used by Lübken. The advantage of the wavelet analysis is that it yields much higher spatial (vertical) resolution, which is theoretically (in an ideal case) the same as for the measured fluctuations profile. In practice, however, it is usually more reasonable to limit the resolution of the analysis (to approximately 30-100 m in the case of rocket measurements in MLT) because of the smoothing properties of the wavelet analysis itself and because of the noisiness of real measurements (see Strelnikov et al., 2003Strelnikov et al., , 2013. In this study, we do not reduce the resolution of the analysis to achieve the most detailed comparison of the turbulence dissipation fields. Also, for the same reason, we interpolate the dissipation fields derived in DNS (ɛ DNS ) to the resolution of the fluctuations data. This makes the ɛ DNS and analysis results ɛ H and ɛ T directly comparable with each other.

Results
In this section, we show the results of the analysis of the potential temperature fluctuations data and compare them with the ɛ DNS -values directly derived in the DNS. First, we show a single profile randomly chosen from the vertical-streamwise cross-section. We note that any profile within the analyzed surfaces shows regions of perfect, good, and strongly biased ɛ-values. Our goal is to find out when the biases occur and quantify how strong these biases are depending on a particular dynamical situation. Next, we compare the entire surfaces of the energy dissipation rates in terms of single values and their statistics. As noted above, the DNS data were scaled to values typical for MLT and the resultant computational domain was between 80 and 90 km altitude. The following discussion will use this altitude range for simplicity.

Profiles
To demonstrate a typical result of the ɛ-derivation we show in Figure 2 profiles of the kinetic energy dissipation rates. The orange profile is directly taken from the DNS data, whereas the blue and green profiles represent the analysis results by using the Heisenberg and Tatarskii spectral models, respectively. It is seen, that in the regions of strong turbulence (ɛ ≳ 1 mW ⋅ kg −1 , above 85 km and around 80 km) both models show values close to the ɛ DNS .
To see more details in the region of a good agreement between ɛ DNS and ɛ H,T we show in Figure 3 a smaller altitude range with the same profiles. It is now seen that the derived energy dissipation rates closely reproduce the general behavior of the ɛ DNS -values directly calculated in DNS. The analysis results, that is, ɛ T and ɛ H , sometimes even coincide with the ɛ DNS -values. The reasons for and implications of the deviations between ɛ DNS and ɛ H,T are discussed in Section 6.
As mentioned in Section 4, when real measurements are analyzed, as a consequence of analysis technique limitations (discussed in Section 6), a smoothing is normally applied. Therefore, to infer the effect of smoothing on the assessment of biases in the estimation of the energy dissipation rates from a single in situ sounding, we show smoothed ɛ-profiles in Figure 4. This plot enlightens several features of the analysis results. First, the general structure of the one-dimensional section of the dissipation field is well reproduced by both ɛ H -and ɛ T -profiles: One can easily recognize major wave-like variations in all three profiles. Herewith the results of the Tatarskii model fit look much closer to the "true", that is, ɛ DNS -values. Second, the high ɛ-values, that is, ɛ ≳ 10 −3 W kg −1 , derived by the spectral model technique based on both models are quite close to the "true" values. Also, both spectral models show results that are close to each other in regions of high energy dissipation rates. In regions of low dissipation, the spectral model analysis results underestimate the amount of energy dissipation. Herewith the Heisenberg model reveals a much stronger bias. At the same time, the Heisenberg results slightly overestimate energy dissipation rates at the peaks of ɛ-profile.

Spectra
In Figure 5, we further demonstrate the performance of the spectral model analysis technique by showing the spectra which yield the energy dissipation rates. The orange line shows a global wavelet spectrum at an altitude of 85.413 km. This altitude is marked by a gray line in Figure 2.
The blue and green lines show the fitted spectra of Heisenberg and Tatarskii models, respectively. The values of energy dissipation rates derived by our analysis are ɛ H = 50 mW kg −1 and ɛ T = 40 mW kg −1 , whereas "true" value calculated in DNS is ɛ DNS = 50 mW kg −1 . We recall, that these ɛ-values are derived from the transition scale l 0 = 2π/k 0 between the inertial-convective and the viscous-diffusive subranges (inner scale) as described in Section 2. The inner scales for the Heisenberg and Tatarskii models are marked by the vertical bold dashed lines in blue and green, respectively. To compare these inner scales with the "true" values inferred from the DNS we show two vertical dashed-dotted lines, which were derived from the ɛ DNS -value. These lines were derived based on the Heisenberg model as and on the Tatarskii model as 0 = 7.06 and are shown in blue and green, respectively. This is an example of a perfect agreement between DNS data and the analysis results. However, already this plot demonstrates how precise (or, in turn, uncertain) are the spectral functions of both models in the dissipation range. One can clearly see that at wavenumbers k ≳ 0.6 cycles m −1 the spectral slopes of both models increasingly deviate from the DNS spectrum. This, however, obviously does not affect the result of the derivation of the energy dissipation rate, ɛ. This is  because the analysis technique only relies on a small part of the spectrum where the transition from the inertial to viscous subrange takes place.
A more detailed analysis of the derived spectra shows that there are different situations of how the DNS spectra are approximated by the model spectra. Figure 6 shows an example when the Tatarskii model with its exponential drop-off in the dissipation range perfectly follows the DNS spectrum. This, however, does not imply coincidence of the energy dissipation rate values ɛ T = 6 × 10 −4 W kg −1 and ɛ DNS = 4 × 10 −4 W kg −1 , even though the difference is not significant. The Heisenberg model, in this case, shows a somewhat opposite situation. The dissipation range slope of k −7 only approximately follows the DNS spectrum and only in the nearby region close to the transition wavenumber (scale). At the same time, the derived energy dissipation rate ɛ H = 7 × 10 −4 W kg −1 is still in an acceptably reasonable agreement with the ɛ DNS -value of 4 × 10 −4 W kg −1 . These spectra correspond to the DNS scaled altitude of 86.367 km. This height is marked in both Figures 2 and 3.
Yet another example of the comparison of DNS with model spectra is shown in Figure 7. In this case, the Tatarskii model demonstrates a somewhat acceptable but far from being a precise approximation of the DNS-spectrum in the dissipation range. At the same time, the derived value of the energy dissipation rate ɛ T = 7 × 10 −2 W kg −1 can be considered as acceptably close to the DNS value of ɛ DNS = 9 × 10 −2 W kg −1 . The Heisenberg model, in turn, follows quite close the DNS spectrum at the beginning of the dissipation range. Whereas the derived value of the energy dissipation rate ɛ H = 2 × 10 −2 W kg −1 is obviously underestimated.
In Figures 8 and 9, we show spectra from the low dissipation part of the profiles shown in Figure 2, that is below 85 km height. In these cases, the approximation of the DNS-spectra by the model-spectra is, like in previous cases, acceptably reasonable. The derived values of the energy dissipation rates are, however, strongly underestimated. These strong biases are discussed in Section 6.

Statistics
After subsequent analysis of every profile of the potential temperature fluctuations in a 2D vertical-streamwise slice of a DNS volume, we reconstruct a surface of the energy dissipation rates.
An example of such a 2D section of the analyzed turbulence field is shown in Figure 10, where panels (a-c) show the "true" ɛ-field, Tatarskii, and Heisenberg model results, respectively.
These figures demonstrate the same features as was inferred from the profile analysis in Section 5.1, but with a stronger statistical basis. Every surface in Figure 10 consists of approximately six thousand profiles or ∼17 million points (single ɛ-values). The main features of the spectral model analysis technique that can be inferred from the comparison of the 2D slices of the "true" and "measured" turbulence fields are as follows.
1. Morphology of the turbulence field, that is, general structure with major features is well reproduced by the analysis regardless of the spectral model used 2. Main regions of strong dissipation are reconstructed quantitatively quite well  3. Analysis technique is not sensitive enough in the regions of weak dissipation, that is, underestimates low ɛ-values 4. Heisenberg model reveals much lower sensitivity to low energy dissipation rate values than the Tatarskii model 5. Heisenberg model tends to overestimate highest ɛ-values 6. Analysis technique is not sensitive enough to resolve the very fine structure of the energy dissipation field Next, in Figure 10d we examine distributions of the energy dissipation rate values from the 2D slices shown in Figures 10a-10c. Histograms in orange, green, and blue show ɛ-distributions for the "true" (DNS), Tatarskii, and Heisenberg analysis results, respectively. Solid lines show Gaussian functions fitted to the respective distributions in the logarithmic domain, that is, represent lognormal distributions of the corresponding energy dissipation rates. Vertical dashed lines mark the median value for each distribution. This figure shows some more details which are not obvious when examining the surface plots shown in Figures 10a-10c. First of all, all three distributions can be described by the Gaussian function acceptably well. Median value inferred from distribution when Tatarskii spectral model applied almost coincides with the median of the "true" DNS distribution. However, both tails of the entire Tatarskii-distribution are slightly expanded relative to the ɛ DNS -distribution. This means, that the highest ɛ T -values are overestimated whereas the lowest values of dissipation rates are underestimated. Distribution of the Heisenberg model results supports the conclusions summarized above and clearly demonstrates that the median ɛ H -value is almost one order of magnitude smaller than the median ɛ DNS .
The statistics have shown so far reflect features of the spectral model analysis technique applied to idealized in situ measurements of well-developed active turbulence. Idealized measurements mean that they are capable of resolving the full range of fluctuations down to the finest scales. By choosing the DNS time t = 11.5 we took for analysis a fully developed active turbulent structure. This implies that the assumptions used in classical turbulence theory are satisfied as much as they can be achieved in these simulations.
In Figure 11, we show another sample of DNS data, taken at a later stage of evolution of the turbulent structure and the analysis results. The DNS time is t = 20.0 meaning that turbulence is already decaying in these data. Even though some classical assumptions of fully developed turbulence most probably do not hold in this case, the key feature for application of the spectral model technique is still present. Namely, at this stage, the decaying turbulence still has prominent inertial and viscous subranges. From an analysis of Figures 11a-11c one can draw the same conclusions as for the case of the developed structure shown above ( Figure 10). However, the histogram plot shown in Figure 11d reveals also some differences if compared with Figure 10d. First, distributions of the results derived using both spectral models are shifted to lower values compared to the developed turbulence case shown in Figure 10d. Second, the distribution width of the Heisenberg model results is significantly narrower than for the developed case and its width is quite close to those of ɛ DNS .

Sensitivity to Instrumental Noise
As noted in the previous section, we applied the spectral model analysis technique to the DNS fluctuations-data assuming there was no instrumental noise. This is what we called idealized measurements. In real measurements, the smallest amplitudes of the measured quantities (e.g., density fluctuations) are usually hidden by instrumental noise. This results in a measured spectrum that only shows a low wavenumber (large scale) part of the viscous subrange. To our knowledge, there are no publications that show spectra measured down to the Kolmogorov scale. This technical imperfection of the in situ measurements motivated us to perform a sensitivity study to assess how experimental limitations affect the analysis results. Figure 12 shows schematics to demonstrate how the instrumental noise affects 1D in situ measurements of turbulence spectra. Bold black curve  shows a spectral function calculated based on Tatarskii model for typical MLT conditions (kinematic viscosity ν = 1 m 2 s −1 ) and turbulence energy dissipation rate ɛ = 1 mW kg −1 .
Black horizontal "tails" to the right of the spectrum show white noise levels from 0.001% to 1.0%. The noise levels are taken as fractions of the maximum amplitude of fluctuations in the spectrum. For example, a noise level of 0.1% means that if measured density fluctuations due to turbulence are at most 2% (e.g., Lehmacher & Lübken, 1995;Lübken, 1992Lübken, , 1997Lübken et al., 1993;Strelnikov et al., 2013), noise flour will hide out all fluctuations smaller than 0.002%. In the spectral domain, these measurements will look like it is shown in Figure 12. The spectra will only be resolved between 10 0 and ∼10 −6 , that is, include six decades of power which is a typical spectral coverage for high-resolution measurements in atmosphere (e.g., Lehmacher & Lübken, 1995;Lübken, 1992Lübken, , 1997Lübken et al., 1993;Söder et al., 2021;Strelnikov et al., 2003Strelnikov et al., , 2013Strelnikov et al., , 2019. The green solid line in Figure 12 shows the part of the spectrum above the noise level of 0.1% which will be fitted by a model. The vertical dashed line shows the inner scale, that is, the visible part of the viscous (viscous-convective) subrange lies between the dashed line and instrumental noise. The shown spectrum is normalized to have its maximum at 10 0 to simplify the estimation of power change between maximum and noise level. It is seen, for example, that an increase of the noise level by factor 10 reduces the visible (resolved by measurements) part of the spectrum by two orders of magnitude. This is because the spectrum is proportional to the square of fluctuations (PSD ∝ Δn 2 ).
In the analyzed DNS data, the large-scale part of turbulence spectra (i.e., to the left of the dashed line in Figure 12) reveals approximately 3-4 decades of power drop and 3.5 decades on average. Note, that it is not necessarily that the inertial (inertial-convective) subrange covers all those large scales. The large-scale (small wavenumber) limit of the inertial subrange does not affect the analysis results and is not discussed in this work. The analysis technique only needs some part of the inertial subrange in the vicinity of the inner sale to be resolved by measurements.
For the sensitivity study, we artificially cut the spectra derived from the DNS fluctuations data below the noise level, as demonstrated in Figure 12 by the green line. Thereby the spectral models were fitted to the "measured" (i.e., DNS) spectra which included inertial-convective subrange and only some part of the viscous-diffusive subrange. By increasing the noise level we shortened the portion of the viscous-diffusive subrange that was used in the fitting process. In this study, we utilized power spectra which covered 8, 6, and 4 orders of magnitude. This approximately corresponds to power drop within the viscous-diffusive subrange of 4.5, 2.5, and 0.5 decades or noise levels of 0.01%, 0.1%, and 1.0%, respectively. Note, that this is not a noise level in terms of a fraction of the dynamical range of an instrument, but a fraction of the largest amplitude of fluctuations produced by turbulence. It is, however, normally possible to relate these quantities in the frame of a defined experiment. Figures 13-15 show the original (i.e., calculated in DNS) and the reconstructed dissipation fields, as well as the related statistical distributions, similar to those shown in Figure 10. Power spectra used for derivation of the ɛ-fields shown in Figures 13-15 were limited to 8, 6, and 4 decades, that is the viscous-diffusive subrange revealed approximately 4.5, 2.5, and 0.5 decades of power change, which is equivalent to noise levels of 0.01%, 0.1%, and 1.0%, respectively. For convenience, hereafter we will refer to these limitations as to spectral coverage, keeping in mind that this describes how much of the viscous subrange is resolved by the measurements.

Developed Turbulence
Results shown in these figures demonstrate the following tendencies: 1. Reduction of spectral coverage (increasing noise level) continuously increases bias in the estimation of ɛ using Tatarskii spectral model 2. Sensitivity of the Tatarskii model to low energy dissipation rates reduces with the reduction of spectral coverage (an increase of noise level)   3. Heisenberg model is less sensitive to the spectral coverage (noise level) within these limits (i.e., demonstrates similar results independent of how much of the viscous subrange is used for the fit) 4. At spectral coverage of 2.5 decades (noise level of 0.1%) both models demonstrate very similar results. This is in accord with the earlier comparisons by Lübken (1992Lübken ( , 1997 and Lübken et al. (1993) 5. For spectral coverage of 0.5 decades (noise level of 1%) median of Heisenberg model results lie closer to the median of the true ɛ-values than the Tatarskii results 6. At the same time, all other features characteristic for an idealized analysis of developed turbulence shown in the previous sections, which do not contradict the 5 listed here items, remain valid

Decaying Turbulence
Next, in Figures 16-18 we show results of the same sensitivity study, but applied to decaying turbulent structures (DNS time t = 20). Interestingly, these results show the same features and lead us to the same conclusions summarized in the previous section. Only a small correction to the last item in that list has to be kept in mind, that the list of the mentioned properties must be extended by the features, characteristic for a decaying structure described at the end of Section 5.3.

Poorly Resolved Viscous Subrange
Further decrease of the spectral range used for the ɛ-derivation gradually increases the negative tendencies of the spectral model analysis technique described above, regardless of a particular model used. The main of them are, that precision of the derived ɛ-values becomes very low and the analysis technique becomes almost insensitive to low energy dissipation rates. Since the large-scale part of the spectra (i.e., down to scale l 0 ) sometimes includes up to four decades of power drop, the spectral coverage of fewer than four decades can completely cut the    viscous-diffusive subrange. In such a case, the fitting process either does not converge or results in a huge fitting error.

Full Spectral Coverage (Low Instrumental Noise)
Statistical basis for analysis of a 2D-slice of the dissipation field discussed in Section 5.3 consists of ∼16.6 and ∼3.4 millions ɛ-values for DNS times 11.5 and 20, respectively. Rigorous derivation of measurement error when applying the spectral model analysis technique to measured spectra of density fluctuations was addressed by Hillert et al. (1994). They showed that the value of ɛ-error (Δɛ H,T ) can be obtained by a proper derivation of the fitting error when applying the least squares technique. However, their error propagation analysis only accounts for the precision of measurements of the tracer and uncertainties in spectral analysis. The fitting errors for our DNS data, are relatively small owing to smooth spectra-a consequence of the idealized measurements. Median fitting errors for both DNS times 11.5 and 20 are 12% and 29% for the Heisenberg and Tatarskii models, respectively. Note, that when spectral models are fitted to turbulent spectra measured in the atmosphere, the fitting errors normally exceed 30% and often reach ∼100%.
Our goal here is to account for the entire scope of possible uncertainties including biases introduced by the spectral models. To assess the distribution of the ɛ-derivation errors we analyzed ratios of the derived to the true values of the energy dissipation rates: ɛ H /ɛ DNS and ɛ T /ɛ DNS . Figure 19 shows these results for active turbulence case (DNS time = 11) in more detail. Bi-dimensional histograms of the two data samples, the derived energy dissipation rates ɛ H,T vs. the ratios ɛ H,T /ɛ DNS are shown in the middle panels of Figures 19a and 19b. The corresponding distributions of ɛ DNS , ɛ H , and ɛ T are shown on the top panels (the same as in Figure 10d).
The bi-dimensional histograms show how the measurement errors (represented by the ratios ɛ H,T /ɛ DNS ) are distributed along with the distributions of the derived ɛ H,T -values (shown in the upper sub-panels). The dashed lines plotted on top of the bi-dimensional histograms show the upper and lower quartiles of the error distributions. That is, the peak of the error distributions for a particular range of ɛ-values lies between the dashed lines. Also, 50% of all the ɛ-values derived within this range lie between the dashed lines and are often referred to as interquartile  Figure 19 show the ratios of ɛ H,T /ɛ DNS = 1, that is where the derived dissipation rates equal the true (ɛ DNS ) value. The right-hand-side panels show histograms of the ratios log 10 (ɛ H,T /ɛ DNS ) for the entire data sets and selection of data around the zero line. The selection was made to mark the region of ɛ-values where analysis yields the most precise results and to see how the distribution of errors in this region looks like.
Thus, it is seen from Figure 19 that the results of analysis using Tatarskii model reveal lowest errors in the range of ɛ-values ∼10 −3 -∼10 −1 W kg −1 . Within this range, 50% of the derived ɛ-values (IQR) have errors lower than half-decade. Whereas for the Heisenberg model the same error is only achieved in the range 10 −2 ≲ ɛ H ≲ 10 −1 W kg −1 . It is also remarkable, that most of the lowest ɛ-values (e.g., all of them beyond the ɛ DNS -distribution) are underestimated whereas the highest values are mostly overestimated.

Limited Spectral Coverage or Dependence on Instrumental Noise
The errors of the energy dissipation rate derivation discussed in the previous section are only relevant for idealized measurements when measured spectra are well resolved down to the smallest scales. To assess the accuracy of the spectral model technique for real measurements we made a series of analyses with artificially reduced resolutions in Section 5.4. From every of those results one can derive the same ratios, that is, ɛ H /ɛ DNS and ɛ T /ɛ DNS , for every derived point (i.e., ɛ-value). In Sections 5.3 and 5.4, we showed results of analysis of eight 2D ɛ-fields for every spectral model, that is, sixteen ɛ-surfaces in total.
Based on the whole statistics of all the derived ɛ-values, we derived median and lower and upper quartiles for the ratios ɛ H,T /ɛ DNS (i.e., the same as Figure 19, but for different resolutions and DNS-times). Thereby we analyzed how many of the derived energy dissipation rates lie close to the true value. It appears that different parts of the ɛ-distribution reveal systematically similar biases for different resolutions and times. To simplify the representation of these results we only consider the median of the ratios ɛ H,T /ɛ DNS . Figure 20 further shows the same curve as the dotted line in the middle panel of Figure 19b (i.e., the median of ɛ-error along the ɛ-distribution) and aims to help in the understanding of the derived statistics. Color-coding (of both line and colorbar) mirrors the ordinate axis. The data were split into ranges of one decade starting from zero and stepping to both positive and negative sides. One special range of 0.0 ± 0.5 decade is additionally marked by white color. Such a curve was made for every instance of our analysis, that is, for different noise levels (i.e., spectral coverage), DNS-times, and spectral models. The error analysis shown in Figure 21 reveals several features. Thus, for example, it is seen, that the right part of the ɛ-distributions (i.e., values to the right side of the median) are more precisely reproduced by both spectral models than low ɛ-values. The best precision is achieved by applying the Tatarskii model to data with low noise levels. Decrease in spectral coverage (i.e., increasing instrumental noise) reduces overall precision of the Tatarskii model results. However, the Tatarskii model shows higher precision than the Heisenberg model for instrumental noise levels above 0.1%, that is, when viscous subrange reveals more than 2.5 decades of power above the noise level. The Heisenberg model, in turn, demonstrates robustness to increasing instrumental noise. If spectral coverage of viscous-convective subrange is decreased to approximately 2.5-2 decades above the noise level, both models demonstrate quite similar results. Although, within a small range of ɛ-values Tatarskii model may reveal slightly lower ɛ-estimates than it will be inferred from the Heisenberg model. At the highest noise level when the viscous-convective subrange reveals ∼0.5 decade of power drop Tatarskii model shows some more underestimates than the Heisenberg model does. At the same time for such noisy data, both models show somewhat   Figure 21 also demonstrates that most of the ɛ DNS -distribution can be approximated by both models with an uncertainty of less than one decade, even when the measured spectra are poorly resolved (i.e., only show half-decade of the viscous-convective subrange).

Discussion
Despite all mentioned imperfections of the spectral model turbulence analysis technique, the analysis results in the regions of moderate to strong dissipation reveal very good agreement with the reference DNS fields. Although the energy dissipation rates are underestimated in the regions of weak turbulence, a general morphology of turbulence in these regions is still reconstructed. That is, if a layer of weak dissipation appears in the DNS data, it also appears in the analysis, though the absolute ɛ-values are smaller.
The results of the assessment of the precision of spectral model turbulence analysis technique shown in previous sections suggest that if measurements allow resolving more than two decades of power for viscous-convective subrange above the noise level, making use of the Tatarskii model yields better overall precision and lower biases at the edges of the actual ɛ-distribution. Also, it better resolves structures in regions where turbulence reveals low dissipation. The higher the spectral resolution of the measurement technique is, the more sensitive is the Tatarskii model to the fine structure of weak dissipation. The best resolved fine structure of turbulence is achieved when the Tatarskii model is applied to the highly resolved spectra, which reveal about six and more decades of power change in the viscous (viscous-diffusive) subrange. However, even in this case, in regions of very weak turbulence analysis underestimates the magnitude of its dissipation considerably. Also, not all fine structure of weak dissipation is reconstructed by the best results of this analysis. The reason for this insensitivity is the limitation of the wavelet spectral analysis technique in the precision of assessment of amplitudes when resolving very fast-changing spectral content. Or, in other words, smoothing properties of the wavelet analysis (e.g., Torrence & Compo, 1998). In this analysis, we applied the Morlet wavelet function of sixths order (e.g., Grossmann & Morlet, 1984) which yields the highest time resolution which is in our case the spatial (altitude) resolution. This represents the main natural limitation of the spectral model turbulence analysis technique. This limitation is due to the width of the wavelet function in the time domain (equivalently spatial domain in our case) leading to that at a given frequency (or wavenumber) the resulting spectral amplitude of a time series under analysis represents an average over a range of the nearest points which is defined by the width of the wavelet function.
Another reason for deviations of the derived energy dissipation rates from the true ɛ-field is the "measurement technique". As noted in Section 2, the measurements are done as a 1D section of the 3D structures. We recall, that the true dissipation field is derived from all three dimensions, that is it accounts for gradients in the fluctuation field perpendicular to the direction of sounding. This can be seen, for example, from Figures 3 and 8 where the good spectral fits yield energy dissipation rates which deviate from the true ɛ DNS -value. This is the reason why energy dissipation rate profiles derived by the spectral model analysis technique shown in Figures 2 and 3 do not exactly reproduce the reference profile ɛ DNS . To address this principal problem in the frame of the spectral model technique it is not only necessary to make 3D soundings, but also to find (either analytically or empirically) a proper 3D spectral function that adequately describes scalar (velocity) spectra in the entire universal range.
The next potential source of uncertainty or biases in the estimation of turbulence energy dissipation rates employing the spectral model technique is the precision of the spectral functions used. The main requirement of these functions is to relate the turbulence kinetic energy dissipation rate with the region of transition from inertial to viscous subranges in wavenumber (or frequency) space as precisely as possible. Whereas it is generally accepted that the inertial (inertial-convective) subrange is precisely described by the k −5/3 power law, there is still no theory that unambiguously defines the spectral function for the viscous (viscous-diffusive) subrange. In fact, there are many suggestions how to describe spectral form in the viscous subrange (e.g., Driscoll & Kennedy, 1981, 1983Gorshkov, 1966;Grant et al., 1962b;Heisenberg, 1948;Hill, 1978;Kovasznay, 1948;Novikov, 1961;Smith & Reynolds, 1991;Tchen, 1973Tchen, , 1975, and many other). However, none of those has received a universally satisfactory confirmation by experiments. All the more uncertain is the approximation of the transition from inertial (inertial-convective) to viscous (viscous-diffusive) subrange in the existing spectral models. This transition is described by interpolation formulas which are not based on physical reasoning but are merely a mathematical convenience.
Since the statistical properties of the viscous subrange are defined by the two physical quantities, ɛ and η, the transition scale l 0 (transition wavenumber k 0 ) must also be defined by these two parameters (e.g., Gurvich et al., 1967;Hinze, 1975;Tatarskii, 1971). For both Heisenberg and Tatarskii models, this dependence is expressed by Equation 3, which states that the transition (inner) scale 0 is proportional to the Kolmogorov scale (see e.g., Gurvich et al., 1967, for a review on this proportionality). The proportionality constant C H,T is of the order 10 as was also noted in early works (e.g., Grant et al., 1962a;Gurvich et al., 1976Gurvich et al., , 1967Hinze, 1975;MacCready, 1962;Pond et al., 1963;Tatarskii, 1971). The range of the suggested values span between 8 and 15 (see e.g., Gurvich et al., 1967). The Kolmogorov scale, in turn, is inversely proportional to 1/4 degree of the energy dissipation rate (η ∝ ɛ (−1/4) ), which makes small changes of l 0 to produce large variations of ɛ. The constants C H,T derived by Lübken (1992) and Lübken et al. (1993), in turn, depend on the constant a 2 or, equivalently, 1 , which are known with a limited precision, as discussed in Sections 2 and 7. The range of a 2 between 2.3 and 3.47, that is, between the lowest possible value (see Section 7) and that one used in our calculations, yields C H between 7.3 and 9.9 and C T between 5.2 and 7.1. This implies an uncertainty of almost four decades for the derivation of ɛ H and one decade for ɛ T . Herewith the lower values of constants C H,T yield lower ɛ. That is, application of lower C H,T -values would introduce an additional negative offset to the derived ɛ H,T -distributions. Making use of the maximum acceptable value for the constant a 2 of 4.02 (see Section 7) will yield ɛ-values which are only twice or half as high as the shown here ɛ-values for the Heisenberg or Tatarskii model, respectively. Taking into account that analysis results yield considerably more underestimates than overestimates, the choice of the constant a 2 = 3.47 looks quite well justified.
After a certain stage of evolution of a turbulence structure, every 2D slice of the DNS volumetric data includes patches of active turbulence and also decaying structures. That is, the turbulence fields derived in this DNS are highly intermittent (e.g., Fritts et al., , 2009b. Detailed comparison of spectra and analysis results for weak and strong, decaying and active turbulences, suggests that the relation between the inner scale l 0 and the energy dissipation rate ɛ given by Equation 3 may be oversimplified. At least, it does not exhibit sufficiently broad universality. Also, the scaling law in wavenumber space for the viscous subrange and, therefore for the transition region, is not precisely described by either of the models in all these considered cases. This fact, however, was already known as a priori (see Section 2) and the spectral models were built upon the assumption of an active and developed turbulence (e.g., Heisenberg, 1948;Hinze, 1975;Tatarskii, 1971). Thus, the better results of this analysis for the developed structures with strong dissipation are somehow expected.

Summary
In this work, we estimated uncertainties and biases in the results of the spectral model turbulence analysis technique applied to in situ measured fluctuations of scalar quantities. Such measurements do only sample fluctuations along one dimension, which forces experimentalists to apply generalized simplifications, for example, to assume isotropy. This, in turn, introduces certain biases in estimated dissipation fields. Uncertainties were determined by the application of the spectral model analysis technique to DNS data, in which ɛ-fields can be rigorously and uniquely determined.
The main results of this study can be summarized as follows.
1. The spectral model technique can reproduce morphology of turbulence field amazingly well and with sufficient details 2. The Tatarskii model reveals high precision of the derived ɛ-values in the range ∼10 −3 -10 −1 W kg −1 if measurements resolve the viscous (viscous-convective) subrange for more than 2 decades of power change, which approximately corresponds to a noise level of 0.1% 3. The Heisenberg model yields a good qualitative picture of the dissipation field, although it is stronger biased than the Tatarskii model A more detailed summary of the uncertainties of the spectral model technique is as follows.
1. This technique robustly detects regions of moderate to strong turbulence with very high precision 2. Kinetic energy dissipation rates derived within such regions reveal uncertainties of less than one order of magnitude 3. At least 50% of those values lie in the 1-sigma interval of their derivation error