Dissipation rate of turbulent kinetic energy in stably stratified sheared flows

Over the years, the problem of dissipation rate of turbulent kinetic energy (TKE) in stable stratification remained unclear because of the practical impossibility to directly measure the process of dissipation that takes place at the smallest scales of turbulent motion. Poor representation of dissipation causes intolerable uncertainties in turbulenceclosure theory and thus in modelling stably stratified turbulent flows. We obtain a theoretical solution to this problem for the whole range of stratifications from neutral to limiting stable; and validate it via (i) direct numerical simulation (DNS) immediately detecting the dissipation rate and (ii) indirect estimates of dissipation rate retrieved via the TKE budget equation from atmospheric measurements of other components of the TKE budget. The proposed formulation of dissipation rate will be of use in any turbulence-closure models employing the TKE budget equation and in problems requiring precise knowledge of the high-frequency part of turbulence spectra in atmospheric chemistry, aerosol science, and microphysics of clouds.


Introduction
Until the present, the dependence of dissipation rate, ε K , of turbulent kinetic energy (TKE), E K , on static stability remained insufficiently understood. This caused principal difficulties in the theory of turbulence energetics and turbulence closure, and intolerable uncertainties in comprehending and modelling stably stratified turbulent flows. Traditionally, the dissipation rate is parameterized in terms of a turbulent length scale, l T , as ε K ∼ E 3/2 K /l T . This solves the problem in neutrally stratified boundary-layer flow, when the only length scale is the distance over the surface, z, so that l T ∼ z. However, in stratified flows, one more length scale appears, namely the Obukhov length scale, L, so that the ratio l T /z becomes an unknown function of z/L. To define this function we combine observational evidence with theoretical analyses. We employ the steady-state TKE budget equation to retrieve data on dissipation versus stability from uncountable data on wind profiles in the moderately stably stratified atmospheric surface layers, supplement this information with our own direct numerical simulation of turbulence in stably stratified Couette flow, and combine the collected empirical knowledge with asymptotic analysis of the TKE equation. The analyses reveal perfect equivalence of our asymptotic formulation of the velocity profile in extremely stable strati-fication and well-known log-linear velocity profile in moderately stable stratifications typical of the atmospheric surface layer -up to the coincidence of empirical dimensionless constants. This very lucky empirical finding yields universal formulation of the dissipation rate versus static stability, valid over the whole range of stratifications from neutral to extremely stable. The formulation is applicable to any stationary and horizontally homogeneous stably stratified sheared flows and can be used within any turbulence-closure model equipped with the TKE budget equation.
For certainty, we consider the dissipation rate of TKE in terms of dry atmosphere, where fluctuation of buoyancy, b = βθ, is proportional to fluctuation of potential temperature, θ ; β = g/T 0 is the buoyancy parameter, g is the gravitational acceleration, and T 0 is a reference value of absolute temperature. Since Kolmogorov (1942), ε K is expressed through the dissipation timescale, t T , or length scale, l T : This formulation is not hypothetical but just defines the scales t T and l T , so that Eq. (1) merely expresses one unknown, ε K , through another, t T or l T . In neutrally stratified boundary-layer flows, the only principal length scale is the height over the surface, z; so that l T is proportional to z, which yields where C l is a dimensionless constant to be determined empirically. Stratification involves the Obukhov length scale: where τ is the absolute value of vertical turbulent flux of momentum τ = (τ, 0), and F z is vertical turbulent flux of potential temperature (Obukhov, 1946). The restraining effect of stable stratification on turbulence is characterized by the dimensionless height, z/L; gradient Richardson number, or flux Richardson number, where ∂U/∂z and ∂ /∂z are vertical gradients of mean wind velocity, U = (U, V ), and mean potential temperature, . Then the dimensionless dissipation rate, ε K z/E 3/2 K , is no longer a constant but depends on stratification (z/L, Ri, or Ri f ). Until recently, practically nothing was known about this dependence beyond the interval of stratifications covered by observations in atmospheric surface layer: 0 < Ri < 0.2, which corresponds to 0 < z/L < 10.
2 Dissipation rate in the steady-state stably stratified sheared flows We consider horizontally homogeneous stationary boundarylayer flow in semi-space z > 0 as an idealized model of atmospheric surface layer. Here, the familiar TKE budget equation expresses the dissipation rate, ε K , through τ , F z , and ∂U/∂z: With increasing static stability, Ri f obviously increases but (because ε T > 0) remains limited, which is why it must tend to a finite limit: Ri f → R ∞ < 1. Atmospheric data and results from direct numerical simulation (DNS) demonstrated below in Figs. 2 and 3 confirm such behaviour and yield quite a certain estimate of R ∞ = 0.2. Then, substituting R ∞ for Ri f in Eq. (5) yields asymptotic expression of the velocity gradient in extremely stable stratification: Here, the von Kármán constant, k, is inserted in numerator and denominator just to highlight consistency of Eq. (7) with the well-known formulation of the velocity gradient in weakly and moderately stable stratifications typical of atmospheric surface layers, obtained in the familiar Monin-Obukhov similarity theory (MOST): where C u = 2 is a well-established dimensionless empirical constant (Monin and Obukhov, 1954;Monin and Yaglom, 1971;Garratt, 1992;Stull, 1997). Originally, Eq. (8) was derived as the first term in the Taylor expansion of the dimensionless velocity gradient, M = kz/τ 1/2 ∂U/∂z, considered in MOST as a universal function of z/L. Subsequently, it was revealed that Eq. (8) with C u = 2 is valid over the whole range of z/L observed in atmospheric surface layers, 0 < z/L < 10, which corresponds to quite small gradient Richardson numbers: 0 < Ri < 0.2 (Monin and Yaglom, 1971). By this means, Eq. (8) with C u = 2 based on massive atmospheric data for moderately stable stratifications yields, as z/L → ∞, precisely the same limit as Eq. (7) with k/R ∞ = 2 resulting from the conventional value of k = 0.4 and new estimate of the critical flux Richardson number: R ∞ = 0.2 obtained from topical DNS and available atmospheric data for maximal stable stratifications (Fig. 2). This lucky coincidence just means that Eq. (8) with C u = k/R ∞ holds true in any stable stratification: Atmos. Chem. Phys., 19, 2489-2496, 2019 www.atmos-chem-phys.net/19/2489/2019/ Then, substituting Eq. (9) for ∂U/∂z into Eqs. (5) and (6) yields the following simple relations linking Ri f with z/L: and exact formulation of the TKE dissipation rate as dependent on static stability: It is worth noting that R ∞ can be derived from wellestablished phenomenological constants of turbulence characterizing the inertial subrange (Katul et al., 2014). The actual value in this case is slightly higher (R ∞ = 0.25) but still within a reasonable range.
To comprehensively validate the above analyses, we performed DNS of stably stratified Couette flow, namely the plain-parallel flow between two horizontal plates separated in the vertical by a distance d, and moving with constant velocity in opposite directions. To assure the accuracy of numerical simulations, we employed two DNS codes: one developed at the Institute of Numerical Mathematics, Russian Academy of Sciences, Moscow State University (hereafter INM-RAS) and another developed at the Institute of Applied Physics, Russian Academy of Sciences (IAP-RAS). Despite using two different codes developed separately and characterized by different spatial and temporal schemes, resolutions, and statistical averaging, the two types of DNSs have shown quite consistent results that can be considered a crossvalidation. For a detailed description of the above numerical models used, see Mortikov (2016), Mortikov et al. (2019), and Druzhinin et al. (2016).
In our DNS total (turbulent + molecular) fluxes of momentum, τ , and potential temperature, F z , are practically equal to the turbulent fluxes elsewhere beyond narrow nearwall sublayers where molecular transports dominate. In the Couette flow, the total fluxes are constant with height, exactly as in the surface-layer flows. Similarly, the flux-profile relations linking τ and F z with vertical gradients of mean velocity, ∂U/∂z, and potential temperature, ∂ /∂z, as well as the budget equations for turbulent energies (in particular Eq. 6), are the same as in the surface-layer flows. The only difference is in the geometry of domains illustrated in Fig. 1.
Following Obukhov (1942), we distinguish between "absolute geometry" characterized by the usual height over the surface, z, and "internal geometry" characterized in Couette flow by specific vertical coordinate, z, dictated by conformal mapping of the Couette-flow domain (0 < z < d) into the semi-space: This coordinate reflects the influences of lower and upper walls on the fluid flow.
In semi-space, the "internal geometry" coincides with "absolute geometry": z = z. Thus, the vertical structure of the Couette flow in terms of z coincides with the vertical structure of the surface-layer flow in terms of z. This allows for showing together the genuine dissipation rate calculated from DNS: ε K = ν (∂u i /∂x k )(∂u i /∂x k ) , where ν is kinematic viscosity, and ε K = τ ∂U/∂z + βF z retrieved from atmospheric observations assuming the steady-state TKE budget.
In Figs. 2-4 we show our DNS data together with atmospheric data on the dissipation rate retrieved from observations in the surface layers indirectly: via the Kolmogorov −5/3 power law from measured spectra of TKE in the inertial subrange (Pearson et al., 2002), and via the steady-state TKE budget Eq. (6) from the measured turbulent fluxes of momentum, τ , and potential temperature, F z , and vertical gradient of mean wind velocity, ∂U/∂z.
In these figures DNS data are shown by bold coloured dots, and atmospheric data by light-grey symbols (Kadantsev et al., 2019). Figure 2 shows flux Richardson number, Ri f = βF z (τ · ∂U/∂z) −1 , versus dimensionless height, z/L, in Couette flow; or versus z/L in the atmospheric surface layer. The black curve is plotted after Eq. (10) taking the conventional value of the von Kármán constant, k = 0.4, and our estimate of the maximal flux Richardson number, R ∞ = 0.2, resulting from the best fit of Eq. (10) to DNS data. Notably, total (turbulent + molecular) fluxes of momentum, τ , and potential temperature, F z , in Couette flow are constant across the flow which assures a very certain specification of Ri f and L, and makes our DNS most suitable for calibrating the theory. We recall that Eqs. (10) and (11) are relevant to the well-developed turbulence regime where molecular transports are negligible, so that turbulent fluxes practically coincide with total fluxes. In our DNS this is true, except for the narrow transition layers dominated by molecular transport near the lower and upper walls: 0 < z < 50ν/τ 1/2 . Data from these layers are indicated by dark-grey points. The light-grey symbols show atmospheric data from the following sources: research observatory Tiksi in eastern Siberia near the Arctic Ocean coast (Grachev et al., 2018), offshore oceanographic platform in the Black Sea (Repina et al., 2009); and acoustic soundings over arid steppe in the Republic of Kalmykia in Southern Russia (Vazaeva et al., 2017). In spite of inevitable heterogeneity, non-stationarity, and other side effects, atmospheric data correlate quite well with DNS data. Figure 3 shows dimensionless dissipation rate, ε K z/τ 3/2 , versus z/L after Eq. (11) and atmospheric data; and ε K z/τ 3/2 versus z/L after DNS of Couette flow. All notations are the same as in Fig. 2. The theoretical curve plotted after Eq. (11) Figure 1. Usual height, z, and vertical coordinate, z, defined by Eq. (12) characterizing "absolute" and "internal" geometry of the domain, respectively. Panel (a) shows semi-space, z > 0, where z = z. Panel (b) shows the layer between two horizontal walls, 0 < z < d, where z coincides with the distance from nearest surface in its close vicinity, but essentially depends on the distances from both surfaces in the central part of the domain. Empirical data used for the calibration are obtained in two series of DNS runs employing INM-RAS code (red dots) and IAP-RAS code (blue dots). Atmospheric data are taken from Arctic coastal observatory Tiksi (light-grey diamonds), Black Sea offshore platform (lightgrey squares), and acoustic soundings in Kalmykia steppe (light-grey stars). Dark-grey dots belong to the very narrow near-surface layer: 0 < z < 50ν/τ 1/2 . The solid black line shows Eq. (10) with the conventional value for the von Kármán constant, k = 0.4, and new empirical value of R ∞ = 0.2 just obtained from the best fit of Eq. (10) to DNS data from elsewhere beyond the layer 0 < z < 50ν/τ 1/2 , where molecular transports are significant and Eq. (10) is not necessarily relevant.
with k = 0.4 and R ∞ = 0.2 is fully consistent with experimental data, except for the narrow transition layer 0 < z < 50ν/τ 1/2 where Eq. (11) is irrelevant. Hence, Fig. 3 justifies the stability dependence of dissipation rate, Eq. (11), and provides additional confirmation to the empirical estimate of R ∞ = 0.2.
In Fig. 2 we consider the flux Richardson number: Ri f = βF z τ ·∂U/∂z , where turbulent fluxes (disregarding molecular contributions in the transition layer) appear in both the numerator and denominator. Hence uncertainties in both fluxes are somehow compensated. This is not the case in Fig. 3 show-ing ε K z/τ 3/2 vs. z/L: the dissipation rate in the numerator is just total dissipation, whereas the momentum flux in the denominator disregards the molecular contribution. This causes the ugly looking but only natural dark-grey points on the left side of Fig. 3. Figure 3. Dimensionless dissipation rate in stable stratification, ε K z/τ 3/2 , versus z/L in Couette flow or ε K z/τ 3/2 versus z/L in the atmospheric surface layer. Empirical data are from the same sources as in Fig. 2. The solid black line shows Eq. (11) with k = 0.4 and R ∞ = 0.2. Dark-grey dots belong to the very narrow near-surface layer: 0 < z < 50ν/τ 1/2 , where molecular transports are significant and Eq. (11) is not relevant.

Turbulent length scales and general criterion of stratification
The concept of the TKE dissipation rate directly relates to the definition of the turbulent timescale, t T ≡ E K /ε K , and length scale, Then Eq. (11) defines l T as a function of z/L: It has the asymptotic limits: where the limits of E K /τ in neutral stratification and in extremely stable stratification are just dimensionless constants. Our DNS yield the following estimates: (E K /τ ) 0 ≈ 4 and (E K /τ ) ∞ ≈ 11. The length scale similar to Eq. (13) was already revealed as inherent to spectra of turbulence in unstably stratified boundary-layer flows (Glazunov, 2014).
We emphasize that l T is the scalar characterizing turbulence as a whole. Contrastingly, turbulent mixing in differ-ent directions is characterized by the mixing-length vector l Ti ≡ E 1/2 Ki t T (i = 1, 2, 3) with generally different streamwise (i = 1), transverse (i = 2), and vertical (i = 3) components. We emphasize principal difference between the scalar length scale and vector mixing length. In literature, the words "turbulent length scale" and "turbulent mixing length" are often used as interchangeable. This causes intolerable confusion because different components of the mixing length differently depend on static stability (Zilitinkevich et al., 2013).
The above analyses are done for the simplest surface-layer (or Couette) flow, where dimensionless height z/L (or z/L) plays the role of criterion quantifying the effect of stratification on turbulence. Luckily, our major result (Eqs. 11 and 13) can be easily extended to a wide range of stratified turbulent flows. We recall that stratified turbulence is characterized, besides E K , by turbulent potential energy (TPE), E P = 1 2 β θ 2 /(∂ /∂z). Hence the effect of stratification on turbulence can be quantified by the "energy Richardson number" defined as In contrast to traditional criteria, such as Ri (Eq. 4), Ri f (Eq. 5), or z/L, the energy Richardson number criterion is valid in heterogeneous and non-stationary flows, for any mechanisms of generation of turbulence (including breaking waves, oscillating grid, etc.) and in flows with complex geometry.
Expressing the dissipation rates of TKE and TPE in the steady state through the dissipation timescale, l T ≡ E 3/2 K /ε K , the budget equations for TKE and TPE become where C P is dimensionless parameter quantifying the difference between the dissipation rates of TKE and TPE (Zilitinkevich et al., 2013). Equations (16) and (17) in combination with Eq. (10) yield the following relations linking Ri E with Ri f or z/L: Figure 4 shows Ri E versus z/L or z/L (like in previous figures) after our DNS and atmospheric observations. The theoretical curve is plotted after Eq. (18) taking k = 0.4, R ∞ = 0.2, and an empirical estimate of the dimensionless parameter C P = 0.62 just obtained from the best fit of Eq. (18) to DNS data. Experimental data reveal the asymptotic limit: Then, using Eq. (18) to express z/L through Ri E , Eq. (11) in terms of Ri E becomes: where ε K(neutral) is dissipation rate in neutral stratification. In the surface layer ε K(neutral) = τ 3/2 /kz; but generally ε K(neutral) depends on concrete energy-generation mechanisms and geometry of flow.
There is an essential advantage to Ri E as criterion of stratification in numerical modelling. Turbulent fluxes are usually calculated through familiar diagnostic down-gradient formulations: τ = −K M ∂U/∂z and F z = −K H ∂ /∂z, where K M is eddy viscosity and K H is eddy conductivity. Then, finitedifference approximation of the gradients causes uncertainties in τ and F z and, hence, in the Obukhov length, L (Eq. 3), flux Richardson number, Ri f (Eq. 5), and gradient Richardson number, Ri (Eq. 4). Contrastingly, TKE and TPE are defined from prognostic budget equations accounting for turbulent diffusion which smooths the energies and assures quite a certain calculation of Ri E .

Concluding remarks
The dissipation rate of TKE, ε K , as dependent on static stability over years remained uncertain because of the impossibility of direct measurement of ε K . Admittedly, ε K can be retrieved via the TKE budget equation from the measured turbulent fluxes, τ and F z , and mean-velocity gradient, ∂U/∂z, and also via the Kolmogorov −5/3 power law from the measured spectra of TKE in the inertial subrange. However, these methods are only justified in stationary and horizontally homogeneous flows and require fully controlled conditions. These provisions, practically unachievable in atmospheric experiments, make estimates of ε K from atmospheric observations rather uncertain. The wide spread of atmospheric data is clearly seen in our figures. Moreover, available atmospheric data cover only weakly to moderately stable stratifications typical of the surface layer. To avoid these difficulties, we performed topical DNS of the steady-state stably stratified turbulent Couette flows up to the strongest attainable stratifications, combined direct data from DNS with data retrieved from atmospheric observations, and employed theoretical analysis to reveal asymptotic behaviour of the mean velocity gradient and the dissipation rate in extremely stable stratification, namely as z/L → ∞ where L is the Obukhov length scale.
By providential coincidence, the formulations happen to be precisely the same in the asymptotic limit z/L → ∞ and in the weakly stable stratifications 0 < z/L < 10 typical of atmospheric surface layer. This yields simple analytical formulations of the dimensionless velocity gradient, kz/τ 1/2 ∂U/∂z, and dissipation rate, kz/τ 3/2 ε K , as uni-