Simulations of Sulfate-Nitrate-Ammonium (SNA) aerosols during the extreme haze events over Northern China in October 2014

: Extreme haze events have occurred frequently over China in recent years. Although many studies have investigated the formation mechanisms associated with PM 2.5 for heavily polluted regions in China based on observational data, adequately predicting peak PM 2.5 concentrations is still challenging for regional air quality models. In this study, we evaluate the performance of one configuration of the Weather Research and Forecasting model coupled with chemistry (WRF-Chem) and use the model to investigate the sensitivity of heterogeneous reactions on simulated peak sulfate, nitrate, and ammonium concentrations in the vicinity of Beijing during four extreme haze episodes in October 2014 over the North China Plain. The highest observed PM 2.5 concentration of 469 μg m -3 occurred in Beijing. Comparisons with observations show that the model reproduced the temporal variability in PM 2.5 with the highest PM 2.5 values on polluted days (defined as days in which observed PM 2.5 is greater than 75 μg m -3 ), but predictions of sulfate, nitrate, and ammonium were too low on days with the highest observed concentrations. Observational data indicate that the sulfur/nitric oxidation rates are strongly correlated with relative humidity during periods of peak PM 2.5 ; however, the model failed to reproduce the highest PM 2.5 concentrations due to missing heterogeneous reactions. As the parameterizations of those reactions is not well established yet, estimates of SO 2 -to-H 2 SO 4 and NO 2 /NO 3 -to-HNO 3 reaction rates that depend on relative humidity were applied which improved the simulation of sulfate, nitrate, and ammonium enhancement on polluted days in terms of both concentrations and partitioning among those species. Sensitivity simulations showed that the extremely high heterogeneous reaction rates and also higher emission rates than those reported in the emission inventory were likely important factors contributing to those peak PM 2.5


Introduction
Anthropogenic PM 2.5 (fine particulate matter with aerodynamic diameters less than 2.5 µm) is known to play a significant role in atmospheric visibility, human health, and climate.Regional haze with extremely high PM 2.5 concentrations (exceeding the WHO standard 10-fold) has become the primary air quality concern in China, especially over the North China Plain (NCP).Severe haze pollution episodes occur frequently over the NCP, especially in Beijing (e.g., Sun et al., 2015), during almost all seasons as described by Wen et al. (2015) for the summer of 2013, Liu et al. (2013), and Yang et al. (2015) for the autumns of 2011 and 2014, and Y. S. Wang et al. (2014) and Han et al. (2015) for the winters of 2013 and 2014.Those haze events provide an opportunity to test the current understanding of sources and the formation mechanisms of PM 2.5 over the NCP as represented in current air quality models.Most air quality modeling studies have focused on winter or summer periods (e.g., L. T. Wang et al., 2014Wang et al., , 2016; Y. X. Wang et al., 2014;Chen et al., 2015;Zhang et al., 2015a, b;Zheng et al., 2015), because emissions resulting from heat production and adverse meteorological conditions during the winter and greater photochemical production rate during the summer play crucial roles in the formation of haze in the two seasons, respectively.Although observations (W.Wang et al., 2014) suggest that the biomass burning associated with autumn harvest in the NCP could be an important source of PM 2.5 , the formation mech-anisms and the performance of models in autumn have not been thoroughly investigated.
Sulfate, nitrate, and ammonium (denoted as SNA) are the predominant inorganic species in PM 2.5 .Observations during the winter of 2013 (e.g., Y. S. Wang et al., 2014) and autumn of 2014 (Yang et al., 2015) show that SNA increases rapidly during the highest haze episodes over the NCP and makes up approximately half of the total PM 2.5 mass.Based on the studies for the winter of 2013 and the autumn of 2014, this rapid SNA increasing is associated with relative humidity (RH).There are strong correlations between RH and sulfur/nitric oxidation ratios, thus the heterogeneous reactions of precursors (sulfur dioxide SO 2 , nitrogen oxides NO x = NO + NO 2 ) are key reactions that lead to the formation of sulfate and nitrate on the surface of particulates (e.g., Li andShao, 2009, 2010;Li et al., 2011; X. F. Wang et al., 2012;Zhao et al., 2013; Y. S. Wang et al., 2014;Yang et al., 2015).However, those reactions are not included in current chemical mechanisms (traditional gas-phase or aqueousphase chemistry) in most air quality models.After identifying this deficiency, Y. X. Wang et al. (2014) and Zhang et al. (2015b) introduced and parameterized the heterogeneous uptake of SO 2 on deliquesced aerosols in the GEOS-Chem model and Zheng et al. (2015) comprehensively evaluated the effects of heterogeneous chemistry in the CMAQ model.Their simulations for the conditions during the 2013 winter showed great improvements when heterogeneous chemistry was included.
Precursor emissions are also important factors in determining SNA concentrations and composition.It is interesting to explore the response of SNA to precursor emission changes due to the uncertainties when applying the emission inventories in the model.Rapid changes related to historical and recent economic developments, or specific political decisions toward emission reductions are often not accounted for by emission inventories.This is particularly the case for the developing countries such as China (Richter et al., 2005).We use the latest emission inventory available for 2010, thus the changes of three precursor emissions (SO 2 , NO x , and NH 3 ) between 2010 and 2014 may not be reflected in our simulations.There is ample evidence to show that SO 2 emissions have decreased rapidly in recent years (especially after 2010-2011) from both bottom-up emission inventories (Lu et al., 2011;Xia et al., 2016) and satellite observations (Wang et al., 2015, Krotkov et al., 2016), driven by wide implementation of flue gas desulfurization in power plants.Around 30 % reductions in SO 2 vertical column densities (VCDs) were reported in eastern China (Krotkov et al., 2016) from 2010 to 2014 and the reduction ratio reaches 50-60 % from 2007-2008 to 2014.For NO x emissions that are primarily from vehicles and industries, inventory estimates based on both bottom-up and satellite remote sensing methods show increases nationwide before 2010 (Lamsal et al., 2011;Zhang et al., 2012;Shi et al., 2014;Krotkov et al., 2016;Xia et al., 2016) and also increases over western China before 2013 (Cui et al., 2016).The Chinese government set a 10 % NO x emission reduction target in the 12th Five-Year Plan (2011)(2012)(2013)(2014)(2015).Satellite observations show that the NO x emissions reached their peak in 2011 over eastern China and then started to decrease during 2012.Compared to 2010, the NO x VCDs in 2014 decreased by 20 % over eastern China (Krotkov et al., 2016).Among the three precursors, NH 3 emissions that are dominated by agricultural sources have the largest uncertainty in terms of the total amount and also seasonal variations.Estimates of NH 3 emissions vary greatly among published papers (Streets et al., 2003;Kim et al., 2006;Dong et al., 2010;Zhang et al., 2010;Huang et al., 2012;Xu et al., 2015).Using the GEOS-Chem model, Wang et al. (2013) and Zhang et al. (2015b) concluded that NH 3 emission plays a critical role in the SNA simulations in the NCP.In addition to the uncertainties in the emission inventory, determining the dynamic variation in emissions due to unpredictable conditions (such as crop biomass burning) during the haze periods is also challenging.
In this study, we parameterized the SNA relevant heterogeneous reactions in WRF-Chem and conducted simulations in the NCP for a haze period in the autumn of 2014.To our best knowledge, this is the first study that SO 2 -NO 2 -NO 3 heterogeneous reactions were taken into account in WRF-Chem.We first evaluate the model results using available surface observations.Then, the relative importance of precursor emissions and the missing heterogeneous chemistry on the simulated results are quantified.We also conducted sensitivity simulations of heterogeneous reaction rates and precursor emissions.The model configuration, observations, and methodology of how heterogeneous reactions are treated are described in Sect. 2. In Sects. 3 and 4, the model evaluation and sensitivity analysis for heterogeneous reaction rates and precursor emissions are presented, respectively.The concluding remarks are given in Sect. 5.
2 Model description, observations, and methodology

WRF-Chem model
The Weather Research and Forecasting (WRF) model coupled with online chemistry (WRF-Chem) is based upon the nonhydrostatic WRF community model (http://www.wrf-model.org/index.php).Details of the WRF-Chem model are described by Grell et al. (2005) and Fast et al. (2006), and there have been many subsequent papers describing recent updates.We used version 3.6.1 in this study and a summary of physical parameterization options is shown in Table 1.The model domain with a 40 km horizontal grid spacing covers most of China and the surrounding region (left panel in Fig. 1) and our interest is the NCP (right panel in Fig. 1).There are 57 vertical levels extending from the surface to 10 hPa.The WRF single-moment 6-class microphysics scheme (Grell and Devenyi, 2002) and the Grell 3-  ) that illustrate the regional pollution on 10 October.
D cumulus parameterization were used to treat clouds and precipitation.The Noah parameterization is used to represent land surface processes and the YSU parameterization is used to represent boundary layer turbulent mixing (Hong et al., 2006).Initial conditions for meteorological variables are obtained from the National Center for Environmental Prediction's (NCEP) Global Forecast System (GFS) analyses that are updated every 6 h.The lateral boundary conditions (LBCs) for the meteorological fields are also provided by the GFS analyses.LBCs for chemistry and aerosol fields are based on prescribed idealized profiles.The simulation started on 1 October 2014 and the first 3 days are treated as a spin-up period and are not used in our analyses.The Carbon Bond Mechanism version Z (CBMZ) and Model for Simulating Aerosol Interactions and Chemistry (MOSAIC) are used as the gas-phase and aerosol chemical mechanisms, respectively, in this study.Aerosol species in MOSAIC are defined as black carbon, organic compounds, sulfate, nitrate, ammonium, sodium, chloride and other inorganic compounds.MOSAIC uses a sectional approach to represent the aerosol size distribution with 4 and 8 size bins available in the public version of the code.Research versions of MOSAIC have used up to 20 size bins in WRF-Chem (Lupascu et al., 2015).In this study, we use 4 size bins with aerosols diameters ranging 0.039-0.1,0.1-1, 1-2.5, and 2.5-10 µm.The Fast-J photolysis scheme is used for photolytic rate calculations.In the standard simulation with the CBMZ-MOSAIC mechanism, the SNA aerosol formation is through oxidation and neutralization/condensation of precursor gases.The sulfate formation starts from the SO 2 to sulfuric acid (H 2 SO 4 ) oxidation, including two pathways in the model, the gas-phase oxidation of SO 2 by hydroxyl radicals (OH), hydrogen peroxide (H 2 O 2 ) and ozone (O 3 ), and aqueous-phase oxidation of SO 2 by H 2 O 2 and O 3 in clouds.The nitrate formation (NO x to nitric acid HNO 3 oxidation) also includes two pathways: the NO 2 oxidation by OH during the daytime and the newly added hydrolysis of dinitrogen pentoxide (N 2 O 5 ) at night (Archer-Nicholls et al., 2014).H 2 SO 4 and HNO 3 are neutralized/condensed mainly by/with NH 3 to form (NH 4 ) 2 SO 4 and NH 4 NO 3 , respectively.As H 2 SO 4 is nonvolatile, the equilibrium surface concentration is assumed to be zero in the model.(NH 4 ) 2 SO 4 is the preferential species in the completion when H 2 SO 4 and HNO 3 are both present and NH 4 NO 3 is formed only if excess NH 3 is available beyond the sulfate requirement.Thus, the amount of NH 3 is a key factor in determining the SNA aerosol formation in the NH 3 -limited environment.

Emissions
The Multi-resolution Emission Inventory for China (MEIC) (Zhang et al., 2009;Lei et al., 2011;He, 2012;Li et al., 2014) for October 2010 is used as the base emission scenario.The original grid spacing of this emission inventory is 0.25 • × 0.25 • and it has been processed to match the model grid spacing (40 km).The gas precursor emissions over China for October are estimated to be 2.03 Tg for SO 2 , 2.07 Tg for NO x (NO 2 + NO), and 0.63 Tg for NH 3 in the base emission scenario and the values over the NCP are 0.60 Tg (9.32 Gmol) SO 2 , 0.63 Tg (13.8 Gmol) NO x , and 0.13 Tg (7.8 Gmol) NH 3 , respectively.The spatial distributions of these three species are shown in Fig. 2. On the molecular basis, NH 3 emissions were less than the sum of 2 × SO 2 and NO x emissions indicating NH 3 -limited conditions over the NCP.
The MEIC-2010 emission inventory has already been applied in other studies (e.g., L. T. Wang et al., 2014Wang et al., , 2016;;Zheng et al., 2015) for simulations over the NCP during the winter of 2013.They found that this inventory provides reasonable estimates of total emissions from cities over the NCP but is subject to uncertainties in the spatial allocations of these emissions over small spatial scales.For our simulation, uncertainties may also arise from two other aspects: the difference between the emission base year and our simulation year, and the monthly allocations (especially for NH 3 ).To address these uncertainties, sensitivity simulations are performed that alter the emission rates as described in Sect.2.5.

Observations
The meteorological data used for all the sites in the NCP are obtained from the National Climate Data Center (NCDC) integrated surface database (http://www.ncdc.noaa.gov/data-access/), shown as red circles in Fig. 1.The following parameters were evaluated: temperature and relative humidity at 2 m (T2 and RH2), wind speed and direction at 10 m (WS10 and WD10), and 24 h accumulated precipitation.Most of the data are of 3 h frequency (instantaneous every 3 h) except for precipitation.Since there are no meteorological data at sites where hourly PM 2.5 species were observed, the hourly data of T2 and RH2 from the Global Data Assimilation System (GDAS) and boundary layer height from the NOAA Air Resources Laboratory (http://www.arl.noaa.gov/index.php)at those PM 2.5 locations were used for evaluation purposes.The meteorological performance is quantified in terms of both site-by-site and also domain-wide overall statistics.The statistical measures calculated include the mean bias (MB), the root mean square error (RMSE), and the correlation (R).
The observed chemical concentrations used in the study are from three data sets: (1) the daily mean concentration of gas phase and PM 2.5 from the Air Pollution Index (API) database in 10 cities in the NCP (shown as black dots in Fig. 1); (2) the average hourly concentrations of gasphase pollutants and PM 2.5 at 34 monitoring sites in Beijing from the China National Environmental Monitoring Center (CNEMC); and (3) the hourly PM 2.5 measured by TEOM (tapered element oscillating microbalance, RP1405F) and 15 min species concentrations (BC, sulfate, nitrate, ammonium) in PM 1 measured in situ by an aerosol chemical speciation monitor (ACSM) at the Beijing Normal University (BNU, blue dot in Fig. 1) from Yang et al. (2015).Details of the BNU instruments can be found in Sun et al. (2013).As the 34 monitoring sites in Beijing fall into eight model grids, the observations within the same grid are averaged and then the averages of the eight grids are compared with the model predictions.There were also five monitoring sites in the grid cell where the BNU site is located.The average SO 2 and NO 2 among the fives sites are also used for oxidation rate calculations.

Heterogeneous reactions
As discussed in Sect.2.1, there are several pathways for SNA formations in the standard version of WRF-Chem.Although aqueous-phase processes within clouds affect sulfate, we found that the precipitation and cloud amounts in October 2014 are very low from both observation and model.Observed averaged 24 h precipitation is less than 0.25 and where the subscript i represents the ith reactant for heterogeneous reactions, d p is the effective diameter of the particles (m), D i is the gas-phase molecular diffusion coefficient for reactant i(m 2 s −1 ), v i is the mean molecular speed of reactant i in the gas phase, γ i is the uptake coefficient for reactant i (dimensionless), and S p is the aerosol surface area per unit volume of air (m 2 m −3 ).
For the most important parameter -uptake coefficients γ i -we used a similar method as in K. Wang et al. (2012) and Zheng et al. (2015).The lower and upper limits are used to present a range of γ values in the laboratory measurements which were applied when RH is lower than 50 % and higher than 90 %, respectively.The values in the 50-90 % RH range are linearly interpolated based on the two limits.The lower and upper limits of NO 2 -and NO 3 -related reactions are based on K. Wang et al. (2012) and those of SO 2related reaction are based on Zheng et al. (2015) (Table 2).We also performed sensitivity tests for those parameters as described in Sect.2.5.

Scenarios
We conducted several experiments aiming to test the model response to different heterogeneous uptake coefficients and different precursor emissions (see Table 3).The baseline emissions used in BASE and HET_BASE simulations are described in Sect.2.2.The heterogeneous reactions parameterization depending on RH is described in Sect.2.4.To determine the appropriate reaction rates, we first separated the impacts of SO 2 and NO 2 -NO 3 relevant heterogeneous reactions on SNA simulations by turning off the NO 2 -NO 3 relevant reactions in different simulations (Fig. S1 in the Supplement).We found that with only the SO 2 heterogeneous reactions simulated nitrate decreases due to competition of SO 2− 4 and NO −  2. When conducting the emission sensitivity scenarios, we took into account the trends of the national emission in recent years and their uncertainties.According to the annual National Environmental Statistical Report, the total amount of SO 2 emissions in China reached the peak in 2006 and subsequently decreased (by ∼ 25 % from 2006 to 2014 nationally).As the MEIC-2010 emissions inventory relied on the annual statistical books in which the data are often 2-3 years older than the actual year, we assumed that the SO 2 emission levels in MEIC-2010 were closer to the previous 2-3 years (2007)(2008).In the NCP where more strict measurements were implemented, larger-than-nationalaverage reduction ratios were expected.For the above reasons, we applied 25 % reduction of SO 2 in the HET_EMIS scenario.Since the NCP is in NH 3 -limited conditions, we increased the NH 3 emission by 30 % to test the sensitivity in HET_EMIS.The emission changes were applied to the whole modeling domain.It should be noted that the NCP is still under NH 3 -sensitive condition after applying the SO 2 and NH 3 emission changes (please refer to the NCP emission totals in Sect.2.2) We also tested three additional scenarios for the 24-25 October period, aiming to improve the simulated peak values of observed SNA aerosols during the highest polluted event.
Although the emission changes in HET_EMIS may reflect some emission trends in recent years and also the uncertainty of the emission inventory, it does not account for dramatic emission increases that are possible during actual polluted events.In addition, the dependent reaction rates will be too low when the simulated RH is also too low.For these two reasons, we increased the emissions and adjusted those uptake coefficients as fixed values in HET_MAX1, HET_MAX1_EMIS1, and HET_MAX2_EMIS2.It should be noted that we did not have any evidence to determine the SO 2 uptake coefficient in those sensitivity simulations due to the limited observations (only one site for a few days).The values here are chosen to best match the observations.More rigorous optimization tests should be conducted in future studies.3 Model evaluation

Meteorology
Table 4 quantifies the performance of the meteorological predictions based on the comparisons with the NCDC data set.
In addition to the averaged statistics in two regions (NCP and the whole study domain), statistics for the three sites (PEK, Baoding, and Shijiazhuang) are also listed.PEK is the only NCDC observational site in Beijing and in our simulation it is one grid cell north of the BNU site where PM 2.5 species were observed.Baoding and Shijiazhuang are two sites generally upwind of Beijing and considered as the major PM 2.5 sources in the region.The temperature and relative humidity at 2 m are overall underestimated in both the NCP (with biases of −0.09 degree for T2 and −7.09 % for RH2) and the whole domain (−0.78 degree and −1.91 %, respectively), while wind speed at 10 m is overestimated with biases of 0.99 m s −1 (27.4 %) over the NCP and 0.78 m s −1 (22.6 %) over China.It should be noted that the intervals of WS10 in NCDC observations are 1 m s −1 and the minimum valid value is 1 m s −1 so that no data are provided when wind speed is less than 1 m s −1 .This may cause some biases in the comparison especially when the simulated wind speeds are smaller than 1 m s −1 .The averaged correlations in the two regions ranged 0.91-0.93 for the 2 m temperature (T2), 0.73-0.79for the 2 m relative humidity (RH2), and 0.52-0.53for the 10 m wind speed (WS10).Compared with another study focusing on the China-NCP that used nested domains in WRF-Chem (Wang et al., 2016), the performance of meteorological simulation in both China and NCP regions are comparable in terms of similar ranges of bias and RMSE, except for the simulated water vapor/RH and wind speed over the NCP.In Wang et al. (2016), the simulated wind speeds in the outer domain (36 km over China) during January 2013 were overestimated by more than 30 % but the overestimation was reduced to less than 10 % within the inner domain (12 km over the NCP).This suggests that nested  domains may help reduce the overprediction in wind speeds.
In Wang et al. (2016), water vapor was usually overestimated by 8-14 % during January 2013, while the RH was underestimated by 2-7 % during October 2014 in our study.Overall, the statistics suggest that the meteorological performance of the model is reasonable for the NCP and the whole domain, while the overestimated wind speed may lead to some negative biases of the chemical species and the underestimated RH may also cause negative biases of RH-dependent reaction coefficients in the simulations.For site-by-site comparisons, the monthly-mean 2 m relative humidity biases at PEK, Baoding, and Shijiazhuang are −11.80,−20.29, and −16.10 %, respectively, which are larger than the NCP average.The biases in the 2 m temperature are also larger at the three sites (from −0.24 to 1.10 degrees).The 10 m wind speed biases are relatively lower at PEK (0.36 m s −1 ) and Baoding (0.37 m s −1 ) than the NCP average, but higher than the average at Shijiazhuang (1.50 m s −1 ).
To clearly show the fluctuations of these meteorological parameters, the time series (in local time) of T2, RH2, WS10, and WD10 at the PEK site are shown in the left panels of Fig. 3. To illustrate the relationship of meteorology and PM 2.5 , the time series of T2, RH2, PBLH, and PM 2.5 at the BNU site are also shown in Fig. 3 (right panels).From the comparisons at both sites, we can see that the T2, WS10 and WD10 simulations in Beijing are reasonable, but the model failed to reproduce high RH on some days.Large low biases of 20-30 % were produced for 13-20, 25, and 29-31 October.The lowest bias occurred on 25 October when the observed RH was between 70 and 100 % and the simulated values were around 50-70 %.The model successfully captured the observed temporal variations of the boundary layer height, including the relatively low boundary layer height during most of the days and also the rapid increases (along with strong winds) on 5, 12, 15, and 26 October.The time series of both the observations and simulations show strong anticorrelations for boundary layer height and PM 2.5 , but positive correlations for boundary layer height and 10 m wind speed, which indicates that boundary layer height and wind patterns are the important meteorological factors that contributed to PM 2.5 accumulation in Beijing.

PM 2.5 and gas-phase pollutants
To evaluate the performance of the transport and chemistry in the original model configuration, the observed 24 h averaged PM 2.5 concentrations at 10 national monitoring sites in the NCP are compared with the BASE simulation in Fig. 4.While Fig. 1 shows that PM 2.5 pollution has regional variations over the NCP, Fig. 4 illustrates four pollutant events with peak values on 9, 19, 24, and 31 October that occurred on almost the same days at those 10 cities.The correlations are most obvious for Beijing, Baoding, and Shijiazhuang.The BASE simulation reproduced the overall PM 2.5 levels Figure 4. Time series of observed 24 h averaged PM 2.5 concentrations (µg m −3 ) at 10 national monitoring system sites in the NCP (locations as black dots in Fig. 1), along with the results from the BASE simulation.and also the four peak events at these sites, but the peak values were too low.The daily mean PM 2.5 values at the 3 aforementioned cities are highest among the 10 cities and bias in the model is also the largest.
Figure 5 shows the time series of observed and simulated hourly PM 2.5 , PM 10 , and four trace gases in the vicinity of Beijing.The observational data are averages among 34 sites.Observed NO 2 and carbon monoxide (CO) show four pol-lution events that are consistent with the high PM 2.5 events.CO is an important tracer since the local ambient concentrations depend mostly on emission rates, transport, and turbulent mixing.Simulated CO level is reasonable for relatively clean days when observed CO is less than 0.7 mg m −3 , but significantly underestimates peak values by 50-70 % for the four polluted events.Since the boundary layer heights and winds are reasonably simulated (Fig. 3), it is possible that the CO emission rates were higher during these events, which were not reflected in the baseline emission scenario.If we assume that CO emissions are underestimated in the baseline emission scenario for the peak days, we should also expect NO x emissions to be underestimated for these days since they both are emitted primarily from vehicles.However, the simulated NO 2 is reasonable, except that it is too low only on 17-18 October.Simulated O 3 is too low during the four peak events based on the averages of the 34 sites (in urban/suburban regions).O 3 formation is more complex as the dependence of O 3 production on NO x and VOCs is significantly different between the so-called NO x -limited regime (in the Beijing urban region) and the VOC-limited regimes (in the suburban region).For SO 2 , the major sources are from combustion (during the heating season) and industry.Severe overestimation (around 3 times the amount) of SO 2 concentrations and underestimation of other pollutants may indicate that the SO 2 emissions from industry are overestimated for the October simulation in Beijing, consistent with new regulations applied in recent years.As shown in Krotkov et al. (2016), 30 % reductions in SO 2 vertical column densities were observed from OMI over eastern China from 2010 to 2014, and the reduction ratio reaches 50-60 % for the pe-riod from 2008 to 2014.The reductions in Beijing are likely larger than the eastern China average since more strict pollution control policy was implemented.Additionally, the overestimation may also indicate that the SO 2 oxidation is too weak/slow in the model.These comparisons not only serve as an evaluation of the simulated chemistry but also provide some insights on how uncertainties in the emission estimates affect the simulated values.

SNA response to heterogeneous reaction rates and precursor emissions
The observations at the BNU site provided not only hourly PM 2.5 concentrations but also concentrations of black carbon, ammonium, nitrate, and sulfate.The latter three are valuable to investigate (1) the impacts of newly added heterogeneous reactions to SNA and PM 2.5 simulation in the model and (2) the relative role of heterogeneous reactions and precursor emissions.
Figure 6 shows the hourly PM 2.5 concentrations as well as the PM 1 BC, sulfate, nitrate, and ammonium concentrations between 15 and 31 October at the BNU site.This period covers three of the highest polluted events when PM 2.5 Figure 6.Time series of observed (OBS) and simulated (BASE, HET_BASE, HET_EMIS) hourly PM 2.5 and PM 1 species (µg m −3 ) including black carbon (BC), sulfate (SO 4 ), nitrate (NO 3 ), and ammonium (NH 4 ).Descriptions of model scenarios are given in Table 3. concentrations started to increase from very clean conditions to larger than 200-400 µg m −3 .The BASE simulation reproduced the 27-31 October event in terms of PM 2.5 and the four species, but it failed to reproduce the 16-21 October and 21-25 October events.The observed BC in PM 1 clearly shows the increase of primary BC emissions during the two events that the model failed to reproduce, indicating a likely underestimation of the primary PM 2.5 emission during the two events.Comparisons of SNA aerosols for the 16-21 October event show that the BASE simulation captured the SNA aerosol species with small low biases on 20 October.During this event, the underestimation of PM 2.5 might be due to organic compounds and other inorganic species.While observed SNA aerosols increased dramatically starting on 24 October and almost doubled (nitrate and ammonium) or tripled (sulfate) within 48 h during the 21-25 Octo-ber event, the BASE simulations severely underestimated the SNA aerosols, especially for sulfate.
Compared with the BASE simulations, the newly added SO 2 heterogeneous reactions in HET_BASE simulations did increase the PM 1 sulfate concentrations with the largest increase occurring on 23 October when sulfate almost doubled from the BASE to HET_BASE simulations.Since the reaction rates are RH dependent, simulated sulfate was overestimated when the RH is too high on 23 October.While the simulated RH is too low (Fig. 3) for the observed peak SNA event on 24-25 October, the sulfate in HET_BASE did not improve too much.The sensitivity simulation results for this 2-day period will be discussed later.As the competition of SO 2− 4 and NO − 3 to form sulfate and nitrate, respectively, in NH 3 -limited conditions, new NO 2 -NO 3 heterogeneous reactions should also be added in the model along with the SO 2 heterogeneous reactions to avoid the nitrate decrease.Compared with the BASE simulation, nitrate in HET_BASE remained nearly the same and ammonium increased slightly.Although HET_BASE changed the SNA ratios, especially the sulfate to nitrate ratio, the total PM 2.5 mass did not increase significantly.This is because Beijing is in NH 3 -limited condition and the SNA mass is highly dependent on the NH 3 emissions.In HET_EMIS, the 30 % NH 3 emission increase leads to noticeable nitrate and ammonium increases when compared to HET_BASE.Moreover, the 25 % SO 2 emission decrease in HET_EMIS leads to slight sulfate decrease.The total PM 2.5 mass increases in HET_EMIS when compared with BASE and HET_BASE, but HET_EMIS still has a large low bias between 24 and 25 October.The mass fraction, depicted as a pie chart in Fig. 7, more clearly illustrates the differences in the SNA ratio changes among the simulations.In this figure, composition concentrations on polluted days (when observed PM 2.5 is larger than 75 µg m −3 ) are averaged between 15 and 31 October.In the BASE simulation, the sulfate and nitrate fractions are significantly lower on polluted days.With the newly added heterogeneous reactions, the SNA fractions from HET_BASE are very close to observations.The sulfur oxidation ratio (SOR = nSO 2− 4 /(nSO 2− 4 + nSO 2 ) (n refers to the molar concentration) and the nitric oxidation ratio (NOR = nNO − 3 /(nNO − 3 + nNO 2 ) are important factors showing the gaseous species oxidation rates and the secondary transformation (Sun et al., 2006(Sun et al., , 2013)).The high fractions of sulfate and nitrate in heavily polluted episodes could be related to the high oxidation rates of SOR and NOR. Figure 8 shows the observed and simulated SOR and NOR for 15-31 October at the BNU site.The colors of those scatter points are associated with PM 2.5 concentrations.As simulated SO 2 were 3 times higher in Beijing compared with observations from the 34 local sites (Fig. 3) possibly due to an overestimation in the SO 2 emissions, the calculated SOR by the simulations would be artificially low due to the SO 2 overestimation.To correct this problem and to check the SOR differences among different simulations, we used the observed  SO 2 when calculating the SOR for the BASE, HET_BASE, and HET_EMIS simulations.The hypothesis is that the environment in Beijing is sulfur rich and the SO 2 overestimation would be corrected by reducing the SO 2 emissions in the model.It should be noted that reducing SO 2 emissions in the model not only reduces SO 2 concentrations but also aerosol sulfate concentrations which may lead to SOR decrease in the three simulations.This uncertainty should be considered when conducting the heterogeneous reaction coefficients tests in the future.From the observations we can see that SOR and NOR have a drastic increase in the 80-100 % RH range and are strongly correlated with high PM 2.5 concentrations.In the BASE simulation without the heterogeneous reactions, the corrected SOR and NOR are within a reasonable range when RH is below 60 %, but the corrected SOR and NOR are too low when RH is between 80-100 and 90-100 %, respectively.SOR and NOR are improved in HET_BASE simulation.In HET_EMIS simulation, the increase of NH 3 emissions and decrease of SO 2 emissions lead to increased NOR and PM 2.5 total mass.
From the 15-31 October period (Figs.6-7), we can see that the BASE simulation generally reproduced the PM 2.5 mass and also SNA species on relatively clean days, and the added heterogeneous reactions in HET_BASE helps to improve the SNA simulations (especially for sulfate) during the polluted events.Nevertheless, there are still large biases between 24 and 25 October.Given that the simulation meteorology is reasonable for these 2 days (Fig. 3), we assume that there might be two possible reasons for the biases: (1) the reaction rates are still too low due to the settings in the upper limit of SO 2 uptake coefficients and also the simulated RH, (2) the increase of precursor emissions in the upwind areas (southwest of Beijing) are not reflected in the model as small fires from autumn biomass burning are not updated in the emission inventory.Fire count from http://rapidfire.sci.gsfc.nasa.gov/cgi-bin/imagery/firemaps.cgi,however, did show intensive fire locations in southern Hebei, which is upwind of Beijing.Another factor is that the underestimations of PM 2.5 in Shijiazhuang and Baoding on 24-25 October are even larger.To examine our assumptions, several different scenarios (details in Sect.2.5) are simulated for these 2 days and the results are shown in Fig. 9.In the HET_MAX1 simulation, the reaction rates for the SO 2 -NO 2 -NO 3 heterogeneous reactions are fixed and triple those of the upper limits in HET_BASE.Without any emission increase, simulated sulfate almost reached the peak values on 24 October but was still underestimated on 25 October.A higher SO 2 heterogeneous reaction rate (7 times the upper limit in HET_BASE) and also doubled NH 3 emissions in HET_MAX2_EMIS2 enable the model to reach the sulfate peak on 25 October.For nitrate, the increase of SO 2 heterogeneous reaction rates in HET_MAX1 lead to lower nitrate concentrations even though the NO 2 -NO 3 heterogeneous reaction rates were also tripled.Only when NO x and NH 3 emissions in HET_MAX1_EMIS1 are increased by 50 % do simulated peak nitrate concentrations become comparable with observations on 24 and 25 October.Therefore, to reach the peak values of SNA aerosols on these 2 days, the sensitivity simulations suggest an increase of SO 2 heterogeneous reaction rates and NO x and NH 3 emissions are essential.For the best simulations (HET_MAX1_EMIS1 on 24 October and HET_MAX2_EMIS2 on 25 October), total PM 2.5 mass was improved by ∼ 100 µg m −3 compared to the BASE simulation.
In the three scenarios for 24-25 October, the SO 2 uptake coefficients were 3 and 7 times the values in HET-BASE.It is hard to justify whether those values are realistic.More rigorous tests should be conducted when additional observations become available.It is also difficult to determine whether the sulfate underestimation is only due to the missing of heterogeneous reactions.Zhang et al. (2015b) also emphasized that aqueous reactions of SO 2 by H 2 O 2 , O 3 , and NO 2 in cloud and on deliquescent aerosols would also help to improve the sulfate simulations significantly in GEOS-Chem.Following the study of Zhang et al. (2015b), we added the relevant aque-objective is to evaluate the capability of one such model, better understand the mechanisms that form sulfate, investigate the uncertainties associated with a set of heterogeneous chemical reactions, and improve the simulations of very high PM 2.5 concentrations during pollution episodes.The evaluations of meteorological parameters in the NCP show that the model is capable of capturing the temporal variations of boundary layer height as well as the low boundary layer heights during four pollution events.The deficiencies in meteorological forecasts were underestimates in relative humidity, especially during the most polluted days, and the overestimates in wind speed.While the default version of the CBMZ-MOSAIC mechanism available in the public version of WRF-Chem was able to simulate the high PM 2.5 concentrations (daily mean up to 200 µg m −3 ) for most of the cities in the NCP, the model severely underestimated the peak values (hourly mean greater than 400 µg m −3 ) in Beijing and at upwind sites (Baoding and Shijiazhuang).PM species observations at the BNU site in Beijing show that the sulfate-nitrate-ammonium aerosols (SNA) increases dramatically during these peak events and the increased sulfuroxidation rates and nitric-oxidation rates are strongly correlated with the high relative humidity (60-90 %) on those days.
The failure of the model to simulate the peak PM 2.5 concentrations is mainly due to the underestimation of SNA and secondary organic compounds.Analyses of the SNA underestimation revealed that missing SO 2 -NO 2 -NO 3 relevant heterogeneous reactions in the current aerosol scheme are likely important in China.Following the methodology in Zheng et al. (2015), the RH-dependent SO 2 , NO 2 , and NO 3 uptake heterogeneous reactions were added to the CBMZ-MOSAIC chemistry scheme.With the newly added reactions, the SNA simulations of the ratios of SNA in PM and the partitioning of sulfate and nitrate aerosols were improved on polluted days.However, there was still a 100 µg m −3 underestimation of SNA aerosols for the 24-25 October period when PM 2.5 concentrations were as high as 400-500 µg m −3 .Two possible explanations are proposed: (1) the RH-dependent reaction rates of those heterogeneous reactions, especially the SO 2 uptake reactions, on those peak days were not high enough either due to the underprediction in RH or the setting of the upper limits of uptake coefficients; (2) comparisons of modeled gas-phase precursors showed the possibility of precursor underestimation (especially NO x ) in the model.Although the two explanations cannot be proved definitively without additional observational evidence, sensitive simulations with increased heterogeneous reaction rates and increased NO x and NH 3 emissions show great improvement of SNA simulations in the model for the 24-25 October peak pollution events.In addition to the heterogeneous reactions, missing aqueous reactions may also play an important role in the sulfate underestimation in the model.
We conclude that RH in the 60-100 % range is a significant factor contributing to peak PM 2.5 values, especially dur-ing the heavily polluted days when sulfur and nitric oxidation rates almost doubled or tripled indicating the rapid heterogeneous reactions.With the underprediction in RH in this range and the corresponding low reaction rates, it is difficult for the model to reproduce the high concentrations of SNA.Data assimilation of meteorological variables, particularly humidity, might be useful from this point of view.Two other concerns should be addressed in future studies.First, the heterogeneous reaction rates applied in this study were based on the available reports in the literature, and the species evaluation was based on only one site in Beijing.Although the comparisons at one site are improved for the polluted days in Beijing, the sensitivity simulations on peak days indicated the current setting of the upper limit reaction rates might be still too low.Evaluations conducted using longer simulation periods of time and additional sites with observed speciated aerosol concentrations are needed to determine the appropriate heterogeneous reaction rates.Second, precursor emissions that are higher than available in the emission inventory due to special conditions might also be another factor contributing to the peak PM 2.5 events.Data assimilation techniques that use observed aerosol concentrations to constrain emission changes might also provide significant improvements in future studies.

Data availability
Data used in this publication can be accessed by contacting the authors Zhiquan Liu (liuz@ucar.edu) and Dan Chen (dchen@ucar.edu).
The Supplement related to this article is available online at doi:10.5194/acp-16-10707-2016-supplement.

Figure 1 .
Figure 1.The model domain and (left panel) and the NCP (right panel).Dots are the observational sites, where red denotes NCDC meteorological sites, blue denotes the Beijing Normal University (BNU) site, black denotes AQI national monitoring sites.Shaded backgrounds are model-simulated PM 2.5 concentrations (µg m −3 ) that illustrate the regional pollution on 10 October.

Figure 2 .
Figure 2. Spatial distribution of SO 2 (left), NO (middle), and NH 3 (right) emissions over the model domain (top panels) and in NCP (bottom panels).Units are in mol km −2 h −1 .

Figure 5 .
Figure 5.Time series of observed averaged pollutant concentrations at 34 sites in Beijing (from local monitoring system) along with the results from the BASE simulation.Units are mg m −3 for CO and µg m −3 for the other panels.

Figure 7 .
Figure7.Pie charts of observed and simulated mass fractions of PM 1 species on relatively clean days (observed PM 2.5 < 75 µg m −3 ) and polluted days (observed PM 2.5 ≥ 75 µg m −3 ) between 15 and 31 October.The units of total PM 1 concentrations (listed in panel titles) are µg m −3 .The species fractions (listed in panel labels) are in percentage (%).

Figure 8 .
Figure8.Observed (OBS) and simulated (BASE, HET_BASE, HET_EMIS) sulfur (SOR) and nitric (NOR) oxidation rates in PM 1 between 15 and 31 October.The x axis is the observed for left panels and simulated relative humidity for the other panels.Colors denote different PM 2.5 concentrations (µg m −3 ).

Table 2 .
Reactions and uptake coefficients added in this study.

Table 4 .
Statistics of meteorological simulations.