Since the 1950s, a countergradient flux term has been added to some K-profile-based first-order PBL schemes, allowing them to simulate the slightly statically stable upper part of the convective boundary layer (CBL) observed in a limited number of aircraft soundings. There is, however, substantial uncertainty in inferring detailed CBL structure, particularly the level of neutral stability (zn), from such a limited number of soundings. In this study, composite profiles of potential temperature are derived from multiyear early afternoon radiosonde data over Beijing, China. The CBLs become slightly stable above zn ~ 0.31–0.33zi, where zi is the CBL depth. These composite profiles are used to evaluate two K-profile PBL schemes, the Yonsei University (YSU) and Shin–Hong (SH) schemes, and to optimize the latter through parameter calibration. In one-dimensional simulations using the WRF Model, YSU simulates a stable CBL above zn ~ 0.24zi, while default SH simulates a thick superadiabatic lower CBL with zn ~ 0.45zi. Experiments with the analytic solution of a K-profile PBL model show that adjusting the countergradient flux profile leads to significant changes in the thermal structure of CBL, informing the calibration of SH. The SH scheme replaces the countergradient heat flux term in its predecessor YSU scheme with a three-layer nonlocal heating profile, with fnl specifying the peak value and specifying the height of this peak value. Increasing fnl to 1.1 lowers zn, but to too low a value, while simultaneously increasing to 0.4 leads to a more appropriate zn ~ 0.36zi. The calibrated SH scheme performs better than YSU and default SH for real CBLs.
In the presence of positive daytime surface heat flux, buoyant turbulence eddies drive the development of the convective boundary layer (CBL). The vertical heat flux in the CBL typically decreases linearly from the surface throughout the CBL and remains positive over most of the CBL, only becoming negative near the top of the CBL with absolute magnitude reaching ~0.2 of the surface heat flux (Driedonks and Tennekes 1984; Chrobok et al. 1992; Sorbjan 2009; Wang et al. 2016). As a result of such a vertical heat flux profile, the mean vertical profile of potential temperature in a typical CBL is characterized by a three-layered structure: a surface layer potential temperature decreasing with height, a near-neutral mixed layer in the middle, and a stably stratified entrainment zone at the top (Chrobok et al. 1992). While large-eddy simulations (LES) can resolve the dominant energy-containing turbulent eddies in the CBL, in operational atmospheric numerical models with much coarser resolutions, turbulent fluxes need to be parameterized using Reynolds averaging via planetary boundary layer (PBL) parameterization schemes. Since only stationary to quasi-stationary and homogeneous flows can satisfy the Reynolds averaging premise, PBL schemes are designed for such idealized conditions and can introduce errors in nonidealized conditions (Arya 2001).
In first-order closure PBL schemes, in direct analogy with molecular diffusion, turbulent transport of heat in the CBL is expressed as being proportional to the vertical gradient of potential temperature (θ):
where Kh is the eddy diffusivity of heat, w is the vertical velocity, and z is elevation. The overbar denotes the grid-resolvable mean state while the prime denotes turbulent perturbations. This is the classic gradient-transport theory, or K-theory (Arya 1999). First-order PBL schemes are widely used in operational and research numerical weather prediction (NWP), climate, and air quality models (Ayotte et al. 1996; Beare et al. 2006; Cuxart et al. 2006; Bosveld et al. 2014), including the NCEP Global Forecast System (GFS) (Han and Pan 2011; Han et al. 2016), the ECMWF Integrated Forecasting System (IFS) (Couvreux et al. 2016; Johnson et al. 2019), the Met Office (UKMO) Unified Model (Brown et al. 2008), and the National Air Quality Forecasting Capability (NAQFC) operational system (Lee et al. 2017). These first-order PBL schemes differ mainly in their approaches to determine the eddy diffusivity. Some schemes calculate eddy diffusivities using the mixing-length theory (Blackadar 1962; Louis 1979; Liu and Carroll 1996; Beare et al. 2006; Steeneveld et al. 2006; Brown et al. 2008), while others specify profiles of eddy diffusivity directly (O’Brien 1970; Troen and Mahrt 1986; Holtslag et al. 1990; Hong and Pan 1996; Lock et al. 2000; Noh et al. 2003; Hong et al. 2006; Sorbjan 2009; Kohler et al. 2011; Shin and Hong 2015).
In K-theory, an assumption of downgradient transport is made. In presence of near-zero gradient of potential temperature in the mixed layer where vertical heat fluxes are still positive, if K-theory is to produce positive fluxes, Kh needs to be infinitely large (or even negative, depending on slight variations around zero of the potential temperature gradient). Such cases violate the original assumption of downgradient transport. To address this, amendment of the K-theory was proposed to parameterize the vertical fluxes in the mixed layer with vanishing or slightly positive gradient of potential temperature (Deardorff 1966). Thus, even though near-zero gradient of potential temperature in the mixed layer is generally accepted, slight variations in the gradient (whether slightly superadiabatic or slightly statically stable) are critical in justifying this amendment to K-theory (Deardorff 1966; Brown 1996), as well as for validation (or invalidation) of other higher-order closure parameterization schemes (Shin and Hong 2011; Wang et al. 2016).
The amendment to K-theory was proposed based on evidence that in the presence of upward heat fluxes the gradient of potential temperature in the upper part of the CBL is slightly positive (i.e., the heat flux is countergradient). To simulate such situations, a countergradient term (γ) was added to the K-theory formulation for the parameterization of turbulent fluxes in CBL (Priestley and Swinbank 1947; Deardorff 1966, 1972):
where Kh and γ are always positive. Adding a certain magnitude of the countergradient term leads to warming of the upper levels in the CBL by nearly 2 K (Deardorff 1972), thus making the upper part of the simulated CBL slightly statically stable. A slightly stable upper part of CBLs has been known to be characteristic of real CBLs since at least 1950s (Stevens 2000). Note that we are referring to local static stability here. Local static stability is typically determined by the local vertical gradient of virtual potential temperature (or the gradient of potential temperature when moisture has small effect on the static stability). Considering the buoyancy of air parcels originating from all possible initial positions (e.g., from the superadiabatic layer near the surface), the entire CBL (including the upper part) is actually unstable in a nonlocal sense (Stull 1988; Arya 2001). Early evidence of a locally slightly stable upper part of the CBL was summarized by Deardorff (1966). Bunker (1956) reported five profiles of potential temperature with positive gradients between ~100–500 m measured by aircraft in the marine boundary layer over Bermuda in winter. Telford and Warner (1964) made profile measurements from aircraft on two days over New South Wales, which disclosed increasing potential temperature with height from ~140 m upward. The average profile of local Richardson number derived from the Great Plains Turbulence Field Program performed in Nebraska from 1 August to 8 September 1953 also implied a positive potential temperature gradient at 100 m and higher (Davidson and Lettau 1957). Later on, more evidence demonstrated the same phenomenon. Two aircraft soundings taken on 25 April 1968 over eastern Colorado also observed a slightly stable CBL (Lenschow 1970). Linear regression of two soundings measured off the eastern coast of Australia on 8 July 1966 also indicated a positive gradient of potential temperature in the lowest kilometer above the ground (Warner 1971). These few measured profiles, even though they provide the early justification for the countergradient amendment of the K-theory, are insufficient in number to draw solid conclusions regarding the general characteristics of CBLs, particularly in terms of detailed CBL structure [e.g., at what altitude the real CBLs start to become slightly stable; referred to as the neutral point, zn by Stevens (2000)]. To answer such questions, a large number of observed CBL profiles with high vertical resolution are needed.
The countergradient correction to the K-theory for turbulent parameterization proposed initially by Ertel (1942), Priestley and Swinbank (1947) and Deardorff (1966) was later widely used in many studies to develop/improve K-profile first-order PBL schemes; examples include Deardorff (1972, 1973), Mailhot and Benoit (1982), Troen and Mahrt (1986), Holtslag and Moeng (1991), Holtslag and Boville (1993), Holtslag et al. (1995), Frech and Mahrt (1995), Sorbjan (2009), Lock et al. (2000) for developing the PBL scheme in the UKMO Unified Model, as well as Hong and Pan (1996) when developing the Medium-Range Forecast (MRF) PBL scheme, and Noh et al. (2003) and Hong et al. (2006) in developing the Yonsei University (YSU) scheme, which has in recent years become one of most widely used PBL schemes (e.g., Hu et al. 2010a; Hu et al. 2013a; Miao et al. 2015; Sun et al. 2016; Zhu et al. 2018; Hu et al. 2019; Yang et al. 2019). The intent of the countergradient flux is to simulate transport by penetrating thermals rising from the unstable surface layer to the upper part of the CBL (Zhou et al. 2018), and it has been shown to play a role in neutralizing/stabilizing the gradients of potential temperature by cooling the lower part of CBL and warming the upper part. Limited tests show that the MRF scheme may overestimate local static stability in the upper CBL due to excessive countergradient fluxes, while the YSU scheme aims to improve the simulated CBL thermal structure (Hong et al. 2006). However, evaluation has only been performed against a limited number of observed soundings.
A new PBL scheme using the YSU treatment of local eddy fluxes (or downgradient fluxes) was developed by Shin and Hong (2015, below abbreviated as SH), in which the countergradient heat flux term was replaced with a nonlocal heat flux profile fitted to LES results, and scale awareness (or horizontal grid spacing dependency) was added to the scheme. To achieve scale awareness, the SH scheme scales both local and nonlocal eddy fluxes according to the normalized grid spacing (, where Δ is actual horizontal model grid spacing and zi is the CBL depth), with the scaling factor gradually approaching 0/1 when approaches 0.02/1. Evaluation against LES benchmarks shows that SH alleviates the tendency of YSU to overestimate the sharpness of the potential temperature profile near the top of the CBL (Shin and Hong 2015). However, the detailed thermal structure in the CBL in terms of vertical gradients of potential temperature profile, including the three-layer structure mentioned earlier, was not thoroughly evaluated, in part due to the scarcity of appropriate observations. For example, operational radiosonde launches at 0000 and 1200 UTC are in the early morning or early evening over North America and eastern Asia, while the mature stage of the CBL occurs during midafternoon. Furthermore, a limited number of aircraft soundings were used in the early days of K-profile first-order PBL scheme development, and there were uncertainties associated with such measurements (Bunker 1956; Telford and Warner 1964).
Starting in 2010, high vertical resolution L-band sounding data have been collected daily at 1300 local time during the rainy season (from June to September) at selected radiosonde sites in China (Guo et al. 2016; W. C. Zhang et al. 2018), including the Beijing site. This study uses 2010–16 L-band radiosonde data from Beijing to investigate the detailed vertical thermal structure of the CBL and to evaluate performance of the YSU and SH first-order PBL schemes within the Weather Research and Forecasting (WRF) Model (Skamarock et al. 2008; Skamarock and Klemp 2008). Viable options for optimization of the SH scheme will be identified, particularly in terms of reproducing the slightly stable upper CBL. The advantage of the scale-aware technique of the SH scheme is not discussed in this study; instead, the performance of the YSU and SH schemes applied beyond the gray zone at is the focus.
The rest of this paper is organized as follows: In section 2, the L-band radiosonde data and their processing method, an analytic solution of a K-profile CBL model, configurations of the WRF simulations used in this study, and relevant details of the SH scheme are described. In section 3, composite profiles of CBL potential temperature over Beijing are presented. The impacts of adjusting countergradient flux profiles on simulated temperature profiles are then demonstrated using the analytic solution of a K-profile CBL model, followed by WRF simulations using YSU, SH, and SH variants with adjusted flux profiles, in both single-column and three-dimensional modes. Finally, section 4 contains a summary of the main findings and some relevant discussion.
2. Data, methods, and simulation experiments
a. L-band radiosonde data and technique to produce composite profiles
In 2010, China shifted from using type 59–701 mechanical radiosondes, which operated at 403 MHz to an L-band (1675 MHz) sounding system for the 120 radiosonde sites managed by the China Meteorological Administration (CMA) (Bian et al. 2011). The L-band system uses a GTS1 digital radiosonde together with a secondary radar to retrieve atmospheric profiles. During the balloon launching process, the sounding data are collected every second, producing data with 2–6 m vertical resolution near the surface (Liu and Chen 2014). The GTS1 radiosonde temperature profiles agree quite well with those of Vaisala RS92 and RS80 radiosondes in the daytime troposphere, having a mean difference of only ~0.2 K (Bian et al. 2011). In addition to the regular daily soundings at 0000 and 1200 UTC [0800 and 2000 local time (LT), UTC = LT − 8 h], additional soundings were launched in the afternoon (~1300 LT) at selected sites during the rainy season (corresponding to the monsoon season from June to September). The Beijing sounding site (world meteorological organization (WMO) station number 54511; located at 39.8°N, 116.467°E) has the greatest number of afternoon soundings since 2010.
Beijing is located in the northern portion of the North China Plains. The Beijing sounding site is located near the southeastern corner of the 5th Ring Road of Beijing (Tian and Lu 2017; Z. Y. Zhang et al. 2018), approximately 60 km away from the Taihang and Yanshan mountains (which are situated to the west and north of Beijing; see Fig. 1). The boundary layer structure at the sounding site thus may experience mountain effects, but outside the region immediately/mostly affected by the mountains, which is estimated to be within ~50 km east of the Taihang Mountain (Hu et al. 2014; Hu et al. 2016). In the presence of southerly prevailing winds during the monsoon season, the Beijing sounding site is in the upwind region of downtown Beijing in an area with relatively flat topography.
In this study, the composite profile of CBL potential temperature over Beijing is produced based on afternoon soundings from 2010 to 2016. From the total datasets of 417 afternoon soundings available from these 7 years, 313 soundings containing CBL features (i.e., a near-neutral mixed layer with a capping inversion) are selected for use in this study. Soundings lacking such features are excluded, as those soundings are likely affected by transient atmospheric processes such as fronts, troughs, and precipitation. The 313 CBL soundings are first normalized using the CBL depth (zi), and are then averaged to obtain composite CBL profiles.
Many methods have previously been used to diagnose CBL top. In the LES community, the CBL top is usually diagnosed as the level of the minimum heat flux. Unfortunately, heat flux profiles are not available from the radiosonde observations used in this study. For sounding data, threshold Richardson number (Guo et al. 2016), and the 1.5-theta-increase method (Nielsen-Gammon et al. 2008) have proven more practical (Hu et al. 2010a; Hu et al. 2010b; Hu et al. 2013b; Miao et al. 2015; Li et al. 2017; Yang et al. 2019). The 1.5-theta-increase method defines the zi as the height where the potential temperature first exceeds the minimum potential temperature within the boundary layer by 1.5 K. An example diagnosis of the CBL top from a L-band radiosonde profile and a WRF-simulated profile is shown in Fig. 2. In this case, for the simulated profile, the 1.5-theta-increase method diagnoses a CBL top 28 m higher than that diagnosed by the YSU algorithm. The YSU scheme diagnoses the CBL top as the lowest altitude for which the bulk Richardson number between the surface and that altitude exceeds zero (simply when virtual potential temperature at that level exceeds surface virtual potential temperature plus thermal excess due to surface buoyancy flux) and thus essentially does not consider wind shear magnitude. A critical Richardson number of 0.25 is used in stable boundary layer (Hong 2010). In this study, the 1.5-theta-increase method is used to diagnose zi for normalization and generating composite profiles for both L-band radiosondes and WRF-simulated profiles.
b. Quasi-steady-state analytical solutions to a K-profile PBL model derived by Stevens (2000)
To examine the uncertainties of K-profile first-order PBL schemes in simulating CBL structures, particularly in terms of vertical temperature gradients, an analytical solution to a K-profile PBL model derived by Stevens (2000) is first examined. If turbulent fluxes are parameterized using (2), with a prescribed Kh = kz(1 − z)2 and a constant γ, the profile of potential temperature in CBL has a quasi-steady analytical solution (Stevens 2000):
where A is defined to be the ratio of the entrainment to the surface flux and k is a proportionality constant. In this study we set A = −0.2 and k = 0.675. Flux profiles using different γ will be examined to show the critical role it plays in simulating the slightly stable upper part of the CBL.
c. Flux profiles in the SH scheme
Except for adding scale-awareness treatment, the other main change to the YSU scheme made by the SH scheme is replacing the countergradient heat flux in YSU (Kh × γ) with a three-layer nonlocal heat flux profile fitted to LES results. The three-layer nonlocal vertical heat flux profile adopted by SH peaks at () and decreases linearly away from the peak value in the boundary layer, with an overlying entrainment layer with negative flux. In contrast, in the YSU scheme the countergradient heat flux profile takes a parabolic shape with maximum at ~0.3zi (Shin and Hong 2015). Given the differences between the local/nonlocal flux profiles extracted from LES simulations and parameterized gradient/countergradient flux profiles from the conventional PBL schemes (Zhou et al. 2018), replacing the countergradient heat flux profile with the LES-fitting nonlocal heat flux profile in SH can lead to model uncertainties. Motivated by the analytical solutions to a K-profile PBL model of Stevens (2000), sensitivity experiments are designed to examine the impact of different nonlocal heat flux profiles in the SH scheme. Two parameters, and fnl, in SH are tuned to adjust the nonlocal heat flux. specifies the normalized height of the surface layer, defined as the lower part of CBL where nonlocal flux increases linearly with height. fnl specifies the ratio of nonlocal heat flux to total heat flux at the top of the surface layer. Note that the surface layer defined here in SH is not necessarily the physical surface layer. Impact of uncertainties in and fnl on simulated CBL structures are examined below in one-dimensional (1D) and three-dimensional (3D) WRF simulations over Beijing and compared with the L-band radiosonde data.
d. One-dimensional (1D) simulations for a clear day (23 July 2010) over Beijing
Sensitivity simulations are first conducted using WRF v4.0 in single-column mode with the YSU and SH schemes, including runs with different values of and fnl in SH. These 1D simulations are initialized from the sounding and soil state extracted from a 3D WRF simulation (described below) at 0000 UTC (0800 LT) 23 July 2010 at the L-band radiosonde site in Beijing and run for 5 h to 1300 LT. On this day, the sky was clear, and mixed CBL characteristics developed in the early afternoon by 1300 LT (Fig. 2). Solar forcing is simulated by the 1D WRF for the specified time and location. Background wind forcing is specified using simulated wind profile extracted from the 3D WRF simulation. 105 model layers are used, with approximately 58 layers below zi to better resolve the CBL, following the approach of our previous study (Hu et al. 2011). Physics schemes include the Dudhia shortwave radiation algorithm (Dudhia 1989), the Rapid Radiative Transfer Model (RRTM) (Mlawer et al. 1997) for longwave radiation, and the Noah land surface model (Chen and Dudhia 2001).
The horizontal grid spacing is assumed to be 4 km. By 0500 UTC (1300 LT), simulated zi is lower than 2 km, thus the normalized grid spacing is larger than 2. As a result, the scale-aware scaling factor in SH for parameterized heat fluxes is approximately 1.0 [see Fig. 2a in Shin and Hong (2015)] such that the scale-aware treatment has essentially no effect. We will be focusing on modification to the nonlocal term in SH and its comparison with the YSU scheme.
e. Three-dimensional (3D) simulations for 14 days with well-developed CBLs over Beijing
To investigate the performance of the YSU and SH schemes and the impact of parameters and fnl in the SH nonlocal flux profile, 3D WRF simulations were conducted for 14 days with mixed CBL features in the summer of 2010. These days included 8, 20, 23, 24, 25, 26, and 27 July; 8, 12, 23, and 24 August; and 5, 9, and 10 September. 3D WRF simulations were initialized at 0000 UTC (0800 LT) each day and run for 30 h. Three one-way nested domains (Fig. 1) were used, with horizontal grid spacings of 27, 9, and 3 km. The CBL depth never exceeded 3 km on these days, and thus the scale-aware treatment has essentially no effect. Each domain has 48 vertical layers extending from the surface to 100 hPa. The lowest 20 model sigma levels are at 1.0, 0.997, 0.994, 0.991, 0.988, 0.985, 0.975, 0.97, 0.96, 0.95, 0.94, 0.93, 0.92, 0.91, 0.895, 0.88, 0.865, 0.85, 0.825, and 0.8 (corresponding to heights of approximately 12, 37, 61, 86, 111, 144, 186, 227, 290, 374, 459, 545, 631, 717, 826, 958, 1092, 1226, and 1409 m above ground). This configuration should be adequate to capture boundary layer structures in terms of temperature and CBL height (Hu et al. 2013c; Miao et al. 2015; Hu and Xue 2016). To the best of our knowledge, detailed gradients of potential temperature in 3D simulations of the CBL have not previously been examined thoroughly in existing literature. The 0.7° × 0.7° European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis is used for initial conditions (for both atmosphere and soil in all three domains) and boundary conditions (for the outer domain) for each selected day, following a common practice used in many studies examining PBL performance (Hu et al. 2010a; Hu et al. 2013a). The boundary conditions of the inner domains are obtained from the outer domain simulation through one-way nesting. A better soil moisture initialization may improve the PBL simulations (Angevine et al. 2014).
All model domains used the Dudhia shortwave radiation algorithm (Dudhia 1989), the Rapid Radiative Transfer Model (RRTM) (Mlawer et al. 1997) for longwave radiation, and the WRF single-moment 6-class (WSM6) microphysics scheme (Hong et al. 2004). The Grell-3 cumulus scheme (Grell and Devenyi 2002) was used on the 27 and 9 km domains, but is turned off for the 3 km domain. In initial tests of the 1D WRF, setting fnl = 1.1 in SH resulted in a slightly stable CBL starting from a very low zn, while setting and fnl = 1.1 in SH resulted in a slightly stable CBL starting from a more appropriate zn. Thus, 3D sensitivity simulations with YSU and SH schemes with default settings, and two additional SH runs with (default) and fnl = 1.1, and with and fnl = 1.1 are performed for the 14 CBL cases selected.
a. Composite profiles of potential temperature in CBL over Beijing
Composite profiles were produced for 14 CBL cases in 2010 and 313 CBL cases during 2010–16 based on the L-band radiosonde data over Beijing. The slightly stable upper part of the CBL starts at 0.32, and 0.31zi in the 14 CBL composite, and the 313 CBL composite, respectively (Fig. 3). When composites were produced using data from individual years between 2010 and 2016, the altitude of this neutral point, zn, varied between 0.28 and 0.36 (not shown); zn derived from these sounding data appears to be lower than zn derived from LESs for other cases. Zhou et al. (2018) reported a zn of 0.4zi and Stevens (2000) reported a zn ~ 0.46zi based on their LES results. Using different definitions of zi may slightly affect zn, which matters because 1.5-theta-increase method is used here, while the altitude of minimal heat flux is used in LES diagnoses. zi diagnosed using the latter method can be lower than zi diagnosed using the former by around 5%, as will be shown later. Thus, mean zn over Beijing during 2010–16 may be 0.33 if diagnosed using the minimal flux method, closer to (but still lower than) that of the LES results. Different meteorological factors, such as capping inversion strength and wind shear, may have modulated the altitude of zn for different cases. For example, when wind shear increases, buoyant turbulence-driven CBL transitions to mechanical turbulence-driven CBL, in which small eddies dominate, CBL will likely become more superadiabatic and zn increases. The actual reasons for the differences warrant future investigating. Nonetheless, the composite profile over Beijing corroborate previous limited observations (Bunker 1956; Telford and Warner 1964; Lenschow 1970), revealing that upper part of CBLs are slightly stable and thus further justifying the countergradient amendment to the original K-profile schemes. In addition, the detailed thermal structure of CBLs in terms of zn height (~0.31–0.33zi) documented by the composite profiles can be used to fine tune the K-profile PBL schemes.
b. Impact of countergradient flux profiles on temperature profile based on analytic K-profile model
Due to atmosphere’s inherent mechanism, vertical heat flux in the CBL linearly decreases from the surface throughout the CBL and becomes negative in the entrainment zone (Driedonks and Tennekes 1984; Wang et al. 2016). We first examine under such a constraint the impact of adjusting countergradient flux profile on simulated temperature profiles using the analytic solution of a K-profile PBL model derived by Stevens (2000). In this model, Kh takes a prescribed cubic form and the entrainment ratio is set as −0.2. In the quasi-steady solution characterized with linear total flux profile, different countergradient flux (Kh × γ) profiles with different γ lead to different profiles of potential temperature following the analytic solution of Eq. (3). As shown in Fig. 4, when countergradient fluxes [called nonlocal fluxes by Stevens (2000)] exceed total fluxes, the local fluxes become negative. Anywhere local fluxes are negative, the CBL becomes statically stable. The crossover point where local flux switches from positive to negative is the neutral point zn. Increasing γ leads to larger countergradient flux, thus lower crossover point.
Entrainment ratio may vary in presence of vertical wind shear (Conzemius and Fedorovich 2006a, b). Varying entrainment ratios also affect zn. In Figs. 4e–h, the analytical solution is shown with a different entrainment ratio of −0.3. Stronger entrainment leads to decreases in zn (Fig. 4), likely due to entrainment-induced turbulence penetrating deeper into the lower CBL, and also leads to a deeper slightly stable upper CBL. For an entrainment ratio of −0.3, γ affects zn the same way as for an entrainment ratio of −0.2.
The analytic solution demonstrates how the countergradient flux profile directly impacts the simulated CBL structure, and therefore provides clues for optimizing countergradient treatment in PBL schemes, including that of SH.
c. CBL structures simulated by YSU, SH, and SH variants in 1D single-column mode
As mentioned earlier, the SH scheme replaces the countergradient heat flux in YSU (Kh × γ) with a nonlocal heat flux profile fitted to LES results. and fnl, together with other parameters, specify the nonlocal CBL heat flux profile in the SH scheme. Effects of and fnl on simulated CBL structures are first examined in 1D WRF simulations and compared with YSU simulation for the Beijing site on 23 July 2010. The YSU scheme simulates a stable profile over most of the CBL with zn ~ 0.24zi, while the SH scheme with default settings ( and fnl = 0.7) simulates a thick superadiabatic layer in the lower half of the CBL with zn ~ 0.45zi (Fig. 5). Given an observed zn of ~0.31–0.33zi based on L-band radiosonde data (Fig. 3), the SH scheme may be simulating too high a zn while the YSU scheme simulates too low a zn. WRF simulations were also conducted for other CBL cases over the U.S. central Great Plains, for which the default SH scheme simulates an even higher zn of ~0.6zi (not shown). These results raise concerns that the default SH scheme may be limited in its ability to simulate the proper CBL thermal structure unless it can be tuned for improved performance.
Simulated flux profiles are examined in Fig. 6; local downgradient flux is computed as , and total flux is defined to be the sum of local and nonlocal fluxes. The default SH scheme (with and fnl = 0.7) specifies the nonlocal flux to be 0.7 times the total flux at the top of the surface layer (at ), with nonlocal flux decreasing linearly from this altitude. The nonlocal fluxes are lower than total fluxes below 0.45zi in this case, and thus the residual local flux switches from positive to negative at 0.45zi (Fig. 6a). Consequently, the model simulates zn ~ 0.45zi (Fig. 5). Note that heights in Fig. 6 are normalized by zi diagnosed using the 1.5-theta-increase method, which are slightly higher (by ~5%) than zi diagnosed from the minimal flux method. If the minimal flux method is used, zn simulated by the default SH becomes ~0.47zi. In this simulated case, the entrainment ratio is approximately ~0.3. If the entrainment ratio were smaller, the slope of the total flux would be larger and the crossover of the nonlocal flux with the total flux would be higher, which would lead to an even higher zn. In other words, if the entrainment ratio were decreased, the SH would simulate a CBL with an even thicker superadiabatic layer in the lower part of the CBL. This situation, where smaller entrainment ratios lead to higher zn, is consistent with the analytic solution (Fig. 4).
Increasing fnl to 1.1 in SH leads to larger nonlocal fluxes in the CBL, exceeding the total flux over a thicker layer in the upper part of the CBL, resulting in negative local fluxes above a shallower surface layer (Fig. 6b). As a result, a more slightly stable profile starting from a very low zn (0.12zi) is simulated. Given the observations (Fig. 3) and LES results (Stevens 2000; Zhou et al. 2018), a zn of 0.12zi is likely too low. Thus, variation of another parameter, , that specifies the normalized height of the surface layer in SH, is considered. By default is set to 0.075, while the LES of Zhou et al. (2018) suggests a deeper superadiabatic lower portion of the CBL (up to 0.4zi). Increasing from default value of 0.075 to 0.4 in SH leads to a more elevated surface layer with negative local flux starting farther aloft (Figs. 6c,d), resulting a higher zn (Fig. 5). When is set to 0.4 and fnl to 1.1, zn is increased to 0.36zi. Meanwhile, temperature in the CBL is lowered by ~0.1 K when is increased from 0.075 to 0.4 (Fig. 5).
These 1D simulations demonstrate that the actual shape of the nonlocal heat flux profile in the SH scheme can be controlled by parameters fnl and and can significantly affect the simulated CBL structures, in particular the vertical gradient of potential temperature, in a manner similar to that predicted by the analytic solution of Stevens (2000). To determine the relative performance of difference schemes for real cases, we next examine 3D simulation results.
d. CBL structures simulated by 3D WRF for 14 cases in 2010
3D WRF simulations are compared at the Beijing sounding site for the 14 cases from 2010 with L-band radiosonde data (Fig. 7), and with National Climatic Data Center (NCDC) global hourly surface data for two cases (Fig. 8). Simulated profiles of potential temperature and water vapor mixing ratios from forecasts using the default SH (with and fnl = 0.7) are shown in Fig. 7. For each 30-h simulation, initialized at 0000 UTC (0800 LT), the next four L-band soundings are available for evaluation: two early in the afternoon (1300 LT), one early in the evening (1900 LT), and one early in the morning (0700 LT). The 3D WRF simulations capture the diurnal variation of PBL structures on these days well. At the scales selected for the x axes in Fig. 7, the profiles simulated by YSU and SH nearly overlap each other and those of YSU are therefore not shown. The simulations also capture the spatial variation of surface temperature (T2) and its day-to-day variation (e.g., higher temperature on 23 July and relatively low temperature on 5 September 2010) (Fig. 8). The difference in T2 simulated by YSU and SH (<0.1 K) is again hardly discernible in Fig. 8. Modification to the SH scheme leads to change of T2 by ~0.1 K. Increasing fnl in SH to 1.1 while keeping at its default value of 0.075 reduces T2 by 0.1 K. Increasing to 0.4 in SH increases T2 by 0.08 K compared to the default SH.
Simulated and observed profiles of potential temperature for two cases, at 1300 LT on two days (23 July and 5 September 2010), are shown in Fig. 9 as zoomed-in plots to better reveal the detailed CBL structures, revealing differences in the finescale structures not discernible in Fig. 7. While all simulations produce smooth profiles, the instantaneous soundings show much greater vertical variation in CBL, likely resulting from transient turbulent processes as well as effects of horizontal nonuniform heat transport. The level of noise in the vertical variation of θ over Beijing is comparable to that of other soundings (e.g., those over Beltsville) (Hu et al. 2012; Hu et al. 2013b, figure not shown). The PBL parameterization instead tries to simulate the mean effects of boundary layer turbulence eddies in the form of an “ensemble average” (Mellor and Yamada 1974; Wyngaard and Coté 1974; Sun and Chang 1986; Nakanishi and Niino 2004). Because of the presence of transient structures, the three-layered CBL structure (i.e., superadiabatic surface layer, near-neutral mixed layer, and upper inversion layer) is harder to clearly define for the observed soundings although the presence of the stable upper CBL is clear in both soundings. The height of the neutral point in the 23 July sounding can be estimated as 0.37 km, while the CBL top can be estimated as 1.3 km (Fig. 9a), although there are significant uncertainties with such estimations. For the 5 September case (Fig. 9b), the observed structure of the lower CBL is even more complicated. Due to the presence of transient structures in the observations, using individual sounding profiles for evaluation purposes is less meaningful, and the use of coarse vertical resolution sounding data is even more problematic. Instead, using composite profiles based on many sounding profiles (such as those shown in Fig. 3) should give more robust results.
Overall, however, the differences between the YSU and SH schemes, and the two variants of the SH scheme, in terms of temperature gradients in CBL shown in Fig. 9 are consistent with the 1D WRF simulations shown in Fig. 5. YSU simulates a slightly stable boundary layer starting from a very low zn (~0.25zi on average), while the default SH gives a higher zn (~0.48zi on average) with overly weak local static stability above zn. Increasing fnl in SH to 1.1 while holding at its default value of 0.075 leads to a much stable layer that is too deep, placing zn at close to 0.1zi. Increasing to 0.4 in SH leads to a higher slightly stable CBL with a zn of 0.34zi on average. Meanwhile the temperature in CBL is lower by ~0.1 K with the increased and fnl compared to the default values of and fnl = 0.7.
To more robustly evaluate the simulated potential temperature profiles, the composite profiles for the 14 CBL cases in 2010 are shown in Fig. 10a. Here the composite profiles from the simulations are produced in the same way as for the observations, using the procedure described in section 2a. The composite of observed sounding profiles is mostly free of the transient structures seen in the individual soundings. Again, the differences between YSU and SH, and two SH variants are consistent with the 1D results. The SH scheme with fnl = 1.1 and simulates a zn of 0.35zi, which shows the best agreement with observations (zn = 0.32zi) among those configurations tested. This SH setting also best reproduces the minimum temperature in the CBL, although it slightly underestimates temperature in the upper CBL—this might be due in part to not considering the aerosol heating effect (Liu et al. 2019). Overall, the simulated composite profile with SH scheme with fnl = 1.1 and is closest to observations among the four experiments.
Given the variation of zn between 0.28 and 0.36zi in different years based on the L-band radiosonde data, we did not try to fine-tune the SH scheme further to exactly match the composite profiles from 2010. Nonetheless, our evaluation for the composite profiles in 2010 illustrates that the YSU scheme simulates too low a zn, and the default SH scheme simulates too high a zn. Adjusting the nonlocal flux profile in the SH scheme by increasing fnl, and can improve the performance of SH in terms of reproducing the detailed thermal structure in CBLs.
Following Shin and Hong (2015), our examinations of the CBL structure focus mainly on the resulting potential temperature (θ) profiles. For dry air, the θ profile determines the static stability. For moister air, virtual potential temperature (θυ), that takes into account of the effect of moisture (q) on the air density, represents the static stability more accurately. Most, if not all, PBL parameterization studies present θ and q profiles individually instead of θυ. Still, we present in Fig. 10b the comparisons of simulated and observed θυ profiles. As can be seen, the shapes of these profiles are very close to those of corresponding θ profiles shown in Fig. 10a, with the main difference being a shift of about 2.4 K toward the right due to the contribution of moisture. Given the very similar vertical structures, not surprisingly, the optimal parameters of the nonlocal heat flux term in the SH scheme would be the same based on θυ. While the observed θυ profile indicates a zn of 0.31, the default SH simulates a zn of 0.51. The two SH schemes with adjusted parameters simulate a zn closer to observation. Overall, The SH scheme with fnl = 1.1 and shows the best agreement with observations in terms of structure of the slightly stable upper CBL and superadiabatic lower CBL manifested in θυ profiles. Note that θυ includes contributions from both θ and moisture. This paper focuses on nonlocal turbulent mixing of θ while that of moisture is not considered. In fact, the treatment of nonlocal mixing of moisture or moment remains uncertain due to different physical processes that are in control (Mahrt 1976, 1991; Lock et al. 2000). Turbulent transport of moisture (and momentum) and their effects on static stability and buoyancy are topics for future studies.
4. Conclusions and discussion
Since the 1950s, a countergradient flux term has been used to amend the classic K-profile first-order PBL schemes based on observations of slightly stable conditions in the upper part of CBLs. These early observations, as documented in the literature, were mostly limited to a few aircraft soundings. Attempts to infer detailed vertical gradients of potential temperature and the neutral stability point (zn) in the CBL from those soundings contain much uncertainty. In this study, composite profiles of potential temperature are derived from multiyear, high vertical resolution L-band radiosonde data taken early in the afternoon over Beijing, China. The CBL over Beijing becomes slightly stable above zn ~ 0.31–0.33zi. These composite profiles are then used to evaluate two K-profile first-order PBL schemes, the YSU and SH schemes. Optimization of the SH scheme through parameter calibration is also proposed.
In 1D WRF simulations for a CBL case over Beijing, the YSU scheme simulates a statically stable profile over most of the CBL with zn ~ 0.24zi, while the SH scheme simulates a thick superadiabatic lower half of the CBL with zn ~ 0.45zi (in other cases, SH simulates an even higher zn of up to ~0.6zi). The uncertainties/biases of simulated zn with different PBL schemes illustrate that even the most recent SH scheme may need further calibration and tuning against a larger observational dataset.
Experiments with the analytic solution of the K-profile PBL model of Stevens (2000) show that adjusting the countergradient flux profile leads to significant changes of thermal structure in the CBL and the associated zn. Increasing countergradient flux leads to lower crossover point of local fluxes becoming negative and consequently a lower zn. These experiments offer insights valuable for calibrating the SH scheme by adjusting the countergradient flux profile used.
The SH scheme replaces the countergradient heat flux profile inherited from the YSU scheme with a three-layer nonlocal heat flux profile, with fnl specifying the peak value of the nonlocal flux and specifying the normalized height of this peak value. Increasing fnl to 1.1 in SH in 1D simulations leads to larger nonlocal fluxes in the CBL, resulting in negative local fluxes above the shallow surface layer. As a result, a more stable upper CBL profile starting from too low a zn(0.12zi) is simulated. Increasing to 0.4 in SH leads to a deeper surface layer and negative local fluxes that start at a higher level, resulting a higher zn. With and fnl to 1.1, zn is elevated to 0.36zi.
The YSU and SH schemes and two SH variants with modified and fnl values are also evaluated in 3D WRF simulations for 14 CBL cases from 2010 and compared to observed L-band radiosonde data. Based on the simulated composite potential temperature profiles, YSU produces a slightly stable boundary layer starting from too low a zn, while SH using its default settings gives too high a zn. The SH scheme with fnl = 1.1 and simulates a zn of 0.35zi, giving the best agreement with observation (zn = 0.32zi) among those configurations tested. This variant of the SH scheme also best reproduces the minimum temperature in the CBL, though it slightly underestimates temperature in the upper CBL. Based on results of this paper, we therefore recommend the use of fnl = 1.1 and for the SH scheme rather than its default values. The optimization approach used in this study can be applied to other PBL schemes as long as they have a similar nonlocal term. We evaluated the SH scheme using the sounding data at one site in this study. Its performance for other locations with different boundary layer structures warrants further investigation.
Note that in addition to modifying the nonlocal heat flux profile from its predecessor YSU, another innovation of the SH scheme is adding scale awareness to both parameterized local and nonlocal turbulence fluxes. With the scare awareness, the parameterized fluxes are scaled according to the normalized horizontal grid spacing (grid spacing normalized by CBL height). The scale-aware aspects of the SH scheme are not discussed in this study. Rather, we focus on the performance of the SH scheme and its variants at grid spacings of 3–4 km that are beyond the gray zone, with , where the scale-aware adjustment to the fluxes has little effect. When gets smaller, more turbulent fluxes can be explicitly resolved so that the parameterized fluxes are scaled back, thus models using the scale-aware SH scheme will perform increasingly more like LESs as the resolution increases. In such a case, the SH-parameterized turbulent fluxes will contribute increasingly less to total Reynolds fluxes and the bias of the default SH scheme discussed in this paper will become less relevant.
Currently the treatments for θ and moisture in YSU and SH are different. Particularly, while the countergradient mixing is considered for heat fluxes, it is not considered for moisture fluxes (Hong et al. 2006). The different treatments for θ and moisture may be due to their different mixing characteristics/processes. It is well realized that the dissipation time scale for moisture fluctuations appears longer than that for temperature fluctuations, and thus moisture is often not well-mixed even when θ is (Mahrt 1991). Evaluation and improvement of the parameterization of fluxes/profiles of moisture, as well as momentum, in mesoscale models also clearly warrants further investigation (Mellado et al. 2017).
This work was supported by the NOAA VORTEX-SE program Grant NA17OAR4590188 and NSF Grant AGS-1917701. We are grateful to Bowen Zhou for discussion. Proofreading by Nate Snook is greatly appreciated. Computations were performed at the Texas Advanced Computing Center (TACC) and earlier simulations were conducted at the San Diego Supercomputer Center (SDSC). The reanalysis dataset was downloaded from https://rda.ucar.edu/, the NCDC data used for model evaluation are downloaded from http://www7.ncdc.noaa.gov/CDO/cdo, and the radiosonde data are archived by the China Metrological Administration (CMA). Model data produced from this study have been archived at University of Oklahoma: http://www.caps.ou.edu/micronet/WRF_1D.html.