Severe winter haze days in the Beijing–Tianjin–Hebei region from 1985 to 2017 and the roles of anthropogenic emissions and meteorology

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 Beijing–Tianjin–Hebei (BTH) from 1985 to 2017 and quantified the roles of changes in anthropogenic emissions and/or meteorological parameters. Observed SWHDs were defined as the days with daily mean PM_2.5 concentration exceeding 150 µg m⁻³, and simulated SWHDs were identified by using the same threshold but with adjustment on the basis of simulation biases. Comparisons between the simulated SWHDs and those obtained from the observed PM_2.5 concentrations and atmospheric visibility showed that the model can capture the spatial and temporal variations in SWHDs in China; the correlation coefficient between the simulated and observed SWHDs is 0.98 at 161 grids in China. From 1985 to 2017, with changes in both anthropogenic emissions and meteorological parameters, the simulated frequency (total severe haze days in winter) and intensity (PM_2.5 concentration averaged over severe haze days in winter) of SWHDs in BTH showed increasing trends of 4.5 d per decade and 13.5 µg m⁻³ per decade, respectively. The simulated frequency exhibited fluctuations from 1985 to 2017, with a sudden decrease from 1992 to 2001 (29 to 10 d) and a rapid growth from 2003 to 2012 (16 to 47 d). 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 and 12.5 d per decade owing to changes in emissions alone and changes in meteorology alone, respectively). The simulated intensity showed a steady increase from 1985 to 2017, which was driven by changes in both emissions and meteorology. Process analysis on all SWHDs during 1985–2017 indicated that transport was the most important process for the formation of SWHDs in BTH with a relative contribution of 65.3 %, followed by chemistry (17.6 %), cloud processes (−7.5 %), dry deposition (−6.4 %), and planetary boundary layer (PBL) mixing (3.2 %). Further examination showed that SWHDs exhibited large interannual variations in frequency and intensity, which were mainly driven by changes in meteorology. The results of this study have important implications for the control of SWHDs in BTH.

1 Introduction

Due to rapid industrialization and urbanization, China has suffered from frequent severe haze episodes in recent years (Cai et al., 2017; K. Li et al., 2018; Zhu et al., 2018), especially during winter (Chen and Wang, 2015; Wu et al., 2017; Wang et al., 2018). The Beijing–Tianjin–Hebei (BTH) region, which is one of the most densely populated economic zones in China, has considerably high levels of PM_2.5 (Jiang et al., 2015; Chen and Wang, 2015). During severe winter haze days (SWHDs), the detected daily mean PM_2.5 concentrations reached 300–600 µg m⁻³ 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; L. Zhang et al., 2015; J. Li et al., 2017). 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; L. Zhang et al., 2015). Regional transport from upwind areas (G. J. Zheng et al., 2015; 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; B. Zhang et al., 2015; Tie et al., 2017). Changes in large-scale circulation, such as the anomalous eastward extension of the Siberian high (Jia et al., 2015), weakened East Asian winter monsoon (Q. Li et al., 2016), 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 to 4 d in the 1960s to 11–16 d 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 d per decade 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 ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$–${\mathrm{NO}}_{\mathrm{3}}^{-}$–${\mathrm{NH}}_{\mathrm{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. By also using the GEOS-Chem, Yang et al. (2016) obtained an increasing trend in wintertime PM_2.5 concentration of 10.5 (±1.5) µg m⁻³per decade averaged over 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 d and found that the annual frequency of such an event averaged over large cites in China (142 cities with populations of more than 1 million) increased from 1960 to 2009 with a trend of 0.68 events per decade. Using the visibility-converted daily aerosol extinction coefficient (AEC) dataset, Li et al. (2018a) examined the long-term trends in extremely high AEC (the 95th percentile in a year) at 272 sites in China and compared these trends with the trends in median AEC (the 50th 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 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 in simulated SWHDs over BTH. The observed datasets, GEOS-Chem model and numerical experiments, and the definition of SWHDs are described in Sect. 2. In Sect. 3, simulated PM_2.5 and SWHDs as well as a model evaluation are presented. In Sect. 4, the simulated changes in SWHDs over China and the BTH region from 1985 to 2015 are presented. The key processes that led to SWHDs in BTH are identified by process analysis in Sect. 5. In Sect. 6, the relative roles of changes in anthropogenic emissions and meteorological parameters in SWHDs in BTH are investigated over the studied period.

2 Methods

2.1 Observed PM_2.5 concentrations

Hourly surface PM_2.5 measurements were mainly obtained from the China National Environmental Monitoring Centre (CNEMC, 2019, http://www.cnemc.cn). The monitoring network expanded from 670 sites in 2013 to 1600 sites in 2018, covering approximately 370 cities nationwide. The 670 sites set up in 2013 have obtained continuous measurements since 2013, and 79 sites are located in BTH. Daily mean PM_2.5 from these 670 sites is used in this study and averaged to the 0.5^∘ latitude × 0.625^∘ longitude MERRA-2 grid (corresponding to 161 model grids) to examine observed SWHDs and for model evaluation. For the BTH region, hourly PM_2.5 concentrations that were observed at the US embassy in Beijing (US embassy, 2019, http://www.stateair.net/web/post/1/1.html) are also used, which are available from 2009 to June 2017.

2.2 Observed visibility

Observed atmospheric visibility data were obtained from the National Climatic Data Center (NCDC) Global Summary of Day (GSOD) database (NCDC, 2019, 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; J. Li et al., 2016). The database includes daily observations from 1985 to 2018 at 346 sites in China. The Beijing site is located at 39.56^∘ N and 116.17^∘ E.

2.3 GEOS-Chem model and numerical experiments

2.3.1 Model description

We simulated PM_2.5 concentrations using the GEOS-Chem model (version 11-01; Harvard University, 2019, http://acmg.seas.harvard.edu/geos/) driven by MERRA-2 assimilated meteorological data (Gelaro et al., 2017) from NASA's Global Modeling and Assimilation Office (GMAO). The nested-grid capacity of GEOS-Chem over Asia (11^∘ S–55^∘ N, 60–150^∘ E) was used with a horizontal resolution of 0.5^∘ latitude by 0.625^∘ longitude and 47 vertical layers up to 0.01 hPa. Boundary conditions for gases and aerosols were updated every 3 h from the coupled global GEOS-Chem simulations performed at a $\mathrm{2}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ horizontal resolution.

The GEOS-Chem model includes fully coupled O₃–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₃ and NO₂ on wet aerosols (Jacob, 2000), hydrolysis of N₂O₅ (Evans and Jacob, 2005), and the uptake of HO₂ 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, black carbon (BC), and organic carbon (OC).

2.3.2 Emissions

For anthropogenic emissions, global SO₂, NO_x, NH₃, and CO emissions followed the EDGAR v4.2 inventory (EC-JRC/PBL, 2011, http://edgar.jrc.ec.europa.eu/overview.php?v=42, last access: 15 August 2019). Global BC and OC emissions were taken from the BOND inventory (Bond et al., 2007). Global non-methane volatile organic compound (NMVOC) emissions followed the RETRO inventory (ftp://rain.ucis.dal.ca/ctm/RETRO, last access: 15 August 2019), where C₂H₆ and C₃H₈ were replaced with emissions from Xiao et al. (2008). In the East Asian domain, MIX 2010 (M. Li et al., 2017) was the baseline anthropogenic inventory, and annual scaling factors were applied for other years during 1985–2017. Over 1985–2009, scaling factors for all species were derived from the EDGAR v4.3 inventory (http://edgar.jrc.ec.europa.eu/overview.php?v=431, last access: 15 August 2019), except those for SO₂ during 1996–1999, which were taken from Lu et al. (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 to 2017. For the SO₂, 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₂). 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 21st century and decrease in recent years as a result of the strictly enforced emission reduction measures during the Eleventh and Twelfth Five-Year Plans in China (Zheng et al., 2018). The NH₃ emissions continued to increase from 1985 to 2017 because of the steady growth in agricultural sources. Overall, relative to 1985, SO₂, NO_x, NH₃, 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).

Figure 1 (a) Variations in total emissions (anthropogenic plus natural emissions) of sulfur dioxide (SO₂, Tg S DJF⁻¹), nitrogen oxide (NO_x, Tg N DJF⁻¹), ammonia (NH₃, Tg N DJF⁻¹), black carbon (BC, Tg C DJF⁻¹), and organic carbon (OC, Tg C DJF⁻¹) over eastern China (105–123^∘ E, 20–45^∘ N) during December–January–February (DJF) from 1985 to 2017 and (b) the percentage changes relative to 1985 values.

2.3.3 Numerical simulations

Daily PM_2.5 concentrations for winters from 1985 to 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 to 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.

1. CTRL is the control simulation with variations in meteorological parameters, anthropogenic emissions, and biomass burning emissions from 1985 to 2017.

2. EMIS is the simulation with changes in anthropogenic and biomass burning emissions over 1985–2017, while the meteorological fields were fixed at the 1985 levels. The aim of this simulation is to quantify the impacts of changes in emissions on SWHDs during 1985–2017.

3. MET is the simulation with changes in meteorological fields over 1985–2017, while anthropogenic and biomass burning emissions were fixed at the 2015 levels. This simulation is set to examine the impacts of changes in meteorological parameters on SWHDs during 1985–2017.

2.4 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⁻³ threshold because the Chinese government issues a “red alert” when the PM_2.5 concentration is forecasted to exceed 150 µg m⁻³ 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) make a scatterplot of 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 Vis $=C_{\mathrm{1}}+C_{\mathrm{2}}\cdot \exp (C_{\mathrm{3}}\cdot {\mathrm{PM}}_{\mathrm{2.5}\mathrm{\_}\mathrm{obs}}$) and obtain the V_t that corresponds to a PM_{2.5_obs} of 150 µg m⁻³. 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 V_t are obtained. Supplement Fig. 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 to 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, M_t 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 (PM_{2.5_mod}) and daily mean observed PM_2.5 (PM_{2.5_obs}) as MB $={\sum }_{i=\mathrm{1}}^N({\mathrm{PM}}_{\mathrm{2.5}\mathrm{\_}\mathrm{mod}}-{\mathrm{PM}}_{\mathrm{2.5}\mathrm{\_}\mathrm{obs}})/N$, where N is total number of all samples; and (3) set M_t=150 µg m⁻³ + MB. Using this approach, M_t is 150 µg m⁻³ if the simulated PM_2.5 has no bias, and M_t is smaller (larger) than 150 µg m⁻³ 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 to 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⁻³, and therefore M_t is 121.9 µg m⁻³. As a result, in Baoding, the number of observed SWHDs is 210 d and that of the simulated SWHDs is 186 d 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 to 2017.

Figure 2 (a) Spatial distributions of simulated (CTRL, shades) and observed (CNEMC, dots) seasonal (DJF) mean surface-layer PM_2.5 concentrations (µg m⁻³) averaged over five winters (2013–2017). (b) Scatterplot of simulated versus observed DJF mean PM_2.5 concentration (µg m⁻³) averaged from 2013 to 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 coefficient between simulated and observed PM_2.5 concentrations. $\mathrm{NMB}\left(\text{normalized mean bias}\right)=({\sum }_{i=\mathrm{1}}^N\left(M_i-O_i\right)/{\sum }_{i=\mathrm{1}}^N\left(O_i\right))\times \mathrm{100}\mathit{\%}$, where O_i and M_i are the observed and simulated PM_2.5 concentrations, respectively. i refers to the ith site, and N is the total number of sites.

3 Simulated severe winter haze and model evaluation

3.1 Surface-layer PM_2.5 concentrations in DJF

Figure 2a presents the simulated seasonal (DJF) mean surface-layer PM_2.5 concentrations averaged from 2013 to 2017 over 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_2.5 concentrations were high over eastern China, covering important economic zones including BTH, the Sichuan Basin (SCB), and the Yangtze River Delta (YRD). The simulated highest concentrations occurred in southern BTH, with seasonal mean PM_2.5 in the range of 135–151 µg m⁻³. 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. Figure 2b shows 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 Simulated (purple solid line, from CTRL simulation) and observed (green solid line, from the US embassy and CNEMC) daily mean surface-layer PM_2.5 concentrations (µg m⁻³) 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⁻³ (red dashed line), the correlation coefficients (R), and NMB values between observed and simulated daily mean PM_2.5 concentrations.

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 US embassy site were used for comparison from 2009 to 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 to 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 low bias 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 19–26 December 2015, there were seven SWHDs in Beijing, and the highest PM_2.5 concentration reached 537 µg m⁻³ on 25 December. However, the model captured only three SWHDs during this period, and the simulated PM_2.5 concentration was 37 µg m⁻³ on this day. Such underpredictions of PM_2.5 concentrations during severe haze events were also reported by previous studies (L. Zhang et al., 2015; Zhang et al., 2014; Wang et al., 2014b) and were attributed to the representation deficiency of the aqueous/heterogeneous processes in the models (Huang et al., 2014; B. Zheng et al., 2015). For Shanghai and Chengdu, the model has high biases with NMBs of 18.6 % and 28.7 %, respectively.

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

3.2 Frequency and intensity of severe winter haze days

Figure 4a shows the simulated frequencies of SWHDs at 161 model grids averaged from 2013 to 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 d, which is more than twice the mean value of 8.7 d for all 161 grids. The highest frequencies of 31–39 d were located in southern BTH, indicating that more than one-third of the wintertime days experienced severe PM_2.5 pollution from 2013 to 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 spatial distributions are captured fairly well by the model. The scatterplot (Fig. 4c) of simulated versus observed frequencies of SWHDs for all model grids in Fig. 4a and b indicates a high correlation coefficient R of 0.98 and an NMB of −12.3 %. Over BTH, the model underestimated the frequency of SWHDs with an NMB of −8.3 %.

Figure 4d shows 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. Note that, for each grid, the simulated intensity was an adjusted value according to the model mean bias (MB) calculated in Sect. 2.4. For example, for the grid at Baoding where the simulated wintertime PM_2.5 has a MB of −28.1 µg m⁻³, the intensity of simulated SWHDs is adjusted as 249.6 µg m⁻³ (221.5 µg m⁻³ minus −28.1 µg m⁻³) over 2013–2017. The simulated mean PM_2.5 concentrations of SWHDs were high over the BTH and Shandong and Henan provinces. The average SWHD PM_2.5 concentration was 228 µg m⁻³ over BTH, which is much higher than the mean value of 195 µg m⁻³ 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 YRD, the Fen-Wei Plain, and northeast China. Accounting for all model grids in Fig. 4d and e, the scatterplot (Fig. 4f) shows an R of 0.58 between simulated and observed intensities of SWHDs, and the simulated results have an NMB of −4.1 %. Over BTH, simulated intensities of SWHDs have a small NMB of −1.7 %.

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). Figure 5 compares 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 d per decade during 1985–2002 and an increasing trend of 1.3 d per decade during 2003–2017; the model partly reproduces these trends with tendencies of −5.2 d per decade from 1985 to 2002 and 3.3 d per decade from 2003 to 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 capture the spatial distributions and historical changes in SWHD frequencies and intensities in China reasonably well.

Figure 5 Temporal evolutions of SWHD frequency (days per DJF) at the Beijing grid, which were derived from observed PM_2.5 (red, 2009–2016), simulated PM_2.5 (gray, 1985–2017), and observed atmospheric visibility (green, 1985–2017).

4 Simulated decadal changes in severe winter haze days

4.1 Decadal changes over China

Figure 6a and b show the spatial distributions of simulated SWHD frequencies and intensities averaged from 1985 to 2017 over China. Simulated SWHDs occurred most frequently in southern BTH, and the intensity showed the highest values in the regions of BTH and Shandong and Henan provinces. For all model grids shown in Fig. 6a, the averaged frequency of SWHDs was 5.9 d, and the mean intensity was 188 µg m⁻³. Figure 6c and d show the simulated temporal changes in SWHD frequency and intensity averaged over all model grids in Fig. 6a from 1985 to 2017. Both the frequency and intensity of SWHDs exhibited increasing trends over the past 3 decades. The simulated mean frequency of SWHDs increased from 1.8 d in 1985 to 4.6 d in 2017, with a large growth rate of 2.6 d per decade. The mean SWHD intensity increased from 178 µg m⁻³ in 1985 to 195 µg m⁻³ in 2017, with a linear trend of 7.1 µg m⁻³ per decade.

Figure 6 Spatial distributions of SWHD (a) frequencies (days) and (b) intensities (µg m⁻³) at 161 model grids averaged from 1985 to 2017. Temporal evolutions of SWHD (c) frequency (days) and (d) intensity (µg m⁻³) in China from 1985 to 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.

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 d per winter from 1985 to 2000, but the value increased rapidly from 2001 to 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 d per decade. Interestingly, negative trends were found over BTH during this period, with a maximum trend of −6.9 d per decade. From 2001 to 2011, the SWHD frequencies increased rapidly on a national scale. The largest increasing trends of more than 20 d per decade 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 Province, with decreasing trends of greater than −45 d per decade.

Figure 7 Linear trends of simulated frequencies (d per decade) 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.

4.2 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 includes 18 grids, a regional SWHD is identified when six 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 Simulated temporal variations in frequency (days, black line), intensity (µg m⁻³, red line), and concentrations of PM_2.5 components (µg m⁻³, bars): ammonia (blue), nitrate (yellow), sulfate (red), black carbon (purple), and organic aerosols (green) of regional SWHDs in BTH from 1985 to 2017. Also shown are linear trends (dashed lines) for frequency and intensity, which are statistically significant above the 95 % confidence level.

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 to 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 d in 1985 to 29 d in 1992, followed by a sudden decrease to 10 d in 2001, a rapid growth to a peak of 47 d in 2012, and a steep decrease thereafter to 15 d in 2017, which coincides with the temporal variations presented in previous studies for haze days (Chen and Wang, 2015), aerosol extinction coefficient (J. Li et al., 2016), and PM_2.5 concentration (Yang et al., 2018) in the North China Plain (NCP). From 1985 to 2017, SWHDs in BTH exhibited large increasing trends of 4.5 d per decade in frequency and 13.5 µg m⁻³ per decade in intensity. In comparison, the seasonal mean PM_2.5 concentrations showed a smaller trend of 10.1 µg m⁻³ per decade (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.6 µg m⁻³ per decade, which was followed by ammonium (4.5 µg m⁻³ per decade). Sulfate and BC made little contribution to the increasing trend, while OC showed the opposite trend of −6.8 µg m⁻³ per decade, which coincided with the emission changes during this period (Fig. 1). Note that these simulated trends of PM_2.5 components are uncertain, which can be influenced by the uncertainties in emission inventories of aerosols and aerosol precursors (Crippa et al., 2018; M. Li et al., 2017; Zheng et al., 2018) as well as by the chemistry scheme in the model (Chen et al., 2019). The simulations in this study did not include secondary organic aerosols; therefore, the concentrations of organic aerosols might have been underestimated. Previous studies also showed that the GEOS-Chem model tends to overestimate nitrate and underestimate sulfate in winter (Wang et al., 2013).

5 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 to 2017. PA has been widely used in previous studies to identify the relative importance of atmospheric processes in specific pollution episodes (Gonçalves 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 the surface to the 11th 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:

$$\begin{array}{}\text{(1)}& \mathit{\%}{\mathrm{PC}}_i={\displaystyle \frac{{\mathrm{PC}}_{\mathrm{SWHD}\mathit{\_}i}-{\mathrm{PC}}_{\mathrm{SM}\mathit{\_}i}}{{\sum }_i^n\mathrm{abs}({\mathrm{PC}}_{\mathrm{SWHD}\mathit{\_}i}-{\mathrm{PC}}_{\mathrm{SM}\mathit{\_}i})}}\cdot \mathrm{100},\end{array}$$

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 to 2017, and PC_{SWHD_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 %.

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 to 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. The %PC_i is calculated following Eq. (1). The pie chart presents the relative contributions made by three wind components: north–south (NS), east–west (EW), and up–down (UD) for the transport difference (PC__{SWHD_tran}−PC__{SM_tran}). Note that carbon here refers to the carbonaceous aerosols, which is the sum of black carbon and organic aerosols.

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 winter from 1985 to 2017 (stripped bars) and over the SWHDs during this period (solid bars). Averaged over all days during winter from 1985 to 2017, the PM_2.5 in BTH was mainly contributed by chemical reactions (10.1 Gg d⁻¹) and local emissions (5.5 Gg d⁻¹). Net transport had an effect of −12.1 Gg d⁻¹, 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, −1.8, and −0.9 Gg d⁻¹, respectively. For all days during winter from 1985 to 2017, the sum of all processes had a very small value of 0.004 Gg d⁻¹. During SWHDs in the selected time period, the outflow of particles decreased greatly, from 12.1 Gg d⁻¹ in mean state to 4.8 Gg d⁻¹, which held more particles in BTH. More secondary aerosols were produced through chemical reactions, with an increase from a mean state of 10.0 to 12.0 Gg d⁻¹. The suppressed PBL height further restricted the PM_2.5 diffusion with reduced aerosol outflow from 0.8 to 0.4 Gg d⁻¹. 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 d⁻¹, leading to the accumulation of PM_2.5 in BTH. The relative contributions of transport, chemistry, cloud 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. It should be noted that although process analysis is a helpful tool to quantify the contribution of each factor, it assumes a linear contribution from each factor. This may not always be true because of the covariation in meteorological parameters.

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.

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 to 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 (Q. Li et al., 2016; 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.

Figure 10 (a) Time series of frequencies (days) of regional SWHD in BTH from three simulations (CTRL, EMIS, MET) for 1985 to 2017. (b–e) Time series of linear trends calculated over different periods for simulated frequencies of the (b) CTRL, (c) EMIS, and (d) MET simulations. 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. The filled color in each square shows the calculated trend value, and those values marked with black borders are significant above the 95 % confidence level.

6 Roles of anthropogenic emissions and meteorological parameters

Figure 10a shows the time series of regional SWHD frequencies in BTH from 1985 to 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 to 2017, with a linear trend of 4.5 d per decade (as described in Sect. 4.2). During 1985–2017, the simulated SWHD frequency exhibited an increasing trend of 5.9 d per decade 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 Fig. 10b–d 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 d per decade) was the largest in the CTRL simulation, the trends in EMIS and MET were −3.3 and −20.1 d per decade, respectively, indicating that the variations in meteorological parameters dominated the decreasing trend over this period. The trend of −3.3 d per decade from 1992 to 2001 in EMIS was consistent with the emission reductions over this period (Fig. 1). During 2003–2012 when the increasing frequency trend (23.6 d per decade) was the largest in the CTRL simulation, the EMIS and MET trends were 27.3 and 12.5 d per decade, 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 to 2001, those from 1984 to 1992 and 2002 to 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 11 (a) Time series of intensities (µg m⁻³) of regional SWHD in BTH from the three simulations (CTRL, EMIS, MET) for 1985–2017. (b–e) Time series of linear trends calculated over different periods for simulated intensities of the (b) CTRL, (c) EMIS, and (d) MET simulations. 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. The filled color in each square shows the calculated trend value, and those marked with black borders are significant above the 95 % confidence level.

Figure 11a shows the time series of the regional SWHD intensities in BTH from 1985 to 2017 based on the CTRL, EMIS, and MET simulations. From 1985 to 2015, the simulated SWHD intensity exhibited increasing trends of 13.5, 5.2, and 8.3 µg m⁻³ per decade 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 to 2001 (Fig. 11b), due to the reasons mentioned above.

Figure 12 Interannual variations in SWHD (a) frequencies (days) and (b) intensities (µg m⁻³) in BTH through a 9-year weighted moving average method. The nine-point high-pass Lanczos filter was used as the weighting coefficient, 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 EMIS simulations (red R).

Additionally, we examine the interannual variations in SWHDs in BTH. The nine-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 to 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 in 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.93 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:

$$\begin{array}{}\text{(2)}& \mathit{\%}{\mathrm{IC}}_i={\displaystyle \frac{\mathrm{abs}(\mathrm{MET}\mathrm{\_}{\mathrm{IA}}_i-\mathrm{MET}\mathrm{\_}{\mathrm{IA}}_{i-\mathrm{1}})}{\mathrm{CTL}\mathrm{\_}{\mathrm{F}}_{i-\mathrm{1}}}}\cdot \mathrm{100},\end{array}$$

where MET_IA_i and MET_IA_i−1 are the interannual anomalies in the MET simulation in years i and i−1, and CTL_F_i−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 to 2013, and 4.6 % 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 ±4.6 % in PM_2.5 concentrations for the SWHDs in BTH interannually. This highlights the importance of variations in meteorological parameters for policy makers when actions are taken to improve air quality in the short term.

7 Conclusions

We investigated the changes in SWHDs (days with observed PM_2.5 concentrations larger than 150 µg m⁻³) from 1985 to 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_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 to 2016, 0.70 % and 18.6 % in Shanghai from 2013 to 2017, and 0.58 % and 28.7 % in Chengdu from 2013 to 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 to 2017, the simulated frequency had an R of 0.98 and NMB of −12.3 %, and the simulated intensity had an R of 0.58 and NMB of −4.1 %.

The Beijing–Tianjin–Hebei region (BTH) had the highest SWHD frequency and intensity based on the observed PM_2.5 concentrations from 2013 to 2017, with a frequency of 21.2 d yr⁻¹ and a mean intensity of 231.6 µg m⁻³. Historically, from 1985 to 2017, with changes in emissions and meteorology, the simulated SWHD frequency in BTH increased with a linear trend of 4.5 d per decade, and the simulated intensity increased by 13.5 µg m⁻³ per decade, which was larger than the increase in the seasonal mean PM_2.5 concentration (10.1 µg m⁻³ per decade). 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 to 2017 indicate that transport (outflow of particles) had the largest change from −12.1 Gg d⁻¹ in the mean state to −4.8 Gg d⁻¹ in SWHDs, which had the largest contribution to SWHD formation. Transport, chemistry, cloud processes, dry deposition, and PBL mixing had relative contributions of 65.3 %, 17.6 %, −7.5 %, −6.4 %, and 3.2 %, respectively, to SWHD formation.

From 1985 to 2017, the simulated frequency of SWHDs in BTH experienced several stages, including an increase from 14 d in 1985 to 29 d in 1992, a sudden drop to 10 d in 2001, a rapid growth to the peak of 47 d in 2012, and a steep decrease thereafter to 15 d in 2017. Sensitivity studies showed that the decrease in frequency from 1992 to 2001 was mainly caused by changes in meteorological parameters. From 2003 to 2012, when the SWHD frequency increased sharply, the simulated frequency trends were 23.6, 27.3, and 12.5 d per decade 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 to 2017; variations in emissions alone and meteorology alone both enhanced the intensity growth with trends of 5.2 and 8.3 µg m⁻³ per decade, 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.