Turbulent Prandtl number and characteristic length scales in stably stratified flows: steady-state analytical solutions

In this study, the stability dependence of turbulent Prandtl number ($Pr_t$) is quantified via a novel and simple analytical approach. Based on the variance and flux budget equations, a hybrid length scale formulation is first proposed and its functional relationships to well-known length scales are established. Next, the ratios of these length scales are utilized to derive an explicit relationship between $Pr_t$ and gradient Richardson number. In addition, theoretical predictions are made for several key turbulence variables (e.g., dissipation rates, normalized fluxes). The results from our proposed approach are compared against other competing formulations as well as published datasets. Overall, the agreement between the different approaches is rather good despite their different theoretical foundations and assumptions.


I. INTRODUCTION
According to the K-theory, based on the celebrated hypothesis of Boussinesq in 1877, turbulent fluxes can be approximated as products of the eddy exchange coefficients (known as the Austausch coefficients in earlier literature) and the mean gradients [45]. Specifically, for incompressible, horizontally homogeneous, boundary layer flows, the along-wind momentum flux (u ′ w ′ ) and the sensible heat flux (w ′ θ ′ ) can be simply written as follows: Here S and Γ denote the vertical gradients of the mean along-wind velocity component and the mean potential temperature, respectively. The eddy viscosity and diffusivity for heat are represented by K M and K H , respectively. In contrast to molecular diffusivities, these eddy exchange coefficients are not intrinsic properties of the fluid [4,64]; rather, they depend on the nature of the turbulent flows (e.g., stability) and position in the flow (e.g., distance from the wall). The ratio of K M and K H is known as the turbulent Prandtl number: This variable is fundamentally different from the molecular Prandtl number: where, ν and α denote kinematic viscosity and thermal diffusivity, respectively. According to a vast amount of literature, P r t is strongly dependent on buoyancy and somewhat weakly dependent on other factors (see below).
For non-buoyant (also called neutral) flows, in this paper, the turbulent Prandtl number is denoted as P r t0 . In the past, for simplicity, a number of studies assumed P r t0 = 1 by invoking the so-called 'Reynolds analogy' hypothesis [55,68,69]. Basically, they implicitly assume that the turbulent transport of momentum and heat are identical. However, this assumption of P r t0 = 1 is not supported by the vast majority of experimental data (see [33] and the references therein). On this issue, Launder [39] commented: "It would also be helpful to dispel the idea that a turbulent Prandtl number of unity was in any sense the "normal" value. We shall see [...] that a value of about 0.7 has a far stronger claim to normality." Perhaps, it is not a mere coincidence that the theoretical study of Yakhot et al. [78] predicted that P r t0 asymptotically approaches 0.7179 in the limit of infinite Re (see also [66]). One of the most cited studies in atmospheric science, by Businger et al. [11], also reported P r t0 = 0.74. According to a review article by Kays [33], for laboratory flows, P r t0 typically falls within the range of 0.7 to 0.9; the most frequent value being equal to 0.85. Most commercial computational fluid dynamics packages (e.g., Fluent, OpenFOAM) assume 0.85 to be the default P r t0 value.
There is some evidence that P r t0 may not be a universal constant; it might weakly depend on P r m , Re, and/or position in the flow. However, there is no general agreement in the literature on this matter (e.g., [3]). Reynolds [56] summarized numerous empirical and semiempirical formulations capturing such dependencies for a wide range of fluids (including air, water, liquid metal) and engineering flows (e.g., pipe flow, jet flow, shear flow). However, to the best of our knowledge, these formulations are yet to be confirmed for high-Re atmospheric flows. In such flows, buoyancy effects have been found to be far more dominant than any other factors.
In atmospheric flows, especially under stably stratified conditions, the value of P r t departs significantly from P r t0 . Over the decades, several empirical formulations have been developed by various research groups (see [40] for a recent review). For example, by regression analysis of aircraft measurements from different field campaigns, Kim and Mahrt [34] proposed: where, Ri g is the gradient Richardson number, commonly used to quantify atmospheric stability. It is defined as follows: Where, g is the gravitational acceleration and Θ 0 represents a reference temperature. The variable β is known as the buoyancy parameter. The so-called Brunt Väisälä frequency is denoted by N . More recently, Anderson [1] conducted rigorous statistical analysis of observational data from the Antarctic. By avoiding the self-correlation issue, he proposed the following empirical relationship for 0.01 < Ri g < 0.25: Clearly, the Ri g -dependence of P r t becomes rather weak as the stability of the flow decreases. In addition to field observational data, laboratory and simulated data were also utilized to quantify the P r t -Ri g relationship. In this regard, a popular semi-empirical formulation by Schumann and Gerz [58] is worth noting: where, R f ∞ is the asymptotic value of the flux Richardson number (R f = Ri g /P r t ) for strongly stratified conditions. Recently, Venayagamoorthy and Stretch [72] used direct numerical simulation (DNS) data and revised the formulation by Eq. (7) as follows: For all practical purposes, the differences between Eq. (7) and Eq. (8) are quite small. In parallel to observational and simulation studies, there have been a handful of attempts to derive the P r t -Ri g formulations from the governing equations with certain assumptions. In the appendices, we have summarized two competing hypotheses by Katul et al. [32] and Zilitinkevich et al. [82]. The readers are also encouraged to peruse the following papers describing other relevant hypotheses: [14], [15], and [29]. In the present study, we report an alternative analytical derivation which leads to a closed-form P r t -Ri g relationship.

II. ANALYTICAL DERIVATIONS
In this section, based on the variance and flux budget equations, we first derive a hybrid length scale (L X ) and establish its relationship with three well-known length scales: the Hunt length scale (L H , [23,24]), the buoyancy length scale (L b , [10,76]), and the Ellison length scale (L E , [16]). Next, the ratios of various length scales (e.g., L b /L E ) are shown to be explicit functions of Ri g and P r t . Equating these functions with one another results in a quadratic equation for P r t . One of the roots of this quadratic equation provides an explicit P r t -Ri g relationship.

A. Budget Equations
The simplified budget equations for turbulent kinetic energy (TKE), variance of temperature (σ 2 θ ), and sensible heat flux (w ′ θ ′ ) can be written as [18,53,75]: where ε and χ θ denote the dissipation rates of TKE and σ 2 θ , respectively. The variance of vertical velocity is σ 2 w . In Eq. (9c), the parameter a p influences the buoyant contribution to the pressure-temperature interaction term; whereas, the last term of this equation is a parameterization of the turbulent-turbulent component of the pressure-temperature interaction. The return-toisotropy time scale is denoted by τ R . Please refer to Appendix 1 for further technical details on the parameterization of pressure-temperature interaction.

B. A Hybrid Length Scale
In analogy to Prandtl's mixing length hypothesis (see [4,50,73]), let us assume that σ w is a characteristic velocity scale for stably stratified flows. Further assume that L X and L X /σ w are characteristic length and time scales, respectively. Then, the eddy diffusivity, the dissipation rates, and turbulent-turbulent component of the pressure-temperature interaction can be re-written as follows: Here the unknown (non-dimensional) coefficients are denoted as c i , where i is an integer. The parameterizations for the dissipation rates (i.e., ε and χ θ ) are further discussed in Section IV.
The Hunt length scale is related to the so-called buoyancy length scale (L b ) as follows: (13) Thus, Eq. (12b) can be re-written as: If we substitute the individual terms of Eq. (9b) by utilizing Eqs. (1b), (2), (10a), and (10c), we get: Simplification of this equation leads to: where L E (= σ θ /Γ) is the Ellison length scale and c E is an unknown (nondimensional) coefficient. We would like to point out that in the appendices of Basu et al. [6,7] we have summarized the characteristics of Hunt, buoyancy, Ellison, Bolgiano, Ozmidov, and several other length scales. For brevity, we do not repeat them here.

C. Ratios of Length Scales
By comparing Eq. (12b) with Eq. (16b), it is rather straightforward to derive: or, where . Using Eq. (13), this equation can be re-written as follows: An alternative expression for can be found if we use Eqs. (1b), (2), (10a), and (10d) to substitute the individual terms of Eq. (9c) as follows: or, where c 5 (= c 1 c 4 ) is an unknown proportionality constant.

D. Derivation of Prandtl Number
By equating Eq. (18) and Eq. (19d), we immediately get the following quadratic equation: Since P r t = P r t0 for neutral conditions (Ri g = 0), via Eq. (20), we find: The roots of Eq. (20) are: where, X = [P r t0 + Ri g + (1 − a p ) c P Ri g ]. Only the larger root is physically meaningful. Eq. (22) includes three unknown parameters (i.e., P r t0 , a p , and c P ). Similarity theory can be used to estimate c P (discussed in the following section). However, P r t0 and a p must be prescribed.
We would like to emphasize that Eq. (22) is a closed form analytical solution for the stability-dependence of P r t . It is derived directly from the budget equations without any additional simplification. Since our derivation makes use of certain length scale ratios (LSRs), we refer to our proposed approach as the LSR formulation.

III. ESTIMATION OF UNKNOWN COEFFICIENTS
For near-neutral conditions, Eqs. (12b) and (16b) simplify to the following expressions, respectively: In order to be consistent with the logarithmic velocity profile in the surface layer, L X should be equal to κz in the surface layer, where κ is the von Kármán constant. Therefore, Numerous studies reported that σ w = c w u * and σ θ = c θ θ * in near-neutral stratified surface layer. The surface friction velocity and temperature scale are denoted by u * and θ * , respectively. Thus, we get: Please note that the non-dimensional velocity gradient, (κzS/u * ), equals to unity according to the logarithmic law of the wall. Whereas, the non-dimensional temperature gradient, (κzΓ/θ * ), equals to P r t0 . By using Eqs. (1a), (10a), and (12b), we can expand the along-wind momentum flux as follows: Thus, the normalized momentum flux can be written as: For neutral condition, R uw simplifies to: Since, σ w = c w u * , we get: Since, c H ≈ 1 cw , the unknown coefficient c 1 is also approximately equal to 1 cw . Typical values of R uw0 are documented in Table I.
In the literature, the most commonly reported values of c w range from 1.25-1.30 [4,27,53,60]. Similarly, c θ values vary approximately from 1.8 to 2.0 [27,60]. In a few publications, somewhat different values were also reported (e.g., [45,74]). In Table I, we have computed c i and other coefficients for a few combinations of P r t0 , c w , and c θ .

A. Energy Dissipation Rate
The energy dissipation rate is commonly parameterized as follows [46]: where q 2 is twice TKE. L M is known as the master length scale and B 1 is a constant coefficient. In this study, following Townsend [70], we use Eq. (10b) as an alternative parameterization for ε which makes use of σ 3 w instead of q 3 . Using Eqs. (12b), and (25), we can re-write this parameterization as follows: If the value of c w is approximately in the range of 1.25-1.30 (refer to Table T1), for small values of Ri g (i.e., weakly stable conditions), we get: It is important to note that Eq. (27b) (with an unknown proportionality constant) was originally proposed by Hunt [24] using heuristic arguments. He hypothesized that the energy dissipation in weakly/moderately stably stratified flows is dictated by mean shear (S) and root-mean-square value of vertical velocity fluctuations (i.e., σ w ) which is the characteristic velocity scale in the direction of S. Later on Schumann and Gerz [58] analyzed various observational and simulation datasets and validated Hunt's parameterization (see their Figure 1). More recently, Basu et al. [7] utilized a database of direct numerical simulations and found: for 0 < Ri g < 0.2. TKE is denoted by e. It is remarkable that the DNS-based empirical formulation of [7] is virtually identical to our analytical prediction, i.e., Eq. (27b). However, we are unable to ascertain the validity of either Eq. (27b) or Eq. (27c) for Ri g > 0.2. We will discuss more on this issue in Section VI. The exact value of B 1 in Eq. (26) is not settled in the literature. Over the years, a number of researchers estimated its value from diverse observational and simulated datasets; see a brief summary in Table II. By combining the analytical results from the present study with the DNS results from Basu et al. [7], we can also estimate B 1 as follows. From Eq. (27c), for 0 < Ri g < 0.2, we can write: 24.0 Janjić [25] 11.9 Cheng et al. [14] 19.3 Basu et al. [7] 25.8 Next, if we assume our proposed length scale (L X ) is equal to the master length scale (L M ), then from Eqs. (12b) and (26), we get: Here we have assumed c H = 0.8 and 1 − Ri g /P r t ≈ 1 for small values of Ri g . Clearly, our estimated value of B 1 agrees reasonably well with some of the published studies; however, it is significantly higher than the widely used value of 16.6. Please note that due to a missing multiplying coefficient of value 2.1, Basu et al. [7] incorrectly reported B 1 = 12.3 instead of 25.8.

B. Dissipation Rate of Temperature Variance
Once again, following Townsend [70], we parameterized the dissipation rate of temperature variance (χ θ ) by Eq. (10c). Combining this equation with Eq. (12b), Eq. (15), and Eqs. (25), we get: For small values of Ri g , we can assume P r t ≈ 0.85. As before, if we also consider c H = 0.8, we arrive at: Almost the same formulation was reported by Basu et al. [6] based on their analysis of a DNS database. For 0 < Ri g < 0.2, they found: In summary of this section, we can state that our analytical formulations of dissipation rates are very reliable for 0 < Ri g < 0.2. However, more research will be needed for their rigorous validation for the very stable regime (i.e., Ri g > 0.2).

A. Turbulent Prandtl Number
Our proposed formulation for the turbulent Prandtl number, Eq. (22), contains 3 unknown coefficients: P r t0 , a p , and c P . Based on the discussion in the Introduction, in this study, we have opted to use P r t0 = 0.85. The value of c P is selected from Table I; it is evident that it should vary within a range of 2.4-4.3 for typical values of c w and c θ . The parameter a p is discussed in Appendix 1.
In Fig. 1, the predictions from our LSR approach are reported for various combinations of a p and c P . In addition to P r t , we have also reported the stabilitydependence of R f . The results are sensitive to a p values for Ri g > 0.1. It is encouraging to see that the predictions are qualitatively in agreement with the published observations. They are also in-line with the predictions from the co-spectral budget (CSB; [32]) and energy-and flux-budget (EFB; [82]) approaches.
We would like to emphasize out that Eq. (22) and Eq. (72b) in Appendix 2 have nearly identical mathematical form despite the fundamental differences in the LSR and CSB approaches. The CSB approach includes prescribed coefficients from Kolmogorov-Obukhov-Corrsin hypotheses and from a parameterization of the pressuretemperature decorrelation (refer to Appendix 2); they are all lumped into a variable called ω CSB in Eq. (72b). However, it does not consider the buoyancy-turbulence interaction term in the sensible heat flux equation. Thus, Eq. (72b) does not include the a p parameter. In contrast, the LSR approach largely depends on c w and c θ coefficients (combined into the c P coefficient) in addition to a p . These coefficients are integral part of surface layer similarity theory for near-neutral conditions. Furthermore, by construction, the CSB approach assumes P r t0 = 1. Whereas, in the case of the LSR approach, P r t0 is assumed to be equal to 0.85.
For very stable condition (i.e., Ri g ≫ 1), Eq. (22) is simplified to: In contrast, Eq. (72b) from the CSB approach leads to: Thus, the CSB approach predicts R f ∞ ≈ 0.25. On the other hand, for a p = 0 and c P = 4.27, R f ∞ equals to 0.19 for the LSR approach. However, for a p = 0.5 and c P = 2.4, R f ∞ increases to 0.46. In the literature (see [16], [20], [70], [79]), R f ∞ has been reported to be within the limits of 0.15 and 0.5; both the LSR-based and CSBbased predictions are in this range.

B. Normalized Variances and Fluxes
In the literature, there is no consensus regarding the exact stability-dependence of a few normalized variables. Different formulations (e.g., [43], [82]) predict different trends. The LSR approach allows us to independently predict some of these ratios without further approximations as elaborated below.

Ratio of Turbulent Potential and Kinetic Energies
We first consider the ratio of the turbulent potential energy (TPE; denoted as e p ) and the vertical component of TKE (i.e., e w ). These variables are commonly written as [43]: where e T = σ 2 θ 2 . By using the definition of the Ellison length scale (L E ), we can re-write e p as follows: Thus, the ratio of e p and e w is simply: By making use of Eq. (18), we can re-write R pw as follows:  [82], and Katul et al. [32] are also shown in these panels for comparison.
In the top panel of Fig. 2, the dependence of R pw on Ri g is shown. Clearly, R pw is strongly influenced by a p for Ri g > 0.2. In contrast, somewhat surprisingly, R pw is not very sensitive to the coefficient c P . In the denominator of R pw , the term (P r t − Rig) appears which strongly depends on c P . It effectively cancels out the influence of c P in the numerator of R pw .

Normalized Momentum Flux
The formulations for R uw and R uw0 are derived earlier in Eqs. (24b) and (24c), respectively. Hence, their ratio becomes: The dependence of the normalized momentum flux on Ri g is shown in the middle panel of Fig. 2. It is marginally sensitive to a p and c P .

Normalized Correlation of w and θ
Similar to the momentum flux expression, the sensible heat flux can be re-written using Eqs. (1b), (2), (10a), and (16b) as follows: Hence, the correlation between w and θ becomes: For neutral condition, we have R wθ0 = − c1cE √ P rt0 . So, the normalized correlation can be written as: Typical values of R wθ0 are documented in Table I. The normalized correlations are plotted in the right panel of Fig. 2. Similar to the normalized momentum flux, this ratio is also very weakly dependent on a p and c P .

Comparison of Different Theoretical Approaches
As documented in Appendix 2, the CSB approach of Katul et al. [32] predicts: where c CSB 0 and c CSB T equal to 0.65 and 0.80, respectively. On the other hand, according to the EFB approach of Zilitinkevich et al. [82], we have (refer to Appendix 3): where, c EF B P is 0.86. Zilitinkevich et al. [82] assumed that the anisotropy parameter A z (discussed in the following section) varies from 0.2 (neutral condition) to 0.03 (strongly stratified condition).
We intercompare Eqs. (35b), (38) and (39) via Fig. 3  (left panel). In comparison to the LSR approach, the CSB approach underestimates R pw by a factor of more than 2. The CSB approach makes an assumption that the temperature spectrum has a flat shape in the buoyancy range (refer to Appendix 2) which is not supported by field observations. We speculate that, as a consequence of this idealization, the CSB approach underestimates the variance of temperature, and in turn, underestimates R pw . The predictions of the EFB approach and the LSR approach agree reasonably well up to R f ≈ 0.15. For higher stability conditions, the EFB predicts a sharp increase in R pw values. This drastic behavior can be attributed to the assumed stability-dependence of A z (see Fig. 6 of [82]).
In the context of normalized momentum fluxes, the CSB and LSR approaches make identical predictions; please compare Eqs. (36) and (76). However, the prediction from the EFB approach include terms involving A z in the numerator [refer to Eq. (85b)]. Thus, owing to the assumed stability-dependence of A z , the EFB approach predicts much higher value of normalized momentum fluxes in comparison to the LSR approach as depicted in the right panel of Fig. 3. Rigorous analyses of observational and simulated data will be needed to (in)validate these predictions.
All the theoretical approaches predict an almost identical relationship for the normalized correlation of w and θ; refer to Eqs. (37c), (77), and (85c). The only difference arises due to the assumed value of P r t0 . The LSR, CSB, and EFB approaches assume P r t0 to be equal to 0.85, 1, and 0.8, respectively.

VI. DISCUSSIONS
In this section, we elaborate on a few limitations of the proposed LSR approach and how to overcome them in a practical manner.

A. Vertical Anisotropy of Turbulence
In this study, we have used Eq. (10b) to parameterize energy dissipation rate (ε). A more common practice would be to use Eq. (26) or its following variant: where, and c * 2 is an unknown coefficient. In Section II, we have implicitly assumed c * 2 A z to be a constant (c 2 ). In the literature, there is some evidence that the anisotropy parameter, A z , may be dependent on Ri g .
Based on observational and simulation data of turbulent air flows, Schumann and Gerz [58] proposed the fol- lowing empirical equation for 0 < Ri g < 1: According to this equation A z is weakly dependent on Ri g ; as a matter of fact, Schumann and Gerz [58] stated "the conclusions do not change much" if A z = 0.22 is used. Based on a DNS database, Basu et al. [7] reported A z to be approximately equal to 0.18 for 0 < Ri g < 0.2. The parameterizations of Canuto et al. [12], Kantha and Clayson [28], and Cheng et al. [15] predict gradual decrease of A z from near-neutral to strongly stratified conditions. Their predicted A  Figure 6) and Cheng et al. [15] (see their Figure 3c), in order to corroborate their respective formulations, do not portray any clear trends. A case in point are the wind tunnel measurements by Ohya [54] which exhibit random fluctuating behavior. Surprisingly, a strongly increasing trend of A z with respect to Ri g was predicted by large-eddy simulation data of [80] (see their Figure 4); this was in direct contradiction to their analytical prediction. Given this diversity in the A z -vs-Ri g relationship, we strongly recommend more research in this arena.
If we utilize Eq. (40a) instead of Eq. (10b), it is straightforward to re-derive all the equations reported in earlier sections. Some of the key equations are given here: Here c * P is an unknown coefficient and can be estimated following the procedure for c P . Furthermore, the quadratic equation for the turbulent Prandtl number be-comes: Ri g P r t + c 5 Ri g = 0. (43)

B. Imbalance of Production and Dissipation of TKE
In Eq. (9a), we have assumed that the production and dissipation of TKE balances exactly. Following Schumann and Gerz [58], we can define their ratio, termed a 'growth factor', as follows: It is likely that under strongly stratified condition, dissipation exceeds production. Thus, G can become less than unity for high values of Ri g . We can re-write Eq. (44) as follows: The key equations will then become: In this case, the quadratic equation for the turbulent Prandtl number becomes: (47) We would like to emphasize that the exact dependence of G on stability is not well studied in the literature. Schumann and Gerz [58] proposed an empirical (exponential decay) equation for G-vs-Ri g based on limited data. We hypothesize that for very stable conditions (Ri g > 1), G should be proportional to Ri −1 g . For practical applications, we propose the following heuristic parameterization for G: Thus, for Ri g < 1, G equals to 1. In other words, production and dissipation of TKE balance each other for weakly and moderately stable condition. However, the balance is lost (i.e., G < 1) for very stable conditions. If Eq. (48) is valid, then according to Eq. (46b), L X will be approximately equal to the buoyancy length scale (L b ) for very stable conditions. Perhaps more interestingly, if Eq. (48) indeed holds, Eq. (47) predicts that P r t should saturate to a constant value for Ri g > 1 . Such a prediction is not in agreement with some of the datasets reported in Fig. (1). However, it is consistent with the findings reported by [35] based on wind tunnel experiments and large-eddy simulations; refer to their Fig. 3.
We would like to emphasize that Eq. (48) is based on a heuristic argument and has not been verified yet. In our future work, we will leverage on Eq. (46e) to extract a reliable formulation for G.

C. Combined Scenario
For the most general case, one should account for the effects of both anisotropy and decay of TKE. In such a combined scenario, both A z and G terms will appear in the aforementioned equations. For example, the length scale equation will read: Similar to Eq. (31), for very stable condition (i.e., Ri g ≫ 1), the Prandtl number equation will be simplified to: Thus, the exact value of R f ∞ depends on a p , A z , G and c * P . Since stability dependencies of a p , A z and G are rather uncertain, empirical parameterizations for the combined terms (e.g., G/A z ) might be more practical for certain applications. High quality data from laboratory experiment (e.g., wind tunnel) and/or direct numerical simulation will be needed to derive such parameterizations.

VII. CONCLUSIONS
In this study, we have analytically derived an explicit relationship between the Prandtl number and the gradient Richardson number. Our derivation is rather simple from a mathematical standpoint and does not make elaborate assumptions beyond variance and sensible heat flux budget equations. Most of the unknown coefficients of the proposed relationship are easily estimated from wellknown surface layer similarity relationships. Our proposed Prandtl number formulation agrees very well with other competing approaches of quite different theoretical foundations and assumptions.
Our original analysis can be easily extended to include the effects of vertical anisotropy. It can also account for an imbalance of production and dissipation of TKE under very stable conditions. We have provided generalized formulations to account for these effects. However, these generalized formulations require stability-dependent formulations for a few parameters (e.g., A z , G) which are not well established in the literature. Currently, we are analyzing wind tunnel measurements and DNS-generated datasets to derive these formulations in a robust manner.
One of the limitations of the present study is that, for simplicity, it omits any discussion of internal gravity waves [61,67]. However, in stable boundary layers, specially under strong stratification, wave-turbulence interactions are extremely important. Thus far, only a handful of analytical studies have looked into such interactions [36,37,65,81]. We hope to further advance our proposed LSR approach along this direction in the future.

APPENDIX 1: PARAMETERIZATION OF THE PRESSURE-TEMPERATURE INTERACTION TERM
In the prognostic equation of sensible heat flux (u ′ i θ ′ ), a pressure-temperature interaction term Π i = − 1 ρ0 θ ′ ∂p ′ ∂xi appears [5,19]. This loss term is significant for atmospheric boundary layer (ABL) flows and requires a reliable parameterization. Using the product rule of calculus, Π i can be decomposed as follows [21,30,64]: Here p and ρ 0 denote pressure and a reference density, respectively. The symbol δ ik represents Kronecker delta. The first term on the right hand side of Eq. (51) represents turbulent diffusion of temperature field by pressure fluctuations and is sometimes neglected under the assumption of isotropy or using scaling argument [64]. As an alternative, in a number of modeling studies, it has been combined with the turbulent transport term, and in turn, the total term is parameterized via K-theory [47].
The second term Φ i = 1 ρ p ′ ∂θ ′ ∂xi is known as the pressure scrambling of the fluctuating temperature field. This term is split into three separate components representing different interactions [14,21]:  [2] 0.75 -Kantha and Clayson [30] 0.70 0.20 Nakanishi [52] 0.65 0.294 The term Φ T T i captures turbulence-turbulence interactions. Following Rotta's celebrated return-to-isotropy hypothesis [57], Monin [49] parameterized this term as follows: where τ R is the return-to-isotropy time scale. In the absence of external forces, this term relaxes turbulence to an isotropic state with zero overall heat flux [71]. Even though Eq. (53) is the most popular in the literature, alternative parameterizations for Φ T T i have been proposed in the past (please refer to [21]).
The mean shear-turbulence interaction is denoted by Φ S i and is parameterized as follows [2,13,21]: In the absence of significant subsidence or under quiescent synoptic condition, the vertical component (i.e., Φ S 3 ) is negligible in the ABL flows since u 3 = w ≈ 0; a comprehensive modeling study by [2] provides supporting results.
The following equation is often used for representing the buoyancy-turbulence interaction [21,38]: Even though this term is known to be important for nonneutral flows, quite interestingly, the well-known parameterizations of Mellor and Yamada [46] disregarded it. Over the years, various studies recommended different sets of values for a s and a p . Some of them are documented in Table III. Additionally, an empirical stabilitydependent formulation for a p was proposed by Wyngaard [75]: However, a p = 1 for Ri g > 1 does not lead to a physically meaningful solution when used in conjunction with Eq. (22). It is trivial to show that the solutions of the quadratic equation lead to two solutions: (i) P r t = P r t0 and (ii) P r t = Ri g . Neither of these solutions are plausible for the strongly stratified regime. In lieu of a realistic stability-dependent parameterization, in this study, we have decided to set a p as a fixed coefficient and have performed simple sensitivity analysis to quantify its influence on the overall predictions. By combining Eqs. (51)(52)(53)(54), the overall pressuretemperature interaction term for the vertical component of sensible heat flux can be simplified as follows: The terms on the right hand side of Eq. (57) are included in the simplified budget equation [i.e., Eq. (9c)] for sensible heat flux. The other terms of Eq. (9c) account for productions due to mean gradient (−σ 2 w Γ) and buoyancy (βσ 2 θ ).

APPENDIX 2: CO-SPECTRAL BUDGET (CSB) APPROACH
In this section, we re-derive the relevant equations of the co-spectral budget (CSB) approach following the footsteps of Katul et al. [32]. Along the way, we point out some of their (implicit) assumptions and differences to our newly proposed LSR approach. During this exercise, we noted certain sign errors in the original derivations of [32]. D. Li [41] confirmed our findings and pointed out additional sign errors in [32]. Fortunately, all these errors cancel out and do not have any effect on the key results. We have communicated our findings to G. G. Katul [31] and he has kindly verified them.
The starting point of the CSB approach is vertical sensible heat and momentum flux budget equations in wavenumber space: Here F θθ is the one dimensional temperature spectrum. k x denotes wavenumber in the along-wind direction. Both these equations assume steady-state condition. They neglect turbulent transport and molecular diffusion terms. Interestingly, the momentum flux equation also neglects the buoyancy term. The pressure-temperature and pressure-velocity interactions are parameterized as follows: Where F wθ and F uw are the cospectra between w-θ and w-u, respectively. τ (k x ) is a relaxation time-scale. A T , are constants which should be prescribed. Katul et al. [32] assumed: A T = A U = 1.8, and c CSB 1T = c CSB 1U = 3/5. Please note that Eq. (59a) does not include the commonly used buoyancy-turbulence interaction term [see Eq. (55) in Appendix 1]. Instead, it includes an unorthodox term which is proportional to the production term.
The production terms are expressed as follows: Here, the one dimensional vertical velocity spectrum is denoted by F ww . Please note that both these equations in [32] contain sign errors as pointed out by [41]. By combining Eqs. (58a), (59a), and (60a), we get: (61a) Similarly, by using Eqs. (58b), (59b), and (60b), we arrive at: Next, Katul et al. [32] assumed that F ww (k x ), F θθ (k x ), and τ (k x ) follow the inertial-range scaling behavior within the range k a ≤ k ≤ ∞ as hypothesized by Kolmogorov-Obukhov-Corrsin: Where N θ is simply half of dissipation rate of temperature variance (χ θ ). c CSB 0 and c CSB T are supposed to be universal constants. Katul et al. [32] assumed c CSB 0 = 0.65 and c CSB T = 0.80. Please note that, for simplicity, they assumed that the inertial-range scaling also holds in the dissipation range.
For low wavenumbers (0 ≤ k x ≤ k a ), the CSB approach assumes flat (i.e., white noise) spectra: Over the decades, several competing hypotheses (e.g., [8,9,44,48,59]) have been put forward to characterize the low wavenumber (aka buoyancy-range) spectra. We would like to point out that none of these hypotheses are in line with the assumption of the CSB approach. Furthermore, in the surface layer, there are ample evidence in the literature (e.g., [26,42]) that temperature spectra follow k −1 x scaling and not k 0 x scaling as assumed by the CSB approach.
In an analogous manner, we get from Eq. (61b) and Eqs. (62)-(63): This equation in [32] contains a sign error. From Eqs. (66a), (66b), and (67), we have: Note that Katul et al. [32] assumed A T = A U . The budget equations of TKE and σ 2 θ can be written as: Unfortunately, a sign error appears in the equation for N θ in Katul et al. [32].
Dividing Eq. (69b) by Eq. (69a) and using the definition of flux Richardson number (R f ), we can write: Hence, and Where, Katul et al. [32] assumed c CSB 0 , c CSB T and c CSB
From Eq. (71b), we can easily derive the following quadratic equation (not reported in previous CSB-related publications): and its roots are: (72b) Only the larger root is physically meaningful.
We would like to point out that Eq. (72b) is cast in a different analytical form than the original CSB formulation in [32] and follow-up studies. For neutral condition (Ri g = 0), according to Eq. (72b), P r t0 equals to 1. Whereas, according Eq. (37) of [32], P r t0 is undetermined for neutral condition.
Using Eqs. (64), (65), and (70), the ratio of turbulent potential and kinetic energies can be derived as follows: For neutral condition, Eq. (69a) simplifies to: Thus, Utilizing this equation in conjunction with Eqs. (64) and (67), after a little algebraic manipulation, we can derive the ratio of normalized momentum flux as: This equation is identical to the prediction by the LSR approach [see Eq. (36)]. Using Eqs. (64), (65), (66a), (66b), and (69b), we can deduce an expression for the normalized correlation between vertical velocity and potential temperature as follows: For neutral condition, by definition Q equals to one. Thus, P r t0 is also unity.

APPENDIX 3: ENERGY-AND FLUX-BUDGET (EFB) APPROACH
Over the past several years, Zilitinkevich and his coworkers have proposed the so-called energy-and fluxbudget (EFB) approach and its several modifications. In this appendix, we briefly discuss some of the salient features of this approach. We follow one of the later versions of the EFB approach as documented by Zilitinkevich et al. [82].
The EFB approach makes use of the steady-state budget equations for both sensible heat and momentum fluxes. As a reminder to the readers, our proposed LSR approach does not utilize the momentum flux equation.
In the case of the sensible heat flux equation, [82] parameterizes the pressure-temperature interaction term as follows: where c EF B θ is an unknown coefficient. In the LSR approach, we use the term a p to denote (1−c EF B θ ). Interestingly, Zilitinkevich et al. [82] neglects the commonly used turbulence-turbulence interactions [i.e., Eq. (53)] in the pressure-temperature interaction term. However, they use this exact term to parameterize the dissipation term (commonly neglected in the literature) of the sensible heat flux equation as follows: Where, τ ε is the dissipation time scale and c EF B F is an unknown coefficient, assumed to be equal to 0.25. Effectively, both the EFB and the LSR approaches use the same form of parameterized sensible heat flux equation. From this equation, with minor algebraic manipulations, [82] derived: Please refer to Eqs. (33a) and (33b) for the definitions of e p and e w , respectively. In the case of the momentum flux equation, Zilitinkevich et al. [82] makes several approximations. They neglect the dissipation term. In addition, they combine the buoyancy and pressure-velocity interaction terms and call it an 'effective dissipation rate'. This combined term is parameterized like a return-to-isotropy term. The resultant momentum flux equation is written as follows: where, c EF B τ is an unknown coefficient, assumed to be equal to 0.2. Thus, the eddy diffusivity can be represented as: Based on Eqs. (79b) and (80b), one can write: Zilitinkevich et al. [82] argued that if P r t → ∞ as Ri g → ∞, then in the limiting case: Even though this equation is only valid for Ri g → ∞, the EFB approach uses c EF B θ as a constant, being equal to 0.105, for all stability conditions. In our proposed LSR approach, the related coefficient is (1 − a p ) and we have also assumed it to be a constant in lieu of a reliable stability-dependent parameterization.
From the budget equations of TKE and variance of potential temperature, along with the definition of flux Richardson number (R f ), [82] derived the following ratios: and, Where c EF B P is an unknown coefficient. Based on available data, [82] assumed it to be equal to 0.86.
By plugging in Eq. (83b) in Eq. (81) and using the definition A z = ew e , one gets the following equation after simplification: For neutral condition (i.e., R f = 0), with the chosen values of c EF B τ and c EF B F , the EFB approach predicts P r t0 = 0.8.
Please note that Eq. (84) requires a parameterization for A z . Zilitinkevich et al. [82] proposed heuristic equations for the redistribution of TKE among various velocity components due to the effects of stratification. Those equations lead to a specific formulation for A z ; please refer to Eq. (50c) of [82]. Using limited data, they further assumed A (Rig =0) z = 0.2 and A (Rig →∞) z = 0.03. In Section VI A we have provided more information on A z .
It is straightforward to derive the following normalized variances and fluxes from the EFB approach (see [43]): In Section V B 4, we have compared these equations against the predictions from the LSR and the CSB approaches.