## Abstract

Based on a priori analysis of large-eddy simulations (LESs) of the convective atmospheric boundary layer, improved turbulent mixing and dissipation length scales are proposed for a turbulence kinetic energy (TKE)-based planetary boundary layer (PBL) scheme. The turbulent mixing length incorporates surface similarity and TKE constraints in the surface layer, and makes adjustments for lateral entrainment effects in the mixed layer. The dissipation length is constructed based on balanced TKE budgets accounting for shear, buoyancy, and turbulent mixing. A nongradient term is added to the TKE flux to correct for nonlocal turbulent mixing of TKE. The improved length scales are implemented into a PBL scheme, and are tested with idealized single-column convective boundary layer (CBL) cases. Results exhibit robust applicability across a broad CBL stability range, and are in good agreement with LES benchmark simulations. It is then implemented into a community atmospheric model and further evaluated with 3D real-case simulations. Results of the new scheme are of comparable quality to three other well-established PBL schemes. Comparisons between simulated and radiosonde-observed profiles show favorable performance of the new scheme on a clear day.

## 1. Introduction

Mesoscale numerical weather prediction (NWP) and global circulation models rely on planetary boundary layer (PBL) schemes to represent the effects of subgrid-scale (SGS) turbulent mixing on the resolved variables. In accordance with the Reynolds-averaged governing equations, conventional PBL schemes parameterize ensemble-averaged SGS turbulent fluxes. Since the 1970s, PBL schemes have evolved from similarity theory-based bulk models to more sophisticated *K*-profile schemes and higher-order closures (LeMone et al. 2019, section 10c). Today, popular community atmospheric models such as the Weather Research and Forecasting (WRF) Model (Skamarock et al. 2008) includes 12 PBL schemes in its latest v4.0 release.

Among WRF’s 12 PBL schemes, 8 are higher-order schemes while the rest are first-order schemes. The order of a scheme refers to the highest order of moments solved by a prognostic equation (Stull 1988, chapter 6). The popularity of higher-order closures stems from the seminal work of Mellor and Yamada (1974). It is generally believed that the uncertainties associated with turbulence closure assumptions can be reduced when constrained by prognostic equations of higher-order moments. On the downside, higher-order closures not only require solving more prognostic equations, but also involve more complicated closure assumptions and a drastic increase in the number of model coefficients to be determined (Holt and Raman 1988; LeMone et al. 2019). Given the pros and cons of higher-order closures, it is not surprising that the 1.5-order turbulence kinetic energy (TKE-1.5)-based schemes have been the most popular among higher-order PBL schemes. They are classified as 1.5 order because the prognostic TKE equation is the only second-order equation retained.

Most TKE-1.5 schemes are constructed upon the gradient-diffusion framework, also known as *K* theory,

where the overline represents the Reynolds-average operator, upper- and lowercase letters represent Reynolds-averaged and turbulent variables, respectively. The vertical turbulent flux $w\varphi \xaf$ is parameterized as a product of the local gradient of resolved variable Φ and an eddy diffusion coefficient *K* (m^{2} s^{−1}). In 1.5-order schemes, *K* is further parameterized as

where $e\xaf\u2261(1/2)u2+\upsilon 2+w2\xaf$ (m^{2} s^{−2}) is the TKE, and *L*_{k} (m) is a turbulent mixing length that must be parameterized or “closed.” The other key length scale of a TKE-1.5 scheme is the dissipation length *L*_{ε} (m), which is used to parameterize viscous dissipation of TKE,

where *ε* (m^{2} s^{−3}) is the dissipation rate. Equation (3) is based on the Richardson–Kolmogorov cascade, where *L*_{ε} is on the order of the integral scale of turbulence (Vassilicos 2015). Note the absence of model coefficients in Eqs. (2) and (3), as they are absorbed into the definitions of *L*_{k} and *L*_{ε}.

In this study, turbulent mixing and dissipation lengths are diagnosed a priori from high-resolution large-eddy simulations (LESs) of a set of convective boundary layer (CBL) cases with different bulk stability as described in section 2. LES is based on the definition of a spatial filter where the large energy-containing eddies are explicitly resolved, while the effects of the smaller eddies on the resolved flow are modeled by a turbulence closure. LES-modeled atmospheric boundary layers (ABL) have been used extensively to construct and improve PBL schemes (see LeMone et al. 2019, section 2c, and references therein).

Based on LES-diagnosed benchmark profiles of *L*_{k} and *L*_{ε}, we evaluate their formulations in some of the popular higher-order schemes, and expose their general limitations in section 3. Improved models of *L*_{k} and *L*_{ε} are then constructed to include more physical insights as well as constraints. Note that this study focuses on the mixing length of heat, while momentum mixing is left for future work. Furthermore, a nongradient TKE flux is added to improve the representation of the turbulent mixing term in the prognostic TKE equation. All improvements are then integrated into an existing 1.5-order TKE-based nonlocal PBL scheme, and implemented in the WRF Model. It is tested in both idealized and real CBL cases in section 4. Its performance is evaluated with benchmark LES and observations, and cross-compared with other well-established PBL schemes.

## 2. Case description and numerical setup

Five idealized CBL cases are simulated and a neutral boundary layer (NBL) case from a previous study (Zhou et al. 2019) is adopted for this work. A list of model parameters is presented in Table 1. The CBL cases have identical model setup as those in Zhou et al. (2019), but with 4 times the horizontal domain size to include more samples (i.e., the organized convective structures) for averaging. Briefly, the simulations are driven by constant surface heat fluxes of 0.2, 0.05, and 0 K m s^{−1} for the *B*, *S*, and *N* cases, respectively. Wind shear across the boundary layer is varied by prescribed geostrophic winds ranging from 1 to 15 m s^{−1}. The five CBL cases cover a broad stability spectrum ranging from case B1 approaching the free convective limit to the shear-dominated case S15. CBL stability is quantified by the dimensionless parameter −*z*_{i}/*L*, which is interpreted as a bulk measure of the relative importance of thermal to mechanical forcings (Lenschow and Stephens 1980), where *z*_{i} (m) is the boundary layer depth, $L\u2261\u2212u*3\Theta 0/\u2061(\kappa gw\theta \xafs)$ (m) is the Obukhov length, $u*$ (m s^{−1}) is the surface friction velocity, *θ* (K) is the potential temperature, Θ_{0} (K) is a reference potential temperature, $w\theta \xafs$ (K m s^{−1}) is the surface heat flux, and *κ* and *g* are the von Kármán and the gravitational acceleration constants. The initial conditions follow Shin and Hong (2013), where an isentropic initial temperature profile is capped by a sharp inversion between 925 and 1075 m, which effectively limits the growth of the boundary layer. The initial wind profile is geostrophic. The Coriolis parameter *f* is 10^{−4} s^{−1}, and the surface roughness length *z*_{0} is 0.1 m. The *B*, *S*, and *N* cases reach quasi-steady state after about 2, 5, and 8 h, and are run for another 2 h with output saved at every 10 min for analysis.

The five CBL cases are performed using (1008, 1008, 328) grid points in the (*x*, *y*, *z*) directions with 10 m horizontal grid spacing. A uniform 4 m vertical grid spacing is adopted below 1.3 km to ensure high resolution for the CBL and is progressively stretched upward to the domain top at 2 km. The NBL is simulated on a (9.216, 9.216, 1.5) km domain with 8 m isotropic grid spacing. All simulations are performed on horizontally periodic domains with Rayleigh damping applied to the top one-quarter of the domain. The 1.5-order TKE closure of Moeng (1984) is used for SGS turbulence. The Advanced Regional Prediction System (ARPS) developed at the Center for Analysis and Prediction of Storms at the University of Oklahoma, is used for the simulations. ARPS is a nonhydrostatic mesoscale and small-scale finite-difference NWP model. It uses a generalized height-based terrain-following coordinate on an Arakawa *C* grid. A mode-splitting time integration scheme is employed (Klemp and Wilhelmson 1978). The LES runs adopt fourth-order advection for the horizontal and vertical directions, with a monotonic flux limiter applied to scalar advection (Xue 2000). More details on ARPS are documented in Xue et al. (2000, 2001).

## 3. Length-scale formulations

For 1.5-order PBL schemes, the eddy diffusion coefficient is formulated according to Eq. (2), in which $e\xaf$ is obtained by solving the Reynolds-averaged prognostic TKE equation,

The Boussinesq approximation for the boundary layer is used here. As PBL schemes assume horizontal homogeneity of SGS turbulence statistics, only terms associated with vertical fluxes are retained. The first term on the right-hand side (rhs) represents buoyancy production, the second and third terms together represent shear production, the fourth is the vertical turbulent mixing term, and the fifth is a pressure correlation term. The turbulence and pressure terms redistribute TKE throughout the depth of the atmospheric column. Therefore, $we\xaf$ and $wp\xaf/\rho 0$ are usually grouped together and parameterized as a turbulent diffusion term (Stull 1988, chapter 6.5),

where *K*_{e} is the eddy diffusivity of TKE. The last term *ε* represents viscous dissipation and is parameterized according to Eq. (3) based on the definition of a dissipation length *L*_{ε}.

### a. Turbulent mixing length

In PBL schemes, turbulent fluxes of momentum, heat, and other scalars are conventionally parameterized based on Eq. (1) with different values of *K*. The eddy viscosity for momentum *K*_{M} and the eddy diffusivity for heat *K*_{H} are linked by a turbulent Prandtl number Pr ≡ *K*_{M}/*K*_{H}, which is approximately 1 for ABL flows. Eddy diffusivities of other scalars including moisture and trace gases are often set equal to *K*_{H} (Li 2019). In this study, we focus on the turbulent mixing of heat, and leave momentum to a separate study. For the rest of the article, the turbulent mixing length *L*_{k} strictly applies to potential temperature.

Parameterization of the vertical heat flux $w\theta \xaf$ for the CBL often requires an additional countergradient correction term *γ* (K m^{−1}) to account for the effects of nonlocal organized convection (see Zhou et al. 2018, and references therein),

Inclusion of *K*_{H}*γ* in Eq. (6) effectively counteracts the downward (negative) gradient-diffusion flux −*K*_{H}∂Θ/∂*z* in the upper part of the CBL, where the temperature stratification is slightly stable (see Hu et al. 2019, and references therein). In this study, Holtslag and Moeng’s (1991) model is adopted,

where $w*\u2261\u2061(g/\Theta 0w\theta \xafszi)1/3$ (m s^{−1}) is the convective velocity scale. Based on Eqs. (2) and (6),

Substituting in the expression of *γ* from Eq. (7), *L*_{k} can then be derived a priori from LES data. In the analysis, Reynolds averaging is replaced by horizontal averaging. The two operators are effectively equivalent when the dominant length scale of the CBL, that is the boundary layer depth (~1 km in the current setup), is much smaller than the domain width (~11 km) (Wyngaard 2004).

The vertical profiles of *L*_{k}/*z*_{i} are presented in Fig. 1a. The parameter *z*_{i} is diagnosed as the height of minimum $w\theta \xaf$. *L*_{k}/*z*_{i} exhibits distinct patterns above and below ~0.4*z*_{i}. In the lower half of the CBL, apparent stability dependence is observed as *L*_{k}/*z*_{i} decreases monotonically with decreasing bulk CBL stability measured by −*z*_{i}/*L*. Above ~0.4*z*_{i}, *L*_{k}/*z*_{i} for the five CBL cases collapse into a single curve that decreases almost linearly with height. In the stably stratified free troposphere above the CBL, *L*_{k} becomes vanishingly small.

The general profiles of *L*_{k} suggest a three-layer treatment that comprises an extended surface layer where *L*_{k} varies with the surface scale *z*/*L*, an upper mixed layer with a universal mixing length that scales with *z*_{i}, and the free-troposphere whose length scale is primarily limited by stable stratification. The surface layer mixing length $Lks$ is based on the work of Redelsperger et al. (2001), and is specifically designed for TKE-based schemes. The value of $Lks$ is determined by Monin–Obukhov (MO) surface-similarity theory, while satisfying the TKE balance equation simultaneously. Compared to purely MO-based length scales, the TKE-constrained mixing length improves the representation of turbulent momentum and heat fluxes for most stability conditions (Drobinski et al. 2006). $Lks$ is defined as

where *A* = 2.79 is a constant, and *ϕ*_{L} is a MO stability function that also depends on the dimensionless shear $\varphi m\u2261(\kappa z/u*)\u2202U/\u2202z$. For unstable conditions (i.e., *z*/*L* < 0),

The coefficient *α* in Eq. (10) represents part of the surface layer TKE parameterized by $u*2$ in the form of $e\xaf=\alpha u*2+\beta w*2$ (Wyngaard and Cote 1974; Panofsky et al. 1977; Garratt 1992). A similar expression $e\xaf=f\u2061(u*,\u2009w*)$ is evaluated with LES data in section 3b, and is exploited to constrain the parameterized dissipation rate.

The profiles of $Lks$ are presented in Fig. 1b, using *α* = 12 in Eq. (10) for the best fit. The suggested value of 4.63 by Redelsperger et al. (2001) overestimates *L*_{k} for the strongly convective cases B1 and B5, and has a minor effect on the rest of the cases (results not shown). This is likely attributed to the absence of the turbulent mixing and pressure correlation terms in Redelsperger et al.’s TKE balance equation (see TKE budgets presented in Fig. 2). The parameterized $Lks$ agree very well with the LES-diagnosed profiles in the surface layer for all CBL cases considered.

In the upper CBL where *L*_{k}/*z*_{i} exhibits universality, the Bougeault–Lacarrere (BouLac; Bougeault and Lacarrere 1989) length scales are adopted as the base scales,

where *l*_{u} and *l*_{d} represent the upward and downward distances that an undiluted air parcel at height *z* with a mean TKE $e\xaf\u2061(z)$ can travel against stratification. They can be derived based on a balance between buoyancy and vertical advection of *w*. In Eq. (13)*l*_{d} is further limited by the distance to ground. The BouLac mixing length $LkBL$ and dissipation length $L\epsilon BL$ are then formulated as

where $CkBL=0.4$, $C\epsilon BL=1.4$ are model constants. Despite its simplicity, the BouLac scheme works quite well for the daytime CBL (Bélair et al. 1999; Cuxart et al. 2000; Shin and Hong 2011; Xie et al. 2012; Huang et al. 2018), and is adopted by the convection-permitting operational NWP Application of Research to Operations at Mesoscale (AROME) model (Seity et al. 2011), as well as the French research-oriented mesoscale model MesoNH (Lac et al. 2018).

As noted by Lenderink and Holtslag (2004), *l*_{u} and *l*_{d} are insensitive to the bulk stability of the CBL. This is clearly seen in Fig. 1c, where *l*_{u} and *l*_{d} for all five CBL cases collapse such that they are well approximated by 1 − *z*/*z*_{i} and *z*/*z*_{i}, respectively. This is because Θ is well mixed inside the CBL due to both bottom-up surface heating and top-down entrainment. For an undiluted parcel, *l*_{u} and *l*_{d} are essentially limited by the CBL capping inversion and the ground. The fact that the diagnosed *L*_{k} in Fig. 1a is also insensitive to bulk CBL stability in the upper half of the boundary layer, makes the BouLac scale a potential candidate to construct *L*_{k}. Here, we only make use of *l*_{u} in the upper CBL, as implied by taking the lower value between *l*_{u} and *l*_{d} in the original BouLac scheme according to Eq. (14).

Clearly, some modifications must be made before adopting *l*_{u} as it consistently overpredicts *L*_{k} in Fig. 1c. This is most likely because Eqs. (12) and (13) assume undiluted parcel ascent and descent and do not account for turbulent entrainment and detrainment processes. In the CBL, rising parcels, as part of the organized thermals, mix with the environment and lose part of their buoyancy along the way. The entrainment processes are represented in eddy-diffusivity mass-flux (EDMF) type of PBL schemes with a bulk plume model (Siebesma et al. 2007). Here, a simple approach is adopted to account for entrainment effects by reducing *l*_{u} through an empirically fitted piecewise linear function,

where $Lkm$ is the proposed turbulent length scale for the mixed layer. Profiles of $Lkm$ are presented as dashed–dotted lines in Fig. 1c, and agree well with the diagnosed *L*_{k} above ~0.4*z*_{i}.

The surface and mixed-layer length scales are combined into a single mixing length $Lksm$ through an averaging operator where

The coefficient *a* is set to −3 based on empirical fitting to the LES-diagnosed *L*_{k}.

In the entrainment zone and free troposphere, approximately 0.8*z*_{i} and above, the turbulent mixing length is mainly limited by stratification. Therefore, the BouLac scale is again adopted to set the entrainment mixing length $Lke=CkBLlu$. Combining $Lksm$ and $Lke$, we arrive at the final form of *L*_{k},

Figure 1d presents the parameterized and the diagnosed profiles of *L*_{k}. In general, Eq. (18) is able to reproduce the diagnosed *L*_{k} in the surface and the upper layers. Some discrepancies are found in the transition region between $Lks$ and $Lkm$, especially with case B5. A more sophisticated averaging operator than Eq. (17) might be able to improve the agreement, but is not further attempted.

Last, note that the parameterized *L*_{k} proposed in the present study is by no means a universal length scale representing turbulent mixing. It is formulated based on the specific parameterization of turbulent heat flux by Eq. (6) as well as the particular form of *γ* in Eq. (7). Adopting a mass flux approach instead of countergradient correction to represent nonlocal heat flux (e.g., Siebesma et al. 2007) for example, or using an alternative formulation of *γ* would all change *L*_{k}. Preliminary tests on the structural and parametric sensitivities of *L*_{k} were performed by replacing Eq. (7) with Troen and Mahrt’s (1986) *γ* formulation (structural sensitivity), and by perturbing the model constant in Eq. (7) (parametric sensitivity). It is found that, although the magnitudes of *L*_{k} do change, the characteristic features of the profiles remain similar to those presented in Fig. 1a (results not shown).

### b. Dissipation length

Compared to the mixing length, the dissipation length *L*_{ε} is associated with more uncertainties and has less available parameterization options. This is primarily due to the difficulties in measuring viscous dissipation rate *ε*, which occurs at the Kolmogorov scale, about 1 mm for the ABL (Kaimal and Finnigan 1994, chapter 2.1). As a result, high- and very high-frequency instruments are required for indirect and direct measurements, respectively. Available field data of *ε* are limited to the surface layer (Piper and Lundquist 2004) and the lower mixed layer (Muñoz-Esparza et al. 2018). Because the dissipation scales are not resolvable even by LES, indirect measures are taken to compute *ε*. One approach is to use the parameterized *ε* by the LES turbulence closure. This is supported by the LES’s capability to resolve part of the inertial subrange, and therefore producing faithful predictions of *ε*. This approach has been commonly used to obtain *ε* from LES and to formulate models of *L*_{ε} accordingly (see, e.g., Moeng and Sullivan 1994, hereafter MS94; Nakanishi 2001; Ito et al. 2015). This study adopts the residual method following Heinze et al. (2015), where *ε* is computed as the residual of the balanced TKE budget in order to reduce uncertainties associated with the numerical schemes. The two methods are almost equivalent as shown by the small budget residual terms (sum of all TKE budget terms including LES parameterized *ε*) presented in Fig. 2. Taking *ε* directly from Moeng’s LES closure to diagnose *L*_{ε} only leads to small changes close to the surface and the top of the CBL for the shear-driven cases (see Figs. 2d,e), and have little effects on the analysis below.

For 1.5-order PBL schemes, a widely used parameterization of *L*_{ε} is proposed by Mellor and Yamada (1974),

where *L*_{ε} represents the height of the “center of TKE” in the vertical grid column. In the original writings of Mellor and Yamada, the idea behind Eq. (19) is to “adopt the simplest-but hopefully general-length scale, characteristic of the extent of the turbulence field.” Equation (19) forms the basis of many Mellor–Yamada (MY)-type PBL schemes such as the MYJ scheme (Janjić 2002), the improved MY scheme of Cheng et al. (2002), and the MYNN schemes (Nakanishi and Niino 2009). Other options include the BouLac scale in Eq. (15) and the Blackadar-type scale in the GBM (Grenier and Bretherton 2001) and UW schemes (Bretherton and Park 2009) where

with *λ*_{0} set proportional to *z*_{i}.

Diagnosed *L*_{ε} based on Eq. (3) for the CBL and NBL cases are presented in Fig. 3a. Normalized *L*_{ε} for all CBL cases exhibits a general three-layer structure that corresponds to the sublayers of the CBL. In the surface layer *L*_{ε} increases rapidly with height. Upon transition to the mixed layer, the vertical gradient of *L*_{ε} becomes much smaller in the presence of vigorous turbulent mixing. In the limiting free convective case of B1, *L*_{ε} is nearly constant with height between roughly 0.3 and 0.9*z*_{i}. At the bottom of the entrainment zone around 0.8*z*_{i}, *L*_{ε} obtains its maximum value and subsequently decreases with height approaching the CBL top. Above the CBL, the current LES resolution is not able to provide reliable information about dissipation in the stratified free troposphere. As a result, the diagnosed *L*_{ε} exhibits large oscillations and is not plotted.

Among the five CBL cases, *L*_{ε}/*z*_{i} decreases in the surface layer and increases in the mixed layer as the bulk CBL stability decreases (i.e., −*z*_{i}/*L* → 0). A similar trend is also observed in Fig. 23 of MS94 for their free and sheared convection cases. In the entrainment layer, *L*_{ε} for different CBL cases converge. This is most likely a result of using the same initial capping inversion profile in all cases. For the neutral case N10, *L*_{ε}/*z*_{i} is smaller than that of B1 at all elevations. The abrupt drop of *L*_{ε} is opposite to what one might have expected if the observed trend of *L*_{ε} as a function of *z*_{i}/*L* for the CBL were extrapolated to neutral conditions at *z*_{i}/*L* = 0. This nonlinear behavior toward neutrality is also documented by Fig. 23 of MS94 for their neutral case. We believe that this is due to the use of a loosely defined TKE in Eq. (3). Strictly, TKE in Eq. (3) is an integral of the energy spectrum over the inertial subrange up to the integral scale that is *L*_{ε}. However, in a PBL scheme, $e\xaf$ used in Eq. (3) to compute *ε* represents TKE of all convective motions that also include those associated with organized CBL structures. The size of these coherent structures is of *O*(*z*_{i}), and is of the same order as the CBL integral scale. Zhou et al. (2019) suggested that TKE associated with the organized CBL motions should not be counted as part of the inertial subrange TKE to compute *ε*. When removed from the total TKE, the computed *L*_{ε}/*z*_{i} exhibit a much more universal profile (results to be presented in a forthcoming study).

Figure 3b presents the modeled *L*_{ε} by the MYNN and BouLac [Eq. (15)] schemes. $L\epsilon M$ of the MYNN scheme is an improved version of the original MY formulation. Based on the three-layer structure of *L*_{ε} observed in Fig. 3a, $L\epsilon M$ combines a MO-based scale for the surface layer, a boundary layer scale [Eq. (19)] for mixed layer, and a modified buoyancy length scale for the stably stratified entrainment layer and free troposphere [see Eqs. (52)–(55) in Nakanishi and Niino 2009]. Unlike the benchmark *L*_{ε} in Fig. 3a, $L\epsilon M/zi$ show little variations with bulk CBL stability in the upper CBL. In fact, $L\epsilon M/zi$ exhibits the *opposite* stability dependence and decreases slightly from free convective to neutral conditions in the upper mixed layer. This is expected from Eq. (19) given the variability of $e\xaf\u2061(z)$ with −*z*_{i}/*L*. Figure 4 presents the averaged $e\xaf\u2061(z)$ of all the simulated cases in the descending order of −*z*_{i}/*L*. With reduced buoyancy and increased shear production, which is most significant in the surface layer, $e\xaf\u2061(z)$ becomes more surface heavy as −*z*_{i}/*L* decreases. As such, Eq. (19)–modeled *L*_{ε} must also decrease monotonically from convective to neutral conditions. For the BouLac scheme [Eq. (15)], the parameterized $L\epsilon BL$ are essentially invariant with −*z*_{i}/*L* as a result of its formulation explained in section 3a. Similarly, the Blackadar scale [Eq. (20)] is also insensitive to the bulk CBL stability (results not shown).

The limited and even opposite variability of the modeled *L*_{ε}/*z*_{i} in Fig. 3b is likely attributed to the model’s overreliance on scaling arguments without accounting for the underlying physical processes that constrain dissipation. MS94 took a different pathway and proposed a process-based model of *ε* with the balanced TKE equation. Under steady-state conditions, the prognostic TKE equation [Eq. (4)] is rearranged, and *ε* is diagnosed as

The key idea of MS94 is to model *ε* with the parameterized rhs forcings in Eq. (21). Buoyancy, turbulence and pressure terms are first considered as they all scale with $w*$ and *z*_{i} in the CBL (Zhou et al. 2019).

Figure 5 presents the normalized profiles of these terms for all CBL cases. The turbulent mixing and pressure transport terms are grouped together as explained in the beginning of section 3. The buoyancy production terms, or effectively $w\theta \xaf$, exhibit the classic linear profiles and are well approximated by

where $RH\u2261\u2212w\theta \xafe/w\theta \xafs$ is the ratio of the maximum (i.e., most negative) entrainment flux to the surface flux of heat, and is set to 0.2 following Holtslag and Moeng (1991). Variation of *R*_{H} with respect to bulk CBL stability (Conzemius and Fedorovich 2006; Liu et al. 2016) is not considered here.

The combined turbulence and pressure terms also collapse reasonably well among all CBL cases except for S15. A plausible reason for such agreement is attributed to the dynamics of the organized motion of the CBL as proposed by Zhou et al. (2019). As the amount of shear further increases in case S15, rightward deviation of the profile become larger, suggesting less TKE from the lower half of the CBL being transported upward. A similar trend is also observed in the direct numerical simulations of sheared CBL in Haghshenas and Mellado (2019; see their Fig. 7a). This is in accordance with the general understanding that shear-generated TKE tends to be dissipated locally. In the neutral limit, turbulence and pressure transport almost vanish as shown in Fig. 2f. Case S15 aside, the profile of the combined turbulence and pressure terms is approximated by a linear function,

Equation (23) is similar to the approximation $w*3/zi\u2061(1.2z/zi\u22120.6)$ by MS94, with slightly different coefficients.

where a step function like *f*(*z*/*z*_{i}) = 1/2{1 − tanh[(*z*/*z*_{i} − 1)/0.025]} is factored in only to reduce the rhs of Eq. (24) to zero above the CBL. The modeled forcing compares well with diagnosed forcings as shown in Fig. 5b, except for case S15. Again, the disagreement is due to the modeled transport terms by Eq. (23). Equation (24) is still adopted for all CBL cases, because as shown in the a posteriori tests in section 4, the relatively poor representation of the turbulent and pressure terms for the strongly sheared case does not affect much the modeled dissipation rates due to the dominance of shear and buoyancy forcings under such conditions. In other words, when Eq. (23) [hence Eq. (24)] starts to fail, the very forcing it represents also becomes less important in the TKE budget.

Vertical profiles of shear production rates are presented in Fig. 6. Case B1 is omitted because under 1 m s^{−1} geostrophic wind, the turbulent momentum flux is small, and shows large wiggles that would require an even larger model domain to obtain a representative profile. Shear production rates are normalized with the surface scale $u*3/\kappa z$ to illustrate the spread among different CBL cases in the surface layer where vertical wind shear is the strongest. With this normalization factor, shear production rates at the top of the CBL are exaggerated, whereas in the original profiles, shear production at the CBL top is of a smaller magnitude than that at the surface as shown in Fig. 2. In the mixed layer in-between, shear production rates are close to zero for all CBL cases. In comparison, the normalized shear production rate for the NBL decreases less rapidly than for the CBL cases in the absence of vigorous convective mixing, and is much greater than that of the CBL in the mixed layer.

Compared to buoyancy-controlled forcings, shear production is more difficult to parameterize, MS94 proposed the following expression,

based on the fact that the momentum flux profile $\tau \u2261uw\xaf2+\upsilon w\xaf2$ in the CBL is well represented by $u*2\u2061(1\u2212z/zi)$ (Zhou et al. 2019), and by adopting a simple approximation of the wind shear by extending surface layer similarity upward (i.e., $\u2202U/\u2202z=u*/\kappa z\varphi m$). Modeled shear production terms by Eq. (25) are plotted as dashed lines in Fig. 6. For the CBL, some agreement is found in the surface layer, where the modeled shear production is mainly controlled by the MO stability function *ϕ*_{m}. Higher up, Eq. (25) overpredicts shear production in the mixed layer, and misses completely the elevated peak at the CBL top. Unfortunately, we do not have a closure for the elevated shear production. Instead, its effect is represented indirectly through a shear length scale as discussed next. In comparison, the agreement between the modeled and diagnosed shear production is quite satisfactory for case N10, with only some slight underestimation.

Substituting Eqs. (24) and (25) into Eq. (21), and using Eq. (3) to parameterize *ε*, after rearrangement we arrive at an expression of *L*_{ε} similar to that of MS94,

where

are the dissipation scales of shear-and buoyancy-related processes.

Equation (26) is incomplete without a parameterization for the stably stratified entrainment zone and the free troposphere. For that purpose, the buoyancy length scale $lb\u2261e\xaf/N$, where *N* (s^{−1}) is the buoyancy frequency, is often adopted by 1.5-order PBL schemes. While *l*_{b} may suffice for stratification-dominated conditions, it is insufficient in the entrainment zone when the effect of shear is significant. In the presence of strong wind shear, turbulence and pressure transport terms in the entrainment zone are small, the dominant TKE balance is between shear production, buoyancy destruction, and dissipation (see Fig. 2). Under such circumstances, the increased shear production of TKE in the entrainment zone tends to be dissipated locally, as shown by reduction of the turbulence and pressure transport terms in the entrainment zone in Fig. 2 (see also Fig. 7a in Haghshenas and Mellado 2019).

In a turbulent boundary layer under neutral conditions, Hunt et al. (1987) suggested that *ε* is controlled by the deformation of small eddies by the large energy-containing ones (see also Hunt and Morrison 2000; Tardu 2017). Therefore, they proposed a shear-dissipation scale $w2\xaf/\u2061(dU/dz)$ away from the wall. In the recent study of CBL entrainment, Haghshenas and Mellado (2019) introduced Δ*U*, the velocity difference across the entrainment zone, into their scaling laws for the lower entrainment sublayer to represent shear effects. Therefore, we propose a mixed entrainment dissipation scale $L\epsilon e$ based on the buoyancy *l*_{b} and shear $ls\u2261e\xaf/S$ length scales,

where $S=\u2061(\u2202U/\u2202z)2+\u2061(\u2202V/\u2202z)2$ (s^{−1}) is the mean shear, *c*_{s} and *c*_{b} are constants. Equation (29) is constructed such that $L\epsilon e$ depends on the smaller of the two length scales, which translates to the relative importance of the shear-and buoyancy-related processes that determine *ε*.

The relative dominance between *l*_{b} and *l*_{s} in Eq. (29) is also affected by the ratio of the model constants *r*_{bs} ≡ *c*_{b}/*c*_{s}. *l*_{s} dominates only when *c*_{s}*l*_{s} < *c*_{b}*l*_{b}, or *l*_{s}/*l*_{b} < *r*_{bs}. Given the definitions of *l*_{s} and *l*_{b}, their ratio $ls/lb=Rig$ is the square root of the gradient Richardson number Ri_{g} ≡ *N*^{2}/*S*^{2}. Therefore, *r*_{bs} can be interpreted as a threshold Ri_{g} that determines the transition from buoyancy-dominated conditions where turbulence and pressure transport are the primary TKE sources, to shear-influenced conditions where shear production becomes a significant contributor to the energetics of the entrainment zone. This transition is found to occur when the flux Richardson number Ri_{f} or Ri_{g}, usually evaluated at the height of the minimum buoyancy flux, falls below 1 (Conzemius and Fedorovich 2006; Pino and De Arellano 2008; Liu et al. 2016; Liu et al. 2018; Haghshenas and Mellado 2019). Therefore, *r*_{bs} is set to 1. The coefficient *c*_{b} = 40 is then set to match the entrainment zone dissipation scale of the MYNN model under free-convective conditions. In the a posteriori tests, however, *r*_{bs} is used as a tuning parameter and an optimal value of 4 is found. This is likely associated with the underprediction of shear production in the entrainment zone. Better representation of momentum fluxes might be able to restore *r*_{bs} to 1, but will be left to future work.

With the new entrainment dissipation scale, the complete model of *L*_{ε} is formulated as

where $L\epsilon e$ is based on Eq. (29) under stable stratification (i.e., ∂Θ/∂*z* > 0), and is set to ∞ otherwise. The modeled *L*_{ε} by Eq. (30), which we refer to as the MS model, is compared to the LES-diagnosed *L*_{ε} in Fig. 7. The MYNN profiles from Fig. 3b are also included for comparison purposes. Unlike the scaling-based MYNN model, the process-based MS model is able to reproduce the observed variability of *L*_{ε} with the stability parameter −*z*_{i}/*L*. As a result, better agreement with the LES benchmark is obtained with the MS model than the MYNN model for all CBL cases, especially the free convective and strongly sheared cases. Given its limited variability, the MYNN model only works well for cases B5 and B10. Also notice that in the surface layer, the MS model performs better than the MYNN model for the three buoyancy cases, and are of comparable quality to the MYNN model for the two shear cases. For the NBL, both the MS and the MYNN models compare well with the LES profile. Figure 8 further illustrates the relative contributions from shear, buoyancy, and entrainment scales to *L*_{ε}. It is no surprise that the shear and buoyancy scales are dominant in the surface and mixed layers, respectively. In this regard, the processes-based MS model can also be viewed as a three-layer model as the MYNN. The improvements are mainly from the additional consideration of physical processes that contribute to the TKE budget in each sublayer of the CBL.

In the a posteriori tests, $L\epsilon s$ and $L\epsilon b$ defined by Eqs. (27) and (28) can lead to overdissipation of TKE during the initial model spinup phase when the TKE is small. Therefore, an additional constraint is imposed to regulate the parameterized dissipation rate,

where

is a parameterized surface-layer TKE $\u2061(es\xaf)$ based on the friction and the convective velocities, and $e1\xaf$ is the TKE at the first vertical grid level. Equation (32) was first used by Wyngaard and Cote (1974) as surface boundary condition for $e\xaf$, where *c*_{u} = 2.3 and *c*_{w} = 0.59 are model constants.^{1}Troen and Mahrt (1986) also adopted Eq. (32) for their PBL scheme using *c*_{u} ≈ 1.0 and *c*_{w} ≈ 0.65. In this study, *c*_{u} and *c*_{w} are set to 1.5 and 0.67 respectively based on a linear regression of LES diagnosed $e\xaf3/2$ against $u*3$ and $w*3$ for all five CBL cases (results not shown). In Eq. (31), the ratio *R* effectively enforces Eq. (32) as a surface boundary condition on $e\xaf$, and therefore constrains (reduces) *ε* during the model spinup phase. *R* drops quickly to unity at the end of the spinup period when the model TKE picks up, and remains nearly constant afterward.

### c. Nongradient mixing of TKE

With LES data, the nongradient mixing of TKE is also examined. The concept is similar to the countergradient correction of $w\theta \xaf$ in Eq. (6) but for turbulent flux of TKE $\u2061(we\xaf)$. Based on the budget analysis of $we\xaf$, Therry and Lacarrère (1983) suggested the inclusion of a correction term *γ*_{e} in Eq. (5) to account for the effect of nonlocal mixing of TKE in the CBL,

where *γ*_{e} (m s^{−2}) is contributed mostly by buoyancy effects. The expression of *γ*_{e} involves further closure assumptions of third-order moments. It is somewhat complicated and not shown here. Shin et al. (2013) used Eq. (33) to diagnose TKE from the YSU scheme. In this study, we do not investigate the physical mechanisms behind nongradient TKE fluxes, but assume Eq. (33) with *K*_{e} = *K*_{H}, i.e., adopting the same mixing length for both heat and TKE as is often done for 1.5-order schemes (Bougeault and Lacarrere 1989; Teixeira and Cheinet 2004), and diagnose *γ*_{e} accordingly.

Figure 9a presents the LES-diagnosed $we\xaf$ and $wp/\rho 0\xaf$. Up to the entrainment zone, $we\xaf$ and $wp/\rho 0\xaf$ are of opposite sign. $we\xaf$ is directed upward for all CBL cases simulated. However, $e\xaf\u2061(z)$ in Fig. 4 show very small vertical gradients in the mixed layer. Close examination reveals that $\u2202e\xaf/\u2202z$ is even positive for case B1 in the lower half of the CBL. Comparing Fig. 9a to Fig. 4 confirms the presence of nongradient vertical fluxes of TKE, especially for the buoyancy-dominated CBLs. It is also interesting to observe that $wp/\rho 0\xaf$ is directed downward and almost compensates for the changes of $we\xaf$ such that their sum remains nearly universal until −*z*_{i}/*L* drops below 5 in case S15.

Vertical profiles of the diagnosed gradient $\u2212Ke\u2202e\xaf/\u2202z$ and nongradient *K*_{e}*γ*_{e} fluxes are presented in Fig. 9b. The nongradient fluxes clearly dominate over their gradient counterparts except for case S15 in the lower half of the CBL. In the mixed layer, vertical fluxes are almost entirely nongradient in nature. Case S15 aside, *K*_{e}*γ*_{e} profiles for the B and S10 cases are close enough to warrant a single curve representation without the need to account for minor variabilities with −*z*_{i}/*L*. A second-order polynomial

is proposed and presented as the gray line in Fig. 9b. For case S15, $\u2202e\xaf/\u2202z$ becomes significant from the surface up to ~0.4*z*_{i} as shown in Fig. 4, and so is the gradient flux. However, over this region, *K*_{e}*γ*_{e} for S15 shows no significant difference from the rest of the cases. Moderate deviations of *K*_{e}*γ*_{e} in S15 are found between 0.4 and 0.8*z*_{i} in the upper mixed layer. In the entrainment zone, nongradient fluxes of all CBL cases again overlap quite well. This means that Eq. (34) should serve as a good model for *K*_{e}*γ*_{e} for the CBL in general, except for the strongly sheared conditions, under which the nongradient $we+wp/\rho 0\xaf$ fluxes would be underestimated in upper mixed layer. The impact of such underestimation is ameliorated to some extent due to the diminishing contribution of the turbulence and pressure transport terms in the TKE budget as −*z*_{i}/*L* decreases toward zero (see Fig. 2). Therefore, Eq. (34) is used as an empirical closure for *K*_{e}*γ*_{e}.

## 4. Evaluation and discussion

All improvements proposed in section 3 are implemented into the BouLac scheme in WRF V3.9 to enable a posteriori evaluation. Since the parameterization of momentum flux is not investigated in this study, Pr_{t} = 1 in the original BouLac scheme is left unchanged. The practice of setting Pr_{t} to constant values is also adopted by the ACM2 (Pleim 2007), the MYJ, and the MYNN2 PBL schemes (the first two sets Pr_{t} to 1, and last uses Pr_{t} = 0.74). Sensitivity tests with Pr_{t} = 0.74 following MYNN2 on the idealized cases show that the vertical profiles presented in Figs. 10 and 11 remain unaffected (results not shown). Previous studies have also suggested that the profiles of the daytime CBL show minor sensitivity to Pr_{t} (Noh et al. 2003; Nielsen-Gammon et al. 2010; Yang et al. 2019). The improved scheme is referred to as the new scheme to differentiate from the original BouLac scheme.

### a. Idealized cases

A set of idealized tests are performed with a single column model, following the same setup as the LES cases described in section 2. Two uniform vertical spacings (100 and 50 m) are adopted to investigate model sensitivity to vertical resolution. Three popular PBL schemes, which are the YSU, the MYNN2, and the BouLac schemes, are also included in the tests for comparison purposes. Profiles at the end of the simulations are compared to the LES benchmark.

Figure 10 presents Θ(*z*) for both the low- and high-resolution runs. In the top row, profiles of the coarse resolution new model essentially overlap with the LES benchmark for all CBL cases considered. In comparison, the two higher-order schemes (BouLac and MYNN2) perform well for the two sheared cases, but underpredict Θ(*z*) for the three buoyancy cases. The YSU scheme, on the other hand, works well for B1, but increasingly overpredicts Θ(*z*) as the bulk CBL stability decreases. The model intercomparison demonstrates the robustness of the new scheme in modeling Θ(*z*) across a wide stability range. Profiles of the 50-m runs presented in the bottom row of Fig. 10 are very close to the 100-m ones, suggesting minor sensitivity to vertical resolution.

Vertical profiles of the prognostic TKE $e\xaf\u2061(z)$ are presented in Fig. 11. The profile of $e\xaf\u2061(z)$ predicted by the new model agrees well with the LES benchmark except for the entrainment zone close to the top of the CBL, where $e\xaf\u2061(z)$ is consistently underestimated. For case S15, the vertical extent of the underprediction is extended farther down to the mixed layer. Doubling the vertical resolution, as shown in the bottom row of Fig. 11, raises $e\xaf\u2061(z)$ of the new model closer to that of the LES benchmark, while the “gap” still remains. A number of issues including possibly insufficient turbulent momentum mixing hence shear production, underestimated TKE mixing and overpredicted dissipation could have contributed to such underestimation in the entrainment zone. Even the surface boundary conditions of TKE [i.e., Eq. (32)] could be involved. By raising *c*_{u} and *c*_{w} in Eq. (32), we could achieve good agreement of the predicted $e\xaf\u2061(z)$ in the mixed layer and the entrainment zone, if some overprediction of $e\xaf\u2061(z)$ is allowed in the surface layer. In this study, we do not attempt to fix this issue with ad hoc adjustments, and shall leave this to a future study.

Overall, $e\xaf\u2061(z)$ of the new model is better than that of the BouLac model for the convective cases, while the two models perform similarly well for the sheared cases. The MYNN2 predicted $e\xaf\u2061(z)$ is close to that of the new model for the *B* cases, except for the entrainment zone where $\u2202e\xaf/\u2202z$ is too mild. For the *S* cases, MYNN2 overpredicts $e\xaf\u2061(z)$ throughout the depth of the CBL. Similar to the findings of Fig. 10, the new scheme demonstrates applicability to all CBL cases simulated.

Next, model sensitivities to the changes in *L*_{k}, *L*_{ε}, and *K*_{e}*γ*_{e} are investigated separately. The BouLac scales [Eqs. (14) and (15)] are adopted as the control for their simplicity, upon which modifications are made one at a time. We focus on cases B10 and S15. The former is representative of all the *B* cases where Θ(*z*) is underpredicted by both the BouLac and MYNN2 schemes as shown in Fig. 10. S15 is chosen because it involves more uncertainties related to the parameterization of turbulence and pressure terms as discussed in section 3. Figure 12 presents the low-resolution sensitivity test results for B10 and S15. In the top row, inclusion of the nongradient TKE flux has the most significant impact on Θ(*z*). Modifications of either *L*_{k} or *L*_{ε} only slightly reduces Θ(*z*) from the control case. However, the relative importance of the three terms in predicting Θ(*z*) differs for B10 and S15. In the former, addition of *K*_{e}*γ*_{e} appears to dominate over both length-scale modifications, while for S15, new formulations of *L*_{k} and *L*_{ε} clearly have a stronger influence.

The simulated $e\xaf\u2061(z)$ in the bottom row of Fig. 12 demonstrates the synergy of combined modifications in the new scheme. In Fig. 12c, none of the three modifications alone is able to match $e\xaf\u2061(z)$ to that of the LES benchmark, an issue that is only resolved when modifications are combined. For S15 in Fig. 12d, *L*_{k} and *K*_{e}*γ*_{e} have stronger influence in the upper mixed layer, while *L*_{ε} regulates the new scheme toward the surface layer.

### b. Semi-idealized case

The new scheme is then tested in a semi-idealized setting with realistic initial conditions, large-scale forcings, and diurnally varying bottom boundary conditions. The classic Wangara day 33 case (Clarke et al. 1971) is selected for this purpose. The model setup is described in detail in Zhou et al. (2014) and not repeated here. LES data of Zhou et al. (2014) are also used as benchmark. An SCM with the new PBL scheme is performed with a 100 m vertical spacing. Refining Δ_{z} to 50 m has little effect on first-order mean profiles, but improves higher-order statistics (results not shown).

Figure 13 presents the vertical profiles of a few selected first-order (top row) and second-order (bottom row) moments. Compared with LES, the SCM is able to reproduce the evolution of the daytime potential temperature Θ(*z*) (Fig. 13a) and horizontal wind speed *U*(*z*) (Fig. 13b). The heat flux $w\theta \xaf\u2061(z)$ is also predicted quite well in Fig. 13c, except for a somewhat deeper region of negative entrainment flux. This is primarily due to the coarse 100-m spacing, and goes away when resolution is refined. The TKE profiles in Fig. 13d reveal the same underestimation issue in the upper mixed layer, as shown in the idealized test results in Fig. 11. Doubling the vertical resolution extends $e\xaf\u2061(z)$ upward and improves the agreement, but small vertical gaps between the SCM and LES profiles toward the CBL top remain, nevertheless.

### c. Real case

Next, the new scheme is put to a real-case test. The case study of Huang et al. (2018) during 6 and 7 July 2015 over northern China is selected for the availability of model setup as well as observation data. Weather conditions were cloudy on 6 and clear on 7 July. Two one-way nested domains with 154 × 154 × 37 grid points are centered at 40.24°N, 116.45°E in Beijing, China. The horizontal grid spacings are 27 and 9 km, respectively. The vertical grid is stretched to include about 15 grid points in the lowest 1 km. The initial and boundary conditions are from the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) 1° × 1° final analysis (FNL). The model is initialized at 12 UTC 5 July and is integrated for 60 h. The same four PBL schemes used in the idealized tests are again adopted for the real case study. The rest of the physical parameterizations follow Huang et al. (2018) exactly. Details can be found in their section 2.1 and are not repeated here. Radiosondes released at 0715, 1315, and 1915 LST at Nanjiao station (39.80°N, 116.47°E) close to the model domain center provide the observations.

Figure 14 presents the observed and simulated Θ(*z*) over Nanjiao station at the three radiosonde launch times on 6 (top row) and 7 July (bottom row). The modeled Θ(*z*) show a general warm bias on 6 July under cloudy conditions, while the overall agreement with observations is much better on 7 July under clear sky. The relative performance among the four PBL schemes is similar for both days. Shortly after sunrise, the MYNN2 scheme produces a colder Θ(*z*) than the other three schemes whose results are closely clustered together (Figs. 14a,d). Modeled Θ(*z*) from different PBL schemes become separated at midday as shown in Figs. 14b and 14e. MYNN2 continues to produce the coldest Θ(*z*) and lowest PBL depth *z*_{i} judging by the height of the capping inversion. YSU produces the warmest Θ(*z*) and the deepest boundary layer. Profiles of Θ(*z*) of the new and the BouLac schemes fall in-between that of the MYNN2 and YSU schemes, with the new scheme predicting slightly warmer Θ(*z*) in the CBL in better agreement with observations. Also notice that the vertical gradient of Θ within the CBL is less steep with the new scheme than with BouLac, suggesting more vigorous vertical mixing in the former. Toward sunset, profiles of all four schemes come close together with comparable *z*_{i} on 6 July in Fig. 14c. On 7 July in Fig. 14f, while Θ(*z*) of the new and the BouLac schemes agree well with observations, the MYNN2 and YSU schemes are both subjected to cold biases up to ~0.5 and ~1 K, respectively. Wind speed profiles exhibit similar patterns among the four schemes. YSU and MYNN2 produce the strongest and the weakest CBL winds, while wind speeds of the BouLac and the new schemes are somewhere in-between (results not shown).

On 7 July, the new scheme clearly outperforms the rest with its superior agreement with observations. However, such agreement should not be overrated since many other factors such as synoptic and land surface forcings could have reduced large-scale biases leading to improved model results. This is most clearly demonstrated in 6 July, where all four PBL schemes suffer from warm biases. The most important conclusion to draw from the real case study is that the performance of the new scheme is comparable to the other three well-established schemes for daytime convective conditions. Overall, the predicted Θ(*z*) of the new scheme is warmer than that of the MYNN2 and BouLac schemes and colder than that of the YSU scheme, same as the general patterns found in the idealized tests.

## 5. Summary and future work

Based on a priori analysis of idealized LES of the CBL, new mixing and dissipation lengths (*L*_{k} and *L*_{ε}) for 1.5-order PBL schemes are proposed. *L*_{k} is constructed based on a combination of Redelsperger et al.’s surface layer scale accounting for both surface similarity and TKE balance, and Bougeault and Lacarrere’s (1989) parcel scale with modifications to correct for the reduction of mixing length due most likely to lateral entrainment. The *L*_{ε} adopts the process-based concept of MS94 to model dissipation based on balanced TKE budgets. Despite the relative simplicity in the representation of the TKE forcings, the modeled *L*_{ε} agrees well with the LES benchmark. Furthermore, unlike conventional scaling-based *L*_{ε}, process-based *L*_{ε} also exhibits the same functional dependence on CBL bulk stability as that diagnosed from LES. In addition, inclusion of a nongradient TKE flux is also proposed to improve the representation of TKE transport based on LES analysis.

The improved length scales are then implemented into the BouLac scheme to form a new scheme, and is evaluated with a posteriori tests. For the idealized tests, profiles of Θ(*z*) and $e\xaf\u2061(z)$ of the new scheme exhibit favorable agreements with the LES benchmark. Most notably, Θ(*z*) is well predicted universally across a stability range from free-convective to near neutral conditions. In comparison, predictions from three other classic PBL schemes are biased at either end of the CBL stability spectrum. Sensitivity studies further reveal the synergy of the two new length scales and the nongradient TKE flux in the new scheme. In the real-case test, a cloudy and a clear-weather CBL are simulated on a 9-km WRF grid. Profiles of temperature and wind speed of the new scheme generally fall in-between the predicted profiles of the well-established MYNN2 and YSU schemes. Profiles of Θ(*z*) of the new scheme compare well to radiosonde observations at sunrise, midday, and sunset, especially on the clear day.

As laid out by Plant and Yano (2016), parameterizations in general should be viewed as a fundamental science problem that demands physical insight. With that in mind, the purpose of this study is not to “invent” yet another scaling-based PBL scheme, but to try to include as many physical processes and constraints in one as possible, in the hope that the new scheme can give the right predictions for the right reasons. This study is a first step toward building a more physically based PBL scheme. It demonstrates the potential improvements of adopting process-based/constrained length scales over scaling-based ones. At this point, the new scheme is still incomplete even for the CBL. Future plans include investigation of turbulent mixing of momentum, TKE, and other scalars. Studies on the dynamics of the entrainment zone, especially on the elevated shear production and its proper parameterization, are also planned.

## Acknowledgments

We thank Prof. Jianning Sun, Prof. Ming Xue, and Prof. Zhe-Min Tan for helpful discussions. This work was supported by the National Key R&D Program of China (Grants 2018YFC1506802 and 2018YFC1507604). Support from the National Natural Science Foundation of China through Grant 41875011 for B. Zhou, Grant 41805011 for Y. Li, and Grant 41975124 for K. Zhu are also gratefully acknowledged. Simulations were performed at the High Performance Computing Center of Nanjing University.

## REFERENCES

*J. Appl. Meteor.*

*Mon. Wea. Rev.*

*J. Climate*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

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

*Geophys. Res. Lett.*

*The Atmospheric Boundary Layer*

*Mon. Wea. Rev.*

*J. Fluid Mech.*

*J. Adv. Model. Earth Syst.*

*Rev. Geophys.*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*Meteor. Atmos. Phys.*

*Eur. J. Mech.*

*Proc. Summer Program*, Stanford University Center for Turbulence Research and NASA, 179–184

*Bound.-Layer Meteor.*

*Atmospheric Boundary Layer Flows: Their Structure and Measurement*

*J. Atmos. Sci.*

*Geosci. Model Dev.*

*A Century of Progress in Atmospheric and Related Sciences: Celebrating the American Meteorological Society Centennial*,

*Meteor. Monogr.*, No. 59, Amer. Meteor. Soc.

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

*Bound.-Layer Meteor.*

*Atmos. Res.*

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

*Adv. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*Bound.-Layer Meteor.*

*J. Meteor. Soc. Japan*

*Mon. Wea. Rev.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*Acta Geophys.*

*J. Atmos. Sci.*

*Theoretical Background and Formulation*. Vol. 1,

*Parameterization of Atmospheric Convection*, Imperial College Press, 515 pp

*J. Appl. Meteor. Climatol.*

*Bound.-Layer Meteor.*

*Mon. Wea. Rev.*

*Bound.-Layer Meteor.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*An Introduction to Boundary Layer Meteorology*

*Int. J. Heat Fluid Flow*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*Annu. Rev. Fluid Mech.*

*J. Atmos. Sci.*

*Bound. -Layer Meteor.*

*J. Geophys. Res.*

*Mon. Wea. Rev.*

*Meteor. Atmos. Phys.*

*Meteor. Atmos. Phys.*

*J. Geophys. Res. Atmos.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

## Footnotes

^{1}

The original expression in Wyngaard and Cote (1974) is formulated as a surface similarity function where $es\xaf/u*=f\u2061(zi/L)$, but is equivalent to Eq. (32) using the relation $\u2212zi/L=\kappa w*3/u*3$.