2.1 Model configuration

A fully coupled online Weather Research and Forecasting with Chemistry model (WRF-Chem v3.7) is used to simulate meteorological fields and concentrations of gases and aerosols simultaneously (Skamarock et al., 2008; Grell et al., 2005). The WRF-Chem model is designed with two domains using 219 (west–east) ×159 (south–north) and 150 (west–east) ×111 (south–north) grid points at the horizontal resolutions of 27 and 9 km, respectively (Fig. 1). The outer domain covers nearly the whole of East Asia, and the inner domain is located in the NCP. In order to minimize the impacts from LBCs (lateral boundary conditions), we only analyze the simulation results from the inner region of the second domain (i.e., BTH), following Chen et al. (2018) and Wu et al. (2012). The vertical dimension is resolved by 29 full sigma levels, with 15 layers located in the bottom 2 km for finer resolution in the PBL; the height of the first layer averaged in BTH is about 30 m.

Figure 1(a) Map of the two nested model domains. (b) Locations of the observations used for model evaluation.

Figure 2Time series of the observed (black dots) and simulated (red dots) hourly 2 m temperature (T₂, K), 2 m relative humidity (RH₂, %), 10 m wind speed (WS₁₀, m s⁻¹), and 10 m wind direction (WD₁₀, ^∘) averaged over the 12 stations from 16 to 29 December 2015.

Download

Meteorological initial and lateral boundary conditions used in the WRF-Chem model are taken from the NCEP (National Center for Environmental Prediction) (Final) Operational Global Analysis data with a spatial resolution of ${\mathrm{1}}^{\circ }\times {\mathrm{1}}^{\circ }$. Four-dimensional data assimilation (FDDA) with the nudging coefficient of $\mathrm{3.0}\times {\mathrm{10}}^{-\mathrm{4}}$ for wind (in and above the PBL), temperature (above the PBL), and water vapor mixing ratio (above the PBL) is adopted to improve the accuracy of simulation results (no analysis nudging is included for the inner domain) (Lo et al., 2008; Otte, 2008; L. Wang et al., 2016; Werner et al., 2016). The forecasts from the MOZART-4 global chemical transport model are processed to provide the chemical initial and boundary conditions for the WRF-Chem model (Emmons et al., 2010).

Anthropogenic emission data are obtained from the MIX Asian emission inventory (http://www.meicmodel.org/dataset-mix.html, last access: 12 August 2019), with a horizontal resolution of 0.25^∘ (M. Li et al., 2017). MIX is developed to support the MICS-Asia III (Model Inter-Comparison Study for Asia Phase III) and the TF HTAP (Task Force on Hemispheric Transport of Air Pollution) projects. This inventory includes SO₂ (sulfur dioxide), NO_x (nitrogen oxides), CO (carbon monoxide), CO₂ (carbon dioxide), NMVOC (non-methane volatile organic compounds), NH₃ (ammonia), BC (black carbon), OC (organic carbon), PM_2.5, and PM₁₀. All of these species are from several sectors, such as agriculture, industry, power, transportation, and residential, and the emission rate of each species for each hour is based on Gao et al. (2015). The biogenic emissions are calculated online using the MEGANv2.04 (Model of Emission of Gases and Aerosol from Nature v2.04) model (Guenther et al., 2006). Biomass-burning emissions are obtained from the GFEDv3 (Global Fire Emissions Database v3) (Randerson et al., 2005). Dust emissions and sea salt emissions are calculated online using algorithms proposed by Shao (2004) and Gong et al. (1997), respectively.

The Carbon-Bond Mechanism version Z (CBMZ) (Zaveri and Peters, 1999) is selected to simulate the gas-phase chemistry, and the eight-bin sectional aerosol module, MOSAIC (Model for Simulating Aerosol Interactions and Chemistry) (Zaveri et al., 2008), with some aqueous chemistry, is used to simulate aerosol evolution. All major aerosol species are considered in the MOSAIC scheme, including sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$), chloride (Cl), sodium (Na), BC, primary organic mass, liquid water, and other inorganic mass (Zaveri et al., 2008). The aerosol size distribution is divided into discrete size bins defined by their lower and upper dry particle diameters (Zhao et al., 2010). In the current CBMZ/MOSAIC scheme, the formation of SOA (secondary organic aerosol) is not included (Zhang et al., 2012; Gao et al., 2016). Aerosol optical properties, including the extinction efficiency, the single scatter albedo, and the asymmetry factor are computed using Mie theory, based on aerosol composition, mixing state, and size distribution (Barnard et al., 2010). The impacts of aerosols on photolysis rates are calculated using the Fast-J photolysis scheme (Wild et al., 2010). Aerosol radiation is simulated using RRTMG (Rapid Radiative Transfer Model for GCMs) for both shortwave (SW) and longwave (LW) radiation (Zhao et al., 2011). More information regarding the parameterizations used in this study can be found in Table 1.

Table 1Parameterizations used in the WRF-Chem model.

Download Print Version | Download XLSX

2.2 Integrated process rate (IPR) analysis

Most air quality models are configured to output only the pollutant concentrations that reflect the combined effects of all physical and chemical processes. Quantitative information on the impacts of individual process is usually unavailable. Process analysis techniques, i.e., integrated process rate (IPR) analysis, can be used in grid-based Eulerian models (e.g., WRF-Chem) to obtain contributions of each physical/chemical process to variations in pollutant concentrations. Eulerian models utilize the numerical technique of operator splitting to solve continuity equations for each species into several simple ordinary differential equations or partial differential equations that only contain the influence of one or two processes (Gipson, 1999).

The IPR analysis method has been fully implemented in Community Multi-scale Air Quality (CMAQ) model, and has been widely applied to study regional photochemical ozone (O₃) pollution (Gonçalves et al., 2009; Khiem et al., 2010; Xing et al., 2017; Tang et al., 2017). Several WRF-Chem model studies have used the IPR analysis to investigate the impacts of physical/chemical process on variations in O₃ concentrations. Gao et al. (2018) investigated the impacts of BC–PBL interactions on O₃ concentrations by analyzing the contributions from photochemistry, vertical mixing, and advection processes. Jiang et al. (2012) calculated the contributions of photochemical reactions and physical processes to O₃ formation using a simplified IPR analysis scheme.

Applying the IPR analysis to diagnose the contributions of each physical or chemical process to variations in aerosol concentrations in the WRF-Chem model is more complex technically; therefore, few studies have utilized the IPR analysis for aerosols. In this study, we developed an improved IPR analysis scheme in the WRF-Chem model to isolate the processes impacting variations in aerosol concentrations into nine different processes, namely advection (TRAN), emission source (EMIS), dry deposition (DYRD), turbulent diffusion (DIFF), sub-grid convection (SGCV), gas-phase chemistry (GASC), cloud chemistry (CLDC), aerosol chemistry (AERC), and wet scavenging (WETP). TRAN includes horizontal and vertical advection, which is highly related to wind and aerosol concentration gradients from upwind regions to downwind areas (Gao et al., 2018). DRYD is based on resistance models for trace gases (Wesely, 1989) and aerosol particles (Ackermann et al., 1998). SGCV refers to the scavenging within the sub-grid wet convective updrafts. CLDC refers to aqueous-phase photolytic and radical chemistry reactions in clouds, including the activation processes. AERC refers to microphysical nucleation, condensation, and coagulation, as well as the mass transfer between the gas phase and condensed phase. WETP contains in-cloud rainout and below-cloud washout during grid-scale precipitation. The contribution of individual processes can be calculated as the difference of aerosol concentrations before and after the corresponding operator.

Based on the principle of mass balance, IPR can be verified by comparing the variations in aerosol concentrations (the concentration at the current time minus the concentration at the previous time) with the sum of the contributions from the nine processes during each time step. As shown in Fig. S1 in the Supplement, the net contributions of all processes match the variations in aerosol concentrations quite well.

2.3 Numerical experiments

Table 2 summarizes the experimental designs. To investigate the contributions of regional transport and local emissions to the PM_2.5 concentrations in the BTH region, four simulations with different anthropogenic emission categories were conducted: (1) CTL – the control simulation with all anthropogenic emissions considered; (2) NoAnth – no anthropogenic emissions are considered in the whole domain; (3) NoBTH_Anth – the same as the CTL, but anthropogenic emissions in the BTH area are excluded; and (4) OnlyBTH_Anth – contrary to the NoBTH_Anth case, anthropogenic emissions are only considered in the BTH region. All the physical and chemical schemes used in these cases are identical. The contributions of regional transport and local emissions to the PM_2.5 concentration in the BTH region can be identified by comparing the simulation results of NoBTH_Anth and NoAnth (i.e., NoBTH_Anth minus NoAnth) and OnlyBTH_Anth and NoAnth (i.e., OnlyBTH_Anth minus NoAnth), respectively.

Table 2Experimental design. “Y” represents yes, and “N” represents no.

Download Print Version | Download XLSX

To quantify the aerosol radiative effects (ARE) on haze pollution, another sensitivity experiment (referred to as the NoARE case) was designed by turning the feedbacks between aerosols and meteorological variables off, including eliminating the aerosol direct effect (ADE) and the aerosol indirect effect (AIE) in the model. The ADE is turned off by removing the mass of aerosol species from the calculation of aerosol optical properties, following Qiu et al. (2017). The AIE is turned off using a prescribed vertically uniform cloud droplet number, which is calculated from the CTL case during the whole simulation period, following Gao et al. (2015) and B. Zhang et al. (2015). The differences between CTL and NoARE (i.e., CTL minus NoARE) represent the impacts of aerosol radiative forcing.

The IPR analysis method is applied to all of the experiments designed. Comparing the contributions of each detailed process between the pollution accumulation stage and the dissipation stage in the CTL case can quantitatively explain the reason for the variation in the PM_2.5 concentrations in the BTH region. Meanwhile, the prominent physical or chemical process responsible for the aerosol radiative impacts on the haze episode can also be investigated by analyzing the IPR analysis method used in the CTL and NoARE cases.

All five simulations are conducted for the period from 13 to 29 December 2015, and the initial 3 days are discarded as the model spin-up to minimize the impacts of initial conditions. Simulation results from the CTL case from 16 to 29 December 2015 are used to evaluate the model performance.

2.4 Observational data

Simulated meteorological parameters in CTL case, including 2 m temperature (T₂), 2 m relative humidity (RH₂), 10 m wind speed (WS₁₀), and 10 m wind direction (WD₁₀), are compared with hourly observations at 12 stations, which are collected from NOAA's National Climatic Data Center (https://gis.ncdc.noaa.gov/maps/ncei/cdo/hourly, last access: 12 August 2019). Due to limited observations of the PBL height in the BTH area, the retrieved PBLH observations at 3 h intervals obtained from the GDAS (Global Data Assimilation System) (https://ready.arl.noaa.gov/READYamet.php, last access: 12 August 2019) in Beijing (39.93^∘ N, 116.28^∘ E) are also used to evaluate the model performance. More detailed information about the GDAS meteorological dataset ($\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$) can be found in Rolph (2013), Kong et al. (2015), and at https://www.ready.noaa.gov/gdas1.php (last access: 12 August 2019). The hourly shortwave downward radiation flux (SWDOWN) at the Xianghe station (39.75^∘ N, 116.96^∘ E) is taken from WRMC-BSRN (World Radiation Monitoring Center-Baseline Surface Radiation Network, http://bsrn.awi.de, last access: 12 August 2019) for the energy budget evaluation. The hourly observed surface-layer PM_2.5 concentrations at the 59 stations are obtained from the CNEMC (China National Environmental Monitoring Center, http://www.cnemc.cn/, last access: 12 August 2019). The daily measurements of the mass concentrations of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, BC, and OC are collected at the Beijing (39.97^∘ N, 116.37^∘ E) and Shijiazhuang (38.03^∘ N, 114.53^∘ E) sites (Huang et al., 2017; Z. Liu et al., 2018). Detailed locations of these observations are shown in Fig. 1b.

3.1 Meteorological parameters

Figure 2 shows the time series of observed and simulated hourly meteorological variables averaged over the 12 stations from 16 to 29 December 2015. Corresponding statistical metrics, including the mean value, the normalized mean bias (NMB), the mean fractional bias (MFB), the mean fractional error (MFE), the index of agreement (IOA), and the correlation coefficient (R) are presented in Table 3. As shown in Fig. 2, simulated T₂, RH₂, WS₁₀, and WD₁₀ agree well with the observational data. For temperature, the WRF-Chem model can perfectly depict its diurnal and daily variations with R and IOA values of 0.90 and 0.94, respectively, but slightly overestimates the low values at night, with a NMB of 1 %. Observed relative humidity can be reasonably reproduced by the model with R and IOA values of 0.73 and 0.82, respectively, but a persistent underestimation is found with a NMB of −12 %. Different surface layer and boundary layer parameterizations may influence the simulated near-surface moisture fluxes, and the settings of these schemes can partially explain the biases of RH₂ between the observations and simulations (Qian et al., 2016). This negative bias of RH₂ can also be simulated by other studies (Zhang et al., 2009; Gao et al., 2015). WRF-Chem can capture the observed low wind speed values from 19 to 23 December and high wind speed values from 16 to 17 and 25 to 27 December. The positive NMB of 28 % probably results from unresolved topographical features in the surface drag parameterization and the coarse resolution used in the nested domain (Yahya et al., 2015; Zheng et al., 2015). For wind direction, the calculated NMB is 1 % and the IOA is 0.65, indicating that the WRF-Chem model can generally reproduce the varied wind direction during the simulation period.

Table 3Statistical metrics between observations and simulations.

^a,b $\stackrel{\mathrm{‾}}{\mathrm{OBS}}$ and $\stackrel{\mathrm{‾}}{\mathrm{SIM}}$ represent the average observations and simulations, respectively. $\stackrel{\mathrm{‾}}{\mathrm{OBS}}=\frac{\mathrm{1}}{\mathrm{n}}\times {\sum }_{i=\mathrm{1}}^{\mathrm{n}}{\mathrm{OBS}}_i$, $\stackrel{\mathrm{‾}}{\mathrm{SIM}}=\frac{\mathrm{1}}{\mathrm{n}}\times {\sum }_{i=\mathrm{1}}^{\mathrm{n}}{\mathrm{SIM}}_i$.
^c NMB is the normalized mean bias, $\mathrm{NMB}=\frac{\mathrm{1}}{\mathrm{n}}\times {\sum }_{i=\mathrm{1}}^{\mathrm{n}}\frac{{\mathrm{SIM}}_i-{\mathrm{OBS}}_i}{{\mathrm{OBS}}_i}\times \mathrm{100}\mathit{\%}$.
^d MFB is the mean fractional bias, $\mathrm{MFB}=\frac{\mathrm{2}}{\mathrm{n}}\times {\sum }_{i=\mathrm{1}}^{\mathrm{n}}\frac{{\mathrm{SIM}}_i-{\mathrm{OBS}}_i}{{\mathrm{SIM}}_i+{\mathrm{OBS}}_i}\times \mathrm{100}\mathit{\%}$.
^e MFE is the mean fractional error, $\mathrm{MFE}=\frac{\mathrm{2}}{\mathrm{n}}\times {\sum }_{i=\mathrm{1}}^{\mathrm{n}}\frac{\left|{\mathrm{SIM}}_i-{\mathrm{OBS}}_i\right|}{{\mathrm{SIM}}_i+{\mathrm{OBS}}_i}\times \mathrm{100}\mathit{\%}$.
^f IOA is the index of agreement, $\mathrm{IOA}=\mathrm{1}-\frac{{\sum }_{i=\mathrm{1}}^{\mathrm{n}}{\left({\mathrm{SIM}}_i-{\mathrm{OBS}}_i\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^{\mathrm{n}}{\left(|{\mathrm{OBS}}_i-\stackrel{\mathrm{‾}}{\mathrm{OBS}}|+|{\mathrm{SIM}}_i-\stackrel{\mathrm{‾}}{\mathrm{SIM}})|\right)}^{\mathrm{2}}}$.
^g R is the correlation coefficient, $R=\frac{{\sum }_i^{\mathrm{n}}\left|\right({\mathrm{OBS}}_i-\stackrel{\mathrm{‾}}{\mathrm{OBS}})\times ({\mathrm{SIM}}_i-\stackrel{\mathrm{‾}}{\mathrm{SIM}}\left)\right|}{\sqrt{{\sum }_i^{\mathrm{n}}({\mathrm{OBS}}_i-\stackrel{\mathrm{‾}}{\mathrm{OBS}}{)}^{\mathrm{2}}+{\sum }_i^{\mathrm{n}}({\mathrm{SIM}}_i-\stackrel{\mathrm{‾}}{\mathrm{SIM}}{)}^{\mathrm{2}}}}$.
In the above OBS_i and SIM_i refer to observations and model predictions, respectively, i refers to a given station, and n is the total number of stations.
^h T₂: temperature at 2 m (K); RH₂: relative humidity at 2 m (%); WS₁₀: wind speed at 10 m (m s⁻¹); WD₁₀: wind direction at 10 m (^∘).

Download Print Version | Download XLSX

Simulated hourly PBLH and SWDOWN are also compared with observations in Fig. 3. It is noted that the PBLH measurements provided by GDAS of NOAA are in 3 h intervals. The simulations in the CTL case agree well with the observations, including capturing the daily maximum at daytime and the low values at night. The correlation coefficients are 0.68 and 0.91 for PBLH and SWDOWN, respectively.

Figure 3Time series of the observed (black dots) and simulated (red lines) hourly planetary boundary layer height (PBLH, m) at the site in Beijing (39.93^∘ N, 116.28^∘ E), and shortwave downward radiation flux (SWDOWN, W m⁻²) at the Xianghe station (39.75^∘ N, 116.96^∘ E) from 16 to 29 December 2015. Notably, PBLH measurements provided by Global Data Assimilation System (GDAS) are in 3 h intervals. All times are converted to China standard time (Beijing time).

Download

3.2 PM_2.5 and its components

Observed hourly surface-layer PM_2.5 concentrations from 16 to 29 December 2015 in the nine cities (Shengyang, Beijing, Xingtai, Hengshui, Baoding, Langfang, Yangquan, Anyang, and Jinan) are compared with the model results from the CTL case (Fig. 4). The statistical metrics are shown in Table 3. Generally, the WRF-Chem model can reasonably reproduce the evolutional characteristics of the observed PM_2.5 concentrations in the nine cities (Rs =0.57–0.90). Both the observed and simulated PM_2.5 concentrations exhibit a growth trend from 16 to 22 and 28 to 29 December, and a decreasing tendency from 23 to 27 December. However, an obvious underestimation is found in Beijing from 25 to 26 December when a maximum hourly concentration of 600 µgm⁻³ was observed. This negative bias is also simulated by previous studies (Chen et al., 2018; Z. Zhang et al., 2018). The possible reasons for the underestimation are as follows: (1) the bias in simulated meteorological conditions (e.g., underestimated RH₂ and overestimated WS₁₀); (2) the missing mechanisms of some gas–aerosol phase partitioning and heterogeneous reactions which may produce secondary inorganic aerosol (X. Huang et al., 2014; Wang et al., 2014); and (3) the lack of SOA simulation in the MOSAIC mechanism (Gao et al., 2016). Generally, the performance statistics of PM_2.5 in almost all cities meet the model performance goal (MFB within ±30 % and MFE ≤50 %) proposed by Boylan and Russel (2006).

Figure 4Time series of the observed (black dots) and simulated (red dots) hourly PM_2.5 concentrations (µg m⁻³) in the nine cities (Shengyang, Beijing, Xingtai, Hengshui, Baoding, Langfang, Yangquan, Anyang, and Jinan) from 16 to 29 December 2015. The “n” in each panel represents the number of observation sites in each city. Beijing time is used for these hourly time series.

Download

Figure 5 compares the simulated and observed surface-layer concentrations of BC, OC, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in Beijing and Shijiazhuang averaged from 16 to 29 December 2015. The WRF-Chem model underestimates the concentrations of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and OC in Beijing (Shijiazhuang) by 19 % (40 %), 14 % (9 %), and 21 % (41 %), respectively, but overestimates the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration by 29 % (44 %). Due to the low reactivity of BC in the atmosphere, the uncertainty in the BC emissions may cause the biases in Beijing (NMB $=+\mathrm{10}$ %) and Shijiazhuang (NMB $=-\mathrm{24}$ %). For OC, the underestimation may result from the lack of SOA in the MOSAIC aerosol module (Qiu et al., 2017). Missing some SO₂ gas-phase and aqueous-phase oxidation mechanisms, as well as heterogeneous chemistry may explain the underestimation of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. It is noted that similar biases of aerosol components were also reported by other WRF-Chem studies (B. Zhang et al., 2015; Qiu et al., 2017).

Figure 5Comparison of observed and simulated surface-layer mass concentrations (µg m⁻³) of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (red), ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (blue), ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (purple), OC (green), and BC (gray) at sites (a) in Beijing (39.97^∘ N, 116.37^∘ E) and (b)  Shijiazhuang (38.03^∘ N, 114.53^∘ E) averaged over the 16–29 December 2015 period. Normalized mean biases (NMBs) are also listed for each species (colored numbers).

Download

4.1 Spatial–temporal evolutions of surface-layer PM_2.5 concentrations

Figure 6a–k show the spatial distributions of the simulated daily mean surface-layer PM_2.5 concentrations from 17 to 28 December 2015. From 17 December, aerosol particles started to accumulate in the near-surface layer in the BTH region under a prevailing southerly air flow. On 20 December, the BTH area was under a uniform pressure field (Fig. S2a). The regional average wind speed was less than 3 m s⁻¹, and the boundary layer became stable, which constrained aerosols within a low mixing layer. Meanwhile, a low-pressure center was situated to the north of the BTH region, where air pollutants from south, southwest, and southeast converged. Consequently, the daily mean PM_2.5 concentration averaged over the BTH area was over 200 µg m⁻³. On 21 December, a weak low-pressure center formed near Bohai Bay and a weak high-pressure center moved to the Shandong Peninsula (Fig. S2b). The synoptic conditions brought more air masses from south to north, and worsened air quality in the BTH region. On 22 December, a weak high-pressure system moved within Inner Mongolia (Fig. S2c), which carried cold air to the BTH region. Meanwhile, the polluted air was also transported back to the BTH, leading to a continuous increase in the PM_2.5 concentration, with the maximum daily mean value exceeding 600 µg m⁻³ (Fig. 6e). Due to the enhanced anticyclone originating from Siberia (Fig. S2d), the accumulation of aerosol particles in the BTH region was terminated by the incursion of a strong cold front from 23 to 27 December. However, frequent transitions between high- and low-pressure systems over the BTH area accompanied by shifting wind directions resulted in a rapid PM_2.5 variation, especially on 24 and 25 December, when a low-pressure system developed northeast of BTH (Fig. S2e). The air mass over BTH was influenced by the pollutants from the south, resulting in a temporary increase in the concentration of PM_2.5 on 25 December. After 27 December, another haze episode gradually formed.

Figure 6(a–k) Spatial distributions of simulated daily PM_2.5 concentrations (shaded, µg m⁻³) and wind vectors (arrows, m s⁻¹). Time series of simulated daily PM_2.5 concentrations averaged over the Beijing–Tianjin–Hebei region are also shown in panel (l).

According to the trends in simulated PM_2.5 concentrations averaged over the BTH region (Fig. 6l), we divide the whole simulation period into three stages: (1) the aerosol accumulation stage (16–22 December, Stage 1), (2) the aerosol dispersion stage (23–27 December, Stage 2), and (3) the formation stage for another haze event (28–29 December, Stage 3). In this paper, we mainly focus on the first two stages to reveal important factors that cause the accumulation and dispersion of particulate matter.

In Stage 1, the daily mean PM_2.5 concentrations averaged over the BTH region increased from 24.2 to 289.8 µg m⁻³, and the average PM_2.5 concentration was 145.6 µg m⁻³ (Fig. 7a), which is close to the “heavily polluted” air quality threshold value (PM_2.5 24 h average concentration >150 µg m⁻³). The WS₁₀ was low (Fig. 7b), especially during the heavily pollution period (20–22 December), and the mean wind speed was 2.3 m s⁻¹, which is less than 3.2 m s⁻¹ (one of the indicators used to define air stagnation by NOAA, https://www.ncdc.noaa.gov/societal-impacts/air-stagnation/overview, last access: 12 August 2019), indicating that the near-surface circulation was insufficient to disperse accumulated air pollutants. The decreased PBLH (from 701.6 to 109.9 m) could compress air pollutants into a shallow layer, resulting in an elevated pollution level. During Stage 2, the PM_2.5 concentration decreased gradually with the increased wind speed and PBLH. The average PM_2.5 concentration during Stage 2 was 107.9 µg m⁻³, which still exceeded the Grade II standard (75 µg m⁻³) defined by the National Ambient Air Quality Standards of China.

Figure 7(a) Contributions of local emissions (red) and regional transport (blue) to the near-surface PM_2.5 concentrations averaged over the Beijing–Tianjin–Hebei region from 16 to 29 December 2015. The absolute contributions (µg m⁻³) are shown using bars, and the percentage contributions (%) are shown using lines. The PM_2.5 concentration and the percentage contributions averaged over each stage are listed at the top of panel (a). Simulated daily 10 m wind speed (WS₁₀; m s⁻¹; black dotted line), specific humidity (g kg⁻¹; green dotted line), and PBLH (m; magenta dotted line) averaged over Beijing–Tianjin–Hebei are also shown in panel (b).

Download

4.2 Contributions of local emissions and regional transport to PM_2.5 concentrations

Previous studies have reported that anthropogenic emissions are the dominant cause of haze events in China (Jiang et al., 2013; Sun et al., 2014; Gu and Liao, 2016; Y. Yang et al., 2016). Emission control measures have been taken to ensure good air quality for major events (e.g., APEC) or to mitigate the severity of coming pollution episodes (Zhou et al., 2018). Other studies, such as Sun et al. (2017) and Wang et al. (2017), have pointed out that regional transport contributed more than 50 % of the particulate concentrations in the BTH region during haze events. This section discusses the contributions of local anthropogenic emission and regional transport to the PM_2.5 concentration in the BTH area, aiming to reveal their relative importance during this haze episode.

As shown in Fig. 7a, the PM_2.5 concentration in BTH during Stage 1 was mainly contributed by the combined effects of local emissions and regional transport. The contributions of local emissions and regional transport to the PM_2.5 concentration were comparable (49 % and 32 %, respectively), especially during the heavy pollution period (20–22 December, 43 % vs. 37 %). During Stage 2, the contributions of regional transport decreased from 30 % to 16 %. The relative high PM_2.5 concentration (107.9 µg m⁻³) was principally caused by the local emissions. On average, the contributions of local emissions and regional transport to the PM_2.5 concentration in Stage 2 were 51 % and 24 %, respectively. The impact of regional transport could be qualitatively expressed by specific humidity, which was treated as an indicator of the origin of air masses (Jia et al., 2008). Air masses from the south were usually warmer and wetter than those from the north; thus, the specific humidity averaged over the BTH region was higher in Stage 1 (1.7 g kg⁻¹) than in Stage 2 (1.4 g kg⁻¹) (Fig. 7b). The evolution of PM_2.5 followed the trend of specific humidity well, with a high correlation coefficient of 0.86.

4.3 Contributions of each physical/chemical process to variations in PM_2.5 concentrations

Figure 8a1–a2 show the diurnal variations of PM_2.5 concentrations averaged over the BTH region during Stage 1 and Stage 2, respectively. The PM_2.5 concentration increased by 43.9 µg m⁻³ (from 136.5 µg m⁻³ at 00:00 LST to 180.4 µg m⁻³ at 23:00 LST) during the period of particulate accumulation (Stage 1), but it decreased by 41.5 µg m⁻³ during the period of particulate elimination (Stage 2).

The hourly PM_2.5 changes induced by each and all physical/chemical processes during Stage 1 and Stage 2 established using the IPR analysis method are shown in Fig. 8b1–b2. During both stages, the dominant sources of surface-layer PM_2.5 were EMIS and AERC, whereas the main sinks were TRAN, DIFF, and DRYD. The maximum positive contribution of EMIS could be found during the rush hours (07:00–08:00 and 16:00–19:00 LST; Fig. S3). The maximum negative contributions of TRAN and DIFF appeared at late night (01:00–05:00 LST) and at noon (11:00–14:00 LST), respectively.

Figure 8(a1–a2) Diurnal variations of PM_2.5 concentrations averaged over Beijing–Tianjin–Hebei during Stage 1 and Stage 2 (purple dotted lines). The colored bars represent different components. The 24 h change in PM_2.5 concentration (23:00 minus 00:00 LST) is also shown in the top-left corner of each panel. (b1–b2) The hourly PM_2.5 changes induced by each physical/chemical process using the IPR analysis method (colored bars). The purple dotted lines represent hourly PM_2.5 changes induced by all processes, also indicating the differences between current and previous-hour PM_2.5 concentrations. (c1–c2) Contributions of each physical/chemical process to 24 h PM_2.5 changes.

Download

To explain the reason for the 24 h PM_2.5 increase during Stage 1 and the 24 h PM_2.5 decrease during Stage 2 (Fig. 8a1–a2), we quantify the contributions of each physical/chemical process to 24 h PM_2.5 changes for both stages (Fig. 8c1–c2), which are calculated by integrating the hourly PM_2.5 changes induced by each process from 00:00 to 23:00 LST (Fig. 8b1–b2). In WRF-Chem, DRYD is intermingled with vertical diffusion, meaning that changes in the column burden during vertical mixing can be attributed to DRYD (Tao et al., 2015). Following Tao et al. (2015), we define vertical mixing (VMIX) as the sum of DIFF and DRYD. As shown in Fig. 8c1–c2, contributions of the AERC, TRAN, and VMIX processes to the 24 h PM_2.5 changes were +29.6 (+17.9) µg m⁻³, −71.8 (−103.6) µg m⁻³, and −177.3 (−221.6) µg m⁻³ for Stage 1 (Stage 2), respectively. Small differences were found for contributions from other processes between Stage 1 and Stage 2 (differences smaller than 5 µg m⁻³). Therefore, the PM_2.5 increase over the BTH region during the haze formation stage was mainly attributed to strong production by the aerosol chemistry process and weak removal by the advection and vertical mixing processes. On the contrary, during haze elimination stage (Stage 2), more aerosols in the BTH area were transported out of the BTH region, dispersed to the upper atmosphere or subsided to the ground. Furthermore, the dry cold air from the north decreased the specific humidity (as shown in Fig. 7b) in the BTH area, leading to weaker production of secondary aerosols by aerosol chemistry process.

5.1 Effects of aerosol radiative forcing on meteorological parameters and PM_2.5 concentrations

Figure 9 illustrates the impacts of aerosols on the downward shortwave radiative flux (SW) at the surface (BOT_SW) and in the atmosphere (ATM_SW), calculated by subtracting the model results of NoARE from those of CTL, during Stage 1, Stage 2, and the whole simulation period. Downward SW at the surface decreased strongly when ARE was considered, especially over high aerosol-loading regions during heavily polluted periods. Generally, the shortwave radiation fluxes at the surface averaged over BTH were reduced by 28 % (23.9 W m⁻²) in Stage 1, 18 % (16.6 W m⁻²) in Stage 2, and 23 % (19.9 W m⁻²) during the whole simulation period. Contrary to the significant negative effects at the surface, as a result of ARE, the downward SW fluxes in the atmosphere averaged over the BTH region were increased by 65 % (19.1 W m⁻²) in Stage 1, 37 % (10.8 W m⁻²) in Stage 2, and 51 % (14.7 W m⁻²) during the whole period.

Figure 9The differences in simulated all-sky radiative forcing (W m⁻²) between the CTL and NoARE cases (CTL minus NoARE) averaged over Stage 1, Stage 2, and the whole simulation period. “BOT_SW” and “ATM_SW” denote the downward shortwave radiative flux at the surface and in the atmosphere, respectively. The calculated differences in the simulated radiative forcing averaged over Beijing–Tianjin–Hebei for each stage are also shown at the bottom of each panel.

The impacts of ARE (including aerosol direct and indirect effects) on meteorological parameters and PM_2.5 concentrations are analyzed in Fig. 10. Because less SW could reach the ground, near-surface temperature decreased over BTH (Fig. 10a), especially during heavy pollution periods, and the largest decrease was up to 2 K. Meanwhile, the increased SW in the atmosphere warmed the upper air. As a result, a more stable atmosphere was expected. It is known that the atmospheric stability can be exactly characterized by the profile of equivalent potential temperature (EPT) (Bolton, 1980; Zhao et al., 2013; J. Yang et al., 2016). If the EPT rises with height, the atmosphere is stable. As shown in Fig. 10b, the EPT decreased in the lower atmosphere (below ∼1000 m) with the largest decrease of 3 K on 22 December, but it increased in the upper atmosphere (above ∼1200 m). The change in the EPT profile indicated that ARE could lead to a more stable atmosphere, which further weakened vertical movement in the BTH (Fig. 10c). As a result of ARE, the PBLH decreased and the relative humidity in the lower atmosphere increased (Fig. 10d). All of the changes in the meteorological variables were beneficial for PM_2.5 accumulation in the lower atmosphere (Fig. 10e). The daily maximum increase in the PM_2.5 concentration was 43.2 µgm⁻³ due to ARE. It was noticed that ARE had a negative impact on the near-surface PM_2.5 concentrations from 23 to 24 December, which could be explained by the fact that absorbing aerosols (i.e., BC) induced anomalous northeasterlies, and then the relatively clean air transported from the northeastern regions to the BTH region (Fig. S4).

Figure 10Time series of differences in (a) temperature (K), (b) equivalent potential temperature (K), (c) vertical wind speed (cm s⁻¹), (d) relative humidity (%), and (e) PM_2.5 concentration (µg m⁻³) between the CTL and NoARE cases (CTL minus NoARE) averaged over the Beijing–Tianjin–Hebei region. The purple and green lines denote the simulated PBLH in the CTL and NoARE cases, respectively. The black line represents the zero contour line.

Download

5.2 Influence mechanism of aerosol radiative effects

As variations in PM_2.5 concentrations are directly caused by physical and chemical processes (Zhu et al., 2015), the IPR method is then used to investigate the detailed influence mechanisms (i.e., the prominent physical or chemical processes responsible for the aerosol radiative impacts on haze episodes). Figure 11a–b show the diurnal variations of PM_2.5 concentrations in the NoARE and CTL cases averaged over the BTH region in Stage 1. A 24 h increase of 39.1 µg m⁻³ was simulated in the NoARE case. When aerosol radiative forcing was considered, the 24 h increase of PM_2.5 concentration was 43.9 µg m⁻³. The enhancement of 4.8 µg m⁻³ (12 %) induced by ARE could be mainly attributed to the contributions of the VMIX, TRAN, and AERC processes, as shown in Fig. 11c. The vertical mixing was strongly restrained by ARE; therefore, fewer particles diffused from the surface to the upper layer, resulting in the accumulation of PM_2.5 in a lower atmospheric boundary layer. The changes induced by ARE in contributions of the VMIX process exhibited positive values in the lower layers and negative values in the upper layers (Fig. S5a). Generally, the VMIX process contributed +22.5 µg m⁻³ to the enhancement in the 24 h PM_2.5 increase (+4.8 µg m⁻³) for Stage 1. The TRAN process, however, contributed −19.6 µg m⁻³. Constrained vertical mixing due to ARE could increase aerosol precursors and water vapor in the thin boundary layer to enhance the formation of secondary particles. Generally, the AERC process contributed +1.2 µg m⁻³. The positive contribution of AERC was mainly distributed over the highly polluted regions in the BTH area (Fig. S5b). Specifically, the average changes in the concentrations of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ during the daytime from 11:00 to 17:00 LST in Stage 1 were −0.5, +1.3, and +0.8 µg m⁻³, respectively. The decreased near-surface temperature caused by ARE may suppress the chemical formation of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. Generally, the total contribution of the VMIX, TRAN, and AERC processes to the change in the 24 h PM_2.5 increase caused by ARE was +4.1 µg m⁻³, and the restrained vertical mixing could be the primary reason for the near-surface PM_2.5 increase when aerosol radiative forcing was considered.

Figure 11Diurnal variations of the near-surface PM_2.5 concentrations in the (a) NoARE and (b) CTL simulations averaged over the Beijing–Tianjin–Hebei region during Stage 1 (purple dotted lines). The colored bars represent different components. The 24 h increase in PM_2.5 concentration (23:00 minus 00:00 LST) is also shown in the top-left corner of each panel. (c) Differences in hourly IPRs caused by aerosol radiative forcing (CTL minus NoARE). The numbers listed in panel (c) represent the contributions of each process to the change in the 24 h PM_2.5 increase caused by aerosol radiative forcing.

Download

Figure 12a shows the vertical profiles of the 24 h increases in the PM_2.5 concentrations (23:00 minus 00:00 LST) averaged over the BTH region during Stage 1 in the CTL and NoARE cases. Below ∼300 m (between L01 and L04), the 24 h increase simulated by CTL was larger than that in NoARE, which could be mainly explained by the fact that the positive contributions of VMIX exceeded the negative contributions of TRAN in the lower atmosphere when aerosol radiative effect was considered (Fig. 12b). However, in the upper layers (from 300 to 2000 m), aerosol radiative forcing weakened the 24 h PM_2.5 increase during Stage 1. When the aerosol radiative effect was considered, less particulate matter, precursors and water vapor were diffused from the surface to the upper layers; therefore, fewer particles were formed in the upper layers. Despite the positive contributions of TRAN, the net contributions of VMIX, TRAN, and AERC to PM_2.5 changes caused by ARE in the upper atmosphere were negative.

Figure 12(a) Vertical profiles of the 24 h increases in PM_2.5 concentrations (23:00 minus 00:00 LST) averaged over Beijing–Tianjin–Hebei during Stage 1 in the CTL and NoARE cases. (b) Vertical profiles of the differences in the 24 h PM_2.5 increases caused by the aerosol radiative effect (CTL minus NoARE, purple dotted line), and the contributions of each physical/chemical process (colored bars).

Download