On dissipation timescales of the basic second-order moments: the effect on the energy and flux budget (EFB) turbulence closure for stably stratified turbulence

The dissipation rates of the basic second-order moments are the key parameters playing a vital role in turbulence modelling and controlling turbulence energetics and spectra and turbulent fluxes of momentum and heat. In this paper, we use the results of direct numerical simulations (DNSs) to evaluate dissipation rates of the basic second-order moments and revise the energy and flux budget (EFB) turbulence closure theory for stably stratified turbulence. We delve into the theoretical implications of this approach and substantiate our closure hypotheses through DNS data. We also show why the concept of down-gradient turbulent transport becomes incomplete when applied to the vertical turbulent flux of potential temperature under stable stratification. We reveal essential feedback between the turbulent kinetic energy (TKE), the vertical turbulent flux of buoyancy, and the turbulent potential energy (TPE), which is responsible for maintaining shear-produced stably stratified turbulence for any Richardson number.

1 Introduction

Turbulence and associated turbulent transport have been studied theoretically, experimentally, observationally, and numerically throughout several decades (see books by Batchelor, 1953; Monin and Yaglom, 1971, 2013; Tennekes and Lumley, 1972; Frisch, 1995; Pope, 2000; Davidson, 2013; Rogachevskii, 2021, and references therein), but some important questions remain. This is particularly true in applications to atmospheric physics and geophysics, where Reynolds and Péclet numbers are extremely large, so the governing equations are strongly nonlinear. The classical Kolmogorov theory (Kolmogorov, 1941a, b, 1942, 1991) has been formulated for neutrally stratified homogeneous and isotropic turbulence.

In atmospheric boundary layers, temperature stratification causes turbulence to become anisotropic and inhomogeneous, making some assumptions underlying Kolmogorov's theory questionable. Numerous alternative turbulence closure theories (see reviews by Weng and Taylor, 2003; Umlauf and Burchard, 2005; Mahrt, 2014) have been formulated using the budget equations for not only turbulent kinetic energy (TKE), but also turbulent potential energy (TPE) (see, for example, Holloway, 1986; Ostrovsky and Troitskaya, 1987; Dalaudier and Sidi, 1987; Hunt et al., 1988; Canuto and Minotti, 1993; Schumann and Gerz, 1995; Hanazaki and Hunt, 1996; Keller and van Atta, 2000; Canuto et al., 2001; Stretch et al., 2001; Cheng et al., 2002; Hanazaki and Hunt, 2004; Rehmann and Hwang, 2005; Umlauf, 2005). The budget equations for all three energies, TKE, TPE, and total turbulent energy (TTE), were considered by Canuto and Minotti (1993), Elperin et al. (2002, 2006), Zilitinkevich et al. (2007), and Canuto et al. (2008).

The energy and flux budget (EFB) turbulence closure theory, which is based on the budget equations for the densities of TKE, TPE, and turbulent fluxes of momentum and heat, has been developed for the stably stratified atmospheric flows (Zilitinkevich et al., 2007, 2008, 2009, 2013, 2019; Kleeorin et al., 2019), surface layers in atmospheric convective turbulence (Rogachevskii et al., 2022), and core of the convective boundary layer (Rogachevskii and Kleeorin, 2024), as well as for passive scalar transport (Kleeorin et al., 2021). The EFB closure theory has shown that strong atmospheric stably stratified turbulence is maintained by large-scale shear (mean wind) for any stratification, and the “critical” Richardson number, considered for many years to be a threshold between the turbulent and laminar states of the flow, actually separates two turbulent regimes: the strong turbulence typical of atmospheric boundary layers and the weak three-dimensional turbulence typical of the free atmosphere and characterised by a strong decrease in the turbulent heat transfer in comparison to the momentum transfer.

Some other turbulent closure models (Mauritsen et al., 2007; Canuto et al., 2008; Sukoriansky and Galperin, 2008; Li et al., 2016) do not employ the concept of the critical Richardson number, so shear-generated turbulent mixing may persist for any stratification. In particular, Mauritsen et al. (2007) have developed a turbulent closure based on the budget equation for TTE (instead of TKE) and different observational findings to take into account the mean flow stability. They used this turbulent closure model to study the turbulent transfer of heat and momentum under very stable stratification. In their model, although the turbulent heat flux tends toward zero beyond a certain stability limit, the turbulent stress stays finite. However, the model by Mauritsen et al. (2007) does not use the budget equation for TPE and the vertical turbulent heat flux.

L'vov et al. (2008) performed detailed analyses of the budget equations for the Reynolds stresses in the turbulent boundary layer (relevant to the strong turbulence regime), explicitly taking into consideration the dissipative effect in the horizontal turbulent heat flux budget equation in contrast to the EFB “effective dissipation approximation”, adopted in the EFB turbulent closure model. However, the theory by L'vov et al. (2008) still contains the critical gradient Richardson number for the existence of the shear-produced turbulence.

Sukoriansky and Galperin (2008) apply a quasi-normal-scale elimination theory that is similar to the renormalization group analysis. Sukoriansky and Galperin (2008) do not use the budget equations for TKE, TPE, and TTE in their analysis. This theory correctly describes the dependence of the turbulent Prandtl number versus the gradient Richardson number and does not employ the concept of the critical gradient Richardson number for the existence of turbulence. However, this approach does not have detailed Richardson number dependences of the other non-dimensional parameters, such as the ratio between TPE and TTE, dimensionless turbulent flux of momentum, or dimensionless vertical turbulent flux of potential temperature. Their background non-stratified shear-produced turbulence is assumed to be isotropic and homogeneous. Canuto et al. (2008) have generalised their original model (see Cheng et al., 2002), introducing the new parameterisation for the buoyancy timescale to accommodate the existence of stably stratified shear-produced turbulence at arbitrary Richardson numbers.

Li et al. (2016) have developed the co-spectral budget (CSB) closure approach, which is formulated in the Fourier space and integrated across all turbulent scales to obtain turbulent characteristics in physical space. The CSB model allows turbulence to exist at any gradient Richardson number. However, the CSB model yields different predictions for the vertical anisotropy versus the Richardson number compared to the EFB theory. All state-of-the-art turbulent closures follow the so-called Kolmogorov hypothesis: all dissipation timescales of turbulent second-order moments are assumed to be proportional to each other, which, at first glance, looks reasonable but is, in fact, hypothetical for stably stratified turbulence.

The present study aims to demonstrate the dependence of dissipation timescales of basic second-order moments on stability through direct numerical simulation (DNS) experiments. The obtained numerical results allow us to modify the EFB turbulence closure theory to account for that dependency. It is worth noting that the DNSs presented here are limited to bulk Richardson numbers (based on the wall velocity and temperature differences and channel height) up to Ri_b=0.11 and Reynolds numbers (based on the wall velocity difference and channel height; see Sect. 3) up to Re=120 000.

This paper is organised as follows. In Sect. 2, we formulate basic budget equations and main assumptions in the framework of the EFB turbulence closure theory. Section 3 describes the setup for DNSs of stably stratified turbulent plane Couette flow to determine the vertical profiles of the dissipation timescales of turbulent second-order moments. In Sect. 4, we formulate the modified EFB turbulence closure theory while considering the dependencies of the dissipation timescales of basic second-order moments on the gradient Richardson number obtained from DNSs. There, we also perform the validation of the modified EFB turbulence closure model, which yields vertical profiles of the basic turbulence parameters (including the turbulent Prandtl number, the ratio of TPE to TKE, and the normalised turbulent heat flux) using the data from the DNS. Finally, in Sect. 5, we discuss the obtained results and draw the conclusions.

2 Problem setting and basic equations

We consider plane-parallel, stably stratified dry-air flow and employ the familiar budget equations underlying turbulence closure theory (e.g. Kleeorin et al., 2021; Zilitinkevich et al., 2013, 2019; Kaimal and Finnigan, 1994; Canuto et al., 2008) for the Reynolds stress, τ_ij=〈u_iu_j〉; the turbulent flux of potential temperature, F_i=〈θu_i〉; and the intensity of potential temperature fluctuations, $E_{\mathit{\theta }}=\langle {\mathit{\theta }}^{\mathrm{2}}\rangle /\mathrm{2}$:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{D}{\mathit{\tau }}_{ij}}{\mathrm{D}t}}+{\displaystyle \frac{\partial }{\partial z}}{\mathrm{\Phi }}_{ij\mathrm{3}}^{\left(\mathit{\tau }\right)}& =-{\mathit{\tau }}_{i\mathrm{3}}{\displaystyle \frac{\partial U_j}{\partial z}}-{\mathit{\tau }}_{j\mathrm{3}}{\displaystyle \frac{\partial U_i}{\partial z}}\\ & -\left[{\mathit{\epsilon }}_{ij}^{\left(\mathit{\tau }\right)}-\mathit{\beta }\left(F_j{\mathit{\delta }}_{i\mathrm{3}}+F_i{\mathit{\delta }}_{j\mathrm{3}}\right)-Q_{ij}\right],\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{D}{F}_i}{\mathrm{D}t}}+{\displaystyle \frac{\partial }{\partial z}}{\mathrm{\Phi }}_i^{\left(F\right)}& =\mathit{\beta }{\mathit{\delta }}_{i\mathrm{3}}〈{\mathit{\theta }}^{\mathrm{2}}〉-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}〈\mathit{\theta }{\displaystyle \frac{\partial p}{\partial x_i}}〉\\ & -{\mathit{\tau }}_{i\mathrm{3}}{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}-F_z{\displaystyle \frac{\partial U_i}{\partial z}}-{\mathit{\epsilon }}_i^{\left(F\right)},\end{array}\text{(3)}& {\displaystyle \frac{\mathrm{D}{E}_{\mathit{\theta }}}{\mathrm{D}t}}+{\displaystyle \frac{\partial }{\partial z}}{\mathrm{\Phi }}^{\left(\mathit{\theta }\right)}=-F_z{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}-{\mathit{\epsilon }}_{\mathit{\theta }}.\end{array}$$

Here, x₁=x and x₂=y are horizontal coordinates; x₃=z is the vertical coordinate; t is the time; $\mathbf{U}=\left(U_{\mathrm{1}},U_{\mathrm{2}},U_{\mathrm{3}}\right)=\left(U,V,W\right)$ is the mean flow velocity; $\mathbf{u}=\left(u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}}\right)=\left(u,v,w\right)$ is the velocity fluctuations; $\mathrm{\Theta }=T{\left(P_{\mathrm{0}}/P\right)}^{\mathrm{1}-\mathrm{1}/\mathit{\gamma }}$ is the mean potential temperature (expressed through absolute temperature, T, and pressure, P); T₀, P₀, and ρ₀ are reference values of temperature, pressure, and density, respectively; $\mathit{\gamma }=c_p/c_v=\mathrm{1.41}$ is the ratio of specific heats; θ and p are fluctuations of potential temperature and pressure; $\mathrm{D}/\mathrm{D}t=\partial /\partial t+U_k\partial /\partial x_k$ is the advective derivative; angle brackets denote the averaging; $\mathit{\beta }=g/T_{\mathrm{0}}$ is the buoyancy parameter; g is the acceleration due to gravity; δ_ij is the unit tensor (δ_ij=1 for i=j and δ_ij=0 for i≠j); and ${\mathrm{\Phi }}_{ij\mathrm{3}}^{\left(\mathit{\tau }\right)}$, ${\mathrm{\Phi }}_i^{\left(F\right)}$, and Φ^(θ) are the third-order moments, which describe turbulent transport of the second-order moments under consideration.

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}{\mathrm{\Phi }}_{ij\mathrm{3}}^{\left(\mathit{\tau }\right)}& =〈u_i{u}_{j}w〉+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left(〈p{u}_i〉{\mathit{\delta }}_{j\mathrm{3}}+〈p{u}_j〉{\mathit{\delta }}_{i\mathrm{3}}\right)\\ & -\mathit{\nu }\left(〈u_i{\displaystyle \frac{\partial u_j}{\partial z}}〉+〈u_j{\displaystyle \frac{\partial u_i}{\partial z}}〉\right),\end{array}\text{(5)}& {\displaystyle }{\mathrm{\Phi }}_i^{\left(F\right)}=〈u_{i}w\mathit{\theta }〉-\mathit{\nu }〈\mathit{\theta }{\displaystyle \frac{\partial u_i}{\partial z}}〉-\mathit{\kappa }〈u_i{\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}〉,\text{(6)}& {\displaystyle }{\mathrm{\Phi }}^{\left(\mathit{\theta }\right)}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}〈{\mathit{\theta }}^{\mathrm{2}}w〉-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{2}}}{\displaystyle \frac{\partial }{\partial z}}〈{\mathit{\theta }}^{\mathrm{2}}〉,\end{array}$$

where Q_ij is the correlation between fluctuations of pressure and strain-rate tensor, which controls the interactions between the Reynolds stress components as follows:

$$\begin{array}{}\text{(7)}& Q_{ij}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}〈p\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)〉.\end{array}$$

Here, ${\mathit{\epsilon }}_{ij}^{\left(\mathit{\tau }\right)}$, ${\mathit{\epsilon }}_i^{\left(F\right)}$, and ε_θ are the dissipation rates of the second-order moments:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\epsilon }}_{ij}^{\left(\mathit{\tau }\right)}=\mathrm{2}\mathit{\nu }〈{\displaystyle \frac{\partial u_i}{\partial z}}{\displaystyle \frac{\partial u_j}{\partial z}}〉,\text{(9)}& {\displaystyle }{\mathit{\epsilon }}_i^{\left(F\right)}=\left(\mathit{\nu }+\mathit{\kappa }\right)〈{\displaystyle \frac{\partial u_i}{\partial z}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}〉,\text{(10)}& {\displaystyle }{\mathit{\epsilon }}_{\mathit{\theta }}=\mathit{\kappa }〈{\left({\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}\right)}^{\mathrm{2}}〉,\end{array}$$

where ν is kinematic viscosity and κ is thermal conductivity.

The budget of TKE components, $E_i=〈u_i^{\mathrm{2}}〉/\mathrm{2}$ ($i=\mathrm{1},\mathrm{2},\mathrm{3}$), is determined by Eq. (1) for i=j, which yields the familiar TKE budget equation:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{D}{E}_K}{\mathrm{D}t}}& +{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}〈u_i^{\mathrm{2}}w〉+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}〈pw〉-{\displaystyle \frac{\mathit{\nu }}{\mathrm{2}}}{\displaystyle \frac{\partial 〈u_i^{\mathrm{2}}〉}{\partial z}}\right)\\ & =-\mathit{\tau }\cdot {\displaystyle \frac{\partial \mathbf{U}}{\partial z}}+\mathit{\beta }{F}_z-{\mathit{\epsilon }}_K,\end{array}\end{array}$$

where $E_K=\sum E_i$ is TKE and ${\mathit{\epsilon }}_K=\sum {\mathit{\epsilon }}_{ii}^{\left(\mathit{\tau }\right)}/\mathrm{2}$ is the TKE dissipation rate. The sum of the term Q_ii (the trace of the tensor Q_ij) is equal to zero because of the incompressibility, constraint on the flow velocity field, $\partial u_i/\partial x_i=\mathrm{0}$; that is, Q_ij only redistributes energy between TKE components.

Likewise, ε_θ is the dissipation rate of the intensity of potential temperature fluctuations, E_θ, and ${\mathit{\epsilon }}_i^{\left(F\right)}$ are the dissipation rates of the three components of the turbulent flux of potential temperature, F_i.

Following Kolmogorov (1941b, 1942), the dissipation rates, ε_K and ε_θ, are taken to be proportional to the dissipating quantities divided by corresponding timescales:

$$\begin{array}{}\text{(12)}& {\mathit{\epsilon }}_K={\displaystyle \frac{E_K}{t_K}},{\mathit{\epsilon }}_{\mathit{\theta }}={\displaystyle \frac{E_{\mathit{\theta }}}{t_{\mathit{\theta }}}},\end{array}$$

where t_K is the TKE dissipation timescale and t_θ is the dissipation timescale of E_θ. Here, the formulation of the dissipation rates is not hypothetical: it merely expresses one unknown (dissipation rate) through another (dissipation timescale).

Figure 1 Profiles of mean flow velocity and mean potential temperature in stably stratified turbulent plane Couette flow. Light-grey dots belong to the viscous sublayer.

In this study, we consider the EFB theory in its simplest, algebraic form, neglecting non-steady terms in all budget equations and neglecting divergence of the fluxes of TKE and TPE (determined by third-order moments). This approach is reasonable because, for example, the characteristic times of variations of the second moments are much larger than the turbulent timescales for large Reynolds and Péclet numbers. We also assume that the terms related to the divergence of the fluxes of TKE and TPE for stably stratified turbulence are much smaller than the rates of production and dissipation in budget Eqs. (3) and (11). In this case, the TKE budget equation, Eq. (11), and the budget equation for E_θ, Eq. (3), become

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathrm{0}=-\mathit{\tau }{\displaystyle \frac{\partial U}{\partial z}}+\mathit{\beta }{F}_z-{\mathit{\epsilon }}_K,\text{(14)}& {\displaystyle }\mathrm{0}=-F_z{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}-{\mathit{\epsilon }}_{\mathit{\theta }}.\end{array}$$

The intensity of the potential temperature fluctuations, E_θ, determines TPE as follows:

$$\begin{array}{}\text{(15)}& E_P={\displaystyle \frac{\mathit{\beta }{E}_{\mathit{\theta }}}{\partial \mathrm{\Theta }/\partial z}},\end{array}$$

so that Eq. (14) becomes

$$\begin{array}{}\text{(16)}& \mathrm{0}=-\mathit{\beta }{F}_z-{\mathit{\epsilon }}_P,\end{array}$$

where ${\mathit{\epsilon }}_P=E_P/t_{\mathit{\theta }}$ is the TPE dissipation time.

The first term on the right-hand side of Eq. (13), $-\mathit{\tau }\partial U/\partial z$, is the rate of the TKE production, while the second term, βF_z, is the buoyancy, which, in stably stratified flow, causes decay of TKE; that is, it results in the conversion of TKE into TPE. The ratio of these terms is the flux Richardson number:

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

and this dimensionless parameter characterises the effect of stratification on turbulence.

Taking into account Eq. (17), the steady-state versions of TKE and TPE budget equations, Eqs. (13) and (14), can be rewritten as

$$\begin{array}{}\text{(18)}& {\displaystyle }{E}_K=\mathit{\tau }{\displaystyle \frac{\partial U}{\partial z}}\left(\mathrm{1}-{\mathit{Ri}}_{\mathrm{f}}\right)t_K,\text{(19)}& {\displaystyle }{E}_P=\mathit{\tau }{\displaystyle \frac{\partial U}{\partial z}}{\mathit{Ri}}_{\mathrm{f}}{t}_{\mathit{\theta }}.\end{array}$$

Thus, the ratio of TPE to TKE is

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{E_P}{E_K}}={\displaystyle \frac{{\mathit{Ri}}_{\mathrm{f}}}{\mathrm{1}-{\mathit{Ri}}_{\mathrm{f}}}}{\displaystyle \frac{t_{\mathit{\theta }}}{t_K}}.\end{array}$$

Zilitinkevich et al. (2013) suggested the following relation, linking Ri_f with another stratification parameter, $z/L$:

$$\begin{array}{}\text{(21)}& {\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}$$

where $L=-{\mathit{\tau }}^{\mathrm{3}/\mathrm{2}}/\mathit{\beta }{F}_z$ is the Obukhov length scale, k=0.4 is the von Kármán constant, and R_∞=0.2 is the maximum value of the flux Richardson number.

On the right-hand side of Eq. (20), there is an unknown ratio between two dissipation timescales, $t_{\mathit{\theta }}/t_K$. The Kolmogorov hypothesis suggests that it is a universal constant. We do not imply that this assumption applies but instead investigate a possible stability dependency of dissipation timescale ratios and improve the EFB turbulence closure model accounting for it. To this end, we perform DNSs of stably stratified turbulent plane Couette flow (see Sect. 3) to measure the dissipation timescales of basic second-order moments and validate the modified EFB turbulence closure model (see Sect. 4).

3 Methods

For our study, we conducted a series of direct numerical simulations of stably stratified turbulent plane Couette flow. This flow occurs between two parallel plates that move relative to each other, producing shear and turbulence, with the plates at different temperatures, thus creating stable stratification. In Couette flow, the total (turbulent plus molecular) vertical fluxes of momentum and potential temperature remain constant, independent of distance from the walls, which, in particular, ensures a very certain fixed value of the Obukhov length scale. Figure 1 illustrates the profiles of mean flow velocity and mean potential temperature. We recall that all our derivations are relevant to the well-developed turbulence regime where molecular transports are negligible compared to turbulent transports, so turbulent fluxes practically coincide with total fluxes. This is the case in our DNS, except for the narrow near-wall viscous turbulent-flow transition layers. Data from these layers, obviously irrelevant to the turbulence regime we consider, are shown by grey points in the figures and are ignored in fitting procedures. In further analysis, we primarily utilise $z/L$ as a stratification parameter instead of Ri or Ri_f because it offers a better dynamic range in our experiments. While Ri remains practically constant in each DNS run and Ri_f is limited in its growth, the parameter $z/L$ is determined by the distance from the walls, thus varying significantly in every DNS run.

Numerical simulation of stably stratified turbulent Couette flow was performed using the unified DNS, LES, and Reynolds-averaged Navier–Stokes equations (RANS) code developed at the Moscow State University (MSU) and the Institute of Numerical Mathematics (INM) of the Russian Academy of Science (see Mortikov, 2016; Mortikov et al., 2019; Bhattacharjee et al., 2022; Debolskiy et al., 2023; Gladskikh et al., 2023; Zasko et al., 2023). The code is designed for high-resolution simulations on modern-day high-performance computing (HPC) systems. The DNS part of the code solves the finite-difference approximation of the incompressible Navier–Stokes system of equations under the Boussinesq approximation. Conservative schemes on the staggered grid (Morinishi et al., 1998; Vasilyev, 2000) of fourth-order accuracy are used in the horizontal direction, while, in the vertical direction, the spatial approximation is restricted to second-order accuracy with near-wall grid resolution refinement sufficient to resolve the near-wall viscous region. The time step used in the simulations was determined by Courant–Friedrichs–Lewy (CFL) restrictions, with CFL maintained at approximately 0.1 in all runs. This corresponds to a value of $u_{\ast }^{\mathrm{2}}\mathrm{\Delta }t/\mathit{\nu }$ on the order of 0.01. The projection method (Brown et al., 2001) is used for the time advancement of momentum equations coupled with the incompressibility condition, while the multigrid method is applied to solve the Poisson equation to ensure that the velocity is divergence-free at each time step. For the Couette flow, periodic boundary conditions are used in the horizontal directions, and no-slip and no-penetration conditions are set on the channel walls for the velocity. The stable stratification is maintained by prescribed Dirichlet boundary conditions on the potential temperature. In all experiments, the value of molecular Prandtl number (ratio of kinematic viscosity and thermal diffusivity of the fluid) was fixed at 0.7 based on its typical value for air (Monin and Yaglom, 1971). The simulations were performed for a wide range of the Reynolds number, Re, defined by the wall velocity difference, channel height, and kinematic viscosity, that is from 40 000 up to 120 000 (see Table 1). All experiments were carried out using the resources of MSU and Center for Scientific Computing (CSC) HPC facilities. For the maximum Re values achieved, the numerical grid consisted of more than 2×10⁸ cells and the calculations used about 10 000 CPU cores.

Table 1 Overview of DNS experiments and key parameters.

For each Reynolds number, we conducted a series of experiments. Beginning with neutral conditions (no imposed gradient of the mean potential temperature), we incrementally increased the bulk Richardson number, which characterises the stable stratification, in each successive experiment. By gradually increasing the stability in each experiment, we were able to cover a wide range of Ri values, extending from neutral to stably stratified states. In each run, the turbulent flow was allowed sufficient time to develop and reach statistical steady-state conditions, which required a spin-up period of at least 15 $H/u_{\ast }$ periods. This ensured that parameters such as the total momentum flux remained constant and that the TKE balance was in a steady state. The fully developed steady state was used as the initial conditions for the higher Ri or Re experiment setups. Additionally, all terms in the second-order moment budget equations (Eqs. 1–3) were evaluated consistently using the finite-difference approximation used, resulting in a negligible residual. This approach enabled us to comprehensively study the characteristics of shear-produced stably stratified turbulence, explicitly resolving all dissipation timescales of turbulent second-order moments.

Figure 2 The ratio of the effective dissipation timescale of τ and the dissipation timescale of TKE $t_{\mathit{\tau }}/t_K$ versus $z/L$. The data used for the calibration are obtained in DNS experiments employing the MSU–INM unified code. Only every sixth data point is presented to increase visibility. For the full dataset, please see Kadantsev and Mortikov (2024). The near-surface layer, essentially affected by molecular viscosity ($\mathrm{0}<z<\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}$), is excluded from the analysis. This sublayer is represented by the dotted light-grey lines. The solid black line shows Eq. (27), with empirical constants $C_{\mathrm{1}}^{\mathit{\tau }K}=\mathrm{0.08}$, $C_{\mathrm{2}}^{\mathit{\tau }K}=\mathrm{0.4}$, and $C_{\mathrm{3}}^{\mathit{\tau }K}=\mathrm{2}$, obtained from the best fit of Eq. (27) to DNS data in the turbulent layer, $z>\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}$.

4 Modified EFB closure model for the steady-state regime of turbulence

In the steady-state, Eq. (1) for the vertical component of the turbulent flux of momentum, τ, becomes

$$\begin{array}{}\text{(22)}& \mathrm{0}=-\mathrm{2}{E}_z{\displaystyle \frac{\partial U}{\partial z}}-\left[{\mathit{\epsilon }}_{\mathit{\tau }}-\mathit{\beta }{F}_x-Q_{\mathrm{13}}\right].\end{array}$$

Following Zilitinkevich et al. (2007, 2013, 2019) we define the sum of all terms in square brackets on the right-hand side of Eq. (22) as the “effective dissipation”:

$$\begin{array}{}\text{(23)}& {\mathit{\epsilon }}_{\mathit{\tau }}^{\left(\mathrm{eff}\right)}={\mathit{\epsilon }}_{\mathit{\tau }}-\mathit{\beta }{F}_x-Q_{\mathrm{13}}\equiv {\displaystyle \frac{\mathit{\tau }}{t_{\mathit{\tau }}}}.\end{array}$$

Thus, Eq. (22) becomes

$$\begin{array}{}\text{(24)}& \mathrm{0}=-\mathrm{2}{E}_z{\displaystyle \frac{\partial U}{\partial z}}-{\displaystyle \frac{\mathit{\tau }}{t_{\mathit{\tau }}}},\end{array}$$

yielding the well-known down-gradient formulation of the vertical turbulent flux of momentum:

$$\begin{array}{}\text{(25)}& \mathit{\tau }=-K_M{\displaystyle \frac{\partial U}{\partial z}},K_M=\mathrm{2}{A}_z{E}_K{t}_{\mathit{\tau }},\end{array}$$

where $A_z\equiv E_z/E_K$ is the vertical share of TKE (the vertical anisotropy parameter). By substituting Eq. (25) into Eq. (18), we obtain

$$\begin{array}{}\text{(26)}& {\left({\displaystyle \frac{\mathit{\tau }}{E_K}}\right)}^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}{A}_z}{\mathrm{1}-{\mathit{Ri}}_{\mathrm{f}}}}{\displaystyle \frac{t_{\mathit{\tau }}}{t_K}}.\end{array}$$

In Eq. (26), all the variables are exactly resolved numerically in DNSs, making a detailed investigation on $t_{\mathit{\tau }}/t_K$ possible. Figure 2 demonstrates that the dissipation timescale ratio, $t_{\mathit{\tau }}/t_K$, is a function of the stratification parameter $z/L$ rather than a constant. We propose to approximate this function with a ratio of two first-order polynomials:

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{t_{\mathit{\tau }}}{t_K}}={\displaystyle \frac{C_{\mathrm{1}}^{\mathit{\tau }K}z/L+C_{\mathrm{2}}^{\mathit{\tau }K}}{z/L+C_{\mathrm{3}}^{\mathit{\tau }K}}}.\end{array}$$

Here, the dimensionless empirical constants are obtained from the best fit of Eq. (27) to DNS bin-averaged data: $C_{\mathrm{1}}^{\mathit{\tau }K}=\mathrm{0.08}$, $C_{\mathrm{2}}^{\mathit{\tau }K}=\mathrm{0.4}$, and $C_{\mathrm{3}}^{\mathit{\tau }K}=\mathrm{2}$. The fitting is done using the rational regression model of Curve Fitting Toolbox version 3.5.13 (R2021a) (2022). The ratio of two first-order polynomials is chosen as a simpler fitting function that could provide monotonicity, reasonable smoothness, and clear asymptotes. The only three adjustable parameters of this approximation correspond to the function value at $z/L=\mathrm{0}$, the $z/L\to \mathrm{\infty }$ limit, and the transition between them.

Figure 3 Comparison of two terms, $\mathit{\beta }〈{\mathit{\theta }}^{\mathrm{2}}〉/\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}〈\mathit{\theta }\frac{\partial p}{\partial z}〉$ and $\mathrm{2}{E}_z\frac{\partial \mathrm{\Theta }}{\partial z}/\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}〈\mathit{\theta }\frac{\partial p}{\partial z}〉$, after the same DNS for stably stratified Couette flow. The solid black line represents the linear dependency of the latter on the former, which turns into Eq. (31) after multiplication by $\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}〈\mathit{\theta }\frac{\partial p}{\partial z}〉$ and simple recombination. The fitting coefficients are C_θ=0.82 and $C_{\mathrm{\nabla }}=-\mathrm{0.80}$.

Proceeding to the vertical flux of potential temperature, F_z, we derive its steady-state budget equation from Eq. (2):

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{\partial }{\partial z}}{\mathrm{\Phi }}_z^{\left(F\right)}=\mathit{\beta }〈{\mathit{\theta }}^{\mathrm{2}}〉-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}〈\mathit{\theta }{\displaystyle \frac{\partial p}{\partial z}}〉-\mathrm{2}{E}_z{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}-{\mathit{\epsilon }}_F.\end{array}$$

DNS modelling has shown that $\frac{\partial }{\partial z}{\mathrm{\Phi }}_z^{\left(F\right)}$ term is of the same order of magnitude as ε_F and of the same sign, so we introduce the “effective dissipation rate”, ${\mathit{\epsilon }}_F^{\left(\mathrm{eff}\right)}$:

$$\begin{array}{}\text{(29)}& {\mathit{\epsilon }}_F^{\left(\mathrm{eff}\right)}={\mathit{\epsilon }}_F+{\displaystyle \frac{\partial }{\partial z}}{\mathrm{\Phi }}_z^{\left(F\right)}\equiv {\displaystyle \frac{F_z}{t_F}}.\end{array}$$

Consequently, Eq. (28) reduces to

$$\begin{array}{}\text{(30)}& \mathrm{0}=\mathit{\beta }〈{\mathit{\theta }}^{\mathrm{2}}〉-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}〈\mathit{\theta }{\displaystyle \frac{\partial p}{\partial z}}〉-\mathrm{2}{E}_z{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}-{\displaystyle \frac{F_z}{t_F}}.\end{array}$$

Traditionally, the pressure term was either assumed to be negligible or declared to be proportional to the β〈θ²〉 term (see Zilitinkevich et al., 2007, 2013). However, our DNS data have shown that it is neither negligible nor proportional to any other term in the budget equation, Eq. (30). Instead, we found that it is well approximated by a linear combination of the production and transport terms of Eq. (30) (see Fig. 3):

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}〈\mathit{\theta }{\displaystyle \frac{\partial p}{\partial z}}〉=C_{\mathit{\theta }}\mathit{\beta }〈{\mathit{\theta }}^{\mathrm{2}}〉+C_{\mathrm{\nabla }}\mathrm{2}{E}_z{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}.\end{array}$$

The dimensionless constants C_θ=0.82 and $C_{\mathrm{\nabla }}=-\mathrm{0.80}$ are obtained from the best fit of Eq. (31) to DNS data.

Figure 4 The ratio of two terms from the square brackets of Eq. (34) versus $z/L$. The data are the same as in Fig. 2. The solid black line shows Eq. (36) with the empirical constant C_Pr=0.65 obtained from the best fit of Eq. (34) to DNS data in the turbulent layer, $z>\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}$.

By substituting Eq. (31) into Eq. (30), we rewrite the budget equation as

$$\begin{array}{}\text{(32)}& \mathrm{0}=\left(\mathrm{1}-C_{\mathit{\theta }}\right)\mathit{\beta }〈{\mathit{\theta }}^{\mathrm{2}}〉-\left(\mathrm{1}+C_{\mathrm{\nabla }}\right)\mathrm{2}{E}_z{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}-{\displaystyle \frac{F_z}{t_F}}.\end{array}$$

Substituting Eq. (15) for 〈θ²〉 into Eq. (32) allows us to express F_z by a familiar temperature gradient expression:

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{F}_z& =-K_H{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}\\ K_H& =\left[\left(\mathrm{1}+C_{\mathrm{\nabla }}\right)-\left(\mathrm{1}-C_{\mathit{\theta }}\right){\displaystyle \frac{E_P}{A_z{E}_K}}\right]\mathrm{2}{A}_z{E}_K{t}_F.\end{array}\end{array}$$

Substituting Eq. (33) into Eq. (14) gives

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{F_z^{\mathrm{2}}}{E_{\mathit{\theta }}{E}_K}}=\mathrm{2}\left[\left(\mathrm{1}+C_{\mathrm{\nabla }}\right)A_z-\left(\mathrm{1}-C_{\mathit{\theta }}\right){\displaystyle \frac{E_P}{E_K}}\right]{\displaystyle \frac{t_F}{t_{\mathit{\theta }}}}.\end{array}$$

Next, the turbulent Prandtl number, defined as ${\mathit{Pr}}_T=K_M/K_H$, is given by

$$\begin{array}{}\text{(35)}& {\mathit{Pr}}_T={\displaystyle \frac{t_{\mathit{\tau }}}{t_F}}/\left[\left(\mathrm{1}+C_{\mathrm{\nabla }}\right)-\left(\mathrm{1}-C_{\mathit{\theta }}\right){\displaystyle \frac{E_P}{A_z{E}_K}}\right].\end{array}$$

Equations (34) and (35) provide us with two constrains on the function in the square brackets. First, the left-hand side of Eq. (34) is non-negative by definition, implying that the same requirement applies for the right-hand side of the equation. Second, the turbulent Prandtl number grows with the increase of the gradient Richardson number – that is, ${\left.{\mathit{Pr}}_T\right|}_{(z/L\to \mathrm{\infty })}\to \mathit{Ri}/R_{\mathrm{\infty }}$ – requiring the function in the square brackets to approach zero under extreme stratification. This leads us to the following approximation (see Fig. 4):

$$\begin{array}{}\text{(36)}& {\displaystyle \frac{\mathrm{1}-C_{\mathit{\theta }}}{\mathrm{1}+C_{\mathrm{\nabla }}}}{\displaystyle \frac{E_P}{A_z{E}_K}}=\mathrm{1}-e^{-C_{\mathit{Pr}}z/L}.\end{array}$$

This function monotonically decreases from 1 to $\mathrm{0}<z/L<\mathrm{\infty }$, satisfying our requirements with C_Pr=0.65. The observed spread of data points might be explained by the simulation time being insufficient to reach a fully statistical steady state for this specific ratio. Although the fully developed steady state was achieved (verified using the standard criterion of stabilised TKE, which showed no significant fluctuations over time), the parameter involving ratios of temperature fluctuations, θ, might require additional time to stabilise. We believe that increasing the experiment time would decrease the spread, but we leave the validation of this hypothesis for future studies.

This leads us to a similar approximation of $t_{\mathit{\tau }}/t_F$ (see Fig. 5):

$$\begin{array}{}\text{(37)}& \begin{array}{rl}{\displaystyle \frac{t_{\mathit{\tau }}}{t_F}}& ={\mathit{Pr}}_T\left(\mathrm{1}+C_{\mathrm{\nabla }}\right)\left[\mathrm{1}-{\displaystyle \frac{\mathrm{1}-C_{\mathit{\theta }}}{\mathrm{1}+C_{\mathrm{\nabla }}}}{\displaystyle \frac{E_P}{A_z{E}_K}}\right]\\ & =C_{\mathrm{1}}^{\mathit{\tau }F}{e}^{-C_{\mathrm{2}}^{\mathit{\tau }F}z/L}.\end{array}\end{array}$$

Figure 5 The ratio of the effective dissipation timescales of τ and F_z, $t_{\mathit{\tau }}/t_F$ versus $z/L$. The data are the same as in Fig. 2. The solid black line shows Eq. (37) with empirical constants $C_{\mathrm{1}}^{\mathit{\tau }F}=\mathrm{0.17}$ and $C_{\mathrm{2}}^{\mathit{\tau }F}=\mathrm{0.62}$ obtained from the best fit of Eq. (37) to DNS data in the turbulent layer, $z>\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}$.

Now, to complete the closure, we need to determine one more dimensionless ratio, $t_{\mathit{\theta }}/t_K$. It is explicitly required for the ratio of turbulent energies, $E_P/E_K$, and consequently for A_z through Eqs. (20) and (36). We approximate it once again with the ratio of two first-order polynomials:

$$\begin{array}{}\text{(38)}& {\displaystyle \frac{t_{\mathit{\theta }}}{t_K}}={\displaystyle \frac{C_{\mathrm{1}}^{\mathit{\theta }K}z/L+C_{\mathrm{2}}^{\mathit{\theta }K}}{z/L+C_{\mathrm{3}}^{\mathit{\theta }K}}}.\end{array}$$

Figure 6 The ratio of the dissipation timescale of 〈θ²〉 and the dissipation timescale of TKE, $t_{\mathit{\theta }}/t_K$ versus $z/L$. The data are the same as in Fig. 2. The solid black line shows Eq. (38) with empirical constants $C_{\mathrm{1}}^{\mathit{\theta }K}=\mathrm{40}$, $C_{\mathrm{1}}^{\mathit{\theta }K}=\mathrm{480}$, and $C_{\mathrm{1}}^{\mathit{\theta }K}=\mathrm{900}$ obtained from the best fit of Eq. (38) to DNS data in the turbulent layer, $z>\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}$.

With Eq. (38), our turbulence closure is now complete, allowing us to proceed with the verification process using quantities not utilised in the fitting procedures. Figure 7 provides empirical evidence that supports the stability dependencies given by Eqs. (20), (26), (27), and (34)–(38). Table 2 summarises the proposed approximations and provides a summary of the resulting turbulent closure.

Table 2 Proposed approximations and resulting revised turbulent parameters of EFB closure.

Figure 7 Validating the closure with quantities not utilised in the fitting procedures. Panel (a) shows the TPE-to-TKE ratio, $E_P/E_K$; panel (b) shows the vertical share of TKE, A_z; panel (c) demonstrates the ratio of dissipation timescales of 〈θ²〉 and F_z; panel (d) shows the turbulent Prandtl number, Pr_T; panel (e) shows the squared dimensionless turbulent flux of momentum, ${\left(\mathit{\tau }/E_K\right)}^{\mathrm{2}}$; and panel (f) shows the squared dimensionless turbulent flux of potential temperature, $F_z^{\mathrm{2}}/E_{\mathit{\theta }}{E}_K$. All quantities are plotted against $z/L$. The solid black lines correspond to theoretical predictions demonstrating acceptable to great agreement with the DNS data in the turbulent layer, $z>\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}$. Empirical data are from the same sources as in Fig. 2. No fitting has been performed for this figure.

For practical reasons, most operational numerical weather prediction and climate models parameterise these dimensionless ratios as functions of the gradient Richardson number rather than $z/L$. This preference arises from the fact that the gradient Richardson number is defined solely by mean quantities – namely, buoyancy and shear productions – which in practice imposes fewer computational restrictions on the model's time step. Since Ri=Pr_TRi_f and both Pr_T and Ri_f are defined as functions of $z/L$ by Eqs. (35) and (21), respectively, we can derive an expression for the gradient Richardson number, Ri, as the function of $z/L$, as shown in Fig. 8:

$$\begin{array}{}\text{(39)}& \mathit{Ri}={\mathit{Ri}}_{\mathrm{f}}{\displaystyle \frac{C_{\mathrm{1}}^{\mathit{\tau }F}}{\mathrm{1}+C_{\mathrm{\nabla }}}}{e}^{-\left(C_{\mathit{Pr}}-C_{\mathrm{2}}^{\mathit{\tau }F}\right)z/L}.\end{array}$$

Figure 8 The resulting approximation of the gradient Richardson number, Ri, after Eq. (39) compared to the exact solution (a) and the relative error of this approximation as a function of gradient Richardson number, Ri (b). The solid black line corresponds to theoretical derivation that shows good agreement with the DNS data in the turbulent layer, $z>\mathrm{50}\mathit{\nu }/{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}$. Empirical data are from the same sources as in Fig. 2. No fitting has been performed for this figure.

5 Concluding remarks

For many years, our understanding of dissipation rates for turbulent second-order moments has been hindered by a lack of direct observations in fully controlled conditions, particularly in a strongly stable stratification. To address this limitation, we conducted topical DNS experiments of stably stratified Couette flows. The main finding of this study is that, contrary to the traditional assumption of them being proportional to a single universal dissipation timescale, the ratios of the dissipation timescales of the basic second-order moments depend on the temperature stratification (e.g. characterised by the gradient Richardson number).

This finding laid the foundation for empirically approximating these ratios with simple universal functions of stability parameters, valid for a wide range of stratifications. Consequently, this allowed us to refine the EFB turbulent closure by accounting for dissipation timescales that are intrinsic to the basic second-order moments. As a result, the revised formulations for eddy viscosity and eddy conductivity reveal greater physical consistency in stratified conditions, thereby enhancing the representation of turbulence in numerical weather prediction and climate modelling.

We have also observed that the dimensionless parameters involving θ fluctuations demonstrate a wider spread of values within and across the DNS experiments, making it more challenging to approximate them with stability functions. This suggests that the stabilisation time for these parameters may be significantly longer than for TKE components.

It is important to note that our DNS experiments were limited to gradient Richardson numbers up to Ri=0.17. Any data reliably indicating different asymptotic values of the timescale dimensionless ratios or demonstrating their different dependency on the temperature stratification would impose the need to readjust the proposed parameterisation.

We deliberately avoided discussing intermittency issues; for that, one needs to determine higher-order two-point (or multi-point) moments. Intermittency is important for small-scale effects, and intermittency implies that higher-order moments of velocity and temperature fields have non-Gaussian statistics. In this study, we focused on larger scales, determining one-point second-order correlation functions that barely touch one-point third-order correlation functions only when it is necessary. However, addressing this topic would be crucial for advancing numerical simulations towards higher stratifications and warrants detailed investigation.

With these considerations in mind, we believe that the most challenging step will be to explicitly explore the transitional region between traditional weakly stratified turbulence and extremely stable stratification, where the behaviour of the turbulent Prandtl number shifts from nearly constant to a linear function with respect to the gradient Richardson number. Investigating this phenomenon would require unprecedented computational resources for DNS or specialised in situ or laboratory experiments.