Turbulence kinetic energy dissipation rate: assessment of radar models from comparisons between 1.3 GHz wind profiler radar (WPR) and DataHawk UAV measurements

The WPR-LQ-7 is a UHF (1.3575 GHz) wind profiler radar used for routine measurements of the lower troposphere at Shigaraki MU Observatory (34.85^∘ N, 136.10^∘ E; Japan) at a vertical resolution of 100 m and a time resolution of 10 min. Following studies carried out with the 46.5 MHz middle and upper atmosphere (MU) radar (Luce et al., 2018), we tested models used to estimate the rate of turbulence kinetic energy (TKE) dissipation ε from the Doppler spectral width in the altitude range ∼ 0.7 to 4.0 km above sea level (a.s.l.). For this purpose, we compared LQ-7-derived ε using processed data available online (http://www.rish.kyoto-u.ac.jp/radar-group/blr/shigaraki/data/, last access: 24 July 2023) with direct estimates of ε (ε_U) from DataHawk UAVs. The statistical results reveal the same trends as reported by Luce et al. (2018) with the MU radar, namely (1) the simple formulation based on dimensional analysis ${\mathit{\epsilon }}_{L_{\mathrm{out}}}={\mathit{\sigma }}^{\mathrm{3}}/L_{\mathrm{out}}$, with L_out∼70 m, provides the best statistical agreement with ε_U. (2) The model ε_N predicting a σ²N law (N is Brunt–Vaïsälä frequency) for stably stratified conditions tends to overestimate for ${\mathit{\epsilon }}_{\mathrm{U}}\lesssim \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³ and to underestimate for ${\mathit{\epsilon }}_{\mathrm{U}}\gtrsim \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³. We also tested a model ε_S predicting a σ²S law (S is the vertical shear of horizontal wind) supposed to be valid for low Richardson numbers ($\mathit{Ri}=N^{\mathrm{2}}/S^{\mathrm{2}}$). From the case study of a turbulent layer produced by a Kelvin–Helmholtz (K–H) instability, we found that ε_S and ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ are both very consistent with ε_U, while ε_N underestimates ε_U in the core of the turbulent layer where N is minimum. We also applied the Thorpe method from data collected from a nearly simultaneous radiosonde and tested an alternative interpretation of the Thorpe length in terms of the Corrsin length scale defined for weakly stratified turbulence. A statistical analysis showed that ε_S also provides better statistical agreement with ε_U and is much less biased than ε_N. Combining estimates of N and shear from DataHawk and radar data, respectively, a rough estimate of the Richardson number at a vertical resolution of 100 m (Ri₁₀₀) was obtained. We performed a statistical analysis on the Ri dependence of the models. The main outcome is that ε_S compares well with ε_U for low Ri₁₀₀ (Ri₁₀₀≲1), while ε_N fails. ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ varies as ε_S with Ri₁₀₀, so that ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ remains the best (and simplest) model in the absence of information on Ri. Also, σ appears to vary as ${\mathit{Ri}}_{\mathrm{100}}^{-\mathrm{1}/\mathrm{2}}$ when Ri₁₀₀≳0.4 and shows a degree of dependence on S₁₀₀ (vertical shear at a vertical resolution of 100 m) otherwise.

1 Introduction

The dissipation rate ε (m² s⁻³ or W kg⁻¹) of turbulence kinetic energy (TKE) is an important variable for assessing the rate of change of TKE with time. This variable appears in a simplified expression of the ensemble-mean TKE budget equation (see, e.g., Stull, 1988, for a complete derivation and for its conditions of validity):

$$\begin{array}{}\text{(1)}& \partial \text{TKE}/\partial t=P-B-\mathit{\epsilon },\end{array}$$

where P is the shear production term and B the buoyancy flux term. Note that advection terms on the left-hand side have been omitted. Under steady-state conditions, the left-hand side is zero, and there exists a balance between the rates of shear production, buoyancy production/dissipation and dissipation of TKE. In practice, ε can be estimated in the atmosphere from VHF stratosphere–troposphere radar and UHF wind profiler measurements of Doppler spectral width, hereafter, noted 2σ_obs (m s⁻¹) (e.g., Hocking, 1983, 1985, 1986, 1999; Fukao et al., 1994; Cohn, 1995, Kurosaki et al., 1996; Bertin et al., 1997; Delage et al., 1997; Naström and Eaton, 1997; Dole et al., 2001; Jacoby-Kaoly et al., 2002; Satheesan and Krishna Murthy, 2002; Naström and Eaton, 2005; Wilson et al., 2005; Kalapureddy et al., 2007; Singh et al., 2008; Dehghan and Hocking, 2011; Kantha and Hocking, 2011; Dehghan et al., 2014; Wilson et al., 2014; Hocking et al., 2016; Li et al., 2016; Kohma et al., 2019; Jaiswal et al., 2020; Chen et al., 2022).

Several models have been proposed to relate 2σ_obs to ε. Some studies have accepted the validity of these models in order to perform statistical analyses of the turbulence characteristics in the troposphere–stratosphere (e.g., Fukao et al., 1994; Kurosaki et al., 1996; Naström and Eaton, 1997; Kalapureddy et al., 2007; Fukao et al., 2011; Chen et al., 2022). Other studies have tested the consistency between models based on spectral width measurement and their consistency with other radar models based on echo power or radial wind velocity variance measurement (e.g., Cohn, 1995; Bertin et al., 1997; Delage et al., 1997; Satheesan and Krishna Murty, 2002; Singh et al., 2008). Yet others assessed the radar estimates from cross-comparisons with indirect estimates based on the Thorpe sorting method applied to potential temperature profiles measured by standard radiosondes (e.g., Clayson and Kantha, 2008; Kantha, 2010; Kantha and Hocking, 2011; Wilson et al., 2014; Li et al., 2016; Kohma et al., 2019; Jaiswal et al., 2020). However, attempts at validations from direct in situ estimates of ε from velocity fluctuation measurements remain very rare. McCaffrey et al. (2017) compared ε estimates derived from a UHF wind profiler and those obtained from sonic anemometer energy spectra made at a 300 m altitude. Shaw and LeMone (2003) and Jacoby-Koaly et al. (2002) evaluated the performance of UHF wind profilers from in situ ε aircraft and/or tower measurements mainly in the convective boundary layer. Dehghan et al. (2014) made ε comparisons between aircraft and the VHF (40.68 MHz) Harrow radar with mixed results.

In addition to being rare, the above-mentioned studies did not aim to test the same radar models. Luce et al. (2018), hereafter denoted L18, assessed different models from comparisons with direct estimates of ε from air speed fluctuation measurements made from highly sensitive Pitot sensors aboard DataHawk UAVs and the VHF 46.5 MHz middle and upper atmosphere (MU) radar observations in the lower troposphere. One of the objectives of the present work is to show that the conclusions obtained from comparisons with the MU radar are also quantitatively valid for the WPR-LQ-7 (Imai et al., 2007), a UHF wind profiler routinely used at the Shigaraki MU Observatory. We also introduce and test another model expected to be valid for weakly stratified or strongly sheared conditions, i.e., low Richardson (Ri) numbers (Hunt et al., 1988; Basu et al., 2021; Basu and Holtslag, 2021) for which the static stability effects can be ignored. Ri is defined as $N^{\mathrm{2}}/S^{\mathrm{2}}$, where $N^{\mathrm{2}}=(g/\mathit{\theta })\mathrm{d}\mathit{\theta }/\mathrm{d}z$ is the squared Brunt–Väisälä frequency (s⁻²), g≈9.81 m s⁻² is the gravitational acceleration, θ (K) is the potential temperature and S (s⁻¹) is the vertical shear of the horizontal wind vector.

Section 2 introduces the expressions for ε used in the present paper with a focus on the newly introduced model in radar studies. Section 3 briefly describes the WPR-LQ-7 and the methods used for the comparisons. Section 4 describes the results for two turbulent layers, one of which was clearly produced by a Kelvin–Helmholtz (K–H; shear flow) instability because of the observation of S-shaped structures specific to this instability in both time–height MU and WPR-LQ-7 echo power cross-sections. The results of comparisons of ε values obtained from the different models applied to the two radars, DataHawk measurements, and a simultaneous radiosonde using the Thorpe sorting method of vertical potential temperature profiles are described for the two layers. Section 5 shows statistics on the consistency between the estimates of ε from the different models and the DataHawks and describes the dependence of the models on Ri. Finally, conclusions are given in Sect. 6.

2 The radar models of ε

2.1 The models tested by L18

The different models and their conditions of application have already been described by L18. Here, the expressions are simply reintroduced. Assuming a vertically pointing radar beam, the first expression is

$$\begin{array}{}\text{(2)}& {\mathit{\epsilon }}_{L_{\mathrm{out}}}={\mathit{\sigma }}^{\mathrm{3}}/L_{\mathrm{out}},\end{array}$$

where σ² is an estimate of the variance 〈w^′2〉 of the vertical wind fluctuations produced by turbulence. L_out has the dimension of a length scale and represents the scale of energy containing turbulent eddies if σ² is an unbiased estimate of 〈w^′2〉. In practice, σ² is obtained after removing the non-turbulent contributions from ${\mathit{\sigma }}_{\mathrm{obs}}^{\mathrm{2}}$ (e.g., Hocking, 1983; Nastrom, 1997; Hocking et al., 2016, and references therein). The practical method used in the present work is described in the Appendix of L18. The dissipation rate is expected to vary as σ³ if σ and L_out are independent or when the typical scale of turbulent eddies exceeds the dimensions of the radar volume, so that L_out would mainly be a function of these dimensions (e.g., Frisch and Clifford, 1974; Labitt, 1979; Doviak and Zrnić, 1984; White et al., 1999).

The second expression (e.g., Hocking, 1983, 1999; Hocking et al., 2016) is

$$\begin{array}{}\text{(3)}& {\mathit{\epsilon }}_N=C_N{\mathit{\sigma }}^{\mathrm{2}}N,\end{array}$$

where C_N is a constant. This expression is expected to be valid for turbulence in a stable stratification (N²>0) whose outer scale is defined by the buoyancy length scale expressed as $L_{\mathrm{B}}=\sqrt{\langle w^{\prime \mathrm{2}}\rangle }/N=\mathit{\sigma }/N$. Equation (3) is thus equivalent to ${\mathit{\epsilon }}_N=C_N{\mathit{\sigma }}^{\mathrm{3}}/L_{\mathrm{B}}$. In a pioneering contribution, Hocking (1983) first derived Eq. (3) from the integration of the transverse 1D spectrum of vertical velocities over the inertial and buoyancy subranges to relate ε to 〈w^′2〉. In its original derivation, the author assumed roughly equal contributions to 〈w^′2〉 from the inertial and buoyancy subranges. More recently, Hocking et al. (2016) proposed a more general expression by introducing a variable factor F, where F is the ratio of the buoyancy contribution to the inertial subrange contribution. This factor can vary from 0.5 to 1. It affects the value of the constant C_N, and Hocking et al. (2016) recommend that a value of (0.5±0.25) be used, which takes into account the variability of F, which is difficult to determine in practice.

The results of comparisons using Eqs. (2) and (3) reported by L18 with DataHawk-derived ε showed that Eq. (2) provides the best overall statistical comparisons for L_out∼ 50–70 m. On the other hand, the analysis of the comparison results with ε_N suggested that N is not a key parameter since the quality of the comparisons appeared to be independent of N.

2.2 The model for strongly sheared or/and weakly stratified flows

While it seems that the conditions of strong shear and weak stratification have not received much attention in the radar community, several studies showed that ε can be written as (Hunt et al., 1988; Schumann and Gerz, 1995)

$$\begin{array}{}\text{(4)}& {\mathit{\epsilon }}_{\mathrm{S}}=C_{\mathrm{S}}{\mathit{\sigma }}^{\mathrm{2}}S.\end{array}$$

Note that we do not make the distinction between σ² and 〈w^′2〉 for simplicity. Equation (4) is equivalent to ${\mathit{\epsilon }}_{\mathrm{S}}=C_{\mathrm{S}}{\mathit{\sigma }}^{\mathrm{3}}/L_{\mathrm{H}}$, where $L_{\mathrm{H}}=\mathit{\sigma }/S$ is the Hunt length scale. Equation (4) can be interpreted as the fact that turbulent eddies are first stretched by shear before being affected by stratification in strongly sheared or weakly stratified flows. This concept was discussed by Hocking and Hamza (1997), but they did not mention the Hunt length scale and did not go further into it. Hunt et al. (1988) suggested that Eq. (4) can be valid up to Ri∼0.5. Schumann and Gerz (1995) even proposed up to Ri∼1 from large eddy simulations. Hunt et al. (1988) proposed C_S≈0.45 for neutral stationary boundary layers. Kaltenbach et al. (1994) found $\mathrm{0.54}<C_{\mathrm{S}}<\mathrm{0.62}$ from large eddy simulations. By using Eq. (1) for a homogeneous shear layer in steady state, i.e., $\mathit{\epsilon }=P-B$, and similarity theory, Basu and Holtslag (2021) re-evaluated the constant C_S and provided a generalization of Eq. (4):

$$\begin{array}{}\text{(5a)}& \begin{array}{rl}{\mathit{\epsilon }}_{\mathrm{S}}^{\prime }& =C_{\mathrm{S}}{\left({\displaystyle \frac{\mathrm{1}-R_{\mathrm{f}}}{\mathit{Ri}}}\right)}^{\mathrm{1}/\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}N=C_{\mathrm{S}}{\left(\mathrm{1}-R_{\mathrm{f}}\right)}^{\mathrm{1}/\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}S\\ & =C_{\mathrm{S}}^{\prime }\left(R_{\mathrm{f}}\right){\mathit{\sigma }}^{\mathrm{2}}S,\end{array}\end{array}$$

with C_S=0.63. R_f is the flux Richardson number. It is related to the turbulent Prandtl number Pr by $R_{\mathrm{f}}=\mathit{Ri}/\mathit{Pr}$. Basu et al. (2021) found from direct numerical simulation (DNS) that Eq. (4) with C_S∼0.60 is valid up to Ri∼0.2 at least. For $\mathrm{0}<\mathit{Ri}\lesssim \mathrm{0.25}$, $C_{\mathrm{S}}^{\prime }\left(R_{\mathrm{f}}\right)$ decreases from 0.63 to 0.60, using the analytical expression (Eq. 22) of Basu and Holtslag (2021) for Pr(Ri). For Ri→0, we have R_f→0, then ${\mathit{\epsilon }}_{\mathrm{S}}^{\prime }\to {\mathit{\epsilon }}_{\mathrm{S}}=\mathrm{0.63}{\mathit{\sigma }}^{\mathrm{2}}S$. For Ri→1, R_f∼0.25, ${\mathit{\epsilon }}_{\mathrm{S}}^{\prime }\to \mathrm{0.5}{\mathit{\sigma }}^{\mathrm{2}}N$. Therefore, Eq. (3) would be quantitatively equivalent to Eq. (5a) for Ri of the order of 1 despite the fact that the two approaches cannot be reconciled because, in essence, there is no contribution from an anisotropic buoyancy subrange in Eq. (5a). Equation (4) removes an inconsistency in Eq. (3), since it wrongly indicates that ε→0 when N→0 for a given σ². If S=0, i.e., if the source of the instability that generates turbulence is removed, then ε=0, which makes more sense.

As discussed by Basu and Holtslag (2021, Sect. 6.2) and Basu and Holtslag (2022, their Appendix 1), the derivation of Eq. (5) does not consider the fact that the steady-state condition (also called “full equilibrium” (FE), Baumert and Peters, 2000) can only be reached for a single value of Richardson number Ri_s, at least for large Reynolds numbers and large shear parameters ST_L, where T_L is the inertial timescale defined as $\text{TKE}/\mathit{\epsilon }$ (see, e.g., Mater and Venayagamoorthy, 2014). For Ri<Ri_s, TKE increases at subcritical Ri, and for Ri>Ri_s, TKE decreases (turbulence decays) at supercritical Ri (e.g Baumert and Peters, 2000). From large eddy simulation (LES) and DNS data, Schumann (1994) and Gerz et al. (1989) reported Ri_S≈0.13 for air, consistent with the value that can be deduced from Fig. 1 of Mater and Venayagamoorthy (2014). Schumann (1994) re-wrote the TKE budget Eq. (1) as $\partial \text{TKE}/\partial t=(G-\mathrm{1})(\mathit{\epsilon }+B)$ where $G=P/(\mathit{\epsilon }+B)$ is called the growth factor. G=1 for FE conditions. By assuming, for simplicity, that G only depends on Ri, Schumann proposed the empirical expression $G\left(\mathit{Ri}\right)=G_{\mathrm{0}}^{(\mathrm{1}-\mathit{Ri}/{\mathit{Ri}}_{\mathrm{s}})}$ with $G_{\mathrm{0}}=\mathrm{1.47}\pm \mathrm{0.13}$ based on wind-tunnel data analysis. By using the same procedure as Basu and Holtslag (2021) from their Eqs. (10) to (12) but starting with $\mathit{\epsilon }=P/G-B$, we get

$$\begin{array}{}\text{(5b)}& {\mathit{\epsilon }}_{\mathrm{S}}^{\prime \prime }={\displaystyle \frac{\mathrm{0.63}}{G^{\mathrm{1}/\mathrm{2}}}}{\left(\mathrm{1}-G{R}_{\mathrm{f}}\right)}^{\mathrm{1}/\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}S.\end{array}$$

Equation (5b) can also be directly obtained by using Eq. (46a) in Eq. (10b) of Basu and Holtslag (2021). For $\mathrm{0}<\mathit{Ri}\lesssim \mathrm{0.25}$, $C_{\mathrm{S}}^{\prime \prime }=C_{\mathrm{S}}/G^{\mathrm{1}/\mathrm{2}}{\left(\mathrm{1}-G{R}_{\mathrm{f}}\right)}^{\mathrm{1}/\mathrm{2}}$ increases from 0.52 to 0.70, i.e., ∼0.60 on average, for Ri_s≈0.13 and G₀=1.47. The Ri dependence of $C_{\mathrm{S}}^{\prime \prime }$ is thus only a small source of dispersion for low Ri values when comparing with other estimates.

From the results of Baumert and Peters (2000) using a “structural equilibrium” approach (i.e., stationarity of ratios of turbulence characteristics) and based on laboratory and LES data (their Fig. 4), we can establish ε_S=0.15σ²S valid for Ri≲0.25. This expression is obtained by combining $L_{\mathrm{H}}/L_{\mathrm{B}}={\mathit{Ri}}^{\mathrm{1}/\mathrm{2}}$, $L_{\mathrm{E}}/L_{\mathrm{O}}=\mathrm{4.2}{\mathit{Ri}}^{\mathrm{3}/\mathrm{4}}$ and $L_{\mathrm{E}}/L_{\mathrm{B}}=\mathrm{1.61}{\mathit{Ri}}^{\mathrm{1}/\mathrm{2}}$, where $L_{\mathrm{E}}=\sqrt{\langle {{\mathit{\theta }}^{\prime }}^{\mathrm{2}}\rangle }/(\mathrm{d}\mathit{\theta }/\mathrm{d}z)$ and $L_{\mathrm{O}}=\sqrt{\mathit{\epsilon }/N^{\mathrm{3}}}$ are the Ellison and Ozmidov length scales, respectively. The constant differs very significantly (by a factor of 3 to 4 less) from the aforementioned estimates. If we use $L_{\mathrm{E}}/L_{\mathrm{O}}=\mathrm{2.4}{\mathit{Ri}}^{\mathrm{3}/\mathrm{4}}$ as proposed by Schumann (1994) for Ri≤0.25, we get ε_S=0.44σ²S with the same $L_{\mathrm{E}}/L_{\mathrm{B}}$ ratio. These expressions are more subject to experimental uncertainties and are thus not considered in this paper. In Appendix A, we propose an alternative derivation of Eq. (4), suggesting $\mathrm{0.45}\le C_{\mathrm{S}}\le \mathrm{0.82}$. We retain the value of C_S=0.63 for the comparisons between the models.

Following the spectral approach proposed by Weinstock (1981), Eq. (4) with C_S equal to C_N≈0.5 can also be obtained from the integration of the 1D Kolmogorov ($-\mathrm{5}/\mathrm{3}$ slope) scalar kinetic energy spectrum over a spherical shell of radius k_H instead of k_B, where k_H(k_B) is the wavenumber corresponding to L_H(L_B). For the context of radar measurements (e.g., Hocking et al., 2016), Eq. (4) can also be obtained from the integration of the 1D transverse vertical velocity spectrum with a $-\mathrm{5}/\mathrm{3}$ slope for large (horizontal) wavenumbers (k>k_H) and a zero slope for (k<k_H), both mathematical developments being equivalent.

Finally, ε_S has the advantage that it can be evaluated entirely from the radar data, since the wind shear S can be estimated at the range and time resolutions of the radar, unlike ε_N which requires N² to be obtained from in situ or radio-acoustic sounding system measurements.

3 The WPR-LQ-7 and methods of comparisons with UAV data

3.1 The WPR-LQ-7

The WPR-LQ-7 is a 1.3575 GHz Doppler radar. It has a phased array antenna composed of seven Luneberg lenses of 800 mm diameter. Its peak output power is 2.8 kW. It can be steered into five directions sequentially (i.e., after fast Fourier transform (FFT) operations), vertical and 14.2^∘ off zenith toward north, east, south and west. The main radar parameters of the WPR-LQ-7 installed at Shigaraki MU Observatory since 2006 are given in Table 1.

Table 1 WPR-LQ-7 parameters in routine observation mode.

The acquisition time for one profile composed of 80 altitudes from 300 m a.g.l. (∼ 684 m a.s.l.) every 100 m in each direction is 59 s after 18 incoherent integrations but for a total of 11.8 s of observations for each direction (due to the intertwining between the directions). The time series are processed by automatic algorithms to remove outliers (e.g., bats, birds, airplanes) and ground clutter as far as possible. Low signals near and below the detection thresholds are removed, and profiles of atmospheric parameters (echo power, radial winds, half-power spectral width, horizontal and vertical winds) averaged over 10 min are made available for routine monitoring (http://www.rish.kyoto-u.ac.jp/radar-group/blr/shigaraki/data/, last access: 24 July 2023). Because of the high data quality control, the 10 min averaged data are used to retrieve ε with the objective to assess the routine data for further analysis. The 59 s resolution data and those collected by the MU radar at a time resolution (sampling) of 24.57 s (∼ 12.3 s) were used to help identify atmospheric structures from height–time signal-to-noise ratio (SNR) or echo power cross-sections, such as convective cells or Kelvin–Helmholtz (K–H) billows. Table 2 shows the acquisition time, the range and transverse resolutions of the WPR-LQ-7 for the altitude range of comparisons, and the range and resolutions of the MU radar for the data used in the present work.

Table 2 Time, range and transverse resolutions of the MU radar and WPR-LQ-7 for the dataset used in the present work. NINCOH refers to the number of incoherent integrations. The range resolution is $\mathrm{\Delta }\mathrm{r}=\mathrm{1}/\mathrm{2}\mathrm{c}\mathit{\tau }$, where c is the light speed and τ is the pulse duration, and the transverse resolution is 2θ₀z, where θ₀ is the half-power half width of the effective (two-way) radar beam and z is the altitude as defined in L18. The time series of MU radar signals are weighted by a Hanning window before FFT calculations.

3.2 The methods of comparison with DataHawk-derived ε

The DataHawk datasets were collected during two field campaigns, called the Shigaraki UAV Radar Experiment (ShUREX), in May–June 2016 and June 2017 at the Shigaraki MU observatory. The DataHawks were flying about 1 km away from the MU radar and the WPR-LQ-7. Kantha et al. (2017) described the instruments and configurations used during a previous ShUREX campaign in June 2015. The processing method used to retrieve ε from Pitot sensor data is not recalled here as it is described in detail by L18 for comparisons with MU radar data. The trajectories of the DataHawks being helicoidal upwards or downwards, pseudo-vertical profiles of ε at a vertical sampling of ∼ 5 m typically were obtained during the ascents and descents of the aircraft, from the ground up to a maximum altitude of ∼ 4.5 km. A total of 36 DataHawk datasets collected during the two campaigns provided 90 full or partial profiles used for the comparisons. Section 4 describes one of these flights with one full ascent (A1) and descent (D2) and one partial ascent (A2) and descent (D1). Three DataHawk flights collected during periods of precipitation contaminating the WPR-LQ-7 returns were rejected. The DataHawk-derived ε profiles were smoothed with a Gaussian window and resampled at the altitude of the radar gates to simulate the radar range resolution. The degraded DataHawk 100 m resolution profiles are hereafter denoted ε_U.

The WPR-LQ-7-derived ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$, ε_N and ε_S profiles were computed at a time resolution of 10 to 30 min, i.e., by averaging up to three consecutive profiles that best correspond to each period of DataHawk ascent or descent. ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ was calculated with L_out=70 m in accordance with the best agreement from comparisons with the MU radar (L18) and the statistics of L_out shown in Sect. 5. The profiles of N² at a vertical resolution of 100 m were estimated from pressure and temperature profiles collected by the DataHawks. These profiles were used to obtain ε_N (Eq. 3). ε_S profiles were computed from Eq. (4) every 10 min from wind shear estimated from radar data and then averaged up to 30 min.

3.3 Estimation of ε from the Thorpe method applied to radiosonde data

The Ozmidov length scale $L_{\mathrm{O}}=\sqrt{\mathit{\epsilon }/N^{\mathrm{3}}}$ is commonly assumed to be proportional to the Thorpe length $L_{\mathrm{T}}\triangleq \langle {\mathrm{d}}^{\prime \mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}$, where d^′ is the Thorpe displacement in the so-called Thorpe layer. Then, we have L_O=cL_T and

$$\begin{array}{}\text{(6)}& {\mathit{\epsilon }}_{\mathrm{T}}=c^{\mathrm{2}}{L}_{\mathrm{T}}^{\mathrm{2}}{N}^{\mathrm{3}}.\end{array}$$

The literature is very divided on the value of c to use $(\mathrm{0.25}<c<\mathrm{4})$ (see Kohma et al., 2019, for a review). Wijesekera and Dillon (1997) showed large temporal variations of $L_{\mathrm{T}}/L_{\mathrm{O}}$ from observations in the ocean. Large temporal variations were also reported from DNS depending on the stage and source of turbulence (e.g., Fritts et al., 2016). An intermediate value of c=1 is sometimes used by default (e.g., Kantha and Hocking, 2011). However, Mater et al. (2013) showed that c∼1 when the turbulent Froude number $\mathit{Fr}=\mathit{\epsilon }/\left(N\text{TKE}\right)$ is near unity (at the transition between shear- and buoyancy-dominated regimes). The basic N² for the Thorpe layers is generally estimated from the sorted potential temperature profile $\left(N_{\mathrm{s}}^{\mathrm{2}}\right)$ or from the rms value of the fluctuations defined as the difference between the measured and sorted profiles $\left(N_{\mathrm{rms}}^{\mathrm{2}}\right)$ (e.g., Smyth et al., 2001; Wilson et al., 2014).

Figure 1 (a) Time–height cross-section of WPR-LQ-7 signal-to-noise ratio (dB) at vertical incidence on 18 June 2017 from 13:30 to 17:30 LT. (b) The corresponding time–height cross-section of MU radar echo power (dB) in (high-resolution) range imaging mode at vertical incidence. “A1”, “D1”, “A2” and “D2” refer to the consecutive ascents and descents of the DataHawk UAV (DH35), emphasized by the red lines. The blue line shows the time–altitude of the radiosonde V6 launched at 14:51 LT from Shigaraki MU Observatory.

Another scale, called the Corrsin length scale, is defined as $L_{\mathrm{c}}=\sqrt{\mathit{\epsilon }/S^{\mathrm{3}}}$. It is the counterpart of the Ozmidov length scale for shear flows under neutral stratification conditions. Similarly, assuming $L_{\mathrm{c}}=c^{\prime }{L}_{\mathrm{T}}$, we can write

$$\begin{array}{}\text{(7)}& {\mathit{\epsilon }}_{T^{\prime }}=c^{\prime \mathrm{2}}{L}_{\mathrm{T}}^{\mathrm{2}}{S}^{\mathrm{3}}.\end{array}$$

Equation (7) is thus a possible alternative to Eq. (6), when the Corrsin length scale is smaller than the Ozmidov length scale. These equations are coherent with the results of Mater et al. (2013), who showed that L_T scales with ${\left(\text{TKE}\right)}^{\mathrm{1}/\mathrm{2}}/S$ in the shear-dominated regime and with ${\left(\text{TKE}\right)}^{\mathrm{1}/\mathrm{2}}/N$ in the buoyancy-dominated regime. Equation (7) can also be justified, and the parameter c^′ can be estimated as follows. The aforementioned ratios $L_{\mathrm{E}}/L_{\mathrm{O}}=\mathrm{4.2}{\mathit{Ri}}^{\mathrm{3}/\mathrm{4}}$ and $L_{\mathrm{E}}/L_{\mathrm{O}}=\mathrm{2.4}{\mathit{Ri}}^{\mathrm{3}/\mathrm{4}}$, found for weakly stratified flows (Ri≤0.25) by Baumert and Peters (2000) and Schumann (1994), respectively, may be representative of $L_{\mathrm{T}}/L_{\mathrm{O}}=\mathrm{1}/c\left(\mathit{Ri}\right)$ for flows free of gravity wave motions. Indeed, L_E=L_T is obtained if $\mathrm{d}\mathit{\theta }/\mathrm{d}z$ is a constant, which implies the absence of gravity waves. By introducing the expressions of $L_{\mathrm{T}}/L_{\mathrm{O}}=\mathrm{1}/c\left(\mathit{Ri}\right)$ into Eq. (6), we obtain Eq. (7), with c^′ (hereafter noted $c_{\mathrm{Sc}}^{\prime }$) equal to $\mathrm{1}/\mathrm{2.4}=\mathrm{0.41}$ or c^′ (hereafter noted $c_{\mathrm{BP}}^{\prime }$) equal to $\mathrm{1}/\mathrm{4.2}=\mathrm{0.24}$. Note that c^′ is a constant, while c depends on Ri. For a shear-dominated regime, from Fig. 3e, f, g and h of Mater and Venayagamoorthy (2014), we can typically deduce $\mathrm{0.25}<c^{\prime }<\mathrm{0.5}$ from DNS and $c^{\prime }\sim \mathrm{0.33}$ from experimental data, which is very consistent with the other values. In Appendix B, we show that $c^{\prime }=\mathrm{0.28}$ can be found from an alternative approach based on the inference of the turbulent Froude number from $L_{\mathrm{E}}/L_{\mathrm{O}}$ for weakly stratified flow condition (Garanaik and Venayagamoorthy, 2019).

Equations (6) and (7) are equivalent if Eq. (5) is written as

$$\begin{array}{}\text{(8)}& {\mathit{\epsilon }}_{\mathrm{T}}=c^{\prime \mathrm{2}}{\mathit{Ri}}^{-\mathrm{3}/\mathrm{2}}{L}_{\mathrm{T}}^{\mathrm{2}}{N}^{\mathrm{3}}.\end{array}$$

Equations (4) and (7) have in common that they are formally independent of N² when Ri≲0.25. If they were both to be confirmed by experimental analysis, they would constitute a coherent whole.

4 Two case studies

Figure 1a shows the time–height cross-section of WPR-LQ-7 SNR (dB) at vertical incidence and a time and range resolution of 59 s and 100 m, respectively, on 18 June 2017 from 13:30 to 17:30 LT and in the altitude range [0.685–7.0 km] a.s.l. (a.s.l. = a.g.l. + 0.385 km). Figure 1b shows the corresponding cross-section of MU radar echo power (dB) at vertical incidence and a time resolution of ∼12 s after doing range imaging with Capon processing (e.g., Luce et al., 2017) in the altitude range [1.275–7.0 km]. Radar echoes from a DataHawk, called “DH35” in reference to the flight numbering, are visible after ∼ 14:30 and before ∼ 15:40 LT on both images. They are the signatures of two ascents (“A1”, “A2”) and two descents (“D1”, “D2”) of a DataHawk. Four red segments emphasize them in Fig. 1b. Incidentally, radar echoes from another DataHawk (DH36) can be noted after 16:30 LT. A Vaisala RS92-SGP radiosonde, called “V6”, was launched at 14:51 LT from the observatory. Its time–height position is indicated by the blue line in Fig. 1b. It roughly coincides with A1.

An approximately 800 m deep enhanced echo power layer with S-shaped structures, signature of Kelvin–Helmholtz billows, is clearly visible in both images in the altitude range [3.0–4.0 km] until ∼ 16:00 LT at least. The layer is denoted by “KHI” on the images. DH35 crossed the layer four times during A1, D1, A2 and D2 between 15:00 and 15:30 LT. DH35 sampled the most obvious case of K–H instability during the entire two campaigns. Since a necessary condition for the development of K–H billows is $\mathit{Ri}<{\mathit{Ri}}_c=\mathrm{0.25}$ at the critical level, their observation suggests that it was fulfilled when sampled by the instruments. Another focus will be given to a turbulent layer between 2100 and 2500 m, sampled twice by DH35 during A1 and D2, even though it is not clearly visible in the radar echo power images (Fig. 1). For this layer, Ri is expected to be ≳1 according to various estimates and layer properties described in Sect. 4.2.

The analysis of these two cases is done to illustrate the differences in the ε behavior of the three different radar models applied to two radars, possibly at different Richardson numbers (Ri≲0.25 and ≳1), compared to the DataHawk-derived ε.

4.1 The K–H layer

4.1.1 Comparisons between DataHawk-derived ε, ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$, ε_N and ε_S

Figure 2a shows the four DataHawk-derived ε profiles during A1, D1, A2 and D2 (dotted black lines) and the profiles of ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$, ε_N and ε_S in the height range [2000–3900] m obtained from the WPR-LQ-7 data (solid red, blue and green lines, respectively) and MU radar data (dashed red, blue and green lines, respectively). Figure 3a, b and c show the same information for the three models but separately. For clarity, because they are virtually identical during D1, A2 and D2, the radar-derived ε profiles are shown for A1 (15:00–15:20 LT) only. The sources of errors and uncertainties on radar estimates are multiple (e.g., Dehghan and Hocking, 2011), and confidence intervals for each individual estimate are difficult or even impossible to establish. However, the consistency between the estimates from successive independent segments of data and concordant temporal evolution as shown in Fig. 4 and Table 3 gives some credence to the significance of the comparisons. Table 3a and the corresponding figure, Fig. 4, show ε, ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$, ε_N and ε_S values averaged over the depth of the K–H layer for A1, D1, A2 and D2. The DataHawk-derived ε values peak in the range of the K–H layer and vary little during the ascents and descents over ∼ 30 min: typically ∼ 2 mW kg⁻¹. During A1, the DataHawk-derived ε profile shows a narrower peak between 3200 and 3600 m. The DataHawk may have sampled a thinner region of the K–H layer (∼400 m), perhaps associated with the edge of a K–H billow. This could also be the case for V6 as the Thorpe analysis suggests a ∼300 m deep layer at the altitude of ∼3.3 km (and an additional thinner layer around the altitude of ∼3.6 km). If we exclude the difference in layer depth during A1, ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ and ε_S estimated from both radars coincide very well with DataHawk-derived ε both in shape and levels during A1, D1, A2 and D2, with very similar variations in time (Table 3a and Fig. 4), indicating that the two radar models are satisfactory and are equivalent in these circumstances. In contrast, the ε_N profiles exhibit the worst agreement with DataHawk-derived ε near the center of the K–H layer where they show a minimum (Figs. 2a and 3b, solid and dashed blue lines). This feature is similar to the one reported by L18 (their Fig. 12) for a turbulent layer generated by a convective instability at a mid-level cloud base. Table 3a and Fig. 4 confirm that 〈ε_N〉 values are lower than the other estimates by a factor of 2 to 3 approximately during A1, D1, A2 and D2 for both radars. This disagreement, occurring repeatedly on the two radars, confirms the inadequacy of the ε_N model for this layer.

Figure 2 (a) DataHawk-derived ε (m² s⁻³) profiles in the height range [2000–3900] m during A1, D1, A2 and D2 of DH35 on 18 June 2017 (dotted black), ε_S(LQ-7) profile (solid red), ${\mathit{\epsilon }}_{L_{\mathrm{out}}}\left(\text{LQ-7}\right)$ profile (solid green), ε_N(LQ-7) profile (solid blue), ε_S(MU) profile (dashed red), ${\mathit{\epsilon }}_{L_{\mathrm{out}}}\left(\text{MU}\right)$ profile (dashed green) and ε_S(MU) profile (dashed blue) derived from radar data between 15:00 and 15:20 LT. Grey dots and lines show ε_T (Eq. 6) with c=1, the depth and altitude of the Thorpe layers. (b) Richardson number profiles estimated from RS92-SGP Vaisala radiosonde V6 data at a vertical resolution of 20 m (black) and 100 m (red).

Figure 3 Same as Fig. 2a but with separate plots for each model: (a) ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$, (b) ε_N, (c) ε_S and (d) ε_T, respectively. Panel (d) shows the results in magenta for ε_T (Eq. 6) with c=1 (solid line), ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }$ (Eq. 7) with ${\mathrm{c}}^{\prime }=\mathrm{1}$ (dashed line), and ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }$ (Eq. 7) with ${\mathrm{c}}^{\prime }=\mathrm{0.41}$ and ${\mathrm{c}}^{\prime }=\mathrm{0.24}$ (dotted line).

Table 3 Mean values of TKE dissipation rates (mW kg⁻¹) according to the different models and instruments for the K–H layer (a) and for the turbulent layer (TL) between 2100 and 2500 m (b).

^a This low value is due to a MU-radar-derived wind shear about 2 times smaller than LQ-7- and balloon-derived wind shear. It is doubtful and will affect L_H and L_C in Table 4.
^b (sorted/rms).

Figure 4 Graphical representation of the averaged estimates of TKE dissipation rates (mW kg⁻¹) (1 mW kg${}^{-\mathrm{1}}={\mathrm{10}}^{-\mathrm{3}}$ m² s⁻³) shown in Table 3 for the K–H layer (3000–3600 m) sampled four times (top) and the layer between 2100 and 2500 m sampled twice during A1 (segment no. 1) and D2 (segment no. 4) (bottom). The horizontal black line at 0.5 mW kg⁻¹ separates the two cases, and each radar-derived value has been successively shifted to the right by 0.025 with respect to segment numbers for clarity.

4.1.2 Comparison between DataHawk-derived ε and ε_T

The altitude and depth of the turbulent layers identified by the Thorpe method from V6 and ε_T (Eq. 6) with c=1 are shown by the dots and the solid vertical segments, respectively, in grey in Fig. 2a and in magenta in Fig. 3d. In Fig. 3d, ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }$ values (Eq. 7) with $c^{\prime }=\mathrm{1}$, $c_{\mathrm{BP}}^{\prime }=\mathrm{0.24}$ and $c_{\mathrm{Sc}}^{\prime }=\mathrm{0.41}$ are also shown for the K–H layer at 3.33 km and the turbulent layer at 2.37 km discussed in Sect. 4.2. $N_{\mathrm{s}}^{\mathrm{2}}$ and $N_{\mathrm{rms}}^{\mathrm{2}}$ for the K–H layer are $\mathrm{7.7}\times {\mathrm{10}}^{-\mathrm{6}}$ and $\mathrm{7.1}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻², respectively, i.e., $\mathrm{7.4}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻² on average. Because L_T=130 m in the K–H layer, we obtain ε_T≈0.35 mW kg⁻¹, which is about 7 times lower than DataHawk-derived ε (2.4 mW kg⁻¹) (Table 3a). The hypothesis that V6 passed through the K–H layer in a region where ε was much lower is not consistent with the low variability (stationarity) of the dissipation rates estimated from DataHawk and radar data for more than 40 min (see Table 3a and Fig. 4). We therefore assume that ε_T must be ≈2.4 mW kg⁻¹. To achieve this condition with Eq. (6), we must have c=2.6.

On the other hand, estimating ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }$ (Eq. 7) requires the retrieval of S, but there is no prescribed method to compute the vertical shear of horizontal wind from balloon data in the Thorpe layers. Here, we estimated S from the difference of the wind vectors at the extremities of the Thorpe layer and from a linear interpolation of the zonal and meridional wind components in the Thorpe layer. We found S=0.013 and 0.010 s⁻¹, respectively, i.e., 0.0115 s⁻¹ on average, so that Ri≈0.055. This value is close to the mean value (〈Ri〉=0.09) obtained at a vertical resolution of 20 m (Fig. 2b). The relevance of ε_T≈0.35 mW kg⁻¹ obtained with c=1 from Eq. (6) can be tested from Eq. (8) with $c=c^{\prime }{\mathit{Ri}}^{-\mathrm{3}/\mathrm{4}}=\mathrm{1}$ using $c_{\mathrm{BP}}^{\prime }=\mathrm{0.24}$ and $c_{\mathrm{Sc}}^{\prime }=\mathrm{0.41}$. We get Ri=0.15 and Ri=0.30, respectively. These values are significantly larger than 0.055. For S=0.0115 s⁻¹ and $c^{\prime }=\mathrm{1}$, we get ${\mathit{\epsilon }}_{{\mathrm{T}}^{\prime }}\approx \mathrm{25.7}$ mW kg⁻¹, i.e., about 11 times larger than DataHawk-derived ε. We must have $c^{\prime }=\mathrm{0.31}$ to be consistent with the DataHawk-derived ε value. $c_{\mathrm{Sc}}^{\prime }=\mathrm{0.41}$ and $c_{\mathrm{BP}}^{\prime }=\mathrm{0.24}$ (and $c^{\prime }=\mathrm{0.28}$ found in Appendix B) reasonably meet the necessary correction, giving credence to the validity of Eq. (7). As a corollary, for Ri=0.055, we would get $c=c^{\prime }{\mathit{Ri}}^{-\mathrm{3}/\mathrm{4}}=\mathrm{2.7}$, i.e., the required value of c for Eq. (6) to be valid. The various estimates of ${\mathit{\epsilon }}_{{\mathrm{T}}^{\prime }}$ for the Thorpe layer are shown in Fig. 3d. Although based on a fragile hypothesis (ε from Thorpe analysis of the radiosonde data is equal to DataHawk-derived ε), Eq. (7) appears to be more adapted than the standard model (Eq. 6). It also has the major advantage over Eq. (6) that c^′ is a true constant, at least when Ri<0.25, although its value remains to be defined more precisely.

4.1.3 Comparison of turbulence scales estimated from radar data

Table 4a shows the Hunt, Corrsin, buoyancy and Ozmidov length scales for the K–H layer calculated from WPR-LQ-7-derived and MU-radar-derived ε, σ and S during A1, D1, A2 and D2. N² is computed from balloon data at the radar resolutions (100 m for the WPR-LQ-7 and 150 m for the MU radar). Once the scales are calculated, they are averaged over the altitude range 3000–3600 m of the K–H layer to compare them with the Thorpe length. All the radar-derived scales reveal the same behaviors between the segments A1, D1, A2, and D2 and do not show substantial differences between the radars, reinforcing the reliability of numerical results and their interpretations. Only the average values for A1, D1, A2 and D2 are discussed. We get 〈L_H〉=46 m, 〈L_C〉=38 m, 〈L_B〉=105 m and 〈L_O〉=131 m. 〈L_H〉 and 〈L_C〉 are substantially smaller than 〈L_B〉 or 〈L_O〉, indicating the latter should not be the scales to consider, as expected from the analysis of Sect. 4.1.1. Depending on the flight segment (A1, D1, A2, D2), $c=\langle L_{\mathrm{O}}\rangle /L_{\mathrm{T}}$ was found between 0.65 and 1.56 and ∼ 1 on average. In contrast, we get $\langle L_{\mathrm{C}}\rangle /L_{\mathrm{T}}=c^{\prime }\approx$ 0.23–0.37. It is smaller than $c_{\mathrm{Sc}}^{\prime }=\mathrm{0.41}$ but close to $c_{\mathrm{BP}}^{\prime }=\mathrm{0.24}$ and the value obtained in Appendix B ($c^{\prime }=\mathrm{0.28}$) and the needed value of 0.31. We obtain ${\left(\langle L_{\mathrm{H}}\rangle /\langle L_{\mathrm{B}}\rangle \right)}^{\mathrm{2}}=\langle \mathit{Ri}\rangle =\mathrm{0.19}$ and ${\left(c^{\prime }/c\right)}^{\mathrm{4}/\mathrm{3}}=\langle \mathit{Ri}\rangle =\mathrm{0.28}$. Both radar estimates are significantly larger than Ri estimated from balloon data with the Thorpe analysis but are close to 〈Ri〉=0.33 estimated from balloon data at the vertical resolution of 100 m (Fig. 2b). The quantitative disagreements result mainly from comparisons between estimates made with different techniques (radar and in situ) and resolutions. Nevertheless, it is interesting to note that these comparisons tend to corroborate the conclusions obtained from in situ measurements alone (Sect. 4.1.2).

Table 4 Mean values of Hunt, Corrsin, buoyancy and Ozmidov length scales for the K–H layer (a) and for the turbulent layer (TL) between 2100 and 2500 m (b).

^* Doubtful (see Table 3).

4.2 The turbulent layer between 2100 and 2500 m

Using the same methods as for the K–H layer, we obtain 〈Ri〉=0.75 (14.5) at a vertical resolution of 20 (100) m from V6 data (Fig. 2b), Ri≈2.0 from Thorpe analysis of V6 in the altitude range [2100–2500 m] (not shown) and 〈Ri〉=4.6 from N² calculated at a vertical resolution of 100 m from V6 data and S calculated from WPR-LQ-7 data during A1 and D2 (Table 4b). Therefore, the Richardson number strongly varies according to the method and data used, but all the estimates are consistent with a Ri value significantly larger than for the K–H layer and likely larger than 1 (Sect. 4.1). Therefore, the weakly stratified condition (Ri<0.25) for which the alternative Eqs. (4) and (7) are valid is likely not verified for this layer. DataHawk-derived ε is about 1 order of magnitude lower than for the K–H layer: ∼ 0.2 mW kg⁻¹ on average (Table 3b). The mean value is less reliable during D2. Many values of DataHawk-derived ε are missing because the algorithm did not detect a $-\mathrm{5}/\mathrm{3}$ subrange in the velocity spectra (see L18). Figure 2a, Table 3b and Fig. 4 show that ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ and ε_N and their mean values derived from both radars and ε_T are close to each other (within a factor of less than ∼ 2) and are very consistent with DataHawk-derived ε. All the radar and DataHawk estimates together show a temporal decrease by a factor of 3 to 4 in the 40 min between A1 and D2 (Fig. 4), giving credence that the agreements between the various estimates during A1 and D2 are not fortuitous. The temporal decrease of ε is consistent with a decaying turbulence when Ri>1. 〈ε_S〉 shows the largest discrepancies with DataHawk-derived ε, perhaps because the model is not valid for large Ri values. 〈L_H〉 and 〈L_C〉 (136 m and 202 m, respectively) exceed 〈L_B〉=60 m and 〈L_O〉=56 m, which are close to L_T (64 m) (Table 4b). Therefore, 〈L_H〉 and 〈L_C〉 should not be the turbulence scales to consider. From the Thorpe analysis, $N^{\mathrm{2}}\approx \mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻² and S≈0.0026 s⁻¹ (Ri≈2). From Eq. (6) with c=1, ε_T≈0.2 mW kg⁻¹ (Table 3b), i.e., very close to the mean value of DataHawk-derived ε or only 2 times lower than the value during A1. It is consistent with the fact that L_T can be assimilated to L_O. From Eq. (7) with Ri=2, $c_{\mathrm{Sc}}^{\prime }=\mathrm{0.41}$ and $c_{\mathrm{BP}}^{\prime }=\mathrm{0.24}$, we obtain ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }\approx \mathrm{0.012}$ and 0.004 mW kg⁻¹, respectively, which is much less than 0.2 mW kg⁻¹. The various estimates of ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }$ are shown in Fig. 3d. As expected ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }$ fails because it is expected to be valid for Ri<0.25 only.

5 Statistical analysis

5.1 Justification of the application of L_out=70 m in Eq. (2)

Figure 5 shows the histogram of log ₁₀(L), where $L=\langle {\mathit{\sigma }}^{\mathrm{2}}{\rangle }^{\mathrm{3}/\mathrm{2}}/{\mathit{\epsilon }}_{\mathrm{U}}$ for $\langle {\mathit{\sigma }}^{\mathrm{2}}{\rangle }^{\mathrm{3}/\mathrm{2}}>\mathrm{0.01}$ as in L18 obtained from the WPR-LQ-7 from data collected during 36 flights (corresponding to 90 profiles). The peak of the distribution has a mean (median) value of 67 m (71 m). These values are almost identical to those obtained from comparisons with MU radar, i.e., 75 m (61 m) (Fig. 7a of L18). This result seems to indicate that the empirical expression ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ with L_out=70 m is not specific to the MU radar but at least to any radar with similar resolution volume (Table 2). However, the acquisition time of the MU radar and WPR-LQ-7 data differs by a factor of ∼ 2.5 (Table 2). It is likely fortunate that the statistical values of ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ (and thus σ) coincide so well.

Figure 5 Histogram of log ₁₀(L_out) for $\langle \mathit{\sigma }{\rangle }^{\mathrm{3}/\mathrm{2}}>\mathrm{0.01}$ as in L18 for MU radar data.

5.2 Statistical evaluation of the models from comparisons with ε_U

Figure 6a, b and c show the scatter plots of log ₁₀(ε_U) vs ${\log }_{\mathrm{10}}\left({\mathit{\epsilon }}_{L_{\mathrm{out}}}\right)$, log ₁₀(ε_S) and log ₁₀(ε_N) with L_out=70 m, the shear estimated from WPR-LQ-7 data and N from DataHawk data at a height resolution of 100 m. Figure 6d and e show the results assuming a constant shear (〈S〉=7.7 m s⁻¹ km⁻¹) and a constant N ($\langle N^{\mathrm{2}}\rangle =\mathrm{6.7}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻²). Of course, Fig. 6d and e differ only in the constant $\mathrm{0.64}/\mathrm{0.5}$.

Figure 6 Scatter plots of log ₁₀(ε_U) vs (a) ${\log }_{\mathrm{10}}\left({\mathit{\epsilon }}_{L_{\mathrm{out}}}\right)$, (b) log ₁₀(ε_S), (c) log ₁₀(ε_N), (d) log ₁₀(ε_S) with constant S and (e) log ₁₀(ε_N) with constant N. The red lines are the result of a line regression (whose slope value is indicated in the panel) for ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³, and R is the correlation coefficient.

The correlation coefficients are fortuitously ∼ 0.66 for all the cases except for ε_N for which the correlation is ∼ 0.60 only. This is an additional clue of the inadequacy of ε_N. The red lines show the results of linear regressions after rejecting dissipation rate values smaller than $\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³ as in L18, even though the quantitative threshold has no reason to be the same since the comparison methods differ. The slope of the regression line between ${\log }_{\mathrm{10}}\left({\mathit{\epsilon }}_{L_{\mathrm{out}}}\right)$ and log ₁₀(ε_U) is ∼ 1.0 (Fig. 6a), confirming the statistical σ³ dependence of ε when no discrimination is made on the conditions under which turbulence occurs. The regression slope obtained with log ₁₀(ε_N) or log ₁₀(ε_S) for S or N equal to a constant (Fig. 6d, e) is 0.60, i.e., close to $\mathrm{2}/\mathrm{3}$, as expected because the two models vary as σ². However, the regression slope between log ₁₀(ε_U) and log ₁₀(ε_N) with measured N (Fig. 6c) is significantly lower than 0.66 (0.50), and the regression slope between log ₁₀(ε_U) and log ₁₀(ε_S) with measured S (Fig. 6b) is significantly larger than 0.66 (0.73). The regression slope between log ₁₀(ε_U) and log ₁₀(ε_N) using MU radar data was 0.55 (L18), i.e., virtually identical to the present case (Fig. 6c). All the regression slopes depend on the quantitative threshold on ε, and, in Fig. 6a, it varies from ∼ 0.9 to ∼ 1.1 for different thresholds, excluding small values. However, all other slope estimates vary in concert, so that the observed trends remain valid for a different threshold.

• 1.

  A regression slope between log ₁₀(ε_U) and log ₁₀(ε_S) that is closer to 1 than the slope between log ₁₀(ε_U) and log ₁₀(ε_N) indicates that ε_S provides estimates more consistent with ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ than ε_N. Figure 7a and b show a comparison between the radar models, i.e., ${\log }_{\mathrm{10}}\left({\mathit{\epsilon }}_{L_{\mathrm{out}}}\right)$ vs log ₁₀(ε_N) and ${\log }_{\mathrm{10}}\left({\mathit{\epsilon }}_{L_{\mathrm{out}}}\right)$ vs log ₁₀(ε_S), respectively, for ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³. The red lines are the results for constant N and constant S. The blue lines show the results of a linear regression. The blue and red lines obviously coincide (slope =0.60) for ε_N. A slope of 0.77 is obtained with ε_S, indicating a greater equivalence between ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ and ε_S, as expected from Fig. 6b. Consequently, our results suggest that ε_S is more relevant than ε_N and should be used instead of ε_N for operational use if the empirical model ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ is not chosen.

• 2.

  We checked that generating a normal random distribution of N or S with a mean and standard deviation similar to the observed distributions produced a regression slope close to $\mathrm{2}/\mathrm{3}$, similar to Fig. 6d and e. Therefore, the observed slopes with measured S and N (Fig. 6b, c) must reveal a statistical dependence of σ with $\mathrm{1}/N$ and S, respectively. The equivalence between ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ and ε_S described in Sect. 4 for the K–H layer implies that σ is simply proportional to S (σ=0.64 L_outS) if L_out is constant. There is a canonical value of L_out (∼70 m), but since L_out is not constant and is unknown (and can vary by 2 orders of magnitude at least; Fig. 5), the correlation between σ and S can only be established for a fixed value of L_out (and for any other variable on which σ depends). Figure 8 shows the same information as Fig. 7 but after dividing by σ² to remove the self-correlation between the variables and to show the relationship between ${\log }_{\mathrm{10}}(\mathit{\sigma }/L_{\mathrm{out}})$ and log ₁₀(0.5N) (Fig. 8a) and between ${\log }_{\mathrm{10}}(\mathit{\sigma }/L_{\mathrm{out}})$ and log ₁₀(0.64S) (Fig. 8b). The two scatter plots show negative and positive correlation coefficients (−0.36 and 0.12, respectively). The correlations are weak but significant according to the P value. If no threshold on ε_U is applied, the correlation coefficients are −0.26 and +0.22. This suggests that σ increases to some extent as S increases and N decreases. It is quite intuitive, but, to our knowledge, this is the first study to suggest and highlight this. The results may reveal a Richardson number dependence. Figure 9a shows scatter plots of σ vs ${\mathit{Ri}}_{\mathrm{100}}^{\mathrm{1}/\mathrm{2}}$, where Ri₁₀₀ (S₁₀₀) now explicitly refers to the Richardson number (shear) calculated at the vertical resolution of 100 m. The red (black) dots show the results without and with a threshold on S₁₀₀ (S₁₀₀>5 m s⁻¹ km⁻¹), respectively. The (negative) correlation coefficient is slightly stronger with the threshold (−0.34 instead of −0.21). The high values of Ri₁₀₀ are mainly associated with a weak shear (S₁₀₀<5 m s⁻¹ km⁻¹) and with the largest variability in σ (Fig. 9a). However, this property may not be significant because the uncertainty on Ri increases as the shear tends to zero. Therefore, we focus on the scatter plot obtained with the threshold on the shear (black dots). It seems to show a linear dependence between log ₁₀(σ) and ${\log }_{\mathrm{10}}\left({\mathit{Ri}}_{\mathrm{100}}^{-\mathrm{1}/\mathrm{2}}\right)$, at least down to ${\log }_{\mathrm{10}}\left({\mathit{Ri}}_{\mathrm{100}}^{\mathrm{1}/\mathrm{2}}\right)\approx -\mathrm{0.2}$, i.e., for Ri₁₀₀<0.4. Attempts of linear regression analysis do not confirm the linear trend, likely due to the strong dispersion and weak correlation. However, the time series obtained from the concatenation of all the profiles of ${\log }_{\mathrm{10}}\left({\mathit{Ri}}_{\mathrm{100}}^{-\mathrm{1}/\mathrm{2}}\right)$ and log ₁₀(σ) after removing their mean values reveal a more obvious dependence between the two variables (Fig. 9b). The curves reveal similar variations and dynamics, especially for records [0−200], compatible with σ² inversely proportional to Ri₁₀₀, at least to a first approximation. For ${\log }_{\mathrm{10}}\left({\mathit{Ri}}_{\mathrm{100}}^{\mathrm{1}/\mathrm{2}}\right)\lesssim -\mathrm{0.2}$, i.e., for low values of Ri₁₀₀(<0.4), log ₁₀(σ) appears to have very little dependence on ${\log }_{\mathrm{10}}\left({\mathit{Ri}}_{\mathrm{100}}^{\mathrm{1}/\mathrm{2}}\right)$. If meaningful, it would be consistent with the fact that N does not play a significant role for low Ri values, as suggested by the ε_S model.

Figure 7 Scatter plots of (a) ${\log }_{\mathrm{10}}\left({\mathit{\epsilon }}_{L_{\mathrm{out}}}\right)$ vs log ₁₀(ε_N) and (b) ${\log }_{\mathrm{10}}\left({\mathit{\epsilon }}_{L_{\mathrm{out}}}\right)$ vs log ₁₀(ε_S) for ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³. The red lines show the results for constant N (a) and constant S (b). The blue lines show the results of a linear regression (whose slope value is indicated in the panel).

Figure 8 Scatter plots of (a) log ₁₀(0.5N) vs ${\log }_{\mathrm{10}}(\mathit{\sigma }/L_{\mathrm{out}})$ and (b) log ₁₀(0.64S) vs ${\log }_{\mathrm{10}}(\mathit{\sigma }/L\mathrm{out})$ for ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³. R is the correlation coefficient, and P is the result of the P test.

Figure 9 (a) Scatter plots of ${\log }_{\mathrm{10}}\left({\mathit{Ri}}_{\mathrm{100}}^{\mathrm{1}/\mathrm{2}}\right)$ vs log ₁₀(σ) for ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³ without the threshold on the shear (grey) and for S₁₀₀>5 m s⁻¹ km⁻¹. (b) The corresponding time series of ${\log }_{\mathrm{10}}\left({\mathit{Ri}}_{\mathrm{100}}^{-\mathrm{1}/\mathrm{2}}\right)$ (grey) and log ₁₀(σ) (black) after subtracting their mean.

Figure 10 shows the scatter plots of ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_{L_{\mathrm{out}}}/{\mathit{\epsilon }}_{\mathrm{U}})$, ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_N/{\mathit{\epsilon }}_{\mathrm{U}})$ and ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_{\mathrm{S}}/{\mathit{\epsilon }}_{\mathrm{U}})$ vs log ₁₀(Ri₁₀₀) applying two thresholds on ε_U: $\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³ in Fig. 10a, b and c and $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³ in Fig. 10d, e and f. The latter is introduced to analyze the dependence of the results on the levels of considered dissipation rates. The red and blue curves show the values averaged in bandwidths of 0.3 from ${\log }_{\mathrm{10}}\left({\mathit{Ri}}_{\mathrm{100}}\right)=-\mathrm{1.7}$. For ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³, the mean curves of ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_{L_{\mathrm{out}}}/{\mathit{\epsilon }}_{\mathrm{U}})$ and ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_{\mathrm{S}}/{\mathit{\epsilon }}_{\mathrm{U}})$ do not reveal a significant dependence on log ₁₀(Ri₁₀₀), at least up to log ₁₀(Ri₁₀₀)∼1, and are almost identical and close to 0 (Fig. 10a, c). Therefore, the applicability of the two models does not seem to depend significantly on the Richardson number on average. For ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³, the curves produced by two models remain close and almost unchanged for log ₁₀(Ri₁₀₀)<0 (Fig. 10d, f). However, the mean values of ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_{\mathrm{S}}/{\mathit{\epsilon }}_{\mathrm{U}})$ now tend to decreases as log ₁₀(Ri₁₀₀) increases. Therefore, when log ₁₀(Ri₁₀₀)>0, ε_S tends to underestimate ε_U when ε_U exceeds ∼ $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³ and inversely when ${\mathit{\epsilon }}_{\mathrm{U}}<\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³. The fact that we experimentally observe that ε_S is not appropriate for large values of Ri₁₀₀ is consistent with the expected domain of applicability of the model, even if it is not clear why it is in this way. For log ₁₀(Ri₁₀₀)<0, ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_N/{\mathit{\epsilon }}_{\mathrm{U}})$ is less than 0 and decreases as log ₁₀(Ri₁₀₀) decreases for both thresholds (Fig. 10b, e). This experimental observation is a confirmation of the inadequacy of ε_N when the Richardson number is low. The results with ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ are difficult to interpret when log ₁₀(Ri₁₀₀)>0. The model is consistent with ε_S when ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³ (Fig. 10c) and seems to be more consistent with ε_N than with ε_S when ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³ (Fig. 10e, f). It may be vain to interpret the properties of this model, since it is only an empirical model for which L_out=70 m only represents a canonical value of a function of multiple variables, including the shear and N.

Figure 10 Scatter plots of log ₁₀(Ri₁₀₀) vs (a, d) ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_{L_{\mathrm{out}}}/{\mathit{\epsilon }}_{\mathrm{U}})$, (b, e) ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_N/{\mathit{\epsilon }}_{\mathrm{U}})$ and (c, f) ${\log }_{\mathrm{10}}({\mathit{\epsilon }}_{\mathrm{S}}/{\mathit{\epsilon }}_{\mathrm{U}})$ for ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³ and ${\mathit{\epsilon }}_{\mathrm{U}}>\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³, respectively.

6 Conclusions

The objective of this work was to test the suitability of TKE dissipation rate models based on Doppler radar spectral width measurements from comparisons with in situ estimates (ε_U) derived from high-resolution Pitot tube measurements aboard DataHawk UAVs. We showed the following:

• 1.

  The models applied to the 46.5 MHz MU VHF radar by L18 produce statistically identical results on the 1.357 GHz WPR-LQ-7.

  • a.

    The empirical model ${\mathit{\epsilon }}_{L_{\mathrm{out}}}={\mathit{\sigma }}^{\mathrm{3}}/L_{\mathrm{out}}$ with L_out∼70 m (as for the MU radar) provides the best statistical agreement with ε_U, at least for $\mathit{\epsilon }\gtrsim \mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³ (Table 2). If L_out really depends on the size of energy containing eddies, it is then independent of radar parameters (assuming σ² is a true indication of 〈w^′2〉 in both radars).

  • b.

    The model ε_N predicting a σ²N law for stably stratified conditions fails to reproduce ε_U. The biases are nearly quantitatively identical to those obtained with the MU radar: ε_N tends to overestimate when ${\mathit{\epsilon }}_{\mathrm{U}}\lesssim \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³ and to underestimate when ${\mathit{\epsilon }}_{\mathrm{U}}\gtrsim \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³.

• 2.

  Applying ε_N to both radars to a turbulent layer attributed to a K–H instability with Ri<0.25 strongly underestimates ε_U in the core of the layer when N² is minimum. On the other hand, in agreement with the statistical results, ε_N provided values consistent with the other estimates in a turbulent layer, likely associated with larger Ri (≳1). These two observations are rather consistent with the domain of validity of ε_N according to the theoretical derivations (Eq. 5), leading to the newly introduced expression of ε_S expected to be valid for weakly stratified or strongly sheared conditions (e.g., Basu et al., 2021).

• 3.

  The application of ε_S to the K–H layer (Ri<0.25) using both radars leads to a good agreement with ε_U. Its application to the turbulent layer associated with larger Ri slightly underestimates ε_U, again in accordance with Eq. (5).

• 4.

  The statistical comparisons between ε_S and ε_U using all data show much better agreement than between ε_N and ε_U, although a bias of the same nature as that observed with ε_N is also noted but to a lesser degree. Empirical ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ remains the most consistent model compared with ε_U. L_out∼70 m is likely a canonical value that results from all the hidden contributions of the various parameters that a most general (and unknown) model should include.

• 5.

  The equivalence between ε_S and ε_U for the K–H layer associated with a low Ri necessarily implies that σ is proportional to S: σ∼0.64 L_outS. For all the layers with the same value of L_out, then σ linearly depends on S. This is a necessary condition if agreement is observed with two models that predict a σ³ and a σ² dependence. For a wide distribution of log ₁₀(L_out) as in Fig. 5, which includes values for all Ri, this linear dependence should be strongly “blurred” because L_out is variable, and Ri is not necessarily low. Moreover, an additional source of dispersion is that the wind shear calculated at the radar resolution S₁₀₀ and at a time resolution of 10–30 min is not necessarily the most effective shear to be considered because S is a scale-dependent parameter (in the same way as Ri). As a result, a very weak, but yet significant, correlation between σ and S₁₀₀ was found (Fig. 8). This weak correlation is responsible for the better agreement obtained with ε_S than with ε_N. However, because the results are based on a limited amount of data and because some degree of coincidence cannot be ruled out, further studies are necessary to analyze the dependence between σ and S₁₀₀, in particular under more suitable conditions (i.e., low Ri).

• 6.

  Reciprocally, σ has a statistical degree of dependence on $\mathrm{1}/N$, as revealed by the poorer statistical agreement between ε_N and ε_U, leading to a regression slope of less than $\mathrm{2}/\mathrm{3}$ (0.50, almost identical to 0.55 obtained from MU radar data by L18) (Fig. 6).

• 7.

  The combination of points (5) and (6) leads to the conclusion that, to some extent, σ depends on ${\mathit{Ri}}_{\mathrm{100}}^{-\mathrm{1}/\mathrm{2}}$, at least for Ri₁₀₀≳0.4 (Fig. 9). This dependence does not seem to be valid for lower Ri₁₀₀ (Fig. 8a), in accordance with the fact that N should not affect turbulence when the Richardson number is low (Eq. 4).

• 8.

  The analysis of the three models, ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$, ε_N and ε_S, with ε_U vs Ri₁₀₀ (Fig. 10) confirms the good agreement between (${\mathit{\epsilon }}_{L_{\mathrm{out}}}$, ε_U) and between (ε_S, ε_U) and the inadequacy of ε_N for Ri₁₀₀≲1. The underestimation of ε_N increases as Ri₁₀₀ decreases. The results for Ri₁₀₀≳1 are more difficult to interpret and more puzzling, but ε_S and ${\mathit{\epsilon }}_{L_{\mathrm{out}}}$ lead to comparable results and do not show substantial bias as a function of Ri₁₀₀. In any case, all results involving large Ri (>1) must be taken with caution because the turbulence may be intermittent. In principle, the interpretation of the results should consider this.

• 9.

  The classical model ${\mathit{\epsilon }}_{\mathrm{T}}=c^{\mathrm{2}}{L}_{\mathrm{T}}^{\mathrm{2}}{N}^{\mathrm{3}}$ (Eq. 6) based on the equivalence between the Thorpe length L_T and the Ozmidov length scale L_O (c=1) fails to reproduce DataHawk-derived ε in the K–H layer for which Ri is expected to be less than 0.25. Although the disagreement can be due to several factors (e.g., an inappropriate choice of c, horizontal inhomogeneity), it can also be due to the fact the model involves the Ozmidov length scale defined for a turbulence affected by the stable stratification. In essence, L_T cannot be related to L_O anymore by a constant if the stratification effects can really be neglected for low Ri. Therefore, an alternative approach using the Corrsin length scale L_C instead of L_O was introduced, leading to ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }=c^{\prime \mathrm{2}}{L}_{\mathrm{T}}^{\mathrm{2}}{S}^{\mathrm{3}}$ (Eq. 7), compatible with studies showing a ${\mathit{Ri}}^{\mathrm{3}/\mathrm{4}}$ dependence of $L_{\mathrm{T}}/L_{\mathrm{O}}$ for Ri<0.25. Contrary to c, c^′ is a true constant (with respect to Ri) for low Ri. It is worth noting that Eqs. (7) and (4) form a coherent pair of models independent of N for a weak stratification or strongly sheared flows. Using values of c^′ deduced from the literature, ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }$ provides estimates consistent with DataHawk- and radar-derived ε (expect ε_N) for the K–H layer. On the other hand, ${\mathit{\epsilon }}_{\mathrm{T}}^{\prime }$ fails for a decaying turbulent layer (Ri>1) as the model is not expected to be valid for Ri>0.25. ε_T with c=1 shows a better agreement with DataHawk- and radar-derived ε (including ε_N), coherent with the fact that the stable stratification should affect the turbulence for large Ri. These results need to be confirmed by statistical studies.

Appendix A

The turbulent Froude number $\mathit{Fr}=\mathit{\epsilon }/Nk$, where k refers to TKE for a simple and standard notation, is often used to characterize turbulent mixing (e.g., Ivey and Imberger, 1991). A strong and weak stratification is associated with Fr<1 and Fr>1, respectively. $T_{\mathrm{L}}=k/\mathit{\epsilon }$ is called the inertial timescale and is a characteristic time of TKE dissipation. The corresponding timescales associated with the turbulence production by the wind shear and with the conversion into potential energy are S⁻¹ and N⁻¹, respectively. When the stratification is weak, i.e., when $N{T}_{\mathrm{L}}={\mathit{Fr}}^{-\mathrm{1}}<\mathrm{1}$, then T_L (dissipation time) and S⁻¹ (production time) should be of the same order, i.e., the shear parameter ST_L=O(1), for stationary turbulence. Several studies (e.g., Mater and Venayagamoorthy, 2014, and references therein) reported a critical value for weakly stratified and stationary flows:

$$\begin{array}{}\text{(A1)}& {\text{ST}}_{\mathrm{Lc}}\approx \mathrm{3.33}.\end{array}$$

By dividing Eq. (A1) by 3.33 NT_Lc, we obtain

$$\begin{array}{}\text{(A2)}& \mathrm{0.3}S/N=\mathrm{1}/N{T}_{\mathrm{Lc}}\left(\iff \mathit{Fr}=\mathrm{0.3}/\sqrt{\mathit{Ri}}\right).\end{array}$$

From the definition of Fr, Eq. (A2) reads

$$\begin{array}{}\text{(A3)}& \mathit{\epsilon }=\mathrm{0.30}kS.\end{array}$$

Assuming isotropy, $k=\mathrm{3}/\mathrm{2}\langle w^{\prime \mathrm{2}}\rangle$, where 〈w^′2〉 is the vertical velocity variance (assumed to be σ² in the paper). We obtain

$$\begin{array}{}\text{(A4)}& \mathit{\epsilon }=\mathrm{0.45}\langle w^{\prime \mathrm{2}}\rangle S.\end{array}$$

For (anisotropic) shear-generated turbulence, $k\approx \mathrm{2}\langle w^{\prime \mathrm{2}}\rangle$, so that

$$\begin{array}{}\text{(A5)}& \mathit{\epsilon }=\mathrm{0.60}\langle w^{\prime \mathrm{2}}\rangle S,\end{array}$$

i.e., virtually Eq. (4) with C_S=0.63. These expressions are valid for Fr>1, i.e., Ri<0.09, according to Eq. (A2).

For $k\approx \mathrm{2.74}\langle w^{\prime \mathrm{2}}\rangle$ when $\mathrm{0}<\mathit{Ri}<\mathrm{0.2}$ (Eq. 28, Basu and Holtslag, 2021), we get

$$\begin{array}{}\text{(A6)}& \mathit{\epsilon }=\mathrm{0.82}\langle w^{\prime \mathrm{2}}\rangle S.\end{array}$$

Appendix B

From Fig. (3) of Garanaik and Venayagamoorthy (2019), showing the turbulent Froude number Fr vs $L_{\mathrm{E}}/L_{\mathrm{O}}$ (or $L_{\mathrm{T}}/L_{\mathrm{O}}$) from DNS, we obtain, for weakly stratified conditions (Fr>1),

$$\begin{array}{}\text{(B1)}& \mathit{Fr}=\mathit{\alpha }{\left(L_{\mathrm{T}}/L_{\mathrm{O}}\right)}^{-\mathrm{2}/\mathrm{3}},\end{array}$$

with α≈0.7. Note that this coefficient is deduced from their Fig. (3) using the linear trend shown by the authors. They did not explicitly refer to this value. By using the definitions $\mathit{Fr}=\mathit{\epsilon }/Nk$ (see Appendix A) and $L_{\mathrm{O}}=\sqrt{\mathit{\epsilon }/N^{\mathrm{3}}}$, Eq. (B1) can be re-written as

$$\begin{array}{}\text{(B2)}& \mathit{\epsilon }={\left(\mathit{Fr}/\mathit{\alpha }\right)}^{\mathrm{3}}{L}_{\mathrm{T}}^{\mathrm{2}}{N}^{\mathrm{3}}.\end{array}$$

For Fr>1 or ${\text{NT}}_{\mathrm{L}}={\mathit{Fr}}^{-\mathrm{1}}<\mathrm{1}$, ST_L≈3.33 (see Appendix 1 and Fig. 1 of Mater and Venayagamoorthy, 2014). Therefore, $\mathit{Fr}\approx \mathrm{0.3}S/N$, so that

$$\begin{array}{}\text{(B3)}& \mathit{\epsilon }={\left(\mathrm{1}/\mathit{\alpha }\right)}^{\mathrm{3}}{L}_{\mathrm{T}}^{\mathrm{2}}{S}^{\mathrm{3}}=c^{\prime \mathrm{2}}{L}_{\mathrm{T}}^{\mathrm{2}}{S}^{\mathrm{3}},\end{array}$$

with $c^{\prime }=\mathrm{0.28}$. This value agrees well with the values reported in the main text.

In addition, by replacing Fr by its definition, Eq. (B2) can be re-written as $\mathit{\epsilon }={\mathit{\alpha }}^{-\mathrm{3}}{L}_{\mathrm{T}}^{\mathrm{2}}{\mathit{\epsilon }}^{\mathrm{3}}/k^{\mathrm{3}}$, so that

$$\begin{array}{}\text{(B4)}& \mathit{\epsilon }\approx k^{\mathrm{3}/\mathrm{2}}/\left(\mathrm{1.5}{L}_{\mathrm{T}}\right).\end{array}$$

Equation (B4) provides a way to relate the isotropic turbulent length scale L_k defined as the scale of the largest eddies weakly affected by the buoyancy and the shear to the Thorpe length: L_k∼1.5 L_T. It also provides an expression of the master length scale L_M defined as $(\mathrm{2}k{)}^{\mathrm{3}/\mathrm{2}}/B_{\mathrm{1}}\mathit{\epsilon }$ (e.g., Mellor and Yamada, 1982) with $\mathrm{11.9}\le B_{\mathrm{1}}\le \mathrm{27.4}$ (Table 2, Basu and Holtslag, 2021). We obtain $L_{\mathrm{M}}=\mathrm{4.24}/B_{\mathrm{1}}{L}_{\mathrm{T}}$.

On the other hand, for strongly stratified conditions (Fr<1), we can write, from Fig. (3) of Garanaik and Venayagamoorthy (2019),

$$\begin{array}{}\text{(B5)}& \mathit{Fr}={\mathit{\alpha }}^{\prime }{\left(L_{\mathrm{T}}/L_{\mathrm{O}}\right)}^{-\mathrm{2}},\end{array}$$

with ${\mathit{\alpha }}^{\prime }\approx \mathrm{0.6}$ according to the figure, or, purposely, $\mathrm{0.66}=\mathrm{2}/\mathrm{3}$. We obtain

$$\mathit{\epsilon }=\mathrm{3}/\mathrm{2}\mathit{Fr}{L}_{\mathrm{T}}^{\mathrm{2}}{N}^{\mathrm{3}},$$

so that

$$\begin{array}{}\text{(B6)}& k=\mathrm{3}/\mathrm{2}{L}_{\mathrm{T}}^{\mathrm{2}}{N}^{\mathrm{2}}.\end{array}$$

For $k=\mathrm{3}/\mathrm{2}\langle w^{\prime \mathrm{2}}\rangle$, then, $\langle w^{\prime \mathrm{2}}\rangle =L_{\mathrm{T}}^{\mathrm{2}}{N}^{\mathrm{2}}$ or

$$\begin{array}{}\text{(B7)}& L_{\mathrm{B}}=\sqrt{\langle w^{\prime \mathrm{2}}\rangle }/N=L_{\mathrm{T}}.\end{array}$$

If $k=\mathit{\beta }\langle w^{\prime \mathrm{2}}\rangle$, then $L_{\mathrm{B}}=\sqrt{\mathrm{1}/\left({\mathit{\alpha }}^{\prime }\mathit{\beta }\right)}{L}_{\mathrm{T}}$. Equation (B5) demonstrates that the equivalence (at least the proportionality) between the buoyancy length scale and the Thorpe length is valid for strongly stratified conditions only (i.e., for conditions for which stability affects the vertical motions before being affected by the wind shear). For a weak stratification, the vertical TKE cannot be fully converted into potential energy because the parcels cannot move vertically over a length of L_B (but L_H only). Therefore, the basic stability does not intervene anymore in the variance of the vertical velocity fluctuations as in the case of a neutral stratification.

Appendix C: List of symbols

α:	constant defined in Eq. (B1)
α′:	constant defined in Eq. (B5)
a.g.l.:	above ground level
a.s.l.:	above sea level
B:	buoyancy production/destruction term in Eq. (1)
β:	coefficient of proportionality between TKE and the variance 〈w′2〉 defined in Appendix B.
c:	constant in Eq. (6)
c′:	constant in Eq. (7)
$c_{\mathrm{Sc}}^{\prime }$:	constant in Eq. (7) for Schumann (1994) model
$c_{\mathrm{BP}}^{\prime }$:	constant in Eq. (7) for Baumert and Peters (2000) model
CN:	constant in Eq. (3)
CS:	constant in Eq. (4)
$C_{\mathrm{S}}^{\prime }$:	pseudo-constant (depending on Rf) in Eq. (5a)
$C_{\mathrm{S}}^{\prime \prime }$:	pseudo-constant (depending on Rf and G) in Eq. (5b)
DH:	DataHawk UAV
DNS:	direct numerical simulation
d′:	Thorpe displacement
ε:	TKE dissipation rate (general term)
${\mathit{\epsilon }}_{L_{\mathrm{out}}}$:	ε estimated from radar data from Eq. (2)
εN:	ε estimated from radar data from Eq. (3)
εS:	ε estimated from radar data from Eq. (4)
${\mathit{\epsilon }}_{\mathrm{S}}^{\prime }$:	ε estimated from radar data from Eq. (5a)
εT:	ε estimated from balloon data from Eq. (6)
${\mathit{\epsilon }}_{{\mathrm{T}}^{\prime }}$:	ε estimated from balloon data from Eq. (7)
εU:	ε estimated from UAV data
F:	fraction of the inertial and buoyancy subrange contributions in Eq. (3)
FE:	full equilibrium
Fr:	Froude number
FFT:	fast Fourier transform
g:	gravitation acceleration
G:	growth factor
G0:	growth factor for Ri=0
k:	an alternative notation for TKE (Appendix)
kB:	buoyancy wavenumber
kH:	Hunt wavenumber
KHI:	Kelvin–Helmholtz instability
L:	$\langle {\mathit{\sigma }}^{\mathrm{2}}{\rangle }^{\mathrm{3}/\mathrm{2}}/{\mathit{\epsilon }}_{\mathrm{U}}$ (Sect. 5.1)
LB:	buoyancy length scale
LC:	Corrsin length scale
LES:	large eddy simulation
LE:	Ellison length scale
LH:	Hunt length scale
LM:	master length scale defined in Appendix B
Lout:	an apparent outer scale of turbulence

LO:	Ozmidov length scale
LT:	Thorpe length
MU:	middle and upper atmosphere (radar)
N:	Brünt–Vaïsälä frequency
Ns:	Brünt–Vaïsälä frequency estimated from sorted θ profile
Nrms:	Brünt–Vaïsälä frequency estimated from the rms θ fluctuations
NINCOH:	number of incoherent integrations
P:	shear production term in Eq. (1)
Pr:	turbulent Prandtl number
R:	correlation coefficient
Rf:	flux Richardson number
Ri:	Richardson number
Ric:	critical Richardson number
Ris:	Richardson number at the FE condition
Ri100:	Richardson number estimated at a vertical resolution of 100 m
S:	vertical shear of horizontal wind
S100:	S estimated at a vertical resolution of 100 m
ShUREX:	Shigaraki UAV Radar EXperiment
σ2:	radar estimate of the variance 〈w′2〉 of the vertical wind fluctuations produced by turbulence
σobs:	half the measured Doppler spectral width
SNR:	signal-to-noise ratio
θ:	potential temperature
θ0:	half-power half width of the two-way radar beam
TKE:	turbulence kinetic energy
TL:	inertial timescale
τ:	radar pulse duration
UAV:	uncrewed aerial vehicle
UHF:	ultrahigh frequency
VHF:	very high frequency
WPR-LQ-7:	LQ-7 wind profiler
z:	altitude (a.g.l.)