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.
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 Km and eddy diffusivity Kh must then be parameterized.
To aid the interpretation of common parameterizations for Km and Kh 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 is the local friction velocity and 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 S2 = (∂u/∂z)2 + (∂υ/∂z)2 and Brunt–Väisälä (BV) frequency N2 = (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 Km,h is parameterized as , where lm,h is the mixing length for momentum or heat. The term Fm,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 Fm,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 Fm,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 Rig > 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 Fm,h, which enables turbulence to remain nonnegligible for Rig → ∞ (e.g., King et al. 2001). Operational models often require even more mixing than these forms of Fm,h provide to avoid the decoupling behavior, leading to much larger Fm,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 Km,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 , where e is the TKE and a commonly employed mixing length is (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 Fm,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 Rif < 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.
thus requiring the mixing lengths for momentum and heat, lm and lh, to be determined. The coefficients Cm and Ch are taken to be constants, with Cm = 0.1 (from the large eddy simulation TKE closure scheme of Deardorff 1980) and Ch = Cm/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 , 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, Km/Kh, 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 (Rig → ∞), to approximate this influence of gravity waves (Lenderink and Holtslag 2004). Note that, by requiring an unbounded increase in Pr with no critical Rig, 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 (Rig → 0) should converge toward that of momentum.
The two limiting behaviors above can be achieved with
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 , the mixing length for heat and the turbulent Prandtl number Km/Kh → ∞. Meanwhile, under neutrally stratified conditions as , we have lh → lm and Km/Kh → Cm/Ch = 0.75, as required.
Finally, we set the dissipation mixing length lε = μlh 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 lm and lh such that . 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 lm and lh 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 lh (and hence lε) already has some dependency on lm.
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 Rig, 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 KmS2 and use (2) to compute Km and Kh, respectively. The mixing lengths lm,h are then obtained assuming the TKE eddy diffusivity parameterization given by (11) and (12). The Rig-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 lm (Fig. 1a) is not particularly sensitive to Rig, 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 lm is stability independent. This would not be the case for first-order models that therefore often incorporate the stability correction function Fm, or mixing lengths defined to be functions of Rig (e.g., Huang et al. 2013).
The mixing length for heat (Fig. 1b), meanwhile, shows a much stronger dependence on Rig. This is somewhat captured by the present lh 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 Rig, 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 Rig) is varied. The resulting TKE predicted by the model shows good agreement with the CASES-99 data and correctly reduces with Rig.
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, Fm,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) = Fm(Rig = 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 mixing-length formulation was used where case C used original operational model constants and had a critical Rig of 1.3, above which turbulence ceased. Case D had no critical Rig 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 , it can be shown that the present model tends toward ϕm/ζ = Bm where Bm = 1 + (Cε/μ)/(Chα2) ≈ 2.04 is just a function of model constants. From (7), this therefore corresponds to a critical flux Richardson number Rif,crit = 1/Bm ≈ 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 Rif,crit < 1, as required. Moreover, since the model has a critical Rif,crit but also lh defined such that the turbulent Prandtl number Km/Kh → ∞ 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 Rif > 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 ζ → ∞, and thus increases with Rig. From the parameterized TKE Eq. (18), the TKE in the stable limit (when ) is with Be = Cm/(ChBmα) ≈ 0.48. Using (3), (11), (12), we can then obtain Rig = (ζ/Bζ)4, where Bζ = (Chα/Bm)1/2κz/(Cmlm). Ultimately this shows that 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, Fm,h directly modifies Km,h as stability correction functions. For the present closure scheme, Fm,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 Fm suggesting substantial mixing.
At large stability, the present model stability functions become and . This is the same scaling, although with different model constants, as the large stability limit of Viterbo et al. (1999), where and . 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.
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, (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 Rig. 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 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 Fm,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 Rig. 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 Rig (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 Rig 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.
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).