Estimation of Turbulence Parameters in the Lower Troposphere from ShUREX (2016–2017) UAV Data

: Turbulence parameters in the lower troposphere (up to ~4.5 km) are estimated from measurements of high-resolution and fast-response cold-wire temperature and Pitot tube velocity from sensors onboard DataHawk Unmanned Aerial Vehicles (UAVs) operated at the Shigaraki Middle and Upper atmosphere (MU) Observatory during two ShUREX (Shigaraki UAV Radar Experiment) campaigns in 2016 and 2017. The practical processing methods used for estimating turbulence kinetic energy dissipation rate ε and temperature structure function parameter C 2 T from one-dimensional wind and temperature frequency spectra are ﬁrst described in detail. Both are based on the identiﬁcation of inertial ( − 5 / 3) subranges in respective spectra. Using a formulation relating ε and C 2 T valid for Kolmogorov turbulence in steady state, the ﬂux Richardson number R f and the mixing e ﬃ ciency χ m are then estimated. The statistical analysis conﬁrms the variability of R f and χ m around ∼ 0.13 − 0.14 and ∼ 0.16 − 0.17, respectively, values close to the canonical values found from some earlier experimental and theoretical studies of both the atmosphere and the oceans. relevance of the interpretation of the inertial subranges in terms of Kolmogorov turbulence is conﬁrmed by assessing the consistency of additional parameters, the Ozmidov length scale L O , the buoyancy Reynolds number Re b , and the gradient Richardson number Ri . Finally, a case study is presented showing altitude di ﬀ erences between the peaks of N 2 , C 2 T and ε , suggesting turbulent stirring at the margin of a stable temperature gradient sheet. The possible contribution of this sheet and layer structure on clear air radar backscattering mechanisms is examined.


Introduction
Turbulence is ubiquitous in the free, stably stratified, atmosphere [1][2][3][4][5] and can be the result of various dynamical instabilities [6]. More generally, a large variety of small-scale structures and dynamics can occur due to complex interactions between background winds and shears, internal gravity waves, convective motions, clouds, etc. These interactions may lead to the observed "sheet and layer" structure [6,7]. This refers to the alternation of steep and thin temperature gradient sheets and deeper layers of weaker stability. Despite these complex dynamics, small-scale turbulence in stratified conditions is often parameterized by the Kolmogorov model valid for locally homogeneous, stationary, isotropic and inertial turbulence. The intensities of the mechanical and temperature turbulence are then described by the turbulence kinetic energy (TKE) dissipation rate ε and temperature structure function parameter C 2 T , respectively. Testing the relevance of these parameters and estimating these parameters accurately are important for correct parameterization of turbulence dissipation and diffusion in theoretical and numerical models of atmospheric dynamics and composition. As different instruments, where g is the acceleration of gravity, and N 2 = g/θ(dθ/dz) is the squared Brunt-Vaïsälä angular frequency expressed in terms of (dry) potential temperature θ. The parameter γ is: where β θ ≈ 3.0 is a universal constant and R f is the flux Richardson number. R f is the ratio of the buoyancy flux to TKE production by shear and is a key turbulence parameter. In literature on oceanic turbulence, a constant value of 0.17 is assumed for R f , even though laboratory experiments, direct numerical simulations and field observations show that R f varies with gradient Richardson number Ri [22][23][24][25]. However, R f is usually observed to vary around a canonical value of~0.17 expected from theory [23]. In literature on atmospheric turbulence, R f is often taken equal to 0.25 [2,3,24]. The direct estimate of ε and use of Equation (1) offer the possibility to evaluate R f in the lower troposphere.
The term R f / 1 − R f is equal to the ratio of the buoyancy flux to the kinetic energy dissipation rate and is commonly called 'mixing (efficiency) coefficient' χ m [25,26]. It is also a key parameter in modelling turbulent flows and appears in the expression of the scalar eddy diffusivity K θ [27]: Based on a part of the dataset collected during ShUREX 2016 and discussed here, a mean value of χ m ∼ 0.16 was found [21], in very close agreement with typical values for oceanic turbulence [25,26] and for turbulence at stratospheric heights [28][29][30]. Kantha and Luce [21] focused on the theoretical derivations of Equation (1) and on the description of this result and its consequences.
The specific goals of this study are: (1) To present the processing methods used for retrieving turbulence parameters ε and C 2 T from T and U data which were also used by Kantha and Luce [21]. The application of Equation (1) requires identification of Kolmogorov turbulence characterized by a −5/3 domain in both one-dimensional (1-D) T and U frequency spectra. (2) To confirm the existence of canonical values of R f and χ m close to those reported in the literature from extended statistics. (3) To analyze the influence of N 2 on these statistics. (4) To assess the relevance of the interpretation of the observed −5/3 domains in terms of Kolmogorov turbulence. (5) To describe a case exhibiting a subtle articulation between the vertical location of the maxima of N 2 , ε and C 2 T . This differs from the case described, for example, by Balsley et al. [16] in that turbulence is not enhanced in a thin stable temperature gradient sheet but at one of its margins, as already described by Dalaudier et al. [7] but using only temperature data.
After a brief description of the experimental set-up (Section 2), we describe in Section 3 the practical method used for identifying bands of the frequency spectra which exhibit a slope consistent with −5/3, a necessary condition for estimating ε and C 2 T . Because a −5/3 domain in 1-D horizontal spectra can also be attributed to non-Kolmogorov regimes especially in strongly stable conditions [13,31,32], this possibility is also tested and discussed. In Section 4, we describe statistical analysis results that tend to confirm a posteriori the hypothesis of a Kolmogorov regime in the selected spectral bands. In Section 5, we describe a case-study and briefly discuss its possible implications to clear air radar backscattering mechanisms, if it is representative of the fine structure of the lower atmosphere. Summary and conclusions are given in Section 6.

Experimental Setup
The sensors on board the DataHawk UAVs were described by Kantha et al. [18] for the ShUREX 2015 campaign and similar configurations were used during ShUREX 2016 and ShUREX 2017. DataHawk is a lightweight and small instrumented fixed wing UAV (about 1.5 m wingspan and about 1 kg mass) developed at the University of Colorado, specifically for in-situ atmospheric measurements (see Figure 1 of Kantha et al. [18]). Its design, the characteristics of ground support components, and some preliminary data collected are described by Balsley et al. [17] and Kantha et al. [18]. Using GPS for navigation, they were usually pre-programmed for moving up and down at a vertical velocity of 2 ms −1 typically, along spiraling ascents and descents in the vicinity of the MU radar. On some occasions, their trajectories were modified in real time, according to the information provided by the MU radar observations, also available in near real time [18]. The DataHawk UAV airspeed ranges from 10 to 20 ms −1 and it can be flown in winds up to 15 ms −1 . Since winds usually increase with height, this was part of the limitation in altitudes of the observations (~up to 4.5 km), in addition to battery life and air traffic regulations.
High-resolution data of the velocity of the flow relative to the DataHawk and temperature were measured by calibrated Pitot-static tube and CWT sensors, respectively. The results were used for estimating turbulence parameters ε and C 2 T . The fast response of the sensors (time constant of 0.5 ms) allowed sampling at an effective rate of 400 Hz and 800 Hz for velocity and temperature measurements, respectively. Pressure and humidity data were also available for each flight at standard (1 Hz) resolution. Most of the flights had long intervals of ascents and descents (denoted by 'A' and 'D' when necessary), separated by horizontal flights of various durations at times. Occasionally, unanticipated blocking of the Pitot tube (used also by the autopilot for flight control) by precipitation produced short time-span (~less than a few tens of seconds) downward motions during ascents. These can be sources of aberrant data points and were manually removed. For the present work, 16 and 23 UAV flights (39 in total) from ShUREX2016 and ShUREX2017 were used. Each flight, hereafter noted FLTXX, where "XX" is the flight number, is composed of one or several ascents and descents between the ground up to variety of altitudes (maximum of 4.5 km). These ascents and descents will be mentioned by an additional 'A' or 'D', respectively, to the flight name.
(a): = −5/3 ± ∆ (b): (a) is satisfied for at least consecutive spectral bands. Then, the corresponding values are selected for use in Equations (6) and (8). After multiple tests, ∆ and were set to 0.25 and 3, respectively. These numerical parameters were chosen in order to fit, as far as possible, the results that would have been obtained from visual inspection of the spectra. In some cases, the criteria may appear too liberal or too restrictive, but they appear to be efficient for rejecting most spectral bands affected by instrumental noise and contaminations as shown below. As the selection procedure is based on the application of thresholds, the results are binary in nature: either or is defined or it is not, depending on the thresholds. When it is not defined, it is either too weak to be selected or the time series are too corrupted by artifacts for retrieving information. (a) "Strong turbulence" Figure 2 shows a case from FLT05A for which the −5/3 domain covered a wide frequency range in both T and U spectra. The spectral slopes for < 2 Hz are rejected because they strongly differ from −5/3, due to a low frequency bump (around 2 Hz) in U spectra and a corresponding dip in the T spectra. For > 2 Hz, the T spectral bands show a tendency close to −5/3 up to about = 40 Hz, but they are rejected because the variability of the estimates exceeds ±0.25. A less stringent threshold would have enabled us to retain more values, but to the detriment of other cases for which the −5/3 tendency is uncertain. Here, spectral bands wider than 0.7 decade would also have provided results that are more satisfying. However, again, it would be to the detriment of cases for

Theoretical Background
The theoretical basis used here for estimating turbulence parameters from spectral analysis was described by Kantha et al. [18]. The most important ones are re-called here. Assuming local isotropy and stationarity of turbulence and using the frozen-turbulence Taylor hypothesis, the theoretical Kolmogorov 1-D power spectral density (spectrum) of the relative wind (airspeed) U is of the form [33,34]: where the coefficient 0.55 holds for motions parallel to the mean streamwise wind U relative the UAV, with U being the mean value of U. Assuming that the calculated spectrumŜ U ( f ) includes a Kolmogorov inertial domain (at least within a frequency band), the spectral data will have the frequency dependenceŜ over the interval of frequency corresponding to the inertial subrange of scales. An estimate of ε can be obtained by calculating the spectral level β U by fitting spectral data, and equating Equations (4) and (5) [35,36]: In the same way, the theoretical Kolmogorov 1-D power spectral density of temperature T is [33]: Using the same procedure as for ε: where β T is the counterpart of β U for temperature spectra.

Practical Methods of Estimating Turbulence Parameters from Pitot and CWT Data
First, T and U frequency spectra are estimated from a variance-conserving, Hanning-weighted Fast Fourier Transform (FFT) using time series of 50 s duration (corresponding to 20,000 and 40,000 points for U and T, respectively). The procedure is applied every 2.5 s, resulting in successive windows overlapping by 95%. Hence consecutive spectral estimates from these time intervals are not independent. However, this has some advantages when comparing profiles (see, e.g., Section 5), and the contribution of remaining isolated aberrant data points in T and U profiles can be more easily detected and rejected.
The selection of frequency bands exhibiting a spectral slope consistent with the presence of a Kolmogorov inertial subrange (~−5/3) is made in a frequency range delimited by 0.1 and 40 Hz (as described below). The lower limit is imposed by the frequency resolution dictated by the time record length (0.02 Hz for 50 s records). The upper limit was selected in order to avoid, as much as possible, artifacts of various origins (mainly motor vibrations) observed in the high frequency range, especially in the airspeed data. These artefacts are not systematic however, and temperature spectra can be free of contaminations at frequencies higher than 40 Hz, as seen later on.
In the frequency range 0.1-40 Hz, 39 overlapping bands of constant width (~0.7 decade, e.g., log 10 (5 (Hz))-log 10 (1 Hz)), shifted by 0.05 decade are used (see Figure 1). Assuming that each band has spectral amplitude in the form β f p , p is estimated by finding the mean spectral amplitude in two spectral sub-bands (noted b1 and b2 in Figure 1) of identical logarithmic width (~0.35 decade). The slope p is given by (log 10 (a2) − log 10 (a1))/(log 10 ( f m2 ) − log 10 ( f m1 )), where a1 and a2 are the mean spectral amplitudes for each of these sub-bands at the mean frequencies f m1 and f m2 . The offset β is obtained by calculating the total power P tot in linear scale, in the band An inertial subrange is inferred when the two criteria below are satisfied: (a): p = −5/3 ± ∆p (b): (a) is satisfied for at least M consecutive spectral bands. Then, the corresponding β values are selected for use in Equations (6) and (8). After multiple tests, ∆p and M were set to 0.25 and 3, respectively. These numerical parameters were chosen in order to fit, as far as possible, the results that would have been obtained from visual inspection of the spectra. In some cases, the criteria may appear too liberal or too restrictive, but they appear to be efficient for rejecting most spectral bands affected by instrumental noise and contaminations as shown below. As the selection procedure is based on the application of thresholds, the results are binary in nature: either C 2 T or ε is defined or it is not, depending on the thresholds. When it is not defined, it is either too weak to be selected or the time series are too corrupted by artifacts for retrieving information. (a) "Strong turbulence" Figure 2 shows a case from FLT05A for which the −5/3 domain covered a wide frequency range in both T and U spectra. The spectral slopes for f max < 2 Hz are rejected because they strongly differ from −5/3, due to a low frequency bump (around 2 Hz) in U spectra and a corresponding dip in the T spectra. For f max > 2 Hz, the T spectral bands show a tendency close to −5/3 up to about f max = 40 Hz, but they are rejected because the variability of the estimates exceeds ±0. 25. A less stringent threshold would have enabled us to retain more values, but to the detriment of other cases for which the −5/3 tendency is uncertain. Here, spectral bands wider than 0.7 decade would also have provided results that are more satisfying. However, again, it would be to the detriment of cases for which the −5/3 domain is much narrower when turbulence is weak. Trade-offs are thus necessary. Here, the non-inclusion of rejected values does not substantially affect the results. For f max > 2 Hz, C 2 T estimates do not vary significantly and the average value is close to the value obtained from all values. By using the additional criterion (b), the final C 2 T estimate is made after rejecting the isolated values at f max = 0.5, 1.58 and 1.99 Hz and by using the groups of 3, 7 and 9 consecutive values, 19 in total. Here, we find C 2 T ∼ 8 × 10 −4 K 2 m −2/3 after applying (b).

Examples of Application
Atmosphere 2019, 10, x FOR PEER REVIEW 6 of 25 which the −5/3 domain is much narrower when turbulence is weak. Trade-offs are thus necessary.
Here, the non-inclusion of rejected values does not substantially affect the results. For > 2 Hz, estimates do not vary significantly and the average value is close to the value obtained from all values. By using the additional criterion (b), the final estimate is made after rejecting the isolated values at = 0.5, 1.58 and 1.99 Hz and by using the groups of 3, 7 and 9 consecutive values, 19 in total. Here, we find ~8 × 10 K m / after applying (b). The spectra shown in Figure 3 and obtained from FLT15A are dominated by instrumental white noise. Both T and U spectra show a slope much steeper than −5/3 at low frequencies and close to 0 at high frequencies. For T (U) spectrum, a −5/3 domain may be observed for 0.7 < < 2 Hz (2 < < 6 Hz). The criterion (b) is met for the U spectrum since there are 5 consecutive values and a very small averaged value of = 10 mWkg is obtained, but not for T spectrum. This is a typical case for which it is difficult to identify with certainty the presence of an inertial domain, because steeper and gentler slopes are observed at low and high frequencies, respectively. The corresponding U spectral bands show a tendency close to −5/3 up to f max ≈ 30 Hz. The decrease of the slope for f max > 30 Hz is due to a narrow and artificial spike just below 40 Hz. The criterion (a) thus permits us to reject the spectral bands corrupted by such contaminations. The average calculated from almost all the values between f max = 2.5 and 30 Hz (= 0.49 mWkg −1 ) provides the estimate of ε for the selected time series.
Atmosphere 2019, 10, 384 7 of 25 (b) "Weak turbulence" The spectra shown in Figure 3 and obtained from FLT15A are dominated by instrumental white noise. Both T and U spectra show a slope much steeper than −5/3 at low frequencies and close to 0 at high frequencies. For T (U) spectrum, a −5/3 domain may be observed for 0.7 < f max < 2 Hz (2 < f max < 6 Hz). The criterion (b) is met for the U spectrum since there are 5 consecutive values and a very small averaged value of ε = 10 −4 mWkg −1 is obtained, but not for T spectrum. This is a typical case for which it is difficult to identify with certainty the presence of an inertial domain, because steeper and gentler slopes are observed at low and high frequencies, respectively. (c) "Strong instrumental contamination" The U spectrum shown in Figure 4 and estimated during FLT15A exhibits obvious instrumental contaminations. The T spectrum was free of artificial spikes. The application of criteria (a) and (b) lead to undefined values of and due to contaminations and instrumental noise, respectively. Figure 4b shows an unusual case, however. The 50 s time series of U and T show dips and steps probably due to a steep bank and sudden loss of altitude. However, these fluctuations are not the cause of the artificial spikes above 10 Hz in the U spectra.

Attempt at Interpretation of the −5/3 Subranges in Terms of Kolmogorov Turbulence
The consistency of the observed −5/3 domain with the assumptions for Kolmogorov inertial range theory can be inferred from the estimation of the Ozmidov length scale . It is defined as: Its inverse defines a wavenumber indicative of the separation between the Kolmogorov inertial and buoyancy subranges. The corresponding Ozmidov frequency can be defined as / (2 ). It is expected to be indicative of the lower limit of the inertial subrange in the U frequency spectrum. Its definition differs from that given by Frehlich et al. [13] and Kit et al. [14] for example, by a factor 2 . Here, this factor is justified by the fact that is not a wavelength, but a scale, and that = (2 ⁄ ) = 1/ . However, the results discussed here are not strongly dependent on this (c) "Strong instrumental contamination" The U spectrum shown in Figure 4 and estimated during FLT15A exhibits obvious instrumental contaminations. The T spectrum was free of artificial spikes. The application of criteria (a) and (b) lead to undefined values of ε and C 2 T due to contaminations and instrumental noise, respectively. Figure 4b shows an unusual case, however. The 50 s time series of U and T show dips and steps probably due to a steep bank and sudden loss of altitude. However, these fluctuations are not the cause of the artificial spikes above 10 Hz in the U spectra.  Hz. The portion of the spectrum exhibiting a −5/3 slope was found at lower frequencies ( 2 < < 6 Hz) in the U spectrum and is within a frequency domain dominated by the instrumental noise. Therefore, it is not consistent with Kolmogorov inertial turbulence. This result raises issues about the interpretation of the −5/3 domain in this case. If a −5/3 slope is not a signature of Kolmogorov turbulence (as discussed in 3.3a), would be overestimated and would be underestimated, reinforcing the conclusion. If it is relevant, another regime must be considered. Stratified turbulence at low Froude number ( = ≪ 1 ⁄ , where is a velocity scale and L the corresponding horizontal velocity scale) consisting of quasi-horizontal turbulent motions may lead to the presence of inertial subranges in horizontal spectra at scales larger than the Ozmidov scale [37][38]. In stable boundary layers, dynamical processes strongly affected by the stable stratification can lead to various regimes producing −5/3 or shallower slopes at scales larger or smaller than the Ozmidov

Attempt at Interpretation of the −5/3 Subranges in Terms of Kolmogorov Turbulence
The consistency of the observed −5/3 domain with the assumptions for Kolmogorov inertial range theory can be inferred from the estimation of the Ozmidov length scale L O . It is defined as: Its inverse k O defines a wavenumber indicative of the separation between the Kolmogorov inertial and buoyancy subranges. The corresponding Ozmidov frequency f Oz can be defined as U/(2πL O ). It is expected to be indicative of the lower limit of the inertial subrange in the U frequency spectrum. Its definition differs from that given by Frehlich et al. [13] and Kit et al. [14] for example, by a factor 2π. Here, this factor is justified by the fact that L O is not a wavelength, but a scale, and that k O = (2π/U) f Oz = 1/L O . However, the results discussed here are not strongly dependent on this factor. f Oz is estimated from ε using the method described in Section 3.2 and from N 2 estimated as follows: for an ascent or descent rate of~2 ms −1 , each temperature data segment corresponds to a vertical segment ∆z of (~2x50 s) = 100 m. The corresponding temperature variation ∆T is estimated after a linear interpolation of the measured temperature profile. Then: where T is the mean temperature within the segment and c p = 1004 Jkg −1 K −1 . The same procedure is used when assessing the relevance of Equation (1) in Section 4. Equation (10) is relevant if |W/W UAV | 1 where W and W UAV are the vertical air and UAV velocities, respectively (see Appendix A). During the UAV flights, this condition was mostly fulfilled in the free atmosphere according to MU radar measurements of vertical air velocities. However, it must be kept in mind that vertical air motions can be a source of errors on N 2 that should increase the dispersion of the statistical distributions of the turbulence parameters depending on N 2 and shown in Section 4.
The values of f Oz are indicated on Figures 2 and 3. For the "Strong turbulence" case ( Figure 2), Hz. The −5/3 behavior observed for higher frequencies is thus consistent with the hypothesis of a Kolmogorov inertial subrange. For the "Weak turbulence" case ( Figure 3 The portion of the spectrum exhibiting a −5/3 slope was found at lower frequencies (2 < f max < 6 Hz) in the U spectrum and f Oz is within a frequency domain dominated by the instrumental noise. Therefore, it is not consistent with Kolmogorov inertial turbulence. This result raises issues about the interpretation of the −5/3 domain in this case. If a −5/3 slope is not a signature of Kolmogorov turbulence (as discussed in 3.3a), ε would be overestimated and f Oz would be underestimated, reinforcing the conclusion. If it is relevant, another regime must be considered. Stratified turbulence at low Froude number (F r = σ U /NL 1, where σ U is a velocity scale and L the corresponding horizontal velocity scale) consisting of quasi-horizontal turbulent motions may lead to the presence of inertial subranges in horizontal spectra at scales larger than the Ozmidov scale [37,38]. In stable boundary layers, dynamical processes strongly affected by the stable stratification can lead to various regimes producing −5/3 or shallower slopes at scales larger or smaller than the Ozmidov scale [14,38]. The identification of −5/3 domain associated with non-Kolmogorov regime may lead to erroneous estimates of ε and C 2 T , unless there is no spectral gap between the two regimes (as in the case reported by Frehlich et al. [13]).
The statistical results described later on in Section 4 do not seem to indicate a posteriori a noticeable misinterpretation with non-Kolmogorov inertial regimes.   are shown at the right side of each plot. For better legibility, only three levels of colors are shown corresponding to three ranges of slopes: p < −5/3 − 0.5 (blue), −5/3 − 0.25 < p < −5/3 + 0.25 (cyan) and p > −5/3 + 0.5 (red), with a gap of 0.25 between each range for a clear distinction. The steep slopes displayed by the blue levels at low frequencies (for f max < 3 Hz) may have an instrumental (unknown) origin and should not be considered. The morphology and intensity of the corresponding spikes strongly depend on the UAV flights (not shown).

Additional Evaluation of the Processing Methods
As expected, the spectral extent of the inertial domains is large when C 2 T and ε values are high, as emphasized by black vertical double arrows in Figure 5. Most of the red levels at the right side of the 2-D plots, labeled 'N', should be understood as the effect of the instrumentation noise. For U spectra and during FLT05A, a spectral contamination by narrow spikes affects the spectra down to 20 Hz. An example was shown in Figure 2 where a spike appears below 40 Hz at the altitude of 1635 m. The contaminations drift toward lower frequencies with time/altitude. The contaminated region is delimited by a thick dashed line and is labeled 'A' in Figure 5. In the top-right-hand corner of the 2-D plot, the blue region (steep slopes) indicates that the spectral band was on the right side of the spike. The criteria used permitted us to avoid this contamination when estimating ε. C 2 T and ε values are low or not defined, when the inertial domains are found to be very narrow: the 'red' domain is largely extended toward the low frequencies. As the CWT data are relatively noisier than the Pitot data, the C 2 T profiles usually exhibit more undefined values than the ε profiles. The values from the band 5-40 Hz are often strongly overestimated due to instrumentation noise (e.g., around 3000 m during descent for ε and above~1500 m for C 2 T during ascent and descent) and due to artifacts (above 1700 m during ascent for ε). FLT05 reveals substantial turbulence at all altitudes. Other flights (as FLT28, see Section 4.3) showed weak turbulence and many more estimates were undetermined after application of the criteria (a) and (b).  T and ε profiles after the selection criteria in black lines, when selecting 1−5 Hz (gray lines) and 5−40 Hz (gray dotted lines). The labels "N" and "A" refer to spectral bands affected by noise and artifacts, respectively (see text for more details).
The frequencies of occurrence of ranges of slopes for T (black) and U (red) spectra are shown in Figure 6. The largest occurrences of −5/3 − 0.25 < p < −5/3 + 0.25 (a qualified inertial domain) are found for 5 < f max < 8 Hz, i.e., for frequency bands of 1-5 Hz to 1.6-8 Hz (~35 % for T spectra, 50 % for U spectra). This frequency range was used by Luce et al. [19] for estimating ε from Pitot data to be compared with MU radar estimates of ε. Above f max ≈ 8 Hz, the frequencies of occurrence slowly decrease down to 20% and 35% at f max = 40 Hz for T and U spectra, respectively, because of instrumentation noise and artifacts. The frequencies of occurrence of p > −5/3 + 0.5 (spectra affected by noise and artifacts) monotonically increase with f max . They reach 60% and 45% at f max = 40 Hz and should still increase for higher frequencies. Even if the frequency of occurrences of p < −5/3 − 0.5 (spectra affected by artifacts apart from the low frequencies) does not exceed~5% (with a slight increase with as f max increases), the quantitative effects can be large because the spikes are generally very energetic (see Figure 5). It is thus important to remove their contribution as much as possible. > −5/3 + 0.5 (cyan), −5/3 − 0.25 < < −5/3 + 0.25 , (blue): < −5/3 − 0.5 , and the corresponding and profiles after the selection criteria in black lines, when selecting 1−5 Hz (gray lines) and 5−40 Hz (gray dotted lines). The labels "N" and "A" refer to spectral bands affected by noise and artifacts, respectively (see text for more details) Figure 6. Frequency of occurrences of ranges of slopes vs for T spectra (black) and U spectra (red) Figure 6. Frequency of occurrences of ranges of slopes vs. f max for T spectra (black) and U spectra (red). (2) can be estimated from Equation (1):

Results of Analyses
where N 2 is estimated using Equation (10). ε and C 2 T values associated with convectively generated turbulence and cloudy air must be, in principle, discarded. This was done when: (1) Relative humidity measured by the humidity sensor on board the UAV exceeded 85%. This threshold is rather arbitrary and may be affected by inaccurate calibration of the humidity sensors. It may lead to discarded data in clear air conditions under some occasions, without impact on the statistical results. In practice, cloud particles can often be observed even when the relative humidity is below but close to 100%, depending on altitudes, and radiosondes [39].
(2) The data were associated with convective boundary layers (CBL). These often extended up to the altitude of 2000-2500 m at midday during sunny periods. They were easily identified from MU radar observations, when CBL top reached the minimum observable altitude, i.e., 1270 m (for example, Figure 23 of Kantha et al. [18]). From temperature and humidity data, CBL is usually weakly stratified and capped by a sharp temperature inversion associated with a relative humidity decrease. In addition, CBL is generally associated with the largest values of ε and C 2 T . All these features helped us identify the vertical extent of CBL.
(3) Turbulence was caused by convective instabilities underneath mid-level cloud base (MCT) [40]. Very few such cases were probed, because the MCT often occurred at altitudes unreachable by the UAVs.
These criteria disqualified 24.7% of the estimates made from the 39 flights. The remaining data are expected to fit the assumptions underlying Equation (1), because the main sources of convection due to thermodynamic effects were excluded. Among the discarded data, an overwhelming majority comes from the CBL, which was sometimes partly cloudy. In the following, we will attribute all the discarded data to CBL for simplicity. Some of the estimates for CBL data were de facto discarded when N 2 was locally negative so that Equations (1) or (11) cannot be applied. The temperature profiles in CBL sometimes revealed suspicious, successive and almost periodic superadiabatic and stable gradients. Their origin is unclear but, for at least some of the observed cases, they could be produced by the possible presence of strong updrafts or downdrafts in the convective cells (see Appendix A). A possible similar effect reported in a review by Gossard [41] was observed from tower measurements during strong buoyancy wave perturbations.
Without CBL data, the distribution of frequency of occurrences of log 10 (γ) obtained from the 39 flights has a nearly log-normal distribution with identical mean and median values of 0.28 for log 10 (γ) and 1.92 for γ (Figure 7). This value is close to 1.95 obtained from ShUREX 2016 data only [21]. The standard deviation σ log 10 (γ) of the distribution is 0.43 and the percentage of values that lie within the band < log 10 (γ) > ± σ log 10 (γ) is 75% (68% for a normal distribution). Therefore, the distribution clearly indicates the existence of a preferential value of γ around 1.92. This value significantly differs from typical values sometimes used by default in the literature of atmospheric studies. γ ≈ 1.0 was used by Gavrilov et al. [3] based on the assumption of R f = 0.25. Using slightly different models, γ ≈ 0.4 was proposed to be used by Hocking [34]. This is smaller than our estimate by a factor~5 and would produce a factor~11 between ε and ε CT2 .

Estimations of , and
We compared to by using < ( ) >= 0.28 with estimated from Equation (10), and < >= 1.47 × 10 rad s (Figure 8a). The constant value < > is the average of for selected height ranges in stratified and clear air conditions. The scatter plots show and vary between ~10 m s and ~10 m s . The regression slope is close to 1 for both measured and averaged (1.11 and 1.17, respectively). It slightly exceeds 1 partly because of a slight asymmetry of the scatter plot for weak values of and (≤ 10 m s ). Rather than a different property of turbulence for weak values, this feature may result from imperfections of the rejection/selection procedure described in Section 3 and should probably not be taken into account. The scatter is more important with < > than with and the correlation coefficients are 0.74 and 0.81, respectively. An increase of the dispersion and a decrease of the correlation when using < > uggest that the estimates from Equation (10) are meaningful and are thus more representative than the average value. The use of an average or standard value cause significant disagreements for individual cases as shown later on, in Section 5.
The histograms of mixing coefficient and flux Richardson number (in scale) cover about two orders of magnitude (Figures 8b,c). Their mean and median values obtained from and < > do not differ significantly: 0.16-0.17 and 0.13-0.15, respectively. As mentioned by [21], these values are close to the canonical values of = 0.2 and = 0.17 reported in the literature on oceanic turbulence [22] and theorized by Ellison [23]. The fact that the upper tail of the distribution of ( ) exceeds 0 (i.e., exceeds 1) indicates that the distribution is affected by statistical estimation errors, in particular on , that make it even wider than it should be. Considering that these distributions are also conditioned by the validity of Equation (1), it is not possible to infer the real variability of and . Interestingly, the spread of the distribution of obtained for oceanic turbulence from different methods and with accurate evaluations of parameters [22] is very similar to that shown in Figure 8c. Including CBL data, the distribution is slightly shifted and skewed to the right (Figure 7). The mean and median values of log 10 (γ) are 0.34 and 0.32, respectively. The distribution is more spread out than the distribution obtained excluding the CBL (σ log 10 (γ) = 0.48). Therefore, CBL data make a significant contribution to the statistics and can be a source of bias when evaluating γ. In the following, we will focus on the statistics made from data outside the CBL. A detailed analysis of the relationship between ε and C 2 T in convective and cloudy regions remains the subject of future studies.

Estimations of ε CT2 , χ m and R f
We compared ε to ε CT2 by using < log 10 (γ) >= 0.28 with N 2 estimated from Equation (10), and < N 2 > = 1.47 × 10 −4 rad 2 s −2 (Figure 8a). The constant value < N 2 > is the average of N 2 for selected height ranges in stratified and clear air conditions. The scatter plots show ε and ε CT2 vary between ∼ 10 −7 m 2 s −3 and ∼ 10 −2 m 2 s −3 . The regression slope is close to 1 for both measured and averaged N 2 (1.11 and 1.17, respectively). It slightly exceeds 1 partly because of a slight asymmetry of the scatter plot for weak values of ε and ε CT2 (≤ 10 −6 m 2 s −3 ). Rather than a different property of turbulence for weak values, this feature may result from imperfections of the rejection/selection procedure described in Section 3 and should probably not be taken into account. The scatter is more important with < N 2 > than with N 2 and the correlation coefficients are 0.74 and 0.81, respectively. An increase of the dispersion and a decrease of the correlation when using < N 2 > suggest that the N 2 estimates from Equation (10) are meaningful and are thus more representative than the average value. The use of an average or standard N 2 value cause significant disagreements for individual cases as shown later on, in Section 5.

Examples of Profiles
We inferred (pseudo) vertical profiles of turbulence parameters and some examples are shown in Figures 9, 10 and 11 for FLT05A1, FLT28A1 and FLT40A1, respectively. We note the Ozmidov scale , ( ) and the eddy diffusivity , ( ) inferred from and , respectively. and ( ) were estimated from Equation (3) using < /(1 − ) >=< >= 0.16. This value is smaller by a factor 2 with respect commonly used formulas assuming = 0.25 [2,24] and by a factor 5 with respect to Hocking [42] and a factor 2.5 with respect to Hocking [43].
The altitude ranges associated with CBL in Figures 9 and 10 show the largest discrepancies between the profiles because Equation (1) should not apply. On the day of FLT40, CBL did not form but cloudy air (relative humidity=100%) was observed up to ~1900 m ( Figure 11). In the height range 1200−1800 m, ( ( ) , ( ) ) strongly underestimate ( , ) by up to two orders of magnitude, whereas there is a good agreement below 1200 m. This observation might indicate that air saturation is not always a sufficient criterion for faulting the validity of Equation (1). The discrepancy in the height range 1200-1800 m is clearly not due to statistical errors and is one of the rare prominent cases not related to CBL at odds with the model. The hypothesis of vanishing temperature turbulence due to strong mixing [44] for stratospheric heights (their Figures 1 and 2) can be discarded because this hypothesis implies enhancements of and at the edges of the mixing layer and a decrease of at the center. These characteristics are not observed here (Figure 11b). Above CBL and outside cloudy regions, Equation (1) with < ( ) >= 0.28 is applicable because the profiles are consistent with each other. The largest discrepancies result mainly from differences between the amplitudes of the peaks (e.g., at the altitude of 1700 m in Figure 9 where the ratio is ≈ 10 for the 3 parameters). In Figure 9, the profiles reveal an alternation of peaks of dissipation rates, whose depths vary from 100 m or less (at 3200 m) up to 600 m around 2600 m (hereafter denoted T1). In the latter case, and reached 10 , and In Figure 10, turbulence was weak above the CBL, and only a few isolated and thin layers were probed (data were missing around 1500 m). The depth of the peaks did not exceed ~100 m typically. and values did not always coincide, but there were likely significant peaks at 2000, 2600, 3700 and 4000 m with and ~10 m s , and  The histograms of mixing coefficient χ m and flux Richardson number R f (in log 10 scale) cover about two orders of magnitude (Figure 8b,c). Their mean and median values obtained from N 2 and < N 2 > do not differ significantly: 0.16-0.17 and 0.13-0.15, respectively. As mentioned by [21], these values are close to the canonical values of χ m = 0.2 and R f = 0.17 reported in the literature on oceanic turbulence [22] and theorized by Ellison [23]. The fact that the upper tail of the distribution of log 10 (χ m ) exceeds 0 (i.e., χ m exceeds 1) indicates that the distribution is affected by statistical estimation errors, in particular on N 2 , that make it even wider than it should be. Considering that these distributions are also conditioned by the validity of Equation (1), it is not possible to infer the real variability of χ m and R f . Interestingly, the spread of the distribution of R f obtained for oceanic turbulence from different methods and with accurate evaluations of parameters [22] is very similar to that shown in Figure 8c

Examples of Profiles
We inferred (pseudo) vertical profiles of turbulence parameters and some examples are shown in Figures 9-11 for FLT05A1, FLT28A1 and FLT40A1, respectively. We note the Ozmidov scale L O , L O(CT2) and the eddy diffusivity K θ , K θ(CT2) inferred from ε and ε CT2 , respectively. K θ and K θ(CT2) were estimated from Equation (3) using < R f / 1 − R f > = < χ m > = 0.16. This value is smaller by a factor 2 with respect commonly used formulas assuming R f = 0.25 [2,24] and by a factor 5 with respect to Hocking [42] and a factor 2.5 with respect to Hocking [43].
The altitude ranges associated with CBL in Figures 9 and 10 show the largest discrepancies between the profiles because Equation (1) should not apply. On the day of FLT40, CBL did not form but cloudy air (relative humidity=100%) was observed up to~1900 m ( Figure 11). In the height range 1200-1800 m, ε CT2 (K θ(CT2) , L O(CT2) ) strongly underestimate ε (K θ , L O ) by up to two orders of magnitude, whereas there is a good agreement below 1200 m. This observation might indicate that air saturation is not always a sufficient criterion for faulting the validity of Equation (1). The discrepancy in the height range 1200-1800 m is clearly not due to statistical errors and is one of the rare prominent cases not related to CBL at odds with the model. The hypothesis of vanishing temperature turbulence due to strong mixing [44] for stratospheric heights (their Figures 1 and 2) can be discarded because this hypothesis implies enhancements of C 2 T and N 2 at the edges of the mixing layer and a decrease of N 2 at the center. These characteristics are not observed here (Figure 11b show the CBL depth. "T1" refers to a turbulent layer discussed in the text. In Figure 11, the noticeable feature is the thick layer of enhanced turbulence parameters between 2200 and 3000 m (hereafter denoted T2). and In T1 and T2, and are both enhanced. This differs from that reported by Bertin et al. [44] for stratospheric layers because was enhanced where was minimum. The difference may suggest an earlier stage of turbulence and/or another model of evolution. A careful inspection of the overall profiles did not reveal events similar to those shown by Bertin et al. [44] in stratified conditions. Figure 10. The same as Figure 9 for FLT28A1. The gray shaded areas show the CBL depth. Above CBL and outside cloudy regions, Equation (1) with < log 10 (γ) >= 0.28 is applicable because the profiles are consistent with each other. The largest discrepancies result mainly from differences between the amplitudes of the peaks (e.g., at the altitude of 1700 m in Figure 9 where the ratio is ≈ 10 for the 3 parameters).
In Figure 9, the profiles reveal an alternation of peaks of dissipation rates, whose depths vary from 100 m or less (at 3200 m) up to 600 m around 2600 m (hereafter denoted T1). In the latter case, ε and ε CT2 reached 10 −3 m 2 s −3 , K θ and K θ(CT2) exceeded 1. show the CBL depth. "T1" refers to a turbulent layer discussed in the text.
In Figure 11, the noticeable feature is the thick layer of enhanced turbulence parameters between 2200 and 3000 m (hereafter denoted T2). and In T1 and T2, and are both enhanced. This differs from that reported by Bertin et al. [44] for stratospheric layers because was enhanced where was minimum. The difference may suggest an earlier stage of turbulence and/or another model of evolution. A careful inspection of the overall profiles did not reveal events similar to those shown by Bertin et al. [44] in stratified conditions. Figure 10. The same as Figure 9 for FLT28A1. The gray shaded areas show the CBL depth. Figure 10. The same as Figure 9 for FLT28A1. The gray shaded areas show the CBL depth.  Figure 11. The same as Figure 9 for FLT40A1. Here, the gray shaded areas does not show CBL range but the height range where air was cloudy (saturated). "T2" refers to a turbulent layer discussed in the text.
These values are consistent with those reported from balloon data collected at stratospheric heights [29,44]. However, 6.7% of values exceeded 1.0 m s (e.g., T1 and T2 of Figures 9 and 11), indicating that some of the observed turbulent events may have had a much more significant impact in terms of diffusion than those reported in the literature [29,44].   Figure 9 for FLT40A1. Here, the gray shaded areas does not show CBL range but the height range where air was cloudy (saturated). "T2" refers to a turbulent layer discussed in the text.
In Figure 10, turbulence was weak above the CBL, and only a few isolated and thin layers were probed (data were missing around 1500 m). The depth of the peaks did not exceed~100 m typically. ε and ε CT2 values did not always coincide, but there were likely significant peaks at 2000, 2600, 3700 and 4000 m with ε and ε CT2~1 0 −5 m 2 s −3 , K θ and K θ(CT2) ∼ 10 −2 − 10 −1 m 2 s −1 and L O and L O(CT2) ∼ 1 − 10 m.
In Figure 11, the noticeable feature is the thick layer of enhanced turbulence parameters between 2200 and 3000 m (hereafter denoted T2). ε and ε CT2 reached 10 −3 m 2 s −3 , K θ and K θ(CT2) In T1 and T2, C 2 T and ε are both enhanced. This differs from that reported by Bertin et al. [44] for stratospheric layers because ε was enhanced where C 2 T was minimum. The difference may suggest an earlier stage of turbulence and/or another model of evolution. A careful inspection of the overall profiles did not reveal events similar to those shown by Bertin et al. [44] in stratified conditions.

Statistics for K θ and L O
< K θ > = 0.03 m 2 s −1 and < L O > ∼ 4 m was obtained for both N 2 and < N 2 > (Figure 12). These values are consistent with those reported from balloon data collected at stratospheric heights [29,44]. However, 6.7% of K θ values exceeded 1.0 m 2 s −1 (e.g., T1 and T2 of Figures 9 and 11), indicating that some of the observed turbulent events may have had a much more significant impact in terms of diffusion than those reported in the literature [29,44]. Figure 13 shows distributions of log 10 (L O ) and log 10 L O(CT2) and a scatter plot of log 10 (L O ) vs. < >= 0.03 m s and < >~4 m was obtained for both and < > (Figure 12).
These values are consistent with those reported from balloon data collected at stratospheric heights [29,44]. However, 6.7% of values exceeded 1.0 m s (e.g., T1 and T2 of Figures 9 and 11), indicating that some of the observed turbulent events may have had a much more significant impact in terms of diffusion than those reported in the literature [29,44].  Finally, Figure 14 shows the distribution of the buoyancy Reynolds number where is the molecular kinematic viscosity and is the Kolmogorov length scale. The average value of ( ) is ~4.0 (i.e., the average value of is ~10,000) and 95% of the estimates fulfill the condition > 200, consistent with isotropic turbulence according to [45] for oceanic turbulence. It is an additional clue that our qualified estimates of are meaningful and related to Kolmogorov turbulence.

Shear and Richardson Number
A clear relationship between enhanced turbulence parameters and both enhanced shear = ( ) ⁄ + ( ) ⁄ / where u and v are the zonal and meridional wind components and small gradient Richardson number = / was reported from UAV data up to altitudes of ~500 m [17]. The estimation of wind shear and from UAV measurements is still under consideration for ShUREX data. However, these parameters can be obtained from RS92G Vaisala radiosondes launched Finally, Figure 14 shows the distribution of the buoyancy Reynolds number Re b = ε/ νN 2 = (L O /η) 4/3 , in log 10 scale, where ν is the molecular kinematic viscosity and η is the Kolmogorov length scale. The average value of log 10 (Re b ) is~4.0 (i.e., the average value of Re b is~10,000) and 95% of the estimates fulfill the condition Re b > 200, consistent with isotropic turbulence according to [45] for oceanic turbulence. It is an additional clue that our qualified estimates of ε are meaningful and related to Kolmogorov turbulence. Finally, Figure 14 shows the distribution of the buoyancy Reynolds number where is the molecular kinematic viscosity and is the Kolmogorov length scale. The average value of ( ) is ~4.0 (i.e., the average value of is ~10,000) and 95% of the estimates fulfill the condition > 200, consistent with isotropic turbulence according to [45] for oceanic turbulence. It is an additional clue that our qualified estimates of are meaningful and related to Kolmogorov turbulence.

Shear and Richardson Number
A clear relationship between enhanced turbulence parameters and both enhanced shear = ( ) ⁄ + ( ) ⁄ / where u and v are the zonal and meridional wind components and small gradient Richardson number = / was reported from UAV data up to altitudes of ~500 m [17]. The estimation of wind shear and from UAV measurements is still under consideration for ShUREX data. However, these parameters can be obtained from RS92G Vaisala radiosondes launched

Shear and Richardson Number
A clear relationship between enhanced turbulence parameters and both enhanced shear S = ((du/dz) 2 + (dv/dz) 2 ) 1/2 where u and v are the zonal and meridional wind components and small gradient Richardson number Ri = N 2 /S 2 was reported from UAV data up to altitudes of~500 m [17]. The estimation of wind shear and Ri from UAV measurements is still under consideration for ShUREX data. However, these parameters can be obtained from RS92G Vaisala radiosondes launched on some occasions from the observatory concurrently with the UAVs. Two radiosondes (hereafter called V06 and V09) were launched around the flights FLT28 and FLT40 described in Section 4.4. V06 was launched about 30 min before FLT28 and drifted horizontally by 5600 m at the maximum altitude of FLT28. V09 flew almost simultaneously to FLT40. V09 drifted horizontally by about 7700 m at the maximum altitude of FLT40.
We first compared N 2 profiles derived from radiosonde and UAV data, hereafter denoted N 2 (V * * ) and N 2 (FLT * * ), respectively, where "**" is the flight number, for checking their consistency (Figures 15  and 16). N 2 (V06) and N 2 (FLT28) compare quite well (Figure 15b). Altitudes of V09 data were shifted downward by 130 m in order to correct a vertical offset between the most prominent peak of N 2 (V09) and N 2 (FLT40) at~2150 m. After this correction, N 2 (V09) and N 2 (FLT40) profiles also show very similar features (Figure 16b). We thus assume that the shear and Ri values estimated from V06 and V09 data are also representative, at least qualitatively, of the dynamical conditions met during FLT28A1 and FLT40A1. and values estimated from V06 and V09 data are also representative, at least qualitatively, of the dynamical conditions met during FLT28A1 and FLT40A1.
For FLT28A1, the profile fluctuates with height with alternating large values ( > 10) and small values ( < 1) in the stratified region above the altitude of 1700 m. It is qualitatively consistent with the sporadic peaks of and (Figure 15d). The narrow peak of and around ~2150 m (Figure 15d) is associated with < 0.25 ( Figure 15c) and a shear maximum of ~20 ( Figure 15a). For FLT40A1, the maxima of (> 10) at 1800, 2200 and 3000 m coincide with the minima of and or even undefined values (Figure 16d). The thick turbulent layer (T2 in Figure 11) with enhanced and (Figure 16d) coincides with a deep minimum of ~< 0.25 ( Figure 16c) and shear maximum of ~20 ms km ( Figure 16a). All these observations are consistent with the commonly accepted criteria on gradient Richardson number for turbulence generated by dynamic shear instabilities and confirm results described by Balsley et al. [17]. It is again an additional clue supporting the hypothesis that we truly selected 3-D Kolmogorov turbulence rather than stratified turbulence.

Analysis of a Data Segment
Here, we mainly describe a case-study selected from FLT05D1 ( Figure 5). The profiles of and show noticeable similarities and differences in the height range 1800-3500 m ( Figure 17). Below the altitude of 2400 m, several peaks of exceed 5 × 10 K m / and generally exceeds the mean value of 1.47 × 10 rad s with peaks up to 8 × 10 rad s . Above 2800 m, where is smaller but exceeds 10 rad s in some altitude ranges, is minimum or even undetermined. Such properties in both the troposphere and stratosphere from high-resolution balloon temperature data were reported by Gavrilov et al. [3]. The peaks labeled A, B and C are closely associated with peaks. Note how A, B and C differ from T1 ( Figure 9) and T2 ( Figure 11). The depths of A, B and C are less than ~100 m whereas T1 and T2 are deep (600-800 m) and not associated with a welldefined enhancement. A careful inspection of the profiles shows that the altitudes of and peaks do not coincide exactly. The differences of altitude, − for A, B and C are −10 m, −9 m and + 35 m, respectively. The other peaks in Figure 17 have not been selected, because the peak association may be quite arbitrary due to multiple consecutive peaks. In order to check the significance of these vertical offsets, 134 additional cases have been selected from the 39 flights. Great care has been taken to select isolated peaks such as A, B and C. A large majority of the selected cases show a positive or negative offset between and peaks of the order of ±(10 − 30) m ( Figure 18). Therefore, the largest temperature fluctuations attributed to Kolmogorov turbulence did not occur at the center of the temperature inversion but at one of its margins, and never at both margins. From the selected cases, it was not possible to associate unambiguously a peak at one margin with another one at the other margin to form a pair of layers of enhanced turbulence at the edges of a more weakly turbulent layer. For FLT28A1, the Ri profile fluctuates with height with alternating large values (Ri > 10) and small values (Ri < 1) in the stratified region above the altitude of 1700 m. It is qualitatively consistent with the sporadic peaks of ε and ε CT2 (Figure 15d). The narrow peak of ε and ε CT2 around~2150 m (Figure 15d) is associated with Ri < 0.25 ( Figure 15c) and a shear maximum of ∼ 20 ms −1 km −1 (Figure 15a). For FLT40A1, the maxima of Ri (> 10) at 1800, 2200 and 3000 m coincide with the minima of ε and ε CT2 or even undefined values (Figure 16d). The thick turbulent layer (T2 in Figure 11) with enhanced ε and ε CT2 (Figure 16d) coincides with a deep minimum of Ri ∼< 0.25 ( Figure 16c) and shear maximum of ∼ 20 ms −1 km −1 (Figure 16a). All these observations are consistent with the commonly accepted criteria on gradient Richardson number for turbulence generated by dynamic shear instabilities and confirm results described by Balsley et al. [17]. It is again an additional clue supporting the hypothesis that we truly selected 3-D Kolmogorov turbulence rather than stratified turbulence.

Analysis of a Data Segment
Here, we mainly describe a case-study selected from FLT05D1 ( Figure 5). The profiles of C 2 T and N 2 show noticeable similarities and differences in the height range 1800-3500 m ( Figure 17). Below the altitude of 2400 m, several peaks of C 2 T exceed 5 × 10 −4 K 2 m −2/3 and N 2 generally exceeds the mean value of 1.47 × 10 −4 rad 2 s −2 with peaks up to 8 × 10 −4 rad 2 s −2 . Above 2800 m, where N 2 is smaller but exceeds 10 −4 rad 2 s −2 in some altitude ranges, C 2 T is minimum or even undetermined. Such properties in both the troposphere and stratosphere from high-resolution balloon temperature data were reported by Gavrilov et al. [3]. The C 2 T peaks labeled A, B and C are closely associated with N 2 peaks. Note how A, B and C differ from T1 ( Figure 9) and T2 ( Figure 11). The depths of A, B and C are less than~100 m whereas T1 and T2 are deep (600-800 m) and not associated with a well-defined N 2 enhancement. A careful inspection of the profiles shows that the altitudes of N 2 and C 2 T peaks do not coincide exactly. The differences of altitude, z CT2max − z N2max for A, B and C are −10 m, −9 m and + 35 m, respectively. The other peaks in Figure 17 have not been selected, because the peak association may be quite arbitrary due to multiple consecutive peaks. In order to check the significance of these vertical offsets, 134 additional cases have been selected from the 39 flights. Great care has been taken to select isolated peaks such as A, B and C. A large majority of the selected cases show a positive or negative offset between N 2 and C 2 T peaks of the order of ±(10 − 30) m ( Figure 18). Therefore, the largest temperature fluctuations attributed to Kolmogorov turbulence did not occur at the center of the temperature inversion but at one of its margins, and never at both margins. From the selected cases, it was not possible to associate unambiguously a C 2 T peak at one margin with another one at the other margin to form a pair of layers of enhanced turbulence at the edges of a more weakly turbulent layer.   Figure 19 focuses on C, in the height range 2550−3000 m. It shows , , with and < > and profiles and the corresponding profile of potential temperature . The peak is wider than the peak and is shifted upward with respect to the peak. and obviously coincide when < > is used and roughly agrees with . However, when using , the agreement is much better. This particular case illustrates the cause of the statistical improvement when estimating from with respect to < > ( Figure 8). The maximum of (or ) above the maximum of and a maximum of between the two is remarkably similar to observations reported by Gossard [41] (Figure 4.5, p. 519), at the bottom of temperature inversions capping the CBL. Thus, mechanisms similar to those observed below capping inversions may occur in the free atmosphere.
The vertical offsets between , , and peaks can be interpreted as follows. On the one hand, peak of (maximum of Kolmogorov turbulence) is found to be far enough from the center of the stable inversion because the maximum of static stability inhibits vertical motion (and thus isotropic turbulence). The largest isotropic billows (maximum of ) can be observed there. They tend to be smaller when "approaching" the stable inversion, i.e., the Kolmogorov domain tends to be suppressed. On the other hand, large requires both enhanced and . When ( ) is maximum, ( ) is reduced so that is also reduced. It follows that is maximum between and peaks. The level of is thus a compromise between the intensity of the dynamic turbulence and the strength of the mean temperature gradient. The potential temperature profile clearly indicates a steep gradient with a laminar appearance at ~2700 m topped by enhanced fluctuations   Figure 19 focuses on C, in the height range 2550−3000 m. It shows , , with and < > and profiles and the corresponding profile of potential temperature . The peak is wider than the peak and is shifted upward with respect to the peak. and obviously coincide when < > is used and roughly agrees with . However, when using , the agreement is much better. This particular case illustrates the cause of the statistical improvement when estimating from with respect to < > ( Figure 8). The maximum of (or ) above the maximum of and a maximum of between the two is remarkably similar to observations reported by Gossard [41] (Figure 4.5, p. 519), at the bottom of temperature inversions capping the CBL. Thus, mechanisms similar to those observed below capping inversions may occur in the free atmosphere.
The vertical offsets between , , and peaks can be interpreted as follows. On the one hand, peak of (maximum of Kolmogorov turbulence) is found to be far enough from the center of the stable inversion because the maximum of static stability inhibits vertical motion (and thus isotropic turbulence). The largest isotropic billows (maximum of ) can be observed there. They tend to be smaller when "approaching" the stable inversion, i.e., the Kolmogorov domain tends to be suppressed. On the other hand, large requires both enhanced and . When ( ) is maximum, ( ) is reduced so that is also reduced. It follows that is maximum between and peaks. The level of is thus a compromise between the intensity of the dynamic turbulence and the strength of the mean temperature gradient. The potential temperature profile clearly indicates a steep gradient with a laminar appearance at ~2700 m topped by enhanced fluctuations  Figure 19 focuses on C, in the height range 2550-3000 m. It shows C 2 T , ε, ε CT2 with N 2 and < N 2 > and L O profiles and the corresponding profile of potential temperature θ. The ε peak is wider than the C 2 T peak and is shifted upward with respect to the C 2 T peak. ε CT2 and C 2 T obviously coincide when < N 2 > is used and ε CT2 roughly agrees with ε. However, when using N 2 , the agreement is much better. This particular case illustrates the cause of the statistical improvement when estimating ε CT2 from N 2 with respect to < N 2 > (Figure 8). The maximum of ε (or ε CT2 ) above the maximum of N 2 and a maximum of C 2 T between the two is remarkably similar to observations reported by Gossard [41] (Figure 4.5, p. 519), at the bottom of temperature inversions capping the CBL. Thus, mechanisms similar to those observed below capping inversions may occur in the free atmosphere. margins, respectively. In practice, the small offset (a few tens of meters or less) would not be enough for distinguishing the difference in altitudes between the VHF and UHF echoes, giving the impression that the same structures are detected. The existence of a small vertical offset between the structures responsible for VHF and UHF radar echoes is an appealing explanation of the paradox, but this explanation requires more thorough investigations, including radar data, for evaluating its relevance and significance. The vertical offsets between N 2 , C 2 T , ε and L O peaks can be interpreted as follows. On the one hand, peak of ε (maximum of Kolmogorov turbulence) is found to be far enough from the center of the stable inversion because the maximum of static stability N 2 inhibits vertical motion (and thus isotropic turbulence). The largest isotropic billows (maximum of L O ) can be observed there. They tend to be smaller when "approaching" the stable inversion, i.e., the Kolmogorov domain tends to be suppressed. On the other hand, large C 2 T requires both enhanced N 2 and ε. When ε (N 2 ) is maximum, N 2 (ε) is reduced so that C 2 T is also reduced. It follows that C 2 T is maximum between ε and N 2 peaks. The level of C 2 T is thus a compromise between the intensity of the dynamic turbulence and the strength of the mean temperature gradient. The potential temperature profile clearly indicates a steep gradient with a laminar appearance at~2700 m topped by enhanced fluctuations around~2730 m (Figure 19b).
These observations do not fit into the scheme of a turbulent sheet associated with enhanced C 2 T [16]. A temperature gradient sheet was defined as a temperature inversion of 3-20 m in depth associated with a temperature increase of 0.2-0.8 K corresponding to a gradient of 30-100 K/km [7]. The temperature inversion at~2700 m (Figure 19b) falls within this definition because the temperature increase is~0.7 K over~14 m (gradient~50 K/km). According to Dalaudier et al. [7], only 29% of the observed temperature gradient sheets in the troposphere and stratosphere were associated with overturning, and thus, potentially enhanced C 2 T . Temperature fluctuations were observed just above or below some sheets but the authors did not explicitly convey turbulence because velocity data were not available at that time. Using identical observation techniques, similar statistics were obtained from data described by Gavrilov et al. [3] (results not published). Our observations are consistent with this type of sheets. They confirm that ongoing mixing (enhanced ε and C 2 T ) is likely not fortuitous at the margins of some temperature sheets and that these sheets may be the result of, or may be reinforced by, turbulent mixing occurring at one of their margins.
This finding can contribute to a longstanding debate about the mechanisms of ground-based clear air radar backscattering. For ST radars operated in the Very High Frequency (VHF) band, echo power is often enhanced in the vertical direction, mainly above the tropopause, but also in the troposphere on some occasions [9]. Partial reflection from thin and horizontally stratified temperature and humidity gradient sheets is often proposed as the primary cause of this phenomenon [9,46]. For radars operated in the UHF band or higher frequencies, the clear air contribution to radar echo power is generally attributed to Bragg backscattering from isotropic turbulence [41]. In stratified conditions, UHF radar echoes also coincide with vertical gradients of temperature and humidity [41]. As a result, time-height cross-sections of VHF and UHF radar echo power do not qualitatively differ at similar range resolutions and both reveal the location of temperature and humidity gradients, when echoes are enhanced. Therefore, a given temperature gradient can be assumed to be either a laminar entity for explaining enhanced VHF radar echoes or a fully turbulent entity for explaining isotropic UHF echoes. Conditions such as those reported in Figure 19 may overcome this paradox. Due to the vertical offset between the temperature sheet and the maximum of C 2 T , the VHF and UHF radar echoes would mainly result from partial reflection from the temperature (and humidity) gradient and scattering from isotropic turbulence at one of its margins, respectively. In practice, the small offset (a few tens of meters or less) would not be enough for distinguishing the difference in altitudes between the VHF and UHF echoes, giving the impression that the same structures are detected. The existence of a small vertical offset between the structures responsible for VHF and UHF radar echoes is an appealing explanation of the paradox, but this explanation requires more thorough investigations, including radar data, for evaluating its relevance and significance.

Summary and Conclusions
In the present study, turbulence parameters in the lower troposphere have been estimated using data from high-resolution Pitot velocity and cold wire temperature sensors on board DataHawk UAVs, collected during ShUREX 2016 and 2017 campaigns at Shigaraki MU observatory. The data processing methods used for these estimations have been described in detail. These methods were designed to retrieve reliable estimates of TKE dissipation rate ε and temperature structure function parameter C 2 T from U and T measurements in the presence of measurement artifacts and instrumentation noise. These methods were based on the identification of the inertial −5/3 subranges in 1-D U and T spectra and on the estimation of their levels. Because a slope of −5/3 in (nearly) horizontal spectra can also be a signature of 2-D stratified turbulence for low Froude numbers or high gradient Richardson numbers [13,32], we considered this possibility by comparing the results with properties expected for Kolmogorov turbulence.
(1) The relationship between ε and C 2 T given by Equation (1), valid for Kolmogorov turbulence generated by shear flow instabilities in stratified conditions, was verified after excluding spurious data, confirming the result described by Kantha and Luce [21]. Equation (1) does not work for turbulence in convective boundary layers (e.g., Figures 9 and 10), and leads to puzzling results in saturated conditions with very good agreements or large discrepancies according to altitude ranges ( Figure 11). The traditional scheme of local diffusion (K-theory) does not work for convective boundary layers [47] and it can be expected that local stability (N 2 ) does not contribute in the relationship between ε and C 2 T . Further investigations are needed for these conditions. The best agreement was obtained with a flux Richardson number R f ∼ 0.13 − 0.15 and a mixing coefficient χ m~0 .16-0.17. These values are very close to the canonical values R f = 0.17 and χ m = 0.2 found for oceanic turbulence [25] and expected from theory [23]. The relatively large dispersion of the distribution of R f and χ m around these values can be partly due to uncertainties on N 2 (Figure 8) but similar dispersion of R f was reported for oceanic turbulence [22].
(2) Large values of ε are associated with small values of Ozmidov frequency, in agreement with the apparent extent of the −5/3 domain down to low frequencies ( Figure 2). Small values of ε, often rejected according to our criteria, were indeed found to be inconsistent with large values of Ozmidov frequency ( Figure 3).
(3) The distribution of turbulent eddy diffusivity K θ values based on ε estimates (Equation 3) fits very well the ranges of values reported from high-resolution radiosonde measurements at stratospheric heights ( Figure 12a) [29,44]. This result seems also to indicate that typical ε values at stratospheric heights do not statistically differ from those observed at tropospheric heights.
(4) The range of values of the buoyancy Reynolds number Re b (95% of the values exceed Re b = 200) is found to be consistent with the existence of isotropic turbulence according to the criteria for oceanic turbulence [45].
(5) Based on concurrent measurements of gradient Richardson numbers Ri from balloon data for a few flights, maxima and minima (or even undefined values) of ε are associated with minima and maxima of Ri, respectively, consistent with the classical scheme of turbulence generated by dynamic shear instabilities, as reported from similar UAV data [17].
All these observations are thus consistent with the Kolmogorov interpretation of the selected −5/3 subranges confirming the relevance of using Equation (1). Finally, we focused on the relationship between turbulence and sheets from a case-study. Turbulence is sometimes observed in close relationship with temperature inversions (consistent with earlier observations [41]). But, on some occasions, peaks of ε, C 2 T and N 2 show vertical offsets, the maximum of C 2 T being confined between the maxima of ε and N 2 (Figures 17 and 18). This characteristic is consistent with the enhancement of inertial temperature turbulence between a region of both enhanced stirring and reduced static stability and a region of both reduced stirring and enhanced static stability associated with a temperature gradient sheet (Figure 19). This feature was suggested from measurements of temperature only [7] and is confirmed here from additional ε measurements. Therefore, some of the temperature gradient sheets can be the result of, or can be reinforced by, turbulent mixing occurring at one of their margins. We speculated the possible impact of this sheet and layer structure on the interpretation of clear air radar echoes by suggesting that VHF radars would be sensitive to mainly the temperature (and humidity) sheet, when echoes are enhanced in the vertical direction and UHF radars to isotropic turbulence at one of the sheet margins. More investigations, including radar observations, will be performed for studying the relevance of this interpretation.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
Estimating the static stability (squared Brunt-Vaïsälä angular frequency) N 2 from time series of temperature data collected by a sensor enabling vertical profiling assumes that the vertical air velocity W air , in absolute terms, is small compared to vertical velocity of the sensor. If this is not the case, the estimates of N 2 are biased because of the temperature variations due to adiabatic compression/expansion of the air parcels during the measurements. Here, we do not consider additional changes due to horizontal advection and diabatic effects such as radiative heating or latent heat release.
For a temperature sensor onboard UAVs, the sensor moves up or down at the vertical velocity W UAV , controlled by the operator and is independent of W air . For a sounding meteorological balloon, the vertical velocity is W Balloon + W air where W Balloon is the ascent rate of the balloon in still air. In practice, N 2 is estimated from time series of duration ∆t. Within this period, the UAV intersected a layer of depth: where z A is the altitude of the bottom of the layer sampled by the UAV at time t and z B is the altitude of the top of the layer sampled by the UAV at time t + ∆t. The top of the layer was located at z B at time t. It moved up or down at z B due to W air . We have: The temperature T B of the air parcel at z B due to adiabatic expansion/compression is: where T B is the temperature of the air parcel when at the altitude z B and Γ a is the adiabatic lapse rate for dry air. Therefore, during ∆t, the sensor measures a temperature variation ∆T = T B − T A . A bulk estimate of N 2 from UAV measurements is: However, the effective value is: i.e., the value for W air = 0. Using (A1), (A2) and (A3), we get: Neglecting the term Γ a 2 W W UAV ∆z: Expression (A5) shows that a fair estimate of N 2 is obtained only if |W| << |W UAV |. For, W W UAV = 1, N 2 UAV = 0 whatever N 2 may be because the UAV is moving vertically with the air parcel that warms or cools at the adiabatic lapse rate. A more detailed analysis is beyond the scope of the paper but the interested reader can easily make quantitative estimates from (A5).