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,

 
wϕ¯=KΦz,
(1)

where the overline represents the Reynolds-average operator, upper- and lowercase letters represent Reynolds-averaged and turbulent variables, respectively. The vertical turbulent flux wϕ¯ is parameterized as a product of the local gradient of resolved variable Φ and an eddy diffusion coefficient K (m2 s−1). In 1.5-order schemes, K is further parameterized as

 
K=e¯Lk,
(2)

where e¯(1/2)u2+υ2+w2¯ (m2 s−2) is the TKE, and Lk (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,

 
ε=e3/2Lε,
(3)

where ε (m2 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 Lk 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 Lk 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 Lk 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 −zi/L, which is interpreted as a bulk measure of the relative importance of thermal to mechanical forcings (Lenschow and Stephens 1980), where zi (m) is the boundary layer depth, Lu*3Θ0/(κgwθ¯s) (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θ¯s (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 z0 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.

Table 1.

List of model parameters and the bulk stability parameter −zi/L averaged from the last 2 h of the simulations. The initial letter in the case name indicates the prescribed heat fluxes, which are 0.20, 0.05, and 0 K m s−1 for the B, S, and N cases, respectively. The numerical value at the end of the case name indicates the prescribed geostrophic wind speeds in m s−1.

List of model parameters and the bulk stability parameter −zi/L averaged from the last 2 h of the simulations. The initial letter in the case name indicates the prescribed heat fluxes, which are 0.20, 0.05, and 0 K m s−1 for the B, S, and N cases, respectively. The numerical value at the end of the case name indicates the prescribed geostrophic wind speeds in m s−1.
List of model parameters and the bulk stability parameter −zi/L averaged from the last 2 h of the simulations. The initial letter in the case name indicates the prescribed heat fluxes, which are 0.20, 0.05, and 0 K m s−1 for the B, S, and N cases, respectively. The numerical value at the end of the case name indicates the prescribed geostrophic wind speeds in m s−1.

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¯ is obtained by solving the Reynolds-averaged prognostic TKE equation,

 
De¯Dt=gΘ0wθ¯uw¯Uzυw¯Vzwe¯z1ρ0wp¯zε.
(4)

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¯ and wp¯/ρ0 are usually grouped together and parameterized as a turbulent diffusion term (Stull 1988, chapter 6.5),

 
we¯+1ρ0wp¯=Kee¯z,
(5)

where Ke 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 KM and the eddy diffusivity for heat KH are linked by a turbulent Prandtl number Pr ≡ KM/KH, which is approximately 1 for ABL flows. Eddy diffusivities of other scalars including moisture and trace gases are often set equal to KH (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 Lk strictly applies to potential temperature.

Parameterization of the vertical heat flux wθ¯ 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),

 
wθ¯=KH(Θzγ).
(6)

Inclusion of KHγ in Eq. (6) effectively counteracts the downward (negative) gradient-diffusion flux −KH∂Θ/∂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,

 
γ=2w*wθ¯sw2¯zi,
(7)

where w*(g/Θ0wθ¯szi)1/3 (m s−1) is the convective velocity scale. Based on Eqs. (2) and (6),

 
Lk=wθ¯e¯(Θzγ).
(8)

Substituting in the expression of γ from Eq. (7), Lk 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 Lk/zi are presented in Fig. 1a. The parameter zi is diagnosed as the height of minimum wθ¯. Lk/zi exhibits distinct patterns above and below ~0.4zi. In the lower half of the CBL, apparent stability dependence is observed as Lk/zi decreases monotonically with decreasing bulk CBL stability measured by −zi/L. Above ~0.4zi, Lk/zi 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, Lk becomes vanishingly small.

Fig. 1.

Normalized vertical profiles of the horizontally and time-averaged turbulent mixing length Lk for the five CBL cases. Solid lines represent LES-diagnosed Lk based on Eq. (8) in all panels. (b) Dashed lines represent Redelsperger’s surface length scale Lks in Eq. (9). (c) Dashed lines represent the BouLac length scales CkBLlu and CkBLld in Eqs. (12) and (13), while the dash–dotted lines represent Lkm in Eq. (16). (d) Dashed lines represent the modeled Lk based on Eq. (18). Note that the LES-diagnosed Lk/zi show large wiggles around 0.8zi, because both the numerator and denominator of Eq. (8) are close to zero at that height. Those presented here have been smoothed between 0.7 and 0.9zi.

Fig. 1.

Normalized vertical profiles of the horizontally and time-averaged turbulent mixing length Lk for the five CBL cases. Solid lines represent LES-diagnosed Lk based on Eq. (8) in all panels. (b) Dashed lines represent Redelsperger’s surface length scale Lks in Eq. (9). (c) Dashed lines represent the BouLac length scales CkBLlu and CkBLld in Eqs. (12) and (13), while the dash–dotted lines represent Lkm in Eq. (16). (d) Dashed lines represent the modeled Lk based on Eq. (18). Note that the LES-diagnosed Lk/zi show large wiggles around 0.8zi, because both the numerator and denominator of Eq. (8) are close to zero at that height. Those presented here have been smoothed between 0.7 and 0.9zi.

The general profiles of Lk suggest a three-layer treatment that comprises an extended surface layer where Lk varies with the surface scale z/L, an upper mixed layer with a universal mixing length that scales with zi, 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

 
Lks=AzϕL(z/L),
(9)

where A = 2.79 is a constant, and ϕL is a MO stability function that also depends on the dimensionless shear ϕm(κz/u*)U/z. For unstable conditions (i.e., z/L < 0),

 
ϕL=ϕm1[1+1α(zL)2/3]1/2,
(10)
 
ϕm=(115zL)1/4.
(11)

The coefficient α in Eq. (10) represents part of the surface layer TKE parameterized by u*2 in the form of e¯=αu*2+βw*2 (Wyngaard and Cote 1974; Panofsky et al. 1977; Garratt 1992). A similar expression e¯=f(u*,w*) 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 Lk 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.

Fig. 2.

Vertical profiles of the horizontally and time-averaged normalized TKE forcing terms for (a)–(e) the five CBL cases and (f) the NBL case. The CBL and NBL forcings are normalized by w*3/zi and u*3/zi, respectively. “Resi” stands for residual in the legends. The dissipation rate ε is parameterized by the TKE-1.5 closure.

Fig. 2.

Vertical profiles of the horizontally and time-averaged normalized TKE forcing terms for (a)–(e) the five CBL cases and (f) the NBL case. The CBL and NBL forcings are normalized by w*3/zi and u*3/zi, respectively. “Resi” stands for residual in the legends. The dissipation rate ε is parameterized by the TKE-1.5 closure.

In the upper CBL where Lk/zi exhibits universality, the Bougeault–Lacarrere (BouLac; Bougeault and Lacarrere 1989) length scales are adopted as the base scales,

 
zz+lugΘ0[Θ(z)Θ(z)]dz=e¯(z),
(12)
 
zldzgΘ0[Θ(z)Θ(z)]dz=e¯(z),
(13)

where lu and ld represent the upward and downward distances that an undiluted air parcel at height z with a mean TKE e¯(z) can travel against stratification. They can be derived based on a balance between buoyancy and vertical advection of w. In Eq. (13)ld is further limited by the distance to ground. The BouLac mixing length LkBL and dissipation length LεBL are then formulated as

 
LkBL=CkBLmin(lu,ld),
(14)
 
LεBL=CεBLluld,
(15)

where CkBL=0.4, Cε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), lu and ld are insensitive to the bulk stability of the CBL. This is clearly seen in Fig. 1c, where lu and ld for all five CBL cases collapse such that they are well approximated by 1 − z/zi and z/zi, 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, lu and ld are essentially limited by the CBL capping inversion and the ground. The fact that the diagnosed Lk 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 Lk. Here, we only make use of lu in the upper CBL, as implied by taking the lower value between lu and ld in the original BouLac scheme according to Eq. (14).

Clearly, some modifications must be made before adopting lu as it consistently overpredicts Lk 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 lu through an empirically fitted piecewise linear function,

 
Lkm=lu×max{min[0.8,1.5×(1.1z/zi)],0},
(16)

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 Lk above ~0.4zi.

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

 
Lksm=[(Lks)a+(Lkm)a]1/a.
(17)

The coefficient a is set to −3 based on empirical fitting to the LES-diagnosed Lk.

In the entrainment zone and free troposphere, approximately 0.8zi 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 Lk,

 
Lk={Lksmz/zi<0.8max(Lksm,Lke)z/zi0.8.
(18)

Figure 1d presents the parameterized and the diagnosed profiles of Lk. In general, Eq. (18) is able to reproduce the diagnosed Lk 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 Lk 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 Lk. Preliminary tests on the structural and parametric sensitivities of Lk 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 Lk 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),

 
Lε0e¯zdz0e¯dz,
(19)

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

 
Lε(1κz+1λ0)1,
(20)

with λ0 set proportional to zi.

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.9zi. At the bottom of the entrainment zone around 0.8zi, 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.

Fig. 3.

Normalized vertical profiles of the horizontally and time-averaged turbulent dissipation length Lε for the CBL and NBL cases. (a) Solid lines represent diagnosed Lε from LES. (b) Solid and dashed lines represent modeled Lε by the MYNN and the BouLac schemes, respectively.

Fig. 3.

Normalized vertical profiles of the horizontally and time-averaged turbulent dissipation length Lε for the CBL and NBL cases. (a) Solid lines represent diagnosed Lε from LES. (b) Solid and dashed lines represent modeled Lε by the MYNN and the BouLac schemes, respectively.

Among the five CBL cases, Lε/zi decreases in the surface layer and increases in the mixed layer as the bulk CBL stability decreases (i.e., −zi/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ε/zi 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 zi/L for the CBL were extrapolated to neutral conditions at zi/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¯ 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(zi), 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ε/zi 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ε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ε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εM/zi show little variations with bulk CBL stability in the upper CBL. In fact, Lε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¯(z) with −zi/L. Figure 4 presents the averaged e¯(z) of all the simulated cases in the descending order of −zi/L. With reduced buoyancy and increased shear production, which is most significant in the surface layer, e¯(z) becomes more surface heavy as −zi/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εBL are essentially invariant with −zi/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).

Fig. 4.

Normalized vertical profiles of the horizontally and time-averaged TKE profiles. TKE is normalized by w*2 and u*2 for the CBL and NBL cases respectively. The NBL TKE is also multiplied by 0.5 for ease of comparison. Angle brackets represent horizontal averaging.

Fig. 4.

Normalized vertical profiles of the horizontally and time-averaged TKE profiles. TKE is normalized by w*2 and u*2 for the CBL and NBL cases respectively. The NBL TKE is also multiplied by 0.5 for ease of comparison. Angle brackets represent horizontal averaging.

The limited and even opposite variability of the modeled Lε/zi 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

 
ε=gΘ0wθ¯uw¯Uzυw¯Vzwe¯z1ρ¯wp¯z.
(21)

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 zi 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θ¯, exhibit the classic linear profiles and are well approximated by

 
gΘ0wθ¯=w*3zi[1(RH+1)zzi],
(22)

where RHwθ¯e/wθ¯s 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 RH with respect to bulk CBL stability (Conzemius and Fedorovich 2006; Liu et al. 2016) is not considered here.

Fig. 5.

(a) Normalized vertical profiles of the horizontally and time-averaged buoyancy production (solid lines) and the combined turbulent mixing and pressure correlation terms (dashed lines) of the TKE budget. (b) The sum of the three forcing terms along with their model [Eq. (24)].

Fig. 5.

(a) Normalized vertical profiles of the horizontally and time-averaged buoyancy production (solid lines) and the combined turbulent mixing and pressure correlation terms (dashed lines) of the TKE budget. (b) The sum of the three forcing terms along with their model [Eq. (24)].

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,

 
we¯z1ρ0wp¯z=w*3zi(0.95zzi0.52).
(23)

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

Adding Eqs. (22) and (23) together, we arrive at

 
gΘ0wθ¯we¯z1ρ0wp¯z=f(zzi)×w*3zi(0.480.25zzi),
(24)

where a step function like f(z/zi) = 1/2{1 − tanh[(z/zi − 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/κ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.

Fig. 6.

Normalized vertical profiles of the horizontally and time-averaged shear production. Solid lines present LES-diagnosed profiles, and dashed lines represent modeled profiles by Eq. (25).

Fig. 6.

Normalized vertical profiles of the horizontally and time-averaged shear production. Solid lines present LES-diagnosed profiles, and dashed lines represent modeled profiles by Eq. (25).

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

 
uw¯Uzυw¯Vz=u*3κz(1zzi)ϕm,
(25)

based on the fact that the momentum flux profile τuw¯2+υw¯2 in the CBL is well represented by u*2(1z/zi) (Zhou et al. 2019), and by adopting a simple approximation of the wind shear by extending surface layer similarity upward (i.e., U/z=u*/κzϕ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,

 
1Lε=1Lεs+1Lεb,
(26)

where

 
Lεs=κzϕm(1z/zi)(e¯u*)3,
(27)
 
Lεb=fzi0.480.25z/zi(e¯w*)3
(28)

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 lbe¯/N, where N (s−1) is the buoyancy frequency, is often adopted by 1.5-order PBL schemes. While lb 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¯/(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εe based on the buoyancy lb and shear lse¯/S length scales,

 
1Lεe=1cblb+1csls,
(29)

where S=(U/z)2+(V/z)2 (s−1) is the mean shear, cs and cb are constants. Equation (29) is constructed such that Lε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 lb and ls in Eq. (29) is also affected by the ratio of the model constants rbscb/cs. ls dominates only when csls < cblb, or ls/lb < rbs. Given the definitions of ls and lb, their ratio ls/lb=Rig is the square root of the gradient Richardson number RigN2/S2. Therefore, rbs can be interpreted as a threshold Rig 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 Rif or Rig, 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, rbs is set to 1. The coefficient cb = 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, rbs 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 rbs to 1, but will be left to future work.

With the new entrainment dissipation scale, the complete model of Lε is formulated as

 
1Lε=1Lεb+1Lεs+1Lεe,
(30)

where Lε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 −zi/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.

Fig. 7.

Normalized vertical profiles of the modeled and LES-diagnosed dissipation length for the CBL and NBL cases.

Fig. 7.

Normalized vertical profiles of the modeled and LES-diagnosed dissipation length for the CBL and NBL cases.

Fig. 8.

Normalized vertical profiles of the modeled dissipation length (MS) and its shear (Shear), buoyancy (Buoy) and entrainment (Entrain) components for case S10. The LES benchmark is also included as the black solid line.

Fig. 8.

Normalized vertical profiles of the modeled dissipation length (MS) and its shear (Shear), buoyancy (Buoy) and entrainment (Entrain) components for case S10. The LES benchmark is also included as the black solid line.

In the a posteriori tests, Lεs and Lε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,

 
ε=e¯3/2LεR,R=(es¯e1¯)3/2,
(31)

where

 
es¯3/2=(cuu*)3+(cww*)3
(32)

is a parameterized surface-layer TKE (es¯) based on the friction and the convective velocities, and e1¯ 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¯, where cu = 2.3 and cw = 0.59 are model constants.1Troen and Mahrt (1986) also adopted Eq. (32) for their PBL scheme using cu ≈ 1.0 and cw ≈ 0.65. In this study, cu and cw are set to 1.5 and 0.67 respectively based on a linear regression of LES diagnosed e¯3/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¯, 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θ¯ in Eq. (6) but for turbulent flux of TKE (we¯). Based on the budget analysis of we¯, 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,

 
we¯+1ρ0wp¯=Ke(e¯zγe),
(33)

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 Ke = KH, 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¯ and wp/ρ0¯. Up to the entrainment zone, we¯ and wp/ρ0¯ are of opposite sign. we¯ is directed upward for all CBL cases simulated. However, e¯(z) in Fig. 4 show very small vertical gradients in the mixed layer. Close examination reveals that e¯/z 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/ρ0¯ is directed downward and almost compensates for the changes of we¯ such that their sum remains nearly universal until −zi/L drops below 5 in case S15.

Fig. 9.

Normalized vertical profiles of the horizontally and time-averaged vertical fluxes of TKE and pressure for the CBL cases. (a) Dashed and dash–dotted lines represent ⟨we⟩ and ⟨wp/ρo⟩, respectively, while the solid lines represent their sum. (b) Dashed and dash–dotted lines represent nongradient and gradient TKE fluxes according to Eq. (33). The gray line represents the parameterized nongradient flux of Eq. (34).

Fig. 9.

Normalized vertical profiles of the horizontally and time-averaged vertical fluxes of TKE and pressure for the CBL cases. (a) Dashed and dash–dotted lines represent ⟨we⟩ and ⟨wp/ρo⟩, respectively, while the solid lines represent their sum. (b) Dashed and dash–dotted lines represent nongradient and gradient TKE fluxes according to Eq. (33). The gray line represents the parameterized nongradient flux of Eq. (34).

Vertical profiles of the diagnosed gradient Kee¯/z and nongradient Keγ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, Keγ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 −zi/L. A second-order polynomial

 
Keγe=0.4w*3zzi(zzi1)
(34)

is proposed and presented as the gray line in Fig. 9b. For case S15, e¯/z becomes significant from the surface up to ~0.4zi as shown in Fig. 4, and so is the gradient flux. However, over this region, Keγe for S15 shows no significant difference from the rest of the cases. Moderate deviations of Keγe in S15 are found between 0.4 and 0.8zi 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 Keγe for the CBL in general, except for the strongly sheared conditions, under which the nongradient we+wp/ρ0¯ 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 −zi/L decreases toward zero (see Fig. 2). Therefore, Eq. (34) is used as an empirical closure for Keγ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, Prt = 1 in the original BouLac scheme is left unchanged. The practice of setting Prt to constant values is also adopted by the ACM2 (Pleim 2007), the MYJ, and the MYNN2 PBL schemes (the first two sets Prt to 1, and last uses Prt = 0.74). Sensitivity tests with Prt = 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 Prt (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.

Fig. 10.

Vertical profiles of potential temperature for the five CBL cases for the (a)–(e) 100-m low-resolution and (f)–(j) 50-m high-resolution results.

Fig. 10.

Vertical profiles of potential temperature for the five CBL cases for the (a)–(e) 100-m low-resolution and (f)–(j) 50-m high-resolution results.

Fig. 11.

As in Fig. 10, but for TKE.

Fig. 11.

As in Fig. 10, but for TKE.

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¯(z) are presented in Fig. 11. The profile of e¯(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¯(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¯(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 cu and cw in Eq. (32), we could achieve good agreement of the predicted e¯(z) in the mixed layer and the entrainment zone, if some overprediction of e¯(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¯(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¯(z) is close to that of the new model for the B cases, except for the entrainment zone where e¯/z is too mild. For the S cases, MYNN2 overpredicts e¯(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 Lk, Lε, and Keγ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 Lk 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 Keγe appears to dominate over both length-scale modifications, while for S15, new formulations of Lk and Lε clearly have a stronger influence.

Fig. 12.

Vertical profiles of (a),(b) potential temperature and (c),(d) TKE for cases (left) B10 and (right) S10. In the legends, “CTRL” stands for the control case, and LK, Lε, and Keγe refer to the sensitivity tests with sole modification of the mixing length, the dissipation length, and the nongradient TKE fluxes, respectively. The x axes in (a) and (b) are set differently for clarity of presentation.

Fig. 12.

Vertical profiles of (a),(b) potential temperature and (c),(d) TKE for cases (left) B10 and (right) S10. In the legends, “CTRL” stands for the control case, and LK, Lε, and Keγe refer to the sensitivity tests with sole modification of the mixing length, the dissipation length, and the nongradient TKE fluxes, respectively. The x axes in (a) and (b) are set differently for clarity of presentation.

The simulated e¯(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¯(z) to that of the LES benchmark, an issue that is only resolved when modifications are combined. For S15 in Fig. 12d, Lk and Keγ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θ¯(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¯(z) upward and improves the agreement, but small vertical gaps between the SCM and LES profiles toward the CBL top remain, nevertheless.

Fig. 13.

Vertical profiles of (a) potential temperature, (b) horizontal wind speed, (c) vertical heat flux, and (d) TKE at five selected times of the Wangara day 33 case. The solid line represent horizontally averaged LES statistics, and dashed lines represent single column results of the new PBL scheme.

Fig. 13.

Vertical profiles of (a) potential temperature, (b) horizontal wind speed, (c) vertical heat flux, and (d) TKE at five selected times of the Wangara day 33 case. The solid line represent horizontally averaged LES statistics, and dashed lines represent single column results of the new PBL scheme.

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

Fig. 14.

Vertical profiles of potential temperature at (left column) 0715, (middle column) 1315, and (right column) 1915 LST (a)–(c) 6 and (d)–(f) 7 Jul 2015 at Nanjiao Station, Beijing, China. Gray dots present radiosonde observations, and solid lines represent 9-km WRF results with different PBL schemes.

Fig. 14.

Vertical profiles of potential temperature at (left column) 0715, (middle column) 1315, and (right column) 1915 LST (a)–(c) 6 and (d)–(f) 7 Jul 2015 at Nanjiao Station, Beijing, China. Gray dots present radiosonde observations, and solid lines represent 9-km WRF results with different PBL schemes.

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 (Lk and Lε) for 1.5-order PBL schemes are proposed. Lk 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¯(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

REFERENCES
Bélair
,
S.
,
J.
Mailhot
,
J. W.
Strapp
, and
J. I.
MacPherson
,
1999
:
An examination of local versus nonlocal aspects of a TKE-based boundary layer scheme in clear convective conditions
.
J. Appl. Meteor.
,
38
,
1499
1518
, https://doi.org/10.1175/1520-0450(1999)038<1499:AEOLVN>2.0.CO;2.
Bougeault
,
P.
, and
P.
Lacarrere
,
1989
:
Parameterization of orography-induced turbulence in a mesobeta-scale model
.
Mon. Wea. Rev.
,
117
,
1872
1890
, https://doi.org/10.1175/1520-0493(1989)117<1872:POOITI>2.0.CO;2.
Bretherton
,
C. S.
, and
S.
Park
,
2009
:
A new moist turbulence parameterization in the Community Atmosphere Model
.
J. Climate
,
22
,
3422
3448
, https://doi.org/10.1175/2008JCLI2556.1.
Cheng
,
Y.
,
V. M.
Canuto
, and
A. M.
Howard
,
2002
:
An improved model for the turbulent PBL
.
J. Atmos. Sci.
,
59
,
1550
1565
, https://doi.org/10.1175/1520-0469(2002)059<1550:AIMFTT>2.0.CO;2.
Clarke
,
R. H.
,
A. J.
Dyer
,
R. R.
Brook
,
D. G.
Reid
, and
A. J.
Troup
,
1971
:
The Wangara experiment: Boundary layer data. CSIRO Tech. Paper 19, 362 pp
.
Conzemius
,
R. J.
, and
E.
Fedorovich
,
2006
:
Dynamics of sheared convective boundary layer entrainment. Part I: Methodological background and large-eddy simulations
.
J. Atmos. Sci.
,
63
,
1151
1178
, https://doi.org/10.1175/JAS3691.1.
Cuxart
,
J.
,
P.
Bougeault
, and
J.-L.
Redelsperger
,
2000
:
A turbulence scheme allowing for mesoscale and large-eddy simulations
.
Quart. J. Roy. Meteor. Soc.
,
126
,
1
30
, https://doi.org/10.1002/qj.49712656202.
Drobinski
,
P.
,
J.-L.
Redelsperger
, and
C.
Pietras
,
2006
:
Evaluation of a planetary boundary layer subgrid-scale model that accounts for near-surface turbulence anisotropy
.
Geophys. Res. Lett.
,
33
,
L23806
, https://doi.org/10.1029/2006GL027062.
Garratt
,
J. R.
,
1992
:
The Atmospheric Boundary Layer
.
Cambridge University Press
,
334
pp.
Grenier
,
H.
, and
C. S.
Bretherton
,
2001
:
A moist PBL parameterization for large-scale models and its application to subtropical cloud-topped marine boundary layers
.
Mon. Wea. Rev.
,
129
,
357
377
, https://doi.org/10.1175/1520-0493(2001)129<0357:AMPPFL>2.0.CO;2.
Haghshenas
,
A.
, and
J. P.
Mellado
,
2019
:
Characterization of wind-shear effects on entrainment in a convective boundary layer
.
J. Fluid Mech.
,
858
,
145
183
, https://doi.org/10.1017/jfm.2018.761.
Heinze
,
R.
,
D.
Mironov
, and
S.
Raasch
,
2015
:
Second-moment budgets in cloud topped boundary layers: A large-eddy simulation study
.
J. Adv. Model. Earth Syst.
,
7
,
510
536
, https://doi.org/10.1002/2014MS000376.
Holt
,
T.
, and
S.
Raman
,
1988
:
A review and comparative evaluation of multilevel boundary layer parameterizations for first-order and turbulent kinetic energy closure schemes
.
Rev. Geophys.
,
26
,
761
780
, https://doi.org/10.1029/RG026i004p00761.
Holtslag
,
A. M.
, and
C.-H.
Moeng
,
1991
:
Eddy diffusivity and countergradient transport in the convective atmospheric boundary layer
.
J. Atmos. Sci.
,
48
,
1690
1698
, https://doi.org/10.1175/1520-0469(1991)048<1690:EDACTI>2.0.CO;2.
Hu
,
X.-M.
,
M.
Xue
, and
X.
Li
,
2019
:
The use of high-resolution sounding data to evaluate and optimize nonlocal PBL schemes for simulating the slightly stable upper convective boundary layer
.
Mon. Wea. Rev.
,
147
,
3825
3841
, https://doi.org/10.1175/MWR-D-19-0085.1.
Huang
,
M.
,
Z.
Gao
,
S.
Miao
, and
F.
Chen
,
2018
:
Sensitivity of urban boundary layer simulation to urban canopy models and PBL schemes in Beijing
.
Meteor. Atmos. Phys.
,
131
,
1235
1248
, https://doi.org/10.1007/S00703-018-0634-1.
Hunt
,
J. C. R.
, and
J. F.
Morrison
,
2000
:
Eddy structure in turbulent boundary layers
.
Eur. J. Mech.
,
19B
,
673
694
, https://doi.org/10.1016/S0997-7546(00)00129-1.
Hunt
,
J. C. R.
,
P.
Spalart
, and
N.
Mansour
,
1987
:
A general form for the dissipation length scale in turbulent shear flows. Proc. Summer Program, Stanford University Center for Turbulence Research and NASA, 179–184
.
Ito
,
J.
,
H.
Niino
,
M.
Nakanishi
, and
C.-H.
Moeng
,
2015
:
An extension of the Mellor-Yamada model to the terra incognita zone for dry convective mixed layers in the free convection regime
.
Bound.-Layer Meteor.
,
157
,
23
43
, https://doi.org/10.1007/s10546-015-0045-5.
Janjić
,
Z. I.
,
2002
:
Nonsingular implementation of the Mellor–Yamada level 2.5 scheme in the NCEP Meso model. NCEP Office Note 437, 61 pp.
, http://www.emc.ncep.noaa.gov/officenotes/newernotes/on437.pdf.
Kaimal
,
J. C.
, and
J. J.
Finnigan
,
1994
:
Atmospheric Boundary Layer Flows: Their Structure and Measurement
.
Oxford University Press
,
304
pp.
Klemp
,
J. B.
, and
R. B.
Wilhelmson
,
1978
:
The simulation of three-dimensional convective storm dynamics
.
J. Atmos. Sci.
,
35
,
1070
1096
, https://doi.org/10.1175/1520-0469(1978)035<1070:TSOTDC>2.0.CO;2.
Lac
,
C.
, and et al
,
2018
:
Overview of the Meso-NH model version 5.4 and its applications
.
Geosci. Model Dev.
,
11
,
1929
1969
, https://doi.org/10.5194/gmd-11-1929-2018.
LeMone
,
M. A.
, and et al
,
2019
:
100 years of progress in boundary layer meteorology. A Century of Progress in Atmospheric and Related Sciences: Celebrating the American Meteorological Society Centennial, Meteor. Monogr., No. 59, Amer. Meteor. Soc.
, https://doi.org/10.1175/AMSMONOGRAPHS-D-18-0013.1.
Lenderink
,
G.
, and
A. M.
Holtslag
,
2004
:
An updated length-scale formulation for turbulent mixing in clear and cloudy boundary layers
.
Quart. J. Roy. Meteor. Soc.
,
130
,
3405
3427
, https://doi.org/10.1256/qj.03.117.
Lenschow
,
D. H.
, and
P. L.
Stephens
,
1980
:
The role of thermals in the convective boundary layer
.
Bound.-Layer Meteor.
,
19
,
509
532
, https://doi.org/10.1007/BF00122351.
Li
,
D.
,
2019
:
Turbulent Prandtl number in the atmospheric boundary layer—Where are we now?
Atmos. Res.
,
216
,
86
105
, https://doi.org/10.1016/j.atmosres.2018.09.015.
Liu
,
C.
,
E.
Fedorovich
, and
J.
Huang
,
2018
:
Revisiting entrainment relationships for shear-free and sheared convective boundary layers through large-eddy simulations
.
Quart. J. Roy. Meteor. Soc.
,
144
,
2182
2195
, https://doi.org/10.1002/qj.3330.
Liu
,
P.
,
J.
Sun
, and
L.
Shen
,
2016
:
Parameterization of sheared entrainment in a well-developed CBL. Part I: Evaluation of the scheme through large-eddy simulations
.
Adv. Atmos. Sci.
,
33
,
1171
1184
, https://doi.org/10.1007/s00376-016-5208-x.
Mellor
,
G.
, and
T.
Yamada
,
1974
:
Hierarchy of turbulence closure models for planetary boundary-layers
.
J. Atmos. Sci.
,
31
,
1791
1806
, https://doi.org/10.1175/1520-0469(1974)031<1791:AHOTCM>2.0.CO;2.
Moeng
,
C.-H.
,
1984
:
A large-eddy-simulation model for the study of planetary boundary-layer turbulence
.
J. Atmos. Sci.
,
41
,
2052
2062
, https://doi.org/10.1175/1520-0469(1984)041<2052:ALESMF>2.0.CO;2.
Moeng
,
C.-H.
, and
P. P.
Sullivan
,
1994
:
A comparison of shear- and buoyancy-driven planetary boundary layer flows
.
J. Atmos. Sci.
,
51
,
999
1022
, https://doi.org/10.1175/1520-0469(1994)051<0999:ACOSAB>2.0.CO;2.
Muñoz-Esparza
,
D.
,
R. D.
Sharman
, and
J. K.
Lundquist
,
2018
:
Turbulence dissipation rate in the atmospheric boundary layer: Observations and WRF mesoscale modeling during the XPIA field campaign
.
Mon. Wea. Rev.
,
146
,
351
371
, https://doi.org/10.1175/MWR-D-17-0186.1.
Nakanishi
,
M.
,
2001
:
Improvement of the Mellor–Yamada turbulence closure model based on large-eddy simulation data
.
Bound.-Layer Meteor.
,
99
,
349
378
, https://doi.org/10.1023/A:1018915827400.
Nakanishi
,
M.
, and
H.
Niino
,
2009
:
Development of an improved turbulence closure model for the atmospheric boundary layer
.
J. Meteor. Soc. Japan
,
87
,
895
912
, https://doi.org/10.2151/jmsj.87.895.
Nielsen-Gammon
,
J. W.
,
X.-M.
Hu
,
F.
Zhang
, and
J. E.
Pleim
,
2010
:
Evaluation of planetary boundary layer scheme sensitivities for the purpose of parameter estimation
.
Mon. Wea. Rev.
,
138
,
3400
3417
, https://doi.org/10.1175/2010MWR3292.1.
Noh
,
Y.
,
W. G.
Cheon
,
S. Y.
Hong
, and
S.
Raasch
,
2003
:
Improvement of the K-profile model for the planetary boundary layer based on large eddy simulation data
.
Bound.-Layer Meteor.
,
107
,
401
427
, https://doi.org/10.1023/A:1022146015946.
Panofsky
,
H. A.
,
H.
Tennekes
,
D. H.
Lenschow
, and
J. C.
Wyngaard
,
1977
:
The characteristics of turbulent velocity components in the surface layer under convective conditions
.
Bound.-Layer Meteor.
,
11
,
355
361
, https://doi.org/10.1007/BF02186086.
Pino
,
D.
, and
J. V.-G.
De Arellano
,
2008
:
Effects of shear in the convective boundary layer: Analysis of the turbulent kinetic energy budget
.
Acta Geophys.
,
56
,
167
193
, https://doi.org/10.2478/s11600-007-0037-z.
Piper
,
M.
, and
J. K.
Lundquist
,
2004
:
Surface layer turbulence measurements during a frontal passage
.
J. Atmos. Sci.
,
61
,
1768
1780
, https://doi.org/10.1175/1520-0469(2004)061<1768:SLTMDA>2.0.CO;2.
Plant
,
R. S.
, and
J.
Yano
, Eds.,
2016
:
Theoretical Background and Formulation. Vol. 1, Parameterization of Atmospheric Convection, Imperial College Press, 515 pp
.
Pleim
,
J. E.
,
2007
:
A combined local and nonlocal closure model for the atmospheric boundary layer. Part I: Model description and testing
.
J. Appl. Meteor. Climatol.
,
46
,
1383
1395
, https://doi.org/10.1175/JAM2539.1.
Redelsperger
,
J. L.
,
F.
Mahe
, and
P.
Carlotti
,
2001
:
A simple and general subgrid model suitable both for surface layer and free-stream turbulence
.
Bound.-Layer Meteor.
,
101
,
375
408
, https://doi.org/10.1023/A:1019206001292.
Seity
,
Y.
,
P.
Brousseau
,
S.
Malardel
,
G.
Hello
,
P.
Bénard
,
F.
Bouttier
,
C.
Lac
, and
V.
Masson
,
2011
:
The AROME-France convective-scale operational model
.
Mon. Wea. Rev.
,
139
,
976
991
, https://doi.org/10.1175/2010MWR3425.1.
Shin
,
H. H.
, and
S.-Y.
Hong
,
2011
:
Intercomparison of planetary boundary-layer parametrizations in the WRF Model for a single day from CASES-99
.
Bound.-Layer Meteor.
,
139
,
261
281
, https://doi.org/10.1007/s10546-010-9583-z.
Shin
,
H. H.
, and
S.-Y.
Hong
,
2013
:
Analysis of resolved and parameterized vertical transports in convective boundary layers at gray-zone resolutions
.
J. Atmos. Sci.
,
70
,
3248
3261
, https://doi.org/10.1175/JAS-D-12-0290.1.
Shin
,
H. H.
,
S.-Y.
Hong
,
Y.
Noh
, and
J.
Dudhia
,
2013
:
Derivation of turbulent kinetic energy from a first-order nonlocal planetary boundary layer parameterization
.
J. Atmos. Sci.
,
70
,
1795
1805
, https://doi.org/10.1175/JAS-D-12-0150.1.
Siebesma
,
A. P.
,
P. M. M.
Soares
, and
J.
Teixeira
,
2007
:
A combined eddy-diffusivity mass-flux approach for the convective boundary layer
.
J. Atmos. Sci.
,
64
,
1230
1248
, https://doi.org/10.1175/JAS3888.1.
Skamarock
,
W. C.
, and et al
,
2008
:
A description of the Advanced Research WRF version 3. NCAR Tech. Note NCAR/TN-475+STR, 113 pp.
, http://doi.org/10.5065/D68S4MVH.
Stull
,
R. B.
,
1988
:
An Introduction to Boundary Layer Meteorology
.
Kluwer Academic
,
666
pp.
Tardu
,
S.
,
2017
:
Near wall dissipation revisited
.
Int. J. Heat Fluid Flow
,
67
,
104
115
, https://doi.org/10.1016/j.ijheatfluidflow.2017.03.006.
Teixeira
,
J.
, and
S.
Cheinet
,
2004
:
A simple mixing length formulation for the eddy-diffusivity parameterization of dry convection
.
Bound.-Layer Meteor.
,
110
,
435
453
, https://doi.org/10.1023/B:BOUN.0000007230.96303.0d.
Therry
,
G.
, and
P.
Lacarrère
,
1983
:
Improving the eddy kinetic energy model for planetary boundary layer description
.
Bound.-Layer Meteor.
,
25
,
63
88
, https://doi.org/10.1007/BF00122098.
Troen
,
I.
, and
L.
Mahrt
,
1986
:
A simple-model of the atmospheric boundary layer; sensitivity to surface evaporation
.
Bound.-Layer Meteor.
,
37
,
129
148
, https://doi.org/10.1007/BF00122760.
Vassilicos
,
J. C.
,
2015
:
Dissipation in turbulent flows
.
Annu. Rev. Fluid Mech.
,
47
,
95
114
, https://doi.org/10.1146/annurev-fluid-010814-014637.
Wyngaard
,
J. C.
,
2004
:
Toward numerical modeling in the “terra incognita.”
J. Atmos. Sci.
,
61
,
1816
1826
, https://doi.org/10.1175/1520-0469(2004)061<1816:TNMITT>2.0.CO;2.
Wyngaard
,
J. C.
, and
O. R.
Cote
,
1974
:
The evolution of a convective planetary boundary layer? A higher-order-closure model study
.
Bound. -Layer Meteor.
,
7
,
289
308
, https://doi.org/10.1007/bf00240833.
Xie
,
B.
,
J. C. H.
Fung
,
A.
Chan
, and
A.
Lau
,
2012
:
Evaluation of nonlocal and local planetary boundary layer schemes in the WRF Model
.
J. Geophys. Res.
,
117
,
D12103
, https://doi.org/10.1029/2011JD017080.
Xue
,
M.
,
2000
:
High-order monotonic numerical diffusion and smoothing
.
Mon. Wea. Rev.
,
128
,
2853
2864
, https://doi.org/10.1175/1520-0493(2000)128<2853:HOMNDA>2.0.CO;2.
Xue
,
M.
,
K. K.
Droegemeier
, and
V.
Wong
,
2000
:
The Advanced Regional Prediction System (ARPS)—A multi-scale nonhydrostatic atmospheric simulation and prediction model. Part I: Model dynamics and verification
.
Meteor. Atmos. Phys.
,
75
,
161
193
, https://doi.org/10.1007/s007030070003.
Xue
,
M.
, and et al
,
2001
:
The Advanced Regional Prediction System (ARPS)—A multi-scale nonhydrostatic atmospheric simulation and prediction tool. Part II: Model physics and applications
.
Meteor. Atmos. Phys.
,
76
,
143
165
, https://doi.org/10.1007/s007030170027.
Yang
,
B.
, and et al
,
2019
:
Parametric and structural sensitivities of turbine-height wind speeds in the boundary layer parameterizations in the Weather Research and Forecasting Model
.
J. Geophys. Res. Atmos.
,
124
,
5951
5969
, https://doi.org/10.1029/2018JD029691.
Zhou
,
B.
,
J. S.
Simon
, and
F. K.
Chow
,
2014
:
The convective boundary layer in the terra incognita
.
J. Atmos. Sci.
,
71
,
2545
2563
, https://doi.org/10.1175/JAS-D-13-0356.1.
Zhou
,
B.
,
S.
Sun
,
K.
Yao
, and
K.
Zhu
,
2018
:
Reexamining the gradient and countergradient representation of the local and nonlocal heat fluxes in the convective boundary layer
.
J. Atmos. Sci.
,
75
,
2317
2336
, https://doi.org/10.1175/JAS-D-17-0198.1.
Zhou
,
B.
,
S.
Sun
,
J.
Sun
, and
K.
Zhu
,
2019
:
The universality of the normalized vertical velocity variance in contrast to the horizontal velocity variance in the convective boundary layer
.
J. Atmos. Sci.
,
76
,
1437
1456
, https://doi.org/10.1175/JAS-D-18-0325.1.
For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

Footnotes

1

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