## Abstract

We present a turbulent kinetic energy (TKE) closure scheme for the stably stratified atmosphere in which the mixing lengths for momentum and heat are not parameterized in the same manner. The key difference is that, while the mixing length for heat tends toward the stability independent mixing length for momentum in neutrally stratified conditions, it tends toward one based on the Brunt–Väisälä time scale and square root of the TKE in the limit of large stability. This enables a unique steady-state solution for TKE to be obtained, which we demonstrate would otherwise be impossible if the mixing lengths were the same. Despite the model’s relative simplicity, it is shown to perform reasonably well with observational data from the 1999 Cooperative Atmosphere–Surface Exchange Study (CASES-99) using commonly employed model constants. Analyzing the scaling behavior of the nondimensional velocity and potential temperature gradients, or of the stability (correction) functions, reveals that for large stability the present model scales in the same manner as the first-order operational scheme of Viterbo et al. Alternatively, it appears as a blend of two cases of the TKE closure scheme of Baas et al. Critically, because a unique steady-state TKE can be obtained, the present model avoids the nonphysical behavior identified in one of the cases of Baas et al.

## 1. Introduction

The stably stratified planetary boundary layer (SBL) forms when the planetary surface is cooler than the air above, frequently occurring at night or in polar regions. The thermal stratification acts to suppress turbulent motions, due to the effort in drawing up heavier (cooler) fluid from below and pulling down lighter (warmer) fluid from above. The degree of thermal stratification is traditionally classified as either weakly stable, associated with moderate winds and sustained turbulence, or very stable, associated with weak winds and weak or intermittent turbulence. The very stable case is also sensitive to other factors such as gravity waves, longwave radiation, and local topography (Nieuwstadt 1984; Mahrt 1999).

Despite the prevalence of the SBL and decades of research, it remains particularly challenging for atmospheric models to realistically represent (Viterbo et al. 1999; Cuxart et al. 2006). Inaccurate representations of the SBL affect the predicted momentum and heat fluxes at the planetary surface, which can go on to have significant large-scale influences in both weather and climate forecasting (Cuxart et al. 2006). Models therefore need to be accurate and robust, both in terms of their operational performance but also in their ability to faithfully reproduce fundamental physical processes and states.

In general, atmospheric turbulence parameterizations use the gradient-flux approach, in which vertical turbulent fluxes are related to the vertical gradient through an eddy viscosity or diffusivity. For momentum (subscript *m*) and heat (subscript *h*) these are

where *u* and *υ* are the mean horizontal wind components, *w* is the vertical wind, *θ*_{υ} is the mean virtual potential temperature, and primes denote fluctuations. The eddy viscosity *K*_{m} and eddy diffusivity *K*_{h} must then be parameterized.

To aid the interpretation of common parameterizations for *K*_{m} and *K*_{h} discussed below, we briefly introduce the local-scaling theory of Nieuwstadt (1984). This theory states that nondimensional quantities within the stable boundary layer are only functions of *ζ* ≡ *z*/Λ, where Λ is the *local* Obukhov length, defined as

where $u*2=\u2061(w\u2032u\u2032\xaf2+w\u2032\upsilon \u2032\xaf2)1/2$ is the local friction velocity and $\theta *=\u2212w\u2032\theta \upsilon \u2032\xaf/u*$ is the local friction potential temperature, both of which vary with height, and *κ* = 0.4 is the von Kármán constant. This local-scaling theory represents a generalization of Monin–Obukhov similarity theory (MOST), wherein MOST is only valid within the near-ground surface layer and uses the *surface* Obukhov length *L*. Essentially, local-scaling theory states that the nondimensional velocity and virtual potential temperature gradients *ϕ*_{m,h} are solely functions of *ζ*, where

Defining the shear rate *S*^{2} = (∂*u*/∂*z*)^{2} + (∂*υ*/∂*z*)^{2} and Brunt–Väisälä (BV) frequency *N*^{2} = (*g*/Θ_{υ})(∂*θ*_{υ}/∂*z*) enables these nondimensional gradient functions to be related to the gradient and flux Richardson numbers as

The turbulent Prandtl number is then given as

First-order closure schemes are commonly used in operational forecast models, wherein *K*_{m,h} is parameterized as $Km,h=lm,h2SFm,h\u2061(Rig)$, where *l*_{m,h} is the mixing length for momentum or heat. The term *F*_{m,h} = 1/(*ϕ*_{m}*ϕ*_{m,h}) is a stability (correction) function that typically reduces with Richardson number to account for the reduction of turbulent motions with increasing stability. These first-order schemes therefore only require solving the prognostic equations for the mean variables. However, experimental data for *ϕ*_{m,h} and therefore *F*_{m,h} show large scatter for high stabilities (see, e.g., Högström 1988; Andreas 2002; Beare et al. 2006), as well as suffering from self-correlation issues (Hicks 1978; Baas et al. 2006). This leads to a wide range of semiempirical forms of *F*_{m,h} suggested in the literature for first-order schemes. Some impose a critical Richardson number of approximately 0.2 above which the flow relaminarizes and turbulence ceases entirely (e.g., Businger et al. 1971; Dyer 1974). While this behavior agrees with MOST in the surface layer, it is contrary to observational data that shows turbulence persisting for Ri_{g} > 1 (Galperin et al. 2007; Huang and Bou-Zeid 2013; Mahrt 2014, and references therein). Moreover, the lack of turbulent mixing can lead to runaway cooling and a decoupling between the near-surface state and that higher up in the SBL (Derbyshire 1999). This motivated so-called sharp forms for *F*_{m,h}, which enables turbulence to remain nonnegligible for Ri_{g} → ∞ (e.g., King et al. 2001). Operational models often require even more mixing than these forms of *F*_{m,h} provide to avoid the decoupling behavior, leading to much larger *F*_{m,h} functions that are tuned based on model performance (e.g., Louis et al. 1982; Beljaars and Holtslag 1991; Viterbo et al. 1999).

So-called 1.5-order closure schemes are more advanced than first-order schemes as they also solve the prognostic turbulent kinetic energy (TKE) equation. This enables *K*_{m,h} to be parameterized on TKE (e.g., Mellor and Yamada 1982; Teixeira and Cheinet 2004, and others) and have been shown to perform well compared to first-order schemes (Cuxart et al. 2006). While popular in research and mesoscale models, these 1.5-order schemes are not typically used in operational global atmospheric models (Cuxart et al. 2006) and have not received as much attention in terms of their high-level scaling behavior (Baas et al. 2008). The parameterization used is often of the form $Km,h=lm,heFm,h$, where *e* is the TKE and a commonly employed mixing length is $lm,h\u221d\tau e$ (Deardorff 1980; Cuxart et al. 2000, 2006; Baas et al. 2008). Here, *τ* is a time scale typically related to the BV frequency in stably stratified flows (Deardorff 1980).

Critically, as will be demonstrated later, this commonly employed mixing-length parameterization does not yield a unique steady-state solution for TKE, when the TKE equation consisting solely of shear production, buoyant destruction and dissipation is considered. This form of the stationary TKE equation is a physically realizable state that has been observed in stably stratified homogeneous sheared turbulence (Gerz et al. 1989; Holt et al. 1992), is assumed in Monin–Obukhov similarity theory and even explicitly exploited to achieve stationary conditions in numerical simulations (Chung and Matheou 2012). Therefore, the inability to gracefully handle this fundamental steady-state behavior represents a troublesome deficiency in these 1.5-order schemes.

A similar, yet distinct, steady-state deficiency has been identified by Baumert and Peters (2000) for the 2.5-order scheme of Mellor and Yamada (1982). This involves, in part, solving prognostic equations for TKE and additionally for a master mixing length *l*. Baumert and Peters (2000) showed that this master mixing-length prognostic equation is essentially the same as the TKE dissipation rate equation, therefore referring to the scheme as a *k*–*ε* closure scheme following engineering nomenclature. At any rate, the form of these two prognostic equations in Mellor and Yamada (1982) yields no steady-state solution for TKE, regardless of mixing-length parameterization. Baumert and Peters (2000) speculated that this issue was a possible reason for ad hoc limiters being applied to the mixing length in future 2.5-order models (e.g., in Galperin et al. 1988). Alternatively, additional terms may be added to the TKE and dissipation equations to represent the transfer of turbulent energy to internal gravity waves, which appears to avoid the issue (Baumert and Peters 2004; Zeng et al. 2020). Rather than adding complexity and computational cost by solving an additional highly parameterized prognostic equation for dissipation, we aim to retain the simplicity of the present 1.5-order closure scheme by only solving for TKE.

The outline of this paper is a follows. The TKE closure model is developed in section 2, with the steady-state deficiency demonstrated and then addressed in section 2b. A brief validation of the model is provided in section 3, although as we are more concerned with the broad properties and steady-state behavior of the model we do not attempt to finely tune the model constants. The scaling behavior of the model is then analyzed in section 4, similar to the analysis by Baas et al. (2008) for their TKE closure scheme, wherein comparisons are made to the first- and 1.5-order closure models in terms of *ϕ*_{m,h} and *F*_{m,h}. Conclusions are then offered in section 5.

## 2. TKE closure model

### a. Model formulation

The prognostic equation for TKE (*e*) assuming horizontally homogeneous conditions can be written as (Stull 1988)

where the terms on the right-hand side correspond to transportation (due to turbulence and pressure diffusion), buoyant production, shear (or mechanical) production, and dissipation, respectively.

The transport term is herein neglected, similar to other SBL modeling studies (e.g., Ellison 1957; Zilitinkevich et al. 2010; Wilson and Venayagamoorthy 2015), as it is often found to be small in the SBL (Nieuwstadt 1984) and typically has negligible impact in models when included in the full TKE prognostic equation (Baas et al. 2008). However, especially during strong intermittent turbulent events, the transport term can become significant (Cuxart et al. 2002) indicating this assumption may restrict us to sustained turbulent conditions such as in weakly and moderately stable cases. From (9), we note that for steady stably stratified turbulence without the transport term, the buoyant destruction cannot exceed the shear production, so that Ri_{f} < 1 and should tend to a constant value for large stability (Monin and Yaglom 1971, their section 7.3; Zilitinkevich et al. 2010).

Following other SBL closure schemes (e.g., Deardorff 1980; Mellor and Yamada 1982; Cuxart et al. 2006; Mauritsen et al. 2007; Baas et al. 2008), the dissipation is parameterized using the Kolmogorov approach, with

where *l*_{ε} is the dissipation length scale and *C*_{ε} a constant. Under neutrally stratified, isotropic, and homogeneous turbulence, *C*_{ε} is often taken to be 0.7; however, these assumptions of isotropy and homogeneity break down with increasing stratification. The value of *C*_{ε} is therefore often reduced as a result (Cuxart et al. 2006); here we use *C*_{ε} = 0.16 from Teixeira and Cheinet (2004), which is similar to the stable-limit value of 0.19 in Deardorff (1980). As with neglecting the TKE transport term above, this dissipation parameterization may not be suitable for very stable and intermittent cases when the turbulence may be anisotropic and inhomogeneous.

Following the eddy-diffusivity approach in (1) and (2), we parameterize the momentum eddy viscosity and heat eddy diffusivity on TKE, with

thus requiring the mixing lengths for momentum and heat, *l*_{m} and *l*_{h}, to be determined. The coefficients *C*_{m} and *C*_{h} are taken to be constants, with *C*_{m} = 0.1 (from the large eddy simulation TKE closure scheme of Deardorff 1980) and *C*_{h} = *C*_{m}/0.75. This will be shown to yield a turbulent Prandtl number in neutrally stratified conditions of 0.75. Finally, under steady-state conditions and with the assumptions made above, (9) then becomes

### b. Mixing-length definitions

For the mixing-length parameterization, a common approach (e.g., Deardorff 1980; Lenderink and Holtslag 2004; Cuxart et al. 2006; Baas et al. 2008) is to put $lm,h,\epsilon =\tau e$, where *τ* is some time scale typically related to the inverse of the BV frequency (Deardorff 1980). However, this is problematic as we see that (13) then results in

which does not yield a unique solution for TKE. This could readily result in atmospheric models predicting the unbounded increase of TKE with time, as was observed in Baas et al. (2008) for particular model constants. We therefore speculate that this issue perhaps contributes to the numerical instabilities often associated with 1.5-order closure schemes (Lenderink et al. 2004). Furthermore, like the similar steady-state issue identified with 2.5-order models in Baumert and Peters (2004), the nonuniqueness may be the reason that additional ad hoc limiters on mixing lengths, or even the TKE, are used to artificially constrain 1.5-order models.

The nonuniqueness of (14) emerges due to the identical mixing-length parameterization employed. To resolve this, we instead use different mixing-length formulations for momentum and heat, a similar approach as done in Teixeira et al. (2004) for unstably stratified flows. Given that the mixing length is a highly conceptual parameter representing the length scale of turbulent mixing, then there is no physical reason that the two must be parameterized in the same manner. We also aim to avoid the use of any ad hoc limiters or stability correction functions on the mixing length or on TKE.

For momentum we simply use the formulation of Blackadar (1962),

which scales as *κz* close to the surface and where *l*_{∞} is the asymptotic turbulent mixing length far from the surface. Typically *l*_{∞} is approximately 40 to 200 m (Cuxart et al. 2006); however, studies using LES (Huang et al. 2013) and field observations (Kim and Mahrt 1992; Sun 2011) have suggested smaller values of approximately 5 to 15 m for stably stratified flows. Here we use *l*_{∞} = 7 m following Huang et al. (2013). Note that while there is no dependence on stability, unlike some other momentum mixing lengths formulations for first-order schemes (Delage 1974; Huang et al. 2013), in 1.5-order closure schemes the eddy viscosity will be indirectly affected by stability due to its dependence on *e* through (11).

To distinguish the parameterization for the mixing length for heat from that of momentum, we can take advantage of the observed increase in the turbulent Prandtl number, *K*_{m}/*K*_{h}, with stability (Ellison 1957; Monin and Yaglom 1971; Kim and Mahrt 1992; Sukoriansky et al. 2006; Venayagamoorthy and Stretch 2010; Huang and Bou-Zeid 2013; Li 2019). This increase is often attributed to momentum being mixed more efficiently than heat due to gravity waves (Lenderink and Holtslag 2004; Anderson 2009). Ideally the gravity waves and turbulence would be parameterized separately, where attempts have been made to do so with higher-order closure schemes (Zilitinkevich 2002) or by adding additional source or sink terms to the TKE and dissipation equations in *k*–*ε* schemes (Baumert and Peters 2004; Zeng et al. 2020). However, this is challenging due to the difficulty in even distinguishing the two motions of turbulence and gravity waves apart from measurement data (Stewart 1969; Jacobitz et al. 2005). Therefore, in the present 1.5-order closure scheme we will simply require the turbulent Prandtl number to become very large for extremely stable situations (Ri_{g} → ∞), to approximate this influence of gravity waves (Lenderink and Holtslag 2004). Note that, by requiring an unbounded increase in Pr with no critical Ri_{g}, this scheme will not obey the stable-limit local-scaling theory of Nieuwstadt (1984), which is based on the *z*-less scaling arguments of (Wyngaard and Coté 1972). Meanwhile, the mixing length for heat under neutrally stratified conditions (Ri_{g} → 0) should converge toward that of momentum.

The two limiting behaviors above can be achieved with

or equivalently

where *τ* = *α*/*N* is a time scale based on the BV frequency and constant *α* = 0.76 (Deardorff 1980; Moeng 1984). We see that, for extremely stable situations as $\tau e/lm\u21920$, the mixing length for heat $lh\u2192\tau e$ and the turbulent Prandtl number *K*_{m}/*K*_{h} → ∞. Meanwhile, under neutrally stratified conditions as $\tau e/lm\u2192\u221e$, we have *l*_{h} → *l*_{m} and *K*_{m}/*K*_{h} → *C*_{m}/*C*_{h} = 0.75, as required.

Finally, we set the dissipation mixing length *l*_{ε} = *μl*_{h} with *C*_{ε}/*μ* = 0.08, a value similar to that used in other studies of the SBL (e.g., Cuxart et al. 2000; Lenderink and Holtslag 2004; Baas et al. 2008). Conventionally the dissipation mixing length is based on that of momentum, albeit in schemes where there is no distinction between *l*_{m} and *l*_{h} such that $l\epsilon \u221d\tau e$. Here, we use the mixing length for heat as it retains this dependency on TKE in the stable limit and is therefore similar to these previous schemes. Using a blend of *l*_{m} and *l*_{h} to define the dissipation mixing length [e.g., Eq. (8) of Teixeira et al. (2004)] does not significantly change the results or conclusions of the present study, presumably because *l*_{h} (and hence *l*_{ε}) already has some dependency on *l*_{m}.

Under the above mixing-length formulation, (13) then becomes

and can now be numerically solved for the steady-state turbulent kinetic energy.

## 3. Model validation with CASES-99

Before analyzing the scaling properties of the model, we provide a short validation that the model agrees reasonably well with field experiments. For this purpose, we use the 1999 Cooperative Atmosphere–Surface Exchange Study (CASES-99; Poulos et al. 2002), a field campaign conducted in southeastern Kansas (37.65°N, 96.74°W; 440 m MSL) in October 1999. Both weakly and very stable conditions were observed, along with other SBL events such as internal gravity waves. A sixth-order polynomial fit to the 60-m main tower wind speed and temperature data is used to determine velocity and potential temperature gradients (Sorbjan and Grachev 2010). The 5-min-averaged measurements are transformed to 1-h averages and bin-averaged for Ri_{g}, as in Wilson and Venayagamoorthy (2015).

Figure 1 shows the mixing lengths for momentum and heat for the CASES-99 data at *z* = 50 m. Here, we assume that the shear production of TKE is equal to *K*_{m}*S*^{2} and use (2) to compute *K*_{m} and *K*_{h}, respectively. The mixing lengths *l*_{m,h} are then obtained assuming the TKE eddy diffusivity parameterization given by (11) and (12). The Ri_{g}-bin averaged mixing lengths from CASES-99 are shown in red, while the blue symbols show the mixing lengths we would obtain using the present formulation given by (15) and (17), where in (17) we use the TKE from CASES-99. We see that *l*_{m} (Fig. 1a) is not particularly sensitive to Ri_{g}, thus justifying the use of a stability independent mixing length for momentum, (15). As mentioned in section 2b, this is because within 1.5-order closure schemes the eddy viscosity remains dependent on stability through *e* in (11), even if *l*_{m} is stability independent. This would not be the case for first-order models that therefore often incorporate the stability correction function *F*_{m}, or mixing lengths defined to be functions of Ri_{g} (e.g., Huang et al. 2013).

The mixing length for heat (Fig. 1b), meanwhile, shows a much stronger dependence on Ri_{g}. This is somewhat captured by the present *l*_{h} formulation [(17), blue line] and could be improved by reducing *α* in the time scale *τ* = *α*/*N* and adjusting other constants accordingly. As we are more interested in the broad properties of the present parameterization, we will continue to use the present model constants that are commonly used in the literature and not attempt to finely tune the constants.

Figure 2a shows the turbulent Prandtl number. The CASES-99 data show a clear increase with Ri_{g}, which is again reasonably captured by the present formulation. These are in agreement with the Pr formulation proposed by Venayagamoorthy and Stretch (2010, gray dashed line). Finally, Fig. 2b shows the TKE from the CASES-99 data along with that obtained from solving the present parameterized steady-state TKE equation of (18). The lapse rate in the model is fixed to the CASES-99 average at *z* = 50 m of *dθ*/*dz* = 40 K km^{−1}, such that the time scale *τ* ≈ 20.1 s, and the shear rate (and thus Ri_{g}) is varied. The resulting TKE predicted by the model shows good agreement with the CASES-99 data and correctly reduces with Ri_{g}.

This analysis can be repeated at different heights where data are available from the CASES-99 main tower and leads to similar results (not shown). Ultimately, the above validation demonstrates that the present TKE model with separately parameterized mixing lengths for momentum and heat captures the essential behavior of the CASES-99 field data. This is despite the model being relatively simplistic, in which only shear production, buoyant destruction and dissipation of the TKE assuming steady-state conditions are considered and standard model constants are employed.

## 4. Scaling behavior

We now look at the scaling behavior of the model in terms of the nondimensional gradients *ϕ*_{m,h} in (4) and (5), and stability functions, *F*_{m,h} = 1/(*ϕ*_{m}*ϕ*_{m,h}), that are often employed in first-order closure schemes. As the parameterized steady-state TKE equation without transport term (18) is independent of *z* we must determine an appropriate vertical length scale. This is achieved by noting that for neutrally stratified surface layer flows *ϕ*_{m}(*ζ* → 0) = 1, which therefore prescribes *z* when *κ* = 0.4 is already specified (Chung and Matheou 2012). This choice of *z* simply guarantees *ϕ*_{m}(*ζ* = 0) = *F*_{m}(Ri_{g} = 0) = 1.

Figure 3 shows the nondimensional velocity gradient *ϕ*_{m} and potential temperature gradient *ϕ*_{h} for the present model. This is accompanied by a selection of first-order closure models often encountered in the literature, introduced in section 1, as well as cases C and D of the 1.5-order TKE closure model of Baas et al. (2008). In Baas et al. (2008), the problematic $lm,h\u221d\tau e$ mixing-length formulation was used where case C used original operational model constants and had a critical Ri_{g} of 1.3, above which turbulence ceased. Case D had no critical Ri_{g} and used model constants that, for their scheme, yielded no solution for TKE in the stable limit.

For momentum (Fig. 3a), the present model appears similar to case C of the TKE closure scheme of Baas et al. (2008). In the limit of large stability, when *ζ* → ∞ and $lh\u2192\tau e$, it can be shown that the present model tends toward *ϕ*_{m}/*ζ* = *B*_{m} where *B*_{m} = 1 + (*C*_{ε}/*μ*)/(*C*_{h}*α*^{2}) ≈ 2.04 is just a function of model constants. From (7), this therefore corresponds to a critical flux Richardson number Ri_{f,crit} = 1/*B*_{m} ≈ 0.49, similar to the critical value of 0.55 for Case C of Baas et al. (2008). Note that, for positive model constants, the present model guarantees that Ri_{f,crit} < 1, as required. Moreover, since the model has a critical Ri_{f,crit} but also *l*_{h} defined such that the turbulent Prandtl number *K*_{m}/*K*_{h} → ∞ in the stable limit, then from (8) there can be no critical gradient Richardson number. The present model limiting behavior for *ϕ*_{m} is less than the model of Businger et al. (1971) and Dyer (1974), which obeys MOST and where in the stable limit *ϕ*_{m}/*ζ* = 4.7. Moreover, it is greater than the high-mixing operational schemes of Beljaars and Holtslag (1991) and Viterbo et al. (1999), where in the stable limit *ϕ*_{m}/*ζ* = 1 and approaches the Ri_{f} > 1 limiting regime shown by the gray triangular area.

The nondimensional temperature gradient *ϕ*_{h} (Fig. 3b), meanwhile, shows a very sharp increase with stability for the present model. This is even more rapid than the nonphysical Case D of Baas et al. (2008); however, the present scheme with separate momentum and heat mixing-length formulations avoids the nonphysical behavior of that particular case. Due to the present model having a critical flux Richardson number but no critical gradient Richardson number, then we see from (6) that in the stable limit as *ζ* → ∞, $\varphi h/\zeta =\u2061(\varphi m/\zeta )2Rig=(1/Rif,crit2)Rig$ and thus increases with Ri_{g}. From the parameterized TKE Eq. (18), the TKE in the stable limit (when $lh\u2192\tau e$) is $estable=Be2\u22c5lm2S2/Rig$ with *B*_{e} = *C*_{m}/(*C*_{h}*B*_{m}*α*) ≈ 0.48. Using (3), (11), (12), we can then obtain Ri_{g} = (*ζ*/*B*_{ζ})^{4}, where *B*_{ζ} = (*C*_{h}*α*/*B*_{m})^{1/2}*κz*/(*C*_{m}*l*_{m}). Ultimately this shows that $\varphi h/\zeta =\zeta 4/\u2061(Rif,crit2B\zeta 4)$ in the stable limit, wherein the large exponent of 4 explains the rapid increase of *ϕ*_{h} with *ζ*.

The stability functions for momentum and heat are shown in Fig. 4. Both of these functions for the present model are much larger than those of the other schemes, which suggests significantly enhanced mixing. However, this interpretation of enhanced mixing comes from the fact that, in first-order models, *F*_{m,h} directly modifies *K*_{m,h} as stability correction functions. For the present closure scheme, *F*_{m,h} is determined diagnostically from (4) and (5) using the TKE derived from solving (18), which only considers shear production, buoyant destruction, and dissipation of TKE. No additional limiters or correction functions are required in specifying the mixing lengths, eddy viscosity or eddy diffusivity. These differences may explain why the steady-state TKE determined by the present model compare reasonably well with observations (Fig. 2b), despite *F*_{m} suggesting substantial mixing.

At large stability, the present model stability functions become $Fm=\u2061(Rif,crit/B\zeta )2Rig\u22121/2\u22480.4331Rig\u22121/2$ and $Fh=(Rif,crit3/B\zeta 2)Rig\u22123/2\u22480.21Rig\u22123/2$. This is the same scaling, although with different model constants, as the large stability limit of Viterbo et al. (1999), where $Fm=0.1Rig\u22121/2$ and $Fh=0.0667Rig\u22123/2$. These functions are based on the Louis–Tiedke–Geleyn (LTG) scheme (Louis et al. 1982), which are ultimately derived from the work of Ellison (1957). As in the current work, Ellison (1957) assumed that the turbulent Prandtl number became unbounded in the stable limit and derived an expression for Pr from the steady-state TKE equation without the transport term (see also Monin and Yaglom 1971, their section 7.4). This suggests that the high-mixing LTG functions are not entirely nonphysical as suggested by some authors (e.g., Baas et al. 2008), but are simply a consequence of assuming that the eddy diffusivity becomes negligible relative to the eddy viscosity in the stable limit.

## 5. Conclusions

TKE closure schemes are a promising avenue for modeling the stably stratified planetary boundary layer; however, they have not been as closely analyzed as the more common first-order approaches. For TKE closure schemes, a commonly employed mixing-length parameterization involves using a time scale and the square root of the TKE, $lm,h,\epsilon =\tau e$ (Deardorff 1980). We show that this common formulation yields no unique solution for TKE when considering the steady-state TKE equation consisting of shear production, buoyant destruction, and dissipation. This deficiency may result in the TKE becoming numerically unstable during steady-state conditions, and may explain the use of artificial, ad hoc limiters to constrain the model.

To avoid this steady-state deficiency, the present TKE closure scheme uses separate mixing lengths for momentum and heat, a similar approach as done in Teixeira et al. (2004) for unstably stratified flows. We emphasize that the present formulation does not actually introduce any new mixing-length parameterizations that have not individually been used in the literature; rather it merges previous parameterizations into a scheme that can faithfully account for steady-state conditions.

We assume that the mixing length for heat tends toward the (stability independent) mixing length for momentum under neutrally stratified conditions, and that the turbulent Prandtl number increases with the gradient Richardson number at large stability. The increase of Pr somewhat approximates the influence of gravity waves, which mix momentum more efficiently than heat (Lenderink and Holtslag 2004; Anderson 2009). These two conditions are satisfied with (17) and results in a scheme with no critical Ri_{g}. The model enables a unique solution for TKE to be obtained from the steady-state TKE prognostic equation without the transport term in (18), as required. The assumptions, notably that the turbulence is sustained, isotropic and homogeneous and neglecting intermittent effects, appear limiting at first. However, comparisons were made with data from the CASES-99 field campaign, which observed both weakly and very stable conditions as well as SBL phenomena such as internal gravity waves (Poulos et al. 2002), and reasonable agreement was found with the model (Figs. 1 and 2).

The present model exhibits some similarities to the TKE closure scheme of Baas et al. (2008), which used the deficient $lm,h=\tau e$ formulation. Figure 3 shows that the present *ϕ*_{m} is similar to their Case C, which was based on model constants using operational values, while the nondimensional potential temperature gradient *ϕ*_{h} increases very rapidly, similar to their Case D. While Case D in Baas et al. (2008) was nonphysical in that there was no unique TKE solution in the stable limit, the present formulation with separately parameterized mixing lengths effectively blends their two cases together and avoids this nonphysical behavior.

Finally, the stability functions *F*_{m,h} of the present scheme are shown to scale in the same manner as the first-order closure Louis–Tiedke–Geleyn (LTG) functions of Louis et al. (1982) [or the revised form in Viterbo et al. (1999)] in the stable limit. This is due to the present scheme and the LTG functions, based on the work of Ellison (1957), both of which assume an unbounded turbulent Prandtl number in the limit of large stability with no critical Ri_{g}. As noted with regards to Case D in Baas et al. (2008), this assumption therefore does not obey the stable limit behavior of the local-scaling theory (Nieuwstadt 1984), which assumes a critical Ri_{g} (*z*-less scaling). However, this suggests that the view that the LTG functions artificially enhance mixing and are nonphysical is instead a consequence of the relaxation of the assumption from local-scaling theory that Ri_{g} must remain finite.

The steady-state results discussed in this paper suggest that a similar mixing-length approach would also produce more realistic results for the stable boundary layer in the context of utilizing a fully prognostic TKE equation to determine the eddy diffusivity and eddy viscosity coefficients.

## Acknowledgments

This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Parts of this research were supported by the NASA MAP Program, the U.S. Department of Energy, Office of Biological and Environmental Research, Earth System Modeling (DE-SC0019242), the ONR Marine Meteorology Program (N000141812432), and by the National Science Foundation Climate Process Team (CPT) Program (AGS 1916619).

## REFERENCES

*Bound.-Layer Meteor.*

*J. Hydrometeor.*

*J. Atmos. Sci.*

*Bound.-Layer Meteor.*

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*Bound.-Layer Meteor.*

*J. Appl. Meteor.*

*J. Geophys. Res.*

*J. Atmos. Sci.*

*J. Fluid Mech.*

*Quart. J. Roy. Meteor. Soc.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*Quart. J. Roy. Meteor. Soc.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*J. Fluid Mech.*

*J. Atmos. Sci.*

*Atmos. Sci. Lett.*

*J. Fluid Mech.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*J. Fluid Mech.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Turbul.*

*Tellus*

*Quart. J. Roy. Meteor. Soc.*

*Quart. J. Roy. Meteor. Soc.*

*Quart. J. Roy. Meteor. Soc.*

*Atmos. Res.*

*Bound.-Layer Meteor.*

*Annu. Rev. Fluid Mech.*

*J. Atmos. Sci.*

*Rev. Geophys. Space Phys.*

*J. Atmos. Sci.*

*Statistical Fluid Mechanics: Mechanics of Turbulence.*Vol. 1. MIT Press, 769 pp

*J. Atmos. Sci.*

*Bull. Amer. Meteor. Soc.*

*Bound.-Layer Meteor.*

*Radio Sci.*

*An Introduction to Boundary Layer Meteorology*. Kluwer Academic Publishers, 670 pp

*Nonlinear Processes Geophys.*

*J. Appl. Meteor. Climatol.*

*Bound.-Layer Meteor.*

*Mon. Wea. Rev.*

*J. Fluid Mech.*

*Quart. J. Roy. Meteor. Soc.*

*J. Atmos. Sci.*

*Quart. J. Roy. Meteor. Soc.*

*k*–

*ε*turbulence model for the stable atmosphere

*J. Atmos. Sci.*

*Quart. J. Roy. Meteor. Soc.*

*Bound.-Layer Meteor.*

Proc. ECMWF Workshop on Boundary Layer Parameterization, Reading, United Kingdom, European Centre for Medium-Range Weather Forecasts, 59–79