Closure scheme for stably stratified turbulence without critical Richardson number

This paper presents a new modification of the second order turbulence closure that removes the critical gradient Richardson number limitation typically found in Mellor-Yamada style models. The mean wind speed and potential temperature profiles are derived for the newly modified model in terms of similarity and structure functions depending on the gradient Richardson number. The derivation is based on a second order boundary layer approximation in neutral to very stable stratification conditions. Some recent closure assumptions for pressure-temperature and heat flux are considered. Variances and covariances of the turbulent fluctuations are also investigated with respect to the gradient Richardson number. The new model predictions are confronted with some well known models. Second-order scheme in the framework of Mellor-Yamada type models employing new heat flux equations. The model does not exhibit a threshold for the gradient Richardson number. Mean wind and temperature profiles, and turbulent fluctuations equations as functions of the gradient Richardson number. Second-order scheme in the framework of Mellor-Yamada type models employing new heat flux equations. The model does not exhibit a threshold for the gradient Richardson number. Mean wind and temperature profiles, and turbulent fluctuations equations as functions of the gradient Richardson number.

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 ≈ 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 ≫ 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 surfacelayer 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 thirdorder 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.

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]) 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 Θ + . Quantities like uw , w 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 0 , respectively. The quantity g j = (0, 0, −g) is the gravity acceleration, is the thermal expansion coefficient, is the kinematic viscosity and the thermal diffusivity. The Coriolis parameter has been neglected as a standard simplification in the surface-layer.

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 , 1 , 2 , 3 and Λ 1 , Λ 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.

Pressure terms
The pressure terms are closed in the sense of "return-toisotropy" proposed by Rotta [30], namely an anisotropic flow tends to isotropy in absence of external forcings. We have, and 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.

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 and

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

Remaining terms
Since there is no isotropic first-order tensor, we have (see e.g. Mellor [24]) and pressure diffusional terms are small (see e.g. Mellor [24])

Boundary layer approximation
We consider the second-order turbulent scheme (1) -(3) with the closure assumptions (4) - (10). We assume where we oriented our coordinate system so that vw = 0 and uw is aligned with the mean wind vector U. Moreover, it could be shown (see Mellor [24]) that 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, and where A 1 , A 2 , B 1 and B 2 are positive coefficients and l represents a length scale such that l → z as z → 0 , where 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. (11)

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 pressuretemperature relaxation time scale [see Eq. (9c) in Canuto et al. [6]]. More precisely, they considered (13) with where t is the turbulent Prandtl number defined as with t0 its neutral value (the constant A ′ 2 is determined so that in the neutral conditions l 2 ∕l = A 2 ) and are the production terms due to shear and buoyancy, respectively. Here, Ri is the gradient Richardson number defined as where N is the Brunt-Väisälä frequency. With the modification given by (15), the equations for u and w change as follows (see Appendix A in Cheng et al. [8] for a more detailed derivation of the new heat flux equations) 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 SN Applied Sciences (2022) 4:214 | https://doi.org/10.1007/s42452-022-05088-8 Research Article

Physical considerations on l 2 ∕l
In order to give a physical insight into the relation (15), we consider the closure assumption (5) where, for simplicity, C 3 = C 4 = 0 . Relation (21) could be rewritten as follows where p is a time scale. Many second-order closure models assume p = ∕c p for stably stratified flows, with, for example, = l 2 ∕q and c p positive constant such that l 2 = c p l (in order to be consistent with the above discussion, c p should be equal to A 2 ). Canuto et al. [6] suggested the use of a buoyancy damping factor in the time scale p 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 p : with c 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 ≫ 1 , the above relation can be generalized as Now, introducing the flux Richardson number the Prandtl number from (16) is expressed as follows For Ri ≫ 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 Relation (24), together with (26), is consistent with (15).

Stability functions
Back to the framework of MY82, we express the turbulent stress uw and the heat flux w in terms of the well-known flux-gradient approximation 2 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]) Note that, from (27) and (28), we can rewrite the turbulent kinetic energy budget equation (second equation in (14)) as follows where = q 3 ∕Λ 1 is the turbulent kinetic energy dissipation. Equilibrium conditions require (P s + P b )∕ = 1 and we end up with Relation (29) will play a role in the subsequent analysis.

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) (30) and (31)

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 Θ as follows Here, u 2 * = −uw , H = −w , * = −H∕u * and L = u 3 * 0 ∕ gH is a length scale first introduced by Obukhov [29] that can be interpreted as a measure of the stability. Here, 0 and 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 ( w > 0 , H < 0 ), positive for stable cases ( w < 0 , H > 0 ) (see for example Garratt [12], Tampieri [34] and Wyngaard [36]). In the following, we are interested in the behavior of m and h . From the definition (33), we can compute where −w = u * * . Consequently, m and h can be expressed in terms of the gradient Richardson number as follows 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 with a limitation on the stability range 0 ≤ 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 parameterizations 3 Consequently, from (34) 1 and (36) we can derive m solving the following relation Given the profile of m , from (34) 2 we can compute h .

Variances and covariances of the turbulent fluctuations
Given the profile of m and 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, u 2 ∕u 2 * , v 2 ∕u 2 * and w 2 ∕u 2 * ; in the horizontal heat flux, u ∕u * * ; in the temperature variance, 2 ∕ 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:

Results
The similarity function 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, m with (35) exhibits a critical Richardson number ( Ri ≈ 0.2 ) like the commonly used log-linear similarity function from Högström [14] (violet). Instead, when (36) is considered, 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 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 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. 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 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 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 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 twotimes the turbulent kinetic energy as function of Ri. All the results exhibit similar behavior, showing a decreasing of the ratio w 2 ∕q 2 , indicative of a higher content of kinetic energy in the horizontal components. Figure 6 shows the ratio 2 ∕ 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.

Discussion
The similarity functions m and 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 z∕l . As already mentioned, the profiles of m and 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 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 m and h , the behaviors described in the relations (38) -(43) is not affected by the choice of the ratio 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 uw on the horizontal heat flux u . 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 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.

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