Assessing the formation and evolution mechanisms of severe haze pollution in the Beijing–Tianjin–Hebei region using process analysis

Fine-particle pollution associated with haze threatens human health, especially in the North China Plain region, where extremely high PM2.5 concentrations are frequently observed during winter. In this study, the Weather Research and Forecasting with Chemistry (WRF-Chem) model coupled with an improved integrated process analysis scheme was used to investigate the formation and evolution mechanisms of a haze event over the Beijing–Tianjin–Hebei (BTH) region in December 2015; this included an examination of the contributions of local emissions and regional transport to the PM2.5 concentration in the BTH area, and the contributions of each detailed physical or chemical process to the variations in the PM2.5 concentration. The mechanisms influencing aerosol radiative forcing (including aerosol direct and indirect effects) were also examined by using process analysis. During the aerosol accumulation stage (16– 22 December, Stage 1), the near-surface PM2.5 concentration in the BTH region increased from 24.2 to 289.8 μgm−3, with the contributions of regional transport increasing from 12 % to 40 %, while the contribution of local emissions decreased from 59 % to 38 %. During the aerosol dispersion stage (23–27 December, Stage 2), the average concentration of PM2.5 was 107.9 μgm−3, which was contributed by local emissions (51 %) and regional transport (24 %). The 24 h change (23:00 minus 00:00 LST) in the near-surface PM2.5 concentration was +43.9 μgm−3 during Stage 1 and −41.5 μgm−3 during Stage 2. The contributions of aerosol chemistry, advection, and vertical mixing to the 24 h change were +29.6 (+17.9) μgm−3, −71.8 (−103.6) μgm−3, and −177.3 (−221.6) μgm−3 during Stage 1 (Stage 2), respectively. Small differences in the contributions of other processes were found between Stage 1 and Stage 2. Therefore, the PM2.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. When aerosol radiative feedback was considered, the 24 h PM2.5 increase was enhanced by 4.8 μgm−3 during Stage 1, which could be mainly attributed to the contributions of the vertical mixing process (+22.5 μgm−3), the advection process (−19.6 μgm−3), and the aerosol chemistry process (+1.2 μgm−3). The restrained vertical mixing was the priPublished by Copernicus Publications on behalf of the European Geosciences Union. 10846 L. Chen et al.: Formation and evolution mechanisms of severe haze pollution mary reason for the enhancement in the near-surface PM2.5 increase when aerosol radiative forcing was considered.


Introduction
Anthropogenic activities associated with rapid industrialization and urbanization have been leading to a sustained increase in the amounts of atmospheric pollutants, especially in quickly developing countries (IPCC, 2013). As one of the largest emission sources of aerosols and their precursors, China has been suffering from serious air pollution for years (Lei et al., 2011;Li et al., 2011;, with severe haze events frequently occurring in winter, especially over large urban agglomerations, such as the North China Plain (NCP) Gao et al., 2015), the Yangtze River Delta (YRD) , and the Sichuan Basin (SCB) Zhang et al., 2019). During severe haze events, the observed maximum hourly surface-layer PM 2.5 (fine particulate matter with an aerodynamic diameter of 2.5 µm or less) concentration can exceed 1000 µg m −3 (Z. Sun et al., 2016;, which can significantly influence visibility , the radiation budget (Steiner et al., 2013), atmospheric circulation , cloud properties (Unger et al., 2009), and human health (Hu et al., 2014;Guo et al., 2017).
Extensive studies have been carried out in recent years to analyze the formation mechanisms of haze episodes in China. Y.  used a synergy of ground-based observations, satellite, and lidar measurements to study a long-lasting, severe haze episode that occurred in eastern China in January 2013, and concluded that stagnant meteorological conditions, which can generally be characterized by a low wind speed, high relative humidity, intense inversion, and a low mixing layer height, were tightly associated with severe haze episodes. Based on National Center for Environmental Prediction (NCEP) reanalysis data, Shu et al. (2017) identified five typical synoptic patterns, and pointed out that each synoptic pattern exerted different impacts on particle pollution over the YRD. By analyzing the simulation results from a large ensemble climate model (MIROC5), K.  investigated the contributions of the anthropogenic influence to severe haze events that occurred over eastern China in January 2013 and December 2015, and found that anthropogenic forcing (i.e., increased emissions of greenhouse gases) could modify the atmospheric circulation pattern, and that these human-induced circulation changes were conducive to the occurrence of severe haze events. B.  used a global 3-D chemical transport model (GEOS-Chem) to quantify the local source contributions to wintertime surface-layer PM 2.5 concentrations over North China from 2013 to 2015, and reported that emissions from residential and industrial sources and transportation contributed most to the high concentrations of atmospheric aerosols in Beijing. Many studies have also reported that the regional transport of aerosols plays an important role in haze episodes (Z. Jiang et al., 2015;N. Li et al., 2018). Z.  reported that the "cross-city clusters transport" outside BTH (Beijing-Tianjin-Hebei) and transport among cities inside the BTH region contributed 20 %-35 % and 26 %-35 % of PM 2.5 concentrations over BTH, respectively. Secondary aerosol formation and their hygroscopic growth were also confirmed to be a large contributor to severe haze episodes (R. J. Han et al., 2015;. The conversion of SO 2 to SO 2− 4 was strongly associated with high relative humidity, and NO − 3 was found to be produced mainly by photochemical and heterogeneous reactions . It is well known that aerosols can scatter and absorb solar radiation to alter the radiative balance of the atmosphere and surface (direct radiative effect), and can serve as cloud condensation nuclei or ice nuclei to affect cloud properties (indirect radiative effect) (Twomey, 1974). These impacts are coupled with atmospheric dynamics to produce a chain of interactions with a large range of meteorological variables that influence both weather and climate (Ramanathan et al., 2001;Huang et al., 2006;Li et al., 2017a;Yang et al., 2017), which will further induce feedbacks on aerosol production, accumulation, and even severe haze pollution (Petaja et al., 2016;Li et al., 2017b;Zhao et al., 2017;Gao et al., 2018;Lou et al., 2019). Based on multi-year measurements (from 2010 to 2016), Huang et al. (2018) found that aerosol radiative effects led to a significant heating in the upper planetary boundary layer (PBL) and a substantial dimming at the surface over North China. This is because high concentrations of light-absorbing aerosols were observed, and the aerosolmeteorology interactions depressed the development of the PBL, and, in turn, aggravated the haze pollution (Su et al., 2018). The light-absorbing aerosols can also amplify haze in the NCP region by decreasing East Asian winter monsoon wind speeds via ocean and cloud feedbacks (Lou et al., 2019). Using the WRF-Chem model, Gao et al. (2015) analyzed the feedbacks between aerosols and meteorological fields over the NCP in January 2013, and found that aerosols caused a significant negative (positive) radiative forcing at the surface (in the atmosphere), resulting in a lower surfacelayer wind speed and lower PBL height (PBLH). The average surface-layer PM 2.5 concentration increased by 10-50 µg m −3 as a result of the more stable atmosphere. By analyzing the observations from a comprehensive field experiment and simulation results from WRF-Chem model, Q.  concluded that the decreased PBLH associated with increased aerosol concentrations could enhance surfacelayer relative humidity by weakening the vertical transport of water vapor, and that the increased relative humidity at the surface accelerated the formation of secondary particulate matter via heterogeneous reactions, leading to an increase of the PM 2.5 concentration by 63 µg m −3 averaged over the NCP from 15 to 21 December 2016.
All of the studies discussed above revealed that the formation of haze episodes was caused by the synergy impacts of local emissions, regional transport, meteorological conditions, and chemical production. Nevertheless, only the net combined effects on the concentrations of pollutants were provided, without the capability to understand and isolate the atmospheric physical and chemical processes involved. The quantitative assessment of the contributions from each detailed physical/chemical process (e.g., vertical mixing process, advection process, emission source process, aerosol chemistry process and cloud chemistry process) is necessary to fully understand the formation and evolution mechanisms of haze episodes (Gonçalves et al., 2009;Xing et al., 2017;Kang et al., 2019). Furthermore, although many previous studies have identified the positive feedback effects of aerosol radiative forcing on particulate accumulation, the detailed influence mechanisms of the forcing-response relationship at each process chain remain largely elusive (i.e., the prominent physical or chemical processes responsible for the aerosol radiative impacts on haze episodes). Since 2013, substantial efforts have been made to improve air quality in China, including emission reduction and energy transition. However, haze events have continued to frequently occur all over the country. For example, a severe, long-lasting, and wide-ranging haze episode was observed in December 2015 over central and eastern China, with the regional average PM 2.5 concentration exceeding 150 µg m −3 . In the BTH region, a red alert for haze (the most serious level) was issued for the period from 20 to 22 December 2015, with the maximum hourly PM 2.5 concentration exceeding 1000 µg m −3 . The formation and evolution mechanisms, and the aerosol radiative feedbacks of this severe haze episode have not yet been fully estimated.
In this study, we develop an improved online integrated process rate (IPR) analysis scheme (i.e., process analysis) in the fully coupled online Weather Research and Forecasting with Chemistry (WRF-Chem) model, to investigate the formation and evolution mechanisms of the severe haze episode that occurred over the NCP from 16 to 29 December 2015. Sensitivity experiments are conducted to examine the contributions of local emissions and regional transport to the PM 2.5 concentrations during the haze episode, while IPR analysis is used to quantify the contributions of each detailed physical/chemical process to the variations in the PM 2.5 concentrations. The effects of aerosol radiative forcing, including the direct and indirect effects, on meteorological parameters and PM 2.5 levels during the haze episode are also quantified, with a special focus on the detailed influence mechanism. We hope that the results from this study may provide a better understanding of the formation mechanisms for severe haze events, and help policy makers design and carry out targeted measures to improve air quality over North China. This paper is arranged as follows. The model configuration, integrated process rate (IPR) analysis (i.e., process analysis), numerical experiments, and observations are presented in Sect. 2. Model evaluation is conducted in Sect. 3. The formation and evolution mechanisms of the haze episode are investigated in Sect. 4. Section 5 provides the impacts of aerosol radiative forcing. Summaries and discussions are presented in Sect. 6.

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 (westeast) ×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.
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 1 • ×1 • . Four-dimensional data assimilation (FDDA) with the nudging coefficient of 3.0 × 10 −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;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. . 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 2 (sulfur dioxide), NO x (nitrogen oxides), CO (carbon monoxide), CO 2 (carbon dioxide), NMVOC (non-methane volatile organic compounds), NH 3 (ammo-  nia), BC (black carbon), OC (organic carbon), PM 2.5 , and PM 10 . 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, MO-SAIC (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 (SO 2− 4 ), nitrate (NO − 3 ), ammonium (NH + 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.

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 3 ) 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 3 concentrations. Gao et al. (2018) investigated the impacts of BC-PBL interactions on O 3 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 3 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 . 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. 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.

Numerical experiments
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. . 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 ana- Table 2. Experimental design. "Y" represents yes, and "N" represents no.

Case
Anthropogenic lyzing 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.

Model evaluation
Accurate representations of observed meteorological fields and pollutant concentrations provide foundations for haze analysis with the WRF-Chem model. Detailed comparisons between observed and simulated meteorological parameters (T 2 , RH 2 , WS 10 , WD 10 , PBLH, and SWDOWN) and pollutant concentrations (PM 2.5 , BC, OC, SO 2− 4 , NO − 3 , and NH + 4 ) are presented in this section. 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 2 , RH 2 , WS 10 , and WD 10 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 2 between the observations and simulations (Qian et al., 2016). This negative bias of RH 2 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.

Meteorological parameters
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 2 : temperature at 2 m (K); RH 2 : relative humidity at 2 m (%); WS 10 : wind speed at 10 m (m s −1 ); WD 10 : wind direction at 10 m ( • ). 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.

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 −3 was observed. This negative bias is also simulated by previous studies . The possible reasons for the underestimation are as follows: (1) the bias in simulated meteorological conditions (e.g., underestimated RH 2 and overestimated WS 10 ); (2) the missing mechanisms of some gas-aerosol phase partitioning and heterogeneous reactions which may produce secondary inorganic aerosol (X. 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 5 compares the simulated and observed surfacelayer concentrations of BC, OC, SO 2− 4 , NO − 3 , and NH + 4 in Beijing and Shijiazhuang averaged from 16 to 29 December 2015. The WRF-Chem model underestimates the con- centrations of SO 2− 4 , NH + 4 , and OC in Beijing (Shijiazhuang) by 19 % (40 %), 14 % (9 %), and 21 % (41 %), respectively, but overestimates the NO − 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 = +10 %) and Shijiazhuang (NMB = −24 %). For OC, the underestimation may result from the lack of SOA in the MOSAIC aerosol module (Qiu et al., 2017). Missing some SO 2 gas-phase and aqueous-phase oxidation mechanisms, as well as heterogeneous chemistry may explain the underestimation of SO 2− 4 . It is noted that similar biases of aerosol components were also reported by other WRF-Chem studies (B. Qiu et al., 2017).

Formation and evolution mechanisms of the haze episode
In this section, we first reproduce the evolution of the severe haze episode, and then investigate the formation and evolution mechanisms, including examining contributions of local emissions and regional transport to the PM 2.5 concentration in the BTH region, and the contributions of each detailed physical/chemical process to the variations in the PM 2.5 concentration. 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 −1 , 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 −3 . 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 −3 (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 highand 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. 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.

Spatial-temporal evolutions of surface-layer PM 2.5 concentrations
In Stage 1, the daily mean PM 2.5 concentrations averaged over the BTH region increased from 24.2 to 289.8 µg m −3 , and the average PM 2.5 concentration was 145.6 µg m −3 (Fig. 7a), which is close to the "heavily polluted" air quality threshold value (PM 2.5 24 h average concentration > 150 µg m −3 ). The WS 10 was low (Fig. 7b), especially during the heavily pollution period (20-22 December), and the mean wind speed was 2.3 m s −1 , which is less than 3.2 m s −1 (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 −3 , which still exceeded the Grade II standard (75 µg m −3 ) defined by the National Ambient Air Quality Standards of China.

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 Sun et al., 2014;Gu and Liao, 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)(21)(22)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 −3 ) 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 −1 ) than in Stage 2 (1.4 g kg −1 ) (Fig. 7b). The evolution of PM 2.5 followed the trend of specific humidity well, with a high correlation coefficient of 0.86.

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 −3 (from 136.5 µg m −3 at 00:00 LST to 180.4 µg m −3 at 23:00 LST) during the period of particulate accumulation (Stage 1), but it decreased by 41.5 µg m −3 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.
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 . 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 −3 , −71.8 (−103.6) µg m −3 , and −177.3 (−221.6) µg m −3 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 −3 ). 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.

Aerosol radiative effects (ARE) on the haze episode
Previous studies have demonstrated that aerosol radiative forcing could increase the near-surface PM 2.5 concentrations by about 12 %-29 % (Gao et al., , 2016Qiu et al., 2017;Zhou et al., 2018). However, the detailed influence mechanisms (i.e., the prominent physical or chemical process responsible for the aerosol radiative impacts on PM 2.5 concentrations) are still unclear. In this section, we examine the effects of aerosol radiative forcing on meteorological parameters and PM 2.5 levels during the haze episode, with a special focus on the detailed influence mechanism using IPR analysis. 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 −2 ) in Stage 1, 18 % (16.6 W m −2 ) in Stage 2, and 23 % (19.9 W m −2 ) 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 −2 ) in Stage 1, 37 % (10.8 W m −2 ) in Stage 2, and 51 % (14.7 W m −2 ) during the whole period. 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;. 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 −3 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). 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.

Influence mechanism of aerosol radiative effects
As variations in PM 2.5 concentrations are directly caused by physical and chemical processes , 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 −3 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 −3 . The enhancement of 4.8 µg m −3 (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 −3 to the enhancement in the 24 h PM 2.5 increase (+4.8 µg m −3 ) for Stage 1. The TRAN process, however, contributed −19.6 µg m −3 . 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 −3 . 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 SO 2− 4 , NO − 3 , and NH + 4 during the daytime from 11:00 to 17:00 LST in Stage 1 were −0.5, +1.3, and +0.8 µg m −3 , respectively. The decreased near-surface temperature caused by ARE may suppress the chemical formation of SO 2− 4 . 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 −3 , 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 9. The differences in simulated all-sky radiative forcing (W m −2 ) 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. 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.

Conclusions and discussions
In this study, an online coupled mesoscale meteorologychemistry model (WRF-Chem) with an improved integrated process rate (IPR) analysis (i.e., process analysis) scheme was applied to investigate the formation and evolution mechanisms of a severe haze episode that occurred in the BTH region from 16 to 29 December 2015. Sensitivity experiments were conducted to examine the contributions of local emissions and regional transport to the PM 2.5 concentrations during the haze event, while IPR analysis was used to quantify the contributions of each physical/chemical process to the variation in PM 2.5 concentration. The impacts of aerosol radiative forcing (including direct and indirect effects) were also quantified, with a special focus on the detailed influence mechanism (i.e., prominent process responsible for the aerosol radiative impacts on the haze event). An integrated comparison between observations and simulations demonstrated good performance for both meteorological and chemical variables, indicating that the WRF-Chem model has the capability to reproduce the haze episode.
Spatial-temporal evolutions of the near-surface PM 2.5 concentration, and the contributions of local emissions and regional transport to the severe haze event in BTH, were first analyzed. During the aerosol accumulation stage (16-22 December, Stage 1), the daily PM 2.5 concentration in the BTH region experienced a consistent increase, with the mean value of 145.6 µg m −3 . The contributions of local emissions and regional transport to the PM 2.5 concentration were comparable (49 % and 32 %, respectively), meaning that the combined effect resulted in the high PM 2.5 concentration in the BTH area. During the aerosol dispersion stage (23-27 December, Stage 2), the average PM 2.5 concentration in BTH was 107.9 µg m −3 . The contributions of local emissions and regional transport were 51 % and 24 %, respectively. Therefore, the relatively high PM 2.5 concentration during Stage 2 was principally caused by local emissions. Over the period Figure 10. Time series of differences in (a) temperature (K), (b) equivalent potential temperature (K), (c) vertical wind speed (cm s −1 ), (d) relative humidity (%), and (e) PM 2.5 concentration (µg m −3 ) 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. from 28 to 29 December (Stage 3), another haze event was formed and developed.
IPR analysis was then used to explain the reason for the PM 2.5 increase during Stage 1 and the decrease during Stage 2, by quantifying the contributions of each physical/chemical process to variations in PM 2.5 concentration. During both stages, the dominant sources were emissions (EMIS) and aerosol chemistry (AERC), whereas the main sinks were turbulent diffusion (DIFF), advection (TRAN), and dry deposition (DRYD). The PM 2.5 concentration increased by 43.9 µg m −3 (23:00 minus 00:00 LST) during Stage 1, but it decreased by 41.5 µg m −3 during Stage 2. The contributions of AERC, TRAN, and VMIX (vertical mixing, the sum of DRYD and DIFF) to the 24 h PM 2.5 changes were +29.6 (+17.9) µg m −3 , −71.8 (−103.6) µg m −3 and −177.3 (−221.6) µg m −3 for Stage 1 (Stage 2), respectively. Small differences in the contributions from other processes were found between Stage 1 and Stage 2. Therefore, the PM 2.5 increase over the BTH region during the haze formation stage was attributed to strong production by aerosol chemistry process and weak removal by advection and vertical mixing processes.
When aerosol radiative forcing was considered, the equivalent potential temperature decreased in the lower layers but increased in the upper layers, leading to a more stable atmosphere. Meanwhile, the decreased PBLH and increased relative humidity were also beneficial for PM 2.5 accumulation. The daily maximum increase of the near-surface PM 2.5 concentration in the BTH region was 43.2 µg m −3 . The IPR method was also used to investigate the detailed influence mechanism of aerosol radiative effects. When aerosol radiative feedback was considered, the 24 h PM 2.5 increase was enhanced by 4.8 µg m −3 (12 %) during Stage 1, which could be mainly attributed to the contributions of VMIX (+22.5 µg m −3 ), TRAN (−19.6 µg m −3 ), and AERC (+1.2 µg m −3 ). The restrained vertical mixing could be the primary reason for near-surface PM 2.5 increase when aerosol radiative forcing was considered.
There are some limitations to this work. The uncertainty of the MIX anthropogenic emission inventory, the lack of secondary organic aerosols, and the missing mechanisms of some heterogeneous reactions may result in large uncertainties in the final simulation results, especially the predicted aerosol chemical compositions, such as SO 2− 4 , NO − 3 , and NH + 4 . The biases in simulated concentrations of SO 2− 4 , NO − 3 , and NH + 4 may have impacts on the contributions of the AERC and CLDC processes to the air pollution variation. Uncertainties should be quantitatively analyzed in future studies. Furthermore, conclusions draw from a case study in the BTH region cannot represent a full view of the underlying mechanisms of haze formation and elimination. A better understanding will be attained by conducting multiplecase simulations in the future. Furthermore, an anomalous northeasterly induced by absorbing aerosols was observed, leading to a decrease in the near-surface PM 2.5 concentrations from 23 to 24 December 2015 in the BTH area, which was different from previous studies that reported that lightabsorbing aerosols could worsen air quality Huang et al., 2018;Gao et al., 2018). More experiments should be designed in the future to examine the changes in atmospheric thermal and atmospheric dynamic caused by absorbing aerosol radiative forcing and their impacts on haze episodes. As Zheng et al. (2018) pointed out, the PM 2.5 concentration in China has been decreasing in recent years; however, this decrease in fine particulate matter could stimulate ozone production (K. Zhu et al., 2019). Multipollutant mixture may be a hot topic in the future, and IPR analysis could be a useful method to provide a quantitative analysis of the formation mechanism of complex air pollution, including figuring out the major physical/chemical process behind these events. Meanwhile, significant differences between model predictions (e.g., O 3 and PM 2.5 ) are found among current multi-scale air quality models (L. , even though the same inputs are used. These different performances can be associated with the differences in model formulations, including parameterizations and numerical methods (Carmichael et al., 2008). In order to acquire a quantitative attribution of the cause of differences between simulation results, a process analysis method should be developed and implemented in these models, and the use of IPR analysis would make it easier to draw conclusions about the fundamental problems that cause the differences between model predictions.
Data availability. Observational datasets and simulation results are available upon request from the corresponding author (hongliao@nuist.edu.cn).