Spatiotemporal distributions of ambient SO2 across China based on satellite retrievals and ground observations: Substantial decrease in human exposure during 2013–2016

Environmental Research 179 (2019) 108795

Contents lists available at ScienceDirect

Environmental Research journal homepage: www.elsevier.com/locate/envres

Spatiotemporal distributions of ambient SO2 across China based on satellite retrievals and ground observations: Substantial decrease in human exposure during 2013–2016

T

Hanyue Zhanga, Baofeng Dia,b, Dongren Liua, Jierui Lia, Yu Zhana,c,d,∗ a

Department of Environmental Science and Engineering, Sichuan University, Chengdu, Sichuan, 610065, China Institute for Disaster Management and Reconstruction, Sichuan University, Chengdu, Sichuan, 610200, China c Sino-German Centre for Water and Health Research, Sichuan University, Chengdu, Sichuan, 610065, China d Medical Big Data Center, Sichuan University, Chengdu, Sichuan, 610041, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: SO2 China Remote sensing Spatiotemporal distribution Human exposure

Multiyear spatiotemporal distributions of daily ambient sulfur dioxide (SO2) are essential for evaluating management effectiveness and assessing human health risk. In this study, we estimate the daily SO2 levels across China on 0.1o grid from 2013 to 2016 by assimilating satellite- and ground-based SO2 observations using the random-forest spatiotemporal kriging (RF-STK) model. The cross-validation R2 is 0.64 and 0.81 for predicting the daily and multiyear averages, respectively. The multiyear population-weighted average of SO2 for China is 28.1 ± 14.0 μg/m3, and the severest SO2 pollution occurs in the northern China (45.1 ± 14.7 μg/m3). The SO2 concentration shows a strong seasonality, i.e., highest in winter (41.6 ± 26.4 μg/m3) and lowest in summer (19.6 ± 8.3 μg/m3). During 2013–2016, the annual SO2 decreases from 34.4 ± 18.2 to 22.7 ± 11.1 μg/m3, and the population% exposed for more than 100 nonattainment days (SO2 > 20 μg/m3) drops from 86% to 48%. While the seasonality of SO2 is mainly determined by the meteorological variation, the substantial decrease attributes to the reduced emissions such as from coal consumption. The effectiveness of SO2 emission reduction varies widely in different prefectures of China. In Shandong province, the SO2 concentration decreases by −45% while the coal consumption increases by 9%. In Shanxi province, the SO2 concentration decreases by −15% while the coal consumption decreases by −3%. The contrasting effectiveness between these two provinces is associated with the much fewer waste gas disposal facilities in Shanxi than Shandong. Stricter regulation is required to further lower the SO2 concentration in order to protect the public health, especially in the northern China.

1. Introduction Sulfur dioxide (SO2), as one of the primary air pollutants in the world, poses risks to human health and leads to environment acidification. Atmospheric SO2 originates from natural sources such as volcanic emission and forest fire, as well as anthropogenic emissions such as combustion of fossil fuel (Lu et al., 2013). SO2 can be oxidized to sulfate aerosols, which are among the major components of fine particulate matter (Philip et al., 2014). SO2 endangers public health, notably damaging respiratory system (Iwasawa et al., 2009; WHO, 2006). Sulfate aerosols contribute to acid deposition, degrade air quality, affect climate, and harm human health (Charlson et al., 1992; Dentener et al., 2006; Lee et al., 2015; van Donkelaar et al., 2008). Considering the hazard effects of SO2 and its derivatives, the World

∗

Health Organization (WHO) set the 24-h air quality guideline (AQG) of 20 μg/m3 for SO2 (WHO, 2006). While SO2 has been properly disposed in the developed countries/regions, SO2 pollution is still severe in many developing countries, such as China and India (Krotkov et al., 2016; Li et al., 2017a). China has taken various measures to prevent and control SO2 emissions, among which the installation of flue gas desulphurization (FGD) device plays an important role (Lu et al., 2010). Recently, Chinese government issued the “Three-year Plan on Defending the Blue Sky” with stricter standards to protect the atmospheric environment (SCC, 2018). Full-coverage and continuous spatiotemporal distributions of daily ambient SO2 levels are essential for evaluating/refining the SO2 management practices and assessing human exposure. Satellite retrievals provide large-scale and regular monitoring of atmospheric SO2 column densities. The instruments of the Tropospheric

Corresponding author. Department of Environmental Science and Engineering, Sichuan University, Chengdu, Sichuan, 610065, China. E-mail address: [email protected] (Y. Zhan).

https://doi.org/10.1016/j.envres.2019.108795 Received 20 May 2019; Received in revised form 31 August 2019; Accepted 3 October 2019 Available online 03 October 2019 0013-9351/ © 2019 Elsevier Inc. All rights reserved.

Environmental Research 179 (2019) 108795

H. Zhang, et al.

Fig. 1. Spatial distribution of SO2 monitoring sites in China during 2013–2016.

Emission Spectrometer (TES), the Atmospheric Infrared Sounder (AIRS), and the Infrared Atmospheric Sounding Interferometer (IASI) relying on infrared (IR) are less sensitive for detecting the SO2 in the planetary boundary layer (PBL), as the temperature differences between ground surface and atmosphere are narrower (Carn et al., 2005; Clarisse et al., 2008; Clerbaux et al., 2008; Fioletov et al., 2016). Instead, sensors in the ultraviolet (UV) spectral range are commonly employed to monitor atmospheric SO2, including the Total Ozone Mapping Spectrometer (TOMS), the Global Ozone Monitoring Experiment (GOME), the Scanning Imaging Absorption spectrometer for Atmospheric Chartography (SCIAMACHY), the GOME-2, the Ozone Monitoring Instrument (OMI), and the Ozone Mapping and Profiler Suite (OMPS), which are summarized in Table S1 (Eisinger and Burrows, 1998; Fioletov et al., 2011; Krueger, 1983; Xu et al., 2010; Yan et al., 2014; Zhang et al., 2016). Among these apparatuses, OMI possesses a long operating period (2004-present), fine spatial (13 × 24 km2) and temporal (daily) resolutions, and high sensitivity (Fioletov et al., 2013). In the official release of OMI, the principal component analysis (PCA) algorithm has recently replaced the band residual difference (BRD) method to obtain more accurate SO2 inversions in the PBL (Li et al., 2013). The OMI SO2 retrievals (OMI-SO2) have been applied to estimate volcanic SO2 fluxes, detect SO2 sources, evaluate SO2 emissions, and analyze SO2 spatiotemporal variations (Fioletov et al., 2011, 2013; Theys et al., 2013; Yan et al., 2014). However, high uncertainties are associated with the satellite retrieved SO2 spatiotemporal patterns. Climatic and topographic conditions would influence satellite inversions. Especially in cloudy weather, radiations are normally prevented from reaching the ground. Satellite instruments and corresponding algorithms introduce uncertainties, although they have been continuously refined. As for SO2, the radiation absorptivity is seriously influenced by ozone (O3) on account of their similar optical properties (Fioletov et al., 2013). The OMI sensor are capable of capturing strong signals within detection limit at high SO2

emission-level regions, but are less effective to quantify low SO2 concentrations (Fioletov et al., 2015, 2016). Most studies thus focus on SO2 point sources such as volcanoes and power plants (Fioletov et al., 2011, 2013, 2015; Theys et al., 2013). Given the limitations of satellite products, OMI-SO2 should be calibrated against ground observations to obtain the spatiotemporal distributions of SO2 with higher accuracy. Various statistical approaches are developed to refine the estimation of atmospheric pollutant levels derived from satellite retrievals. The generalized additive model (GAM), the geographically weighted regression (GWR), and the land use regression (LUR) are used in deriving PM2.5 and NO2 concentrations based on satellite aerosol optical depth (AOD) and NO2 column densities (Liu et al., 2009; Song et al., 2014; Vienneau et al., 2013). In addition, machine learning algorithms gradually become popular in various fields for solving predictive problems. Machine learning models can handle complex relationships (e.g., nonlinearity and interactions) between independent and dependent variables, showing superior performance than traditional statistical models (Hastie et al., 2009). Based on machine learning and geostatistics, we have developed a random-forest spatiotemporal kriging (RF-STK) model to reconstruct the spatiotemporal distributions of daily ambient NO2 concentrations (Zhan et al., 2018b). In this study, we refine the original RF-STK model to estimate the spatiotemporal distributions of daily ambient SO2 concentrations across China on 0.1° grid (98,341 cells) for the period of 2013–2016. The spatial and temporal configurations are consistent with our previous study, which belongs to a series of air pollution researches (Zhan et al., 2017, 2018a, 2018b). The influential factors on the SO2 levels are characterized by the statistical measures of variable importance. Based on the prediction results, we analyze the spatial patterns and temporal trends of SO2, as well as assess the human exposure levels according to the WHO AQG. This study delineates the multiyear spatiotemporal distributions of daily SO2 across China during 2013–2016 in order to facilitate the air quality management and human health risk 2

Environmental Research 179 (2019) 108795

H. Zhang, et al.

the area-weighted average method. An additional set of predictor variables are generated by operating the spatial convolutions with Gaussian kernels on the emission inventories, land uses, NDVI, road density, population density, and elevation, which reflect the pollutant dispersion or neighboring effect (Goodfellow et al., 2016). Moreover, day of year (DOY) and year are introduced as surrogate variables to indicate intra- and inter-annual variabilities. The specific information of all the independent variables is listed in Table S2.

assessments. 2. Materials and methods 2.1. SO2 ground monitoring The daily average ambient SO2 concentrations during 2013–2016 are summarized from the hourly concentrations monitored using the ultraviolet fluorescence method at the air quality monitoring network in China (EPAROC, 2015; EPDHK, 2015; MEPC, 2015; Zhao et al., 2016). All the state-managed sites (1657 in total) are used in this study, with 744, 1022, 1612, and 1604 sites for 2013, 2014, 2015, and 2016, respectively (Fig. 1). Most of the monitoring sites are distributed in urban areas and are relatively sparse in remote areas such as Tibet and Xinjiang Uygur Autonomous Regions. After data screening, a total of 1.67 million records of daily average SO2 concentrations are obtained for training the RF-STK model.

2.4. RF-STK model The RF-STK model is comprised of a random forest (RF) submodel and a spatiotemporal kriging (STK) submodel. The RF submodel combines the SO2 ground observation, OMI-SO2, and environment factors to estimate the ambient SO2 levels, which is implemented in python with the package scikit-learn (Pedregosa et al., 2011). The SO2 ground observation is the dependent or response variable, while the OMI-SO2 and environment factors are used as the independent or predictor variables. The training data contain the information of both dependent and independent variables at all the sites, while the prediction data include information of only independent variables for all the grid cells. The RF submodel is built with 500 regression trees, each of which is trained with the bootstrapped samples from the training data (Breiman, 2001). The SO2 levels are estimated as the average predictions of all the trees based on the prediction data. The variable importance (IM) and partial dependence plots are employed to interpret the RF submodel. IM reflects the variable contributions to SO2, i.e., the influences of predictor variables on the SO2 levels. The predictor variables are divided into ten categories for simpler interpretations (Table S2). The IM of each category is the sum of the IMs of all the individual predictor variables in that category. The partial dependence plots delineate the relationship between each predictor variable and SO2 levels while keeping the effects of all the other predictor variables at their averages. The STK submodel, implemented with the package gstat, is employed to interpolate the out-of-bag (OOB) residuals of RF, utilizing the spatiotemporal information remaining in residuals (Graler et al., 2016). Note that the OOB residuals are screened by the three-sigma (3σ) rule to remove outliers. The final SO2 concentrations are estimated as the sum of the predictions made by the two submodels.

2.2. OMI SO2 retrievals The SO2 vertical column densities in the PBL are extracted from the level-3 OMI SO2 products (OMSO2e) released by the Goddard Earth Sciences Data and Information Services Center (GES DISC) of the National Aeronautics and Space Administration (NASA) (GES DISC, 2015). OMI is aboard the polar sun-synchronous Aura satellite, which was launched by the NASA's Earth Observing System (EOS) on July 15, 2004. OMI provides daily and global information of aerosols and trace gases with transmit local time around 13:45 p.m. The column density of atmospheric SO2 is derived based on the OMI bands at the wavelength of 310.8–314.4 nm. The OMSO2e, with a regular latitude-longitude grid at a resolution of 0.25o, contains the best pixel data for SO2 in PBL, where SO2 is most closely related to human activities. The OMI-SO2 are presented in the Dobson unit (DU; 1 DU = 2.69 × 1016 molecules cm−2), which is defined as the thickness of pure gas layer under the condition of standard temperature and pressure (OMI Team, 2012). Note that the Tropospheric Monitoring Instrument (TROPOMI) aboard on the Sentinel-5 Precursor (launched at October 13, 2017) retrieves SO2 at higher resolution (7 × 7 km2), which are valuable for future work deriving ambient SO2 since 2018 (Theys et al., 2017). The OMISO2 are resampled to the predefined 0.1° grid using the area-weighted average method, which are then processed with the temporal convolution with Gaussian kernels for filtering noises and filling data gaps (Goodfellow et al., 2016). The OMI-SO2 coverage rate of each grid cell is calculated by dividing the number of days with OMI-SO2 data by the total number of days over a specified period (e.g., a season or a year), and the national OMI-SO2 coverage rate is the average over all the grid cells.

2.5. Cross-validation The ten-fold cross-validation is conducted to assess the prediction performance of the RF-STK model. With stratification by monitoring sites, the whole training data are approximately equally divided into ten parts at random. In each round, nine parts are for training a new testing RF-STK model, and then the remaining part is for validating the predictions made by this testing model. After ten rounds, the complete predictions corresponding to the observations in training data are obtained. The daily predictions and observations are summarized to different levels (i.e., monthly, seasonal, annual, and multiyear levels) for calculating the predictive performance across multiple temporal scales. The commonly used statistical indices, including coefficient of determination (R2), root mean square error (RMSE), relative prediction error (RPE), and slope, are employed to evaluate the prediction performance of the RF-STK model.

2.3. Environmental factors Environmental factors, including weather conditions (CMA, 2015), Planetary Boundary Layer Height (PBLH) (GMAO, 2015), emission inventories (Li et al., 2017b), land uses (Jun et al., 2014), Normalized Difference Vegetation Index (NDVI) (Didan et al., 2015), road density (OpenStreetMap contributors, 2015), population density (CIESIN, 2016), and elevation (Jarvis et al., 2008), are compiled from miscellaneous sources. These covariates are selected based on the empirical knowledge of their influences on ambient SO2 and the data availability (Li et al., 2019; Zhan et al., 2018b). The daily weather conditions are interpolated by utilizing the co-kriging interpolation with elevation (Deutsch and Journel, 1998). The emission inventories for the years of 2013 and 2015 are linearly interpolated from the available data of 2012, 2014, and 2016. The annual NDVI are aggregated from the 8-day NDVI through the maximum value compositing approach (Holben, 1986). The PBLH, emission inventories, land uses, NDVI, road density, population density, and elevation are assigned to the grid by applying

2.6. Spatiotemporal distribution analysis and human exposure assessment The RF-STK model trained by the complete dataset of training data predicts daily SO2 concentrations on 0.1o grid in China from 2013 to 2016. The prediction results are presented as the population-weighted average SO2 concentrations at different spatial and temporal levels using Eq. (S1). The four-year average population-weighted SO2 concentrations for the whole nation and the 34 prefectures are calculated to indicate distribution characters. The monthly, seasonal, and annual 3

Environmental Research 179 (2019) 108795

H. Zhang, et al.

The OMI-SO2 generally manifested a normal distribution, with the mean ± standard deviation of 0.03 ± 0.40 DU (Table S4). The average OMI-SO2 were the highest in Shandong province (0.34 ± 0.84 DU) and the lowest in Fujian province (−0.01 ± 3.7 DU). The OMISO2 was higher in fall (0.05 ± 0.41 DU) than in spring (0.02 ± 0.38 DU), which was inconsistent with the seasonality of ground observations. The OMI-SO2 showed a decreasing trend with annual OMI-SO2 of 0.05 ± 0.44 and 0.02 ± 0.36 DU for the year of 2013 and 2016, respectively. The correlation between the OMI-SO2 and ground concentrations was extremely low at daily level (R = 0.17; Table S5), as their relationship was affected by various atmospheric conditions such as vertical profiles. The atmospheric conditions became less influential when evaluating the correlation between the OMI-SO2 and ground concentrations at annual level (R = 0.67). In addition, the four-year nationwide coverage rate of the OMI-SO2 was 29%, with the highest coverage in summer (35%) and the lowest in winter (30%; Fig. S2). After the temporal convolution, the data gaps in the OMI-SO2 were filled, and the correlation between the daily OMI-SO2 and ground observations was strengthened (R = 0.45).

Table 1 Site-based ten-fold cross-validation results of the RF-STK model. Metrica 2

R RMSE RPE Slope MFB MFE MNB MNE N

Daily

Monthly

Seasonal

Annual

Multiyear

0.64 19.5 73% 0.67 0.12 0.44 0.46 0.70 1,672,984

0.74 14.0 52% 0.77 0.09 0.33 0.26 0.46 56,835

0.78 12.0 44% 0.80 0.08 0.29 0.21 0.38 19,364

0.79 9.4 34% 0.80 0.06 0.24 0.13 0.29 4982

0.81 7.9 30% 0.79 0.05 0.22 0.11 0.26 1657

a R2: coefficient of determination; RMSE: root mean square error (μg/m3); RPE: relative prediction error; MFB: mean fractional bias; MFE: mean fractional error; MNB: mean normalized bias; MNE: mean normalized error; and N: number of samples.

average population-weighted SO2 levels are computed to investigate the intra- and inter-annual trends. The seasonal trend decomposition using loess (STL) algorithm is used to obtain long-term deseasonalized trends (Cleveland, 1990). In addition, the predicted SO2 distributions are correlated with the coal consumption and the number of FGD facilities, which are retrieved from the China Energy Statistical Yearbook and the China Statistical Yearbook on Environment, respectively (NBSC, 2017a; NBSC, 2017b). The annual human exposure is evaluated through associating daily SO2 concentrations with the WHO 24-h Interim Target-1 (IT-1) of 125 μg/m3, Interim Target-2 (IT-2) of 50 μg/ m3, and AQG of 20 μg/m3 (WHO, 2006). We evaluate the human exposures under different levels subject to the constant cumulative human exposure (i.e., exposure intensity × duration). The areas with dense population (> 100 people/km2), which is an empirical reference value distinguishing between rural and urban areas (Zhan et al., 2018b), are annotated with the number of nonattainment days (SO2 > 20 μg/m3) in order to show the spatial distributions of the population at high SO2 exposure risk.

3.2. Prediction performance The RF-STK model shows a good prediction performance in the cross-validation, with R2 = 0.64, RMSE = 19.5 μg/m3, RPE = 73%, and slope = 0.67 in predicting the daily average concentrations (Table 1). In order to obtain a parsimonious RF-STK model, we performed the backward variable selection by removing the least important variable one at a time (Zhan et al., 2018a, 2018b). Whereas the predictions made by the reduced RF-STK model were dubious for the remote areas such as the Tibetan Plateau due to erroneous removals of important predictors. We therefore choose the RF-STK model with all the predictor variables for the present research, as accurate predictions are more important than parsimonious models. Better prediction performance is achieved at monthly (R2 = 0.74), seasonal (R2 = 0.78), annual (R2 = 0.79), and multiyear (R2 = 0.81) levels. The mean normalized error (MNE) shows no significant difference among the years and the seasons (except for summer), indicating the robust performance of RF-STK over the time (Table S6). Relatively higher MNE for summer is mainly due to the lower signal-to-noise ratio of OMI-SO2 under the situation of low SO2 concentrations. The RF-STK model presents better performance than other approaches to reconstructing the SO2 spatiotemporal patterns. Compared to the low correlation between the daily OMI-SO2 and the ground SO2 observations, the RF-STK model largely improves the accuracy of mapping daily SO2. The predictions by the RF-STK model are featured by the detailed characterization of spatially heterogeneous distributions when compared to the kriging interpolations (Figs. S3 and S4). The RFSTK model moreover shows superior performance to random forest, kriging interpolation, and regression kriging (Li et al., 2019; Zhan et al., 2018b). However, the RF-STK model is not applicable for long-term back-extrapolation (i.e., making predictions for the period earlier than the observations), as spatially and temporally neighbored observations are necessary for running the STK submodel. Considering the importance of high-resolution data to local-scale characterization and acute human exposure assessment, we plan to refine the current simulation resolution (i.e., 0.1o and 1 day) to 1 km and an hour in the future. The RF-STK model can also be employed to derive the spatiotemporal distributions of other air pollutants (e.g., PM2.5 and NO2), which are regularly monitored at ground-based sites and through spaceborne remote sensing.

3. Results and discussion 3.1. Basic statistics of ground- and space-based SO2 observations The daily SO2 concentrations observed at the monitoring sites ranged from 0 to 1553.2 μg/m3, with the mean ± standard deviation of 26.6 ± 32.1 μg/m3 and the median of 16.7 μg/m3, following a lognormal distribution (Table S3). The average diurnal cycle of the observations showed that the ground concentrations peaked around 10am (Fig. S1). The SO2 concentrations at 14pm when the OMI sensor transmitted around were highly correlated with the daily average concentrations (Spearman correlation coefficient, R = 0.86). The multiyear average SO2 levels were the highest in Shanxi province (62.0 ± 65.7 μg/m3) and the lowest in Hainan province (5.2 ± 3.5 μg/m3). The highest SO2 concentration was observed in winter (42.2 ± 47.5 μg/m3), followed by spring (24.5 ± 23.9 μg/m3), fall (24.1 ± 27.6 μg/m3), and summer (16.1 ± 15.1 μg/m3). The annual average SO2 concentrations of China decreased from 33.8 ± 39.8 to 21.8 ± 26.9 μg/m3 during 2013–2016. Ground observations are essential for improving satellite retrievals and supporting air quality management. Short-term (e.g., daily) observations help improve the accuracy of satellite-based estimation of air pollution distributions through data assimilation (van Donkelaar et al., 2012). Long-term (e.g., multiyear) observations provide reference values for setting air quality standards suitable for the status quo. In addition, the observed SO2 distributions are affected by the environmental factors both spatially and temporally, providing essential information for analyzing the influences of meteorological conditions and anthropogenic emissions on the spatiotemporal variations of SO2 concentrations (Luvsan et al., 2012).

3.3. Variable contribution Among all the individual predictor variables in the RF-STK model, the OMI-SO2 is the most important predictor variable (IM = 14.0%; Fig. 2 and Table S7). The partial dependence plot shows that the SO2 4

Environmental Research 179 (2019) 108795

H. Zhang, et al.

Fig. 2. Importance (IM; %) of predictor variable categories in the RF-STK model. All the predictor variables are categorized into ten groups (Table S2). The IM of each individual variable is detailed in Table S7. Table 2 Population-weighted average of SO2 concentrations and temporal trends in all the prefectures/nation during 2013–2016. Prefecturea

Concentrationb

Trendc

%d

Prefecture

Concentration

Trend

%

Shanxi Hebei Shandong Ningxia Tianjin Liaoning Henan Inner Mongolia Shaanxi Jiangsu Qinghai Guizhou Anhui Gansu Hubei Jilin Hunan

56.7 50.8 46.5 41.3 38.2 37.6 37.5 29.7 28.7 28.0 27.6 27.3 27.1 26.9 25.8 24.5 23.7

−2.9 −8.8 −9.2 −3.2 −9.7 −3.8 −5.1 −3.8 −4.3 −4.0 −2.3 −4.6 −2.7 −2.1 −5.0 −1.7 −3.0

−14.5 −40.3 −45.1 −21.3 −53.8 −26.5 −33.5 −31.9 −36.0 −34.7 −21.4 −39.5 −25.1 −20.9 −44.6 −18.8 −31.2

Heilongjiang Jiangxi Chongqing Sichuan Beijing Shanghai Xinjiang Yunnan Guangxi Zhejiang Guangdong Tibet Fujian Hong Kong Taiwan Hainan Nation

23.6 ± 9.1 23.5 ± 6.0 23.1 ± 4.3 22.3 ± 4.6 21.1 ± 4.9 20.9 ± 2.8 20.9 ± 6.4 20.4 ± 5.9 19.5 ± 5.0 19.5 ± 6.9 16.9 ± 5.1 16.6 ± 3.8 13.5 ± 3.3 12.2 ± 1.4 9.1 ± 2.7 8.8 ± 3.1 28.1 ± 14.0

−1.5 −2.4 −5.3 −4.5 −5.4 −3.1 −2.0 −2.1 −1.7 −2.6 −2.1 −0.1 −1.5 −1.0 −0.5 −0.5 −3.9

−17.7 −26.8 −48.7 −44.1 −53.3 −36.3 −23.5 −26.2 −22.3 −32.0 −29.9 −1.9 −26.5 −20.6 −14.3 −14.9 −34.1

a b c d

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

13.2 12.8 14.8 18.6 4.4 14.7 9.6 12.1 10.3 5.4 7.3 5.7 5.8 7.3 5.0 7.3 4.6

Macao is not included in the Table. Multiyear average SO2 concentrations (mean ± standard deviation; μg/m3). Annual absolute deseasonalized trends (μg/m3). All trends are significant: P < 0.01. Four-year relative trends (%).

emissions of SO2. In warm season, SO2 dispersion is favored by enhanced atmospheric turbulence and convection, and high relative humidity promotes the conversion of SO2 to SO42−, resulting in lower SO2 concentrations (Cheng et al., 2016; Yang et al., 2017). In cold season, the increased coal consumption for heating leads to intensified SO2 emission, while low PBLH causes condensed pollutant density near ground surface (Wang et al., 2014; Zhan et al., 2018b). With negative radiance forcing, sulfate aerosols formed from SO2 in turn affect the weather conditions in different aspects, including reduction of solar radiance, decrease of atmospheric temperature, and shift of monsoon rain band (Charlson et al., 1992; Dey and Tripathi, 2008; Xu, 2001). In addition, the emissions of atmospheric pollutants are highly correlated with each other and collectively indicate the spatiotemporal variations of SO2 emissions (Fig. S7). The emission inventory of SO2 (IM = 4.9%) is enriched with those of organic carbon (OC; IM = 6.1%) and black carbon (BC; IM = 2.3%) for estimating the SO2 concentrations. This suggests that the actual SO2 emissions from residential activities, which are main emission sources of OC and BC, might be underestimated in

predictions are nonlinearly associated with the OMI-SO2 (Fig. S5). The predicted SO2 concentration increases approximately proportionally with the elevation of OMI-SO2 at higher level. The predictions are nevertheless insensitive to the OMI-SO2 variations when the SO2 concentrations are low. In addition, the OMI-SO2 shows higher correlation with the SO2 observations at the monitoring sites with severer pollution (Fig. S6). Both the partial dependence plot and the correlation trend demonstrate that the OMI-SO2 is capable of identifying highly polluted areas but is inadequate to distinguish variations in relatively clean areas which account for the major parts of China. It is thus crucial to include the covariates in the model to obtain the spatiotemporal distributions of ambient SO2 levels at high resolution. Among all the predictor variable categories, the weather category accounts for the largest portion of the overall variable importance (IM = 37.7%), followed by the category of emission inventory (IM = 25.9%). The weather conditions, such as temperature (IM = 12.4%), PBLH (IM = 5.3%), and relative humidity (IM = 4.7%), are associated with the atmospheric dispersion, dissipation and 5

Environmental Research 179 (2019) 108795

H. Zhang, et al.

Fig. 3. Spatial distributions of annual average ambient SO2 concentrations in (a) 2013, (b) 2014, (c) 2015, and (d) 2016.

SO2 pollution emerged in winter but was not a primary concern in the other seasons. The national total amount of SO2 emission was approximately 25% higher in winter than in summer, and the relative difference in the ambient SO2 concentrations between the two seasons was about 110%. Therefore, while the intensified anthropogenic emissions contributed to the elevated SO2 pollution in winter, the stagnant atmospheric conditions played a more important role. The predicted SO2 concentrations exhibited wider seasonal variation than the OMI-SO2 did (Fig. S10), as the shrinkage in cold season and the expansion in warm season of PBL caused the condensation and dilution of ambient SO2, respectively. The population-weighted average SO2 of China linearly reduced from 34.4 ± 18.2 in 2013 to 22.7 ± 11.1 μg/m3 in 2016, at the decreasing rate of −3.9 μg/m3/year (Figs. 3 and 5). The SO2 levels showed significantly decreasing trends in all the prefectures (P < 0.01; Table 2), which were also reflected in the temporal trend of OMI-SO2 (Fig. S8). The largest absolute decrease of SO2 occurred in the northern China, especially in Tianjin (−9.7 μg/m3/year), Shandong (−9.2 μg/ m3/year), and Hebei (−8.8 μg/m3/year). The southwestern China suffered severe SO2 pollution in 2013 (37.2 ± 7.4 μg/m3), after when the concentrations rapidly decreased to 19.0 ± 4.6 μg/m3 that was lower than the national average in 2016. These three municipalities, including Tianjin (−54%), Beijing (−53%), and Chongqing (−49%), showed the rapidest relative decrease. In contrast, the SO2 concentrations decreased relatively slower in Tibet (−1.9%), Taiwan (−14.3%), and Shanxi (−15%). The decreased coal consumption and enforced FGD installation contributed to the decline of SO2 emission and ambient concentration (van der A et al., 2017). In China, the air pollution prevention and control measures have been put forward to reduce SO2, including the Air Pollution Prevention and Control Action Plan, as well as the 12th and 13th Five Year Plans (SCC, 2011; SCC, 2013; SCC, 2017). Pollution prevention intends to avoid SO2 production in industrial and residential activities by for example replacing coal with cleaner energies. Pollution control refers to the actions of reducing SO2 emissions to the atmosphere, e.g., FGD installation. Consequently, both the coal consumption

the inventory (Table S8). The variable category of road density has low impact on the SO2 concentrations (IM = 1.9%), and the correlation between the distance to roads and the SO2 concentrations is low (R = 0.10), both of which reflect the fact that transportation contributes only 1.8% of the anthropogenic SO2 emissions (Table S8). 3.4. Spatiotemporal distributions of estimated SO2 concentrations The population-weighted average of SO2 concentration for China was 28.1 ± 14.0 μg/m3 during 2013–2016, and the spatial distribution exhibited high heterogeneity (Table 2 and Fig. 3). The predicted spatial pattern was highly similar to the OMI-SO2 distribution (Fig. S8). The northern China (45.1 ± 14.7 μg/m3) was the severest SO2 polluted region in China and throughout the world (Krotkov et al., 2016), especially in Shanxi (56.7 ± 5.8 μg/m3), Hebei (50.8 ± 4.9 μg/m3), and Shandong (46.5 ± 4.3 μg/m3) provinces. In contrast, the SO2 pollution was the lowest in Hainan (8.8 ± 3.1 μg/m3), Taiwan (9.1 ± 2.7 μg/m3), and Fujian (13.5 ± 3.3 μg/m3) provinces. Given the short lifetime of SO2 in the atmosphere (4–48 h), transport of SO2 emitted from industrial and residential activities (such as coal consumption) tends to influence the distribution of ambient SO2 at local rather than regional/provincial scale (Fioletov et al., 2015; Lee et al., 2011). The multiyear averages of SO2 concentrations for all the prefectures were correlated with those of coal consumptions (R = 0.51; Table S9). In addition, the spatial distributions of SO2 emission and concentrations were similar to each other, except for the Yangtze River Delta and the Pearl River Delta, where the emissions might be overestimated (Fig. S9). The predicted SO2 concentrations showed a strong seasonality across China during 2013–2016 (Fig. 4). The multiyear populationweighted averages of SO2 concentrations were the highest in winter (41.6 ± 26.4 μg/m3), followed by spring (25.9 ± 12.2 μg/m3), fall (25.7 ± 12.2 μg/m3), and summer (19.6 ± 8.3 μg/m3). Over all the seasons, the SO2 pollution was consistently higher in the northern China than in the other regions of China. In the southwestern China, as well as several major cities in northwestern and northeastern China, the 6

Environmental Research 179 (2019) 108795

H. Zhang, et al.

Fig. 4. Spatial distributions of average ambient SO2 concentrations in (a) spring, (b) summer, (c) fall, and (d) winter during 2013–2016.

sulfate aerosols, which are relatively stable and have more extensive impact on the environment through long-range transport (Perry et al., 1999). Most of the sulfate aerosols are finally removed from atmosphere through acid deposition (Chin and Jacob, 1996). In the previous studies, the air pollutant observations obtained from ground monitoring were used to investigate the long-range transport of sulfur and the classification of aerosols (Perry et al., 1999; Wang et al., 2006). The SO2 distributions derived by the RF-STK model supply the comprehensive information of the important precursor to sulfate aerosol. Spatial heterogeneity of SO2 distributions can be further characterized by using the geographically weighted regression based on both site-based SO2 observations and satellite-based SO2 predictions (Brunsdon et al., 1996; Hu et al., 2013). In addition, the estimated daily SO2 variations reveal the temporal evolution of plume dispersions. 3.5. Human exposure Fig. 5. Monthly deseasonalized trends of SO2 concentrations in the northern China and whole nation. The shaded bands represent the point-wise 99% confidence intervals.

On the basis of the ambient SO2 concentration predictions and the population distributions, the main results on the cumulative human exposure during 2013–2016 are summarized as follows (Fig. 6):

and SO2 emission decreased in most of the prefectures. There existed an interesting phenomenon in Shandong and Shanxi provinces, with the severest SO2 pollution and the largest coal consumption. Shandong experienced a large relative decrease in SO2 concentrations (−45%) and a substantial increase in coal consumption (9%). Comparatively, Shanxi exhibited a small relative reduction of SO2 (−15%) and a small decline of coal consumption (−3%) (Table S9). This may be attributed to that the number of FGD facilities in Shandong was approximately ten times of that in Shanxi, largely decreasing the ambient SO2 in Shandong. Stricter regulations on waste gas emission are required to further reduce the SO2 concentrations across China, especially in the northern China. The predicted SO2 distributions provide not only basic data for assessing and improving environmental management practices, but also important information for characterizing emission dispersions and evaluating transport models. Atmospheric SO2 can be converted to

• Exceedance of IT-1 (> 125 μg/m ): The population% exposed for > 16 nonattainment days decreased from 10% to 0%; • Exceedance of IT-2 (> 50 μg/m ): The population% exposed for > 40 nonattainment days decreased from 45% to 3%; • Exceedance of AQG (> 20 μg/m ): The population% exposed 3

3

3

for > 100 nonattainment days decreased from 86% to 48%;

where the exposure durations were set subject to the constant constrain on the exceeding exposure intensity over time (i.e., intensity × duration ≡ 2000 μg/m3 day)). The estimate uncertainties mainly originated from the residential air exchange rate, personal daily activities, and sub-grid variation in SO2 concentration. The results showed that China has generally achieved IT-1 and IT-2 but required more efforts to meet AQG. In particular, the population under high exposure to SO2 in 2016 were mostly distributed in the northern China, 7

Environmental Research 179 (2019) 108795

H. Zhang, et al.

Fig. 6. Human exposure assessments of daily ambient SO2 for (a) 2013, (b) 2014, (c) 2015, and (d) 2016. The daily Air Quality Guideline (AQG), Interim Target-1 (IT-1), and Interim Target-2 (IT-2) proposed by the World Health Organization (WHO) are 20, 50, and 125 μg/m3, respectively.

Declaration of competing interest

such as the Fenwei Plain located in Shanxi and Shaanxi provinces (Figs. 3 and S11). Under the current status, the “Three-year Plan on Defending the Blue Sky” issued in 2018 brought the Fenwei Plain into the key regions urgently requiring air pollution prevention and control (SCC, 2018).

The authors have no conflict of interest to disclose. Acknowledgements This study was supported by the National Natural Science Foundation of China (21607127), the Fundamental Research Funds for the Central Universities (YJ201765), and the Sichuan “1000 Plan” Young Scholar Program.

4. Conclusions The daily ambient SO2 concentrations are derived across China (0.1o grid) during 2013–2016 based on the satellite OMI retrievals and extensive ground observations. The SO2 concentrations predicted by the RF-STK model well capture the spatiotemporal variations of the ground observations. The ambient SO2 concentrations substantially decreased over these four years, showing the overall effectiveness of SO2 emission reduction. More efforts should be made to further reduce SO2 emission in order to achieve the AQG (20 μg/m3) all over China. Particularly, the regulation should be enforced in the Fenwei Plain located in Shanxi and Shaanxi provinces, where the SO2 pollution remained the severest. The comprehensive dataset of ambient SO2 provides important support for the air quality management and environmental epidemiological analyses. In the future, the TROPOMI SO2 product will be used to derive ambient SO2 distribution at higher spatial resolution.

Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.envres.2019.108795. References Breiman, L., 2001. Random forests. Mach. Learn. 45, 5–32. Brunsdon, C., et al., 1996. Geographically weighted regression: a method for exploring spatial nonstationarity. Geogr. Anal. 28, 281–298. Carn, S.A., et al., 2005. Quantifying tropospheric volcanic emissions with AIRS: the 2002 eruption of Mt. Etna (Italy). Geophys. Res. Lett. 32. Charlson, R.J., et al., 1992. Climate forcing by anthropogenic aerosols. Science 255, 423–430. Cheng, Y., et al., 2016. Reactive nitrogen chemistry in aerosol water as a source of sulfate during haze events in China. Sci. Adv. 2, e1601530. Chin, M., Jacob, D.J., 1996. Anthropogenic and natural contributions to tropospheric sulfate: a global model analysis. J. Geophys. Res. Atmos. 101, 18691–18699. CIESIN, 2016. Gridded Population of the World, Version 4 (GPWv4): Population Count. NASA Socioeconomic Data and Applications Center (SEDAC), Palisades, NY. Clarisse, L., et al., 2008. Tracking and quantifying volcanic SO2 with IASI, the september 2007 eruption at jebel at tair. Atmos. Chem. Phys. 8, 7723–7734.

Funding This study was supported by the National Natural Science Foundation of China (21607127), the Fundamental Research Funds for the Central Universities (YJ201765), and the Sichuan “1000 Plan” Young Scholar Program. 8

Environmental Research 179 (2019) 108795

H. Zhang, et al.

HomePage.aspx?id=N2017110016&name=YCXME&floor=1, Accessed date: 30 March 2018. NBSC, 2017b. China statistical Yearbook on environment. http://www.stats.gov.cn/ztjc/ ztsj/hjtjzl/, Accessed date: 30 March 2018. OMI Team, 2012. Ozone monitoring instrument (OMI) data user's guide. https:// docserver.gesdisc.eosdis.nasa.gov/repository/Mission/OMI/3.3_ ScienceDataProductDocumentation/3.3.2_ProductRequirements_Designs/README. OMI_DUG.pdf, Accessed date: 10 April 2018. OpenStreetMap contributors, 2015. Planet OSM. http://planet.openstreetmap.org, Accessed date: 10 April 2018. Pedregosa, F., et al., 2011. Scikit-learn: machine learning in Python. J. Mach. Learn. Res. 12, 2825–2830. Perry, K.D., et al., 1999. Long-range transport of anthropogenic aerosols to the national oceanic and atmospheric administration baseline station at mauna loa observatory, Hawaii. J. Geophys. Res. Atmos. 104, 18521–18533. Philip, S., et al., 2014. Global chemical composition of ambient fine particulate matter for exposure assessment. Environ. Sci. Technol. 48, 13060–13068. SCC, 2011. Notice of the State Council on printing and distributing the comprehensive work program of 12th Five-Year Plan for energy saving and emission reduction. http://www.gov.cn/zhengce/content/2011-09/07/content_1384.htm, Accessed date: 10 April 2018. SCC, 2013. Notice of the state council on printing and distributing the air pollution prevention and control action plan. http://www.gov.cn/zhengce/content/2013-09/ 13/content_4561.htm, Accessed date: 10 April 2018. SCC, 2017. Notice of the State Council on printing and distributing the comprehensive work program of 13th Five-Year Plan for energy saving and emission reduction. http://www.gov.cn/zhengce/content/2017-01/05/content_5156789.htm, Accessed date: 10 April 2018. SCC, 2018. Notice of the state council on printing and distributing the three-year plan on defending the Blue Sky. http://www.gov.cn/zhengce/content/2018-07/03/content_ 5303158.htm, Accessed date: 11 September 2018. Song, W.Z., et al., 2014. A satellite-based geographically weighted regression model for regional PM2.5 estimation over the Pearl River Delta region in China. Remote Sens. Environ. 154, 1–7. Theys, N., et al., 2013. Volcanic SO2 fluxes derived from satellite data: a survey using OMI, GOME-2, IASI and MODIS. Atmos. Chem. Phys. 13, 5945–5968. Theys, N., et al., 2017. Sulfur dioxide retrievals from TROPOMI onboard Sentinel-5 Precursor: algorithm theoretical basis. Atmos. Meas. Tech. 10, 119–153. van der A, R.J., et al., 2017. Cleaning up the air: effectiveness of air quality policy for SO2 and NOX emissions in China. Atmos. Chem. Phys. 17, 1775–1789. van Donkelaar, A., et al., 2008. Analysis of aircraft and satellite measurements from the Intercontinental Chemical Transport Experiment (INTEX-B) to quantify long-range transport of East Asian sulfur to Canada. Atmos. Chem. Phys. 8, 2999–3014. van Donkelaar, A., et al., 2012. Improving the accuracy of daily satellite-derived groundlevel fine aerosol concentration estimates for North America. Environ. Sci. Technol. 46, 11971–11978. Vienneau, D., et al., 2013. Western European land use regression incorporating satelliteand ground-based measurements of NO2 and PM10. Environ. Sci. Technol. 47, 13555–13564. Wang, Y., et al., 2006. The ion chemistry, seasonal cycle, and sources of PM2.5 and TSP aerosol in Shanghai. Atmos. Environ. 40, 2935–2952. Wang, Y., et al., 2014. Spatial and temporal variations of six criteria air pollutants in 31 provincial capital cities in China during 2013-2014. Environ. Int. 73, 413–422. WHO, 2006. Air quality guidelines: global update 2005. Particulate matter, ozone, nitrogen dioxide and sulfur dioxide. http://www.euro.who.int/__data/assets/pdf_file/ 0005/78638/E90038.pdf, Accessed date: 4 October 2018. Xu, Q., 2001. Abrupt change of the mid-summer climate in central east China by the influence of atmospheric pollution. Atmos. Environ. 35, 5029–5040. Xu, X., et al., 2010. Temporal-spatial variations of tropospheric SO2 over China using SCIAMACHY satellite observations. In: 18th International Conference on Geoinformatics, pp. 1–5 2010. Yan, H.H., et al., 2014. SO2 columns over China: temporal and spatial variations using OMI and GOME-2 observations. In: 35th International Symposium on Remote Sensing of Environment (Isrse35) Vol. 17, 012027. Yang, X., et al., 2017. The impact of anthropogenic emissions and meteorological conditions on the spatial variation of ambient SO2 concentrations: a panel study of 113 Chinese cities. Sci. Total Environ. 584–585, 318–328. Zhan, Y., et al., 2017. Spatiotemporal prediction of continuous daily PM2.5 concentrations across China using a spatially explicit machine learning algorithm. Atmos. Environ. 155, 129–139. Zhan, Y., et al., 2018a. Spatiotemporal prediction of daily ambient ozone levels across China using random forest for human exposure assessment. Environ. Pollut. 233, 464–473. Zhan, Y., et al., 2018b. Satellite-based estimates of daily NO2 exposure in China using hybrid random forest and spatiotemporal kriging model. Environ. Sci. Technol. 52, 4180–4189. Zhang, Y., et al., 2016. Continuation of long-term global SO2 pollution monitoring from OMI to OMPS. Atmos. Meas. Tech. Discuss. 1–21 2016. Zhao, S., et al., 2016. Annual and diurnal variations of gaseous and particulate pollutants in 31 provincial capital cities based on in situ air quality monitoring data from China National Environmental Monitoring Center. Environ. Int. 86, 92–106.

Clerbaux, C., et al., 2008. Measurements of SO2 profiles in volcanic plumes from the NASA tropospheric emission spectrometer (TES). Geophys. Res. Lett. 35. Cleveland, R.B., 1990. STL : a seasonal-trend decomposition procedure based on loess. J. Off. Stat. 6, 3–33. CMA, 2015. China meteorology data. http://data.cma.cn/, Accessed date: 4 October 2018. Dentener, F., et al., 2006. Nitrogen and sulfur deposition on regional and global scales: a multimodel evaluation. Glob. Biogeochem. Cycles 20, n/a–n/a. Deutsch, C.V., Journel, A.G., 1998. GSLIB Geostatistical Software Library and User's Guide. Oxford University Press, New York. Dey, S., Tripathi, S.N., 2008. Aerosol direct radiative effects over Kanpur in the IndoGangetic basin, northern India: long-term (2001–2005) observations and implications to regional climate. J. Geophys. Res. 113. Didan, K., et al., 2015. MODIS Vegetation Index User's Guide (MOD13 Series). (Collection 6) Version 3.00. Eisinger, M., Burrows, J.P., 1998. Tropospheric sulfur dioxide observed by the ERS-2 GOME instrument. Geophys. Res. Lett. 25, 4177–4180. EPAROC, 2015. Taiwan air quality monitoring network. http://taqm.epa.gov.tw, Accessed date: 10 April 2018. EPDHK, 2015. Hong Kong air quality monitoring data. http://epic.epd.gov.hk/EPICDI/ air/station/, Accessed date: 10 April 2018. Fioletov, V.E., et al., 2011. Estimation of SO2 emissions using OMI retrievals. Geophys. Res. Lett. 38 n/a-n/a. Fioletov, V.E., et al., 2013. Application of OMI, SCIAMACHY, and GOME-2 satellite SO2 retrievals for detection of large emission sources. J. Geophys. Res. Atmos. 118 (11) 399-11,418. Fioletov, V.E., et al., 2015. Lifetimes and emissions of SO2 from point sources estimated from OMI. Geophys. Res. Lett. 42, 1969–1976. Fioletov, V.E., et al., 2016. A global catalogue of large SO2 sources and emissions derived from the Ozone Monitoring Instrument. Atmos. Chem. Phys. 16, 11497–11519. GES DISC, 2015. OMSO2e: OMI/Aura Sulfur Dioxide (SO2) Total Column L3 1 day Best Pixel in 0.25 degree x 0.25 degree V3. https://disc.gsfc.nasa.gov/datasets/OMSO2e_ V003/summary?keywords=OMI, Accessed date: 10 April 2018. GMAO, 2015. MERRA-2 tavg1_2d_flx_Nx: 2d,1-Hourly, Time-Averaged, Single-Level, Assimilation, Surface Flux Diagnostics V5.12.4. Goddard Earth Sciences Data and Information Services Center (GES DISC), Greenbelt, MD, USA. Goodfellow, I., et al., 2016. Deep Learning. MIT press. Graler, B., et al., 2016. Spatio-Temporal Interpolation using gstat. RFID J. 8, 204–218. Hastie, T., et al., 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer-Verlag, New York. Holben, B.N., 1986. Characteristics of maximum-value composite images from temporal avhrr data. Int. J. Remote Sens. 7, 1417–1434. Hu, X., et al., 2013. Estimating ground-level PM2.5 concentrations in the southeastern U.S. using geographically weighted regression. Environ. Res. 121, 1–10. Iwasawa, S., et al., 2009. Effects of SO2 on respiratory system of adult Miyakejima resident 2 years after returning to the island. J. Occup. Health. 51, 38–47. Jarvis, A., et al., 2008. Hole-filled seamless SRTM data V4.1. International Centre for Tropical Agriculture (ciat). http://srtm.csi.cgiar.org, Accessed date: 10 April 2018. Jun, C., et al., 2014. China: open access to Earth land-cover map. Nature 514, 434. Krotkov, N.A., et al., 2016. Aura OMI observations of regional SO2 and NO2 pollution changes from 2005 to 2015. Atmos. Chem. Phys. 16, 4605–4629. Krueger, A.J., 1983. Sighting of el chichón sulfur dioxide clouds with the nimbus 7 total ozone mapping spectrometer. Science 220, 1377–1379. Lee, C., et al., 2011. SO2 emissions and lifetimes: estimates from inverse modeling using in situ and global, space-based (SCIAMACHY and OMI) observations. J. Geophys. Res. Atmos. 116. Lee, C.J., et al., 2015. Response of global particulate-matter-related mortality to changes in local precursor emissions. Environ. Sci. Technol. 49, 4335–4344. Li, C., et al., 2013. A fast and sensitive new satellite SO2 retrieval algorithm based on principal component analysis: application to the ozone monitoring instrument. Geophys. Res. Lett. 40, 6314–6318. Li, C., et al., 2017a. India is overtaking China as the world's largest emitter of anthropogenic sulfur dioxide. Sci. Rep. 7, 14304. Li, M., et al., 2017b. MIX: a mosaic Asian anthropogenic emission inventory under the international collaboration framework of the MICS-Asia and HTAP. Atmos. Chem. Phys. 17, 935–963. Li, R., et al., 2019. Satellite-based prediction of daily SO2 exposure across China using a high-quality random forest-spatiotemporal Kriging (RF-STK) model for health risk assessment. Atmos. Environ. 208, 10–19. Liu, Y., et al., 2009. Estimating regional spatial and temporal variability of PM2.5 concentrations using satellite data, meteorology, and land use information. Environ. Health Perspect. 117, 886–892. Lu, Z., et al., 2010. Sulfur dioxide emissions in China and sulfur trends in East Asia since 2000. Atmos. Chem. Phys. 10, 6311–6331. Lu, Z., et al., 2013. Ozone monitoring instrument observations of interannual increases in SO2 emissions from Indian coal-fired power plants during 2005-2012. Environ. Sci. Technol. 47, 13993–14000. Luvsan, M.E., et al., 2012. The influence of emission sources and meteorological conditions on SO2 pollution in Mongolia. Atmos. Environ. 61, 542–549. MEPC, 2015. Air quality daily report for China. http://datacenter.mee.gov.cn/websjzx/ queryIndex.vm, Accessed date: 10 April 2018. NBSC, 2017a. China energy statistical Yearbook. http://tongji.cnki.net/kns55/Navi/

9