Supplementary material for the source apportionment of primary PM2.5 in an aerosol pollution event over Beijing-Tianjin-Hebei region using WRF-Chem, China

Key Laboratory of Meteorological Disaster, Ministry of Education (KLME), Joint International Research Laboratory of Climate and Environment Change (ILCEC), Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters, Key Laboratory for Aerosol-Cloud-Precipitation of China Meteorological Administration, Nanjing University of Information Science & Technology, Nanjing, China Qianshan Meteorological Bureau, Qianshan Country, Anqing City, Anhui Province, China State Key Laboratory of Atmospheric Boundary Layer Physics and Atmospheric Chemistry (LAPC), Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China


INTRODUCTION
In recent years, severe aerosol pollution events have occurred frequently in the Beijing-Tianjin-Hebei (BTH) region, impacting residents' health, life, and productivity (Pope et al., 2003;Parrish and Zhu, 2009;Brook and Rajagopanlan, 2010;Ding et al., 2013;Huang et al., 2014;Zhang et al., 2015).Researchers within China and around the world have investigated the distribution of aerosol concentrations and the causes of such serious pollution in this region (Okajima et al., 2001;Zhao et al., 2009;He et al., 2011;Liu et al., 2013;Li et al., 2014).
Atmospheric chemistry models can be used in the quantitative analysis of regional pollution source distribution and the characterization of local pollution.Such information can be provided to the government to shape policies that reduce and control pollution.Sensitivity analysis tests are widely used to assess the roles of various pollution sources (An et al., 2007;Streets et al., 2007).For example, Streets et al. (2007) used Community Multiscale Air Quality (CMAQ) to point out that about 34% of Beijing's fine particulate matter (PM 2.5 ) pollution originated outside Beijing.Analytical methods developed using online aerosol source apportionment in atmospheric chemistry models can reduce calculations while providing more accurate quantitative estimates of the contribution of source regions to overall pollution levels within a target region.Li et al. (2008) and Wu et al. (2011) used the Nested Air Quality Prediction Modeling System (NAQPMS) to quantify the contributions of Beijing and its surrounding areas to transportation-related pollutants within the city and determined that the surrounding regions were the greatest regional contributors to such pollution in Beijing.Long-distance transport of pollutants also significantly affects concentrations (Kevin et al., 1997;Wang et al., 2008;Liu et al., 2010).Wang et al. (2014) found that pollutants traveling across the urban agglomeration of the BTH region from outside the region contributed just as much to overall pollution levels as did local pollutants, and during the period of pollution case, the influence of cross-border transport was more significant.
Aerosols in the atmosphere can be classified as primary or secondary aerosols.Primary aerosols are directly released by natural and anthropogenic emissions, while secondary aerosols are produced by chemical reactions occurred in the atmosphere.Secondary aerosols are contingent upon the presence of primary aerosols and the relevant gaseous precursors in the atmosphere.Because of this relationship and the differences of source regions and distribution of primary and secondary aerosol, the most effective way to reduce aerosol pollution is to reduce primary emissions (Fan et al., 2005).In addition, some researches has also indicated that meteorological conditions and synoptic condition, including low wind speeds, a temperature inversion structure, low boundary layer height, an equalizing pressure field, and the area being behind high-pressure conditions, are the main factors of high pollutant concentrations in this area (Ji et al., 2012;Wang et al., 2012;Zhang et al., 2013;Wang et al., 2014).Moreover, the special terrain of the region surrounded on three sides by tall mountains in the BTH region is also a major contributor (Huang et al., 2012;Wang et al., 2014;Wu et al., 2014).So far, studies on heavy pollution in the BTH region have mainly focused on the polluted area along with surrounding areas, and little research has examined the connections between two different polluted areas.In this study, the Weather Research and Forecasting model coupled to Chemistry (WRF-Chem model) was used to simulate a severe aerosol pollution episode that occurred in the BTH region in the fall of 2015.We analyzed the causes of this pollution and quantified the contributions of primaryaerosols from all the geographical sources regions in model domain by the primary PM 2.5 source apportionment method.Understanding the contributions of different source regions to pollution within the BTH will provide a valid scientific basis for conservation-oriented policies and governance of these source regions.

Model Description and Modification
The WRF-Chem model is an online mesoscale atmospheric dynamics-chemical coupling model developed by the U.S. National Oceanic and Atmospheric Administration (NOAA) and other collaborators.This model system comprises meteorological (Skamarock et al., 2008) and air quality components (Grell et al., 2005) with relatively comprehensive chemical and physical processes.An advantage of this model is that all physical quantities in the chemical and meteorological modules and the flux transmission scheme are fully coupled online model, which enables air quality to simulate at the same time as the meteorological model runs.
We used version 3.4 of the WRF-Chem model and developed and integrated a primary PM 2.5 tagging method into the model in order to identify the primary PM 2.5 contributions of different source regions to BTH during this episode.The tagging method is a mass-balance technique similar to particulate source apportionment technology (PSAT) that can be applied for this purpose to any region (Ying and Kleeman, 2006;Wagstrom et al., 2008).It tracked and calculated the concentrations of PM 2.5 originating in different source regions (the independent variable) within the larger simulated region.The tracers in this tagging method account for the effects of physical processes (e.g., advection, vertical mixing, convection), as the normal simulation does, but do not disrupt the normal calculations.This case study was limited in that we only tagged primary PM 2.5 pollution, and the PM 2.5 source-apportionment is not complete, the secondary aerosols were not included.
In this method, the simulation domain was divided into 17 source regions (Fig. 1(a)).And the primary PM 2.5 identified with each source region was tracked as an independent variable.For each time step, the concentration of primary PM 2.5 with identification by i is determined using the following equation: where C i (t) is the primary PM 2.5 concentration with identification of i at the time of t, ∆t is the time step, ∆C i is the change in concentration of primary PM 2.5 at this time step.∆C i is determined by the following equation: where ∆Chem i is the concentration changes caused by chemical reactions.In aerosol chemistry, primary PM 2.5 doesn't participate in the chemical calculation which means the ∆Chem i is 0. ∆Phy i is the sum of the physical processes (including advection, vertical mixing, dry and wet deposition, and convection).Each tracked variable undertakes all the relevant calculations of the physical processes as the normal simulation do, but does not perturb the normal model calculations, so the concentration changes caused by physical processes will be obtained directly.∆Emis i is the gross emission production within region i, determined by: (x, y, z) (x, y, z) i

Emis
:inside region i Emis 0 :outside region i For any grid in domain, if it is in the region i, then the ∆Emis i (x, y, z) equals ∆Emis(x, y, z), otherwise, it equals zero.
Compared with classic sensitivity analysis, this primary PM 2.5 -tagging model provides a more accurate measurement of primary PM 2.5 concentrations in all related regions, while reducing calculation errors (Li et al., 2008).The primary PM 2.5 -tagging by source regions method can also be applied to other pollutants.For example, Gao et al. (2016) used this method to analyze ozone sources over the Yangtze River Delta.

Model Settings and Basic Data
Fig. 1(a) shows the model domains, which included two nested domains on a lambert projection with horizontal resolutions of 36 and 12 km and grids of 90 × 90 and 108 × 108, respectively.The parent domain, D1, supported reasonable meteorological and chemical boundary conditions for the nest domain, D2.Vertically, there were 27 sigma levels between the surface pressure and 50 hPa.The simulation was conducted between 08:00 on October 27, 2015 Local Standard Time (LST) and 08:00 on November 8, 2015 (LST), and the first two days were used as the spin-up time.The carbon-bond mechanism Z (CBM-Z) (Zaveri and Peters, 1999), a gas-phase physical mechanism based on CBM-IV and containing 134 reactions for 55 species, and MOSAIC (4 bins) (Zaveri et al., 2008) aerosol chemistry (divided into four aerosol bins: 0.039-0.156µm, 0.156-0.625µm, 0.625-2.5 µm, and 2.5-10 µm) were applied to the WRF-Chem model.The remaining model parameter options are shown in Table 1.

MODEL EVALUATION
The observed and simulated meteorological data on sealevel pressure (SLP), temperature (T), wind speed (Wspd), and wind direction (Wdir) and the PM 2.5 data from these four cities and eight cities are shown in Figs. 2 and 3, respectively, and all cities statistical results are summarized in Table 3.The Mean Bias (MB), Normalized Mean Bias (NMB), Normalized Mean Error (NME), the Root-Mean-Square Error (RMSE), and correlation coefficient (r) are displayed.We found that the WRF-Chem model accurately simulated trends and peak distribution, as well as the PM 2.5 levels at all stations except SJZ.As shown in Table 3, the correlation coefficient r of Wspd was low (0.46), and RMSE was too high (131.62)at SJZ. Possible causes of these faulty simulations are stronger and higher uncertainty of emission sources, and worse capability in simulation of wind direction and wind speed (in Fig. 2(k)) under complex terrain around SJZ.
Fig. 3 shows that the model successfully reproduced PM 2.5 concentrations during this period.From these PM 2.5 levels, we can see that severe pollution was present on November 3 and 5, with the most serious pollution (PM 2.5 concentration exceeding 300 µg m -3 ) in the region appearing at around 00:00 (LST) on November 5. Based on the PM 2.5 concentration levels and evolution in the BTH region, we divided the pollution process into four stages: stage 1, from October 29 to 30; stage 2, from October 31 to 12:00 LST on November 2; stage 3, from November 3 to 06:00 LST on November 5; and stage 4, from 12:00 LST on November 5 to 12:00 LST on November 7, representing clean, light, heavy, and then clean pollution periods, respectively (We discuss this in more detail in the following part.).

Description of Synoptic Conditions and Simulated PM 2.5
Distributions Meng and Cheng (2002) pointed out that pollutants can accumulate in light horizontal winds on the local scale more easily when Beijing is under uniform pressure type of surface synoptic conditions.The condition of Beijing in front of high-pressure systems has been found to facilitate the spread of pollutants (Miao et al., 2015).Fig. 4 presents surface synoptic conditions, which varied throughout the study period.The wind was blowing northwest when the BTH region was under high-pressure conditions at 02:00 LST on October 30 showed in Fig. 4(a).This high pressure then weakened near the central region of Mongolia, but another high pressures then covered a wide area of southeastern China, leading to a weak pressure gradient and low wind speeds in the BTH region at 02:00 LST on November 2 presented in Fig. 4(b).In Fig. 4(c), the BTH region was in saddle field with low wind speeds at 02:00 on November 4. Finally, strong high-pressure conditions near Mongolia and Inner Mongolia and low pressures in southern China led to a dense pressure gradient and gale at 02:00 LST on November 6 in Fig. 4(d).
Figs. 5(a)-5(d) show the average spatial distribution of surface PM 2.5 and wind fields during the four different stages.From Figs. 5 and Fig. 1(b), we can see that the terrain of the Taihang Mountains (northeast-southwest) and Yanshan Mountain (east-west) served as important demarcation lines for wind conditions and pollution concentrations.The altitude of most of the Taihang Mountains is over 1,200 m, and that of the Yanshan Mountains ranges between 600 and 1,500 m.Thus, this special terrain exerts a significant effect on the BTH region.The terrain weakens the north winds, while contaminants easily accumulate at the foot of the mountain under the control of southerly wind (Miao et al., 2015).
During stage 1 (Fig. 5(a)), the northern BTH region was experiencing a northwest gale, which improved air quality.The pollutants at the foot of the two mountains in the BTH region could not be eliminated completely.However, owing to low wind speeds in southern BTH region, the average PM 2.5 concentration in this region was almost under 60 µg m -3 , within Grade II of the Chinese National Ambient Air Quality Standards.In stage 2 (Fig. 5(b)), when the wind in the south turned southwest wind, light levels of pollutants were beneficial to transport into the BTH region.Moreover, the mountainous terrain weakened the speed of the north wind, which in turn aggravated pollution levels at the foot of the mountain and formed a pollution zone from southwest to northeast (including in Shijiazhuang, Baoding, and Tianjing), similar to findings by Wu et al. (2014) and Ren et al. (2004).It is worth noting that northwestern BJ was still in a clean state owing to greater wind speeds in the North.Shortly afterwards, North China was dominated by southerly winds and wind in south region greater than before (stage 3, Fig. 5(c)), which facilitated the longdistance transport of pollutants from the southern regions.The wind speed over BJ and TJ was low, which allowed pollutants to accumulate after being transported to these areas.Furthermore, BJ was not protected by the northwest wind as the wind shifted south in the northern mountains, and a high level of pollution moved into BJ.The maximum observation concentration of PM 2.5 was 304 µg m -3 , occurring at 20:00 LST on October 4 in BJ.This was also the point at which PM 2.5 reached its widest scope and largest concentration in such a special terrain structure.Finally, when the wind over most of North China shifted northeast (stage 4, Fig. 5(d)), fresh air was driven into the BTH region, and PM 2.5 concentrations decreased quickly to below 40 µg m -3 , within Grade I of the Chinese National Ambient Air Quality Standards.

Primary PM 2.5 Contributions of Different Geographical Source Regions to the BTH Region
Pollution source apportionment is an effective method of studying the source regions of air pollutants.Fig. 6(a) illustrates the primary PM 2.5 contributions of different geographical source regions to the BTH region.Fig. 6(b) presents the average spatial distribution of the primary-tototal PM 2.5 ratio within the study period.We can see that the percentage of primary PM 2.5 almost exceeds 50% in most of domain 1 and nearly surpasses 60% in the BTH region.Because the high ratio of primary PM 2.5 and the precursors of secondary PM 2.5 close related to primary PM 2.5 in source regions, further source apportionment of primary PM 2.5 will help characterize total PM 2.5 .As seen in Fig. 6(a), in stage 1, primary PM 2.5 originated mainly in the local region.The wind on the northern mountains was stronger than the southern.The upstream region with clean air had almost no primary PM 2.5 contribution to the downstream region.In stage 2, when the wind speed become small on the northern side of the two mountains and wind shifted southwest on the southern side, a light pollution zone at the foot of the mountains from SJZ to TJ was created by both wind fields and the terrain.Moreover, northwestern BJ is clean due to it was still experiencing northwest winds, so only southeastern BJ occurred light air pollution.Fig. 7 shows the average primary PM 2.5 contributions to the BTH region by several geographical source regions.We can clearly Table 3. Comparisons of statistical parameters between simulated and observed variables in eight representative cities of BTH region (Unit: PM 2.5 (µg m -3 ); P (pa); T (°C); Wspd (m s -1 ); Wdir (°)).see that BJ was mostly affected by primary PM 2.5 originating locally and from HB, accounting for 68.7% and 23.4%, respectively, of total primary PM 2.5 contributions.In TJ, primary PM 2.5 originated mainly in TJ itself and HB, accounting for 31.2% and 54.2% of the total, respectively.The primary PM 2.5 in HB was mainly from HB itself, accounting for 69.4% of the total.So we can know that the local contribution of PM 2.5 to the BTH region accounted for 90.7% of total primary PM 2.5 before the heavy pollution event (stage 1 and stage 2), and HB, TJ, and BJ made up 53.6%, 21.0%, and 16.1%, respectively.In stage 3, when the wind in the northern mountains shifted southerly, pollutant concentrations rapidly increased in BJ.Southwest winds shifted south, and a gale weakened into soft winds in the middle of the BTH region, which allowing for greater concentrations and a wider range of pollution.In this stage, BJ, TJ, and HB were most affected by the regions also shown in Fig. 7. BJ and HB were the regions that most significantly contributed to primary PM 2.5 in BJ, accounting for 25.0% and 26.7% of the total primary PM 2.5 , respectively.Pollution in TJ was contributed primarily by TJ itself, as well as HB and SD, accounting for 30.5%, 20.9%, and 15.0% of the total, respectively.HB was mainly affected by pollution from HB, HN, and SD, accounting for 38.2%, 16.9%, and 10.2% of the total, respectively.So the BTH region contributions by BTH, Shandong, Henan, Jiangsu and Anhui accounted for 50.7%, 11.6%, 11.8%, 9.2% and 6.7%, respectively.AH, JS, SD and HN all had an obvious effect on the BTH region, and the primary PM 2.5 concentration contribution in S-HB (Shijiazhuang, Xingtai, Handan, Hengshui, Cangzhou, Baoding, and Langfang) was greater than in N-HB (Zhangjiakou, Tangshan, Chengde and Qinhuangdao).Finally, pollution in the BTH region cleared out quickly owing to the shifting of the northeast gale in stage 4.

Sites
Through all the study period, by combining the pollution processes with the synoptic conditions map, we determined that the pollution phenomenon began with the appearance of high-pressure conditions in southern China and decreasing high-pressure in northern China in stage 2, while an isobaric synoptic condition in the BTH region led to the most severe pollution and long-distance transport of pollutants from JS, AH, and HN to the BTH region in stage 3.So we will discuss stage 3 in detail to reveal the reason of the severe pollution.

Long-Distance Transport during Heavy Pollution Periods in Stage 3 PBL Characteristics
Lyu et al. (2012) indicated that Planetary Boundary Layer Height (PBLH) determines the efficiency of pollutant diffusion.Higher PBLH is beneficial to pollutant diffusion, while lower PBLH allows pollution to accumulate in the boundary layer (Guo et al., 2011;Yang et al., 2011;Ye et al., 2011).As shown in Fig. 8, only 200 m of the PBLH and horizontal distribution of wind divergence were drew in stage 3, we can see that the PBLH of BJ and SJZ was below 200 m.Compared Fig. 5(c) with Fig. 8, which reveal high PM 2.5 concentrations correspond to low PBLH.And we also can clearly see a large wind convergence region (the parallelogram area in Fig. 8) caused by high wind speed in the south and low speed in the region of the southwestnortheast direction (from BJ to SJZ).Therefore PM 2.5 pollutants easily accumulated in this parallelogram area, which is similar to the study by Yang et al. (2016).
To help characterize the vertical structure, we used the profiles shown in Figs.9(a) and 9(b) (from the black straight lines in Fig. 8).Fig. 9(a) (from the black southeast-northwest straight line in Fig. 8, including BJ and TJ) reveals that pollutants were continuously transported to the BTH region by the southerly wind, but the tall mountains blocked their  8, including SJZ and the midpoint between TJ and BJ) shows that SJZ, BJ, and TJ lie in a large valley.And Pollutants were transported to the valley by the southerly wind and did not spread out easily because of the terrain.There is also a temperature inversion and a low PBLH in this valley.The height of the PM 2.5 pollution in the BTH region is higher than the other regions.
In general, a PBLH near 200 m with a temperature inversion was not conducive to the vertical diffusion of pollutants, which also further aggravated surface pollution levels in the BTH region.However, the BTH region is surrounded on three sides by tall mountains, which lead to vertical airflow.Therefore, the pollutants transported by the southerly wind were unable to spread horizontally, forced them to extend nearly 3 km upraised by terrain.

Long-Distance Transport Implicated by Primary PM 2.5 Tagging
Using the tagging method, Fig. 10 shows the distributions of average primary PM 2.5 contributions by four major source regions (JS, AH, HN, and SD) in stages 3a (on November 3) and 3b (from November 4 to 06:00 LST on November 5).In stage 3a (Figs.10(a)-10(d)), JS, AH, HN and SD were behind high-pressure conditions, so this four source regions were controlled by the southerly wind.There was almost no contribution from JS on BTH, although the primary PM 2.5 concentration in JS reached about 100 µg m -3 .AH only exerted slight influence on S-HB.HN and SD are near BTH, so pollutants from HN and SD were therefore transported to BTH more easily than from JS and AH.Therefore, HN affected the whole BTH region, and SD mainly affected TJ and S-HB via the southerly wind.
As high-pressure conditions moved east (stage 3b, Figs.10(e)-10(h)), the site of gale also moved east so that JS experienced a southeast gale.JS significantly contributed pollutants to BTH, especially in southeastern HB and TJ.The primary PM 2.5 from AH and SD almost covered the entire BTH region.The wind speed in HN decreased after high-pressure conditions moved east, but HN continued to contribute pollutants to the BTH region.SD has a significantly contribution to the most BTH region.
Fig. 11 presents the average primary PM 2.5 contributions by different source regions in stages 3a and 3b.In stage 3a, primary PM 2.5 concentrations in the BTH region mainly originated locally and from HN, SX&SN.SD only contributed greatly to the TJ region (average primary PM 2.5 concentration 9.3 µg m -3 , accounting for 18.3%), while HN contributed significantly to the BTH region throughout stage 3a.During this period, JS exerted almost no effect on BTH.In stage 3b, the contributions in most significant source regions were greater than stage 3a, although the direct line distance from JS to BTH is about 850 km, but JS had a significant impact on BTH, especially in TJ, where the mean value reached 14.5 µg m -3 (accouting for 18.3%).The contributions by SD, HN, AH, and JS to S-HB were nearly twice those to N-HB during stage 3b, thus, the major polluted area in HB was S-HB.
In conclusion, under the influence of a high-pressure system, the JS, AH, HN and SD region contributed significant levels of primary PM 2.5 to the BTH region via the southerly wind.So a source in the Yangtze River Delta city group region played a significant role in the BTH region.Although JS only affected the BTH region about one day, its primary PM 2.5 contribution could not be ignored.The higher primary PM 2.5 concentration in HB illustrated that HB, especially S-HB, played a key role in high primary PM 2.5 events over the BTH region.

CONCLUSION
In this study, the WRF-Chem model coupled with an online primary aerosol tagging method was used to simulate   an episode of high concentration of PM 2.5 that occurred in the BTH region in late October and early November 2015.The primary PM 2.5 contributions by different geographical source regions to the BTH region were quantified.Using the model, we evaluated both meteorological variables and PM 2.5 concentrations with confidence, and our conclusions are summarized below.
The simulation showed that pollution in the BTH region was mainly affected by synoptic conditions, the terrain, and boundary layer characteristics.The distribution and concentrations of pollutants differed under various surface synoptic conditions, with pollution primarily originating locally from October 29 to 30, locally and in surrounding areas from October 31 to 12:00 LST on November 2, and in a wide geographic area (and brought to the BTH region via long-distance transport ) from November 3 to 06:00 LST on November 5.
The clear contribution and spatial distribution of average primary PM 2.5 from different geographical source regions showed the trajectory of primary PM 2.5 transported to the BTH region.To a certain extent, the primary PM 2.5 with a high percentage within total PM 2.5 can represent the characteristics of total PM 2.5 .The primary PM 2.5 source appointment results clearly presented that before the heavy pollution event, the BTH region was mainly originated locally accounted for 90.7% of the total primary PM 2.5 .When the BTH region was under isobaric synoptic condition and the pollution was most serious throughout the four stages, contributed by locally and other polluted areas, especially SD, HN, AH, and JS, with average primary PM 2.5 contributions of 50.7%, 11.6%, 11.8%, 6.7% and 9.2%, respectively.Primary PM 2.5 from local (the BTH region) and non-local regions contributed nearly equally during the heavy pollution.These accurate contribution results showed that the BTH region was affected by not only the locally and the surrounding cities but also the areas from the regional scale or even larger.
Vertically, during periods of heavy pollution, a PBLH near 200 m with a temperature inversion and the terrain of the BTH region surrounded on three sides by tall mountains were not conducive to the vertical and horizontal diffusion of pollutants.However, the vertical airflow caused by the terrain forced the pollutants to upraise and extend nearly 3 km via the southerly wind.

Fig. 5 .
Fig. 5. a-d) Average spatial distribution of surface PM 2.5 and wind fields during four stages.

Fig. 6 .
Fig. 6. a) Primary PM 2.5 contributions of different geographical source regions to the BTH region (including N-HB and S-HB), b) average spatial distribution of the primary-to-total PM 2.5 ratio.

Fig. 7 .
Fig. 7. Average primary PM 2.5 contributions by different geographical source regions to the BTH region during stages 2 and 3.

Fig. 8 .
Fig. 8. Distribution of average PBLH and wind divergence in stage 3 (blue contours: boundary layer height of 200 m, black straight line: the positions of two vertical profiles used in Fig. 9, yellow shaded: horizontal distribution of wind divergence).

Fig. 9 .
Fig. 9. a) Vertical profile of temperature and PM 2.5 concentration with southeast-northwest boundary layer height, b) vertical profile of temperature and PM 2.5 concentration with southwest-northeast boundary layer height (black line: isotherm, blue line: boundary layer height, black areas: terrain , and the vertical wind was enlarged 200 times).

Fig. 10 .
Fig. 10.Spatial distribution of average primary PM 2.5 contributions by four major source regions during stage 3a: a) JS, b) AH, c) and HN, d) SD and stage 3b: e) JS, f) AH, g) HN, and h) SD.

Fig. 11 .
Fig. 11.Contributions of average primary PM 2.5 by different source regions to BJ, TJ, and HB (S-HB and N-HB) in stages 3a and 3b.

Table 1 .
Physical parameters used in the modeling study.