Assessment of precision of spectral model turbulence analysis technique using DNS-data

,


Introduction
Turbulence measurements in atmosphere and ocean comprise remote sensing techniques and vast of in situ methods.The most detailed picture of turbulence dissipation or intensity fields can only be acquired by in situ measurements.In situ measurement techniques in turn, utilize different principles depending on altitude (depth) region and consequently its accessibility means.
Different physical quantities inside turbulent flows when measured with a sufficient precision reveal fluctuations around a mean background value.A spectral analysis of these fluctuations shows that they are distributed in a continuous wavenumber space and might obey a mathematical law called spectrum function.In the case of the velocity fluctuations one may find spectrum functions which are the Fourier transform of correlation functions (see e.g., Hinze, 1975).Similar functions can also be applied to describe spectral distribution of other, scalar quantities ϑ, also referred to as tracers.In general case these functions are three dimensional (3D) in space and include time dependency.The most measurement techniques, however, do measure a quasi-instant one-dimensional (1D) cross section of this 3D-spectrum.This technical limitation can normally be circumvented by assuming an isotropy of the spectral distribution of the measured fluctuations allowing for 3D-to-1D transform of the spectrum functions.Also, the measurements must be performed during a short time period so that the time dependence can be neglected.Lübken (1992) introduced a spectral model method for derivation of turbulence energy dissipation rate, ε, based on a theory of spectral distribution of a scalar quantity in turbulence field (see e.g., Tatarskii, 1971;Hinze, 1975;A. M. Obukhov, 1988).This technique was successfully applied to fluctuations of neutral air density measured in mesosphere by sounding rockets (e.g., Lübken, 1992Lübken, , 1997;;Lübken et al., 1993Lübken et al., , 2002;;Strelnikov et al., 2003Strelnikov et al., , 2013Strelnikov et al., , 2017Strelnikov et al., , 2019;;Szewczyk et al., 2013) and to velocity fluctuations measured by stratospheric balloons Theuerkauf et al. (2011); Haack et al. (2014); Schneider et al. (2015Schneider et al. ( , 2017)); Söder et al. (2019Söder et al. ( , 2020)).However, the underlying theory and therefore the ε-derivation technique are based on assumptions which might introduce some uncertainties, which are not quantified yet.
In 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 et al. (1997); Arendt (1998); Andreassen et al. (1998); 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., 2009b,a).Highly resolved velocity field produced by such DNS allows for precise and detailed derivation of kinetic energy dissipation rate with 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 details in e.g., Fritts et al. (2009b,a).
More recent DNS by Fritts et al. (2013) and Fritts & Wang (2013) studied multiscale dynamics (MSD) 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 fine structure of the velocity and dissipation fields and its 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 validation of such data analysis techniques and for quantification of their precision.This paper 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 analysis of in situ measurements.For detailed discussion and comparison of those assumptions and gained results the reader is referred to e.g., Reid (1960); Tatarskii (1971); Hinze (1975) and to number of more focused works that address specific topics of analysis techniques or review articles cited in our manuscript.
The paper is structured as follows.In the next section the spectral model analysis technique is briefly described and main equations are summarized.The DNS data itself and how the analysis is applied to these data are described in Sec. 3 and 4, respectively.The results of this analysis are described in Sec. 5 and critically discussed in Sec. 6.
In Sec. 7 we summarize the main results.
2 Spectral model technique 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).
Here we shortly summarize 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 spectral model technique.
In general case spectra of scalar field at high wavenumbers beyond the inertial subrange additionally depend on scalar properties described by the dimensionless numbers Sc or Pr.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 a Sc-(Pr-) dependent spectral model (e.g., Hill, 1978;Driscoll & Kennedy, 1985).We do not consider the Sc (Pr) dependencies in this work but only treat cases where Sc (Pr) value is close to unity, which covers 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 which describes both inertial-convective and viscous-diffusive subranges (e.g., Heisenberg, 1948;Novikov, 1961;Grant et al., 1962a;Tatarskii, 1971;Driscoll & Kennedy, 1985;Smith & Reynolds, 1991).The spectral model technique aimed at 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 viscousdiffusive 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 (Lübken et al., 1993): 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 Sec.6, and Γ is gamma function.
The key feature of the adapted models is that they explicitly include l 0 (ε) dependence in the form: (3) -4-manuscript submitted to JGR: Atmospheres where η is Kolmogorov scale and the dimensionless constant C is model dependent.
There are different approaches how to derive the constant C. Thus e.g., A. Obukhov (1949) defined the inner scale l 0 as intersection of asymptotic extensions of the structure functions (which can be related to the spectrum) in inertial and viscous subranges.
Lübken (1992), Lübken et al. (1993), andLübken (1997) applied the spectral model technique using models of Heisenberg (1948) and Tatarskii (1971), i.e.Eq. 1 and 2, to relative fluctuations of neutral air density measured in mesosphere.Based on a limited set of data Lübken et al. (1993) and Lübken (1997) showed that application of these models reveals values of the derived energy dissipation rates which are close to each other.

DNS data
In this work we make use of the DNS by Fritts et al. (2013) and Fritts & Wang (2013) where they studied spanwise-and domain-averaged turbulence evolutions and statistics which yields knowledge on evolution of turbulent patches as whole, as well as their morphological and dynamical properties.In particular, Fritts et al. (2013) studied influences of FS orientation and character on GWs, instability, and turbulence evolutions arising in these flows.
Fig. 1 shows an example of 2D slices taken from 3D fields obtained by Fritts et al. (2013).The dimensions of the shown surfaces are normalized to the vertical size of the simulation domain.For a typical GW-breakdown scenario in mesosphere this vertical size will be 3 to 15 km.The shown 2D-fields are tilted at an angle of ∼ 5 • for consistency and comparability with figures in Fritts et al. (2009bFritts et al. ( ,a, 2013)), Fritts &Wang (2013), andFritts et al. (2017).
As in their previous studies Fritts et al. (2013) and Fritts & Wang (2013) solve the nonlinear Navier-Stokes equations subject to the Boussinesq approximation in a Cartesian domain aligned along 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 P r = 1; a sufficiently high value of Reynolds number Re = λ 2 z /νT b = 2 × 10 5 appropriate for a GW in the mesosphere having λ z ∼3 to 15 km.
The dissipation data calculated in these DNS are directly derived from the gradients of the velocity fluctuations (see e.g., Landau & Lifshitz, 1987) which results that the dimensions of the ε DN S -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 verticalstreamwise surfaces, i.e. the data to be analyzed in this work is shown in Fig. 1.These data were taken at DNS time of t = 11.5 T b when the structures were in its 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.0T 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 Sec. 6.
In situ measurements (either from rockets, aircraft, or balloons) do only measure a single profile across the 2D-field shown in Fig. 1a or Fig. 1b.Such a profile is a subject for further analysis using the described in Sec. 2 spectral model technique.

Analysis approach
For an incompressible flow (i.e. for motions significantly slower than 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.
This implies that by analyzing the potential temperature fluctuations derived in these 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 (Fig. 1b) and applying Lübken's spectral model analysis technique (Sec.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 (Fig. 1c).
As mentioned in Sec. 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 S ε = λ 2 z /T 3 b (Fritts & Wang, 2013;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 MLT region.
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 profile of the energy dissipation rate calculated in the DNS, ε DN S .
As noted by Fritts et al. (2017), their DNS studies show that a single (or even several sporadic) ε-profile(s) cannot adequately characterize turbulence field in terms of their mean ot highest values.Therefore, it makes more sense to obtain some statistics by analyzing vertical-streamwise cross sections, similar to those shown in Fig. 1b, by subsequently deriving ε-profiles and, thereby constructing ε H -and ε T -surfaces for compar-ison with the ε DN S -surface (Fig. 1c).This will also yield a statistical basis for assessment of biases introduced by the fluctuation data analysis technique.
The exact analysis technique is described in detail by e.g., 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 Sec. 2, but utilizes wavelet spectral analysis technique instead of the Fourier transform originally used by Lübken.Advantage of the wavelet analysis is that it yields much higher spatial (vertical) resolution, theoretically (in 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 to 100 m in case of rocket measurements in MLT) because of smoothing properties of the wavelet analysis itself and because of noisiness of real measurements (see Strelnikov et al., 2003Strelnikov et al., , 2013, for details), for details).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 (ε DN S ) to the resolution of fluctuations data.This makes the ε DN S and analysis results ε H and ε T to be directly comparable with each other.

Results
In this section we show the results of analysis of the potential temperature fluctuations data and compare them with the ε DN S -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 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 Fig. 2 profiles of the kinetic energy dissipation rates.The blue profile is directly taken from the DNS data whereas orange 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 ε DN S .In the region where DNS reveals low ε-values analysis results show different deviations.Mean ratios of the derived-to-DNS ε-values are ε T /ε DN S =1.07 and ε H /ε DN S =1.14 for Tatarskii and Heisenberg models, respectively.
To see more details in the region of a good agreement between ε DN S and ε H,T we show in Fig. 3 a smaller altitude range with the same profiles.It is now seen that the derived energy dissipation rates closely reproduce general behavior of the ε DN S -values directly calculated in DNS.The analysis results, i.e. ε T and ε H , sometimes even coincide with the ε DN S -values.The reasons for and implications of the deviations between ε DN S and ε H,T are discussed in Sec. 6.
As mentioned in Sec. 4, when real measurements are analyzed, as a consequence of analysis technique limitations (discussed in Sec.6), a smoothing is normally applied.
Therefore, to infer the effect of smoothing on the assessment of biases in estimation of the energy dissipation rates from a single in situ sounding, we show smoothed ε-profiles in Fig. 4.This plot enlightens several features of the analysis results.First, general structure of the 1D section of 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 ) models in orange and green, respectively.Vertical dashed-dotted lines show the inner scales derived from the DNS data (ε DN S -value) based on the Heisenberg model (l H 0 = 9.9(ν 3 /ε DN S ) 1/4 ) and Tatarskii (l T 0 = 7.06(ν 3 /ε DN S ) 1/4 ) model in orange and green, respectively.
The orange 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 ε DN S =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      Sec. 2. The inner scales for the Heisenberg and Tatarskii models are marked by the vertical bold dashed lines in orange 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 ε DN S -value.These lines were derived based on the Heisenberg model as l H 0 = 9.9(ν 3 /ε DN S ) 1/4 and on the Tatarskii model as l T 0 = 7.06(ν 3 /ε DN S ) 1/4 , and are shown in orange and green, respectively.This is an example of 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 the spectral slopes of the both models increasingly deviate from the DNS spectrum.This, however, obviously does not affect the result of 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.Fig. 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 ε DN S =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 ε DN S -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 Fig. 2 and 3.
Yet another example of the comparison of DNS with model spectra is shown in Fig. 7.
In this case the Tatarskii model demonstrates a somewhat acceptable but far from being 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 ε DN S =9•10 −2 W kg −1 .The Heisenberg model, in turn, follows quite close the DNS spectrum in 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 Fig. 8 and 9 we show spectra from the low dissipation part of the profiles shown in Fig. 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 Sec. 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 Fig. 10, where panels a, b, and 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 Sec.5.1, but with stronger statistical basis.Every surface in Fig. 10 consist of approximately six thousands profiles or ∼17 millions points (single ε-values).The main fea-   tures 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.
• Morphology of the turbulence field, i.e. general structure with major features is well reproduced by the analysis regardless of spectral model used.
• Main regions of strong dissipation are reconstructed quantitatively quite well.
• Analysis technique is not sensitive enough in the regions of weak dissipation, i.e. underestimates low ε-values.
• Heisenberg model reveals much lower sensitivity to low energy dissipation rate values than the Tatarskii model.The statistics shown so far reflects features of the spectral model analysis technique applied to idealized in situ measurements of well developed active turbulence.Idealized measurements means that they are capable of resolving full range of fluctuations down to 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 it can be achieved in these simulations.
In Fig. 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 the spectral model technique is still present.Namely, at this stage the decaying turbulence still has a prominent inertial and the viscous subranges.From analysis of Fig. 11a-c one can draw the same conclusions as for the case of the developed structure shown above (Fig. 10).However, the histogram plot shown in Fig. 11d reveals also some differences if compared with Fig. 10d.First, distributions of the results derived using both spectral models are shifted to lower values compared to the developed turbulence case shown in Fig. 10d.Second, distribution width of the Heisenberg model results is significantly narrower than for the developed case and its width is quite close to those of ε DN S .

Sensitivity to instrumental noise
As noted in the previous section, we applied the spectral model analysis technique to the DNS fluctuations-data assuming there were 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 which only shows a low wavenumber (large scale) part of the viscous subrange.To our knowledge there are no publications which show spectra measured down to Kolmogorov scale.This technical imperfection of the in situ measurements motivated us to perform a sensitivity study to asses how experimental limitations affect the analysis results.Black horizontal "tails" to the right of the spectrum show white noise levels from 0.001-1.0%.The noise levels are taken as fractions of maximum amplitude of fluctuations in spectrum.For example, noise level of 0.1% means that if measured density fluctuations due to turbulence are at most 2% (e.g., Lübken, 1992Lübken, , 1997;;Lübken et al., 1993;G. Lehmacher & Lübken, 1995;Strelnikov et al., 2013), noise flour will hide out all fluctuations smaller than 0.002%.In spectral domain these measurements will look like it is shown in Fig. 12.The spectra will only be resolved between 10 0 and ∼10 −6 , i.e. include six decades of power which is a typical spectral coverage for high resolution measurements in atmosphere (e.g., Lübken, 1992Lübken, , 1997;;Lübken et al., 1993;G. Lehmacher & Lübken, 1995;Strelnikov et al., 2003Strelnikov et al., , 2013Strelnikov et al., , 2019;;Söder et al., 2021).Green solid line in Fig. 12 shows the part of the spectrum above the noise level of 0.1 % which will be fitted by a model.Vertical dashed line shows the inner scale, i.e., 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 estimation of power change between maximum and noise level.It is seen, e.g., that an increase of the noise level by factor 10 reduces visible (resolved by measurements) part of spectrum by two orders of magnitude.This is because the spectrum is proportional to the square of fluctuations (P SD ∝ ∆n 2 ).
In the analyzed DNS data the large-scale part of turbulence spectra (i.e. to the left of the dashed line in Fig. 12) reveal approximately 3 to 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 Fig. 12 by the green line.Thereby the spectral models were fitted to the "measured" (i.e.DNS) spectra which included inertialconvective 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 to noise levels of 0.01, 0.1, and 1.0 %, respectively.Note, that this is not a noise level in terms of fraction of dynamical range of instrument, but a fraction of largest amplitude of fluctuations produced by turbulence.It is, however, normally possible to relate these quantities in the frame of a defined experiment.

Developed turbulence
Fig. 13, 14, and 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 Fig. 10.Power spectra used for derivation of the ε-fields shown in Fig. 13, 14, and 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 this limitations as to spectral coverage, keeping in mind that this describes how much of the viscous subrange is resolved by the measurements.
Results shown in these figures demonstrate the following tendencies: • Reduction of spectral coverage (increasing noise level) continuously increases bias in estimation of ε using Tatarskii spectral model.
• Sensitivity of the Tatarskii model to low energy dissipation rates reduces with the reduction of spectral coverage (increase of noise level).
• 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).
• For spectral coverage of 0.5 decade (noise level of 1 %) median of Heisenberg model results lies closer to the median of the true ε-values than the Tatarskii results.
• At the same time, all other features characteristic for an idealized analysis of a developed turbulence shown in the previous sections, which do not contradict 5 listed here items, remain valid.

Decaying turbulence
Next, in Fig. 16, 17, and 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 in the end of Sec.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 less 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 Sec.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 tech-   nique 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 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 model, 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 introduces by the spectral models.To assess distribution of the ε-derivation errors we analyzed ratios of the derived to the true values of the energy dissipation rates: ε H /ε DN S and ε T /ε DN S .Fig. 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 versus the ratios ε H,T /ε DN S are shown in the middle panels of Figs.19a     and b.The corresponding distributions of ε DN S , ε H , and ε T are shown on the top panels (the same as in Fig. 10d).
The bi-dimensional histograms show how the measurement errors (represented by the ratios ε H,T /ε DN S ) are distributed along 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, 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 range (IQR).The dotted lines plotted on the bi-dimensional histograms show medians of the measurement errors (i.e., of the ratios ε H,T /ε DN S ).The   horizontal solid zero lines in the middle panels of Fig. 19 show the ratios of ε H,T /ε DN S = 1, that is where the derived dissipation rates equal the true (ε DN S ) value.The right-handside panels show histograms of the ratios log 10 (ε H,T /ε DN S ) for the entire data sets and for selection of data around the zero line.The selection was made to mark region of εvalues where analysis yields most precise results and to see how the distribution of errors in this region looks like.Thus, it is seen from Fig. 19 that the results of analysis using Tatarskii model reveal lowest errors in the range of ε-values ∼10 −3 to ∼10 −1 W kg −1 .Within this range 50% of the derived ε-values (IQR) have error 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 ε DN Sdistribution) 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 an idealized measurements when measured spectra are well resolved down to 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 Sec.5.4.From every of those results one can derive the same ratios, i.e. ε H /ε DN S and ε T /ε DN S , for every derived point (i.e., ε-value).In Sec.5.3 and 5.4 we showed results of analysis of eight 2D ε-fields for every spectral model, i.e. 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 /ε DN S (i.e., the same as Fig. 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 is appeared, that different parts of the εdistribution reveal systematically similar biases for different resolutions and times.
To simplify representation of these results we only consider the median of the ratios ε H,T /ε DN S .Fig. 20 further shows the same curve as the dotted line in the middle panel of Fig. 19b (i.e., median of ε-error along the ε-distribution) and aims to help in understanding of the derived statistics.Color-coding (of both line and colorbar) mirrors the ordinate axis.The data were split in 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, i.e. for different noise levels (i.e., spectral coverage), DNS-times, and spectral models.

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 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 precision of spectral model turbulence analysis technique shown in previous sections suggest that if measurements allow to resolve more than two decades of power for viscous-convective subrange above 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 measurement technique is, the more sensitive is the Tatarskii model to 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 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 limitation of the wavelet spectral analysis technique in 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 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 range of the nearest points which is defined by the width of the wavelet function.
Another reason of deviations of the derived energy dissipation rates from the true ε-field is the "measurement technique".As noted in Sec. 2, the measurements are done as a one dimensional 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 fluctuation field perpendicular to the direction of sounding.This can be seen, e.g. from Fig. 3 and 8 where the good spectral fits yield energy dissipation rates which deviate from the true ε DN S -value.This is the reason why energy dissipation rate profiles derived by the spectral model analysis technique shown in Fig. 2 and 3 do not exactly reproduce the reference profile ε DN S .To address this principal problem in 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 which adequately describes scalar (velocity) spectra in the entire universal range.
The next potential source of uncertainty or biases in estimation of turbulence energy dissipation rates by means of the spectral model technique is the precision of the spectral functions used.The main requirement to 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 which 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., Heisenberg, 1948;Kovasznay, 1948;Novikov, 1961;Grant et al., 1962b;Gorshkov, 1966;Tchen, 1973Tchen, , 1975;;Hill, 1978;Driscoll & Kennedy, 1981, 1983, 1985;Smith & Reynolds, 1991, 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 a physical reasoning but they 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., A. Gurvich et al., 1967;Tatarskii, 1971;Hinze, 1975).For both Heisenberg and Tatarskii models this dependence is expressed by Eq. 3, which states that the transition (inner) scale l H,T 0 is proportional to the Kolmogorov scale (see e.g., A. Gurvich et al., 1967, for a review on this proportionality).The proportionality constant C H,T is of the order ten as was also noted in early works (e.g., MacCready Jr., 1962;Grant et al., 1962a;Pond et al., 1963;A. Gurvich et al., 1967;Tatarskii, 1971;Hinze, 1975;A. S. Gurvich et al., 1976).The range of the suggested values span between 8 and 15 (see e.g., A. 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, C 1 ϑ , which are known with a limited precision, as discussed in Sec. 2 and 7.The range of a 2 between 2.3 and 3.47, i.e. between the lowest possible value (see Sec. 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 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 Sec. 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 these DNS are highly intermittent (e.g., Fritts et al., 2009bFritts et al., , 2013)).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 Eq. 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 obviously not precisely described by either of models in all these considered cases.This fact, however, was already known a priori (see Sec. 2) and moreover, the spectral models were build upon assumption of active developed turbulence (e.g., Heisenberg, 1948;Tatarskii, 1971;Hinze, 1975).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 results of 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, e.g. to assume isotropy.This, in turn introduces certain biases in estimated dissipation fields.Uncertainties were determined by 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.
• The spectral model technique can reproduce morphology of turbulence field amazingly well and with sufficient details.
• The Tatarskii model reveals high precision of the derived ε-values in the range ∼10 −3 to 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 noise level of 0.1 %.
• The Heisenberg model yields a good qualitative picture of the dissipation field, although it is stronger biased than the Tatarskii model.
-25-manuscript submitted to JGR: Atmospheres Some more detailed summary of the uncertainties of the spectral model technique are as follows.
• This technique robustly detects regions of moderate to strong turbulence with very high precision.
• Kinetic energy dissipation rates derived within such regions reveal uncertainties of less than one order of magnitude.
• At least 50 % of those values lie in 1-sigma interval of their derivation error.
• The minimum spectral coverage needed to reliably apply spectral model technique only requires half decade of power drop within viscous (viscous-convective) subrange (which corresponds to noise level of 1 % of fluctuations' amplitude).
• If the viscous (viscous-convective) subrange is resolved to reveal at least two decades of power drop, the models of Heisenberg (1948) and Tatarskii (1971) demonstrate similar results and relatively high precision.
• If the viscous (viscous-convective) subrange is resolved within more than two decades of power change, the model of Tatarskii (1971) shows more accurate results and reveals relatively high sensitivity to low ε-values.
• The spectral model of Heisenberg (1948), on the other hand, is almost insensitive to the quality of measured spectra (i.e., reveals near the same accuracy regardless of how much of the viscous subrange is resolved by measurements).
Specifically for MLT, that is taking into account the applied scaling of the dimensionless DNS data we can additionally highlight several features.
• Low values of energy dissipation rates, i.e. ε 1 mW•kg −1 are mostly underestimated, meaning that the true ε-value can exceed the measured ones.
• If the derived energy dissipation rates lie in the range between ∼ 2•10 −5 W kg −1 and ∼ 1 W kg −1 , their value does not deviate from the true ε-value by more than one decade with probability of 50 %.
With all the uncertainties critically discussed above, the spectral model analysis technique of in situ measurements reproduces the ε-reference fields not only amazingly well, but also in much more details compared to other techniques available for atmospheric or oceanographic turbulence soundings.
Appendix A: Uncertainties of constants used in spectral functions Eq. 6 and 7 show that the constants f a and a 2 are explicitly used to derive the constant C which connects the inner scale l 0 and the energy dissipation rate ε.The constant f a was introduced by Lübken (1992) to make it possible to apply the same formulae for both energy (i.e.velocity) and scalar spectra.For energy and scalar spectra f a takes values of 1 and 2, respectively.The constant a 2 is somewhat worse defined.It appears from derivation of the Obukhov-Corrsin law for the inertial subrange when comparing different derivation approaches.Constant a 2 , in particular can be related to the Obukhov-Corrsin constant C 1 ϑ as (see e.g., Tatarskii et al., 1992): Since in the inertial-convective subrange the 3D-spectrum has the same form as the 1Dspectrum, it must be distinguished between the Obukhov-Corrsin constants for these cases, with C 1 ϑ replaced by a different constant C ϑ for 3D spectrum.Isotropy implies that they are related as (e.g., Hill, 1978;Sreenivasan, 1996): Potential temperature fluctuations.(c) Kinetic energy dissipation rate.

Figure 2 :Figure 3 :
Figure2: Example of vertical profiles of the derived energy dissipation rates.Blue profile shows the DNS data, whereas the orange and green profiles show the analysis results using the Heisenberg and Tatarskii spectral models, respectively.gray bold horizontal lines mark altitudes where power spectra are taken from for demonstration in section 5.2.

Figure 5 :
Figure 5: Example of power spectra which yield the ε-profiles shown in Fig. 2 taken at an altitude of 85.413 km.Blue, orange, and green lines show the DNS, Heisenberg, and Tatarskii data.Bold vertical dashed lines show the inner scales (l 0 = 2π/k 0 ) derived from the fit of the Heisenberg (l H 0 4e-04 W/kg Heisenberg fit; = 7e-04 W/kg ±16% Tatarskii fit; = 6e-04 W/kg

Figure 6 :
Figure 6: Same as Fig. 5, but for an altitude of 86.367 km.

Figure 7 :
Figure 7: Same as Fig. 5, but for an altitude of 86.470 km.

Figure 8 :
Figure 8: Same as Fig. 5, but for an altitude of 83.109 km.

Figure 9 :
Figure 9: Same as Fig. 5, but for an altitude of 82.073 km.

Figure 10 :
Figure 10: 2D fields of the kinetic energy dissipation rates.Panel (a) shows the "true" DNS data used as reference (the same as Fig. 1c).DNS-time=11.5,i.e. for well developed turbulence.Panels (b) and (c) show analysis results using Tatarskii and Heisenberg models, respectively.Lighter colors correspond to higher ε-values.Panel (d) shows the same data as in panels (a), (b), and (c), but as histograms of ε-distributions and fitted PDFs.Vertical dashed lines show medians of corresponding data sets (here the DNS and Tatarskii results almost coincide and are hardly distinguishable).

•
Heisenberg model tends to overestimate highest ε-values.• Analysis technique is not sensitive enough to resolve very fine structure of the energy dissipation field.Next, in Fig. 10d we examine distributions of the energy dissipation rate values from the 2D slices shown in Fig. 10a, b, and c.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 logarithmic domain, i.e. represent lognormal distributions of the corresponding energy dissipation rates.Vertical dashed lines mark median value for each distribution.This figure shows some more details which are not obvious when examining the surface plots shown

Fig. 12 Figure 12 :
Fig. 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 .

Figure 13 :
Figure 13: Same as Fig. 10, but for noised spectra with 4.5 decades of visible power within the viscous subrange (noise level of 0.01 %).

Figure 14 :
Figure 14: Same as Fig. 10, but for noised spectra with 2.5 decades of visible power within the viscous subrange (noise level of 0.1 %).

Figure 15 :
Figure 15: Same as Fig. 10, but for noised spectra with 0.5 decades of visible power within the viscous subrange (noise level of 1 %).

Figure 16 :
Figure 16: Same as Fig. 11, but for noised spectra with 4.5 decades of visible power within the viscous subrange (noise level of 0.01 %).

Fig. 21
Fig. 21 shows compilation of these analysis results, where eight curves like that one in Fig. 20 are shown for every spectral model.Abscissa in Fig. 21 shows energy dissipation rates in logarithmic scale, log 10 (ε), and the orange curves schematically show the PDFs of the ε DN S -distributions.Upper and lower panels in Fig. 21a and 21b show results for active and decaying turbulence (DNS times 11.5 and 20), respectively.Left and right panels (Fig. 21a and 21b) show results for the Heisenberg and Tatarskii spectral model, respectively.Reddish colors show regions where the ratios ε H,T /ε DN S are greater than unity, that is the derived values ε H and ε T are overestimated.Blueish colors show regions where the derived energy dissipation rates are underestimated.Gray color marks region outside the derived range of values.The error analysis shown in Fig. 21 reveals several features.Thus, e.g., it is clearly 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.De-

Figure 17 :
Figure 17: Same as Fig. 11, but for noised spectra with 2.5 decades of visible power within the viscous subrange (noise level of 0.1 %).

Figure 18 :
Figure 18: Same as Fig. 11, but for noised spectra with 0.5 decades of visible power within the viscous subrange (noise level of 1 %).

Figure 19 :
Figure19: Distributions of the derived energy dissipation rates and of the errors (represented by the ratios ε H,T /ε DN S ) in logarithmic scale.The left (a) and right (b) subfigures are for the Heisenberg and Tatarskii model, respectively.Middle panels: Bi-dimensional histograms of the derived energy dissipation rates (ε H,T ) and of the measurement errors (ε H,T /ε DN S ).Solid horizontal lines show the ratio ε H,T /ε DN S = 1.The dashed lines show upper and lower quartiles of the respective ratio ε H,T /ε DN S , i.e. the area between dashed lines shows the interquartile range for the ε-derivation error.The dotted line shows the median error.Upper panels: Distributions of ε H , ε T , and ε DN S in blue, green, and orange, respectively.Right-hand-side panels: Distributions of the ratios ε H,T /ε DN S .Red color in the mid-and right-panels shows selection of data with errors within one decade around the zero-line: -0.5<log 10 (ε H,T /ε DN S ) <0.5.Red histograms show distributions of errors for the selection of data.The black dotted lines show the lower and upper quartiles for red histograms.

Figure 20 :
Figure 20: Example of the derived ε-error represented by the ratio in logarithmic scale: log 10 (ε T /ε DN S ).Color-coding of both line and colorbar mirrors the ordinate axis.White color shows range of ε H,T -values where their error falls within half-decade interval: −0.5 < log 10 (ε H,T /ε DN S ) < 0.5.Red and blue colors show regions where derived ε H,T are over-and underestimated, respectively.

Figure 21 :
Figure 21: Ratios of the derived to true energy dissipation rates in logarithmic scale: log 10 (ε H,T /ε DN S ) shown by colors as a function of ε H,T -value (abscissa) and spectral resolution (ordinate).White color shows range of ε H,T -values where their error falls within half-decade interval: −0.5 < log 10 (ε H,T /ε DN S ) < 0.5.Red and blue colors show regions where derived ε H,T are over-and underestimated, respectively.Orange Gaussians schematically show PDF(ε DN S ).