Parameterizing the Energy Dissipation Rate in Stably Stratified Flows

We use a database of direct numerical simulations to derive parameterizations for energy dissipation rate in stably stratified flows. We show that shear-based formulations are more appropriate for stable boundary layers than commonly used buoyancy-based formulations. As part of the derivations, we explore several length scales of turbulence and investigate their dependence on local stability.


I. INTRODUCTION
Energy dissipation rate is a key variable for characterizing turbulence [1]. It is a sink term in the prognostic equation of turbulent kinetic energy (TKE; e): ∂e ∂t + ADV = BN C + SHR + T RP + P RC − ε, (1) where, ε is the mean energy dissipation rate. The terms ADV , BN C, SHR, T RP , and P RC refer to advection, buoyancy production (or destruction), shear production, transport, and pressure correlation terms, respectively. Energy dissipation rate also appears in the celebrated "-5/3 law" of Kolmogorov [2] and Obukhov [3,4]: where, E(κ) and κ denote the energy spectrum and wavenumber, respectively. In field campaigns or laboratory experiments, direct estimation of ε has always been a challenging task as it involves measurements of nine components of the strain rate tensor. Thus, several approximations (e.g., isotropy, Taylor's hypothesis) have been utilized and a number of indirect measurement techniques (e.g., scintillometers, lidars) have been developed over the years. In parallel, a significant effort has been made to correlate ε with easily measurable meteorological variables. For example, several flux-based and gradient-based similarity hypotheses have been proposed [e.g., [5][6][7][8].
In addition, a handful of papers also attempted to establish relationships between ε and either the vertical * sukanta.basu@gmail.com † drpinghe@umich.edu ‡ awdemarc@ncsu.edu velocity variance (σ 2 w ) or TKE (e). One of the first relationships was proposed by Chen [9]. By utilizing the Kolmogorov-Obukhov spectrum (i.e., Eq. 2) with certain assumptions, he derived: Since this derivation is only valid in the inertial-range of turbulence, a band-pass filtering of vertical velocity measurements was recommended prior to computing σ w . A few years later, Weinstock [10] revisited the work of [9] and again made use of Eq. 2, albeit with different assumptions (see Appendix 2 for details). He arrived at the following equation: where, N is the so-called Brunt-Väisäla frequency. Using observational data from stratosphere, Weinstock [10] demonstrated the superiority of Eq. 4 over Eq. 3. In a recent empirical study, by analyzing measurements from the CASES-99 field campaign, Bocquet et al. [11] proposed to use ε as a proxy for σ 2 w . In the present work, we quantify the relationship between ε and e (as well as between ε and σ w ) by using turbulence data generated by direct numerical simulation (DNS). To this end, we first compute several wellknown "outer" length scales (e.g., buoyancy length scale and Ozmidov scale), normalize them appropriately, and explore their dependence on local stability. Next, we investigate the inter-relationships of certain (normalized) outer length scales (OLS) which portray qualitatively similar stability-dependence. By analytically expanding these relationships, we arrive at two ε-e and two ε-σ w formulations; only the shear-based formulations portray quasi-universal scaling.
The organization of this paper is as follows. In Sect. 2, we describe our DNS runs and subsequent data analyses. Simulated results pertaining to various length scales are included in Sect. 3. The ε-e and ε-σ w formulations are derived in Sect. 4. A few concluding remarks, including implications of our results for atmospheric modeling, are made in Sect. 5. In order to enhance the readability of the paper, either a heuristic or an analytical derivation of all the length scales is provided in Appendix 1. Given the importance of Eq. 4, its derivation is also summarized in Appendix 2. Last, in Appendix 3, we elaborate on the normalization of various variables which are essential for the post-processing of DNS-generated data.

II. DIRECT NUMERICAL SIMULATION
Over the past decade, due to the increasing abundance of high-performance computing resources, several studies probed different types of stratified flows by using DNS [e.g., [12][13][14][15][16][17][18][19]. These studies provided valuable insights into the dynamical and statistical properties of these flows (e.g., intermittency, structure parameters). In the present study, we use a DNS database which was previously generated by using a massively parallel DNS code, called HERCULES [20], for the parameterization of optical turbulence [21].
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 turbulent kinetic energy and its mean dissipation are computed as follows (using Einstein's summation no-tation): In these 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.

III. LENGTH SCALES
In this section, we discuss various length scales of turbulence. To enhance the readability of the paper, we do not elaborate on their derivations or physical interpretations here; for such details, the readers are requested to peruse Appendix 1.
From the DNS-generated data, we first calculate the integral length scale (L) and Kolmogorov length scale (η). They are defined as [22,23]: In Fig. 1, normalized values of L and η are plotted against the gradient Richardson number (Ri g = N 2 /S 2 ); where, S is the magnitude of wind shear. It is evident that the simulations with larger cooling rates result in smaller L as would be physically expected. In contrast, η marginally increases with higher stability due to lower ε. The ratio of L to η decreases from about 100 to 20 as stability is increased from weakly stable condition to strongly stable condition.
Next, we compute four outer length scales: Ozmidov (L OZ ), Corrsin (L C ), buoyancy (L b ), and Hunt (L H ). They are defined as [24][25][26][27][28][29][30][31]: Please note that, in the literature, L b and L H have also been defined as σ w /N and σ w /S, respectively. Both L OZ and L C are functions of ε, a microscale variable. In contrast, L b and L H only depend on macroscale variables. Both shear and buoyancy prefer to deform the larger eddies compared to the smaller ones [15,[32][33][34]. Eddies which are smaller than L C or L H are not affected by shear. Similarly, buoyancy does not influence the eddies of size less than L OZ or L b . In other words, the eddies can be assumed to be isotropic if they are smaller than all these OLSs.
Since L changes across the simulations, all 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.
Both normalized L OZ and L b decrease monotonically with Ri g ; however, the slopes are quite different. The length scales L C and L H barely exhibit any sensitivity to Ri g (except for Ri g > 0.1). Even for weakly-stable condition, these length scales are less than 25% of L.
Based on the expressions of the OLSs (i.e., Eqs. 7) and the definition of the gradient Richardson number, we can write: Thus, for Ri g < 1, one expects L C < L OZ and L H < L b . Such relationships are 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.
Owing to their similar scaling behaviors, L b /L against L OZ /L are plotted in Fig. 3 (left panel). Once again, all the simulated data collapse nicely in a quasi-universal (nonlinear) curve. Since in a double-logarithmic representation (not shown) this curve is linear, we can write: where, m is an unknown power-law exponent. By using L b ≡ e 1/2 /N and the definitions of L OZ and L, we arrive at: Further simplification leads to: ε = eN ; please note that the exponent m cancels out in the process. Instead of e 1/2 , if we utilize σ w in the definitions of L b and L, we get: ε = σ 2 w N . This equation is identical to Eq. 4 which was derived by Weinstock [10]. His derivation, based on inertial-range scaling, is summarized in Appendix 2.
In the right panel of Fig. 3, we plot L H /L versus L C /L. Both these normalized length scales have limited ranges; nonetheless, they are proportional to one another. Like Eq. 9, we can write in this case: where, n is an unknown power-law exponent. The expansion of this equation leads to either ε = eS or ε = σ 2 w S, depending on the definition of L H and 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.

IV. PARAMETERIZING THE ENERGY DISSIPATION RATE
Earlier in Fig. 3, we have plotted normalized OLS values against one another. It is plausible that the apparent data collapse is simply due to self-correlation as same variables (i.e., L, N , and S) appear in both abscissa and ordinate. To further probe into this problematic issue, we produce Fig. 4. Here, we basically plot normalized ε as functions of normalized eN , eS, σ 2 w N , and σ 2 w S, respectively. These plots have completely independent abscissa and ordinate terms and do not suffer from self-correlation. Please note that the appearance of Re b and Ri b in these figures are due to the normalization of variables in DNS; Appendix 3 provides further details. Throughout the paper, the subscript "n" is used to denote a normalized variable.
It is clear that the plots in the left panel of Fig. 4, which involve N , do not show any universal scaling. For low CR values, normalized ε values do not go to zero; this behavior is physically realistic. One cannot expect ε to go to zero for neutral condition (i.e., N → 0). With increasing cooling rates, the curves seem to converge to an asymptotic curve which passes through the origin. As e or σ w continually reduces with increasing stability, one does expect ε to approach zero.
In a seminal paper, Deardorff [35] proposed a parameterization for ε which for strongly stratified condition approaches to 0.25eN . In Fig. 4 (top-left panel), we overlaid this line ε = 0.25eN on the DNS-generated data. Clearly, it only overlaps with the simulated data at the strongly stratified region. If ε = 0.25eN is used in the definition of L OZ , after simplification, one gets L OZ = L b /2. The line L b = 2L OZ is drawn in Fig. 3. As would be anticipated, it only overlaps with the simulated data when the OLS values are the smallest (signifying strongly stable

condition).
Compared to the left panels, the right panels of Fig. 4 portray very different scaling characteristics. All the data collapse on quasi-universal curves remarkably; especially, for the ε ≈ eS case. The slopes of the regression lines, estimated via conventional least-squares approach and bootstrapping [36,37], are shown on these plots. Essentially, we have found: We note that our estimated coefficient 0.63 is within the range of values reported by Schumann and Gerz [38] from laboratory experiments and large-eddy simulations (please refer to their Fig. 1).
In summary, neither ε = eN nor ε = σ 2 w N are appropriate parameterizations for weakly or moderately stratified conditions; they may provide reasonable predictions for very stable conditions. In contrast, the shear-based parameterizations should be applicable from a wide range of stability conditions, from near-neutral to at least Ri g ≈ 0.2. Since within the stable boundary layer, Ri g rarely exceeds 0.2 [see 39], we believe Eq. 12a or Eq. 12b will suffice for most practical boundary layer applications. However, for free atmosphere, where Ri g can exceed O(1), Deardorff's parameterization (i.e., ε = 0.25eN ) might be a more viable option. Unfortunately, we cannot verify this speculation using our existing DNS dataset.

V. CONCLUDING REMARKS
The boundary-layer community almost always utilizes buoyancy-based energy dissipation rate parameterizations for numerical modeling studies. Only a few studies [e.g., 40,41] have explored the possibilities of combining buoyancy-based and shear-based parameterizations. In this paper, we demonstrated that shear-based parameterizations are much more appropriate for stable boundary layer flows as long as Ri g does not exceed 0.2. Our DNSbased results are in complete agreement with the theoretical work (supported by numerical results) of Hunt et al. [28]. They concluded: "...when the Richardson number is less than half, it is the mean shear ... (rather than the buoyancy forces) which is the dominant factor that determines the spatial velocity correlation functions and hence the length scales which determine the energy dissipation or rate of energy transfer from large to small scales." We sincerely hope that our community will begin to appreciate this incredible insight of Hunt and his co-workers in future modeling studies (including large-eddy simulations).
Before closing, we would like to emphasize that the proposed shear-based parameterizations are only applicable away from the surface. Near the surface, due to the blocking effect [see 28,29], L C or L H cannot be a representative length scale. They should be properly combined with an explicit parameterization involving height above ground (e.g., the harmonic mean of 0.4z and L H ).

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 thoughtprovoking 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.  [22] and Pope [23] provided a heuristic derivation of the integral length scale. 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 (L) can be approximated as e 3/2 /ε.
In the literature, the autocorrelation function of the longitudinal velocity series is commonly used to derive an estimate of the integral length scale (L 11 ). The relationship between L and L 11 is discussed by Pope [23].
b. Kolmogorov Length Scale: Pope [23] paraphrased the first similarity hypothesis of Kolmogorov [2] as (the mathematical notations were changed by us for consistency): "In every turbulent flow at sufficiently high Reynolds number, the statistics of the smallscale motions (l ≪ L) have a universal form that is uniquely determined by ν and ε." Based on ν and ε, the following length scale can be formulated using dimensional analysis: η ≡ ν 3 ε 1/4 . At this scale, TKE is converted into heat by the action of molecular viscosity.
c. Ozmidov Length Scale: Dougherty [25] and Ozmidov [26] independently proposed this length scale. Here, we briefly summarize the derivation of Ozmidov [26]. Based on Kolmogorov [2], the first-order moment of the velocity increment (∆u) in the vertical direction (z) can be written as: where, the overlines denote ensemble averaging. Using this equation, the vertical gradient of longitudinal velocity component can be approximated as: Similar equation can be written for the vertical gradient of the lateral velocity component ( ∂v ∂z ). Thus, the magnitude of wind shear (S) can be written as: By definition, Ri g = N 2 /S 2 . Thus, Ozmidov [26] assumed that for a certain critical Ri g (which is assumed to be an unknown constant), ∆z becomes the representative outer length scale (L OZ ). Thus, Eq. 16 can be re-written as: The unknown proportionality constant is a function of the critical Ri g and is assumed to be on the order of one. d. Buoyancy Length Scale: The following heuristic derivation is based on Brost and Wyngaard [27] and Wyngaard [31]. In an order-of-magnitude analysis, the inertia term of the Navier-Stokes equations, can be written as: where, L s , T s , and U s represent certain length, time, and velocity scales, respectively. In a similar manner, the buoyancy term can be approximated as: where, Θ 0 and θ ′ denote a reference temperature and temperature fluctuations, respectively. Equating the inertia and the buoyancy terms, we get: For stably stratified flows, either e 1/2 or σ w can be used as an appropriate velocity scale. Accordingly, the length scale (L s ) can be approximated as e 1/2 N or σw N . In the literature, this length scale is commonly known as the buoyancy length scale (L b ).
e. Corrsin Length Scale: The derivation of Corrsin [24] leverages on a characteristic spectral time-scale, T s (κ), which is representative of the inertial-range. Based on dimensional argument, Onsager [43] proposed: In order to guarantee local isotropy in the inertial-range, Corrsin [24] hypothesized that T s (κ) must be much smaller than the time-scale associated with mean shear (S). In other words, Using the "-5/3 law" of Kolmogorov [2] and Obukhov [3,4], this equation can be re-written as: If we assume that for a specific wavenumber κ = 1/L C , the equality holds in Eq. 23, then we get: From this equation, we can estimate L C as defined earlier in Eq. 7b. f. Hunt Length Scale: Hunt [28] hypothesized that in stratified shear flows, ε is controlled by mean shear (S) and σ w . From dimensional analysis, it follows that: The associated length scale, L H , is assumed to be on the order of σ w /S.
The starting point of Weinstock's derivation was "-5/3 law" of Kolmogorov [2] and Obukhov [3,4]. He integrated this equation in the wavenumber space and set the upper integration limit to infinity. The lower integration limit was fixed at the buoyancy wavenumber (κ b ). Furthermore, he assumed that the eddies are isotropic for wavenumbers larger than κ b (i.e., in the inertial and viscous ranges). His derivation can be summarized as: Weinstock [10] assumed that κ b can be parameterized by N σw (basically, the inverse of the buoyancy length scale L b ). By plugging this parameterization into Eq. 26 and simplifying, we get:

APPENDIX 3: NORMALIZATION OF DNS VARIABLES
In DNS, the relevant variables are normalized as follows: After differentiation, we get: (29d) 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 and TKE can be normalized as: Following the above normalization approach, we can also derive the following relationship for the energy dissipation rate: In order to expand ε = eN , we use Eq. 32, Eq. 33d, and Eq. 34 as follows: In a similar manner, ε = eS can be re-written as: