1 Introduction

The molecular dissipation of temperature fluctuations (\(\chi\)) is an important variable for characterizing turbulent mixing in various environmental flows. It is frequently used in micrometeorology (e.g., [66]) and atmospheric optics (e.g., [43]). Furthermore, any higher-order closure model requires solving a prognostic equation or a diagnostic parameterization for ensemble-averaged \(\overline{\chi }\) (refer to [15, 36, 40, 65]) .

Over the years, a number of studies focused on the correlation between turbulent kinetic energy dissipation rate (\(\varepsilon\)) and \(\chi\) (e.g., [2, 5,6,7, 29, 69]). In addition, some papers reported on the probability density function, spatio-temporal intermittency and anomalous scaling of \(\chi\) (e.g., [6, 52, 54]). Often, \(\chi\) has been found to be more intermittent (commonly quantified by the multifractal scaling exponents) and non-Gaussian than \(\varepsilon\) (e.g., [50, 54]).

Most of these previous studies primarily focused on the instantaneous, localized traits of the dissipation fields. Instead, we are interested to better quantify their spatially averaged characteristics. Towards this goal, we first investigate the statistical properties of several length scales which can be formulated based on \(\overline{\chi }\) and other relevant variables. Based on these findings, we then derive simple parameterizations for \(\overline{\chi }\) and temperature structure parameter (\(C_T^2\)). For all the analyses, we utilize a direct numerical simulation (DNS) database of stably stratified flows which is discussed in the following section.

2 Direct numerical simulation

Recently, for the parameterization of optical turbulence, He and Basu [31] created a DNS database using a massively parallel DNS code, called HERCULES [30]. The DNS runs were conducted by solving the normalized Navier–Stokes and temperature equations, as shown in Eqs. 1–3 (using Einstein’s summation notation):

$${\frac{\partial u_{n,i}}{\partial x_{n,i}}=0,}$$
(1)
$${\frac{\partial u_{n,i}}{\partial t_n}+\frac{\partial u_{n,i} u_{n,j}}{\partial x_{n,j}}=-\frac{\partial p_n}{\partial x_{n,i}} +\frac{1}{Re_b} \frac{\partial }{\partial x_{n,j}} \frac{\partial u_{n,i}}{\partial x_{n,j}} +{\varDelta }P {\delta }_{i1}+Ri_b {\theta_n } {\delta }_{i3}, } $$
(2)
$${\frac{\partial \theta_n }{\partial t_n}+\frac{\partial \theta_n u_{n,j}}{\partial x_{n,j}}=\frac{1}{Re_b Pr} \frac{\partial }{\partial x_{n,j}} \frac{\partial \theta_n }{\partial x_{n,j}},} $$
(3)

where \(u_n\) and \(x_n\) are the normalized velocity and coordinate vectors, respectively, with the subscript  i  denoting the ith vector component, tn is the normalized time, pn is the normalized pressure, \({\varDelta }P\) is the streamwise pressure gradient driving the flow, and \(\theta_n\) is the normalized potential temperature. The bulk Richardson number is denoted by:

$$\begin{aligned} Ri_b = \frac{\left( \varTheta _{top}-\varTheta _{bot}\right) g h}{U_b^2 \varTheta _{top}}, \end{aligned}$$
(4)

where g denotes the gravitational acceleration, and \(\varTheta _{top}\) and \(\varTheta _{bot}\) represent potential temperature at the top and the bottom of the channel, respectively. \(Pr = \nu /k = 0.7\) is the Prandtl number with k being the thermal diffusivity, and \(Re_b = \frac{U_b h}{\nu }\) is the bulk Reynolds number with h, \(U_b\), and \(\nu\) being the channel height, the bulk (averaged) velocity in the channel, and the kinematic viscosity, respectively. The bulk Reynolds number was fixed at 20,000 for all the simulations.

The computational domain size for all the DNS runs was \(L_x \times L_y \times L_z = 18 h \times 10 h \times h\). The domain was discretized by \(2304 \times 2048 \times 288\) grid points in streamwise, spanwise, and wall-normal directions, respectively. A total of five simulations were performed with gradual increase in the temperature difference between the top and bottom walls (effectively by increasing \(Ri_b\)) to mimic the nighttime cooling of the land-surface. The normalized cooling rates (CR), \(\partial Ri_b/ \partial T_n\), ranged from \(1\times 10^{-3}\) to \(5\times 10^{-3}\); where, \(T_n\) is a non-dimensional time (\(=tU_b/h\)).

All the simulations used fully developed neutrally stratified flows (\(Ri_b = 0\)) as initial conditions and evolved for up to \(T_n = 100\). The simulation results were output every 10 non-dimensional time. To avoid spin-up issues, in the present study, we only use data for the last five output files (i.e., \(60 \le T_n \le 100\)). Furthermore, we only consider data from the region \(0.1 h\le z \le 0.5 h\) to discard any blocking effect of the surface or avoid any laminarization in the upper part of the open channel. Vertical profiles and some basic statistics from these simulations have been documented in “Appendix 3”.

The mean dissipation of turbulent kinetic energy and temperature fluctuations are computed as follows:

$$\begin{aligned} \overline{\varepsilon }&= \nu \overline{\left( \frac{\partial u_i'}{\partial x_j} \frac{\partial u_i'}{\partial x_j}\right) }, \end{aligned}$$
(5a)
$$\begin{aligned} \overline{\chi }&= 2k \overline{\left( \frac{\partial \theta '}{\partial x_j} \frac{\partial \theta '}{\partial x_j} \right) }. \end{aligned}$$
(5b)

In the above equations, and in the rest of the paper, the “overbar” notation is used to denote mean quantities. Horizontal (planar) averaging operation is performed for all the cases. The “prime” symbol is used to represent the fluctuation of a variable with respect to its planar averaged value.

In a recent paper, Basu et al. [10] utilized this DNS database to derive parameterizations for \(\overline{\varepsilon }\). In the present work, the focus is placed on \(\overline{\chi }\).

3 Integral length scales

From the DNS-generated data, we first calculate two different integral length scales as follows:

$$\begin{aligned} {\mathcal {L}}\equiv & {} \frac{\overline{e}^{3/2}}{\overline{\varepsilon }}, \end{aligned}$$
(6a)
$$\begin{aligned} {\mathcal {L}}_{\theta }\equiv & {} \frac{\overline{e}^{1/2}\sigma _\theta ^2}{\overline{\chi }}, \end{aligned}$$
(6b)

where \(\overline{e}\) and \(\sigma _\theta ^2\) denote turbulent kinetic energy (TKE) and the variance of temperature, respectively.

Based on the original ideas of Taylor [58], both [49, 59] provided a heuristic derivation of \({\mathcal {L}}\). Given TKE (\(\overline{e}\)) and mean energy dissipation rate (\(\overline{\varepsilon }\)), an associated integral time scale can be approximated as \(\overline{e}/\overline{\varepsilon }\). One can further assume \(\sqrt{\overline{e}}\) to be the corresponding velocity scale. Thus, an integral length scale can be approximated as \(\overline{e}^{3/2}/\overline{\varepsilon }\). Using dimensional arguments, an analogous length scale \({\mathcal {L}}_{\theta }\) can be formulated based on temperature fluctuations [1, 68].

In the top-panels of Fig. 1, normalized values of \({\mathcal {L}}\) and \({\mathcal {L}}_\theta\) are plotted against the gradient Richardson number (\(Ri_g = N^2/S^2\)); where, N is the Brunt–Väisäla frequency and S is the magnitude of wind shear. In these plots, we have marked four specific points based on the data from DNS run with imposed cooling rate of 10\(^{-3}\) to better understand the effects of height and stability on the integral length scales. The points \(p_1\) and \(p_2\) represent data from \(z/h = 0.1\) and \(z/h = 0.5\), respectively at non-dimensional time (\(T_n\)) of 60. Similarly, \(q_1\) and \(q_2\) are associated with data from \(z/h = 0.1\) and \(z/h = 0.5\), respectively at non-dimensional time (\(T_n\)) of 100.

Fig. 1
figure 1

Top panel: integral length scales as functions of gradient Richardson number. Both the length scales are normalized by the height of the open channel (h). Bottom-left panel: scatter plot of \({\mathcal {L}}\) versus \({\mathcal {L}}_\theta\). Bottom-right panel: normalized \(\overline{\chi }\) as a function of normalized \(\overline{\varepsilon }\), \(\overline{e}\), and \(\sigma _\theta ^2\). Please refer to Eq. 7. Simulated data from five different DNS runs are represented by different colored symbols in these plots. In the legends, CR represents normalized cooling rates

Full size image

Physically, one would expect the integral scales to be increasing with height as long as the eddies feel the presence of the surface (near-neutral or weakly stable condition). For very stable conditions, the eddies no longer feel the presence of the surface. In the atmospheric boundary layer literature, such a situation is known as the z-less condition [28, 64]. Under the influence of strong stability, the integral length scales become more-or-less independent of the height above the surface.

From Fig. 1, it is clear that the integral length scales increase with height and they slowly decrease with time in all the simulations due to the increasing stability effects. Simulations with higher cooling rates have smaller integral length scales. Some of these runs (e.g., \(CR = 5\times 10^{-3}\)) exhibit z-less behavior due to strong stability effects.

Given the similar trends of normalized \({\mathcal {L}}\) and \({\mathcal {L}}_\theta\), they are plotted against each other in the bottom-left panel of Fig. 1. There is (approximately) a linear relationship between these length scales. If \({\mathcal {L}} \propto {\mathcal {L}}_{\theta }\), it is straightforward to derive from Eqs. 6:

$$\begin{aligned} \frac{\overline{\chi }}{\overline{\varepsilon }} \propto \frac{\sigma _\theta ^2}{\overline{e}}. \end{aligned}$$
(7)

This relationship was first reported by Béguier et al. [11] for shear flow turbulence. In a follow-up study, Elghobashi and Launder [22] hypothesized that the similarity of the generation processes of TKE and scalar variance is at the root of this intriguing relationship. In contrast to shear flows, they did not find Eq. 7 to hold for thermal mixed layer.

In the bottom-right panel of Fig. 1, we demonstrate the approximate validity of Eq. 7. Linear least-square regression with bootstrapping [21, 42] is used to estimate the slope of the fitted line. Given that the collapse of the data points is quite reasonable, the relationship \(\overline{\chi } = 1.74 \overline{\varepsilon } \frac{\sigma _\theta ^2}{\overline{e}}\) might be useful for practical applications.

Please note that the appearance of the Prandtl number (Pr) in Fig. 1 (bottom-right panel) is due to the normalization of variables in DNS; “Appendix 2” provides further details. Throughout the paper, the subscript “n” is used to denote a normalized variable.

4 Outer length scales

Both shear and buoyancy prefer to deform larger eddies compared to smaller ones [16, 34, 39, 53]. Turbulent eddies are not affected by shear and buoyancy if they are smaller than the outer length scales (OLSs). Ozmidov (\(L_{OZ}\)) and Corrsin (\(L_C\)) length scales are the most commonly used OLSs in the literature. They are defined as [17, 20, 47]:

$$\begin{aligned} L_{OZ}\equiv & {} \left( \frac{\overline{\varepsilon }}{N^3}\right) ^{1/2}, \end{aligned}$$
(8a)
$$\begin{aligned} L_C\equiv & {} \left( \frac{\overline{\varepsilon }}{S^3}\right) ^{1/2}. \end{aligned}$$
(8b)

Eddies which are smaller than \(L_{OZ}\) are not affected by buoyancy; similarly, shear does not influence the eddies of size less than \(L_{C}\). In other words, the eddies can be assumed to be isotropic if they are smaller than both \(L_{OZ}\) and \(L_{C}\).

Fig. 2
figure 2

Ozmidov (left panel) and Corrsin (right panel) length scales as functions of gradient Richardson numbers. These length scales are normalized by the integral length scale (\({\mathcal {L}}\)). Simulated data from five different DNS runs are represented by different colored symbols in these plots. In the legends, CR represents normalized cooling rates

Full size image

Since \({\mathcal {L}}\) changes across the simulations, the OLS values are normalized by corresponding \({\mathcal {L}}\) values and plotted as functions of \(Ri_g\) in Fig. 2. The collapse of the data from different runs, on to seemingly universal curves, is remarkable for all the cases except for \(Ri_g > 0.2\). We would like to mention that similar scaling behavior was not found if other normalization factors (e.g., h) are used.

Normalized \(L_{OZ}\) decreases monotonically with \(Ri_g\). In contrast, normalized \(L_C\) barely exhibits any sensitivity to \(Ri_g\) (except for \(Ri_g > 0.1\)). Even for weakly-stable condition, it is less than 20% of \({\mathcal {L}}\). Based on the expressions of \(L_{OZ}\), \(L_C\) and \(Ri_g\), we can write:

$$\begin{aligned} \frac{L_C}{L_{OZ}} = \left( \frac{N}{S}\right) ^{3/2} = Ri_g^{3/4}. \end{aligned}$$
(9)

Thus, for \(Ri_g < 1\), one expects \(L_C < L_{OZ}\); this relationship is fully supported by Fig. 2. In comparison to the buoyancy effects, the shear effects are felt at smaller length scales for the entire stability range considered in the present study.

Dissipation rate of turbulent kinetic energy is used in the definitions for both \(L_{OZ}\) and \(L_{C}\). However, it is also possible to formulate OLSs based on the dissipation rate of temperature fluctuations as follows:

$$\begin{aligned} L_1\equiv & {} \left( \frac{g}{\varTheta _0}\right) ^{-1/4} \overline{\chi }^{1/2} \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^{-5/4}, \end{aligned}$$
(10a)
$$\begin{aligned} L_2\equiv & {} \overline{\varepsilon }^{-1/4} \overline{\chi }^{3/4} \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^{-3/2}, \end{aligned}$$
(10b)
$$\begin{aligned} L_3\equiv & {} \overline{\chi }^{1/2} \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^{-1} S^{-1/2}, \end{aligned}$$
(10c)
$$\begin{aligned} L_4\equiv & {} \left( \frac{g}{\varTheta _0}\right) \overline{\chi }^{1/2} S^{-5/2}, \end{aligned}$$
(10d)

where \((\partial {\overline{\theta }}/\partial z)\) is the vertical gradient of mean potential temperature and \(\varTheta _0\) is a reference potential temperature. These length scales were proposed by Panchev based on dimensional analysis [41, 48]. Characteristics of yet another OLS proposed by Bolgiano [13, 14] and Obukhov [46] is discussed separately in “Appendix 1”.

Fig. 3
figure 3

\(\chi\)-based length scales as functions of gradient Richardson numbers. These length scales are normalized by the integral length scale. Simulated data from five different DNS runs are represented by different colored symbols in these plots. In the legends, CR represents normalized cooling rates

Full size image

In Fig. 3, the \(\overline{\chi }\)-based length scale formulations are plotted against \(Ri_g\). Similar to \(L_{OZ}\), the normalized \(L_1\) monotonically decrease with increasing \(Ri_g\). Whereas, the normalized \(L_4\) increase with \(Ri_g\) in an unphysical manner. It is quite evident that both the normalized \(L_2\) and \(L_3\) scales behave very similar to \(L_C\) (see right panel of Fig. 2).

Fig. 4
figure 4

Variation of the normalized \(\overline{\chi }\)-based length scales against the normalized Ozmidov length scale (top-left panel) and the normalized Corrsin length scale (top-right and bottom-left panels). Bottom-right panel: normalized \(\overline{\chi }\) as a function of normalized \(\overline{\varepsilon }\), S, and \((\partial \overline{\theta }/\partial z)\). Please refer to Eq. 11. Simulated data from five different DNS runs are represented by different colored symbols in these plots. In the legends, CR represents normalized cooling rates

Full size image

Given the trends in Fig. 3, we plotted a few inter-relationships of OLSs in Fig. 4. In each case, the collapse of DNS-based data on a single curve is excellent. In the case of \(L_1\)-versus-\(L_{OZ}\) plot, the curve is nonlinear. However, in the case of \(L_2\)-versus-\(L_{C}\) and \(L_3\)-versus-\(L_{C}\) plots, the data fall on more-or-less straight lines. The regressed slopes are reported in the legends of these plots.

If we assume \(L_2 \equiv L_C\), based on Eq. 8b and Eq. 10b, it is trivial to arrive at:

$$\begin{aligned} \frac{\overline{\chi }}{\overline{\varepsilon }} \approx \frac{\left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2}{S^2}. \end{aligned}$$
(11)

Interestingly, the assumption of \(L_3 \equiv L_C\) also leads to the same equation. As a matter of fact, this equation can be derived from the budget equations of TKE and temperature variance with certain assumptions as elaborated below. Assuming steady-state condition, horizontal homogeneity, and neglecting the secondary terms (e.g., transport), we can write:

$$\begin{aligned} \overline{\varepsilon } = -\overline{u'w'}\left( \frac{\partial \overline{u}}{\partial z}\right) -\overline{v'w'}\left( \frac{\partial \overline{v}}{\partial z}\right) + \left( \frac{g}{\varTheta _0}\right) \overline{w'\theta '}, \end{aligned}$$
(12a)

and

$$\begin{aligned} \overline{\chi } = -\,2 \overline{w'\theta '}\left( \frac{\partial \overline{\theta }}{\partial z}\right) . \end{aligned}$$
(12b)

If we apply K-theory, these equations can be further simplified to:

$$\begin{aligned} \overline{\varepsilon } = K_M S^2 - K_H \left( \frac{g}{\varTheta _0}\right) \left( \frac{\partial \overline{\theta }}{\partial z}\right) , \end{aligned}$$
(13a)

and

$$\begin{aligned} \overline{\chi } = 2 K_H \left( \frac{\partial \overline{\theta }}{\partial z}\right) ^2, \end{aligned}$$
(13b)

where \(K_M\) and \(K_H\) are eddy viscosity and diffusivity, respectively. By utilizing the definitions of \(Ri_g\) and turbulent Prandtl number (\(Pr_T = K_M/K_H\)), we can deduce from Eq. 13a and Eq. 13b:

$$\begin{aligned} \frac{\overline{\chi }}{\overline{\varepsilon }} =\frac{2}{\left( Pr_T - Ri_g\right) } \frac{\left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2}{S^2}. \end{aligned}$$
(14)

Anderson [4] conducted a rigorous statistical analysis of the observational data collected at the British Antarctic Survey’s Halley station on the the Antarctic. He avoided the self-correlation issue and proposed the following empirical relationship for \(0.01< Ri_g < 0.25\):

$$ Pr_T^{-1} = (0.84 \pm 0.03) Ri_g^{-0.105 \pm 0.012}.$$
(15)

Clearly, the \(Ri_g\)-dependence of the Prandtl number is rather weak for small values of \(Ri_g\). Similar findings were reported in other experimental and modeling studies (e.g., [35, 38, 55]).

In the bottom-right panel of Fig. 4, we have plotted Eq. 11 in a normalized form. The slope of the fitted line is 2.35. For \(Ri_g = 0.2\), according to Eq. 15, \(Pr_T \approx 1\). Thus, the ratio \(2/(Pr_T - Ri_g)\) is approximately 2.48. When \(Ri_g\) equals to 0.1, \(Pr_T \approx 0.93\) following Eq. 15. In this case, the ratio \(2/(Pr_T - Ri_g)\) is close to 2.40. These values are not far from the estimated slope of 2.35 in Fig. 4 (bottom-right panel). In other words, our DNS-based results are in-line with past observational studies.

5 Structure parameter of temperature (\(C_T^2\))

Using the DNS database of the current study, Basu et al. [10] recently found that \(\overline{\varepsilon } = 0.23 \overline{e} S\) and \(\overline{\varepsilon } = 0.63 \sigma _w^2 S\) for \(0< Ri_g < 0.2\). If we insert these formulations in Eq. 11, we get:

$$\begin{aligned} \overline{\chi } \approx \left( \frac{\overline{e}}{S}\right) \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2, \end{aligned}$$
(16a)

and

$$\begin{aligned} \overline{\chi } \approx \left( \frac{\sigma _w^2}{S}\right) \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2. \end{aligned}$$
(16b)
Fig. 5
figure 5

Top panels: normalized \(\overline{\chi }\) as a function of normalized \(\overline{e}\), \(\sigma _w^2\), S, and \((\partial \overline{\theta }/\partial z)\). Please refer to Eqs. 16. Bottom panels: normalized \(\overline{\varepsilon }^{-1/3}\overline{\chi }\) as a function of normalized \(\overline{e}\), \(\sigma _w^2\), S, and \((\partial \overline{\theta }/\partial z)\). Please refer to Eqs. 17. Simulated data from five different DNS runs are represented by different colored symbols in these plots. In the legends, CR represents normalized cooling rates

Full size image

The top panels of Fig. 5 strongly support the validity of these formulations. The proportionality constants in these equations are found to be equal to 0.55 and 1.47, respectively.

By definition, \(C_T^2 \approx \overline{\varepsilon }^{-1/3} \overline{\chi }\). The proportionality constant is usually taken equal to 1.6 [31, 67]. Thus, we can write:

$$\begin{aligned} C_T^2 \approx \overline{\varepsilon }^{-1/3} \overline{\chi } \approx \left( \frac{\overline{e}^{1/2}}{S}\right) ^{4/3} \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2, \end{aligned}$$
(17a)

and

$$\begin{aligned} C_T^2 \approx \overline{\varepsilon }^{-1/3} \overline{\chi } \approx \left( \frac{\sigma _w}{S}\right) ^{4/3} \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2. \end{aligned}$$
(17b)

In the bottom panels of Fig. 5, we establish that these equations (especially Eq. 17b) nicely hold for our DNS-generated data.

Based on theoretical and numerical work, Hunt et al. [32, 33] proposed the shear-based length scales, \(L_H \equiv \left( \frac{\overline{e}^{1/2}}{S}\right)\) and \(L_H \equiv \left( \frac{\sigma _w}{S}\right)\), as the characteristic length scales for \(0< Ri_g < 0.5\). Thus, we can re-write Eqs. 17 as:

$$\begin{aligned} C_T^2 \approx L_H^{4/3} \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2. \end{aligned}$$
(18)

A very similar equation was proposed by Tatarskii more than 50 years ago [56, 57], albeit with an OLS which needs to be prescribed. In the literature, several empirical parameterizations were proposed for this unknown length scale [8, 18, 19, 62]. In this study, based on DNS-generated data, we demonstrate that the outer length scale in Tatarskii’s equation should be equal to \(L_H\) for \(0< Ri_g < 0.2\).

At this point, we point out an interesting relationship that one can further deduce from our findings. If we compare Eq. 7 against Eq. 11, we get:

$$\begin{aligned} \frac{\sigma _\theta ^2}{\overline{e}} \approx \frac{\left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2}{S^2}. \end{aligned}$$
(19)

Equivalently, one can write:

$$\begin{aligned} \frac{\sigma _\theta }{\left( \frac{\partial \overline{\theta }}{\partial z} \right) } \approx \frac{\overline{e}^{1/2}}{S}, \end{aligned}$$
(20a)

or,

$$\begin{aligned} L_E \equiv L_H. \end{aligned}$$
(20b)

where \(L_E\) is a length scale proposed by Ellison [23]. The dependence of \(L_E\) on \(Ri_g\) is documented in the left panel of Fig. 6. In the right panel, we show the one-to-one relationship between \(L_E\) and \(L_H\). With the exception of a few data points from the simulation with the strongest cooling rate, it is clear that these length scales are linearly related to each other. Thus, the following equation can be used as a viable alternative to Eq. 18:

$$\begin{aligned} C_T^2 \approx L_E^{4/3} \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2. \end{aligned}$$
(21)
Fig. 6
figure 6

Left panel: normalized Ellison length scale as a function of gradient Richardson number. Right panel: normalized Ellison length scale as a function of normalized Hunt length scale. Simulated data from five different DNS runs are represented by different colored symbols in these plots. In the legends, CR represents normalized cooling rates

Full size image

In the literature, several studies have demonstrated the similarities between the so-called Thorpe scale (\(L_T\); [60, 61]) and \(L_E\) using observed and simulated data (e.g., [34, 39]). A simple heuristic derivation was also provided by Gavrilov et al. [27]. Thus, it is plausible to replace \(L_E\) with \(L_T\) in Eq. 21:

$$\begin{aligned} C_T^2 \approx L_T^{4/3} \left( \frac{\partial \overline{\theta }}{\partial z} \right) ^2. \end{aligned}$$
(22)

This equation was proposed by Basu [8] and was validated using observational data from a field campaign over Mauna Kea, Hawaii.

In summary, we conjecture that Eqs. 18,  21, and  22 are all valid parameterizations for \(C_T^2\) as long as \(Ri_g\) does not exceed 0.2. For larger values of \(Ri_g\), a different length scale might be more appropriate; our present DNS runs cannot shed light on such a strong stability regime.

6 Concluding remarks

In this study, we analyze DNS-generated data to characterize several integral and outer length scales. From these results, we propose simple parameterizations for \(\chi\) and \(C_T^2\) when gradient Richardson number is less than 0.2. In the continuously turbulent atmospheric stable boundary layer, \(Ri_g\) is usually less than 0.2 [26, 45]. Thus, the proposed parameterizations should be suitable for certain practical boundary layer problems. However, they will have limited applications for intermittently stable conditions.

In closing, we would like to emphasize the importance of Eq. 19. To the best of our knowledge, it was first reported by Fulachier and Dumas [25] from boundary layer experiments over a slightly heated plate. In a latter study, Fulachier and Antonia [24] found this formulation to hold for various other types of flows. They even concluded:

It seems therefore reasonable, from both mathematical and physical points of view, to seek a relationship, not between momentum and heat fluxes, as in the case with the Reynolds analogy, but preferably between the turbulent kinetic energy and the temperature variance.

To the best of our knowledge, Eq. 19 is not used in atmospheric boundary layer studies. Since our findings are in agreement, we strongly endorse the assertion of Fulachier and Antonia [24] and advocate further research on this equation.