Simulated impacts of vertical distributions of black carbon aerosol on meteorology and PM_2.5 concentrations in Beijing during severe haze events

Vertical profiles of black carbon (BC) play a critical role in modifying the meteorological conditions such as temperature, planetary boundary layer height (PBLH), and regional circulation, which influence surface layer concentrations of PM_2.5 (particulate matter with a diameter of 2.5 µm or less; the surface layer covers from 0 to 79.5 m). However, BC vertical profiles in current models usually have large uncertainties. In this study, by using measurements of BC vertical profiles in Beijing collected by King Air 350 aircraft and the Weather Research and Forecasting with Chemistry model (WRF-Chem) coupled with an improved integrated process (IPR) analysis scheme, we investigated the direct radiative effect (DRE) of BC with different vertical profiles on meteorology and PM_2.5 concentrations in Beijing during two severe haze events (11–12 and 16–19 December 2016). Compared with measurements in Beijing, the model overestimated BC concentrations by 87.4 % at the surface and underestimated BC mass by 14.9 % at altitudes of 300–900 m as averaged over the two pollution events. The BC DRE with the default vertical profiles from the model heated the air around 300 m altitude, but the warming would be stronger when BC vertical profiles were modified for each day using the observed data during the two severe haze events. Accordingly, compared to the simulation with the default vertical profiles of BC, PBLH was reduced further by 24.7 m (6.7 %) and 6.4 m (3.8 %) in Beijing in the first and second haze events, respectively, with the modified vertical profiles, and hence the surface layer PM_2.5 concentrations were higher by 9.3 µg m⁻³ (4.1 %) and 5.5 µg m⁻³ (3.0 %) over central Beijing, owing to increased positive contributions of vertical mixing and chemical processes. Furthermore, we quantified by sensitivity experiments the roles of BC vertical profiles with six exponential decline functions ($C\left(h\right)=C_{\mathrm{0}}\times e^{-h/\mathrm{hs}}$ and hs​​​​​​​ = 0.35, 0.48, 0.53, 0.79, 0.82, and 0.96) parameterized on the basis of the observations. A larger hs means less BC at the surface and more BC in the upper atmosphere, resulting in less solar radiation reaching the ground and consequently a stronger cooling at the surface (+0.21 with hs of 0.35 vs. −0.13^∘ with hs of 0.96). Our results indicate that it is very important to have accurate vertical profiles of BC in simulations of meteorology and PM_2.5 concentrations during haze events.

1 Introduction

With the rapid economic development and large increases in fossil energy consumption, haze pollution has become one of the most serious challenges in China, especially in the Beijing–Tianjin–Hebei (BTH) region (H. Wang et al., 2015; Zhang et al., 2019). In 2014 and 2015, the numbers of extremely serious PM_2.5 (particulate matter with an aerodynamic equivalent diameter of 2.5 µm or less) pollution days (with daily mean PM_2.5 > 150 µg m⁻³) in Beijing reached 45 and 54, respectively (He et al., 2017). The real-time hourly average concentration of PM_2.5 in Beijing even reached 1000 µg m⁻³ during the severe haze events in January 2013, far exceeding the Chinese Ambient Air Quality Grade I Standards (35 µg m⁻³ for daily mean PM_2.5) (Liu et al., 2017). With the implementation of the toughest-ever clean air policy since 2013, the observed annual mean PM_2.5 concentrations averaged over 74 cities in China fell from 61.8 µg m⁻³ in 2013 to 42.0 µg m⁻³ in 2017 (Zhang et al., 2016; Wang et al., 2017; K. Li et al., 2019; Zhang et al., 2019). However, severe haze events still occurred in Beijing during the COVID-19 lockdown period (January–February 2020) (Huang et al., 2020; Zhu et al., 2020). Therefore, understanding the mechanisms responsible for the occurrence of severe haze is important for air quality management planning.

BC, an important component of PM_2.5, is emitted mainly from the incomplete combustion of fossil fuel, biofuel, and biomass burning. BC particles can strongly absorb solar radiation in the atmosphere, which alters the Earth's radiation balance (Bond et al., 2013; Huang et al., 2015; Hu et al., 2020). In recent years, researchers have found that the radiative effect of BC significantly affects the structure of the planetary boundary layer (PBL) during severe haze pollution events (Ding et al., 2016; Huang et al., 2018; Wang et al., 2018; Liu et al., 2019). By using the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem), Ding et al. (2016) illustrated that BC suppressed the development of PBL by heating the air in the upper PBL and reducing the solar radiation at the surface in Beijing in December 2013. This process was defined as the “dome effect” of BC by Ding et al. (2016). This dome effect was also found over the Indian Ocean (Wilcox et al., 2016). BC can also change the land–sea thermal contrast and induce circulation anomalies during severe haze events (Gao et al., 2016b; Qiu et al., 2017; Q. Ding et al., 2019; Chen et al., 2021). By using the WRF-Chem model, Q. Ding et al. (2019) showed that, during a haze event in December 2013, the direct radiative effect (DRE) of BC enhanced advection between land and sea by causing a cooling (−1.0^∘) in air temperature over land and a warming (+0.8^∘) in air temperature over sea, which transported moist air from the sea to the Yangtze River delta region. By using the WRF-Chem model, Qiu et al. (2017) and Chen et al. (2021) also reported that the radiative effect of BC induced strong anomalous northeasterly winds from the sea during a haze event in the North China Plain (NCP) in February 2014.

BC can influence concentrations of PM_2.5 during haze events because of its impact on PBL and other meteorological fields (Gao et al., 2016b; Wilcox et al., 2016; Miao et al., 2017; Qiu et al., 2017; Gao et al., 2018; Wang et al., 2019; Chen et al., 2021). Gao et al. (2016b) used the WRF-Chem model to simulate the haze event that occurred in the NCP in January 2010 and found a maximum increase in PM_2.5 of 14.4 µg m⁻³ (5.1 %) due to the DRE of BC. Qiu et al. (2017) also analysed the impact of BC on surface layer PM_2.5 during haze pollution in NCP in February of 2014 by using the WRF-Chem model and found that the average PM_2.5 concentration increased by 2.1 µg m⁻³ (1.0 %) owing to the DRE of BC. Chen et al. (2021) analysed, by using the WRF-Chem model, the DRE from the ageing of BC and its impact on PM_2.5 concentration over the BTH region during a haze event in February 2014. They found that the near-surface PM_2.5 concentration average over BTH increased by 9.6 µg m⁻³ (7.0 %) due to the ageing of BC.

So far few studies examined the impacts of vertical distributions of BC aerosol on meteorology and PM_2.5 concentrations. Wang et al. (2018) examined the role of BC at different altitudes in influencing PBL height (PBLH) by considering a single column using WRF-Chem version 3.6.1. They divided the height from 150 to 2250 m evenly into seven layers and increased BC concentrations from 0 to 30 µg m⁻³ with an increment of 2 µg m⁻³ at one of the layers, with the BC concentrations at the other layers fixed to 0 µg m⁻³. Model results showed that the near-surface BC could increase PBLH by 0 %–4 %, while BC aloft would decrease PBLH by 2 %–16 % due to the warming of the atmosphere by BC. Current chemistry–climate models were reported not to be able to represent the BC vertical profiles accurately, so sensitivity studies were carried out to adjust vertical profiles of BC in the model by changing the vertical resolution, aerosol microphysical scheme, and emission height (Wang et al., 2019; Yang et al., 2019; Watson-Parris et al., 2019).

In recent years, measurements of BC vertical distribution have been conducted by aircraft during the severe haze events in Beijing, using a single-particle soot photometer (SP2) (Zhao et al., 2018; Tian et al., 2019; Zhao et al., 2019; Tian et al., 2020; Liu et al., 2020). During the period of severe pollution from 11 to 19 December 2016, Zhao et al. (2019) collected BC vertical profiles over Beijing by Air350 aircraft and found that the vertical profiles can be classified into two types. The first type was characterized by decreases in BC concentration with altitude, which was the case mainly controlled by local emissions. The second type had maximum BC concentration around 900 hPa, which was mainly affected by regional transport from the polluted south-southwest region. Generally, the first type occurred more frequently than the second type during haze events in Beijing. Observations of vertical profiles of BC in severe haze events over Beijing in 2018 by a King Air 350 aircraft by S. Ding et al. (2019) also obtained the same types of profiles.

In this work, we use the BC vertical profiles observed during two severe haze events (11–12 and 16–19 December 2016) over Beijing and the online-coupled WRF-Chem model to investigate the DRE of BC vertical profiles on meteorology and PM_2.5 concentrations. Compared with previous studies that examined the impact of BC on meteorology and PM_2.5, our study is the first to pay attention to the role of BC vertical profile as well as the underlying mechanism. The description of the model, observational datasets, and numerical experiments are presented in Sect. 2. Section 3 evaluates simulated meteorological and chemical variables by comparing with observations. Section 4 compares the DRE of BC with the original and modified vertical profiles, and Sect. 5 discusses the role of BC vertical profiles in influencing meteorological parameters and PM_2.5 concentrations. The conclusions of this study are given in Sect. 6.

2 Method

2.1 Model configuration

A fully coupled online Weather Research and Forecasting with Chemistry model (WRF-Chem version 3.7.1) (Grell et al., 2005; Chapman et al., 2009) was employed to simulate the two severe haze events in Beijing from 7 to 20 December 2016, and the initial 4 d are spin-up. This model adopts Lambert projection and two nested domains with grid resolutions of 30 km (domain 01) and 10 km (domain 02). Figure 1 shows that the outer domain covers most of China with 100 (west–east) × 100 (south–north) grid cells, and the second domain covers the BTH region with 58 (west–east) × 76 (south–north) grid cells. The number of vertical layers is 29, with the first 15 layers below 2 km for finer resolution in the PBL. Meteorological initial and boundary conditions in this model were derived from National Centers for Environmental Prediction (NCEP) Final (FNL) Operational Model Global Tropospheric Analyses (ds083.2) with a spatial resolution of 1^∘ × 1^∘. MOZART-4 (Model for Ozone And Related chemical Tracers-4) simulation results provided the initial and lateral boundary conditions for the concentrations of chemical species in our model (Emmons et al., 2010).

Figure 1 (a) Two nested domains with grid resolutions of 30 km (d01) and 10 km (d02). (b) The BC vertical profiles were modified for the blue box, which covers all of Beijing.

Anthropogenic emission data in the year 2016 were obtained from the MEIC inventory with a spatial resolution of 0.25^∘ × 0.25^∘ (Zheng et al., 2018). This inventory includes sulfur dioxide (SO₂), nitrogen oxides (NO_x), carbon monoxide (CO), non-methane volatile organic compounds (NMVOCs), ammonia (NH₃), BC, organic carbon (OC), PM_2.5, PM₁₀, and carbon dioxide (CO₂), which were categorized into agriculture, industry, residence, transport, and power generation sectors (Li et al., 2017). The biogenic emissions were calculated online using MEGAN (Model of Emissions of Gases and Aerosol from Nature), including isoprene, terpene, and other substances emitted by plants (Guenther et al., 2006). Biomass burning emissions were taken from the Fire INventory from NCAR (FINNv1.5), which provides daily emissions at a horizontal resolution of ∼ 1 km² (Wiedinmyer et al., 2011).

Table 1 Physical and chemical options for WRF-Chem.

The parameterization schemes of physical and chemical processes of the WRF-Chem model adopted in the study are summarized in Table 1. The Carbon-Bond Mechanism version Z (CBMZ) is chosen to simulate the gas-phase chemistry. The aerosol scheme is the Model for Simulating Aerosol Interactions and Chemistry (MOSAIC), which includes sulfate, nitrate, ammonium, chloride, sodium, BC, OC, and other inorganic aerosol, and the aerosol particles are divided into eight particle size segments. However, the formation of secondary organic aerosol is not considered in this scheme (Zhang et al., 2012; Gao et al., 2016a). In MOSAIC, the aerosol particles are assumed to be an internal mixture, and aerosol optical properties are calculated by the volume averaging mixing method (Barnard et al., 2010; Stelson, 1990). The choice for photolysis schemes is the Fast-J photolysis scheme.

2.2 Integrated process rate (IPR) analysis

The IPR analysis has been widely applied to illustrate the impacts of each physical and chemical process on the variations in O₃ concentrations (Zhang and Rao, 1999; Jiang et al., 2012; Gao et al., 2017, 2018). The improved IPR analysis method developed by Chen et al. (2019) in the WRF-Chem model is used in this work to quantitatively analyse the contributions of physical and chemical processes to PM_2.5 concentrations, including the contributions from the sub-grid convection (CONV), vertical mixing (VMIX), chemistry (CHEM), regional transport (TRA), wet scavenging (WET), emission source (EMI), and other processes (OTHER). CONV refers to the transport within the sub-grid wet convective updrafts, downdrafts, and precipitation (Chen et al., 2019), and VMIX is affected by atmospheric turbulence and vertical distribution of PM_2.5 concentrations (Zhang and Rao, 1999; Gao et al., 2018). CHEM represents PM_2.5 production and loss including gas-phase, cloud, and aerosol chemistry. TRA is caused by advection, which is highly related to wind and horizontal distribution of PM_2.5 concentrations (Gao et al., 2018; Chen et al., 2019). WET represents the wet removal processes of aerosols by in-cloud scavenging and below-cloud washout. EMI is controlled by emission source. OTHER represents the processes other than the above six processes in the model. The NET is the sum of all physical and chemical processes, which matches the variations in PM_2.5 concentrations. It is worth noting that each IPR variable is an accumulated value which is the sum of each time step.

2.3 Observational data

To evaluate the model performance in simulating near-surface meteorological fields, the observed hourly temperature at 2 m (T2​​​​​​​), relative humidity at 2 m (RH2), wind speed at 10 m (WS10), and wind direction at 10 m (WD10) at Beijing Capital International Airport station (40.08^∘ N, 116.58^∘ E) are collected from NOAA's National Climatic Data Center (http://gis.ncdc.noaa.gov/maps/ncei/cdo/hourly, last access: 3 February 2022). Due to the limited observations of planetary boundary layer heights (PBLHs), shortwave downward radiation flux (SWDOWN), and total cloud cover in Beijing, the reanalysis data of 3-hourly PBLH, SWDOWN, and total cloud cover for Beijing from the Global Data Assimilation System (GDAS) with a spatial resolution of 1^∘ × 1^∘ (http://ready.arl.noaa.gov/READYamet.php, last access: 3 February 2022) were used for model evaluation. More details about the GDAS dataset can be found in Rolph (2013) and Kong et al. (2015). The radiosonde data (temperature and relative humidity profiles) in Beijing were obtained from the University of Wyoming, Department of Atmospheric Science (http://weather.uwyo.edu/upperair/sounding.html, last access: 3 February 2022). Hourly concentrations of PM_2.5, CO, NO₂, SO₂, and O₃ at Beijing station were obtained from the China National Environmental Monitoring Center (CNEMC, http://www.cnemc.cn/, last access: 3 February 2022), which were used to evaluate the model performance in simulating pollutants at the surface. Aerosol optical depth (AOD) at 550 nm over China retrieved from the MODIS (Moderate Resolution Imaging Spectroradiometer) satellite was used to evaluate the horizontal distribution of simulated optical properties of aerosols in this study (https://ladsweb.modaps.eosdis.nasa.gov/, last access: 3 February 2022). The MYD03 (Level-1A) product with 1 km spatial resolution from the Aqua platform and the MOD03 (Level-1A) product with 1 km spatial resolution from the Terra platform were used in this study. The values of daily aerosol optical depth (AOD) at 500 and 675 nm in Beijing were obtained from the AERONET dataset (https://aeronet.gsfc.nasa.gov/, last access: 3 February 2022).

The vertical profiles of BC mass concentrations in Beijing were collected by King Air 350 aircraft using SP2 from 11–12 and 16–19 December 2016. The aircraft departed from Shahe (∼ 20 km to central Beijing) (Fig. 1) at 12:00–13:00 local time (LT) and returned around 15:00–16:00 LT, which avoided the possible diurnal variation in the PBL among flights. Most flights could reach 2.5 km. Zhao et al. (2019) reported that these vertical profiles of BC could be expressed as an exponential decline function $C\left(h\right)=C_{\mathrm{0}}{e}^{-\frac{h}{hs}}$ except 11 December 2016, where C(h) (µg m⁻³) is BC concentration at altitude h (km), C₀ (µg m⁻³) is BC concentration at the surface, and each hs value is calculated for each flight of BC vertical profile using nonlinear regression (Table S1 in the Supplement). Tian et al. (2019) observed a regional transport of pollution in Beijing from 10 to 12 December 2016 using SP2, and they found a different vertical structure of BC from that of Zhao et al. (2019), with the BC concentration at the altitudes of 400–900 m being 1.5 times higher than the near-surface BC concentration on 11 December 2016. More detailed information about the King Air 350 aircraft dataset can be found in Zhao et al. (2019), S. Ding et al. (2019), and Tian et al. (2019).

2.4 Numerical experiments

To compare the DRE of BC with the original and corrected vertical profiles and quantify the role of BC vertical profiles in influencing meteorological conditions and air pollutants, we performed the following numerical experiments as summarized in Table 2.

1. CTRL. This is the control simulation with the direct and indirect radiative effects of all aerosols (BC, OC, sulfate, nitrate, ammonium, Na⁺, Cl⁻, and other inorganics – OIN) included for the time period of 11–20 December 2016. The vertical profiles of BC were the default ones simulated by the model.

2. NoBCrad. This is the same as the CTRL simulation, except that the DRE of BC was turned off.

3. VerBC_obs. This is the same as the CTRL simulation, except that the BC vertical profiles in the model were modified according to the observed ones. The specific method will be discussed below.

4. VerBC_hs1-6. This is the same as the CTRL simulation, except that the vertical profiles of BC in the model were modified according to the exponential decline function ($C\left(h\right)=C_{\mathrm{0}}\times e^{-h/\mathrm{hs}}$). The values of hs in VerBC_hs1 to VerBC_hs6 were 0.35 to 0.96 (from small to large), respectively.

5. VerBC_RT. This is the same as the VerBC_hs1-6 simulations, except that the BC vertical profiles in the model were modified according to the observed transport BC vertical profile on 11 December 2016, which was affected by regional transport.

Table 2 Numerical experiments. Y indicates “on”, and N indicates “off”.

^∗ The values of hs in VerBC_hs1, VerBC_hs2, VerBC_hs3, VerBC_hs4, VerBC_hs5, and VerBC_hs6 are 0.35, 0.48, 0.53, 0.79, 0.82, and 0.96, respectively.

In the case of NoBCrad, the BC DRE was turned off by setting the BC mass concentration equal to zero when calculating the optical properties of BC, following the studies of Qiu et al. (2017) and Chen et al. (2021). In the VerBC_obs experiment, we modified the simulated BC vertical profile online using the observed BC vertical profile on the corresponding day. Firstly, we interpolated the observed BC concentrations to the height of each layer in the model as C_{obs_int(i)}. Each layer in the model has a top height and a bottom height, and we selected the middle height of this layer for interpolation. Secondly, we used C_{obs_int(i)} to calculate the BC mass column burden in each layer (M_{obs_int(i)}) in the model, and P_{obs_int(i)} is the percentage of BC mass column burden in each layer to the total BC mass column burden (Fig. S1) calculated by

$$\begin{array}{}\text{(1)}& M_{\mathrm{obs}\mathrm{\_}\mathrm{int}\left(i\right)}=C_{\mathrm{obs}\mathrm{\_}\mathrm{int}\left(i\right)}\cdot \left(H_{\mathrm{sim}\mathrm{\_}\mathrm{top}\left(i\right)}-H_{\mathrm{sim}\mathrm{\_}\mathrm{bot}\left(i\right)}\right),\end{array}$$

where H_{sim_top(i)} is the top height of layer i and H_{sim_bot(i)}is the bottom height of layer i.

$$\begin{array}{}\text{(2)}& P_{\mathrm{obs}\mathrm{\_}\mathrm{int}\left(i\right)}={\displaystyle \frac{M_{\mathrm{obs}\mathrm{\_}\mathrm{int}\left(i\right)}}{\sum _{i=\mathrm{1}}^{\mathrm{16}}{M}_{\mathrm{obs}\mathrm{\_}\mathrm{int}\left(i\right)}}}\cdot \mathrm{100}\mathit{\%}\end{array}$$

In the VerBC_obs simulation, the simulated BC mass column burden was redistributed to each layer below 2.5 km according to the calculated P_{obs_int(i)}. These procedures ensure that the modification of the BC vertical profile for each day in the model by using the observed data does not change the total BC mass column burden. Since the aircraft measured only BC concentrations below 2.5 km, we modify the BC profile up to the 16th model layer (about 2.5 km in Beijing).

In the experiment of VerBC_hs1∼6, we also used the above method to modify the BC vertical profile by an exponential decline function, which is $C\left(h\right)=C_{\mathrm{0}}\times e^{-h/\mathrm{hs}}$. However, in cases of VerBC_hs1∼6, we modified for the dates of 12 and 16–19 December. On 11 December, BC did not show an exponential decline with height due to the regional transport. In the simulation of VerBC_RT, the method and setting were the same as VerBC_hs1∼6, except that the BC vertical profile in the model was modified according to the observed one on 11 December 2016. In the VerBC_obs, VerBC_hs1∼6, and VerBC_RT cases, the modifications of BC vertical profiles were performed only when the direct radiative forcing of BC was calculated. All other physical and chemical processes in these experiments still used the original BC vertical profiles simulated by the model. The BC vertical profiles were only modified in the blue square shown in Fig. 1a.

3 Model evaluation

3.1 Near-surface air pollutants and BC vertical profiles

Figure 2a–i show the horizontal distributions of simulated surface layer PM_2.5 concentrations at 14:00 LT from 11 to 19 December 2016. In Beijing, high PM_2.5 concentrations of 234.1 and 165.9 µg m⁻³ occurred on 11 and 12 December, respectively. The severe pollution in Beijing on 11 December was caused by regional transport from the heavily polluted southern area under a prevailing southerly air flow (Tian et al., 2019). From 16 December, PM_2.5 started to accumulate in eastern China, and the concentrations of PM_2.5 reached the highest value of 217.8 µg m⁻³ on 19 December in Beijing. The daily PM_2.5 concentrations (Fig. 2j) in Beijing had low values from 13–15 December 2016. The model results for Beijing in this paper are the averages over the region of the blue square shown in Fig. 1a unless stated otherwise. The severe pollution from 16–19 December 2016 was mainly affected by local emissions. We are mainly focused on the two heavy pollution incidents (11–12 and 16–19 December 2016) in the following sections.

Figure 2 (a–i) Simulated spatial distributions of surface layer PM_2.5 concentrations (µg m⁻³) and winds (m s⁻¹) at 850 hPa at 14:00 LT from 11 to 19 December 2016. Black and blue squares in each panel denote the regions of Beijing–Tianjin–Hebei and Beijing, respectively. (j) Time series of simulated mean daily PM_2.5 concentration from 11 to 19 December 2016 averaged over Beijing (blue square) and BTH (black square).

Figure 3 Time series of the observed (black dots) and simulated (red lines) hourly concentrations of PM_2.5 (µg m⁻³), NO₂ (ppbv), O₃ (ppbv), CO (ppmv), and SO₂ (ppbv) in Beijing from 11 to 19 December 2016. The observations and simulations in Beijing were averaged over 12 observational sites and corresponding grid points, respectively.

Table 3 Statistical metrics for temperature at 2 m (T2; ^∘), relative humidity at 2 m (RH2; %), wind speed at 10 m (WS10; m s⁻¹), wind direction at 10 m (WD10, ^∘), PBLH (m), SWDOWN (W m⁻²), PM_2.5 (µg m⁻³), SO₂ (ppbv), NO₂ (ppbv), CO (ppmv), and O₃ (ppbv).

^a,b SIM and OBS represent the averaged model results and observations in Beijing from 11 to 19 December 2016.
^c R is the correlation coefficient, which is calculated between hourly observations and simulations in Beijing from 11 to 19 December 2016, $R=\frac{{\sum }_{i=\mathrm{1}}^n\mathrm{|}({\mathrm{OBS}}_i-\mathrm{OBS})\cdot ({\mathrm{SIM}}_i-\mathrm{SIM})\mathrm{|}}{\sqrt{\sum _{i=\mathrm{1}}^n({\mathrm{OBS}}_i-\mathrm{OBS}{)}^{\mathrm{2}}+{\sum }_{i=\mathrm{1}}^n({\mathrm{SIM}}_i-\mathrm{SIM}{)}^{\mathrm{2}}}}$, where OBS_i and SIM_i are the hourly observed and simulated data in Beijing and n is the total number of hours.
^d MB is the mean bias, $\mathrm{MB}=\frac{\mathrm{1}}{n}\cdot {\sum }_{i=\mathrm{1}}^n{\mathrm{SIM}}_{i-{\mathrm{OBS}}_i}$.
^e NMB is the normalized mean bias, $\mathrm{NMB}=\frac{\mathrm{1}}{n}\cdot {\sum }_{i=\mathrm{1}}^n\frac{{\mathrm{SIM}}_i-{\mathrm{OBS}}_i}{{\mathrm{OBS}}_i}\cdot \mathrm{100}\mathit{\%}$.
^f MFB is the mean fraction bias, $\mathrm{MFB}=\frac{\mathrm{2}}{n}\cdot {\sum }_{i=\mathrm{1}}^n\frac{{\mathrm{SIM}}_i-{\mathrm{OBS}}_i}{{\mathrm{SIM}}_i+{\mathrm{OBS}}_i}\cdot \mathrm{100}\mathit{\%}$.

Results from the CTRL simulation were compared with the observed hourly surface concentrations of PM_2.5, NO₂, O₃, CO, and SO₂ from 11–19 December 2016 in Beijing in Fig. 3. The observed maximum PM_2.5 concentration of 219.5 µg m⁻³ occurred on 18 December, far exceeding the national air quality standard for daily PM_2.5 of 75 µg m⁻³ (Wang et al., 2017). The correlation coefficient (R), mean bias (MB), normalized mean bias (NMB), and mean fraction bias (MFB) are summarized in Table 3. The model can reasonably reproduce the temporal variations in PM_2.5, NO₂, O₃, and CO; the correlation coefficients between simulated and observed hourly concentrations are 0.77, 0.78, 0.66, and 0.73, respectively. The correlation coefficient for SO₂ is lower (0.38). Gao et al. (2016b) explained that the WRF-Chem model cannot represent the SO₂ concentration and its change with time well due to the uncertainty in SO₂ emissions and missing heterogeneous oxidation. Compared with observations, the model overestimates the concentrations of PM_2.5 and NO₂ in Beijing with the MBs and NMBs of 13.2 µg m⁻³ and 10.0 % and 8.5 ppbv and 21.6 %, respectively, and underestimates the concentrations of O₃ (−0.1 ppbv, −1.2 %) and CO (−0.1 ppmv, −4.9 %). It should be noted that the model performance in simulating PM_2.5, SO₂, CO, and O₃ is better during the two haze events than on clean days. For hourly PM_2.5, for example, the MBs (NMBs) are 29.1 µg m⁻³ (82.5 %) on clean days and 6.3 µg m⁻³ (3.5 %) during the two haze events. The possible reasons for the overall overestimation of PM_2.5 are as follows: (1) the model biases in underestimating WS10 and daytime PBLH and (2) the uncertainties in anthropogenic emission data (e.g. the overestimation in the BC emissions) (Qiu et al., 2017; Chen et al., 2021). Overall, the model can capture the two severe pollution events in Beijing from 11–19 December 2016 fairly well.

Figure 4 Observed (black line), simulated (red line), and modified (blue line) BC vertical profiles in Beijing on 11–12 and 16–19 December 2016. The time of observation is indicated on top of each panel. The model results are 2 h averages around the observation time.

Because of the lack of measured BC vertical profiles from 13–15 December 2016 in Beijing, Fig. 4 compares only the simulated vertical profiles of BC with observations for the two polluted events (11–12 and 16–19 December 2016). Observed mass concentrations of BC decreased exponentially with altitude on all days except for 11 December when regional transport of pollution dominated. On 11 December, the observed maximum mass concentration of BC (7.0 µg m⁻³) occurred at 850 m altitude, which was much higher than the surface layer concentration of 4.7 µg m⁻³. Compared with the observed vertical profiles of BC, the model can represent the decreases in BC mass concentration with height on 12 and 16–19 December well but cannot reproduce the vertical profile on 11 December. Possible reasons for the model's failing to represent the BC vertical profile on 11 December are as follows: (1) the model cannot capture the wind at high altitudes and does not reproduce the high-altitude BC concentrations in the surrounding areas of Beijing, and (2) the model underestimates the daily maximum PBLH in Beijing, which inhibits the upward transport of surface layer BC. Averaged over the two pollution events, the simulated BC mass concentration was overestimated by 87.4 % on the ground and underestimated by 33.1 % at an altitude of 1000 m compared with the observations in Beijing. The inaccuracy of the vertical distribution of BC would lead to inaccurate representation of the interactions between BC and PBL, especially in heavily polluted events.

3.2 Meteorological parameters

The first haze event started on 11 December when southeasterlies transported polluted air from southern BTH to Beijing (Fig. 2a). Although the southeasterlies turned into northeasterlies in Beijing on 12 December, PM_2.5 concentrations were still high because of the high relative humidity (63.2 %) that was conductive to the formation of secondary aerosols. With the relatively high wind speed of 3.6 m s⁻¹ and low relative humidity of 37.2 % in Beijing from 13–15 December, the haze pollution gradually disappeared (Fig. 2c–e). From 16 to 19 December, PM_2.5 began to accumulate again with unfavourable diffusion conditions (WS10 of 1.4 m s⁻¹) and enhanced formation of secondary aerosols under high relative humidity of 67.1 % (J. Li et al., 2019; Dai et al., 2021). Throughout the simulated period of 11–19 December 2016, Beijing had no precipitation and was partly cloudy (Fig. 5g–h). Figure 5 shows the hourly simulated and observed T2, RH2, WS10, WD10, PBLH, SWDOWN, precipitation, and total cloud cover in Beijing from 11 to 19 December 2016. The statistical metrics are summarized in Table 3. In the two severe haze events, the observed maximum RH2 on each day exceeded 70.0 %, which accelerated the formation of secondary aerosols (Sun et al., 2006; Wang et al., 2014). Compared with the observations, the model can represent the temporal variation in T2 and RH2 well with correlation coefficients of 0.77 and 0.75, respectively, but slightly overestimates T2 with a MB of 0.1^∘ and underestimates RH2 with a MB of 3.4 %. For WS10, observations and simulated results both show low wind speed with the mean values of 1.5 and 1.4 m s⁻¹ in Beijing during the two periods of haze events. Such a meteorological condition was very beneficial to the accumulation of near-surface pollution. The WRF-Chem model also captures the high values of WS10 from 14 to 15 December. For wind direction at 10 m, the NMB is −9.0 % and the R is 0.45, which indicates that the model can simulate the change of wind direction during the period of heavy pollution. For PBL, the reanalysis PBL was 118.7 m during the two severe haze events, compared to 287.5 m during the clean period. The model can represent the change of PBLH in Beijing from 11 to 19 December 2016 with R of 0.72. The model overestimates PBLH by 30.9 m (17.7 %) in Beijing as averaged over 11–19 December 2016. The overestimation is mainly in the hours of 00:00–08:00 and 17:00–23:00 LT. It is noted that the reanalysis PBLH values provided by GDAS of NOAA were mostly 50–60 m at 00:00–08:00 and 17:00–23:00 LT in Beijing, far below the simulated mean value of 154.5 m in these hours. There might be biases in the PBLH from GDAS. Several previous studies showed that the values of observed PBLH from lidar measurements were about 200 m at night during haze events (Wang et al., 2012; Luan et al., 2018; Chu et al., 2019). The simulated SWDOWN in CTRL experiment agrees well with the observations with R and MB of 0.76 and −14.9 W m⁻². Due to the limitation of the model outputs, the model provides only information of whether there is cloud in the grid or not. The model can reproduce the presence of cloud from 11–19 December 2016 well. Both observations and model results show no precipitation in the studied time period.

Figure 5 Comparisons of simulated meteorological parameters from CTRL simulation with measurements. The black dots and red lines are the observed (reanalysis) and simulated hourly data of T2 (^∘), RH2 (%), precipitation (mm), and 3-hourly data of PBL height (m), SWDOWN (W m⁻²), total cloud cover (%), 6-hourly data of WS10 (m s⁻¹), and daily data of WD10 (^∘) in Beijing from 11 to 19 December 2016. PBLH, SWDOWN, and total cloud cover are taken from GDAS. The WRF-Chem model output shows only whether a grid has cloud (Y) or no cloud (N).

The simulated and observed vertical profiles of temperature in Beijing from 11–19 December 2016 are shown in Fig. S2. The observed vertical temperature profiles are available only at 08:00 and 20:00 LT. During the two severe pollution events, strong temperature inversions below 1500 m were observed in Beijing, which inhibited vertical mixing and caused the accumulation of pollutants near the ground. The model captures these temperature inversions well but overestimates the inversion layer height on 11 December and underestimates the inversion layer height from 16 to 19 December. The inaccuracy of the simulated inversion layer height may be due to the fact that the model cannot correctly represent the vertical profiles of BC (Fig. 4).

3.3 AOD and AAOD

AOD (AAOD) is the measure of aerosols (absorbing aerosols) distributed within a column of air from the surface to the top of the atmosphere (Khor et al., 2014). Figure S3 shows the horizontal spatial distributions of observed and simulated AOD at 550 nm over the NCP averaged over 11–19 December 2016. The model can generally reproduce the horizontal distribution of AOD, with a spatial correlation coefficient of 0.89. However, the model underestimates AOD over the NCP region. Many previous studies have shown that MODIS retrieval tends to overestimate AOD over NCP (Li et al., 2016; Qiu et al., 2017). We also compared the simulated hourly AAOD at 550 nm with AERONET AAOD at Beijing and Xianghe stations in Fig. 6. The correlation coefficient between simulations and observations is 0.85. Compared with AERONET AAOD, simulated AAOD values at Beijing and Xianghe are overestimated by 0.02 (33.3 %) and 0.02 (39.9 %), respectively.

Figure 6 Comparison of simulated absorption aerosol optical depth (AAOD) at 550 nm with observations at Beijing (39.98^∘ N, 116.38^∘ E) and Xianghe (39.75^∘ N, 116.96^∘ E) stations from 11 to 19 December 2016.

4 A comparison of BC DRE with the original and modified vertical profiles

As shown in Fig. 4, the model does not represent the vertical distribution of BC concentrations well during the two heavily polluted events, especially on 11 December. So, in this section, we examine the differences in the BC DRE on meteorology and concentrations of pollutants with the original and modified vertical profiles.

4.1 Direct radiative effect of BC on meteorology

Figure 7 shows the atmospheric temperature and PBLH simulated from the CTRL simulation and their changes caused by BC DRE with the original profiles (CTRL minus NoBCrad) and modified profiles (VerBC_obs minus NoBCrad), over Beijing from 11 to 19 December 2016. Light-absorbing BC heated the air at around 300 m on 11 and 16–19 December, regardless of the original or modified BC vertical profiles (Fig. 7b and c). With the original and modified BC profiles, the maximum warming effects in the PBL were 0.8 and 0.9^∘, respectively, at 14:00 LT on 18 December. Although BC concentration was the highest at the surface, the largest increase in temperature occurred in the upper layers because of the stronger shortwave absorption efficiency of BC at higher altitude (Ding et al., 2016; Wang et al., 2018). The warming at around 300 m resulted in a more stable stratification, thereby weakening convective motions (Gao et al., 2018). The largest reductions in PBLH were 133.8 m (28.0 %) at 14:00 LT on 12 December and 141.2 m (59.0 %) at 14:00 LT on 18 December in Beijing with the original and modified BC vertical profiles, respectively. On 11 December when regional transport of pollution dominated, relative to the simulation with the original BC profile, simulated air temperature with the modified profile was lower by about 0.5^∘ within the PBL (Fig. 7d), which was caused by the observed maximum mass concentration of BC around 850 m altitude (Fig. 4a). Correspondingly, the maximum reduction in PBLH of 74.2 m was also simulated on 11 December. On 16–19 December when local emissions dominated, compared to the effects of the original BC profiles, the air temperatures at around 300 m were all higher with modified BC. The largest difference of +0.1^∘ was simulated in the PBL on 18 December (Fig. 7d).

Figure 7 (a) Simulated hourly vertical profiles of temperature (contour) and PBLH (red solid line) over Beijing at local time (LT) from 11 to 19 December 2016. (b–d) Time series of changes in vertical temperature (contour) and PBLH (ΔPBLH; red solid line) induced by BC DRE with the original (b; CTRL minus NoBCrad) and modified vertical profiles (c; VerBC_obs minus NoBCrad) and the difference between the effects of the modified and original BC profiles (d; VerBC_obs minus CTRL) over the Beijing region from 11 to 19 December 2016.

Figure 8 shows the spatial distributions of changes in T2, sea-level pressure (SLP), and wind at 10 m caused by BC DRE with the original and modified vertical profiles. BC DRE with both the original and modified vertical profiles produced anomalous northeasterlies in eastern BTH during the two haze events. The mechanism of such changes is that BC DRE induced a strong warming over the Bohai Sea in the east of BTH with a maximum warming of about 1.8^∘, resulting in an anomalous low pressure here and consequently anomalous northeasterlies in eastern BTH (Fig. 8c–d). The similar changes in winds caused by the heating effect of BC were also reported in previous studies (Gao et al., 2016b; Qiu et al., 2017; Chen et al., 2021).

Figure 8 The changes in T2, SLP, and wind at 10 m induced by BC DRE with the original (a, c; CTRL minus NoBCrad) and modified vertical profiles (b, d; VerBC_obs minus NoBCrad) averaged over the period of 12:00–18:00 LT of the two haze events, respectively. (c–d) The northeasterlies in the east of BTH are denoted in red.

4.2 Direct radiative effect of BC on PM_2.5 concentration

By altering the meteorological conditions, BC exerts feedback onto PM_2.5 concentrations. Figure 9 illustrates the impacts of BC DRE with the original and modified vertical profiles on surface layer PM_2.5 as well as the differences in simulated surface layer PM_2.5 between VerBC_obs and CTRL (VerBC_obs minus CTRL) in Beijing during the two haze events. Because of the differences in BC-induced changes in air temperature, wind field, and PBLH, changes in surface layer PM_2.5 concentrations in northern and southern Beijing were different. In the first haze event of 11–12 December, although PBLH was reduced in northern Beijing due to BC DRE, enhanced northerlies brought in relatively clean air to northern Beijing, leading to decreases in surface layer PM_2.5 concentrations with maximum values of 12.5 µg m⁻³ (9.4 %) and 10.6 µg m⁻³ (8.0 %) in this region with the original and modified BC vertical profiles, respectively. Nevertheless, PM_2.5 concentrations increased by up to 17.8 µg m⁻³ (6.6 %) and 24.0 µg m⁻³ (9.3 %) in southern Beijing due to the BC effect with the original and modified vertical profiles, respectively. In the second haze event of 16–19 December, the surface layer PM_2.5 concentrations increased in most areas of Beijing with both vertical profiles. Compared to the simulation with the original profiles, the modified profiles of BC led to larger increases in PM_2.5 concentrations over Beijing, and the maximum differences in PM_2.5 were simulated over central Beijing, which were 9.3 µg m⁻³ (3.6 %) and 5.5 µg m⁻³ (3.1 %) in the first and second haze events, respectively.

Figure 9 The spatial distribution of changes in near-surface PM_2.5 concentrations induced by BC DRE with the original (CTRL minus NoBCrad; a1, a2) and modified vertical profiles (VerBC_obs minus NoBCrad; b1, b2) and the difference between VerBC_obs and CTRL (VerBC_obs minus CTRL; c1, c2) over Beijing averaged over the period of 12:00–18:00 LT of the two haze events. Panels (a1)–(c1) represent the first pollution event of 11–12 December 2016, and (a2)–(c2) represent the second pollution event of 16–19 December 2016. (a3–c3) The daily 6 h contributions of each physical and chemical process (coloured bars, each of which is calculated as the concentration at 18:00 LT minus that at 12:00 LT) to the change in PM_2.5 in Beijing from 11 to 19 December 2016. The black dotted line represents the 6 h net contribution to PM_2.5 change by summing over all processes.

To explain the changes in surface layer PM_2.5 concentrations in Beijing due to BC effects, we carried out process analysis for PM_2.5 for 12:00–18:00 LT of each day when the DRE of BC is the largest (Fig. 9a3, b3, and c3). From 11 to 19 December 2016, VMIX had a dominant positive contribution to changes in PM_2.5 concentration, which reached the maximum contributions of 32.4 and 33.9 µg m⁻³ on December 18 with the original and modified BC vertical profiles, respectively. The vertical mixing was strongly restrained by PBLH; therefore, the decreases in PBLH caused accumulation of PM_2.5 in the lower layers. Meanwhile, CHEM contributed 4.8 and 6.1 µg m⁻³ to PM_2.5 changes because more aerosol precursors restrained in the boundary layer led to the formation of secondary particles. TRA was the major process that had negative contribution to the changes in PM_2.5, which can be explained by the enhanced northerlies in the central part of the NCP due to BC effects as shown in Fig. 8. Relative to the case with the original BC vertical profiles, VMIX and CHEM contributions increased largely with modified profiles, with increases of 8.6 µg m⁻³ (6.5 %) and 7.7 µg m⁻³ (26.8 %), respectively, as averaged over the two haze events, reflecting the further decreases in PBLH (Fig. 7d).

Figure 10 Vertical profiles of the 6 h contributions of physical and chemical processes (coloured bars; each is calculated as the concentration at 18:00 LT minus that at 12:00 LT) to the changes in PM_2.5 induced by BC DRE with the original (CTRL minus NoBCrad; a1, b1) and modified vertical profiles (VerBC_obs minus NoBCrad; a2, b2) and the difference between the original and modified BC profiles (VerBC_obs minus CTRL; a3, b3) over Beijing. The red dotted lines represent the 6 h net contributions to PM_2.5 changes by summing over all processes.

Figure 10 shows the vertical profiles of the contributions of physical–chemical processes to changes in PM_2.5 over Beijing due to BC DRE with the original (CTRL minus NoBCrad; Fig. 10a1 and b1) and modified vertical profiles (VerBC_obs minus NoBCrad; Fig. 10a2 and b2) in the two haze events. In the first haze event of 11–12 December when regional transport of pollution dominated, the NET contribution to PM_2.5 was positive below 256 m, because the positive contribution of VMIX was larger than the negative contribution of TRA. However, in the upper layers (from 256 to 1555 m), the contributions of VMIX and CHEM became negative with both the original and modified vertical profiles, which can be explained by the decreases in PBLH inhibiting the transport of low-layer pollutants to the upper layer. Compared to the original BC vertical profiles, the modified BC vertical profiles increased PM_2.5 in all the vertical layers below 2080 m, in which the positive contribution between 256–757 m was caused by TRA. These results agree with the observed high concentrations of BC at altitudes of 600–1500 m on 11 December (Fig. 4a). In the second haze event, the NET contribution to PM_2.5 was positive below 127 m and negative at 127–504 m. However, the effects of BC on PM_2.5 were small above 504 m because BC concentrations decreased rapidly with altitude.

4.3 Model performance in simulating meteorology and PM_2.5 with the original and modified BC vertical profiles

It is of interest to compare the performance of CTRL (with the original BC vertical profiles) with that of VerBC_obs (with modified BC vertical profiles) in simulating meteorological parameters and PM_2.5 during the two haze events. Figure S4 shows the comparisons between observed T2, RH2, WS10, WD10, and PBLH and the simulated values from the CTRL and VerBC_obs simulations in Beijing in the two haze events (11–12 and 16–19 December 2016). Relative to the CTRL simulation with the original BC vertical profiles, the VerBC_obs simulation with modified BC vertical profiles has better performance in simulating T2, WS10, and PBLH except for WD10 and RH2 in the first pollution event. While the MBs of T2, RH2, WS10, WD10, and PBLH are 0.2^∘, 0.0 %, −0.4 m s⁻¹, 3.8^∘, and 45.0 m in CTRL, they are 0.0^∘, 2.2 %, −0.1 m s⁻¹, 10.9^∘, and 29.0 m in VerBC_obs, respectively (Table S3). In the second pollution event, the positive bias in PBLH (MB = 43.7 m, NMB = 42.9 %) in CTRL is reduced to 33.9 m and 33.3 % in VerBC_obs.

Table S4 shows the statistical comparison between observed hourly surface layer PM_2.5 and the model results from CTRL and VerBC_obs in Beijing for each day during the two haze events. The model with modified BC vertical profiles can enhance the capability in simulating the temporal variation in PM_2.5 for each day; the correlation coefficient between simulated hourly concentrations and hourly observations in each day of the studied period increased from 0.04–0.84 in CTRL to 0.24–0.93 in VerBC_obs.

5 Roles of BC vertical profiles

BC has higher light-absorbing efficiency at higher altitudes (Ding et al., 2016; Wang et al., 2018). As described in Sect. 2.3, the observed vertical profiles of BC on heavily polluted days (12 and 16–19 December) can be parameterized as exponential decline functions using nonlinear regression ($C\left(h\right)=C_{\mathrm{0}}\times e^{-h/\mathrm{hs}}$) with hs values of 0.35, 0.48, 0.53, 0.79, 0.82, and 0.96, and the profiles affected by regional transport had high concentrations of BC at high altitudes. We conducted seven sensitivity experiments which applied six exponential functions and one observed transport-dominated vertical profile, as described in Sect. 2.4, to examine the roles of BC vertical profiles in influencing meteorological conditions and PM_2.5 during severe haze events. In these sensitivity experiments, we only modify the BC vertical profiles for the dates of 12 and 16–19 December. In the function of $C\left(h\right)=C_{\mathrm{0}}\times e^{-h/\mathrm{hs}}$, a larger hs means less BC at the surface and more BC in the upper atmosphere, as shown in Fig. 11.

Figure 11 Vertical profiles of BC concentrations parameterized as six exponential functions for 12 and 16–19 December 2016.

5.1 Impacts of BC vertical profiles on meteorology

Figure 12 shows the simulated changes in atmospheric temperature induced by BC DRE with exponential functions (VerBC_hs1-6 minus NoBCrad) and with the transport-dominated vertical profile (VerBC_RT minus NoBCrad). BC had a significant warming effect at altitudes of 256–421 m from 12:00 to 18:00 LT (Fig. 7). Generally, with the value of hs gradually increasing, the BC-induced warming in the afternoon around 300 m became smaller, which can be explained by the highest mass fraction of BC at the altitudes of 256–421 m to total BC column burden in VerBC_hs1 case (31.7 %) and the lowest percentage in the VerBC_hs6 case (21.7 %) among the six sensitivity experiments (Fig. S1). The maximum warming around 300 m was 0.42^∘ in the VerBC_hs1 case and 0.19^∘ in the VerBC_hs6 case. It should be noted that BC led to a significant cooling effect at the surface (below 80 m) when hs values were 0.79, 0.82, and 0.96, with changes in temperature by −0.08, −0.09, and −0.13^∘, respectively. Because more BC mass was assigned to high altitudes (above 1000 m) with higher hs, less solar radiation could reach the ground (Fig. S5). These results are consistent with those found in previous modelling and observational studies (Cappa et al., 2012; Ferrero et al., 2014; Ding et al., 2016; Wang et al., 2018). Meanwhile, in the case of VerBC_RT, BC also had a cooling effect of −0.30^∘ at the surface (Fig. 12g). In current regional air quality models, the uncertainties in BC profiles could influence the capability of a model to simulate a cooling effect of BC on surface air temperature (Wang et al., 2019).

Figure 12 Time series of changes in vertical temperature induced by BC DRE with six exponential functions (VerBC_hs1-6 minus NoBCrad) and one transport-dominated vertical profile (VerBC_RT minus NoBCrad) averaged over 12 and 16–19 December 2016.

Figure 13 (a) Variation in ΔT_BC caused by BC DRE with increasing hs values averaged for 12 and 16–19 December. The black dashed line is the linear fit. (b) Time series of changes in PBLH in Beijing caused by different BC vertical profiles averaged for 12 and 16–19 December 2016.

We further use the difference in temperature between the upper PBL (T_H; 256–421 m) and the ground (T_L; 0–127 m) ($\mathrm{\Delta }{T}_{\mathrm{BC}}=T_{\mathrm{H}}-T_{\mathrm{L}}$) averaged over 12:00–18:00 LT of 12 and 16–19 December to quantify temperature inversion caused by BC DRE. BC aerosol leads to cooling at the surface and warming in the upper PBL, both of which weaken the convection in the boundary layer and consequently reduce the PBLH (Ding et al., 2016; Wang et al., 2018; Chen et al., 2021). In our study, with hs values increasing from 0.35 to 0.96, ΔT_BC increased from 0.17 to 0.42^∘, and the ΔT_BC value was 0.51^∘ in the VerBC_RT case (Fig. 13a). The larger ΔT_BC indicates stronger cooling at the surface. Such temperature inversion at 12:00–18:00 LT resulted in more stable stratification and further inhibited the development of PBL. The cooling at the surface also reduced sensible heat flux from the surface (Fig. S5), suppressing vertical turbulence and hence reducing PBLH (Wilcox et al., 2016). As a result, the reductions in PBLH were larger with higher hs (Fig. 13b). The minimum decrease in PBLH was −31.9 m (−14.3 %) with a hs value of 0.35, and the maximum decrease was 48.9 m (22.0 %) with a hs value of 0.96, as averaged over the period of 12:00–18:00 LT on 12 and 16–19 December. In the case of VerBC_RT, the mean PBLH was reduced by 56.9 m (25.6 %) during the period of 12:00–18:00 LT.

Figure 14 (a) Time series of the changes in surface layer PM_2.5 in Beijing caused by BC with six exponential functions (VerBC_hs1-6 minus NoBCrad) and one observed transport-dominated vertical profile (VerBC_RT minus NoBCrad) averaged for 12 and 16–19 December 2016. (b–d) The hourly contributions of each physical and chemical process to PM_2.5 changes caused by BC DRE with two exponential functions (hs = 0.35 and 0.96) and one transport-dominated vertical profile.

5.2 Impacts of BC vertical profiles on PM_2.5 concentration

Figure 14a shows the changes in surface layer PM_2.5 concentration caused by BC DRE with six exponential functions (VerBC_hs1-6 minus NoBCrad) and the transport-dominated vertical profile (VerBC_RT minus NoBCrad) averaged over 12 and 16–19 December 2016. From 00:00 to 11:00 LT, the surface layer PM_2.5 exhibited larger BC-induced decreases with a higher value of hs. This can be explained by the negative contribution of the TRA process that increased during the period of 00:00–05:00 LT (Fig. 14b and c) when the hs value changed from 0.35 to 0.96. The surface layer PM_2.5 concentration was reduced by up to 9.1 µg m⁻³ (6.2 %) and 12.6 µg m⁻³ (8.6 %) at 05:00 LT with hs values of 0.35 and 0.96, respectively. Compared to the NoBCrad case, the surface layer PM_2.5 concentrations were reduced by up to 13.8 µg m⁻³ (9.4 %) at 05:00 LT due to BC DRE in the VerBC_RT case. From 12:00 to 18:00 LT, the BC-induced increase in surface layer PM_2.5 concentrations was larger as hs values were higher; relative to NoBCrad simulation, the mean PM_2.5 concentrations were increased by 5.5 µg m⁻³ (3.4 %) and 7.9 µg m⁻³ (4.9 %) with the hs values of 0.35 and 0.96, respectively. Because the PBL was suppressed by BC DRE from 12:00 to 15:00 LT, the contributions of VMIX and CHEM to surface layer PM_2.5 were positive and larger in magnitude than the negative contribution of TRA. The NET of all processes was negative from 16:00 to 18:00 LT due to the continuous growth of a negative contribution of TRA. The negative contribution of the TRA process from 12:00 to 18:00 LT can be explained by the enhanced northerlies in the central part of BTH caused by BC DRE, which transported cleaner air masses into Beijing (Fig. S6). From 19:00 to 23:00 LT, the surface layer PM_2.5 concentrations were decreased by BC DRE, which can be explained by the dominant negative contribution of TRA from 19:00 to 21:00 LT. At 22:00 LT, the reduction in surface layer PM_2.5 was 7.5 µg m⁻³ (4.0 %) when the hs value was 0.35 and 6.6 µg m⁻³ (3.5 %) when hs was 0.96.

6 Conclusions

In this study, a fully coupled online WRF-Chem model with an improved integrated process rate (IPR) analysis scheme is employed to investigate the direct radiative effects (DRE) of BC vertical profiles on meteorology and PM_2.5 concentrations during two severe haze events (11–12 and 16–19 December 2016). Sensitivity experiments are conducted to compare the DRE of BC with the original and corrected vertical profiles and to quantify the role of BC vertical profiles in influencing meteorological conditions and air pollutants.

Compared to the measured vertical profiles of BC in Beijing, the default vertical profiles of BC from the WRF-Chem model can capture the decreases in BC mass concentration with altitude on 12 and 16–19 December when local emissions dominated, but they cannot reproduce the observed maximum mass concentration of BC around 850 m altitude on 11 December when regional transport of pollutants dominated. Averaged over the two severe pollution events, the model overestimated BC mass concentration by 87.4 % at the surface but underestimated BC by 33.1 % at 1000 m altitude compared with the observations in Beijing.

We carried out simulations with both the default original BC vertical profiles and the modified vertical profiles using the observations (keep the column burden of BC from WRF-Chem but distribute BC mass vertically according to the observed fractions of BC in individual model layers for each day). Comparing the simulation with original BC profiles, the warming by BC DRE around 300 m altitude was stronger with the modified profiles. Accordingly, the BC-induced reductions in PBLH in Beijing averaged over the two severe haze events were 43.4 m (18.4 %) and 55.4 m (23.5 %), respectively, with the original and modified profiles. As a result, relative to the simulation with the original profiles, the modified profiles of BC led to larger increases in PM_2.5 concentrations by BC DRE, because VMIX and CHEM had the dominant positive contributions to the changes in surface layer PM_2.5 due to the reductions in PBLH.

Seven sensitivity experiments were further carried out to understand the roles of BC vertical profiles. In six assumed exponential functions ($C\left(h\right)=C_{\mathrm{0}}\times e^{-h/\mathrm{hs}}$) with hs values of 0.35, 0.48, 0.53, 0.79, 0.82, and 0.96, a larger hs means less BC at the surface and more BC in the upper atmosphere. As a result, with higher hs, the surface had less warming or larger cooling, leading to stronger temperature inversion and hence larger BC-induced increases in PM_2.5.

Results from our study highlight the importance of accurate representation of BC vertical profiles in models, which alter the radiation balance, BC-PBL interaction, and hence the simulated PM_2.5 concentrations. Due to the limitation of observational data, this study was focused on the DRE of BC vertical profiles on meteorology and PM_2.5 concentration in Beijing during severe haze events. However, the results from this study should be generally important for understanding severe haze for urban areas.

There are channels for further improvement in near-future research. We distribute BC mass vertically according to the observed fractions of BC in individual model layers for each day without considering the hourly variations in BC vertical profiles due to the lack of data. Such assumed distribution of BC based on observation may not be consistent with the dynamical (winds, temperature, etc.) and chemical processes of the atmosphere. Further efforts are needed to examine the roles of BC vertical profiles in coupled chemistry–weather models.