1 Introduction

Atmospheric turbulence is one of the key factors affecting the weather, climate, air quality as well as many engineering applications like wind turbines (farms) design, wind effects on buildings and structures, etc. Its understanding and modeling is thus a crucial issue in climate, weather and local environmental processes predictions. The low resolution models based on first-order closures served well for many years, at the beginning of computer based numerical prediction models. The increased computational power becoming available in past decades made possible to develop and use many high resolution models that also brought higher standards and requirements to parametrizations used within the underlying mathematical models. At that stage, the first-order turbulence closures became largely insufficient and obsolete, giving the rise to second and even higher order closures.

The problem of higher order turbulence closures is a subject of intensive investigation since at least the beginning of 1970’s (see e.g. Mellor [24]; Launder et al. [19]; Lumley [22]; André et al. [1]). This effort led soon to a hierarchy of turbulence closures proposed and summarized by Mellor and Yamada [25]. These parametrizations offered a deeper insight into the problem and became a basis for a number of usable strategies for numerical simulations of atmospheric flows. One of the most popular models was presented in Mellor and Yamada [26]; hereafter MY82. This approach was widely adopted and further developed in the following decades (Galperin et al. [10]; Canuto [5]; Shih and Shabbir [31]). Since the very beginning it was clear that this model, as any closure model, includes number of compromises and simplifications leading to various limitations. One of the most frequently discussed limitations of this kind of turbulence closure schemes is the prediction of stably stratified flows. The Mellor Yamada model implies the existence of a finite critical gradient Richardson number, beyond which the turbulence ceases and can not exist anymore (Kantha and Clayson [15, 16]; Nakanishi [28]; Cheng et al. [7]; Sukoriansky et al. [33]; Galperin et al. [11]; Canuto et al. [6]; Kantha and Carniel [17]; Bretherton and Park [2]; Zilitinkevich et al. [37]). The existence of a critical value of gradient Richardson number (\(Ri_{cr}\approx 0.19\) for the MY82 model) contradicts the observations of stably stratified flows where turbulence persists even at stratifications with \(Ri > Ri_{cr}\), i.e. turbulence survives at \(Ri \gg 1\).

The problem of critical Richardson number was discussed and addressed by a number of authors (among others Cheng et al. [7], Canuto et al. [6], Kantha and Carniel [17], Zilitinkevich et al. [37], Li et al. [20]). All these efforts were focused either on the increase of the \(Ri_{cr}\) to an acceptable level (e.g. model of Cheng et al. [7] increases \({Ri}_{cr}\) to O(1), which is a fair estimation in geophysical flows) or even into complete removal of this barrier (e.g. models by Canuto et al. [6], Zilitinkevich et al. [37]).

The work presented in this paper aims to propose a kind of modified second-order closure (in the sense of MY82 model), that does not have any critical Richardson number limitations. The modification is based on variation of the length scale associated with the return-to-isotropy, following the ideas from Cheng et al. [8], Canuto et al. [6], incorporating an asymptotic limit to the dynamic length scale according to Nakanishi [28]. The dimensionless mean wind and temperature gradient profiles in the surface-layer as well as the normalized variances and covariances of the turbulent fluctuations are derived using the similarity functions in the framework of the Monin-Obukhov Similarity Theory [27]. The results are put into context and direct comparison with some previous relevant works in this area. The present analysis extends the previous work by Caggio et al. [4] (see also Caggio and Bodnár [3]) mainly in terms of the investigation of second-order quantitites.

The paper is organized as follows. In Sect. 2 we introduce the second-order scheme for the turbulent variables, together with the parameterizations of the third-order terms, in the boundary layer approximations. Next, we discuss the new heat flux equations. Then, we recover the framework of the MY82 model and we compute the appropriate quantities in order to derive the mean wind speed and potential temperature profiles in terms of similarity functions and the variances and covariances of the turbulent fluctuations. Section 3 and 4 are devoted to the presentation and discussion of our results, and to the conclusions.

2 Second-order scheme

In the following, we introduce the general second-order scheme for the turbulent variables (see e.g. Garrat [12], Tampieri [34], Wyngaard [36])

$$\frac{{\partial \overline{{u_{i} u_{j} }} }}{{\partial t}} + \frac{\partial }{{\partial x_{k} }}\left[ {U_{k} \overline{{u_{i} u_{j} }} + \overline{{u_{k} u_{i} u_{j} }} - \nu \frac{\partial }{{\partial x_{k} }}\overline{{u_{i} u_{j} }} } \right] + \frac{1}{{\rho _{0} }}\left( {\frac{\partial }{{\partial x_{j} }}\overline{{pu_{i} }} + \frac{\partial }{{\partial x_{i} }}\overline{{pu_{j} }} } \right) = - \overline{{u_{k} u_{i} }} \frac{{\partial U_{j} }}{{\partial x_{k} }} - \overline{{u_{k} u_{j} }} \frac{{\partial U_{i} }}{{\partial x_{k} }} - \frac{1}{{\theta _{0} }}\left( {g_{j} \overline{{u_{i} \theta }} + g_{i} \overline{{u_{j} \theta }} } \right) + \frac{1}{{\rho _{0} }}\overline{{p\left( {\frac{{\partial u_{i} }}{{\partial x_{j} }} + \frac{{\partial u_{j} }}{{\partial x_{i} }}} \right)}} {\text{ }} - 2\nu \overline{{\frac{{\partial u_{i} }}{{\partial x_{k} }}\frac{{\partial u_{j} }}{{\partial x_{k} }}}} ,$$
(1)
$$\frac{{\partial \overline{{u_{j} \theta }} }}{{\partial t}} + \frac{\partial }{{\partial x_{k} }}\left[ {U_{k} \overline{{\theta u_{j} }} + \overline{{u_{k} u_{j} \theta }} } \right. - \left. {\alpha \overline{{u_{j} \frac{{\partial \theta }}{{\partial x_{k} }}}} - \nu \overline{{\theta \frac{{\partial u_{j} }}{{\partial x_{k} }}}} } \right] + \frac{1}{{\rho _{0} }}\left( {\frac{\partial }{{\partial x_{j} }}\overline{{p\theta }} } \right) $$$$\begin{aligned} =-\overline{u_{j}u_{k}}\frac{\partial \Theta }{\partial x_{k}}-\overline{\theta u_{k}}\frac{\partial U_{j}}{\partial x_{k}}-\frac{1}{\theta _{0}}g_{j}\overline{\theta ^{2}}+\frac{1}{\rho _{0}}\left( \overline{p\frac{\partial \theta }{\partial x_{j}}}\right) \nonumber \\ -\left( \alpha +\nu \right) \overline{\frac{\partial u_{j}}{\partial x_{k}}\frac{\partial \theta }{\partial x_{k}}}, \end{aligned}$$
(2)
$$\begin{aligned} \frac{\partial \overline{\theta ^{2}}}{\partial t}+\frac{\partial }{\partial x_{k}}\left[ U_{k}\overline{\theta ^{2}}+\overline{u_{k}\theta ^{2}}-\alpha \frac{\partial \overline{\theta ^{2}}}{\partial x_{k}}\right] =-2\overline{u_{k}\theta }\frac{\partial \Theta }{\partial x_{k}}\nonumber \\ -2\alpha \overline{\frac{\partial \theta }{\partial x_{k}}\frac{\partial \theta }{\partial x_{k}}}. \end{aligned}$$
(3)

Here and hereafter, capital and tiny letters denote the mean physical quantity and its turbulent fluctuation respectively in the sense of Reynolds decomposition. The ensemble averages are marked by overbar. In this sense we split the wind speed in \(U + u\) and the potential temperature in \(\Theta + \theta\). Quantities like \(\overline{uw}\), \(\overline{w\theta }\) or analogous, are interpreted as turbulent stress and heat fluxes, respectively. Pressure fluctuations and base state profile for the density are denoted by p and \(\rho _0\), respectively. The quantity \(g_{j} = (0, 0, -g)\) is the gravity acceleration, \(\beta\) is the thermal expansion coefficient, \(\nu\) is the kinematic viscosity and \(\alpha\) the thermal diffusivity. The Coriolis parameter has been neglected as a standard simplification in the surface-layer.

2.1 Closure assumptions

In order to close the system (1) – (3), we need to express the third-order terms as functions of second-order quantities. In the following, \(l_{1},l_{2}\), \(\lambda _{1},\lambda _{2}, \lambda _{3}\) and \(\Lambda _{1},\Lambda _{2}\) will be length scales, \(C_{1},C_{2},C_{3}\) and \(C_{4}\) positive coefficients and the quantity q will represents the square-root of two-times the turbulent kinetic energy.

2.1.1 Pressure terms

The pressure terms are closed in the sense of “return-to-isotropy” proposed by Rotta [30], namely an anisotropic flow tends to isotropy in absence of external forcings. We have,

$$\begin{aligned} \overline{{p\left( {\frac{{\partial u_{i} }}{{\partial x_{j} }} + \frac{{\partial u_{j} }}{{\partial x_{i} }}} \right)}} = & - \frac{q}{{3l_{1} }}\left( {\overline{{u_{i} u_{j} }} - \frac{{\delta _{{ij}} }}{3}q^{2} } \right) + C_{1} q^{2} \left( {\frac{{\partial U_{i} }}{{\partial x_{j} }} + \frac{{\partial U_{j} }}{{\partial x_{i} }}} \right) \\ + C_{2} \beta \left( {g_{i} \overline{{u_{j} \theta }} + g_{j} \overline{{u_{i} \theta }} - \frac{2}{3}\delta _{{ij}} g_{k} \overline{{u_{k} \theta }} } \right) \\ \end{aligned}$$
(4)

and

$$\begin{aligned} \overline{p\frac{\partial \theta }{\partial x_{j}}}=-\frac{q}{3l_{2}}\overline{u_{j}\theta }+C_{3}\beta g_{j}\overline{\theta ^{2}}+C_{4}\overline{u_{k}\theta }\frac{\partial U_{i}}{\partial x_{k}}. \end{aligned}$$
(5)

More precisely, Rotta [30] suggested a closure for the terms (4) and (5) with \(C_2=C_3=C_4=0\). The extension that concerns the presence of the heat fluxes and the temperature variance is possible to find in Yamada [35] and Nakanishi [28] (see also Denby [9]) with different values for \(C_2\), \(C_3\) and \(C_4\). In the following analysis we will keep this extension.

2.1.2 Third-order covariances

The third-order covariances of velocity components and temperature are expressed in terms of the flux-gradient approximation (see also relation (27) below), namely

$$\begin{aligned} \overline{{u_{k} u_{i} u_{j} }} = & - q\lambda _{1} \left( {\frac{{\partial \overline{{u_{i} u_{j} }} }}{{\partial x_{k} }} + \frac{{\partial \overline{{u_{i} u_{k} }} }}{{\partial x_{j} }} + \frac{{\partial \overline{{u_{j} u_{k} }} }}{{\partial x_{i} }}} \right),\;\;\overline{{u_{k} u_{j} \theta }} \\ & = - q\lambda _{2} \left( {\frac{{\partial \overline{{u_{k} \theta }} }}{{\partial x_{j} }} + \frac{{\partial \overline{{u_{j} \theta }} }}{{\partial x_{k} }}} \right) \\ \end{aligned}$$
(6)

and

$$\begin{aligned} \overline{u_{k}\theta ^{2}}=-q\lambda _{3}\frac{\partial \overline{\theta ^{2}}}{\partial x_{k}}. \end{aligned}$$
(7)

2.1.3 Dissipative terms

The dissipative terms are expressed according to the Kolmogorov hypothesis of small-scale isotropy (see Kolmogorov [18]) as

$$\begin{aligned} 2\nu \overline{\frac{\partial u_{i}}{\partial x_{k}}\frac{\partial u_{j}}{\partial x_{k}}}=\frac{2}{3}\frac{q^{3}}{\Lambda _{1}}\delta _{ij}, \ \ \ \ 2\alpha \overline{\frac{\partial \theta }{\partial x_{k}}\frac{\partial \theta }{\partial x_{k}}}=2\frac{q}{\Lambda _{2}}\overline{\theta ^{2}}. \end{aligned}$$
(8)

2.1.4 Remaining terms

Since there is no isotropic first-order tensor, we have (see e.g. Mellor [24])

$$\begin{aligned} \left( \alpha +\nu \right) \overline{\frac{\partial u_{j}}{\partial x_{k}}\frac{\partial \theta }{\partial x_{k}}}=0 \end{aligned}$$
(9)

and pressure diffusional terms are small (see e.g. Mellor [24])

$$\begin{aligned} \overline{pu_{i}}=\overline{p\theta }=0. \end{aligned}$$
(10)

2.2 Boundary layer approximation

We consider the second-order turbulent scheme (1) – (3) with the closure assumptions (4) – (10). We assume horizontally homogeneous conditions in the steady state and, for high Reynolds’ number, we neglect diffusion terms. The system of equations reads as followsFootnote 1

$$\begin{aligned} \overline{u^{2}}=\gamma _{1}q^{2}+2C_{2}\frac{l_{1}}{q}\beta g\overline{w\theta }-6\frac{l_{1}}{q}\overline{uw}\frac{\partial U}{\partial z}, \ \ \ \ \overline{v^{2}}=\gamma _{1}q^{2}\nonumber \\ +2C_{2}\frac{l_{1}}{q}\beta g\overline{w\theta }, \ \ \ \ \overline{w^{2}}=\gamma _{1}q^{2}+2(3-2C_{2})\frac{l_{1}}{q}\beta g\overline{w\theta }, \end{aligned}$$
(11)
$$\begin{aligned} \overline{uw}=3\frac{l_{1}}{q}\left[ -\left( \overline{w^{2}}-C_{1}q^{2}\right) \frac{\partial U}{\partial z}+(1-C_{2})\beta g\overline{u\theta }\right] , \end{aligned}$$
(12)
$$\begin{aligned} \overline{u\theta }=3\frac{l_{2}}{q}\left[ -\overline{uw}\frac{\partial \Theta }{\partial z}-(1-C_{4})\overline{w\theta }\frac{\partial U}{\partial z}\right] , \ \ \ \ \overline{w\theta }\nonumber \\ =3\frac{l_{2}}{q}\left[ -\overline{w^{2}}\frac{\partial \Theta }{\partial z}+(1-C_{3})\beta g\overline{\theta ^{2}}\right] , \end{aligned}$$
(13)
$$\begin{aligned} \overline{\theta ^{2}}=-\frac{\Lambda _{2}}{q}\overline{w\theta }\frac{\partial \Theta }{\partial z}, \ \ \ \ \frac{q^{3}}{\Lambda _{1}}=\left( -\overline{uw}\frac{\partial U}{\partial z}+\beta g\overline{w\theta }\right) , \end{aligned}$$
(14)

where we oriented our coordinate system so that \(\overline{vw}=0\) and \(\overline{uw}\) is aligned with the mean wind vector U. Moreover, it could be shown (see Mellor [24]) that \(\partial V/\partial z=\overline{uv}=\overline{v\theta }=0\).

Note that the second equation in (14) is the budget for the turbulent kinetic energy under equilibrium conditions. Accordingly, the production of turbulent kinetic energy by shear and buoyancy is totally balanced by the dissipation of turbulence. This is consistent with the approximation previously discussed.

In the system above,

$$\begin{aligned} \gamma _{1}=\frac{1}{3}-2\frac{A_{1}}{B_{1}} \end{aligned}$$

and

$$\begin{aligned} (l_{1}, l_{2}, \Lambda _{1}, \Lambda _{2})=(A_{1}, A_{2}, B_{1}, B_{2})l, \end{aligned}$$

where \(A_{1}\), \(A_{2}\), \(B_{1}\) and \(B_{2}\) are positive coefficients and l represents a length scale such that \(l \rightarrow \kappa z\) as \(z \rightarrow 0\), where \(\kappa\) is the von Kármán constant and can be prescribed or solved from a prognostic equation. Physical considerations about l will be, in the following, the key point of our analysis.

2.3 Length scale and new heat flux equations

Under the boundary layer approximation discussed in the previous section, the system (11) – (14) is presented in the framework of the MY82 scheme. In a very recent paper, Cheng et al. [8] focused in a modification of the heat flux equations proposing a new closure for the return-to-isotropy term of the pressure-temperature correlation (5). This closure was previously proposed by Canuto et al. [6] when they examined the wavenumber spectra of the pressure-temperature relaxation time scale [see Eq. (9c) in Canuto et al. [6]]. More precisely, they considered (13)

$$\begin{aligned} \overline{u\theta }=3\frac{l_{2}}{q}\left[ -\overline{uw}\frac{\partial \Theta }{\partial z}-(1-C_{4})\overline{w\theta }\frac{\partial U}{\partial z}\right] , \ \ \ \ \overline{w\theta }\\ =3\frac{l_{2}}{q}\left[ -\overline{w^{2}}\frac{\partial \Theta }{\partial z}+(1-C_{3})\beta g\overline{\theta ^{2}}\right] \end{aligned}$$

with

$$\begin{aligned} \frac{l_2}{l} = \frac{A_2'}{(1 + \sigma _t)}, \ \ A_2' = A_2 (1 + \sigma _{t0}), \end{aligned}$$
(15)

where \(\sigma _t\) is the turbulent Prandtl number defined as

$$\begin{aligned} \sigma _t = \frac{-P_s {Ri}}{P_b} \end{aligned}$$
(16)

with \(\sigma _{t0}\) its neutral value (the constant \(A_2'\) is determined so that in the neutral conditions \(l_2/l = A_2\)) and

$$\begin{aligned} P_s = -\overline{uw}\frac{\partial U}{\partial z}, \ \ P_b = \beta g\overline{w\theta } \end{aligned}$$
(17)

are the production terms due to shear and buoyancy, respectively. Here, Ri is the gradient Richardson number defined as

$$\begin{aligned} {Ri} = \frac{N^2}{S^2}, \ \ N^2 = \beta g \frac{\partial \Theta }{\partial z}, \ \ S^2 = \left( \frac{\partial U}{\partial z}\right) ^2, \end{aligned}$$
(18)

where N is the Brunt-Väisälä frequency. With the modification given by (15), the equations for \(\overline{u \theta }\) and \(\overline{w \theta }\) change as follows (see Appendix A in Cheng et al. [8] for a more detailed derivation of the new heat flux equations)

$$\begin{aligned} \overline{u\theta }=-3\frac{A_2'l}{q}\left[ (1-C_{4})\overline{w\theta }\frac{\partial U}{\partial z}\right] , \ \ \overline{w\theta }\nonumber \\ =-3\frac{A_2'l}{q}\left[ \overline{w^{2}}\frac{\partial \Theta }{\partial z}-(1-C_{3})\beta g\overline{\theta ^{2}}\right] - \left( \overline{uw}\frac{\partial \Theta }{\partial z}\right) /\frac{\partial U}{\partial z}. \end{aligned}$$
(19)

Consequently, in the following analysis we will consider the system (11), (12), (14), (19) with values of the coefficients from Cheng et al. [8], namely

$$\begin{aligned} A_1 = 0.92, \ \ A_2' = 1.332, \ \ B_1 = 16.6, \ \ B_2 = 10.1, \ \ C_1 = 0.08, \ \ C_2 = 0.25, \ \ C_3 = 0.22, \ \ C_4 = 0. \end{aligned}$$
(20)

2.3.1 Physical considerations on \(l_2/l\)

In order to give a physical insight into the relation (15), we consider the closure assumption (5)

$$\begin{aligned} \overline{p\frac{\partial \theta }{\partial x_{j}}}=-\frac{q}{3l_{2}}\overline{u_{j}\theta }, \end{aligned}$$
(21)

where, for simplicity, \(C_3=C_4=0\). Relation (21) could be rewritten as follows

$$\begin{aligned} \overline{p\frac{\partial \theta }{\partial x_{j}}}= - \frac{\overline{u_{j}\theta }}{\tau _{p\theta }} , \end{aligned}$$
(22)

where \(\tau _{p\theta }\) is a time scale. Many second-order closure models assume \(\tau _{p\theta }=\tau /c_{p\theta }\) for stably stratified flows, with, for example, \(\tau = l_2 / q\) and \(c_{p\theta }\) positive constant such that \(l_2 = c_{p\theta } l\) (in order to be consistent with the above discussion, \(c_{p\theta }\) should be equal to \(A_2\)). Canuto et al. [6] suggested the use of a buoyancy damping factor in the time scale \(\tau _{p\theta }\) in order to reduce the effect of eddies that, working against the gravity, lose kinetic energy that is converted to potential energy. This translates in the following formulation for \(\tau _{p\theta }\):

$$\begin{aligned} \tau _{p\theta } = \left( \frac{1}{c_{p\theta }} \right) \frac{\tau }{1+c_{\tau }\tau ^2N^2}, \end{aligned}$$
(23)

with \(c_{\tau }\) positive constant. Moreover, Canuto et al. [6] (see Section 5, relation (14e)) also showed that, including shear and buoyancy, for high gradient Richardson number, \({Ri}\gg 1\), the above relation can be generalized as

$$\begin{aligned} \tau _{p\theta } = \left( \frac{1}{c_{p\theta }} \right) \frac{\tau }{1+{Ri}}. \end{aligned}$$
(24)

Now, introducing the flux Richardson number

$$\begin{aligned} R_{f}=\frac{g}{\theta _{0}}\overline{w\theta }/\overline{uw}\frac{\partial U}{\partial z}, \end{aligned}$$
(25)

the Prandtl number from (16) is expressed as follows

$$\begin{aligned} \sigma _t = \frac{{Ri}}{R_f}. \end{aligned}$$

For \({Ri}\gg 1\), \(R_f\) tends to a constant value that we denote \(R_{fc}\) (see, for example, [8], Fig. 2b) and, consequently, the Prandtl number increases linearly with Ri (see, for example, [8], Fig.2a), namely

$$\begin{aligned} \sigma _t \approx \frac{{Ri}}{R_{fc}}. \end{aligned}$$
(26)

Relation (24), together with (26), is consistent with (15).

2.4 Stability functions

Back to the framework of MY82, we express the turbulent stress \(\overline{uw}\) and the heat flux \(\overline{w\theta }\) in terms of the well-known flux-gradient approximationFootnote 2

$$\begin{aligned} \overline{uw}=-K_{m}\frac{\partial U}{\partial z}=-lqS_{m}\frac{\partial U}{\partial z}, \ \ \ \ \overline{w\theta }=-K_{h}\frac{\partial \Theta }{\partial z}=-lqS_{h}\frac{\partial \Theta }{\partial z}. \end{aligned}$$
(27)

Here, \(K_{m}\) and \(K_{h}\) are the eddy diffusivities, functions of l, q and non-dimensional stability functions \(S_{m}\) and \(S_{h}\). Moreover, we define the non-dimensional gradients (see Mellor and Yamada [26])

$$\begin{aligned} G_{m}=\frac{l^{2}}{q^{2}}\left( \frac{\partial U}{\partial z}\right) ^{2},\ \ \ \ G_{h}=-\frac{l^{2}}{q^{2}}\frac{g}{\theta _{0}}\frac{\partial \Theta }{\partial z}. \end{aligned}$$
(28)

Note that, from (27) and (28), we can rewrite the turbulent kinetic energy budget equation (second equation in (14)) as follows

$$\begin{aligned} \frac{P_{s}+P_{b}}{\varepsilon }=B_{1}\left( S_{m}G_{m}+S_{h}G_{h}\right) , \end{aligned}$$

where \(\varepsilon =q^{3} / \Lambda _{1}\) is the turbulent kinetic energy dissipation. Equilibrium conditions require \((P_{s}+P_{b})/\varepsilon =1\) and we end up with

$$\begin{aligned} B_{1}\left( S_{m}G_{m}+S_{h}G_{h}\right) =1. \end{aligned}$$
(29)

Relation (29) will play a role in the subsequent analysis.

2.4.1 Derivation of \(S_m\) and \(S_h\)

In order to derive the dimensionless mean wind and temperature gradient profiles (see (33) in Section 2.5 below), we need to express the stability functions \(S_{m}\) and \(S_{h}\), as well as the non-dimensional gradient \(G_m\) and \(G_h\) in terms of the gradient Richardson number. A first step is to derive the stability functions \(S_{m}\) and \(S_{h}\) in terms of the non-dimensional gradient \(G_h\). From (11), (12) , (14), (19), (27) and (28), we have (for a detailed derivation see Appendix A below)

$$\begin{aligned} S_h = C / \left[ 1-G_h (A+B)\right] \end{aligned}$$
(30)

and

$$\begin{aligned} S_m = \left[ c' - (c'(A+B) -C(a' + b')) G_h \right] / \left[ 1-G_h (A+B)\right] \end{aligned}$$
(31)

where

$$\begin{aligned} A = a - a', \ \ B = b - b', \ \ C = c - c'\\ a' = 9 A_1 A_2' (1-C_2)(1-C_4), \ \ b' \\ = 6 A_1^2 (3 - 2C_2), \ \ c' = 3A_1 (\gamma _1 - C_1)\\ a = 3 A_2' B_2 (1-C_3), \ \ b \\ = 6 A_1 A_2' (3-2C_2), \ \ c = 3A_2' \gamma _1. \end{aligned}$$

Now, an expression for the non-dimensional gradient \(G_h\) as a function of the gradient Richardson number (see (18)\(_1\)) is obtained from the turbulent kinetic energy balance (29). Indeed, from (29), we have

$$\begin{aligned} S_m + S_h \frac{G_h}{G_m} - \frac{1}{B_1 G_m} = 0 \end{aligned}$$

and using the definition of the gradient Richardson number in (18), we have

$$\begin{aligned} S_m - S_h{Ri} - \frac{1}{B_1 G_m} = 0. \end{aligned}$$

Now, setting

$$\begin{aligned} s_0 = C, \ \ s_2 = c', \ \ s_3 = c' (A+B)-C(a' + b'), \ \ d_1 = A+B, \end{aligned}$$

from the definition (30) and (31), we can derive

$$\begin{aligned} \frac{s_2-s_3 G_h}{1-d_1 G_h} - {Ri} \frac{s_0}{1-d_1 G_h} - \frac{1}{B_1 G_m} = 0 \end{aligned}$$

and, after some algebra, we end up with

$$\begin{aligned} c_2 G_h^2 + c_1 G_h +c_0 = 0, \end{aligned}$$
(32)

with \(c_0 = - {Ri}\), \(c_1 = (B_1 s_0 + d_1){Ri} - B_1 s_2\) and \(c_2 = B_1 s_3.\) Given \(G_h\) from (32), we derive \(G_m\) in terms of the gradient Richardson number from (29).

2.5 Similarity functions

In the framework of the Monin-Obukhov Similarity Theory we define the non-dimensional vertical gradients of the mean wind speed U and mean potential temperature \(\Theta\) as follows

$$\begin{aligned} \frac{\kappa z}{u_{*}}\frac{\partial U}{\partial z}\equiv \phi _{m}\left( \frac{z}{L}\right) , \ \ \ \ \frac{\kappa zu_{*}}{H}\frac{\partial \Theta }{\partial z}=-\frac{\kappa z}{\theta _{*}}\frac{\partial \Theta }{\partial z}\equiv \phi _{h}\left( \frac{z}{L}\right) . \end{aligned}$$
(33)

Here, \(u_{*}^{2}=- \overline{uw}\), \(H=-\overline{w\theta }\), \(\theta _{*}=-H/u_{*}\) and \(L=u_{*}^{3}\theta _{0}/\kappa gH\) is a length scale first introduced by Obukhov [29] that can be interpreted as a measure of the stability. Here, \(\theta _0\) and \(\kappa\) are the base state profile for the potential temperature and the von Kármán constant, respectively. The sign of the heat flux determines the sign of L, negative for unstable cases (\(\overline{w\theta }>0\), \(H<0\)), positive for stable cases (\(\overline{w\theta }<0\), \(H>0\)) (see for example Garratt [12], Tampieri [34] and Wyngaard [36]). In the following, we are interested in the behavior of \(\phi _m\) and \(\phi _h\). From the definition (33), we can compute

$$\begin{aligned} \frac{\kappa z}{u_*} \frac{\partial U}{\partial z} = \frac{\kappa z}{l} \left( \frac{l}{q} \frac{\partial U}{\partial z} \right) ^{1/2} \left( \frac{lq \frac{\partial U}{\partial z}}{-\overline{uw}} \right) ^{1/2} = \frac{G_m^{1/4}}{S_m^{1/2}} \frac{\kappa z}{l}, \end{aligned}$$

where \(- \overline{uw} = u_{*}^{2}\). Similarly,

$$\begin{aligned} \frac{\kappa z}{\theta _{*}}\frac{\partial \Theta }{\partial z} =\frac{\kappa z}{u_*} \frac{\partial U}{\partial z} \left( \frac{- \overline{uw}}{lq \frac{\partial U}{\partial z}} \right) \left( \frac{lq \frac{\partial \Theta }{\partial z}}{- \overline{w \theta }} \right) =\phi _m \frac{S_m}{S_h}, \end{aligned}$$

where \(- \overline{w\theta } = u_{*}\theta _{*}\). Consequently, \(\phi _m\) and \(\phi _h\) can be expressed in terms of the gradient Richardson number as follows

$$\begin{aligned} \phi _m = \frac{G_m^{1/4}}{S_m^{1/2}} \frac{\kappa z}{l}, \ \ \ \ \phi _h = \phi _m \frac{S_m}{S_h}. \ \ \ \ \end{aligned}$$
(34)

In Eq. (34) the length scale l must be parameterized. In near neutral conditions it can be assumed to be proportional to the distance from the surface z. A common extension to stable conditions is

$$\begin{aligned} l=\kappa z/(1+\widetilde{\alpha } z/L), \ \ \widetilde{\alpha }=2.7, \end{aligned}$$
(35)

with a limitation on the stability range \(0 \le z/L < 1\) (see Cheng et al. [8], relation (11e)). To exploit the behaviour of the present theory as the stability range is extended we refer again to Nakanishi [28], relation (39), who suggests a limit for l that shall not become smaller that the value at \(z/L=1\). Thus, the following formula is used in order to interpolate between the different parameterizationsFootnote 3

$$\begin{aligned} \frac{\kappa z}{l} = \frac{\widetilde{\beta } (1+\alpha '\frac{z}{L})}{\widetilde{\beta } + \alpha ' \frac{z}{L}} = \frac{\widetilde{\beta } (1+\alpha '\phi _m {Ri}\frac{S_h}{S_m})}{\widetilde{\beta } + \alpha ' \phi _m {Ri}\frac{S_h}{S_m}}, \ \ \left( \alpha ' = \frac{\widetilde{\alpha }}{1-\frac{1}{\widetilde{\beta }}}, \ \widetilde{\alpha }= 2.7, \ \widetilde{\beta }=3.7\right) . \end{aligned}$$
(36)

Consequently, from (34)\(_1\) and (36) we can derive \(\phi _m\) solving the following relation

$$\begin{aligned} b_2 \phi _m^2 + b_1 \phi _m + b_0 = 0, \end{aligned}$$
(37)

where

$$\begin{aligned} b_2 = \alpha ' {Ri} \frac{S_h}{S_m}, \ \ \ \ b_1 = \widetilde{\beta } \left( 1 - \alpha ' {Ri} \frac{S_h}{S_m} \frac{G_m^{1/4}}{S_m^{1/2}}\right) , \ \ \ \ b_0 = - \widetilde{\beta } \frac{G_m^{1/4}}{S_m^{1/2}}. \end{aligned}$$

Given the profile of \(\phi _m\), from (34)\(_2\) we can compute \(\phi _h\).

2.6 Variances and covariances of the turbulent fluctuations

Given the profile of \(\phi _m\) and \(\phi _h\), from (11), (12) , (14), (19) and (33), it is possible to derive the normalized variances and covariances of the turbulent fluctuations in terms of the gradient Richardson number. In particular, we are interested in the three-components of the velocity variances, \(\overline{u^{2}} / u_{*}^{2}\), \(\overline{v^{2}} / u_{*}^{2}\) and \(\overline{w^{2}} / u_{*}^{2}\); in the horizontal heat flux, \(\overline{u\theta } / u_{*}\theta _{*}\); in the temperature variance, \(\overline{\theta ^{2}} / \theta _{*}^2\); and in the kinetic energy \(q^{2} / u_{*}^{2}\). A derivation of these quantities is discussed in the Appendix B below. Here, we list the obtained results:

$$\begin{aligned} \frac{\overline{u^{2}}}{u_{*}^{2}}= \left( \frac{G_m^{1/4}}{S_m^{1/2}} \right) ^{2/3} \left[ B_1 \left( 1-Ri \frac{S_h}{S_m} \right) \right] ^{2/3} \nonumber \\ \left[ \gamma _1 + 2A_1 \frac{\left( 3-C_2 Ri \frac{S_h}{S_m}\right) }{B_1 \left( 1- Ri \frac{S_h}{S_m} \right) } \right] , \end{aligned}$$
(38)
$$\begin{aligned} \frac{\overline{v^{2}}}{u_{*}^{2}}= \left( \frac{G_m^{1/4}}{S_m^{1/2}} \right) ^{2/3} \left[ B_1 \left( 1-Ri \frac{S_h}{S_m} \right) \right] ^{2/3}\nonumber \\ \left[ \gamma _1 - 2A_1C_2 \frac{ Ri \frac{S_h}{S_m}}{B_1 \left( 1- Ri \frac{S_h}{S_m} \right) } \right] , \end{aligned}$$
(39)
$$\begin{aligned} \frac{\overline{w^{2}}}{u_{*}^{2}}= \left( \frac{G_m^{1/4}}{S_m^{1/2}} \right) ^{2/3} \left[ B_1 \left( 1-Ri \frac{S_h}{S_m} \right) \right] ^{2/3}\nonumber \\ \left[ \gamma _1 - 2A_1(3-2C_2) \frac{ Ri \frac{S_h}{S_m}}{B_1 \left( 1- Ri \frac{S_h}{S_m} \right) } \right] , \end{aligned}$$
(40)
$$\begin{aligned} \frac{\overline{u\theta }}{u_{*}\theta _{*}}=3A_{2}'(1-C_{4}) \left( \frac{G_m^{1/4}}{S_m^{1/2}} \right) ^{2/3} \left[ B_{1}\left( 1-Ri \frac{S_h}{S_m}\right) \right] ^{-1/3}, \end{aligned}$$
(41)
$$\begin{aligned} \frac{\overline{\theta ^{2}}}{\theta _{*}^{2}} =B_2 \left( \frac{G_m^{1/4}}{S_m^{1/2}} \right) ^{2/3} \left[ B_{1}\left( 1-Ri \frac{S_h}{S_m}\right) \right] ^{-1/3} \frac{S_m}{S_h}, \end{aligned}$$
(42)
$$\begin{aligned} \frac{q^{2}}{u_{*}^{2}}=\left[ B_{1}\frac{G_m^{1/4}}{S_m^{1/2}}\left( 1-Ri \frac{S_h}{S_m}\right) \right] ^{2/3}. \end{aligned}$$
(43)

3 Results

The similarity function \(\phi _m\) computed from Eq. (34)\(_1\) with parameterization (35) and (36) is reported in Fig. 1 (open circles and full squares, respectively) and compared with results from literature. It results that in the near neutral and weak stability range (\(Ri<0.1\)) all the similarity functions are equivalent, while they become quite different as stability increases. In particular, \(\phi _m\) with (35) exhibits a critical Richardson number (\(Ri \approx 0.2\)) like the commonly used log-linear similarity function from Högström [14] (violet). Instead, when (36) is considered, \(\phi _m\) exists for all \(Ri > 0\), so there is no threshold, and levels off for \(Ri > 0.5\). This is consistent with the similarity functions suggested by Sorbjan [32] (green) and by Gryanik et al. [13] (red), that are based on observations and thus are limited in the observed Ri range, and with the theoretical curve by Zilitinkevich et al. [37] (yellow), that is based on a different closure and also does not exhibit a critical Richardson number. The similarity function \(\phi _h\) computed from Eq. (34)\(_2\) with parameterization (35) and (36) is reported in Fig. 2 and compared with other functions from the literature. Similarly to \(\phi _m\), when (35) is used, a critical Ri is observed as in the log-linear relation from Högström [14], whilst no threshold for Ri occurs when (36) is considered, in agreement with Gryanik et al. [13] and, in particular, with the theoretical curve by Zilitinkevich et al. [37].

In Figs. 4 – 7 we report some of the similarity relations presented above concerning the variances and covariances of the turbulent fluctuations, namely the vertical velocity variance, the temperature variance and the turbulent kinetic energy from the relations (40), (42) and (43), respectively, and compared with some literature results; in particular, with the MY82 model, the Energy- and Flux-Budget model (EFB) developed by Zilitinkevich et al. [37] and the Cospectral Budget (CSB) model presented in Li et al. [20]. The EFB approach solves the budgets of turbulent momentum and heat fluxes and turbulent kinetic and potential energies, while the Cospectral Budget (CSB) approach is formulated in wavenumber space and integrated across all turbulent scales to obtain flow variables in physical space. Both approaches are recent advances in the study of stably stratified atmospheric flows and unlike the MY82 model allow turbulence to exist at any gradient Richardson number. However, as we will discuss below, a revision of the MY82 model allows also to avoid the threshold value for the gradient Richardson number.Footnote 4 The empirical curve suggested by Mauritsen and Svensson [23] is also presented in Fig. 6 and 7. This curve was obtained by identifying the weakly and very stable regimes and interpolating between them, from six different data sets, and can thus be interpreted as a summary of the field observations.

Figure 4 shows the ratio \(\overline{w^2}/u_*^2\). In the neutral and weak stability conditions, our result shows a good agreement with the MY82 model and the EFB model, and increasing slightly as stability increases. Similar increase is visible in the MY82 model with bigger values of the ratio \(\overline{w^2}/u_*^2\). For very stable conditions, our curve shows similar behavior as the MY82 model and the CSB model stabilizing on a constant value. Note that the MY82 behaves differently depending on the choice of the coefficient \(C_2\). While for \(C_2=0.3\) the model presents a divergent behavior, for \(C_2=1\) can encompass arbitrary large gradient Richardson numbers. The rational behind this choice will be better clarified in the next section. Note also that, while our result, together with the result of MY82, shows an increase of the ratio \(\overline{w^2}/u_*^2\) when stability increases, the EFB model shows a decrease of the same ratio. Figure 5 shows the ratio between the vertical velocity variance and two-times the turbulent kinetic energy as function of Ri. All the results exhibit similar behavior, showing a decreasing of the ratio \(\overline{w^2}/q^2\), indicative of a higher content of kinetic energy in the horizontal components. Figure 6 shows the ratio \(\overline{\theta ^{2}} / \theta _{*}^{2}\). In neutral conditions our result presents a good agreement with most of the results present in the literature. As the stability increases, our curve follows the the MY82 and EFB model. Concerning the turbulent kinetic energy (see Fig. 7), our result shows a good agreement with some of the results present in literature for the whole stability range. In particular, for very stable conditions our curve follows the MY82 and EFB model.

Fig. 1
figure 1

Blue squares: similarity function from Eq. (34) as function of Ri with the scale length parameterized as in Eq. (36). Open dots: the same, with the scale length from Eq. (35). Coloured continuous lines: expression from various authors

Full size image
Fig. 2
figure 2

Blue squares: similarity function from Eq. (34) as function of Ri with the scale length parameterized as in Eq. (36). Open dots: the same, with the scale length from Eq. (35). Coloured continuous lines: expression from various authors

Full size image
Fig. 3
figure 3

Ratio \(l/(\kappa z)\) as a function of Ri from Eqs. (35) and (36)

Full size image
Fig. 4
figure 4

Blue squares: vertical velocity variance from Eq. (40) as function of Ri. Coloured lines: expression from various authors

Full size image
Fig. 5
figure 5

Blue squares: ratio between the vertical velocity variance and two-times the turbulent kinetic energy from Eq. (40) and Eq. (43) as function of Ri. Coloured lines: expression from various authors

Full size image
Fig. 6
figure 6

Blue squares: temperature variance from Eq. (42) as function of Ri. Coloured lines: expression from various authors

Full size image
Fig. 7
figure 7

Blue squares: turbulent kinetic energy from Eq. (43) as function of Ri. Coloured lines: expression from various authors

Full size image

4 Discussion and conclusions

4.1 Discussion

The similarity functions \(\phi _m\) and \(\phi _h\) derived in the present analysis are the result of the combination between the second-order turbulent scheme discussed in Section 2.2 with the new heat flux equations proposed by Cheng et al. [8], together with the MY82 framework discussed in Section 2.4. and the formula (36) for the ratio \(\kappa z / l\). As already mentioned, the profiles of \(\phi _m\) and \(\phi _h\) do not present any threshold for the gradient Richardson number. In this, a crucial role is played by the parameterization proposed for the ratio \(\kappa z / l\). Indeed, we have noticed that the present approach, if the scale length is allowed to decrease without limit, as in Eq. (35), exhibits a singular behaviour at a finite value of Ri (see Fig. 3), and thus does not give a solution for larger Ri. This behaviour highligths the key role of the lenght scale parameterization.

Differently from \(\phi _m\) and \(\phi _h\), the behaviors described in the relations (38) –(43) is not affected by the choice of the ratio \(\kappa z / l\), but only by the parameterizations introduced in Section 2.1 and Section 2.3. Concerning the MY82 model, we would like to clarify the difference between the results due to the choice of the coefficient \(C_2\). As already mentioned, for \(C_2 = 1\) the model can encompasses arbitrary large gradient Richardson numbers, while it presents a divergent behavior for \(C_2 = 0.3\). This is due to the fact that the model predicts a threshold value for Ri and consequently a degeneration of turbulence (see Li et al. [20], Section 2 and Fig. 2). As discussed in Li et al. [20], Section 3, the rational behind the choice of \(C_2 = 1\) is that it eliminates the dependence of the turbulent momentum flux \(\overline{uw}\) on the horizontal heat flux \(\overline{u\theta }\). This revision allows the MY82 model to show results for any Ri. Similarly, the modification (15) proposed by Cheng et al. [8] modifies the equations of the heat fluxes in a way that a dependence on the momentum flux is not present in the equation of the horizontal heat flux, see (19)\(_{1}\).

Apart for the numerical values, the behavior described by the relation (42), (43) and the ratio between (40) and (43) is similar to the results of the literature presented here. Regarding the vertical velocity variance, (40), the difference in the behavior of our result, together with the result of MY82, with respect to the EFB model could be potentially related to the different parameterizations of the pressure terms. Indeed, while in our approach, as in the MY82 model, the coupling between the pressure fluctuations and the gradient of the velocity or temperature fluctuations is expressed in terms of “return-to-isotropy” in the sense of Rotta [30], the EFB model [37] parametrizes the pressure terms depending on the stability. This choice seems more suitable for stably stratified flows.

4.2 Conclusions

The analysis presented in this work concerned the modification of the second-order turbulence closure model that avoids the threshold for the gradient Richardson number through the introduction of new heat fluxes equations. In term of similarity and structure functions, the mean wind speed and potential temperature profiles have been derived for the new model together with the second-order quantities as a function of the gradient Richardson number. The predictions of the new model have been confronted with well know results from the literature.

Our findings are in a reasonable agreement with the results from the literature presented here. However, one important discrepancy is present in the profiles of the vertical velocity variance due to possibly, as already mentioned, the different parameterizations of the pressure terms. Consequently, it could be of some interest to investigate how the relations proposed in the EFB model could, potentially, modify the present results. In other words, it could be matter of future research to implement the pressure-terms parameterizations introduced by Zilitinkevich et al. [37] in the framework of the MY82 scheme with the new heat fluxes equations proposed by Cheng et al. [8]. Particular attention should also be given to the expression of the ratio \(\kappa z / l\) where different alternatives could be suggested.

We conclude that the present work could be seen as an extension, and perhaps an improvement, of the analysis developed by Cheng et al. [8], at least regarding the analysis of the similarity functions that do not present a critical value for the gradient Richardson number and the derivation of the equations describing the behavior of the variances and covariances of the turbulent fluctuations in terms of the gradient Richardson number.