On the Dissipation Rate of Temperature Fluctuations in Stably Stratified Flows

In this study, we explore several integral and outer length scales of turbulence which can be formulated by using the dissipation of temperature fluctuations ($\chi$) and other relevant variables. Our analyses directly lead to simple yet non-trivial parameterizations for both $\chi$ and the structure parameter of temperature ($C_T^2$). For our purposes, we make use of high-fidelity data from direct numerical simulations of stratified channel flows.


I. INTRODUCTION
The molecular dissipation of temperature fluctuations (χ) is an important variable for characterizing turbulent mixing in various environmental flows. It is frequently used in micrometeorology [e.g., 55] and atmospheric optics [e.g., 35]. Furthermore, any higher-order closure model requires solving a prognostic equation or a diagnostic parameterization for χ [refer to 13,32,54].
Over the years, a number of studies focused on the correlation between turbulent kinetic energy dissipation rate (ε) and χ [e.g., 2, 4-6, 25, 58]. In addition, some papers reported on the probability density function, spatiotemporal intermittency and anomalous scaling of χ [e.g., 5,43,45]. Often, χ has been found to be more intermittent and non-Gaussian than ε.
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 χ and other relevant variables. Based on these findings, we then derive simple parameterizations for χ and temperature structure parameter (C 2 T ). For all the analyses, we utilize a direct numerical simulation (DNS) database of stably stratified flows which is discussed in the following section.

II. DIRECT NUMERICAL SIMULATION
Recently, for the parameterization of optical turbulence, He and Basu [27] created a DNS database using a massively parallel DNS code, called HERCULES [26]. The computational domain size for all the DNS runs was L x × L y × L z = 18h × 10h × h, where h is the height of the open channel. The domain was discretized by 2304×2048×288 grid points in streamwise, spanwise, and wall-normal directions, respectively. The bulk Reynolds number, Re b = U b h ν , for all the simulations was fixed at 20000; where, U b and ν denote the bulk (averaged) velocity in the channel and kinematic viscosity, respectively. The bulk Richardson number was calculated as: ; where, Θ top and Θ bot represent potential temperature at the top and the bottom of the channel, respectively. The gravitational acceleration is denoted by g. A total of five simulations were performed with gradual decrease in the temperature of the bottom wall (effectively by increasing Ri b ) to mimic the nighttime cooling of the land-surface. The normalized cooling rates (CR), Ri b /T n , ranged from 1 × 10 −3 to 5 × 10 −3 ; where, T n is a non-dimensional time (= tU b /h). Since we were considering stably stratified flows in the atmosphere, the Prandtl number, P r = ν/k was assumed to be equal to 0.7 with k being the thermal diffusivity.
All the simulations used fully developed neutrally stratified flows as initial conditions and evolved for upto 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 ≤ T n ≤ 100). Furthermore, we only consider data from the region 0.1h ≤ z ≤ 0.5h to discard any blocking effect of the surface or avoid any laminarization in the upper part of the open channel.
The mean dissipation of turbulent kinetic energy and temperature fluctuations are computed as follows (using Einstein's summation notation): 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. [8] utilized this DNS database to derive parameterizations for ε. In the present work, the focus is placed on χ.

III. INTEGRAL LENGTH SCALES
From the DNS-generated data, we first calculate two different integral length scales as follows: where, e and σ 2 θ denote turbulent kinetic energy (TKE) and the variance of temperature, respectively.
Based on the original ideas of Taylor [49], both Tennekes and Lumley [50] and Pope [41] provided a heuristic derivation of L. Given TKE (e) and mean energy dissipation rate (ε), an associated integral time scale can be approximated as e/ε. One can further assume √ e to be the corresponding velocity scale. Thus, an integral length scale can be approximated as e 3/2 /ε. Using dimensional arguments, an analogous length scale L θ can be formulated based on temperature fluctuations [1,57].
In the top-panel of Fig. 1, normalized values of L and L θ 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. It is evident that the simulations with larger cooling rates result in smaller integral length scales as would be physically expected.
Given the similar trends of normalized L and L θ , they are plotted against each other in the bottom-left panel of Fig. 1. These length scales are significantly (but not perfectly) correlated. If L ≈ L θ , it is straightforward to derive from Eqs. 2: This relationship was first reported by Béguier et al. [9] for shear flow turbulence. In a follow-up study, Elghobashi and Launder [20] 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. 3 to hold for thermal mixed layer.
In the bottom-right panel of Fig. 1, we demonstrate the approximate validity of Eq. 3. Linear least-square regression with bootstrapping [19,34] is used to estimate the slope of the fitted line. Given that the collapse of the data points is quite reasonable, the relationship χ = 0.87ε σ 2 θ e might be useful for practical applications. Please note that the appearance of the Prandtl number (P r) 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.

IV. OUTER LENGTH SCALES
Both shear and buoyancy prefer to deform larger eddies compared to smaller ones [14,30,31,44]. 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 [15,18,39]: 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 . Since L changes across the simulations, the OLS values are normalized by corresponding 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 L. Based on the expressions of L OZ , L C and Ri g , we can write: 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 θ . Please refer to Eq. 3. Simulated data from five different DNS runs are represented by different colored symbols in these plots. In the legends, CR represents normalized cooling rates.
possible to formulate OLSs based on the dissipation rate of temperature fluctuations as follows: where, (∂θ/∂z) is the vertical gradient of mean potential temperature and Θ 0 is a reference potential temperature. These length scales were proposed by Panchev based on dimensional analysis [33,40]. Characteristics of yet another OLS proposed by Bolgiano [11,12] and Obukhov [38] is discussed separately in Appendix 1.
In Fig. 3, the χ-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).
Given the trends in Fig. 3, we plotted a few interrelationships 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 -vs-L OZ plot, the curve is nonlinear. However, in the case of L 2 -vs-L C and L 3 -vs-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 ≡ L C , based on Eq. 4b and Eq. 6b, it is trivial to arrive at: Interestingly, the assumption of L 3 ≡ 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: and If we apply K-theory, these equations can be further simplified to: and where, K M and K H are eddy viscosity and diffusivity, respectively. By utilizing the definitions of gradient Richardson number (Ri g ) and turbulent Prandtl number (P r T = K M /K H ), we can deduce from Eq. 9a and Eq. 9b: Thus, if P r T = c+ Ri g , where c is a constant, then, Eq. 7 holds. In the bottom-right panel of Fig. 4 Using the DNS database of the current study, Basu et al. [8] recently found that ε = 0.23eS and ε = 0.63σ 2 w S for 0 < Ri g < 0.2. If we insert these formulations in Eq. 7, we get: and 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.28 and 0.74, respectively. By definition, C 2 T ≈ ε −1/3 χ. The proportionality constant is usually taken equal to 1.6 [27,56]. Thus, we can write: and In the bottom panels of Fig. 5, we establish that these equations (especially Eq. 12b) nicely hold for our DNSgenerated data. Based on theoretical and numerical work, Hunt et al. [28,29] proposed the shear-based length scales, L H ≡ e 1/2 S and L H ≡ σw S , as the characteristic length scales for 0 < Ri g < O(0.5). Thus, we can re-write Eqs. 12 as: A very similar equation was proposed by Tatarskii more than 50 years ago [47,48], albeit with an OLS which needs to be prescribed. In the literature, several empirical parameterizations were proposed for this unknown length scale [7,16,17,53]. In this study, based on DNSgenerated data, we demonstrate that the outer length scale in Tatarskii's equation should be equal to L H for 0 < Ri g < O(0.2). At this point, we point out an interesting relationship that one can further deduce from our findings. If we compare Eq. 3 against Eq. 7, we get: Equivalently, one can write: or, where, L E is a length scale proposed by [21]. 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. 13: In the literature, several studies have demonstrated the similarities between the so-called Thorpe scale [L T ; 51,52] and L E using observed and simulated data [e.g., 30,31]. A simple heuristic derivation was also provided by Gavrilov et al. [24]. Thus, it is plausible to replace L E with L T in Eq. 16: This equation was proposed by Basu [7] and was validated using observational data from a field campaign over Mauna Kea, Hawaii. In summary, we conjecture that Eqs. 13, 16, and 17 are all valid parameterizations for C 2 T as long as Ri g does not exceed O(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.

VI. CONCLUDING REMARKS
In this study, we analyze DNS-generated data to characterize several integral and outer length scales. From w , S, and (∂θ/∂z). Please refer to Eqs. 11. Bottom panels: normalized ε −1/3 χ as a function of normalized e, σ 2 w , S, and (∂θ/∂z). Please refer to Eqs. 12. Simulated data from five different DNS runs are represented by different colored symbols in these plots. In the legends, CR represents normalized cooling rates.
these results, we propose simple parameterizations for χ and C 2 T when gradient Richardson number is less than O(0.2). In the atmospheric stable boundary layer, Ri g is usually less than 0.2 [37]. Thus, the proposed parameterizations should be applicable for practical boundary layer problems.
In closing, we would like to emphasize the importance of Eq. 14. To the best of our knowledge, it was first reported by Fulachier and Dumas [23] from boundary layer experiments over a slightly heated plate. In a latter study, Fulachier and Antonia [22] 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. 14 is not used in atmospheric boundary layer studies. Since our findings are in agreement, we strongly endorse the assertion of Fulachier and Antonia [22] and advocate further research on this equation.

DATA AND CODE AVAILABILITY
The DNS code (HERCULES) is available from: https://github.com/friedenhe/HERCULES. Upon acceptance of the manuscript, all the analysis codes and processed data will be made publicly available via zenodo.org. Given the sheer size of the raw DNS dataset, it will not be uploaded on to any repository; however, it will be available upon request from the authors.

ACKNOWLEDGMENTS
The first author thanks Bert Holtslag for having interesting discussions on this topic. The authors acknowledge computational resources obtained from the Department of Defense Supercomputing Resource Center (DSRC) for the direct numerical simulations. The views expressed in this paper do not reflect official policy or position by the U.S Air Force or the U.S. Government.

APPENDIX 1: BOLGIANO-OBUKHOV LENGTH SCALE
Bolgiano [11,12] and Obukhov [38] independently proposed a buoyancy-range scaling and the following OLS based on theoretical arguments: Several laboratory and numerical studies [e.g., 10,36] reported the existence of Bolgiano scaling in unstable condition. However, studies involving stably stratified conditions are rather limited [3,42] and somewhat inconclusive in terms of the existence of this OLS. In Fig. 7, we show the traits of L BO as a function of Ri g . Similar to L OZ and L 1 , this length scale also shows a decreasing trend with increasing stability. However, the relationship between L BO and L OZ is nonlinear. As a consequence, we were unable to derive any simple formulation involving ε, χ, and other variables.

APPENDIX 2: NORMALIZATION OF DNS VARIABLES
In DNS, the relevant variables are normalized as follows: After differentiation, we get: The gradient Richardson number can be expanded as: Using the definition of Ri b (see Sect. 2), we re-write Ri g as follows: Similarly, N 2 can be written as: The velocity variances, TKE, and temperature variance can be normalized as: Following the above normalization approach, we can also derive the following relationships for the dissipation rate of TKE and variance of temperature fluctuations: We can combine Eqs. 24d, 24e, 25a, and 25b, we can re-write Eq. 3 as follows: χ n = ν k σ 2 θn e n ε n = P r σ 2 θn e n ε n .
In a similar fashion, we can utilize Eqs. 20c, 20d, 24d, and 25b to re-write Eq. 11a as follows: