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 turbulence-closure 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.

1 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 ${\mathit{\epsilon }}_{\mathrm{K}}\sim E_{\mathrm{K}}^{\mathrm{3}/\mathrm{2}}/l_{\mathrm{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 stratification 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, θ; $\mathit{\beta }=g/T_{\mathrm{0}}$ is the buoyancy parameter, g is the gravitational acceleration, and T₀ is a reference value of absolute temperature. Since Kolmogorov (1942), ε_K is expressed through the dissipation timescale, t_T, or length scale, l_T:

$$\begin{array}{}\text{(1)}& {\mathit{\epsilon }}_{\mathrm{K}}={\displaystyle \frac{E_{\mathrm{K}}}{t_{\mathrm{T}}}}={\displaystyle \frac{E_{\mathrm{K}}^{\mathrm{3}/\mathrm{2}}}{l_{\mathrm{T}}}}.\end{array}$$

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

$$\begin{array}{}\text{(2)}& l_{\mathrm{T}}=C_{l}z,{\mathit{\epsilon }}_{\mathrm{K}}={\displaystyle \frac{E_{\mathrm{K}}^{\mathrm{3}/\mathrm{2}}}{C_{l}z}},\end{array}$$

where C_l is a dimensionless constant to be determined empirically.

Stratification involves the Obukhov length scale:

$$\begin{array}{}\text{(3)}& L={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{3}/\mathrm{2}}}{-\mathit{\beta }{F}_z}},\end{array}$$

where τ is the absolute value of vertical turbulent flux of momentum $\mathit{\tau }=\left(\mathit{\tau },\mathrm{0}\right)$, 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,

$$\begin{array}{}\text{(4)}& \mathit{Ri}={\displaystyle \frac{\mathit{\beta }\partial \mathrm{\Theta }/\partial z}{{\left(\partial U/\partial z\right)}^{\mathrm{2}}}};\end{array}$$

or flux Richardson number,

$$\begin{array}{}\text{(5)}& {\mathit{Ri}}_{\mathrm{f}}={\displaystyle \frac{\mathit{\beta }{F}_z}{\mathit{\tau }\cdot \partial U/\partial z}},\end{array}$$

where $\partial U/\partial z$ and $\partial \mathrm{\Theta }/\partial z$ are vertical gradients of mean wind velocity, $U=\left(U,V\right)$, and mean potential temperature, Θ. Then the dimensionless dissipation rate, ${\mathit{\epsilon }}_{\mathrm{K}}z/E_{\mathrm{K}}^{\mathrm{3}/\mathrm{2}}$, 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: $\mathrm{0}<\mathit{Ri}<\mathrm{0.2}$, which corresponds to $\mathrm{0}<z/L<\mathrm{10}$.

2 Dissipation rate in the steady-state stably stratified sheared flows

We consider horizontally homogeneous stationary boundary-layer 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 $\partial U/\partial z$:

$$\begin{array}{}\text{(6)}& {\mathit{\epsilon }}_{\mathrm{K}}=\mathit{\tau }{\displaystyle \frac{\partial U}{\partial z}}+\mathit{\beta }{F}_z=\mathit{\tau }{\displaystyle \frac{\partial U}{\partial z}}\left(\mathrm{1}-{\mathit{Ri}}_{\mathrm{f}}\right).\end{array}$$

With increasing static stability, Ri_f obviously increases but (because ε_T>0) remains limited, which is why it must tend to a finite limit: ${\mathit{Ri}}_{\mathrm{f}}\to R_{\mathrm{\infty }}<\mathrm{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:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial U}{\partial z}}\to {\displaystyle \frac{\mathrm{1}}{R_{\mathrm{\infty }}}}{\displaystyle \frac{{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}}{L}}\equiv {\displaystyle \frac{{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}}{L}}\left({\displaystyle \frac{k}{R_{\mathrm{\infty }}}}{\displaystyle \frac{z}{L}}\right)\mathrm{as}{\displaystyle \frac{z}{L}}\to \mathrm{\infty }.\end{array}$$

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):

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial U}{\partial z}}={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}}{kz}}\left(\mathrm{1}+C_u{\displaystyle \frac{z}{L}}\right),\end{array}$$

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, ${\mathrm{\Phi }}_M=\left(kz/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}\right)\partial U/\partial 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, $\mathrm{0}<z/L<\mathrm{10}$, which corresponds to quite small gradient Richardson numbers: $\mathrm{0}<\mathit{Ri}<\mathrm{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\to \mathrm{\infty }$, precisely the same limit as Eq. (7) with $k/R_{\mathrm{\infty }}=\mathrm{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_{\mathrm{\infty }}$ holds true in any stable stratification:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial U}{\partial z}}={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}}{kz}}\left(\mathrm{1}+{\displaystyle \frac{k}{R_{\mathrm{\infty }}}}{\displaystyle \frac{z}{L}}\right)\mathrm{as}\mathrm{0}<z/L<\mathrm{\infty }.\end{array}$$

Then, substituting Eq. (9) for $\partial U/\partial z$ into Eqs. (5) and (6) yields the following simple relations linking Ri_f with z∕L:

$$\begin{array}{}\text{(10)}& {\mathit{Ri}}_{\mathrm{f}}={\displaystyle \frac{kz/L}{\mathrm{1}+k{R}_{\mathrm{\infty }}^{-\mathrm{1}}z/L}},{\displaystyle \frac{z}{L}}={\displaystyle \frac{R_{\mathrm{\infty }}}{k}}{\displaystyle \frac{{\mathit{Ri}}_{\mathrm{f}}}{R_{\mathrm{\infty }}-{\mathit{Ri}}_{\mathrm{f}}}},\end{array}$$

and exact formulation of the TKE dissipation rate as dependent on static stability:

$$\begin{array}{}\text{(11)}& {\mathit{\epsilon }}_{\mathrm{K}}={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{3}/\mathrm{2}}}{kz}}\left[\mathrm{1}+k\left(R_{\mathrm{\infty }}^{-\mathrm{1}}-\mathrm{1}\right){\displaystyle \frac{z}{L}}\right]={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{3}/\mathrm{2}}}{kz}}{\displaystyle \frac{\mathrm{1}-{\mathit{Ri}}_{\mathrm{f}}}{\mathrm{1}-{\mathit{Ri}}_{\mathrm{f}}/R_{\mathrm{\infty }}}}.\end{array}$$

It is worth noting that R_∞ can be derived from well-established 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 cross-validation. For a detailed description of the above numerical models used, see Mortikov (2016), Mortikov et al. (2019), and Druzhinin et al. (2016).

Figure 1 Usual height, z, and vertical coordinate, $\stackrel{\mathrm{̃}}{z}$, defined by Eq. (12) characterizing “absolute” and “internal” geometry of the domain, respectively. Panel (a) shows semi-space, z>0, where $\stackrel{\mathrm{̃}}{z}=z$. Panel (b) shows the layer between two horizontal walls, $\mathrm{0}<z<d$, where $\stackrel{\mathrm{̃}}{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.

In our DNS total (turbulent + molecular) fluxes of momentum, τ, and potential temperature, F_z, are practically equal to the turbulent fluxes elsewhere beyond narrow near-wall 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, $\partial U/\partial z$, and potential temperature, $\partial \mathrm{\Theta }/\partial 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, $\stackrel{\mathrm{̃}}{z}$, dictated by conformal mapping of the Couette-flow domain $\left(\mathrm{0}<z<d\right)$ into the semi-space:

$$\begin{array}{}\text{(12)}& {\displaystyle }\stackrel{\mathrm{̃}}{z}={\displaystyle \frac{d}{\mathit{\pi }}}\sin {\displaystyle \frac{\mathit{\pi }z}{d}}\text{in Couette flow}\left(\mathrm{0}<z<d\right).\end{array}$$

This coordinate reflects the influences of lower and upper walls on the fluid flow.

In semi-space, the “internal geometry” coincides with “absolute geometry”: $\stackrel{\mathrm{̃}}{z}=z$. Thus, the vertical structure of the Couette flow in terms of $\stackrel{\mathrm{̃}}{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: ${\mathit{\epsilon }}_{\mathrm{K}}=\mathit{\nu }\langle (\partial u_i/\partial x_k)(\partial u_i/\partial x_k)\rangle$, where ν is kinematic viscosity, and ${\mathit{\epsilon }}_{\mathrm{K}}=\mathit{\tau }\partial U/\partial z+\mathit{\beta }{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 $-\mathrm{5}/\mathrm{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, $\partial U/\partial 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 Flux Richardson number, Ri_f, in stable stratification versus $\stackrel{\mathrm{̃}}{z}/L$ in Couette flow or versus z∕L in the atmospheric surface layer. 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 (light-grey squares), and acoustic soundings in Kalmykia steppe (light-grey stars). Dark-grey dots belong to the very narrow near-surface layer: $\mathrm{0}<\stackrel{\mathrm{̃}}{z}<\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{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 $\mathrm{0}<\stackrel{\mathrm{̃}}{z}<\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}$, where molecular transports are significant and Eq. (10) is not necessarily relevant.

Figure 2 shows flux Richardson number, ${\mathit{Ri}}_{\mathrm{f}}=\mathit{\beta }{F}_z{\left(\mathit{\tau }\cdot \partial U/\partial z\right)}^{-\mathrm{1}}$, versus dimensionless height, $\stackrel{\mathrm{̃}}{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: $\mathrm{0}<\stackrel{\mathrm{̃}}{z}<\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{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 Dimensionless dissipation rate in stable stratification, ${\mathit{\epsilon }}_{\mathrm{K}}\stackrel{\mathrm{̃}}{z}/{\mathit{\tau }}^{\mathrm{3}/\mathrm{2}}$, versus $\stackrel{\mathrm{̃}}{z}/L$ in Couette flow or ${\mathit{\epsilon }}_{\mathrm{K}}z/{\mathit{\tau }}^{\mathrm{3}/\mathrm{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: $\mathrm{0}<\stackrel{\mathrm{̃}}{z}<\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}$, where molecular transports are significant and Eq. (11) is not relevant.

Figure 4 The energy Richardson number, ${\mathit{Ri}}_{\mathrm{E}}=E_{\mathrm{P}}/E_{\mathrm{K}}$, versus $\stackrel{\mathrm{̃}}{z}/L$ in Couette flow or versus z∕L in the atmospheric surface layer. Empirical data are from the same sources as in Figs. 2 and 3. The solid black line shows Eq. (18) with k=0.4, R_∞=0.2, and C_P=0.62.

Figure 3 shows dimensionless dissipation rate, ${\mathit{\epsilon }}_{\mathrm{K}}z/{\mathit{\tau }}^{\mathrm{3}/\mathrm{2}}$, versus z∕L after Eq. (11) and atmospheric data; and ${\mathit{\epsilon }}_{\mathrm{K}}\stackrel{\mathrm{̃}}{z}/{\mathit{\tau }}^{\mathrm{3}/\mathrm{2}}$ versus $\stackrel{\mathrm{̃}}{z}/L$ after DNS of Couette flow. All notations are the same as in Fig. 2. The theoretical curve plotted after Eq. (11) with k=0.4 and R_∞=0.2 is fully consistent with experimental data, except for the narrow transition layer $\mathrm{0}<\stackrel{\mathrm{̃}}{z}<\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{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: ${\mathit{Ri}}_{\mathrm{f}}=\frac{\mathit{\beta }{F}_z}{\mathit{\tau }\cdot \partial U/\partial 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 showing ${\mathit{\epsilon }}_{\mathrm{K}}z/{\mathit{\tau }}^{\mathrm{3}/\mathrm{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.

3 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_{\mathrm{T}}\equiv E_{\mathrm{K}}/{\mathit{\epsilon }}_{\mathrm{K}}$, and length scale, $l_{\mathrm{T}}\equiv E_{\mathrm{K}}^{\mathrm{1}/\mathrm{2}}{t}_{\mathrm{T}}=E_{\mathrm{K}}^{\mathrm{3}/\mathrm{2}}/{\mathit{\epsilon }}_{\mathrm{K}}$. Then Eq. (11) defines l_T as a function of z∕L:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{l}_{\mathrm{T}}\equiv E_{\mathrm{K}}^{\mathrm{1}/\mathrm{2}}{t}_{\mathrm{T}}={\displaystyle \frac{E_{\mathrm{K}}^{\mathrm{3}/\mathrm{2}}}{{\mathit{\epsilon }}_{\mathrm{K}}}}=kz{\left({\displaystyle \frac{E_{\mathrm{K}}}{\mathit{\tau }}}\right)}^{\mathrm{3}/\mathrm{2}}\\ \text{(13)}& {\displaystyle }& {\displaystyle }{\left[\mathrm{1}+k\left(R_{\mathrm{\infty }}^{-\mathrm{1}}-\mathrm{1}\right){\displaystyle \frac{z}{L}}\right]}^{-\mathrm{1}}.\end{array}$$

It has the asymptotic limits:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{l}_{\mathrm{T}}\to k{\left(E_{\mathrm{K}}/\mathit{\tau }\right)}_{\mathrm{0}}^{\mathrm{3}/\mathrm{2}}z\sim z\mathrm{as}z/L\to \mathrm{0}\\ \text{(14)}& {\displaystyle }& {\displaystyle }\mathrm{and}{l}_{\mathrm{T}}\to {\displaystyle \frac{R_{\mathrm{\infty }}}{\mathrm{1}-R_{\mathrm{\infty }}}}{\left({\displaystyle \frac{E_{\mathrm{K}}}{\mathit{\tau }}}\right)}_{\mathrm{\infty }}^{\mathrm{3}/\mathrm{2}}L\sim L\mathrm{as}z/L\to \mathrm{\infty },\end{array}$$

where the limits of E_K∕τ in neutral stratification and in extremely stable stratification are just dimensionless constants. Our DNS yield the following estimates: ${\left(E_{\mathrm{K}}/\mathit{\tau }\right)}_{\mathrm{0}}\approx \mathrm{4}$ and ${\left(E_{\mathrm{K}}/\mathit{\tau }\right)}_{\mathrm{\infty }}\approx \mathrm{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 different directions is characterized by the mixing-length vector $l_{\mathrm{T}i}\equiv E_{\mathrm{K}i}^{\mathrm{1}/\mathrm{2}}{t}_{\mathrm{T}}$ $\left(i=\mathrm{1},\mathrm{2},\mathrm{3}\right)$ 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 $\stackrel{\mathrm{̃}}{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_{\mathrm{P}}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\beta }\langle {\mathit{\theta }}^{\mathrm{2}}\rangle /(\partial \mathrm{\Theta }/\partial z)$. Hence the effect of stratification on turbulence can be quantified by the “energy Richardson number” defined as

$$\begin{array}{}\text{(15)}& {\mathit{Ri}}_{\mathrm{E}}={\displaystyle \frac{E_{\mathrm{P}}}{E_{\mathrm{K}}}}.\end{array}$$

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_{\mathrm{T}}\equiv E_{\mathrm{K}}^{\mathrm{3}/\mathrm{2}}/{\mathit{\epsilon }}_{\mathrm{K}}$, the budget equations for TKE and TPE become

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{K}}=t_{\mathrm{T}}\left(-\mathit{\tau }\cdot \partial U/\partial z+\mathit{\beta }{F}_z\right),\text{(17)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{P}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\beta }\langle {\mathit{\theta }}^{\mathrm{2}}\rangle }{\partial \mathrm{\Theta }/\partial z}}=-C_{\mathrm{P}}{t}_{\mathrm{T}}\mathit{\beta }{F}_z,\end{array}$$

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:

$$\begin{array}{}\text{(18)}& {\mathit{Ri}}_{\mathrm{E}}={\displaystyle \frac{C_{\mathrm{P}}}{{\mathit{Ri}}_{\mathrm{f}}^{-\mathrm{1}}-\mathrm{1}}}={\displaystyle \frac{C_{\mathrm{P}}kz/L}{\mathrm{1}+\left(R_{\mathrm{\infty }}^{-\mathrm{1}}-\mathrm{1}\right)kz/L}}.\end{array}$$

Figure 4 shows Ri_E versus z∕L or $\stackrel{\mathrm{̃}}{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:

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathit{Ri}}_{\mathrm{E}}\to R_{E\mathrm{\infty }}={\displaystyle \frac{C_{\mathrm{P}}}{R_{\mathrm{\infty }}^{-\mathrm{1}}-\mathrm{1}}}=\mathrm{0.155}\mathrm{as}z/L\to \mathrm{\infty },\end{array}$$

Then, using Eq. (18) to express z∕L through Ri_E, Eq. (11) in terms of Ri_E becomes:

$$\begin{array}{}\text{(20)}& {\mathit{\epsilon }}_{\mathrm{K}}={\mathit{\epsilon }}_{\mathrm{K}\left(\mathrm{neutral}\right)}{\left(\mathrm{1}-{\displaystyle \frac{{\mathit{Ri}}_{\mathrm{E}}}{R_{E\mathrm{\infty }}}}\right)}^{-\mathrm{1}},\end{array}$$

where ε_K(neutral) is dissipation rate in neutral stratification. In the surface layer ${\mathit{\epsilon }}_{\mathrm{K}\left(\mathrm{neutral}\right)}={\mathit{\tau }}^{\mathrm{3}/\mathrm{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: $\mathit{\tau }=-K_{\mathrm{M}}\partial U/\partial z$ and $F_z=-K_{\mathrm{H}}\partial \mathrm{\Theta }/\partial z$, where K_M is eddy viscosity and K_H is eddy conductivity. Then, finite-difference 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.

4 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, $\partial U/\partial z$, and also via the Kolmogorov $-\mathrm{5}/\mathrm{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\to \mathrm{\infty }$ where L is the Obukhov length scale.

By providential coincidence, the formulations happen to be precisely the same in the asymptotic limit $z/L\to \mathrm{\infty }$ and in the weakly stable stratifications $\mathrm{0}<z/L<\mathrm{10}$ typical of atmospheric surface layer. This yields simple analytical formulations of the dimensionless velocity gradient, $\left(kz/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}\right)\partial U/\partial z$, and dissipation rate, $\left(kz/{\mathit{\tau }}^{\mathrm{3}/\mathrm{2}}\right){\mathit{\epsilon }}_{\mathrm{K}}$, as universal functions of z∕L (Eqs. 9 and 11) across the whole range of stratifications from neutral to extremely stable.

Universal analytical formulation of $\left(kz/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}\right)\partial U/\partial z$ versus z∕L yields the single-valued relations linking z∕L as the criterion of stratification in the surface-layer flow or $\stackrel{\mathrm{̃}}{z}/L$ as the same criterion in Couette flow with alternative criterions: flux Richardson number, Ri_f (Eq. 5), and the newly introduced “energy Richardson number”, Ri_E (Eq. 13), applicable to any turbulent regimes. This opens up prospects for extending the obtained dependence of dissipation rate on static stability to any stably stratified turbulent flows.