The Stability Functions and Realizability of the Turbulent Kinetic Energy–Scalar Variance Closure for Moist Atmospheric Boundary-Layer Turbulence

The problem of realizability of the second-order turbulence closure models (parametrization schemes) is addressed through the consideration of the so-called “stability functions”. The emphasis is on the turbulence kinetic energy–scalar variance (TKESV) closure scheme that carries prognostic transport equations for the turbulence kinetic energy (TKE) and for the variances and covariance of two quasi-conservative scalars suitable for describing moist atmospheric boundary-layer turbulence. Stability functions appear within the framework of truncated closure schemes, where (i) the Reynolds-stress and scalar-flux equations (and, within the framework of one-equation TKE schemes, also equations for scalar variances and covariance) are reduced to the diagnostic algebraic formulations by neglecting the substantial derivatives and the third-order transport terms, and (ii) simplified linear parametrizations of the pressure-scrambling terms are used. The stability functions are ill-behaved (tend to infinity or become negative) over a certain range of governing parameters, e.g., mean velocity shear and buoyancy gradient. Using the approach of Helfand and Labraga (J Atmos Sci 45:113–132, 1988), we develop regularized stability functions for the TKESV scheme that reveal no pathological behaviour over their entire parameter space. The physical meaning of the regularization procedure and its relation to non-linear parametrizations of the pressure-scrambling terms and to weak non-equilibrium hypothesis are discussed. Finally, realizability of turbulence closures is considered within a more general framework of the moments problem of the probability theory.


Introduction
A turbulence parametrization scheme is an integral part of a physical parametrization package of any model of atmospheric circulation, including numerical weather prediction (NWP) and climate models. The importance of turbulence parametrization schemes will likely further increase as the resolution of the atmospheric models is refined. As the mesh size becomes small, quasi-organized flow structures, such as deep convective plumes that are chiefly responsible for non-local convective transport of momentum and scalars, are increasingly resolved. The focus of parametrizations of the subgrid-scale processes is then shifted towards motions at smaller scales and towards other issues, as, e.g., the anisotropy of turbulence near the surface and in stably-stratified regions of the flow, and the interaction between boundary-layer turbulence and shallow clouds. Turbulence closure schemes based on truncated second-order moment equations are viable tools for describing these features.
The so-called one-equation turbulence closure schemes have been very popular in geophysical applications over several decades. Those schemes carry the turbulence kinetic energy (TKE) equation, including the time rate-of-change (or substantial derivative if the TKE advection is included) and the third-order turbulent transport terms. All other second-moment equations, viz., the equations for the Reynolds stress, for the fluxes of scalar quantities (e.g., temperature and humidity), and for the scalar variances and covariance, are reduced to diagnostic algebraic expressions. The closure schemes, which carry prognostic transport equations not only for the TKE but also for the variances and covariance of scalar quantities, represent a natural step beyond one-equation TKE schemes. Examples are the schemes presented by Mellor and Yamada (1974, their level-3 closure), Cheng et al. (2002), Kenjereš andHanjalić (2002) andZeli et al. (2019) for the dry atmosphere, and by Nakanishi and Niino (2004) and Mironov and Machulskaya (2017) (see also Machulskaya and Mironov 2013) for the moist atmosphere. The latter scheme, termed the TKE-scalar variance (TKESV) scheme, receives major attention in the present study. An obvious advantage of the TKESV and similar schemes over the TKE schemes is a more consistent treatment of the flow energies. In atmospheric flows, where the buoyancy stratification is either unstable or stable but virtually never neutral, the TKE and the turbulence potential energy (TPE), characterized by the variances and covariances of scalar quantities, are equally important and should be treated in a similar way (see Mironov 2009, for discussion). An improved treatment of scalar (co)variances helps to improve the performance of atmospheric models in several respects. These are, for example, a consistent treatment of counter-gradient fluxes of scalars (which is not possible with the one-equation TKE schemes) and an improved representation of the fractional cloudiness within the framework of statistical cloud schemes.
In the present study, the problem of realizability of truncated second-order closure schemes for moist atmospheric boundary-layer turbulence is addressed through consideration of the so-called stability functions (Mellor and Yamada 1974). Stability functions appear within the framework of truncated closure schemes in the expressions for the Reynolds stress and scalar fluxes if (i) the Reynolds-stress and scalar-flux equations (and, within the framework of one-equation TKE schemes, also equations for the scalar variances and covariances) are reduced to the diagnostic algebraic formulations by neglecting the substantial derivatives and the third-order transport terms, and (ii) linear parametrizations (in the second-order moments involved) of the pressure-scrambling terms are used. The stability functions are ill-behaved over a certain range of governing parameters, e.g., mean velocity shear and buoyancy gradient. The ill-behaved stability functions yield physically meaningless values of the Reynolds stress and scalar fluxes, e.g., negative velocity variances or infinite scalar fluxes.
An appealing way to tackle the problem of ill behaviour of truncated closures was proposed by Helfand and Labraga (1988), see also du Vachat (1989) and Helfand and Labraga (1989). These authors analyzed the pathological behaviour of stability functions of the one-equation TKE scheme. Using rather plausible physical arguments, they developed regularized functions that reveal no pathological behaviour over their entire parameter space. In the present study, we extend the approach of Helfand and Labraga to the TKESV scheme. We develop regularized stability functions that are well-behaved over their entire parameter space and cause no pathological behaviour of the Reynolds stress and scalar fluxes. We discuss the physical meaning of the regularization procedure and examine its relation to the non-linear, more physically sound parametrizations of the pressure-scrambling terms and to weak nonequilibrium hypothesis often used to develop turbulence closures.
In what follows, a well-accepted terminology of Yamada (1974, 1982;Yamada 1977) is used. The level 3 refers to the TKESV and similar schemes that carry transport equations for the TKE and for the scalar (co)variances. The level 2 refers to the schemes, where all second-moment equations are reduced to the diagnostic algebraic formulations. The level 2.5 refers to one-equation TKE schemes.
An outline of the salient features of the TKESV scheme pertinent to the present analysis is given in Sect. 2. Section 3 presents the diagnostic algebraic equations for the Reynolds stress and scalar fluxes used within the framework of the TKESV scheme. Those equations are obtained with the aid of the boundary-layer approximation and are one-dimensional (the flow variables depend on the vertical coordinate only). In Sect. 4, the stability functions are introduced. The regularization procedure of Helfand and Labraga (1988) is presented in Sect. 5, and the regularized, well-behaved stability functions of the TKESV scheme are developed in Sect. 6. In Sect. 7, the relation of the regularization procedure to the nonlinear parametrizations of the pressure-scrambling terms is discussed. The relation of the Helfand and Labraga (1988) regularization procedure to the so-called weak non-equilibrium hypothesis (Rodi 1976;Gibson and Launder 1976;Hanjalić and Launder 2011) is discussed in Sect. 8. In Sect. 9, realizability of turbulence closures is considered within a more general framework of the moments problem of the probability theory. Conclusions are presented in Sect. 10. In order to keep the discussion in the main body of the text more concise, we include three Appendices. Turbulence potential energy is discussed in Appendix 1. The ill behaviour of stability functions in the shear-free flow is demonstrated analytically in Appendix 2. Appendix 3 outlines the derivation of regularized stability functions of the TKESV scheme in the case of shear-free flow.

Outline of the TKESV Scheme
In this section, a brief outline of the TKE-scalar variance turbulence parametrization scheme is given, where the focus is on the aspects pertinent to the present analysis. A detailed description of the scheme is given in Mironov and Machulskaya (2017). The governing equations are presented in their full three-dimensional form. The one-dimensional equations (where the flow variables depend on the vertical coordinate only) obtained with the aid of the boundary-layer approximation are used in the subsequent text.
The TKESV scheme is formulated in terms of variables that are approximately conserved for phase changes in the absence of precipitation (Betts 1973;Deardorff 1976). These are the total water specific humidity q t and the liquid water potential temperature θ l defined as where q is the water vapour specific humidity (the mass of water vapour per unit mass of moist air), q l is the liquid water specific humidity (the mass of cloud water per unit mass of moist air), c p is the specific heat of air at constant pressure, L v is the latent heat of vaporization, θ is the potential temperature related to the absolute temperature T through θ = T (P 0 /P) R d /c p , P andP 0 being the atmospheric pressure and its reference value, respectively, and R d being the gas constant for dry air. No supersaturation is assumed, so that q l = q t − q s if q t > q s , q s being the saturation specific humidity, and q l = 0 otherwise. Clearly, q t and θ l reduce to q and θ , respectively, in unsaturated conditions. In what follows, the total water specific humidity and the liquid water potential temperature are also referred to as simply humidity and temperature, respectively. The TKESV scheme uses prognostic transport equations for the TKE, e ≡ 1 2 u 2 i , for the scalar variances θ 2 l and q 2 t , and for the scalar covariance θ l q t . Here, u i is the velocity, and u i u j is the Reynolds stress. 1 The angle brackets denote mean quantities, and a prime denotes a turbulent fluctuation. The Einstein summation convention for repeated indices is adopted. The Boussinesq approximation is used, and the molecular diffusion terms in the second-moment budget equations are neglected (a good approximation for the atmospheric flows, where the Reynolds number is very high). The equations for the TKE and for the temperature-humidity covariance read Replacing q t with θ l in Eq. 3 yields the temperature-variance equation, and replacing θ l with q t yields the humidity-variance equation. In Eqs. 2 and 3, t is the time, x i are the space coordinates, β i = −g i α T is the buoyancy parameter, g i is the acceleration due to gravity, α T is the thermal expansion coefficient, ρ is the density, and p is the kinematic pressure (deviation of pressure from the hydrostatically balanced pressure divided by the reference density ρ r ). The last terms on the right-hand sides (r.h.s.) of Eqs. 2 and 3 are the molecular destruction (dissipation) rates of the TKE and of the scalar covariance. The quantity θ v is the virtual potential temperature that is defined with due regard for the water loading effect (Lilly 1968;Bannon 2007), where R v is the gas constant for water vapour. Equations 2 and 3 are not closed. Parametrizations (closure assumptions) are required for the buoyancy flux in the TKE equation (the term incorporating θ v ), for the dissipation rates, and for the third-order turbulent transport terms. The buoyancy terms are discussed at the end of this section. The dissipation rates of the TKE and of the scalar covariance (similar for the scalar variances) are parametrized in terms of the TKE dissipation time scale τ through the following algebraic expressions, where R τ is the dissipation time scale ratio. The TKE dissipation time scale is determined diagnostically through e and the turbulence length scale, see Eq. 11 below. For the thirdorder transport terms, either the simplest isotropic down-gradient parametrizations, or more advanced anisotropic parametrizations are used (see Mironov and Machulskaya 2017, for details). The specific form of those parametrizations are not required for the present analysis.
A key point is that the third-order transport terms along with the substantial derivatives (which reduce to the time rate-of-change in the one-dimensional case) are retained in the equations for the TKE and for the scalar variances and covariance. Within the framework of the TKESV scheme, the transport equations for the Reynolds stress u i u j and for the scalar fluxes u i θ l and u i q t are reduced (truncated) to the diagnostic algebraic expressions by neglecting the substantial derivatives and the third-order transport of u i u j , u i θ l , and u i q t , and the pressure transport of u i u j . In atmospheric flows, the Coriolis terms in the second-moment equations are typically small and can safely be neglected. The resulting algebraic equations are recast in terms of the departure-from-isotropy tensor The diagnostic algebraic equation for a i j (and hence for u i u j ) reads − 4 3 eS i j + a ik S jk + a jk S ik − 2 3 δ i j a km S km + a ik W jk + a jk W ik where are the symmetric and the antisymmetric parts, respectively, of the mean-velocity gradient tensor, and d i j = i j − 2 3 δ i j is the deviatoric part of the Reynolds-stress dissipation tensor. In order to obtain Eq. 7, a decomposition of the pressure gradient-velocity covariance u i ∂ p /∂ x j + u j ∂ p /∂ x i into pressure-strain (traceless) and pressure transport (diffusion) is used. The traceless part does not alter the TKE but acts to redistribute the energy between the components, thus driving turbulence towards the isotropic state (hence the term "pressure scrambling"). The non-uniqueness of decomposition of the pressure gradient-velocity covariance into pressure scrambling and pressure diffusion and the associated issues are briefly discussed in Mironov and Machulskaya (2017). The diagnostic algebraic equation for the temperature flux u i θ l reads where θi is the temperature-flux dissipation rate. Equation for the humidity flux u i q t is obtained from Eq. 8 by replacing θ l with q t . The following parametrizations of the pressure-scrambling terms in Eqs. 7 and 8 are used that are linear in the second-order moments, The expression for the pressure-scrambling term q t ∂ p /∂ x i in the humidity-flux equation is obtained from Eq. 10 by replacing θ l with q t . It is significant that the use of linear formulations of the pressure-scrambling terms yields a linear equation system for the Reynolds-stress and scalar-flux components. The following estimates of dimensionless constants are used: C u t = 1.8, C u s1 = 4/5, C u s2 = 3/5, C u s3 = 13/15, C u b = 3/10, C θ t = 5.0, C θ s1 = 3/5, C θ s2 = 1, and C θ b = 1/3. The constants C q t , C q s1 , C q s2 , and C q b in the humidity-flux equation are set equal to the respective constants in the temperature-flux equation. A detailed consideration of disposable constants and parameters of the TKESV scheme is presented in Mironov and Machulskaya (2017).
The deviatoric part of the Reynolds-stress dissipation rate d i j and the temperature-flux dissipation rate θi (similar for the humidity-flux dissipation rate qi ) are incorporated into the return-to-isotropy parts of the respective pressure-scrambling terms, the first terms on the r.h.s. of Eqs. 9 and 10, respectively.
The TKE dissipation time scale is expressed in terms of the TKE e and the turbulence length scale l as where C is a dimensionless constant. The turbulence length scale is computed through a Blackadar-type interpolation formula (Blackadar 1962) which includes a correction term due to stable density stratification (e.g., Stull 1973;Zeman and Tennekes 1977), or through a more advanced non-local formulation (e.g., Bougeault and Lacarrère 1989), see Mironov and Machulskaya (2017) for details. The buoyancy terms (the terms with β i ) in Eqs. 2 and 7-10 are parametrized with due regard for the presence of cloud condensate. The problem essentially amounts to modelling the virtual potential-temperature flux u i θ v and the scalarvirtual-potential-temperature covariances, θ l θ v and q t θ v , in terms of other moments that include fluctuations of u i , θ l , and q t (but not of q l that enters the definition of θ v ). The following formulation is used within the framework of the TKESV scheme, where a generic variable f stands for u i , θ l , or q t . In Eq. 12, , T l being the liquid water absolute temperature, is computed from the The variableR satisfies 0 ≤R ≤ 1. Equation 12 interpolates between the "dry" limit (R = 0), where a given grid box of a host atmospheric model is cloud-free (but the water vapour may be present), and the "wet" limit (R = 1), where the entire atmospheric-model grid box is saturated (i.e., q l = q t − q s > 0 at each point of the grid box). The interpolation variableR is related to the fractional cloud coverĈ which is determined with the aid of a statistical cloud scheme. In the stratus and stratocumulus regimes,R may be set equal toĈ. In the cumulus regime, highly localized cumulus clouds account for much of the buoyancy terms in the Reynolds-stress, scalar-flux and TKE equations, although the fractional cloud cover is typically small. Then,R should be (considerably) larger thanĈ.
A comprehensive discussion of various issues related to the parametrization of clouds and cloud-turbulence interaction in atmospheric models is beyond the scope of the present paper. Readers are referred to Sommeria and Deardorff (1977), Golaz et al. (2002), Tompkins (2003Tompkins ( , 2005, Machulskaya (2015), and Larson (2017). Essential for the present discussion is the fact that the covariances u i θ v , θ l θ v and q t θ v are expressed as linear combinations of u i θ l and u i q t , θ 2 l and θ l q t , and q 2 t and θ l q t , respectively, using Eq. 12.

Algebraic Reynolds-Stress and Scalar-Flux Equations
A simplification typically made in geophysical applications is the so-called boundary-layer approximation. This approximation is fairly accurate for large-scale and mesoscale atmospheric models whose grid-box aspect ratio (the ratio of the horizontal grid size to the vertical grid size) is large. Within the framework of the boundary-layer approximation, all derivatives in the x 1 and x 2 horizontal directions in the second-moment equations are neglected and the grid-box mean vertical velocity u 3 is zero in the second-moment equations (but not in the equations for the mean fields). The equations become one-dimensional, i.e., the flow variables depend on the x 3 vertical coordinate only. Substituting parametrizations of the pressure-scrambling terms (9) and (10) into Eqs. 7 and 8, respectively, and applying the boundary-layer approximation, we obtain the following algebraic equations for the departure-from-isotropy tensor (and hence for the Reynolds stress) and for the scalar fluxes, The estimates of dimensionless constants are given in Sect. 2. Within the framework of the boundary-layer approximation, , and other components of S i j and W i j are zero. The vertical x 3 axis is taken to be aligned with the vector of gravity, so that the only non-zero component of the buoyancy parameter is β 3 , which is denoted by β in what follows.
Equations 13-15 constitute a system of 12 linear equations for the components of a i j (six independent components), of u i θ l (three independent components), and of u i q t (three independent components). This linear system is readily solved, yielding explicit expressions for the Reynolds stress and scalar fluxes. Unfortunately, the solution to Eqs. 13-15 is nonrealizable. That is, it yields physically meaningless values of a i j (and hence of u i u j ), u i θ l , and u i q t , such as negative velocity variances or infinite scalar fluxes, over a certain range of governing parameters.
Note that the term "realizability" is used in a somewhat loose sense throughout much of the paper. Basically efforts are mounted to make sure that the quantities, which must be nonnegative and finite by definition (e.g., velocity and scalar variances), remain non-negative and finite. In Sect. 9, the notion of realizability is defined in a mathematically more rigorous way within the framework of the moments problem of the probability theory.

The Stability Functions
The system of linear equations 13-15 yields explicit expressions for the Reynolds-stress and scalar-flux components. The solutions for u 3 u 1 and u 3 u 2 (components of the vertical flux of horizontal momentum) have the following down-gradient form, The solutions for the vertical scalar fluxes contain both the down-gradient terms and the non-gradient terms that describe the generation of scalar fluxes by buoyancy, The quantities F M , F H1 , and F H2 are referred to as the stability functions (Mellor and Yamada 1974). It should be emphasized at once that the stability functions are merely the notation introduced to represent the formulations of fluxes in a compact form. Basically, F M , F H1 , and F H2 reflect the parametrization assumptions invoked to arrive at the algebraic Reynoldsstress and scalar-flux equations, first of all, the parametrizations of the pressure-scrambling terms.
The functions F M and F H1 that appear in the down-gradient terms in Eqs. 16-18 depend on three dimensionless governing parameters. These are the squared dimensionless shear (τ S) 2 , the squared dimensionless buoyancy frequency (τ N ) 2 , and the potential to kinetic energy ratio P/e. The stability function F H2 that appears in the non-gradient terms in Eqs. 17 and 18 depends on (τ S) 2 and (τ N ) 2 but does not depend on P/e. Here, S 2 , N 2 , and P are given by where A, P, Q andR are defined in Sect. 2. The thermodynamic factors I θ and I q weigh the relative contributions of the two scalars to N 2 and P. The quantity P characterizes potential energy of the turbulent flow and may be referred to as the turbulence potential energy. Other definitions of TPE have also been used in the analyses of atmospheric turbulence. A brief discussion is given in Appendix 1. Instead of (τ S) 2 and (τ N ) 2 , the stability functions can be presented in terms of (τ S) 2 [or (τ N ) 2 ] and the gradient Richardson number Ri = N 2 /S 2 . The stability functions in (16)-(18) that result from the solution of Eqs. 13-15 are illbehaved over a certain range of governing parameters. Figure 1 shows F M and F H1 as dependent on Ri and (τ S) 2 at a fixed value of P/e. As seen from the figure, the stability functions have physically meaningless values, i.e., they tend to infinity or become negative over a part of their parameter space. The scalar-flux stability function F H2 (not shown) reveals qualitatively similar ill behaviour as F H1 . The expressions of F M , F H1 , and F H2 are rather lengthy and are not presented here. However, those expressions are drastically simplified in the case of shear-free (and rotation-free) flow, and the ill behaviour of the stability functions can be easily demonstrated analytically. The analysis is presented in Appendix 2. The full linear algebraic system of equations for the Reynolds-stress and scalar-flux components (12 equations for the TKESV scheme) was solved using symbolic algebra software. [The resulting expressions can be received from the first author upon request.] In order to obtain the Reynolds-stress and scalar-flux formulations that remain in force over the entire parameter space, the stability functions should be modified (regularized) in one way or the other. A straightforward approach is to simply impose the upper and lower bounds on the stability functions themselves. Alternatively, constraints may be imposed on the arguments, e.g., on (τ S) 2 and (τ N ) 2 , so that the stability functions cannot become infinite or negative (see, e.g., Mellor and Yamada 1982;Hassid and Galperin 1994). These methods applied within the framework of a one-equation TKE scheme may lead to either the model blow-up or to physically implausible solutions (Helfand and Labraga 1988). One more method is to replace the ill-behaved stability functions with more simple functions that reveal no pathological behaviour. For example, the stability functions of the algebraic turbulence parametrization scheme (level 2), where all second-moment equations, including the TKE equation, are reduced to the diagnostic algebraic formulations, are well-behaved over their entire parameter space. The use of the algebraic-scheme stability functions within the one-equation TKE scheme yields regular solutions. However, an arbitrary replacement of stability functions makes the resulting scheme internally inconsistent. Furthermore, with the level-2 stability functions, the adjustment of the TKE towards the production-dissipation equilibrium state would occur on too short a time scale. In other words, turbulence would respond too quickly to changes in the shear and buoyancy forcing. A more detailed discussion of the above (and some other) regularization methods is given in Helfand and Labraga (1988), see also Hanjalić and Launder (2011) and Lazeroms et al. (2013Lazeroms et al. ( , 2015.

Regularization Procedure of Helfand and Labraga
An appealing way to tackle the problem was proposed by Helfand and Labraga (1988), see also du Vachat (1989) and Helfand and Labraga (1989). These authors analyzed the stability functions of the level-2.5 scheme (one-equation TKE scheme). They considered dry atmosphere, where θ is the only thermodynamic variable that affects buoyancy. As the subsequent analysis of the TKESV scheme is for the moist air, we consider the moist level-2.5 scheme, using thermodynamic variables θ l and q t . Within the framework of the level-2.5 scheme, the scalar-variance and scalar-covariance transport equations are truncated to the diagnostic expressions reflecting the steady-state local balance between the mean-gradient production and the dissipation. The vertical scalar fluxes are given by the following down-gradient formulations, where F 2.5 H is the level-2.5 scalar-flux stability function. The momentum-flux components u 3 u 1 and u 3 u 2 are given by Eqs. 16, where F M is replaced with the level-2.5 stability function F 2.5 M . The level-2.5 stability functions depend on (τ S) 2 and (τ N ) 2 , or, alternatively, on (τ S) 2 [or (τ N ) 2 ] and Ri. The level-2.5 functions suffer from the same deficiencies as the level-3 functions, i.e., they are ill-behaved over a part of their parameter space (Appendix 2). Helfand and Labraga (1988) recast the functions F 2.5 M Ri, τ 2 S 2 and F 2.5 H Ri, τ 2 S 2 in terms of Ri and the ratio e/e e of the actual TKE e = e 2.5 of the level-2.5 scheme (hence the subscript) to the equilibrium value of the TKE e e = e 2 of the level-2 scheme. In order to change the variables Ri, τ 2 S 2 to the variables (Ri, e/e e ), the equality τ 2 S 2 = (τ /τ e ) 2 τ 2 e S 2 = (e e /e) τ 2 e S 2 is used, where τ e is the equilibrium TKE dissipation time scale, i.e., τ computed with e e (see Eq. 11). The value of τ 2 e S 2 for a given Ri is found by solving the equilibrium TKE equation of the level-2 scheme, which is a quadratic equation in τ 2 e S 2 . Since the stability functions F 2 M (Ri) and F 2 H (Ri) are well-behaved over their entire parameter space (any value of Ri), a physically meaningful solution can always be found. Helfand and Labraga (1988) noticed that F 2.5 M and F 2.5 H become pathological in the case of growing turbulence, e/e e = e 2.5 /e 2 < 1, where the TKE production (locally) dominates over the TKE dissipation. Figure 2 illustrates the pathological behaviour of the level-2.5 stability functions on the Ri × e 2.5 /e 2 plane. Note that F 2.5 M and F 2.5 H are computed with the estimates of dimensionless constants given in Sect. 2, which are somewhat different from the estimates used by Helfand and Labraga (1988). This difference does not appreciably affect the result; Fig. 2 of the present paper and Figs. 8 and 9 of Helfand and Labraga (1988) are very similar both qualitatively and quantitatively [within a factor of 2 1/2 C since the expression for fluxes are cast in terms of l(2e) 1/2 in op. cit. and in terms of τ e in the present paper]. As seen from Fig. 2, F 2.5 M and F 2.5 H are well-behaved in the case of equilibrium or decaying turbulence, e/e e ≥ 1, but become pathological at e/e e < 1 (a brief explanation of this behaviour in the regime of shear-free convection is given in Appendix 2).
The failure of the level-2.5 scheme is attributed to the closure assumptions used to develop the scheme, namely, (i) the neglect of the substantial derivative and the third-order transport (diffusion) of the Reynolds stress, scalar fluxes, and scalar variances and covariance, and (ii) the use of simplified linear parametrizations of the pressure-scrambling terms in the Reynolds-stress and scalar-flux equations. These closure assumptions are valid for the flows where turbulence is not far from the stationary, homogeneous, and isotropic state. Helfand and Labraga (1988) argue that it is precisely in the regime of growing TKE, where e/e e < 1, that the turbulence becomes strongly anisotropic (also non-stationary and non-homogeneous) and the above assumptions are too crude, leading to the ill behaviour of the stability functions and to malfunctioning of the turbulence closure scheme. In order to remedy the situation, Helfand and Labraga proposed to account, in an approximate way, for the effects of time rate-of-change, advection and diffusion in those truncated second-moment equations in which these effects have been neglected. In other words, the idea is to restore to the second-moment equations some information that has been lost because of truncation. The following approximation was proposed, Here, M stands for (component of) the Reynolds stress, (component of) the scalar flux, and (within the framework of the level-2.5 scheme) the scalar variances and covariance, d/dt is the substantial derivative that coincides with ∂/∂t within the framework of the boundarylayer approximation, D is the diffusion rate of the respective second-order moment, and G is its production (generation) rate. The production rate G M includes the terms due to mean velocity shear and buoyancy (in the Reynolds-stress and scalar-flux equations), and due to mean scalar gradients (in the scalar-flux, scalar-variance, and scalar-covariance equations). It should be mentioned that, within the framework of the Helfand and Labraga approach, the production rate is defined with due regard for the rapid parts of the pressure-scrambling terms (all but the first terms on the r.h.s. of Eqs. 9 and 10). Equation 24 states that the ratio of the tendency of the second-order moment M due to d/dt and D to the production rate of M is the same as the respective ratio for the TKE. Some comments on the validity of Eq. 24 are in order. Equation 24 is obviously exact in the regime of production-destruction equilibrium. In this regime, the production rate of any moment M, including the TKE, is balanced by the destruction rate of M due to return-toisotropy (M stands for the Reynolds stress or scalar flux) or viscous dissipation (M stands for the scalar variance, scalar covariance, or TKE). In the regime of production-destruction equilibrium, both the r.h.s. and the left-hand side (l.h.s.) of Eq. 24 are zero. Equation 24 is also valid asymptotically if the production rate of the moment M strongly dominates over its destruction rate, and the same holds for the TKE. A transport equation for M can be written in the following form, where M denotes the destruction rate of the moment M. The r.h.s. of Eq. 25 approaches one as the ratio M /G M tends to zero. If the ratio /G e for the TKE also tends to zero, both the r.h.s. and the l.h.s. of Eq. 24 approach one. In the real-world flows, this never holds exactly. However, M /G M 1 may hold approximately following a rapid increase of forcing (by mean scalar gradient, mean velocity shear, or buoyancy), i.e., precisely in the regime of growing turbulence.
Using Eq. 24, the truncated, algebraic equations of the level-2.5 scheme can be modified to approximately account for (mimic) the effects of the time rate-of-change and turbulent diffusion. To this end, approximations dM/dt + D M = α e G M , where α e = (de/dt + D e ) /G e is known from the solution of the TKE equation, are added to the truncated algebraic equations. In substance, these modifications amount to multiplying the production terms in the Reynolds-stress and scalar-flux equations (all but the first terms on the l.h.s. of Eqs. 13-15) and in the scalar-variance and scalar-covariance equations (the first two terms on the r.h.s. of Eq. 3) by a factor 1 − α e . With due regard to the definition of α e , the one-dimensional TKE equation is written as The resulting system of the second-moment equations remains algebraic and linear in the Reynolds stress and scalar fluxes. Its solution (see Helfand and Labraga 1988, for details) yields the down-gradient formulations of the vertical momentum and scalar fluxes, Eqs. 16 and 23, with the following expressions of the stability functions of the regularized level-2.5 scheme,   where the subscript "r" indicates the modified (regularized) stability functions. The regularized functions given by Eq. 27 are used in the regime of growing turbulence, where the actual value of TKE e 2.5 of the level-2.5 scheme is less than the equilibrium TKE e 2 of the level-2 scheme. In the regime of equilibrium or decaying turbulence, the original, non-modified stability functions of the level-2.5 scheme are used that are well-behaved at e 2.5 /e 2 ≥ 1 and pose no problem. Note that matching of regularized and non-regularized stability functions is smooth since F 2.5 M = F 2.5 Mr = F 2 M at e 2.5 /e 2 = 1 (similar for the scalar-flux stability functions). Recall that the stability functions F 2 M and F 2 H of the level-2 algebraic scheme are well-behaved over their entire parameter space. Figure 3 illustrates the regularized stability functions of the level-2.5 scheme that reveal no pathological behaviour.
A note on the critical Richardson number is in order. The critical Richardson number Ri cr is defined here as the value of Ri beyond which the shear-generated turbulence cannot be maintained against the destructive action of gravity, i.e., the TKE (as well as the TPE) is zero at Ri ≥ Ri cr . The Richardson number is the local quantity, and any turbulence scheme that accounts for non-local effects in space (third-order transport) and/or time (nonstationarity) should be able to maintain turbulence at supercritical values of Ri. Clearly, there is no critical Richardson number within the framework of the level-2.5 scheme, where the effects of non-stationarity and third-order transport in the TKE equation are accounted for. This is different for the level-2 scheme, where the non-local effects are neglected, the turbulence energy is always in the equilibrium state, and a finite value of Ri cr does exist. With the estimates of dimensionless constants given in Sect. 2, Ri cr = 0.96 (cf. Cheng et al. 2002). As Ri approaches Ri cr , the equilibrium TKE e 2 of the level-2 scheme approaches zero and so do the stability functions F 2 M and F 2 H . At Ri ≥ Ri cr , a non-trivial solution of the level-2 system does not exist. In this part of the domain, e 2 = 0, e 2.5 /e 2 ≥ 1, and the original, non-regularized stability functions are applied (which are well-behaved there). Thin white stripes near the right margins of the plots in Figs. 2 and 3 merely indicate that F 2.5 M and F 2.5 H cannot be shown on the Ri × e 2.5 /e 2 plane beyond the critical Richardson number of the level-2 scheme.

Regularized Stability Functions of the TKESV Scheme
The regularization method of Helfand and Labraga (1988) is now applied to the TKESV (level-3) scheme. Using Eq. 24, the Reynolds-stress and scalar-flux equations are modified, and the resulting system of algebraic equations for a i j , u i θ l and u i q t is solved. In the general case of stratified shear flow, the derivations are cumbersome and are omitted here. The derivations for the shear-free flow are presented in Appendix 3. The following expressions of the stability functions in Eqs. 16-18 are obtained, where "3" refers to the level-3 (i.e., to the TKESV) scheme, and "r" indicates the regularized stability functions. The scheme referred to as "2.p" utilizes the equilibrium TKE equation but carries transport equations for the scalar variances and scalar covariance. That is, the TKE is determined diagnostically using the steady-state production-dissipation balance, whereas the scalar variances and covariance (and hence the TPE) are determined prognostically with due regard for the third-order turbulent transport. 2 The equilibrium TKE equation of the level-2.p scheme is cubic in e e = e 2.p (cf. the level-2 scheme, where the equilibrium TKE equation is quadratic in e e = e 2 ). Figure 4 shows the stability functions of the level-2.p scheme as dependent on Ri and s 2 = P/ e τ 2 S 2 2 , where s 2 is the dimensionless TPE normalized so that it does not contain the TKE. The momentum-flux stability function F 2.p M is well-behaved over its entire parameter space. The scalar-flux functions F 2.p H1 and F 2.p H2 are also well-behaved over most of their parameter space, but they become pathological with strong static stability (large values of Ri) where the TKE is small. For most practical purposes, a simple clipping should be sufficient to remedy the situation. 2.p H2 plots shows the equilibrium value of s 2 as function of Ri, i.e., the TPE of the level-2 scheme A more physically appealing remedy is based on the observation that the ill behaviour is caused by the factor that appears in the denominator of the expressions of both scalar-flux stability functions.
Here, C * is a dimensionless constant (a combination of model constants given in Sect. 2), and the second expression on the r.h.s. is obtained using (11). As the static stability becomes strong, the equilibrium e decreases, τ increases, and B eventually becomes negative, leading to the ill behaviour of F 2.p H1 and F 2.p H2 . However, this happens only if the turbulence length scale l is taken to be independent of the density (buoyancy) stratification. This assumption is in fact unrealistic where the stratification is stable. In the TKESV scheme, a stability correction to the Blackadar-type interpolation formula for the turbulence length scale is used (e.g. Stull 1973;Zeman and Tennekes 1977). In the limit of strongly stable stratification, the length scale behaves as where C N is a dimensionless constant of order 1. Note that an advanced non-local formulation of the length scale proposed by Bougeault and Lacarrère (1989) also accounts for the effect of stable density stratification on l; it reduces to Eq. 30 in a particular case of height-constant N . With due regard for (30), the expression of B becomes As the static stability increases, the second term on the r.h.s. of Eq. 31 decreases, B remains positive, and no pathological behaviour of the scalar-flux stability functions takes place. Thus, it is necessary to use a stability-dependent formulation of the turbulence length (time) scale in order to ensure that the scalar-flux stability functions of the level-2.p scheme are well-behaved over their entire parameter space. Figure 5 shows the original, non-modified stability functions F 3 M and F 3 H1 of the level-3 scheme (F 3 H2 is not shown) on the Ri × e 3 /e 2.p plane at a fixed value of s 2 = 0.01. Again, the stability functions are well-behaved in the case of equilibrium or decaying turbulence, but become pathological in the case of growing turbulence. Although the ill behaviour is only seen over a part of growing-turbulence domain, we apply the regularized functions at e 3 /e 2.p < 1. In order to apply the regularized stability functions over a part of the growingturbulence domain, additional objective criteria are required to separate out the "dangerous" part of the e 3 /e 2.p < 1 domain (for example, the realizability constraints considered in Sect. 9 can be used). Note also that matching F 3 M and F 3 H1 with F 3 Mr and F 3 Hr , respectively, at any value different from e 3 /e 2.p = 1 would result in discontinuous stability functions. Figure 6 shows the regularized stability functions of the level-3 scheme. As can be seen, the stability functions reveal no pathological behaviour (and are continuous). Note that the regularized, well-behaved expressions are obtained for all components of the Reynolds stress and scalar fluxes, not only for u 3 u 1 , u 3 u 2 , u 3 θ l , and u 3 q t . Those expressions are not presented here.
It is interesting to draw analogies between Eqs. 27 and 28. As Eq. 27 shows, the regularized level-2.5 stability functions are expressed in terms of the level-2 stability functions. For the level-2.5 scheme, the level-2 scheme is the nearest lower-level scheme, whose stability functions are well-behaved. Likewise, for the level-3 scheme, the level-2.p scheme is the nearest lower-level scheme with well-behaved stability functions (provided the stabilitydependent formulation of the length/time scale is used). The level-2 and the level-2.p schemes differ significantly in the way they treat the scalar variances and covariance. However, both schemes utilize a steady-state equilibrium TKE equation, suggesting that it is the use of the TKE transport equation in combination with the simplified algebraic formulations of the Reynolds stress and scalar fluxes that causes ill behaviour of truncated turbulence closure schemes.
Note an important factor e 3 /e 2.p 1/2 in the expressions (28) for the regularized level-3 stability functions. It is this factor that takes care of a gradual transition of the TKE towards the production-dissipation equilibrium. Dropping the factor e 3 /e 2.p 1/2 in Eq. 28 would result in a too quick response of turbulence to changes in the forcing and in overestimation of momentum and scalar fluxes in the regime of growing TKE. The factor (e 2.5 /e 2 ) 1/2 in Eq. 27 has a similar effect on the behaviour of the level-2.5 scheme. Neither the level-3 scheme nor the level-2.p scheme has a critical Richardson number. Within the framework of the level-3 scheme, the effects of non-stationarity and third-order transport are accounted for in both the TKE and the scalar (co)variance equations. Within the In the case of growing turbulence, e 3 /e 2.p < 1, the stability functions are computed from Eq. 28 framework of the level-2.p scheme, the TKE transport equation is truncated to the diagnostic expression reflecting the steady-state local balance between production and dissipation. However, the TKE can be maintained at an arbitrary large value of Ri at the expense of the TPE, whose local value can be different from zero due to the effects of non-stationarity and third-order transport of scalar (co)variances.
The various truncated second-order closure schemes are summarized pictorially in Fig. 7 on the e/e e × P/P e plane, where P e is the equilibrium TPE resulting from the productiondestruction balance of the truncated equations for the scalar variances and scalar covariance. Within the framework of the level-2 scheme, all second-moment equations are reduced to diagnostic algebraic expressions. That is, all second-order moments are in the state of local production-destruction equilibrium. The level-2 scheme corresponds to a single point (e/e e = 1, P/P e = 1) on the e/e e × P/P e plane. A horizontal line P/P e = 1 depicts the level-2.5 scheme, where the TPE is in the state of local production-destruction equilibrium, whereas the TKE can take on any value different from e e due to the effects of time rateof-change, advection and diffusion. A vertical line e/e e = 1 corresponds to the level-2.p scheme, where the TKE is in the production-destruction equilibrium, whereas the TPE can

Relation to Non-linear Parametrizations of the Pressure-Scrambling Terms
The idea of the regularization method of Helfand and Labraga (1988) is to restore to the second-moment equations some information that has been lost because of truncation. Equation 24 is used to approximately account for the effects of the time rate-of-change and turbulent diffusion in the Reynolds-stress and scalar-flux (and, in the case of level-2.5 scheme, in the scalar-variance and scalar-covariance) equations. An alternative a posteriori interpretation can be suggested, however. As discussed above, the use of Eq. 24 amounts to multiplying the production terms in the truncated equations by a factor 1 − α e , where α e = (de/dt + D e ) /G e depends on time t and vertical coordinate x 3 . Consider, for example, the buoyancy production term in the temperature-flux equation 14, where C θ b is a constant that stems from a linear parametrization of the buoyancy contribution to the pressure-scrambling term. With the estimate of C θ b = 1/3, the isotropic turbulence constraint is satisfied. Recall that the use of simplified linear parametrizations of the pressurescrambling terms is identified as one reason for failure of truncated second-order schemes. If Eq. 24 is utilized to develop regularized, well-behaved expressions for fluxes, the buoyancy production term in the temperature-flux equation becomes where C θ b * = C θ b + α e 1 − C θ b is no longer a constant, but a function of t and x 3 . It varies from C θ b * = C θ b where α e = 0 (which recovers the original, non-regularized formulation) to C θ b * = 1 where α e = 1. Then, the Helfand and Labraga regularization procedure can be viewed as a parametrization assumption that is in effect similar to the use of non-linear parametrizations of the pressure-scrambling effects.
A turbulence closure proposed by Craft et al. (1996) (see also Mironov 2001) is an illustrative example of a closure based on non-linear parametrizations of the pressure-scrambling terms. The closure of Craft et al. (1996) incorporates a non-linear formulation, where C θ b * is a function of the departure-from-isotropy tensor. It satisfies both the isotropic constraint and the so-called two-component limit (TCL) constraint of strongly anisotropic turbulence. The TCL is the limit that turbulence approaches when the velocity component (both mean velocity and turbulent fluctuation) in one direction vanishes. This occurs, for example, near the rigid boundary, where the velocity component normal to the boundary vanishes as the boundary is approached. It also occurs in stably stratified layers, where the velocity component aligned with the vector of gravity vanishes as the stratification becomes strong. In order to satisfy the TCL constraint, C θ b * should tend to one as the velocity component aligned with the vector of gravity tends to zero. Otherwise, spurious generation of the temperature flux by the buoyancy forces takes place. Clearly, it is impossible to satisfy both the isotropic and the TCL constraints if linear parametrizations of the pressure-scrambling terms with a constant C θ b * are used. The formulation (33) is not the same as the TCL formulation. For example, α e → 1 and hence C θ b * → 1 may occur in the case of quickly growing nearly isotropic turbulence, which is clearly not the TCL. However, if the quickly growing turbulence is nearly two-component, Eq. 33 would affect the buoyancy terms in the scalar-flux equations in a similar way as the TCL formulation.

Relation to Weak Non-equilibrium Hypothesis
The so-called weak non-equilibrium hypothesis (Rodi 1976;Gibson and Launder 1976;Hanjalić and Launder 2011) is often used to develop turbulence closure schemes of wider applicability than the schemes, where the time rate-of-change, advection, and third-order transport of the Reynolds stress and scalar fluxes are entirely neglected. Within the framework of the weak non-equilibrium hypothesis, the above effects are accounted for in an approximate way. In this section, the relation between the Helfand and Labraga (1988) approach and the weak non-equilibrium hypothesis is discussed.
The weak non-equilibrium hypothesis (Rodi 1976) states that during the flow evolution the sum of the substantial derivative (equal to ∂/∂t within the framework of the boundary-layer approximation) and the diffusion of the dimensionless departure-from-isotropy tensor a i j /e is zero. This yields which states that the ratio of the tendency of the Reynolds-stress component due to d/dt and D to that component itself is the same as the respective ratio for the TKE. This is at variance with the Helfand and Labraga approach, where a parametrization assumption is made about the ratio of the Reynolds-stress component tendency to its production rate, see Eq. 24. The application of the weak non-equilibrium hypothesis to the scalar (temperature) flux (Gibson and Launder 1976) assumes that the normalized scalar flux (the u i -θ correlation coefficient) θ u i / θ 2 e 1/2 remains unchanged as the flow evolves. This yields Consider a simplified case of shear-free temperature-stratified flow, where the potential temperature is the only thermodynamic variable that affects buoyancy. Using Eqs. 34 and 35, we obtain the following equations for a 33 and u 3 θ , where the quantities γ e = de/dt + D e e and γ θ = d θ 2 /dt + D θθ θ 2 have the dimensions of reciprocal time and are taken to be known from the solutions of the TKE and the θ 2 equations. If γ e > 0, which corresponds to the regime of growing TKE, the use of the weak non-equilibrium hypothesis results in the reduction of the Reynolds-stress relaxation (return-to-isotropy) time scale by a factor of 1 + γ e τ /C u t −1 (as compared to the non-modified system). The reduction of the scalar-flux relaxation time scale is by a factor of 1 + (γ e + γ θ ) τ /2C θ t −1 (provided γ e + γ θ is positive). The application of the Helfand and Labraga procedure effectively results in the reduction of the Reynolds-stress and scalar-flux production rates by a factor of 1 − α e , where α e = (de/dt + D e ) /G e is known from the solution of the TKE equation. A comparison of Eqs. 36 and 37 with Eqs. 59 and 60, respectively, suggests that in the regime of growing turbulence the two approaches have similar effect on the Reynolds stress and scalar fluxes. If γ e = γ θ and in addition C θ t = C u t , the two approaches are equivalent. Both approaches yield Eqs. 59-62 where 1 − α e = 1 + γ e τ /C u t −1 . It should be pointed out, however, that the weak non-equilibrium hypothesis is usually applied not only in the regime of growing turbulence, but also during the turbulence decay, where γ e < 0 and/or γ e + γ θ < 0. As seen from Eqs. 36 and 37, large negative γ e and/or γ e + γ θ may yield infinite or negative relaxation time scale(s), i.e., a spurious result. Then, further modifications of the Reynolds-stress and scalar-flux formulations are necessary to avoid their ill behaviour (see, e.g., Lazeroms et al. 2013Lazeroms et al. , 2015Zeli et al. 2019).

Realizability of Turbulence Closures and the Problem of Moments
Until now, we used the notion of realizability in a somewhat loose sense, being concentrated mostly on the physically meaningful boundedness of the quantities in question in general and on the non-negativeness of the non-negatively defined quantities in particular. Whereas the non-negativeness is a well-defined criterion, the physically meaningful boundedness had a rather diffuse meaning. Observing continuous unbounded growth of the stability functions, we were unable to argue without any further knowledge where exactly they still have physically meaningful values and where already not. The notion of realizability can, however, be defined mathematically more precisely as the assurance of the existence of the probability distribution of a random variable with prescribed statistical moments (du Vachat 1977). This is also known as the moments problem (Shohat and Tamarkin 1943). The solution of the moments problem can be represented as a precise and complete set of constraints that the modelled statistical moments should obey. If the modelled statistical moments do not obey those necessary constraints, then there cannot exist any probability distribution or, equivalently, any random variable possessing those moments. In such a case, one may justifiably say that the model produces an unphysical solution and is non-realizable. As far as the stability functions are concerned, the solution to the moments problem does not tell us what their values should be. However, it tells what the values of stability functions may not take on.
A general solution of the moments problem is deduced from the condition that the probability density function should be non-negative. For a multi-variate random variable, it consists in the requirement that a matrix of a special kind containing all statistical moments should be positive semi-definite. It is not easy to present that matrix in a generic and compact way. By way of illustration, we denote the components of the multi-variate random variable as {a, b, ..., z} and present the first row of the matrix which should contain all moments in the ascending order, starting with the zeroth-order moment. That is 1 a b · · · z a 2 ab · · · az b 2 · · · bz · · · z 2 a 3 a 2 b · · · The same sequence as in (39) constitutes the first column. The element in the ith row and the jth column of the matrix is determined by multiplying the (not averaged) quantities standing in the first row of the jth column and in the first column of the ith row and by subsequent averaging. Note that in our application the variables are centered (fluctuations about the mean values are considered). Hence a = 0, and the same holds for all other variables. In accordance with the Sylvester's criterion (see, e.g., Shafarevich and Remizov 2013), positive semi-definiteness of a matrix is equivalent to non-negativeness of all its principal minors (the determinants of smaller quadratic matrices obtained by crossing out all possible subsets of columns and rows with the same number). This requirement leads to an infinite chain of inequalities relating the moments of all orders to each other. If, however, positive definiteness of the matrix in question rather than its positive semi-definiteness were required (e.g., if we would not allow variances to be exactly zero), then it would be sufficient to demand that all leading principal minors be positive. Recall that the leading principal minors are the determinants of the smaller quadratic matrices (of all orders), which are merely the upper left corners of the original matrix. It may seem that the order of the variables in the matrix matters and different variables will not have equal rights. This is not the case, however. Crossing out all columns and rows except for m first columns and rows (m = 1, . . . ) is in fact equivalent, in the case of positive definiteness, to crossing out all possible sub-sets of columns and rows. This can readily be seen from the following example. Consider the second-order and the third-order leading principal minors of the matrix ⎛ ⎜ ⎜ ⎝ 1 0 0 · · · 0 a 2 ab · · · 0 ab b 2 · · · · · · · · · · · · · · · ⎞ ⎟ ⎟ ⎠ .
The requirements a 2 > 0 (second-order leading principal minor is positive) and a 2 b 2 − ab 2 > 0 (third-order leading principal minor is positive) immediately yield b 2 > 0. However, this is no longer true if a 2 is allowed to be zero. Indeed, if a 2 = 0 and ab = 0, a 2 b 2 − ab 2 ≥ 0 is satisfied at any value of b 2 . This point was apparently overlooked in Schumann (1977), where it was pointed out that non-negative leading principal minors are sufficient for the matrix to be positive semi-definite. Now we consider moist turbulence and limit ourselves to the second-order moments. For five random fields (three components of u i , θ l and q t ), the matrix has the following form ⎛ where the row and the column containing only one and zeros (the variables are centered) are omitted. The requirement of non-negative definiteness of the matrix (41) generates five inequalities from the first-order minors, ten inequalities from the second-order minors, ten inequalities from the third-order minors, five inequalities from the fourth-order minors, and one inequality expressing the non-negativeness of the determinant of the entire matrix. Those inequalities impose powerful constraints on the second-order moments and hence on the stability functions. Consideration of all 21 constraints is rather cumbersome. By way of illustration, we analyze a simplified case of dry, shear-free flow considered in Appendices 2 and 3. Note though that the white area in Fig. 5 where the stability functions are ill-behaved widens as Ri decreases, suggesting that it is this seemingly simple shear-free convection case that is potentially the most dangerous.
In the dry, shear-free flow, u 1 q t = u 2 q t = u 3 q t = 0, q 2 t = 0, q t θ l = 0, u 1 u 2 = u 1 u 3 = u 2 u 3 = 0, u 1 θ l = u 2 θ l = 0, and θ l = θ . The matrix (41) is simplified to give ⎛ The requirement of non-negative definiteness of the matrix (42) yields the following inequalities, which state that the variances are non-negative and the magnitude of the temperature-flux correlation coefficient does not exceed 1. In what follows, we derive the constraints these inequalities impose on the second-order moments of the regularized TKESV scheme.
For the level-3 schemes, the temperature variance θ 2 is an independent prognostic variable that may be regarded as an external parameter for the algebraic equation system for the second-order moments. Non-negativeness of θ 2 must be ensured by its transport equation. Since in shear-free convection regime u 3 θ ≥ 0, Eq. 59 shows that the vertical-velocity variance is non-negative as well.
Substituting the expression of u 2 3 from (59) and (67) into (44), then eliminating u 3 θ with the help of (61) and (67), we obtain the following constraint, Equation 45 shows that, although the stability functions of the regularized TKESV scheme are well-behaved, the scheme can still be non-realizable at certain values of TPE and/or TKE. In other words, realizability in strict mathematical sense imposes additional, more rigid constraints as compared to the requirements that some quantities are non-negative and finite.
Although it may seem different at a glance, the realizability constraint (45) is in fact very mild as far as practical applications are concerned. The above inequality may be violated where either e is (very) small as compared to e e , or P is (very) small as compared to e (or both conditions hold). This is rarely encountered in practice. Numerical experiments with the TKESV scheme (results not shown) indicate that the inequality (45) is virtually never violated. The violation may occur over a certain adjustment period if TPE is poorly initialized, e.g., if e and e e are both positive but P is set equal to zero. However, mean-gradient production and vertical diffusion will ensure a quick build-up of the scalar variances, and hence of the TPE, if the TKE is non-zero. Still there is no guarantee that a turbulence scheme will never fall into the non-realizable state, but it is highly desirable that the scheme is able to handle any theoretically allowable values of the mean and turbulence quantities. We therefore propose to apply a simple clipping to the scalar (co)variances, or to the scalar fluxes, so that the realizability requirements are satisfied.
The realizability constraints given by the first two members of Eq. 43 state that the horizontal velocity variances are non-negative. Using Eq. 62 with due regard for (60), (61) and (67), we obtain after some algebra This relation demonstrates the usefulness of the realizability considerations in developing turbulence models. The realizability constraints may help limit the allowable values of the model constants. With the values of C u t and C u b used within the TKESV scheme (see Sect. 2), Eq. 46 is satisfied.
The above analysis is limited to the second-order turbulence moments. Principle minors containing higher-order moments yield further powerful constraints that help limit disposable parameters of turbulence models and help control the model behaviour during the integration. Illustrative examples are the relations between skewness and kurtosis of fluctuating fields (e.g., André et al. 1976;Gryanik et al. 2005). Consideration of realizability constraints on high-order turbulence moments is beyond the scope of the present paper.
The problem of realizability of truncated second-order closure schemes for moist atmospheric boundary-layer turbulence is addressed through the consideration of the so-called stability functions. The emphasis is on the TKE-scalar variance closure schemes (the level-3 schemes in the nomenclature of Mellor and Yamada 1974) that carry prognostic transport equations for the TKE and for the (co)variances of scalar quantities (e.g., liquid water potential temperature and total water specific humidity). The stability functions were introduced Yamada 1974, 1982;Yamada 1977) to represent the formulations of the Reynolds stress and scalar fluxes in a compact form. They appear within the framework of truncated closure schemes, where (i) the Reynolds-stress and scalar-flux equations (and, within the framework of one-equation TKE schemes, also equations for scalar variances and covariance) are reduced to the diagnostic algebraic formulations by neglecting the substantial derivatives and the third-order transport terms, and (ii) linear parametrizations (in the second-order moments involved) of the pressure-scrambling terms are used. The stability functions are ill-behaved over a certain range of governing parameters, e.g., mean velocity shear and buoyancy gradient. They tend to infinity or become negative in the regime of growing turbulence, where the actual value of the TKE is smaller than the equilibrium TKE corresponding to the steadystate production-dissipation balance. The ill-behaved stability functions inevitably lead to malfunctioning of turbulence closure schemes.
Regularized stability functions for the TKESV closure scheme (Mironov and Machulskaya 2017) are developed that reveal no pathological behaviour over their entire parameter space (provided a stability-dependent formulation of the turbulence length/time scale is used). To this end, the approach (regularization procedure) of Helfand and Labraga (1988) is invoked. The approach, originally developed for the one-equation TKE schemes, is extended to the TKESV scheme. The idea is to account, in an approximate way, for the effects of time rate-ofchange, advection, and diffusion in those truncated second-moment equations in which these effects have been neglected. That is, some information that has been lost because of truncation is restored to the second-moment equations. The physical meaning of the regularization procedure and its relation to the non-linear, more physically sound parametrizations of the pressure-scrambling terms are discussed. It is shown that the Helfand and Labraga approach can be viewed as a parametrization assumption that is in effect similar to the use of non-linear parametrizations of the pressure-scrambling effects.
The relation between the Helfand and Labraga (1988) approach and the weak nonequilibrium hypothesis (Rodi 1976;Gibson and Launder 1976;Hanjalić and Launder 2011) often used to develop turbulence closure schemes is considered. It is demonstrated that, in the regime of growing turbulence, the use of the Helfand and Labraga procedure effectively results in the reduction of the Reynolds-stress and scalar-flux production rates, whereas the use of the weak non-equilibrium hypothesis results in the reduction of the Reynolds-stress and scalar-flux relaxation (return-to-isotropy) time scales. Where the TKE is smaller than its equilibrium value, the two approaches have similar effect on the Reynolds stress and scalar fluxes.
Finally, realizability of turbulence closures is considered within a more general and mathematically more rigorous framework of the moments problem of the probability theory. Realizability in strict mathematical sense imposes additional, more rigid constraints as compared to the requirements that some quantities handled by a turbulence scheme are merely non-negative and finite. The realizability constraints help control the model behaviour by imposing limits on the values of turbulence moments computed by a closure scheme. They can also limit disposable parameters of turbulence schemes and are, therefore, of considerable help for the model development purposes.
The results from the present study are useful for numerical weather prediction, climate modelling, and related applications. Numerical weather prediction and climate models require physically sound and, importantly, robust turbulence parametrization schemes. Those schemes should be applied around the clock and (within global models) over the entire atmosphere, they should be able to handle any atmospheric turbulence regime. Unless regularized stability functions are applied, the TKESV scheme, as well as similar rather advanced turbulence parametrization schemes, cannot be utilised. Ill-behaved stability functions will inevitably lead to an atmospheric model blow-up. The use of regularized stability functions developed herein makes the schemes much more robust.
The TKESV closure scheme of Mironov and Machulskaya (2017) with the regularized, well-behaved stability functions reported in the present study is implemented into the numerical weather prediction model ICON (ICOsahedral Nonhydrostatic modelling framework). Testing through numerical experiments is underway at the German Weather Service; results will be reported later.
The quantity P defined by Eqs. 21 and 22 naturally appears in the equations for the Reynolds stress and scalar fluxes. It is characteristic of the potential energy of turbulent flow and can, therefore, be referred to as the turbulence potential energy (TPE). By way of illustration, consider the temperature-stratified fluid, where potential temperature θ is the only thermodynamic variable that affects buoyancy. Then, q l = 0, q t = 0, θ l = θ , I θ = 1, N 2 = β∂ θ /∂ x 3 , and the definition of potential energy (21) is simplified to give The definition of TPE is not unique, however. For example, the following definition of TPE has been used in the analyses of atmospheric turbulence (see e.g. Zilitinkevich et al. 2007;Mauritsen et al. 2007), A transport equation for P N can be readily derived using the transport equation for buoyancy (potential temperature) variance. Assuming that N 2 varies slowly in space and time, and using the boundary-layer approximation, we obtain In the boundary-layer approximation, the transport equation for the TKE reads The buoyancy flux β u 3 θ appears with the opposite signs in Eqs. 49 and 50. It describes the rate of conversion of the TPE to the TKE that occurs in unstable density (buoyancy) stratification, and vice versa where the stratification is stable. The use of TPE defined through Eq. 48 is advantageous as it makes the TPE transport equation particularly convenient and facilitates the analysis of turbulence energetics. For example, Eqs. 49 and 50 can be added, leading to the total (TPE + TKE) turbulence energy equation that has been used by some researchers in the analyses of stably-stratified turbulent flows (e.g. Zilitinkevich et al. 2007Zilitinkevich et al. , 2009Zilitinkevich et al. , 2013Mauritsen et al. 2007).
Note, however, that the TPE defined through Eq. 48 is a well-defined quantity (nonnegative and finite) only if the flow is stably stratified (N 2 > 0). The quantity P N bears a close analogy to the available potential energy defined as a part of the total potential energy of the stratified flow that can be converted to kinetic energy (see, e.g., Vallis 2006, for a comprehensive discussion). Where the buoyancy stratification is neutral (N 2 = 0) or unstable (N 2 < 0), P N is infinite or negative and is not really convenient to use. For the Earth's atmosphere, where the buoyancy stratification is due to both temperature and humidity, and the thermodynamics is strongly complicated by phase changes, a TPE transport equation as simple and elegant as Eq. 49 is difficult to derive. It is therefore advantageous to work with the variance and covariance equations for the scalar quantities that can be derived from the first principles in a fairly straightforward way. The quantity P defined by Eqs. 21 and 22, or by Eq. 47 in the case of temperature-stratified flow, is merely a diagnostic quantity that naturally appears in the equations for the Reynolds stress and scalar fluxes. It has the physical meaning or turbulence potential energy but a separate transport equation for P is actually not needed.
where u 3 θ v is expressed through u 3 θ l and u 3 q t using Eq. 12. The vertical scalar fluxes are given by where N 2 and P are given by Eqs. 20-22. In the shear-free limit, F H2 is merely a constant and F H1 is a function of τ 2 N 2 and P/e. There is no problem with the stability function F H1 in the case of stable buoyancy stratification, where N 2 is positive. If the stratification is unstable, N 2 is negative, and the expression in the denominator (first set of square brackets) on the r.h.s. of Eq. 55 may approach zero or become negative, leading to physically meaningless values of F H1 . A similar problem with the stability functions is encountered within the framework of a one-equation closure scheme that carries prognostic transport equation for the TKE, whereas all other second-moment equations, including the scalar-variance and scalar-covariance equations, are reduced to diagnostic algebraic expressions. Consider a temperature-stratified flow where θ l = θ . Within the framework of the one-equation TKE scheme, the following downgradient formulation for the vertical potential-temperature flux holds, As Eq. 58 suggests, F H may become infinite or negative in convective conditions, where N 2 < 0. The above considerations help explain why the stability functions are well-behaved in the case of equilibrium or decaying turbulence but may become pathological in the case of growing turbulence (see Sect. 5 and 6). Recall that no ill behaviour is encountered if the TKE is in the state of production-dissipation equilibrium (level-2 and level-2.p schemes). That is, the value of τ 2 N 2 corresponding to the equilibrium TKE e e does not lead to ill behaviour of F H1 and F H given by Eqs. 55 and 58, respectively. As Eq. 11 suggests, the larger the TKE, the smaller the dissipation time scale. Then, in the case of decaying turbulence, where e > e e by definition, the value of τ 2 N 2 is smaller in magnitude than in the equilibrium case, and no pathology is encountered. As the term "growing" implies, the TKE is small as compared to its equilibrium value. Then, τ 2 N 2 in Eqs. 55 and 58 may become large in magnitude (N 2 is negative in convective conditions), leading to infinite or negative stability functions. Now substituting τ (1 − α e ) for τ e in Eq. 66 and using Eq. 67 to express 1 − α e through the equilibrium and non-equilibrium TKE, e e = e 2.p and e = e 3 , respectively, yields Eqs. 28b and 28c for the regularized stability functions of the level-3 (TKESV) scheme, where the stability functions of the level-2.p scheme given by Eq. 64 are well-behaved.
A consideration of the general case of stratified shear flow is principally the same as of the shear-free flow, but the derivations are more lengthy.