Severe winter haze days in the Beijing-Tianjin-Hebei region from 1985-2017 and the roles of anthropogenic emissions and meteorological parameters

We applied a global 3-D chemical transport model (GEOS-Chem) to examine the variations in the frequency and intensity in severe winter haze days (SWHDs) in BTH from 1985-2017 and quantified the roles of changes in anthropogenic emissions and/or meteorological parameters. Comparisons between the simulated SWHDs and those obtained from the observed PM2.5 concentrations and atmospheric visibility showed that the model can capture the spatial and temporal 15 variations of SWHDs in China; the correlation coefficient between the simulated and observed SWHDs is 0.98 at 161 grids in China. From 1985-2017, with changes in both anthropogenic emissions and meteorological parameters, the simulated frequency (total severe haze days in winter) and intensity (PM2.5 concentration averaged over severe haze days in winter) of SWHDs in BTH showed increasing trends of 4.5 days decade and 13.7 μg m decade, respectively. The simulated frequency exhibited fluctuations from 1985-2017, with a sudden decrease from 1992-2001 (29 days to 10 days) and a rapid 20 growth from 2003-2012 (16 days to 47 days). The sensitivity simulations indicated that variations in meteorological parameters played a dominant role during 1992-2001, while variations in both emissions and meteorological parameters were important for the simulated frequency trend during 2003-2012 (simulated trends were 27.3 days decade and 12.5 days decade owing to changes in emissions alone and changes in meteorology alone, respectively). The simulated intensity showed a steady increase from 1985-2017, which was driven by changes in both emissions and meteorology. The results of 25 this study have important implications for the control of SWHDs in BTH.


2018).
The Beijing-Tianjin-Hebei (BTH) region, which is one of the most densely populated economic zones in China, has considerably high level of PM 2.5 (Jiang, 2015;Chen and Wang, 2015).During severe winter haze days (SWHDs), the detected daily mean PM 2.5 concentrations reached 300-600 µg m -3 in BTH (Wang et al., 2014a;Li et al., 2018b).Exposure to these fine particulate matters detrimentally affects human health, causing problems from respiratory illnesses to cardiovascular diseases and premature death (Lelieveld et al., 2015;Zhang et al., 2017).High PM 2.5 levels during SWHDs also severely deteriorate ambient visibility, which endangers ground and air traffic and consequently disrupts economic activities (Wang et al., 2015;Gao et al., 2015a).Therefore, understanding the long-term variations in SWHDs is essential for air quality planning.
Previous studies on severe winter haze were mainly focused on a single episode, which reported the physical and chemical processes that led to SWHDs.Severe winter haze episodes usually occur under stagnant weather conditions and high anthropogenic emissions (Zhang et al., 2014;Zhang et al., 2015b;Li et al., 2017a).Stagnant weather conditions, characterized by weak surface wind speed, enhanced temperature inversion in the lower troposphere, and depressed planetary boundary layer height (PBLH) restrain the dispersion of pollutants (Liu et al., 2013;Zhao et al., 2013b).High relative humidity (RH) during SWHDs facilitates the rapid growth of secondary aerosols (Sun et al., 2014;Zhang et al., 2015b).Regional transport from upwind areas ( Zheng et al., 2015b;Sun et al., 2016;Ma et al., 2017) and the feedbacks between aerosols and meteorological variables are also important drivers for the formation of severe winter haze events (Gao et al., 2015b;Zhang et al., 2015a;Tie et al., 2017).The changes in large-scale circulation, such as the anomalous eastward extension of the Siberian high (Jia et al., 2015), weakened East Asian winter monsoon (Li et al., 2016b), and anomalous warm southerlies at 850 hPa (Cai et al., 2017) were also found to be conducive to the formation and maintenance of severe winter haze in BTH.
Some studies examined the changes in air quality in China during the past decades, but these studies were generally focused on the monthly or seasonal mean haze days.Because of the lack of long-term measurements of PM 2.5 concentrations in China, observed atmospheric visibility was often used as a proxy for PM 2.5 .Haze days were defined in previous studies as the days with atmospheric visibility of less than 10 km and relative humidity of less than 90 % (Chang et al., 2009;Yang et al., 2016).Ding and Liu (2014) examined haze days at 553 sites in China and reported that annual haze days averaged over these sites increased from 2-4 days in the 1960s to 11-16 days after 2005.It was also reported that haze days mainly increased in economically developed eastern China and decreased in northeastern and northwestern China.Yang et al. (2016), using visibility data, found that the winter haze days averaged over eastern China (196 sites, 105-122.5°E, 20-45° N) showed a large increasing trend of 2.6 days decade -1 during 1980-2012.Historical changes in air quality in China were also examined through modeling studies.Using the GEOS-Chem model, Jeong and Park (2013) reported a 73 % increase in simulated DJF mean SO 4 2--NO 3 --NH 4 + (SNA) concentration between 1985-1989 and 2002-2006 over East Asia (90-145° E, 20-45° N) and largely attributed this increase to changes in anthropogenic emissions.Also by using the GEOS-Chem, Yang percentile in a year) at 272 sites in China and compared these trends with the trends in median AEC (the 50 th percentile in a year).They showed that the extreme trends exceeded the median trends at most sites in the 1980s, the extreme trends were larger than the median trends in the south (south of 33° N) but smaller in the north (north of 33° N) in the 1990s, and then, most sites exhibited a faster increase in the median trends in the 2000s.These studies, performed through statistical methods based on visibility datasets, could not provide information about the historical changes in aerosol components in SWHDs and associated chemical/physical processes (transport or diffusion, chemical production or loss, wet and dry deposition).
Additionally, no quantitative information has ever been given regarding the separate roles of changes in anthropogenic emissions and meteorological parameters in historical changes of SWHDs in China.
In this study, we applied the GEOS-Chem model to simulate SWHDs over China for 33 winters from 1985 to 2017 and to quantify the roles of changes in anthropogenic emissions and meteorological parameters in the variations of simulated

Observed visibility
Observed atmospheric visibility data were obtained from the National Climatic Data Center (NCDC) Global Summary of Day (GSOD) database (https://www7.ncdc.noaa.gov/CDO/cdoselect.cmd),which has been widely used in previous studies to examine the haze pollution trend in China (Che et al., 2007;Chang et al., 2009;Deng et al., 2012;Yang et al., 2016;Li et al., 2016a The GEOS-Chem model includes fully coupled O 3 -NO x -hydrocarbon and aerosol chemistry with more than 80 species and 300 reactions (Bey et al., 2001;Park et al., 2004).The PM 2.5 components simulated in GEOS-Chem include sulfate (Park et al., 2004), nitrate (Pye et al., 2009), ammonium, black carbon and primary organic carbon (Park et al., 2003), mineral dust (Fairlie et al., 2007), and sea salt (Alexander et al. (2005).Aerosol thermodynamic equilibrium is computed by the ISORROPIA II package, which calculates the gas-aerosol partitioning of the sulfate-nitrate-ammonium system (Fountoukis and Nenes, 2007).Heterogeneous reactions of aerosols include irreversible absorption of NO 3 and NO 2 on wet aerosols (Jacob, 2000), hydrolysis of N 2 O 5 (Evans and Jacob, 2005), and the uptake of HO 2 by aerosols (Thornton et al., 2008).Wet deposition, including below-cloud washout, in-cloud rainout and scavenging in moist convective updrafts follows the scheme of Liu et al. (2001).Dry deposition is computed based on the resistance-in-series model from Wesely (1989) with a number of modifications (Wang et al., 1998).Considering that sea salt is not a dominant aerosol species in China and that the concentrations of mineral dust aerosols are generally low in winter (Duan et al., 2006;Zhao et al., 2013a), we calculated the PM 2.5 concentration in this study as the sum of the simulated mass of sulfate, nitrate, ammonium, BC, and OC.

Emissions
For anthropogenic emissions, global SO 2 , NO x , NH 3 , and CO emissions followed the EDGAR v4.(2011).For 2011-2017, scaling factors for all species followed the study of Zheng et al. (2018).Global biomass burning emissions were taken from GFED4 for 1997-2016 with a monthly temporal resolution (van der Werf et al., 2017).For the years out of the available range, emissions of the closest available year were used.For example, the biomass burning emissions of 1997 were used for 1985-1996.
Figure 1 shows the variations in wintertime emissions (both anthropogenic and natural emissions) of aerosols and their precursors over eastern China (105-122.5°E, 20-45° N) from 1985-2017.For the SO 2 , NO x , BC and OC emissions, all are a two-peak type with one peak occurring in the 1990s and another in approximately 2010 (2005 for SO 2 ).The relatively low emissions around 2000 were due to a slowdown in economic growth and a decline in fuel consumption (Streets et al., 2000;Hao et al., 2002;Streets and Aunan, 2005;Ohara et al., 2007).Emissions exhibited drastic increases in the beginning of the 21 st century and decreases in recent years as a result of the strictly enforced emission reduction measures during the 11 th and 12 th Five-Year Plan in China (Zheng et al., 2018).The NH 3 emissions continued to increase from 1985-2017 because of the steady growth in agricultural sources.Overall, relative to 1985, SO 2 , NO x , NH 3 , BC and OC emissions in 2017 changed by 11 %, 88 %, 81 %, 18 %, and -24 %, respectively.The variations in emissions used in this work are consistent with those from the studies of Jeong and Park (2013) and Yang et al. (2016).

Numerical simulations
Daily PM 2.5 concentrations for winters from 1985-2017 were simulated using the GEOS-Chem model driven by MERRA-2 meteorological fields.The following three simulations were conducted to examine the changes in SWHDs in China from 1985-2017 and to identify the relative roles of changes in anthropogenic emissions and meteorological parameters.The winter (December, January, February, DJF) of a specific year includes December of this year and January and February of the following year.

Definition of severe winter haze
To understand the decadal changes in SWHDs, we rely on the GEOS-Chem simulation and use the observed PM 2.5 concentration and atmospheric visibility for model evaluation.SWHDs must be defined for each of the datasets of observed PM 2.5 , observed atmospheric visibility, and simulated PM 2.5 .For observed PM 2.5 , following the study of Cai et al. ( 2017), a severe winter haze day is identified when the daily mean observed PM 2.5 concentration exceeds the 150 µg m -3 threshold because the Chinese government issues a "red alert" when the PM 2.5 concentration is forecasted to exceed 150 µg m -3 for 72 consecutive hours.
With respect to the definition of a severe winter haze day based on observed atmospheric visibility, a severe haze day is defined as a day with visibility of less than a visibility threshold (V t ) and relative humidity of less than 90 %.A statistical approach is used to obtain V t for each site of interest: 1) select the time period with available observations of both PM 2.5 and atmospheric visibility; 2) scatterplot the daily atmospheric visibility (Vis) versus daily mean PM 2.5 concentration (PM 2.5_obs ) for all samples in the time period selected; and 3) perform an exponential fit as obtain the V t that corresponds to a PM 2.5_obs of 150 µg m -3 .Notably, atmospheric visibility was observed manually before 2013 and has been observed automatically since 2013.To account for the discrepancies in data before and after 2013, two sets of Vt are obtained.Figure S1 provides an example of how to compute V t for Beijing.Observations of PM 2.5 and atmospheric visibility are available for 2009-2016.Using the statistical approach described above, V t values are calculated as 6.5 km over the manual period of 2009-2012 and 4.5 km over the automatic period of 2013-2016.Then, the V t values obtained from 2009-2012 were used to obtain SWHDs for the manually observed period of 1985-2012.
To identify SWHDs based on simulated PM 2.5 , a severe winter haze day is defined as that with a simulated daily mean PM 2.5 larger than a threshold concentration (M t ).Considering the biases in the model results caused by the representation of chemical and physical processes or the horizontal and vertical resolution in the model, Ct is also obtained from a statistical approach for each of the model grids with corresponding PM 2.5 observations: 1) select the time period with both observed and simulated PM 2.5 available; 2) for all samples selected in 1), calculate the mean bias (MB) between the daily simulated PM 2.5 ( !.!_!"# ) and daily mean observed PM 2.5 (PM 2.5_obs ) as  = , where N is total number of all samples; and 3) set M t = 150 µg m -3 + MB.Using this approach, M t is 150 µg m -3 if the simulated PM 2.5 has no bias, and M t is smaller (larger) than 150 µg m -3 if the simulated PM 2.5 concentrations have low (high) biases.Take the grid at Baoding (115.6°E, 39° N), one of the most polluted cities in BTH, as an example.From 2013-2017, when both simulated and observed PM 2.5 data are available, the simulated wintertime PM 2.5 from GEOS-Chem has a mean bias of -28.1 µg m -3 , and therefore, M t is 121.9 µg m -3 .As a result, in Baoding, the number of observed SWHDs is 210 days and that of the simulated SWHDs is 186 days during 2013-2017.Figure S2 displays the calculated M t values at 161 grids using the above statistical approach.These values were used to identify the SWHDs at each grid from 1985-2017.China.These years are chosen depending on the available years with measured PM 2.5 so that the model results can be evaluated by using observations.Simulated PM concentration reached 537 µg m -3 on December 25.However, the model captured only 3 SWHDs during this period, and the simulated PM 2.5 concentration was 37 µg m -3 on this day.Such under predictions of PM 2.5 concentrations during severe haze events were also reported by previous studies (Zhang et al., 2015b;Zhang et al., 2014;Wang et al., 2014b) and were attributed to the representation deficiency of the aqueous/heterogeneous processes in the models (Huang, 2014;Zheng et al., 2015a).For Shanghai and Chengdu, the model has high biases with NMBs of 18.6 % and 28.7 %, respectively.Over 1985-2017, while the SWHD intensity showed steady growth over the whole period, the changes in SWHD frequency exhibited three distinct stages.The frequency was quite stable, with a value of approximately 3.5 days per winter from 1985-2000, but the value increased rapidly from 2001-2011 and then decreased sharply during 2012-2017; this change can also be seen clearly in Fig. 7, which presents the spatial distributions of linear trends in SWHD frequencies during these three time periods.During 1985-2000, most grids showed no prominent trends, except for some grids in the Henan and Shanxi provinces as well as the Yangtze River Valley, where growth rates were in the range of 0.4-7.5 days decade -1 .Interestingly, negative trends were found over BTH during this period, with a maximum trend of -6.9 days decade -1 .From 2001-2011, the SWHD frequencies increased rapidly on a national scale.The largest increasing trends of more than 20 days decade -1 were centered on the North China Plain (NCP) and Fen-Wei Plain.Benefiting from the emission reduction actions, the frequencies in most grids began to experience a decline in 2012.Most prominent improvements occurred in the southern BTH and SCB and Shaanxi provinces, with decreasing trends of greater than -45 days decade -1 .

Decadal changes over BTH
As shown above, BTH merits special attention because the region has high DJF mean PM 2.5 concentrations, and this region has experienced the most frequently and intensely occurring severe winter haze over the years.To represent the SWHDs over BTH, a regional SWHD is defined here as the day when severe winter haze (defined in Sect.2.4) occurs simultaneously in more than one third of the grids in a region.Therefore, for the BTH region that include 18 grids, a regional SWHD is identified when 6 grids or more report severe winter haze at a certain time.The intensity of this regional SWHD is calculated as the mean PM 2.5 concentration of the grids with SWHDs on that day.
Figure 8 shows the temporal evolution of the simulated frequency and intensity, as well as the concentrations of regional SWHD aerosol components in BTH from 1985-2017.Differing from the mean values averaged over 161 grids shown in Fig. 6c, the frequency of regional SWHD in BTH exhibited a much higher level before 2000, indicating serious air pollution by that time.Specifically, the frequency increased from 14 days in 1985 to 29 days in 1992, followed by a sudden decrease to 10 days in 2001, a rapid growth to a peak of 47 days in 2012 and a steep decrease thereafter to 15 days in 2017, which coincides with the temporal variations presented in previous studies for haze days (Chen and Wang, 2015), aerosol extinction coefficient (Li et al., 2016a) and PM 2.5 concentration (Yang et al., 2018) in the North China Plain (NCP).From 1985-2017, SWHDs in BTH exhibited large increasing trends of 4.5 days decade -1 in frequency and 13.7 µg m -3 decade -1 in intensity.In comparison, the seasonal mean PM 2.5 concentrations showed a smaller trend of 10.1 µg m -3 decade -1 (Supplementary Fig. S3), implying a more serious problem with severe winter haze pollution.When considering each species, nitrate aerosols exhibited the largest increasing trend of 15.7 µg m -3 decade -1 , which was followed by ammonium (4.5 µg m -3 decade -1 ).Sulfate and BC made little contribution to the increasing trend, while OC showed the opposite trend of -6.8 µg m -3 decade -1 , which coincided with the emission changes during this period (Fig. 1).Notably, the simulations in this study do not include secondary organic aerosols; therefore, the concentrations of organic aerosols may have been underestimated.In eastern China, the GEOS-Chem model overestimates nitrate and underestimates sulfate in winter (Wang et al., 2013).

Key processes for the occurrence of SWHDs over BTH
Concentrations of aerosols are jointly determined by emissions, transport, chemical reactions and deposition.To further isolate and understand the roles of each process in the formation of SWHDs over BTH, a process analysis (PA) was performed for the CTRL simulation from 1985-2017.PA has been widely used in previous studies to identify the relative importance of atmospheric processes in specific pollution episodes (Goncalves et al., 2009) or annual (Zhang et al., 2009) to decadal simulations (Mu and Liao, 2014;Lou et al., 2015) of air pollutants.
For each aerosol species, six processes were diagnosed at every time step, consisting of primary species emissions, horizontal and vertical transport, chemical reactions, cloud processes, dry deposition, and PBL mixing, which altogether determined the mass balance of aerosols.Here, the cloud processes include the scavenging of soluble tracers and mixing due to cloud convections.PBL mixing refers to the mass flux brought by turbulence within the planetary boundary layer.PA was carried out for the BTH box (from surface to the 11 th vertical model layer that corresponds to approximately 850 hPa) on the basis of daily values.The contribution of each process to the formation of SWHDs was calculated by the following equation: where n is the number of processes (n = 6 as described above), PC SM_i is the daily mean mass flux of process i from 1985-2017, and PC SWH_i is the daily mean mass flux of process i averaged over SWHDs during 1985-2017.Note that the sum of abs(%PC i ) is 100 %. and local emissions (5.5 Gg day -1 ).Net transport had an effect of -12.1 Gg day -1 , exporting the particles to the downwind areas as well as to the upper troposphere.The contributions of turbulent diffusion (PBL mixing), dry deposition, and cloud processes were small, with values of -0.8 Gg day -1 , -1.8 Gg day -1 , and -0.9 Gg day -1 , respectively.For all days during winters from 1985-2017, the sum of all processes had a very small value of 0.004 Gg day -1 .During SWHDs in the selected time period, the outflow of particles decreased greatly, from 12.1 Gg day -1 in mean state to 4.8 Gg day -1 , which held more particles in BTH.More secondary aerosols were produced through chemical reactions, with an increase from a mean state of processes, dry deposition and PBL mixing to the SWHDs were 65.3 %, 17.6 %, -7.5 %, -6.4 % and 3.2 % (Table 1), respectively, indicating that transport was the most important factor that facilitated the occurrence of SWHDs in BTH.
Notably, the relative contribution of emissions is zero because the emissions in the model were monthly values, which were the same for either polluted or clean days.During the SWHDs, the relative contribution of transport is the largest for all aerosol species except sulfate (Table 1).For sulfate aerosols, chemical reactions had a large relative contribution of 49.3 %, which is comparable to the transport contribution (39.5 %).This finding highlighted the importance of local chemical production of sulfate during SWHDs.With regard to the carbonaceous aerosols, transport was the dominant process with a relative contribution of 83.2 % since carbonaceous aerosols are assumed to be chemically inert in the model.For nitrate and ammonium, transport had relative contributions of 62.4 % and 61.3 %, respectively.
To further explore the wind component (zonal, meridional or vertical) that contributed most to the transport process, mass fluxes brought by winds in three directions (east-west, north-south, up-down) were also determined for the BTH box from 1985-2017.Then, relative contributions were calculated following Eq.( 1) except that n is the number for the three directions.
The north-south (NS) wind was calculated to have a relative contribution of 93 % and hence a dominant role (pie chart in Fig. 9).This result agrees with previous studies (Li et al., 2016b;Cai et al., 2017), which found that anomalous southerlies in the lower troposphere favor the formation and accumulation of PM 2.5 during severe haze events.

Roles of anthropogenic emissions and meteorological parameters
Figure 10a shows the time series of regional SWHDs frequencies in BTH from 1985-2017 based on three simulations (CTRL, EMIS, and MET).In the CTRL simulation, with changes in both emissions and meteorology, the SWHD frequency in BTH increased from 1985-2017, with a linear trend of 4.5 days decade -1 (as described in Sect.4.2).During 1985-2017, the simulated SWHD frequency exhibited an increasing trend of 5.9 days decade -1 in EMIS with changes in emissions alone but no significant trend was seen in the MET simulation with variations in meteorological parameters alone.Considering that the frequency trend is dependent on the selected simulation years, we show in decade -1 , respectively, indicating that the variations in meteorological parameters dominated the decreasing trend over this period.The trend of -3.3 days decade -1 from 1992-2001 in EMIS was consistent with the emission reductions over this period (Fig. 1).During 2003-2012 when the increasing frequency trend (23.6 days decade -1 ) was the largest in the CTRL simulation, the EMIS and MET trends were 27.3 days decade -1 and 12.5 days decade -1, respectively, reflecting that both the increases in emissions and changes in meteorological parameters contributed to the increasing trend in CTRL, although the role of meteorology was smaller compared to that of emissions.Notably, the decrease in frequency over the 1990s in MET was mainly caused by the interdecadal variations in atmospheric circulations (Chen and Wang, 2015).The researchers found that compared to the meteorological conditions from 1993-2001, those from 1984-1992 and 2002-2010 were more favorable for the occurrence of severe winter haze events over North China because of the weaker wind speeds and increased moisture in the lower troposphere, weaker East Asian trough in the midtroposphere, and greater northward shift of the East Asian jet stream in the upper troposphere.
Figure 11a shows the time series of the regional SWHD intensities in BTH from 1985-2017 based on the CTRL, EMIS, and MET simulations.From 1985-2015, the simulated SWHD intensity exhibited increasing trends of 13.7 µg m -3 decade -1 , 6.2 µg m -3 decade -1 , and 8.0 µg m -3 decade -1 due to changes in both emissions and meteorology, changes in emissions alone, and changes in meteorological parameters alone.Therefore, both anthropogenic emissions and meteorological parameters were major factors that led to high PM 2.5 concentrations in the SWHDs in BTH.The intensity in CTRL exhibited an overall increasing trend except for a small decrease from 1992-2001 (Fig. 11b), due to the reasons mentioned above.
Additionally, we examine the interannual variations in SWHDs in BTH.The 9-point weighted running mean is employed for the frequencies in Fig. 10a and intensities in Fig. 11a, respectively, to remove the fluctuations in periods of more than 9 years while reserving the interannual anomalies.The weighting coefficients are based on the Lanczos filter (Duchon, 1979).
Figure 12 shows the interannual variations in frequencies and intensities from the three simulations (CTRL, EMIS, and MET) from 1989-2013.The interannual frequency variation in MET generally resembles that from the CTRL simulation, with a high correlation coefficient R of 0.80 between the simulations, while the R between the interannual variations of CTRL and EMIS is only 0.05 (Fig. 12a).Similarly, for intensity, the interannual variation from the MET exhibited the same peaks and troughs as that in the CTRL simulation, with a high R of 0.90 between the simulations (Fig. 12b).Therefore, the interannual SWHD variations in BTH were mainly driven by variations in meteorological parameters.To quantify the interannual variation caused by the meteorological variations, the relative interannual change (%IC) was calculated by the following equation: where MET_IA i and MET_IA i-1 are the interannual anomalies in the MET simulation in years i and i- This highlights the importance of variations in meteorological parameters for policy makers when actions are taken to improve air quality in the short term.

Conclusions
We investigated the changes in SWHDs (days with observed PM 2.5 concentrations larger than 150 µg m -3 ) from 1985-2017 and quantified the respective roles of changes in anthropogenic emissions and meteorological parameters in the simulated SWHD changes using the nested-grid version of the GEOS-Chem model.
The simulated SWHDs were evaluated using the observed PM  alone and meteorology alone both enhanced the intensity growth with trends of 6.2 µg m -3 decade -1 and 8.0 µg m -3 decade -1 , respectively.These results suggest that meteorological parameters must be considered for the control of SWHDs.
Note that a number of factors contribute to the uncertainties in our results.We obtained long-term SWHDs by using the observed atmospheric visibility, which were influenced by the changes in the way of observation or relocation of observation site, especially before 2013 when observations were carried out manually.The uncertainties with emissions inventories of aerosols and aerosol precursors might lead to uncertainties in our simulated decadal changes in aerosol concentrations.
Furthermore, we only accounted for anthropogenic aerosols in our simulations; natural aerosols, such as secondary organic aerosol, mineral dust, and biomass burning aerosols, contributed to the observed changes in SWHDs and should be considered for in model simulation.These issues need to be examined in future studies of severe haze trend over eastern China.
Tables Table 1.Relative contributions of each of the five processes (transport, chemistry, cloud, dry, and PBL mixing) to SWHD formation in BTH (unit: %).The process analyses are carried out for both PM 2.5 and its components, following Eq.( 1) described in Sect. 5.
et al. (2016) obtained an increasing trend in wintertime PM 2.5 concentration of 10.5 (±1.5) µg m -3 decade -1 averaged over Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2019-306Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 30 April 2019 c Author(s) 2019.CC BY 4.0 License.eastern China (105-122.5°E, 20-45° N) during 1985-2005 and found that the changes in anthropogenic emissions and meteorological parameters over the studied period contributed 83 % and 17 % to the increasing trend in PM 2.5 , respectively.To date, only a few studies have examined historical changes in SWHDs in China using observed atmospheric visibility.Fu et al. (2013) defined an extremely poor visibility event as a case with atmospheric visibility of less than 5 km for at least 3 days and found that the annual frequency of such event averaged over large cites in China (142 cities with populations of more than one million) increased from 1960 to 2009 with a trend of 0.68 events decade -1 .Using the visibility-converted daily aerosol extinction coefficient (AEC) dataset, Li et al. (2018a) examined the long-term trends in extremely high AEC (the 95 th Figure 2a presents the simulated seasonal (DJF) mean surface-layer PM 2.5 concentrations averaged from 2013-2017 over

Figure
Figure4ashows the simulated frequencies of SWHDs at 161 model grids averaged from 2013-2017 over China.Consistent with the spatial distribution of seasonal mean PM 2.5 in Fig.2a, simulated SWHDs occurred most frequently in BTH.The frequency averaged over BTH was 19.4 days, which is more than twice the mean value of 8.7 days for all 161 grids.The highest frequencies of 31-39 days were located in southern BTH, indicating that more than one third of the wintertime days experienced severe PM 2.5 pollution from 2013-2017.For comparison, the observed frequencies of SWHDs, averaged over the same years, are presented in Fig.4b.The SWHDs were observed most frequently in BTH, and the horizontal

Figure
Figure4dshows the spatial distribution of simulated SWHD intensities.For each grid, the intensity is calculated as the average of the daily mean PM 2.5 concentrations over all SWHDs during 2013-2017.The simulated mean PM 2.5 concentrations of SWHDs were high over the BTH, Shandong and Henan provinces and the Sichuan Basin.The average SWHD PM 2.5 concentration was 232 µg m -3 over BTH, which is much higher than the mean value of 207 µg m -3 for all grids.The observed intensities of SWHDs over the same period are also displayed in Fig.4e.The observed high PM 2.5 values of SWHDs are centered over BTH, which is generally reflected in the simulation.The model results have low biases in the Shanxi province and Inner Mongolia and high biases in the Sichuan Basin and Shandong and Henan provinces.Accounting for all model grids in Fig.4d and 4e, the scatterplot (Fig.4f) shows that simulated intensities of SWHDs have an NMB of 1.3 %.Over BTH, simulated intensities of SWHDs have a small NMB of 0.3 %.We also evaluate the model's performance to reproduce the long-term variation in SWHD frequency.To evaluate the historical changes, the SWHDs estimated from the observed atmospheric visibility are used.Beijing is the only grid with observed PM 2.5 concentrations available before 2013 for establishing the visibility threshold for identifying visibility-based SWHDs over 1985-2012 (as described in Sect.2.4).Figure5compares the temporal evolution of the SWHD frequency at the Beijing grid obtained from the GEOS-Chem simulation, observed PM 2.5 (as described in Sect.3.1), and observed atmospheric visibility (as described in Sect.2.2).During 2009-2016, when all three datasets are available, frequencies obtained from the simulation and atmospheric visibility capture the interannual variation (peaks and troughs) in observed frequency well.The model underestimates the frequency with an NMB of -4.9 %, and the visibility result has an NMB of -6.6 %.Over the entire 1985-2017 period, when both observed visibility and simulation are available, the correlation coefficient R is 0.48.The frequency obtained from the atmospheric visibility exhibited a decreasing trend of -4.2 days decade -1 during 1985-2002 and an increasing trend of 1.3 days decade -1 during 2003-2017; the model partly reproduces these trends with tendencies of -5.2 days decade -1 from 1985-2002 and 3.3 days decade -1 from 2003-2017.However, a distinct low bias was found in the simulation before 2003, which may be attributed to the inconsistency problem in the atmospheric visibility dataset obtained from the NCDC.In general, the CTRL simulation can reasonably well capture the spatial distributions and historical changes in SWHD frequencies and intensities in China.

Figure 9
Figure 9 presents the average daily mass fluxes of PM 2.5 and its components in the BTH box from each of the six processes over all days during winters from 1985-2017 (stripped bars) and over the SWHDs during this period (solid bars).Averaged over all days during winters from 1985-2017, the PM 2.5 in BTH was mainly contributed by chemical reactions (10.1 Gg day -1 ) Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2019-306Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 30 April 2019 c Author(s) 2019.CC BY 4.0 License.10.0 Gg day -1 to 12.0 Gg day -1 .The suppressed PBL height further restricted the PM 2.5 diffusion with reduced aerosol outflow from 0.8 Gg day -1 to 0.4 Gg day -1 .Because of the high PM 2.5 concentrations during the SWHDs, cloud processes and dry deposition removed more aerosols compared to the mean state.Overall, the sum of all processes had a large positive value of 8.1 Gg day -1 , leading to the accumulation of PM 2.5 in BTH.The relative contributions of transport, chemistry, cloud Fig. 10b-10d a more comprehensive analysis of frequency trends over different selected periods, with a minimum duration of 10 years.Here, the x-axis indicates the start year, and the y-axis indicates the number of years since the start year during which period the trend is calculated.As shown by the trends in the 10 years since the start year, the CTRL frequency experienced a moderate increase in the late 1980s, a significant decrease beginning in 1991, and a large increase beginning in 1997.Similar patterns are also found in the EMIS and MET simulations, although different magnitudes are simulated.During 1992-2001 when the decreasing trend (-20.1 days decade -1 ) was the largest in the CTRL simulation, the trends in EMIS and MET were -3.3 days decade -1 and -20.1 days Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2019-306Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 30 April 2019 c Author(s) 2019.CC BY 4.0 License.
2.5 concentrations and atmospheric visibility.The GEOS-Chem model captured the daily variations in PM 2.5 concentrations fairly well.The correlation coefficient (R) and normalized mean bias (NMB) between the simulated and observed wintertime daily PM 2.5 concentrations were 0.61 and -9.2 % in Beijing from 2009-2016, 0.70 and 18.6 % in Shanghai from 2013-2017, and 0.58 and 28.7 % in Chengdu from 2013-2017.The model also captured the SWHD frequency (total severe haze days during winter) and intensity (PM 2.5 concentration averaged over severe haze days during winter); compared to the observed SWHD frequency and intensity at 161 grids in China from 2013-2017, the simulated frequency had an R of 0.98 and NMB of -12.3 %, and the simulated intensity had an NMB of 1.3 %.The Beijing-Tianjin-Hebei region (BTH) had the highest SWHD frequency and intensity based on the observed PM 2.5 concentrations from 2013-2017, with a frequency of 21.2 days yr -1 and a mean intensity of 231.6 µg m -3 .Historically, from 1985-2017, with changes in emissions and meteorology, the simulated SWHD frequency in BTH increased with a linear trend of 4.5 days decade -1 , and the simulated intensity increased by 13.7 µg m -3 decade -1 , which was larger than the increase in the seasonal mean PM 2.5 concentration (10.1 µg m -3 decade -1

Figure 2 .
Figure 2. (a) Spatial distributions of simulated (CTRL, shades) and observed (CNEMC, dots) seasonal (DJF) mean surface-layer PM 2.5 concentrations (µg m -3 ) averaged over 5 winters (2013-2017).(b) Scatterplot of simulated versus observed DJF mean PM 2.5 concentration (µg m -3 ) averaged from 2013-2017 for 161 grids in (a), where the green grids are 18 grids located in the BTH (Beijing-Tianjin-Hebei) region.Also shown in (b) are the y=x line (dashed), linear fit (solid line and equation) and values of r and NMB.Here, r is the correlation 5

Figure 3 .
Figure 3. Simulated (purple solid line, from CTRL simulation) and observed (green solid line, from the U.S. Embassy and CNEMC) daily mean surface-layer PM 2.5 concentrations (µg m -3 ) for grids of (a) Beijing, (b) Shanghai and (c) Chengdu over the period when observations are available.Also shown are the threshold of PM 2.5 = 150 µg m -3 (red dashed line), the correlation coefficients (R) and NMB values between observed and simulated daily mean PM 2.5 concentrations.

Figure 4 .
Figure 4. Spatial distributions of (a) simulated and (b) observed frequencies (days) of SWHDs averaged from 2013-2017 at 161 model grids in China and (c) scatterplot of simulated versus observed results in (a) and (b); (d-f) are the same as (a-c) except that the values are SWHD intensities (µg m -3 ).Also shown in the scatterplots are y=x lines (dashed), linear fits (solid line and equation), correlation coefficients (R) and NMB values.

Figure 6 .
Figure 6.Spatial distributions of SWHD (a) frequencies (days) and (b) intensities (µg m -3 ) at 161 model grids averaged from 1985-2017.Temporal evolutions of SWHD (c) frequency (days) and (d) intensity (µg m -3 ) in China from 1985-2017.In (c) and (d), solid lines are the averages of all samples (161 model grids), and the shades indicate the range between the first and third quartiles of the samples.

Figure 7 .
Figure 7. Linear trends of simulated frequencies (days decade -1 ) of SWHDs at 161 model grids over the three periods of (a) 1985-2000, (b) 2001-2011 and (c) 2012-2017.Large dots indicate statistical significance above the 95 % confidence level, while smaller ones are not statistically significant.

Figure 9 .
Figure 9. Daily mean mass fluxes of aerosols by physical/chemical processes (bottom axis: emission, transport, chemistry, cloud processes, dry deposition, diffusion and the sums) in the BTH box.Two kinds of averages are given here: averaged over all days during winters from 1985-2017 (PC SM , left stripped bars) and averaged over all SWHDs during this period (PC SWHD , right solid bars).Black dots give the relative contributions (%PC, right axis) of each of the six processes i (unit: %) for the SWHD formation when compared with mean values.

Figure 12 .
Figure 12.Interannual variations in SWHD (a) frequencies (days) and (b) intensities (µg m -3 ) in BTH through a nine-year weighted moving average method.The 9-point high-pass Lanczos filter was used as the weighting coefficients, which effectively removes fluctuations with periods of more than 9 years and reserves the interannual anomalies.Also shown are the correlation coefficients that were calculated between the interannual signals of the CTRL and MET simulations (green R) and between the results of the CTRL and 5
2.5 concentrations were high over Eastern China, covering important economic zones including BTH, the Sichuan Basin and the Yangtze River Delta (YRD).Simulated highest concentrations occurred in southern BTH, with seasonal mean PM 2.5 in the range of 135-151 µg m -3 .For model evaluation, observed PM 2.5 concentrations from the CNEMC dataset (averaged into 161 model grids as described in Sect.2.1) are shown in Fig.2a.The simulated spatial patterns of PM 2.5 agree well with the measurements, with a high correlation coefficient (R) of 0.82 between the observed and simulated seasonal mean PM 2.5 values.Figure2bshows the scatterplot of simulated versus observed seasonal mean PM 2.5 concentrations at all 161 model grids shown in Fig.2a.The model overestimates PM 2.
5 concentrations with a normalized mean bias (NMB) of +8.8 % for all grids and an NMB of +4.3 % for BTH. Figure 3 compares the time series of simulated and observed daily mean surface-layer PM 2.5 concentrations at the Beijing, Shanghai and Chengdu grids, which represent the three most polluted regions of BTH, YRD and the Sichuan basin, respectively.For Beijing (116.25°E,40° N), observations at the U.S. embassy site were used for comparison from 2009-2017.For Shanghai (121.25°E,31° N) and Chengdu (103.75°E,30.5°N), observations from the gridded CNEMC dataset were available from 2013-2017.The model generally captures the daily variations (peaks and troughs) in the observed PM 2.5 concentrations, with R values of 0.61, 0.70 and 0.58 for Beijing, Shanghai and Chengdu, respectively.The model has a lowbias in Beijing with an NMB of -9.2 % and is unable to predict the maximum PM 2.5 concentration in some cases.For example, observations show that during December 19-26 th , 2015, there were 7 SWHDs in Beijing, and the highest PM 2.5 1,is the original frequency or intensity in the CTRL simulation in year i-1.The calculated mean %IC was 26 % for frequency from 1989-2013, and 5.3 % for intensity over the same period.Here, the %IC indicates that on average, meteorological variations alone can lead to a change of ±26 % in frequency and ±5.3 % in PM 2.5 concentrations for the SWHDs in BTH interannually.Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2019-306Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 30 April 2019 c Author(s) 2019.CC BY 4.0 License.
).Further process analysis was carried out to identify the key processes for SWHD occurrences in BTH.The comparisons of the six processes (emission, transport, chemistry, cloud processes, dry deposition, and PBL mixing) averaged over the SWHDs with those averaged over all days during winters from 1985-2017 indicate that transport (outflow of particles) had the largest change from -12.1 Gg day -1 in the mean state to -4.8 Gg day -1 in SWHDs, which had the largest contribution to SWHD formation.Transport, chemistry, cloud processes, dry simulated frequency trends were23.6,27.3, and 12.5 days decade -1 owing to the changes in both emissions and meteorology, emissions alone, and meteorology alone, respectively, highlighting the contributions of both emissions and meteorological parameters.Interestingly, the simulated SWHD intensity in BTH increased steadily from 1985-2017; variations in emissions Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2019-306Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 30 April 2019 c Author(s) 2019.CC BY 4.0 License.