Modeling of aerosol property evolution during winter haze episodes over a megacity cluster in northern China: roles of regional transport and heterogeneous reactions of SO₂

Regional transport and heterogeneous reactions have been shown to play crucial roles in haze formation over a megacity cluster centered on Beijing. In this study, the updated Nested Air Quality Prediction Model System (NAQPMS) and the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model were employed to investigate the evolution of aerosols – in terms of the number concentration, size distribution, and degree of aging – in Beijing during six haze episodes between 15 November and 15 December 2016, as part of the Air Pollution and Human Health–Beijing (APHH-Beijing) winter campaign of 2016. The model exhibited reasonable performance not only with respect to mass concentrations of PM_2.5 and its components in Beijing but also regarding the number concentration, size distribution, and degree of aging. We revealed that regional transport played a non-negligible role in haze episodes, with contributions of 14 %–31 % to the surface PM_2.5 mass concentration. The contribution of regional transport to secondary inorganic aerosols was larger than that to primary aerosols (30 %–63 % and 3 %–12 %, respectively). The chemical transformation of SO₂ along the transport pathway from source regions to Beijing was the major source of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. We also found that sulfate formed outside Beijing from SO₂ emitted in Beijing; this sulfate was then blown back to Beijing and considerably influenced haze formation. Along the transport pathway, aerosols underwent aging, which altered the mass ratio of the coating of black carbon to black carbon (R_BC) and the size distribution of number concentrations. During the episodes, the geometric mean diameter (GMD) increased from less than 100 nm at the initial site to approximately 120 nm at the final site (Beijing), and the R_BC increased from 2–4 to 4–8. During haze episodes with high humidity, the average contributions of gas and aqueous chemistry, heterogeneous chemistry, and primary emission to sulfate were comparable. However, their relative contributions varied with pollution levels. Primary emissions had the greatest impact under light to moderate pollution levels, whereas heterogeneous chemistry had a stronger effect under high pollution levels.

1 Introduction

In past decades, a megacity cluster in China that is centered on Beijing and includes 28 cities (272 500 km², with a population of 191.7 million people) has experienced frequent severe and persistent haze episodes (Zhao et al., 2013; Sun et al., 2014, 2016). PM_2.5 levels exceeding 500 µg m⁻³ have often been reported. The adverse effects of PM_2.5 on visibility, climate, and particularly human health have drawn widespread public attention (Hyslop, 2009; Chen et al., 2018; Yang et al., 2017a, b; Anderson et al., 2010). Although the PM_2.5 concentration in Beijing decreased by 35 % over a recent 5-year (2013–2017) span, benefiting from the implementation of the Atmospheric Pollution Prevention and Control Action Plan, the PM_2.5 level in 2017 still reached 58 µg m⁻³, which is 1.7 times the safe level of 35 µg m⁻³ recommended by the World Health Organization (http://www.bjepb.gov.cn/bjhrb/index/index.html, last access: 13 June 2019). Thus, understanding the mechanism behind haze episodes in this megacity cluster is an urgent task for policymakers.

Observations have revealed that haze episodes in this megacity cluster are mainly caused by the rapid formation of secondary inorganic aerosol (SIA) species, including sulfate, nitrate, and ammonium (Huang et al., 2014; Zheng et al., 2015; Han et al., 2016). The SIA mass fraction in PM_2.5 can be up to 55 % on severely polluted days, which is 2.5 times that observed on clear days (Ma et al., 2017). G. Tang et al. (2016) proposed that local chemical transformation associated with humidity dominated the rapid formation of SIA in Beijing. Yang et al. (2015) argued that local chemical conversion would not be fully able to explain the observed rapid formation of SIA over a short period of time. Using a ceilometer and in situ observation data, Zhu et al. (2016) and Ma et al. (2017) further proposed that regional transport was the major cause of the initial haze stage and that local chemistry, particularly heterogeneous chemistry, dominated the subsequent later rise in Beijing. This result differs from the findings of modeling studies (Timmermans et al., 2017; Li and Han, 2016; J. Li et al., 2017), in which regional transport was identified as the dominant factor during haze episodes in the megacity cluster. Comprehensive observations of the physicochemical properties (e.g., mixing state, number concentration, and size distribution) of aerosols can provide more insights into the accuracy of regional transport and chemistry assessment. Black carbon (BC) is usually more thickly coated by SIA and organic aerosols in transported and aged air masses than in fresh particles, as indicated by the higher fractal dimension (Y. Y. Wang et al., 2017), the larger coating fraction (ratio of variation in BC mass equivalent diameter to initial BC diameter, $\mathrm{\Delta }{D}_{\text{me}}/D_{\text{me},\mathrm{0}}$) (Peng et al., 2016), and the higher mass ratio of the coating of BC to BC (R_BC) (Wang et al., 2019). Massoli et al. (2015) and Y. Y. Wang et al. (2017) reported that the R_BC exceeded 10 in remote sites after BC had undergone long-term transport. This value was much higher than that reported in an urban area with high fresh-particle emissions, where the R_BC was generally less than 1.5 (Liu et al., 2017). The geometric mean diameter (GMD) of PM_2.5 also changed significantly due to the impact of regional transport. In haze episodes in Beijing, the GMD increased to 120 nm in regionally transported air masses, which was twice that observed under clean conditions (Ma et al., 2017). Thus, investigating the evolution of aerosol properties other than mass concentration during regional transport is useful for assessing the roles of regional and local chemistry. Such investigations are rarely conducted using the current three-dimensional chemical transport models. The current models generally account for only a part of the observed ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations (Y. Wang et al., 2014). Heterogeneous chemistry is considered critical to improving model performance (Zheng et al., 2015; Cheng et al., 2016; Li et al., 2018), and the treatment of heterogeneous chemistry is likely another source of modeling uncertainty. Studies found that uptake coefficients (γ) were dynamically related to aerosol liquid water and aerosol mixing states (coating thickness) (Riemer et al., 2009; Morgan et al., 2015). Although these parameters are obtained in some observation studies (Bian et al., 2014; Zhang et al., 2016), the parameterization of heterogeneous chemistry has rarely been linked to them in previous studies (Morgan et al., 2015; Zheng et al., 2015; Li et al., 2018).

From 15 November to 15 December 2016, a field campaign was carried out in Beijing within the framework of the UK–China Air Pollution and Human Health (APHH) project. Details regarding this project are given in Shi et al. (2018). Aerosol properties such as the size distribution, number concentration, and mixing states were simultaneously measured in China. APHH-Beijing aimed to explore the sources and processes affecting urban atmospheric pollution in Beijing. In this study, we used the Nested Air Quality Prediction Model System (NAQPMS) to simulate aerosol properties during the campaign period as a part of the APHH research. To improve model performance, the NAQPMS was updated by incorporating an advanced particle microphysics (APM) module that explicitly accounts for the microphysical process (Chen et al., 2014) and a new heterogeneous chemistry scheme (Li et al., 2018). The Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model was also employed to explicate the evolution of aerosol properties (e.g., mixing state, number concentration, and size distribution). Detailed analysis of the transport of precursors or secondary products, and heterogeneous reactions was mainly focused on sulfate, as recent studies have indicated that sulfate is a key driver for severe haze events (Huang et al., 2014; Zheng et al., 2015). Crucially, the effects of the regional transport and heterogeneous chemistry of SO₂ on aerosol properties were quantified. To our knowledge, this is the first study to distinguish the contributions of the transport of SIA alone and its precursors to PM_2.5 in Beijing as well as combining trajectories with the evolution of microphysical properties. We believe that this study will help to understand the causes of haze in this megacity cluster.

2 Model description and methodology

2.1 Model description

The Nested Air Quality Prediction Model System (NAQPMS) developed by the Institute of Atmospheric Physics, Chinese Academy of Sciences (IAP/CAS) is a three-dimensional Eulerian terrain-following chemical transport model. The Weather Research and Forecasting model (version 3.6.1; WRFv3.6.1), driven by Final Analysis (FNL) data from the National Centers for Environmental Prediction (NCEP), provides the meteorology field for the NAQPMS. The NAQPMS includes emissions, horizontal and vertical advection and diffusion, dry and wet deposition, and chemical (including gas, aqueous, and heterogeneous) reaction processes (Wang et al., 2001; Li et al., 2012, 2018). In the model, aqueous chemistry only occurs in cloud water, including the oxidation reactions of S(IV) by O₃, H₂O₂, methyl hydroperoxide, peroxyacetic acid, and oxidation catalysis of transition metal ions (Fe³⁺ and Mn²⁺). Heterogeneous chemistry reactions occur on aerosol's surface with the aqueous layer and are related to aerosol liquid water. For sulfate, assuming that the upper limit of ${\mathit{\gamma }}_{{\mathrm{SO}}_{\mathrm{2}}}$ does not exceed the uptake coefficient on dust surfaces, the upper limit of ${\mathit{\gamma }}_{{\mathrm{SO}}_{\mathrm{2}}}$ is $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ for aerosol water content (AWC) >  300 µg m⁻³ and its lower limit is $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ if the AWC < 10 µg m⁻³. More details can be found in Li et al. (2018). It also incorporates online source tagging, process analysis, an online WRF coupler, and other techniques (Wu et al., 2017; Z. Wang et al., 2014). The Carbon Bond Mechanism version Z (CBMZ) is used for gas-phase chemistry mechanisms. The thermodynamic model ISORROPIA (version 1.7) is used to calculate the composition and phase state of an ${\mathrm{NH}}_{\mathrm{4}}^{+}$–${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$–${\mathrm{NO}}_{\mathrm{3}}^{-}$–Cl⁻–Na⁺–H₂O inorganic aerosol system (Nenes et al., 1998). Six secondary organic aerosols (SOAs) are managed using a two-product module. Further details regarding the NAQPMS can be found in J. Li et al. (2013, 2014, 2017), and numerous subsequent papers have been published describing recent updates.

To accurately describe aerosol properties (e.g., number concentration, size distribution, and mixing states), an advanced multi-type, multicomponent, size-resolved advanced particle microphysics (APM) module is coupled to the NAQPMS (Chen et al., 2014). The APM module explicitly describes microphysical aerosol processes, including nucleation, condensation, evaporation, coagulation, thermodynamic equilibrium with local humidity, hygroscopic growth, and dry and wet deposition (Yu and Luo, 2009), and it has already been applied in the global GEOS-Chem model (Ma and Yu, 2014). In the updated NAQPMS, 40 sectional bins covering 0.0012–12 µm were used to represent secondary particle distribution (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and secondary organic aerosols) (Chen et al., 2014). The size distribution of BC and primary organic aerosol was represented using 28 section bins. Other primary particles such as dust and sea salt were represented using four bins. The coating of secondary species on primary particles (sea salt, BC, OC, and dust) was explicitly simulated using a scheme that dynamically calculates the aerosol aging time with an hourly resolution on the basis of aerosol microphysics. The mixing state is assumed to be semi-external. Primary particles coated with SIA or SOA are considered to be core–shell mixtures whereas nucleated secondary particle is internally mixed (Chen et al., 2014). The mass concentrations of coating species were also tracked in the model. Chen et al. (2017) employed the updated NAQPMS and revealed that the daytime aging time of BC in Beijing can be less than 2 h in winter. This is much less than the fixed aging timescale of 1.2 d that has been stipulated in previous studies (Liu et al., 2009) but is close to observed levels (2–4 h) (Peng et al., 2016). Li et al. (2018) further developed a heterogeneous chemical scheme based on mixing states to reproduce the chemical transformation of gaseous precursors on aerosol surfaces, which largely altered the sizes and hygroscopicity of particles. Heterogeneous chemistry includes the oxidation of S(IV) on the aqueous layer of aerosols and it is parameterized according to the scheme of Li et al. (2018). Comparison with long-term observations has proven that the updated NAQPMS can successfully estimate aerosol mass and the number concentration, size distribution, mixing states, and BC aging time in China (J. Li et al., 2017, 2018; Chen et al., 2014, 2017).

Distinguishing the contributions of the transport of SIA itself and its precursors to PM_2.5 is always difficult (Sun et al., 2014; J. Li et al., 2014, 2017; Ying et al., 2014). These contributions have generally been named regional transport in studies; this leads to ambiguity in regional transport. In this study, we divided the secondary species (e.g., SIA) in the ith receptor region into four parts: (1) SIA locally produced from the ith region locally emitted precursors (LC); (2) SIA chemically formed in other regions from the ith region locally emitted precursors (LTC); (3) SIA chemically formed along the transport pathway to the ith receptor region from precursors emitted in the jth source region (RTC); and (4) SIA produced in the jth region from precursors emitted in the jth source region (RLC).

Figure 1 (a) Simulation domains: “do2” represents the second domain, and “d03” represents the third domain. (b) Primary PM_2.5 emission rates of the innermost domain and locations of observation sites (black dots). (c) Tracer tagging regions that are described in Table 1.

An online tracer-tagging module in the NAQPMS was used to resolve the contributions from LC, LTC, RTC, and RLC. The module is capable of tracing both the emission regions of precursors and the formation regions of secondary aerosols. First, the mass contribution from the locations in which SIA was formed, referred to as C₂, was tagged. The mass contribution from precursors emitted in different locations, referred to as C₁, was then tagged. More technical details can be found in Li et al. (2014) and Wu et al. (2017). The following equation can be employed to calculate the degree of chemical conversion during transport (TC):

$$\begin{array}{}\text{(1)}& \text{TC}=\sum _{i=\mathrm{1}}^n(C_{\mathrm{1}i}-C_{\mathrm{2}i}\times {\text{CC}}_i),\end{array}$$

where i represents region, n is the total number of regions, and n is 10 in this study. C_1i refers to the absolute mass concentration transported to the receptor site, produced by precursors emitted in region i; C_2i refers to the absolute mass concentration formed in region i and transported to receptor site; CC_i refers to the local contribution ratio of precursors in region i; C_2i×CC_i refers to the absolute mass transported to the receptor site and generated at region i by the chemical conversion of precursors released at region i. For Beijing, i=1, refers to LC; i1, refers to RLC; is the mass concentration generated in all regions except I by the chemical conversion of precursors released at region I and finally transported to the receptor site. For Beijing, i=1, refers to LTC; i1, refers to RTC.

In this study, 10 regions, according to administrative division, are selected for source tagging (Fig. 1c), six of which – Chengde, Zhangjiakou, and Qinhuangdao (NHB); Beijing (BJ); Tianjin (LT); Hengshui, Xingtai, and Handan (SHB); Baoding and Shijiazhuang (WHB); and Tangshan, Langfang, and Cangzhou (EHB) – are parts of the Beijing–Tianjin–Hebei (BTH) region. Henan (HN), Shandong (SD), Shanxi (SX), and “Other regions” (OT) are areas outside of the BTH region.

Table 1 Source-tagging regions and primary PM_2.5 emissions from 15 November to 15 December 2016 in this study^a.

^a Regions are shown in Fig. 1c; ^b GDP unit in 2016 is Chinese Yuan (CNY) (http://www.tjcn.org/tjgb/, last access: 13 June 2019); ^c PM_2.5 emissions data are obtained from the 2016 Multi-resolution Emission Inventory for China (MEIC) at a $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ resolution.

2.2 Model configuration

The simulation was conducted from 10 November to 15 December 2016, and the first 5 d were set aside as a spin-up period. The three nested model domains are shown in Fig. 1a. The horizontal resolutions of the domains were 27, 9, and 3 km from the coarsest to the innermost domain. The first level of the NAQPMS was approximately 20 m in height, and there were about 17 layers under 2 km.

To quantitatively assess the contribution of primary emissions, traditional chemistry reactions (gas-phase and aqueous chemistry reactions), and heterogeneous formation of sulfate, three sensitivity simulations were conducted. The baseline scenario (Base) involved all heterogeneous reactions and considered primary sulfate emissions; its results were used for model verification and source apportionment analysis. Control 1 (C1) involved all heterogeneous reactions but did not consider primary sulfate emissions. Compared with Base, Control 2 (C2) excluded the heterogeneous reactions of SO₂. Base-C2 represents the effect of heterogeneous reactions on sulfate. Base-C1 represents the contribution of primary sulfate emissions.

The HYSPLIT model was used to analyze the trajectories of air masses (Draxler and Hess, 1998). The calculated trajectories are helpful for resolving the evolution of aerosol properties along the transport pathway by extracting the simulated results by the NAQPMS along trajectories. In this study, the same meteorology data (obtained hourly data of the third domain) used in the NAQPMS were employed to perform trajectory analysis; this avoided the errors caused by inconsistency between the two models (the NAQPMS and HYSPLIT).

2.3 Emission inventory

The anthropogenic emissions were obtained from the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ Multi–resolution Emission Inventory for China (MEIC), and the base year was 2016 for the BTH region (http://www.meicmodel.org/publications.html, last access: 12 November 2019). In addition, observation data collected at sites within the BTH region were used to update the MEIC on the basis of their latitude and longitude information. Biomass burning emissions were taken from the Fire Inventory from NCAR (National Center for Atmospheric Research) (Wiedinmyer et al., 2011). Primary sulfate was assumed to constitute 5 % of SO₂ emissions in the original MEIC inventory. Cao et al. (2014), Wang et al. (2009), Zheng et al. (2013), and Ma et al. (2015) discovered that primary sulfate comprised large amounts of primary PM_2.5 from industrial, power, and residential emissions in the main form of (NH₄)₂SO₄ via in situ measurement of source profiles. Thus, we took 40 %, 6 %, and 15 % of primary PM_2.5 from industrial, power, and residential emissions, respectively, as primary sulfate emissions in the inventory. Figure 1b displays the hourly primary PM_2.5 emission rate.

Figure 2 Comparison between the simulated and observed hourly concentrations of PM_2.5 for different sites. The black lines refer to observations (Obs.) and the red lines are simulation results (Sim.); the light blue shading denotes the six episodes identified; the green lines refer to 75 µg m⁻³, as a criterion for determining if pollution is present or not. FAC refers to the fraction of the simulation within 0.5 and 2 times the observation.

2.4 Observations

The surface meteorological parameters were obtained from the China Meteorological Administration, whereas the vertical profiles of meteorological parameters were obtained from the University of Wyoming (http://weather.uwyo.edu/upperair/sounding.html, last access: 13 June 2019). Observations of PM_2.5, SO₂, NO₂, and O₃ concentrations were obtained from the China National Environmental Monitoring Center (http://www.cnemc.cn/, last access: 13 June 2019). Aerosol components (including organic matter – OM, sulfate, nitrate, and ammonium) were measured in situ in Beijing using an Aerodyne high-resolution time-of-flight aerosol mass spectrometer. Details regarding this instrument can be found in Sun et al. (2015). A seven-wavelength Aethalometer (AE33, Magee Scientific Corp.) was used to measure BC in Beijing (Xie et al., 2019). The OC ∕ EC in aerosol was measured using a field semi-online OC ∕ EC analyzer from Sunset Laboratory Inc. (USA) with a PM_2.5 cyclone inlet at Tianjin and Langfang (Gao et al., 2016). Two ambient ion monitors (AIM; Model URG 9000D, URG Corporation) were used to measure hourly concentrations of water-soluble ions in PM_2.5, including ${\mathrm{NH}}_{\mathrm{4}}^{+}$, Na⁺, K⁺, Ca²⁺, Mg²⁺, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and Cl⁻ at TJ and Langfang (Gao et al., 2016). The particle number size distributions at ground level were obtained using a scanning mobility particle sizer (SMPS) with a time resolution of 5 min. Details regarding this instrumentation can be found in Du et al. (2017). All data in this study are presented in Beijing local time (UTC+8 h).

3 Model validation

3.1 PM_2.5 mass and number concentrations and degree of aging

The time series of simulated and observed PM_2.5 in different cities in the BTH region from 15 November to 15 December 2016, are illustrated in Fig. 2. During the study period, six regional haze episodes were identified, namely, 15–20 November (Ep1), 23–26 November (Ep2), 28–30 November (Ep3), 2–4 December (Ep4), 6–8 December (Ep5), and 10–12 December (Ep6). The PM_2.5 mass concentration frequently exceeded 200 µg m⁻³ and the average concentration reached 120 µg m⁻³ during episodes. Haze usually formed in several hours; for example, the increasing rate of PM_2.5 reached 200 µg m⁻³ h⁻¹ and lasted approximately 12 h in Tangshan. These observed haze patterns were generally reproduced by the NAQPMS. The correlation coefficient (R) between the observation and simulation in most cities was 0.6–0.8, and 60 %–80 % of simulation results were within a factor of 2 of the observation. The mean fractional bias (MFB) and mean fractional error (MFE) ranged from −0.07 to 0.7, meeting the criteria of a MFB ≤ 0.6 and a MFE ≤ 0.75 (Boylan and Russell, 2006). The simulation did, however, underestimate PM_2.5 in Beijing and Baoding during Ep2. This was caused by the failure of the mineral aerosol transport simulation. Compared with other cities in the cluster, Beijing and Baoding are closer to the Gobi Desert, a major dust source in East Asia; thus, they are more easily affected by dust storm transport. Pan et al. (2019) found a pronounced peak in the size distribution at 4–5 µm for Ep2 in Beijing. The concentration of Ca²⁺ was 7 times the campaign average (Fig. S1).

Figure 3 Comparison between the simulated (red) and observed (solid black) hourly components including sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), ammonia (${\mathrm{NH}}_{\mathrm{4}}^{+}$), black carbon (BC), and organic matter (OM) at (a) Beijing, (b) Tianjin, and (c) Langfang. Blue lines refer to sulfate produced by gas and aqueous chemistry reactions.

The aerosol components in Beijing, Langfang, and Baoding are compared in Fig. 3. In general, the simulation largely reproduced the variation in primary and secondary aerosols. In particular, the rapid increase in SIA during Ep1, Ep2, and Ep4 was captured by the simulation. Interestingly, the NAQPMS underestimated the sulfate concentration in Beijing during Ep2 and Ep4, but the nitrate and ammoniate concentrations during these two episodes were successfully reproduced. This was related to the transport of mineral dust (Ep2) and local emissions (Ep4). As discussed in the last paragraph, Beijing had high mineral loadings for Ep2, which provided a favorable medium for the chemical transformation of anthropogenic SO₂ into sulfate in the form of CaSO₄ or MgSO₄ (Z. Wang et al., 2018, 2017). Underestimation of the sulfate concentration in Ep4 may have been caused by local emissions in Beijing. As illustrated in Fig. 3, the simulation failed to reproduce the sharp increase in both sulfate and BC in Beijing during this episode. This is different from the case of Ep2, in which sulfate was underestimated but BC was favorably reproduced. Wang et al. (2009) and Ma et al. (2015) found that sulfate accounted for about 40 % and 6.6 % of primary PM_2.5 emissions from industry and power plants, respectively, which also emit a large amount of BC. This sharp increase in BC was a local-scale episode. In Langfang, a site located 50–60 km from Beijing, both the observed and simulated BC concentration increased slowly to 20 µg m⁻³, which is much less than that in Beijing (45 µg m⁻³). The monthly emissions employed in this study made it difficult to capture these short-term local-scale emission changes. The simulated SO₂ concentrations are compared with the observations in Fig. S2, and the normalized mean bias (NMB) values were less than 40 %.

Figure 4 (a) Particle size distribution in Beijing at ground level. (b) Comparison of the geometric mean diameter (GMD) for particles in the range of 16–600 nm between observations and simulated values for Beijing. The black solid line and the red solid line represent the observations (Obs.) and simulation (Sim.) results, respectively.

The number size distribution is critical for examining aerosol evolution during haze episodes (Ma et al., 2017). In this study, both the simulation and observations revealed a rapid increase in the GMD from 50 to approximately 120 nm during the initial stages of episodes in Beijing (Fig. 4). The observed mean number concentration of aerosols (dN∕dlog D_p) showed a unimodal distribution and was mainly concentrated in the Aitken mode (25–100 nm) and accumulation mode (100–1000 nm). The highest concentration was approximately 1.8×10⁴ cm⁻³ at a 100 nm diameter. These patterns were favorably reproduced by the simulation. The simulated number concentrations were underestimated in the 10–60 nm range by 20 %–30 % and overestimated in the 80–150 nm range by 20 %. This indicated that the model needs to be improved regarding its treatment of new particle formation and the volatility of primary organic aerosols.

Herein, the mass ratio of the coating of BC to BC (R_BC) was used as an indicator of the degree of aging – which is a measure that has been widely used in previous studies (Oshima et al., 2009; Collier et al., 2018). Figure 11 shows that the mean simulated R_BC values in Beijing were 4.5 and 5.0 for the entire study period and for during pollution episodes, respectively, which are extremely close to the observations (∼ 5.0 and 5.1) (Wang et al., 2019). The high performance of the model in terms of mass and number concentrations, compositions, and the degree of aging of aerosols gives us confidence for analyzing aerosol evolution during transport in the megacity cluster centered on Beijing.

3.2 Meteorology

The simulated wind direction and speed coincided with the observations for the haze episodes. In particular, the model captured low wind speeds, and the moments when the wind shifted direction were well reproduced (Fig. S3). Regarding relative humidity and temperature, WRF produced high R values (0.68–0.93) and low NMB values (−0.51 to 0.44; Table S1). In particular, the high relative humidity during Ep1 was well reproduced. Inversion layers were present during the initial stage of haze formation (Fig. S4). The height of the inversion layers varied among episodes. During Ep1 and Ep6, strongly elevated inversion layers were present between 1 and 2 km, whereas the inversion layers were close to the surface during other episodes. Temperature inversion is favorable for pollution accumulation, and the model reproduced this feature well. In summary, the high performance of the meteorological simulation gave us confidence for PM_2.5 simulation.

4 Results and discussion

4.1 Source apportionment of surface PM_2.5

The simulated spatial distribution of average surface PM_2.5 levels and the wind vector during the six haze episodes are shown in Fig. 5. In general, two types of patterns were observed. The first pattern corresponded to Ep1, Ep4, and Ep6 and reflected that a highly polluted belt with PM_2.5 over 200 µg m⁻³ extended from the southwest to the northeast along the Taihang mountain range. In the second pattern (Ep2, Ep3, and Ep5), the PM_2.5 level of 150–200 µg m⁻³ was concentrated in three northern cities (Beijing, Tianjin, and Tangshan). In the other cities, the PM_2.5 mass concentrations ranged from 75 to 115 µg m⁻³, indicating a light pollution level according to the Technical Regulation on Ambient Air Quality Index (on trial).

Figure 5 Spatial distribution of simulated average surface PM_2.5 (µg m⁻³) and wind (m s⁻¹) over the BTH area. (a) Average of the whole study period. (b)–(g) Episode average of Ep1–Ep6 (identified in the text). Solid circles represent observations and have the same color bar as the simulations.

Figure 6 The contribution of regional transport and local emissions to the average (a) total PM_2.5; (b) secondary inorganic aerosol (SIA); and (c) primary inorganic aerosol (PIA), including BC and primary inorganic PM_2.5, over the BTH area. The numbers above the pie charts represent average concentrations (µg m⁻³) of certain species in certain cities.

Figure 6 shows the contributions of regional transport and local emissions to average PM_2.5, primary inorganic aerosol (PIA), including BC and inorganic primary PM_2.5), and SIA levels in different cities during the study period. The contribution from local emissions was higher than that of regional transport to the PM_2.5 mass concentration in all cities, except Heng Shui, Cangzhou, Langfang, and Qinhuangdao; the magnitude of local emission contributions was 49 %–80 %. The principle reason for this was the accumulation of local PIA emissions. In most cities, 64 %–93 % of PIA originated from local emissions (Fig. 6c). In contrast to PIA, the SIA contribution was dominated by regional transport of emissions from other cities (50 %–87 %). Even the emissions of cities outside the city cluster (e.g., Henan, Shanxi, and Shandong) were transported to the megacity cluster, traveling 500–1000 km. In Beijing, the local contribution to total PM_2.5 and PIA was 74 % and 94 %, respectively, whereas regional transport from other cities was the major source of SIA, contributing 51 %. The difference in source apportionment between PIA and SIA was related to the emission of PIA and formation mechanisms of SIA. Regarding PIA, the inversion layer and weak winds during stable weather conditions prevented PIA transport and resulted in local-scale accumulation of PIA emissions. The regional transport provided sufficient time (1–3 d) and aerosol surface for chemical transformation of precursors to SIA (Li et al., 2015; J. Li et al., 2017). This also indicates that regional controls on precursors would be the most efficient way to decreasing the SIA concentration in this megacity cluster. Our results agree favorably with the observed impact of regional emission controls during Asia-Pacific Economic Cooperation China 2014. During this gathering, the SIA concentration in Beijing decreased to a greater degree than the PIA concentration due to regional controls (Sun et al., 2016).

Figure 7 (a) Source contribution of PM_2.5 in Beijing; the pie charts represent the average status of each episode. (b) The relative contribution of different regions to SIA, OM, and PIA in Beijing at the surface layer during each episode (shaded). The concentrations are also shown (blue line).

The source apportionment in haze episodes in Beijing is illustrated in Fig. 7. Regional transport contributed 14 %–31 % to the surface PM_2.5 mass concentration during the six episodes. The highest contribution of regional transport occurred in Ep1 and Ep5 (29 % and 31 % of the total PM_2.5, respectively). In Ep1 and Ep5, the contribution of the SIA originating from regional transport reached 53 % and 63 %, respectively. Interestingly, the regionally transported SIA had different source regions in Ep1 and Ep5. In Ep5, SX, WHB, and NHB were the dominant source regions, whereas the source regions for Ep1 were more diverse. This indicates the complexity of regional transport in this megacity cluster. Compared with the episodes in November 2015, the effects of regional transport of PM_2.5 and SIA mass concentrations were weaker in this study, which may be related to the weather system and emission controls in 2016 (J. Li et al., 2017). Therefore, more studies on regional transport should be conducted to help further understand regional haze formation mechanisms. In other episodes (Ep2, Ep3, Ep4, and Ep6), regional transport of surface PM_2.5, PIA, OM (organic matter; sum of primary organic aerosols, POA, and secondary organic aerosols ,SOA), and SIA mass concentrations were in the range of 14 %–23 %, 3 %–12 %, 3 %–14 %, and 30 %–51 %, respectively. Back trajectories and emissions source regions can be connected (Fig. S5). In Ep6, for example, the air mass mainly originated from Shandong, traveled through SHB and WHB, and finally arrived in Beijing. Furthermore, the height of the trajectory within WHB was low, meaning that the contribution of WHB should have been big, which agrees with results shown in Fig. 7b; additionally, source apportionment results showed that WHB contributed 24 % to SIA in Beijing during Ep6.

Figure 8 (a) Relative contribution of regionally transported SIA under different pollution levels in Beijing during the whole study period. (b) Variation of wind direction under different pollution levels in Beijing during the whole study period.

Figure 8 presents the relative contribution of regionally transported SIA and variation of wind direction under different pollution levels in Beijing. The source regions varied considerably under different pollution levels. When Beijing is controlled by strong northerly winds, NHB and SX are the main source regions, contributing up to 30 % and 19 % to SIA, respectively, resulting in clean conditions (SIA < 50 µg m⁻³). When Beijing is mainly affected by southerly winds (southeast, south, and southwest), WHB, EHB, and SD become the main source regions, contributing 27 %, 13 % and 15 % to SIA, respectively. Strong emissions of source regions lead to heavier pollution in Beijing. When Beijing is dominated by weak southeasterly winds, the contributions from distant regions like HN and SD increase. Continuous transport and accumulation lead to severe pollution (SIA > 150 µg m⁻³). This indicates that wider regional emission control is necessary to reduce severe pollution.

4.2 The impact of regional transport of sulfate and its precursors in Beijing

Quantifying the impact of regional transport of sulfate and its precursors is a crucial task. Sun et al. (2014) considered sulfate formed outside Beijing as regionally transported sulfate, and they estimated that its contribution reached 75 % during winter haze episodes. By tagging emission regions of precursors in models and ignoring where secondary aerosols were formed, J. Li et al. (2017) and Timmermans et al. (2017) estimated the contribution of transport to be 40 %–50 %. These estimated contributions of regional transport are different in physical meaning, which may confuse policymakers. In this study, we divided the sulfate mass concentration in Beijing into four parts, LC, LTC, RLC, and RTC as described in Sect. 2.1. The regional transport defined by Sun et al. (2014) was LTC + RLC + RTC, whereas in the studies by J. Li et al. (2017) and Timmermans et al. (2017), it was RLC + RTC. In this study, we employed RLC + RTC as regional transport.

Figure 9 (a) Sources of secondary sulfate in Beijing. (b) Variation of secondary sulfate sources with the surface sulfate concentration in Beijing for the whole study period. LC refers to sulfate locally produced from SO₂ emitted in Beijing; LTC refers to sulfate chemically formed in regions outside of Beijing from SO₂ emitted in Beijing; RTC is sulfate chemically formed along the transport pathway to Beijing from SO₂ emitted in source regions outside Beijing; RLC is sulfate produced in regions outside Beijing from locally emitted SO₂ and then transported to Beijing. PSO₂ refers to particulate sulfate. Sections of the pie charts where the values are not specified represent values of less than 3 %.

Figure 9a shows the contributions of LC, LTC, RLC, and RTC to the daily average sulfate concentration in Beijing during the study period. RTC and LC were the dominant sources of sulfate, contributing 71 %–89 % in total. The contributions of RTC ranged from 29 % in Ep6 to 59 % in Ep2, and contributions of LC were 30 %–42 %. RTC dominated the regional transport over the whole period, which indicates that chemical conversions along the transport pathway of SO₂ were critical to haze formation. Notably, the LTC contribution was comparable to that of LC in Ep3, Ep4, and Ep6. This suggests that the SO₂ emitted in Beijing was blown away and formed sulfate outside Beijing, and recirculation of air masses can be supported by HYSPLIT trajectories (Fig. S6). Take trajectories at 23:00 on 12 December (LST), for example, air masses were blown away from Beijing by southwesterly winds, moved through Chengde, Tianjin, and Langfang, and finally traveled back to Beijing. These formed sulfates may have been blown back to Beijing under specific weather conditions and were previously considered regional transport. The contribution of LTC also largely explains the difference in estimated regional transport contributions between Sun et al. (2014) and J. Li et al. (2017). In the present study, LTC + RLC + RTC accounted for 58 %–70 % of the sulfate concentration in the six episodes, which is relatively similar to the estimation (75 %) of Sun et al. (2014), which was based on the observed hourly increase rate of the local sulfate concentration.

In the initial and subsequent pollution stages, LC, LTC, and RTC showed different patterns in Beijing. In Ep1, local contributions dominated before the sulfate concentration increased rapidly (15 and 16 November). In particular, sulfate blown back to Beijing from its local emissions (LTC) made a larger contribution (35 %) than RTC (25 %). In the rapid rising phase of sulfate (17 and 18 November), the contribution of RTC increased from 25 % to 47 %. LC was also significant and its contribution increased considerably from 37 % to 41 %. These two parts (LC and RTC) explained the rapid formation of sulfate in Beijing. This suggests that the joint control of local and regional SO₂ emissions is essential to prevent the rapid formation of haze in this region, which is receiving considerable attention and eliciting widespread interest among researchers and policymakers (Sun et al., 2014; Ma et al., 2017; J. Li et al., 2017). This feature is also reflected in Fig. 9b. Under clear conditions (sulfate < 20 µg m⁻³), the local contributions (LC and LTC) were positively correlated with the sulfate mass concentration. In total, they contributed 40 %–60 % of the sulfate mass concentration. The ratio of LC to LTC was approximately 2 : 1. Under moderate sulfate levels (20 µg m⁻³ < sulfate < 35 µg m⁻³), the local contribution was lower – particularly the LTC – leading to a ratio of LC to LTC of approximately 8. The contribution of sulfate formed along the regional transport pathway (RTC) significantly increased from 40 % to 65 %. Under heavy pollution levels (> 35 µg m⁻³), the LC contributed up to 50 % due to extremely stable boundary layers. Our results are consistent with those of Ma et al. (2017), in which regional transport and local heterogeneous chemistry were qualitatively discovered to make high contributions to initial and subsequent pollution stages.

4.3 Evolution of aerosol properties in Beijing during haze episodes

Aerosol properties such as the particle size and the degree of aging can change dramatically on haze days due to fresh emissions, subsequent chemical conversions, and regional transport, which strongly affect regional radiation and climate (Cappa et al., 2012). As illustrated in Fig. 4b, the GMD of aerosols in Beijing increased remarkably to approximately 120 nm during the six haze episodes, compared with 50 nm under clean conditions. Two stages were identified: an initial rising stage and a sustained increase stage. In the initial stage, the GMD of aerosols increased by 50–60 nm over several hours, and the GMD then remained at 100–120 nm for several days in the subsequent elevated pollution stage. This GMD increase during the initial stage was mainly caused by the increase of accumulation-mode particles with diameters of 100–1000 nm and Aitken-mode particles (Fig. 10). Under clean conditions (SIA < 50 µg m⁻³), the average contributions of the three modes (nucleation, Aitken, and accumulation modes) to the number concentration were comparable, although the number of nucleation-mode particles decreased with SIA concentration. Under light to moderate pollution conditions (50 < SIA < 150 µg m⁻³), the proportion of accumulation-mode particles significantly increased from 35 % to 60 %, whereas the proportion of Aitken-mode particles slowly decreased. As discussed in previous sections, regional transport played a dominant role during the initial stage. This indicates that condensation, coagulation, and chemical transformation along the transport pathway increased the number of particles with a diameter in the range of 100–1000 nm. Finally, the contributions of Aitken-mode and accumulation-mode particles remained stable under heavy pollution conditions (SIA > 150 µg m⁻³).

Figure 10 Variation of number concentration fraction of particles with SIA in Beijing during the whole study period.

Figure 11 (a) Average and standard variation of the mass ratio of the coating of BC to BC (R_BC) during different episodes and pollution levels. (b) Diurnal cycles of R_BC under different pollution levels. (c) Temporal variation of R_BC during study period. “Ave” refers to the average over the whole study period.

Aging processes play a critical role in the growth of particles during haze episodes. According to observations, a significant coating of secondary components was found on BC in the study period (Wang et al., 2019). Figure 11 presents a time series of the simulated R_BC, which is a favorable indicator of the degree of aging (Oshima et al., 2009; Collier et al., 2018). A higher R_BC indicates that BC has undergone more ageing. In this study, the simulated R_BC was 2–10, with an average value of 4.5. Under polluted conditions, the R_BC was higher than under clean conditions, with an average value of 5.0. R_BC in Beijing even exceeded 10.0 during some extremely severe pollution events, which is close to observations of remote sites (J. Wang et al., 2017; Massoli et al., 2015) and aged particles (Cappa et al., 2012). Urban aerosols usually have a lower R_BC owing to fresh emissions, and the high R_BC values in this study indicate that Beijing aerosol particles were more aged during the haze episodes. On clean days, R_BC values ranged from 2 to 5, with an average of 2.8. This is similar to the R_BC of vehicle emissions (< 3) (Liu et al., 2017), indicating that Beijing is affected by local emission. Vehicle emissions contributed 70 % of BC in downtown Beijing in 2016 after strict controls on coal burning had been implemented (Kebin He, personal communication, 2016).

Figure 12 Variation of aerosol properties along transport pathways. Panels (a–f) refer to Ep1–Ep6. The red lines refer to 24 h backward trajectories at an altitude of 100 m. Aerosol properties include the geometric mean diameter (GMD; in nanometers – the red numbers beside the solid blocks), the mass ratio of the coating of BC to BC (R_BC; an indicator of the degree of aging – represented by the black numbers beside the solid blocks), the regional sources of BC (shown by the pie charts – red represents regional transport and gray is the local contribution). Solid black triangles denote the end points of back trajectories, referred to as T₀. Solid black circles denote points along trajectories per 6 h. T₋₆, T₋₁₂, T₋₁₈, and T₋₂₄ denote 6, 12, 18, and 24 h before arriving at the end point. End times of back trajectories are before pollution peaks at 21:00 on 18 November, 22:00 on 25 November, 16:00 on 29 November, 22:00 on 3 December, 00:00 on 8 December, and 22:00 on 11 December (LST), respectively. Sections of the pie charts where the values are not specified represent values of less than 6 %.

Figure 13 The contribution of different formation pathways for sulfate in Beijing. (a) The daily average. The blue line shows relative humidity in Beijing. The pie charts show the average contribution of different pathways during each episode. (b) The relationship between the sulfate concentration and different formation pathways of sulfate during Ep1. “Het” refers to the contribution of heterogenous chemistry of SO₂, and “Emis” refers to the contribution of primary sulfate. Sections of the pie charts where the values are not specified represent values of less than 3 %.

Figure 12 shows the evolution of R_BC, GMD, and regional sources of BC along the transport pathway from the source region to Beijing during the six haze episodes. Aerosol properties including number concentration along the transport pathway per 6 h are shown in Table S2. The transport pathway was calculated using the HYSPLIT model. The figure clearly shows that the aerosol properties changed considerably along the transport pathway. In Ep1, the GMD of aerosols was only 97 nm at the initial site of the 24 h back trajectory (T₋₂₄). At a larger transport distance, the diameters of aerosol particles were markedly increased to 128 nm in the middle (T₋₁₂) and 134 nm at the final site (T₀) of the back trajectory. R_BC increased from 3.6 at T₋₂₄ to 8.7 in Beijing (T₀) as BC was coated during the transport. This indicates that BC underwent considerable aging and increased in size while moving along the transport pathway. Similar characteristics were discovered for Ep3–Ep6. In Ep3, Ep4, Ep5, and Ep6, the GMD in Beijing (T₀) was 126, 117, 124, and 116 nm, respectively, compared with 96, 95, 99, and 111 nm mid-transport (T₋₁₂). R_BC also increased to between 4.6 and 7.6. Ep2 was an exception, as the GMD (106 nm) and R_BC (3.8) at the final end point (Beijing, T₀) were lower than those 6 h prior (T₋₆). Regional transport contributed 95 % of BC at T₋₆, whereas local emissions accounted for 87 % of BC at T₀. The number concentration was smaller at T₋₆ than at T₀. Therefore, we conclude that regional transport of aged aerosols led to a high GMD at T₋₆, and that the addition of locally emitted fresh air caused a high number concentration but low GMD at T₀. In clean areas, such as at T₋₂₄ in Ep5, the R_BC was higher than 10 and the GMD was considerably smaller.

4.4 Impact of heterogeneous chemistry on sulfate mass concentration

Current models generally account for a part of the observed ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in China (Y. Wang et al., 2014). Heterogeneous chemistry on aerosol surfaces under high relative humidity has been considered a potential missing source of sulfate formation (Cheng et al., 2016; Zheng et al., 2015; G. Li et al., 2017; M. J. Tang et al., 2016). Li et al. (2018) developed a simple parameterization of heterogeneous chemistry and discovered that SO₂ uptake on aerosols partly closed the gap between simulation results and observations. In their study, uptake coefficients were dependent on the aerosol core and shell species, shell thickness, and the amount of aerosol liquid water. Zheng et al. (2013) and Yang (2014) measured local source profiles, and they reported that primary sulfate from industry and power plants accounted for a large fraction of PIA. Researchers also found that household coal burning can emit a certain amount of sulfate (Zhang et al., 2018).

In this study, we examined the contributions of gas (SO₂+OH) and aqueous chemistry, heterogeneous chemistry, and primary sulfate emissions to the sulfate mass concentration in Beijing (Fig. 13). In Ep1, under high relative humidity, the contribution of heterogeneous chemistry was 33 %. Primary emissions mostly exerted an effect under light to moderate pollution levels (sulfate < 20 µg m⁻³), whereas heterogeneous chemistry played a crucial role under high pollution levels (sulfate > 30 µg m⁻³). The contributions of gas and aqueous chemistry were largely consistent under all pollution conditions (∼ 30 %). This indicates that high relative humidity and aerosol loading accelerated the SO₂ chemical transformation. Interestingly, the contribution of heterogeneous chemistry was markedly higher when the sulfate mass concentration exceeded the threshold of 20 µg m⁻³. Under high relative humidity and mass concentration conditions, a higher aerosol surface area resulting from hygroscopic growth and haze particles under high RH (60 %–80 %) with an aqueous shell due to phase transition provided a favorable media for heterogeneous reactions (Tie et al., 2017; Sun et al., 2018). The aforementioned threshold is relatively similar to that during the haze episodes in the winter of 2013 (Li et al., 2018). For policymakers, implementing measures to prevent the sulfate concentration from exceeding this threshold is essential. Such measures would be effective for avoiding extremely high sulfate levels. In other episodes, heterogeneous chemistry was depressed because of the low relative humidity (< 70 %). Gas and aqueous chemistry and primary emissions contributed 35 %–40 % and 58 %–61 %, respectively. It should be noted that failure of the model to simulate mineral dust led to an underestimation of the sulfate level in Ep2. The interaction between SO₂ and alkaline dust can contribute considerably to the sulfate concentration.

5 Conclusions

The contributions of regional transport to haze episodes over a megacity cluster centered on Beijing have been under debate in recent decades. Investigating the evolution of aerosol properties along the transport pathway may provide more information regarding how researchers can improve the accuracy of regional transport and chemistry impact assessments. To address one of the aims of the APHH 2016 winter campaign, we employed a Eulerian chemical transport model (NAQPMS) and a Lagrangian trajectory model (HYSPLIT) to assess the evolution of aerosols – in terms of the number concentration, size distribution, and degree of aging – in Beijing during six haze episodes between 15 November and 15 December 2016. The transport of sulfate and its precursors was also quantitatively investigated.

The results demonstrated that regional transport contributed 14 %–31 % to the surface PM_2.5 mass concentration in Beijing during the six episodes, with a monthly average contribution of 26 %. Regarding aerosol components, 30 %–62 % of the SIA in Beijing was regionally transported, whereas little PIA (< 10 %) was contributed by emissions from other regions. Source regions differed between episodes. During high pollution periods, WHB, SD, and EHB were the main source regions of SIA regionally transported to Beijing, whereas NHB and SX made greater contributions under clean and light pollution conditions. This indicates the complexity of regional transport in this megacity cluster.

The chemical transformation of SO₂ along the transport pathway from source regions outside of Beijing to Beijing (RTC) was the major source of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ regional transport and was more critical than the transport of sulfate formed in source regions outside Beijing (RLC). Compared with sulfate that was chemically transformed from SO₂ emitted in Beijing and then blown back to Beijing (LTC), the contribution of sulfate produced in Beijing from SO₂ emitted in Beijing (LC) was generally greater. However, RTC markedly increased in some episodes and explained the rapid formation of sulfate in Beijing. This suggests that the joint control of local and regional SO₂ emissions is essential to reduce the rapid formation of haze in this region.

Aerosols became considerably aged during transport in haze episodes, which altered the R_BC and the size distribution of number concentrations. During haze episodes, the GMD increased from less than 100 nm at the initial site to approximately 120 nm at the final site (Beijing), and the R_BC increased from 2–4 to 4–8. The number of accumulation-mode particles with a diameter between 100 and 1000 nm increased considerably more than the number of particles of different modes. The R_BC in Beijing during the episodes was higher than that in other urban regions (Collier et al., 2018). At high pollution levels, the R_BC was close to that observed in remote regions (J. Wang et al., 2017), which indicates that BC in Beijing under polluted conditions was more aged and thus more likely to affect radiation and circulation.

Contributions from different pathways to sulfate in Beijing were also examined. In episodes with high humidity (Ep1), the average contributions of gas and aqueous chemistry, heterogeneous chemistry, and primary sulfate were comparable; however, their relative contributions varied with pollution levels. Under light to moderate pollution levels, primary emissions mostly had an effect; however, under high pollution levels during Ep1, heterogeneous chemistry played a relatively crucial role. In other episodes (Ep2, Ep3, Ep4, Ep5, and Ep6), gas and aqueous chemistry and primary emissions contributed 35 %–40 % and 58 %–61 %, respectively.