Long-Term Dynamics of Atmospheric Sulfur Dioxide in Urban and Rural Regions of China: Urbanization and Policy Impacts

Abstract

High levels of sulfur dioxide (SO₂) due to human activities pose a serious air pollution issue in China, especially in urban agglomerations. However, limited research has investigated the impact of anthropogenic emissions on higher SO₂ concentrations in urban regions compared to rural areas in China. Here, we analyzed the trends in SO₂ concentrations from 1980 to 2021 in China using the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) dataset. SO₂ column concentrations from the Copernicus Atmosphere Monitoring Service (CAMS) and the Ozone Monitoring Instrument (OMI) during the years 2007–2021 were also examined for validation and comparison purposes. Eight representative areas, including four urban regions (Pearl River Delta [PRD], Beijing-Tianjin-Hebei [BTH], Yangtze River Delta [YRD], and Sichuan Basin [SCB]) and four rural regions (Northeast Region [NER], Mongolian Region [MR], West Region [WR], and Tibetan Plateau Region [TR]) were selected for the analysis. Overall, a significant but fluctuating increase in SO₂ concentrations over China was observed during 1980–2021. During 1980–1997 and 2000–2010, there was an increase in SO₂ concentration, while during 1997–2000 and 2010–2021, a decreasing trend was observed. The average increase in SO₂ concentration was approximately 16 times higher in urban regions than in the rural background. We also found that SO₂ dynamics were highly associated with expansion of urban areas, population density, and gross domestic product. Nonetheless, since 2007, SO₂ concentrations have exhibited a downward trend, which is mainly attributed to the air pollution policies implemented by the Chinese government. Our findings highlight the need for further studies on the impact of SO₂ on regional climate change in China.

1. Introduction

Sulfur dioxide (SO₂), a main atmospheric gas pollutant, plays a key role in the formation of sulfate aerosols and contributes to acid rain [1,2,3,4]. Thus, the oxidation of SO₂ products (sulfuric acid and sulfate) can adversely affect the urban environment, air quality, the Earth’s climate system, and human health [5,6,7,8,9]. Anthropogenic emissions, particularly those from fossil fuel combustion, are a major source of SO₂ in economically developing countries [10,11,12,13,14]. SO₂ emissions are the primary cause of acid deposition in North America, Europe, and China [15,16], and SO₂ is also a key precursor of secondary sulfate aerosols in ambient air [17,18,19,20]. Over the last two decades, the distribution of SO₂ emissions has changed dramatically on a global scale, with major changes detected over the Middle East, East Asia, and the Indian subcontinent [12,17,21,22,23]. Volcanic eruptions can also disperse large amounts of volcanic ash and gases such as SO₂ into the atmosphere, which can influence radiative transfer, climate, ecosystems, agriculture, and aviation, as well as air quality and health [24,25].

Urbanization and economic growth in China have resulted in a significant increase in SO₂ emissions over recent decades [26,27]. In addition, chemical processes and meteorological conditions play significant roles in SO₂ concentrations in the atmosphere over China [24,25]. In this regard, anthropogenic SO₂ concentration data from various sources such as the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) [23,28,29], Copernicus Atmosphere Monitoring Service (CAMS) [17,30] and Ozone Monitoring Instrument (OMI) [8,17] have been widely used to assess the spatiotemporal variation of SO₂ concentrations in different regions [14,31,32]. Of note, the real SO₂ emission trends in developed countries, including different regions of China, differ largely from the projections. Previous studies have concluded that India is overtaking China as the world’s largest emitter of anthropogenic SO₂ [17,33].

Previous studies have enhanced our current understanding about the variations in SO₂ emissions in China. An early study estimated that China emits approximately 15 million metric tons of SO₂ annually [34], while China was the largest SO₂ emitter in East Asia between the years 1975–1987, presenting a significant increasing trend in SO₂ emissions during this period [18]. SO₂ emissions increased over Asia by 119% during 1980–2003, which can be attributed to large increases in SO₂ emissions in China associated with its rapid growth in total energy consumption [35]. During the years 1949–2015, an increasing trend in SO₂ emissions was identified in China [26], although 2006 was a turning point. In 2015, the SO₂ emissions in China were approximately 27 million metric tons [26]. In addition, different spatio-temporal dependencies of SO₂ emissions were observed over most parts of eastern China during 1998–2014 [25], by using monthly total emissions for Eastern China provided by Peking University emissions inventory (PKU; http://inventory.pku.edu.cn/, accessed on 1 October 2023). These inventories include six major air pollutants, with a spatial resolution of 10 km. On the other hand, SO₂ columns from different satellite sensors, including OMI, were also analyzed over China, revealing negative trends in OMI SO₂ concentrations over 90% of the power plant sites in Chinese megacities during 2007–2014 [36], while the SO₂ concentrations decreased by 80.4% from the years 2007 to 2019 in Northern China [37]. The changes in carbonaceous aerosols were in good agreement with the changes in Aerosol Optical Depth (AOD) and SO₂ over India and China during the years 1996–2010 [38]. Anthropogenic aerosols have recently decreased in some parts of China due to Chinese emissions control policies [28]. Using an approach of adoption curve modeling, there is a projection for the reduction potential of SO₂ emissions from the years 2020 to 2050 [39]. Meteorological conditions and changes in emission rate may also explain various distributions of atmospheric SO₂ loading during different seasons in China [12,24,25].

However, important issues remain to be addressed to further advance our understanding about spatio-temporal variations and trends of SO₂ concentrations in China. First, how different are the long-term (1980–2021) SO₂ trends in urban regions compared to those in rural regions in China? Second, besides previous studies analyzing the impact of meteorological conditions on anthropogenic SO₂ changes, a significant environmental issue that warrants discussion is the extent to which human activity is responsible for SO₂ trends in different regions of China and, in particular, how urbanization, gross domestic product (GDP), and population affect SO₂ variability over urban regions compared to rural areas of China. Third, the impact of the air pollution mitigation policies implemented by the Chinese government on regional SO₂ trends across China is lacking.

Therefore, in view of the above issues, in this study, we present a comprehensive long-term trend analysis of SO₂ concentrations over four main urban and four main rural regions in China during 1980–2021, using SO₂ retrievals from the MERRA-2, CAMS, and OMI to better assess anthropogenic-related SO₂ changes. To examine the reliability of the computed SO₂ trends, SO₂ concentrations from ground air pollution stations were used for comparison. Urbanization, GDP, and population data were used to investigate the changes in SO₂ concentrations in response to human activities in the selected regions.

2. Materials and Methods

2.1. Regions of Interest

China is a global hotspot for high concentrations of aerosol particles [40,41]. In this study, eight representative regions of China were selected (Figure 1) to analyze the long-term dynamics of SO₂ concentrations via satellite observations and reanalysis. The four regions of interest (ROIs) represent urban areas with high population density and loading of anthropogenic emissions [25,42]; these regions are the Yangtze River Delta (YRD, 117°E_122°E, 28°N_33°N), Beijing-Tianjin-Hebei (BTH, 113°E_119°E, 36°N_42°N), Pearl River Delta (PRD, 108°E_114°E, 20°N_26°N), and the Sichuan Basin (SCB, 103°E_108°E, 26°N_31°N) [43,44,45]. Meanwhile, the four other examined ROIs representing rural regions in this study are largely covered by forests or desert (Figure 1), and are arid areas with less population density and low GDP (Figure S1); these are the Northeast Region (NER, 118°E_129 ° E, 49 °N_53 °N), Mongolian Region (MR, 97°E_111°E, 38°N_42°N), West Region (WR, 81°E_89°E, 37°N_43°N), and Tibetan Plateau Region (TR, 80°E_98°E, 28°N_33°N) [46,47,48].

2.2. Database

2.2.1. SO₂ Data

The Chinese National Environmental Monitoring Center (CNEMC) provides ground observations of SO₂ concentrations. Many ground stations were established in China in 2012. In this study, 100 sites with complete ground-based SO₂ concentrations data (2013–2021) were used (Figure 1) to evaluate the MERRA-2 SO₂ products, based on monthly averages. The SO₂ observation sites were also used in the trend analysis. Further details regarding the Chinese observational data can be found at (http://www.cnemc.cn/, accessed on 3 March 2020) (

Table 1. Data used.

Data Type	Parameters	Period	Date Sources
SO2 Data	SO2 concentrations	1980–2021	Monthly MERRA-2 reanalysis dataset at 0.5° × 0.625° spatial resolution https://disc.gsfc.nasa.gov, accessed on 15 March 2023
	Total column SO2	2007–2021	Monthly CAMS reanalysis dataset at 0.75° × 0.75° spatial resolution https://ecmwf.int/, accessed on 15 May 2023
	Total column SO2	2007–2021	Monthly OMI satellite dataset at 0.25° × 0.13° spatial resolution https://disc.gsfc.nasa.gov, accessed on 15 May 2023
	SO2 concentrations	2013–2021	Monthly ground observation http://www.cnemc.cn/, accessed on 15 May 2023
Driving factors	Population count	1990–2020	5-year intervals of GPW V 3 and 4 http://sedac.ciesin.columbia.edu, accessed on 15 March 2023
	LULC	1990–2020	5-year intervals of six primary land use categories http://www.resdc.cn/, accessed on 15 March 2023
	GDP	1990–2020	5-year intervals of grid map using Chinese governmental data https://www.resdc.cn/, accessed on 15 March 2023

).

MERRA-2 is the second version of the atmospheric reanalysis of the modern satellite era and uses the Atmospheric Data Assimilation System (ADAS) and the Goddard Earth Observing System Model, Version 5 (GEOS-5). MERRA-2 utilizes a radiatively coupled version of the Goddard Chemistry, Aerosol, Radiation and Transport (GOCART) aerosol module [49,50,51,52] to simulate 15 externally mixed aerosol mass mixing ratio tracers, including SO₂ concentrations. In this study, SO₂ concentrations, based on monthly means from MERRA-2, were used for our assessment during the 1980–2021 period (Table 1). Based on the SO₂ ground-based observations from 2013 to 2021, we evaluated the uncertainty of MERRA-2 SO₂ products in China based on SO₂ ground-based observations from 2013 to 2021. The results showed a significant overall agreement between MERRA-2 SO₂ and ground-based SO₂ across China, with an acceptable bias for annual mean data, and correlation coefficient of R = 0.76 [root-mean square error (RMSE): 8.3 µg m⁻³; mean absolute error (MAE): 8.9 µg m⁻³; fractional gross error (FGE): 1.39; mean fractional error (MFE): 45.3%; and relative mean bias (RMB): 1.11]. For more details on the comparison and validation processes, please refer to Text S1 and Figures S1 and S2, in the Supplementary Material.

In addition to MERRA-2, this study utilized the total column SO₂ concentration over China from the CAMS model reanalysis initiated by the European Center for Medium-range Weather Forecasts (ECMWF) (https://ecmwf.int/, accessed on 15 December 2022). CAMS combines observation data assimilation and forecasting systems to analyze daily products and 5-d forecasts of atmospheric composition, including forecasts of SO₂ in near real time [17,53,54]. Monthly means of the SO₂ total column datasets (at 0.75° × 0.75° longitude–latitude grid of resolution) (Table 1), provided by CAMS during the period 2007–2021, were used in the current study to improve the assessment, including annual and seasonal variabilities and inter-comparisons with the MERRA-2 dataset.

The OMI sensor is onboard the Earth Observing System (EOS) Aura satellite [55]. In the spectral zoom-in measurement mode, OMI enables the retrieval of aerosols, OClO, HCHO, ozone, and particularly SO₂ observations at a spatial resolution of 25 km × 13 km at the nadir [8], with near global coverage [55]. The OMI columnar SO₂ is given in Dobson units (DU), with 1 DU = 2.69 × 10¹⁶ molecules·cm⁻² [8,17]. The single retrieval error of SO₂ was estimated to be about 70–150% [21]. However, a previous study demonstrated the validity of using OMI SO₂ data in China [8]. In this study, the OMI monthly mean SO₂ total column datasets from 2007 to 2021 were used from (https://disc.sci.gsfc.nasa.gov/, accessed on 15 May 2023).

2.2.2. Key Driving Factors

Driving factors associated with human activity play a key role in regional-scale SO₂ changes. More specifically, population has been recognized to be one of the main factors influencing changes in anthropogenic emissions. In addition, urban-land expansion presents a critical influence on understanding changes in urban emission rates [56,57]. GDP is used to describe the economic development of a country, which in China often drives construction activities and increases fossil fuel consumption and anthropogenic emissions [58,59].

In this study, population, urbanization, and GDP were selected as the main driving factors for analyzing the dynamics of SO₂ in China. Population data were taken from the Gridded Population of the World version 3 (GPWv3) from 1990 to 2000, and version 4 (GPWv4) from 2000 to 2020, provided by the NASA Socioeconomic Data and Applications Center. Both versions were at 2.5-arcminute resolution, with 5-year intervals [25], and further information related to the population are available at (http://sedac.ciesin.columbia.edu/, accessed on 15 March 2023) (Table 1). Urban-land expansion is presented as an urban area percentage, which refers to the ratio between the urban and total land areas, serving as a measure of the degree of urbanization [60].

Land-use data in the ROIs were provided by the Chinese Academy of Sciences (CAS) through the National Resources and Environment Database, including a series of five years of temporal resolution from late 1990 to 2019. There are six primary land use categories (http://www.resdc.cn/, accessed on 15 March 2023). In this study, artificial surface classes were extracted and used from 1990 to 2020 to determine the influence of urbanization. The data were highly recommended as being of high quality and have been used in previous studies [61,62].

The GDP data used in this study were provided by the CAS covering the period from 1990 to 2019. The datasets were bilinearly interpolated to a grid map using Chinese governmental data [58,63].

2.3. Trend and Regression Analysis

Trend analysis was implemented using monthly timeseries of SO₂ concentrations from MERRA-2, and the total column from CAMS and OMI over the examined ROIs, using the least-squares linear regression technique, which is often used for long-term analyses [29,64]. Before the linear regression analysis, the MERRA-2, CAMS, and OMI retrievals were bilinearly interpolated to spatial grids with 1° × 1° resolution. For each grid, the annual cycle of monthly averages was removed [65,66] using annual average data. The Student’s t test was applied, and the 95% confidence level (p < 0.05) was set as the cut-off for statistical significance [67,68].

To better understand the effects of the total population (TP), GDP, and urbanization percentage (UP) on SO₂ concentrations in each ROI, MERRA-2 SO₂ concentrations obtained during the years 1990–2021 were averaged based on 3-year means. For example, the year 2000 value represents the average of the values in 1999, 2000, and 2001, with five-year temporal intervals. Subsequently, the averaged MERRA-2 SO₂ concentrations were linked to population, urbanization, and GDP for 1990, 1995, 2000, 2005, 2010, 2015, and 2020. The strength of the relationships was measured using the Pearson correlation coefficient (R). An ordinary least squares regression model (OLS) [59,69,70,71,72] was used in this study. The OLS method is given in Equation (1)

$$y={\beta }_{ }{x}_{ }+{\beta }_0+ϵ,$$

(1)

where y is the extracted mean SO₂ concentration of an ROI as the dependent variable; and $\beta \mathrm{x}$ is the regression coefficient of the independent variables; for example, TP, GDP, and UP. ${\beta }_0$ is the value of the intercept, and ε refers to the standard error. Geographically weighted regression (GWR) models explain the spatial variation in the strength of the relationships between factors across selected geolocation regions [59,73]. Therefore, the GWR model was used in a form similar to that of the OLS model. Nevertheless, the selected parameters differ according to the geographical location. The GWR is given in Equation (2) as:

$$y_j={\beta }_o\left(u_j,v_j\right)+{\sum }_{i=1}^p{\beta }_i\left(u_j,v_j\right)x_{ij}+{\epsilon }_j$$

(2)

where ${\beta }_o\left(u_j,v_j\right)$ is the intercept value for position (j) at coordinates $\left(u_j,v_j\right)$; and ${\beta }_i\left(u_j,v_j\right)$ is the estimated independent variable of local parameters, such as TP, GDP, and UP of (x_i) at location (j). In addition to the localization of the parameters, spatial autocorrelation was accurately included in the GWR model. A defining feature of GWR is that the selected parameters are allowed to vary across spaces to measure non-static spatial relationships. Using different bandwidths and a representation of the matrix for each period, the estimated value of the selected variable can be expressed by Equation (3)

$${\beta }_f\left(u_j,v_j\right)={\left(X^{T}W\left(u_j,v_j\right)X\right)}^{-1}{X}^{T}W\left(u_j,v_j\right)y$$

(3)

where T is the matrix transpose operation; ${\beta }_f\left(u_j,v_j\right)$ is the weight of the specific location of the variable; and $W\left(u_j,v_j\right)$ is a (number) by (number) spatial weight matrix at location j. All the calculations of GWR were performed in ArcGIS software version 10.8 via a model (GTWR Beta 1.0) developed by [74,75]. To examine the relative contributions (%) of the selected variables to SO₂ concentration changes, the method of Lindeman, Merenda, and Gold (LMG) [12,28,76] was used. The LMG can be calculated as.

$$LMG(x_{i })=\frac{1}{p!} {\sum }_{\mathrm{r} \mathrm{p}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}^{ }seq{R}^2\left({\{x}_{i }\right\}\left\{r\right\})$$

(4)

the rth (r = 1, 2, …, p!) represented by r, and $seq{R}^2\left({\{x}_i\right\}\left\{r\right\})$ denotes the sequential of sum squares for the r of x_i in the ordering of the r in the rth.

3. Results and Discussion

3.1. Spatial Patterns of SO₂ Concentrations

Figure 2 displays the spatial patterns of the annual mean SO₂ calculated from the MERRA-2, CAMS, and OMI data over China during 1980–2021, 2007–2021, and 2007–2021, respectively. The three datasets showed similar patterns of SO₂ over China despite the different time periods, with the highest values observed over the Eastern China plains, near large urban agglomerations such as PRD, YRD, BTH, and SCB. Due to different units in SO₂ concentrations between the three datasets, only a qualitative visualized comparison can be extracted from the spatial distribution maps (Figure 2a–c). Figure 2d shows the seasonal mean SO₂ concentrations (in μg m⁻³) over the eight selected regions in China, as analyzed from MERRA–2 during 1980–2021.

High SO₂ loading was mainly observed in the urban regions of China. The YRD and BTH regions located in eastern China are highly dominated by anthropogenic emissions [43,45,77,78] and exhibited a very high amount of SO₂. In contrast, the SO₂ concentrations were the lowest in the TR and WR remote areas, as expected. The multi-year average of MERRA-2 SO₂ concentrations is the highest in YRD (28.9 ± 8.5 μg m⁻³) followed by BTH (20.9 ± 9.1 μg m⁻³), PRD (14.9 ± 6.1 μg m⁻³), and SCB (12.4 ± 6.8 μg m⁻³). Industrialization and urbanization play a key role in the increase in SO₂ emissions in eastern China [8]. SO₂ from areas of high anthropogenic and industrial emissions tends to affect the SO₂ concentration distribution at a local rather than regional scale [79,80]. Over the four selected rural regions, the annual average SO₂ concentrations were highest in NER (3.4 ± 2.6 μg m⁻³) and lowest in TR (0.26 ± 0.2 μg m⁻³). Note that the multi-year average of MERRA-2 SO₂ concentrations over the four selected rural regions is about 11 times lower than that over the corresponding urban regions.

A comparison of the three products showed that the mean spatial distribution of the MERRA-2 SO₂ concentrations agreed well with the results of CAMS and OMI in most regions of China, with relatively small differences. On a seasonal timescale, high SO₂ loading was mainly observed in winter (December, January, and February; DJF) and autumn (September, October, and November; SON), whereas the concentration of MERRA–2 SO₂ was relatively low in spring (March, April, and May; MAM) and summer (June, July, and August; JJA) (Figure 3). Previous studies have shown that approximately 35% of the annual average of SO₂ concentrations over East China occur during the winter, whereas the summer season contributes only 15% to the annual SO₂ loading [24]. Previous study [8] also reported highest SO₂ concentrations in winter in China from 2013 to 2016.

The current results of the highest seasonal mean SO₂ distribution in winter, followed by autumn, were also similar among the three SO₂ datasets (Figure 3). Particulate matter (PM), including sulfate aerosols produced by human activities in winter and autumn, is the main contributor to high aerosol loading in eastern China [81]. The higher concentrations of SO₂ in autumn and winter can be attributed mainly to the increased demand for coal consumption, resulting in an increase in SO₂ loading along with increased black carbon, NO_x and organic aerosols [8,27,31,37]. Moreover, lower humidity and cold temperatures, attributed to the winter East Asian monsoon, jointly weaken wet deposition and oxidation rates, extending the lifetime of SO₂ molecules [24], thus favoring the SO₂ accumulation. SO₂ loading is relatively low in spring and summer because of seasonal variations in atmospheric turbulence driven by solar radiation, which makes wet removal effective during the monsoon period [12,82,83].

3.2. Trends in SO₂ Concentrations

To better assess the temporal variation and trends of MERRA-2 SO₂ concentrations, normalized annual anomalies were examined during the period 1980–2021 over the selected urban (YRD, BTH, PRD, and SCB) and rural (NER, MR, WR, and TR) regions (Figure 4). The results show that SO₂ concentrations have experienced relatively similar and significant (p < 0.05) positive trends in both urban and rural regions during 1980–2021. Significant increasing trends were observed in YRD (0.71 µg m⁻³ yr⁻¹), BTH (0.5 µg m⁻³ yr⁻¹), PRD (0.52 µg m⁻³ yr⁻¹), and SCB (0.35 µg m⁻³ yr⁻¹). The percentage of relative change in SO₂ concentrations between 1980 and 2021 exceeded 100% over the selected urban regions (

Table 2. Trend and relative change percentage values from the MERRA-2 SO₂ (µg m⁻³ yr⁻¹), CAMS SO₂ (mg m⁻³ yr⁻¹), and OMI SO₂ (DU yr⁻¹) datasets in the eight selected regions during different periods. Bold text refers to significant trends (p < 0.05).

	Urban Regions				Rural Regions
Regions/Data and Periods	YRD	BTH	PRD	SCB	NER	MR	WR	TR
	Trend/Change	Trend/Change	Trend/Change	Trend/Change	Trend/Change	Trend/Change	Trend/Change	Trend/Change
MERRA-2 1980–2021	0.71	0.50	0.52	0.35	0.04	0.06	0.01	0.012
	137%	142%	232%	157%	38%	176%	33%	256%
MERRA-2 1980–1997	0.57	0.44	0.37	0.22	0.01	0.03	0.007	0.008
	50%	65%	89%	45%	5%	66%	120 %	190%
MERRA-2 1997–2001	−0.92	−0.91	−0.38	−0.33	−0.10	−0.01	−0.01	0.03
	−25%	−28%	−11%	−18%	−40%	−6%	−13%	80%
MERRA-2 2001–2010	2.42	1.68	1.52	1.23	0.21	0.16	0.04	0.014
	98%	89%	119%	129%	74%	73%	40%	34%
MERRA-2 2007–2021	−0.04	−0.05	−0.01	−0.03	−0.02	0.02	−0.008	0.003
	−7%	−4%	−10%	−8%	−14.2%	5%	−15%	7%
CAMS 2007–2021	−0.15	0.07	−0.05	−0.06	−0.04	−0.008	−0.01	−0.004
	−14%	13%	−10%	−13%	−24%	−5%	−5%	−6%
OMI 2007–2021	−0.024	−0.03	−0.002	−0.02	0.001	−0.001	−0.001	−0.001
	−51%	−65%	−18%	−54%	−4%	−20%	−6%	−11%

), with the highest percentage increase (232%) detected in PRD. Over the four selected rural regions, statistically significant, but weaker, increasing trends were observed during 1980–2021 with rates of 0.04 µg m⁻³ yr⁻¹ in MER, 0.06 µg m⁻³ yr⁻¹ in MR, 0.01 µg m⁻³ yr⁻¹ in WR, and 0.012 µg m⁻³ yr⁻¹ in TR (Figure 4).

However, the relative change in SO₂ concentrations in TR was particularly high (256%) from 1980 to 2021 (Table 2), likely due to low SO₂ concentrations in this remote region. Notably, the average of the increasing trends of MERRA-2 SO₂ over the four urban regions was approximately 16 times higher than the average SO₂ trends over the four rural regions, which is attributed to rapid economic growth, increase in energy consumption, and anthropogenic emissions in urban compared to rural regions [12,26,35]. For a more precise investigation of the long-term temporal changes in SO₂ concentrations, sliding trend analyses were applied to calculate the annual trends over urban and rural regions during different periods from 1980 to 2011 using spans of 10 years. As shown in Figure 5, the MERRA-2 SO₂ concentrations showed different trends throughout the study period. Therefore, the long-term variations and trends in SO₂ can be divided into three sub-periods: 1980–1997, 1997–2001, and 2001–2010.

From 1980 to 1997, MERRA-2 SO₂ concentrations increased over the four selected urban regions (Table 2) at rates of 0.57 µg m⁻³ yr⁻¹ in YRD, 0.44 µg m⁻³ yr⁻¹ in BTH, 0.37 µg m⁻³ yr⁻¹ in PRD, and 0.22 µg m⁻³ yr⁻¹ in SCB. Furthermore, it is worth noting that the SO₂ levels over the unpolluted and economically underdeveloped regions were also raised during 1980–1997, as SO₂ concentrations significantly increased over the four selected rural regions of China at rates of 0.01, 0.03, 0.007, and 0.008 µg m⁻³ yr⁻¹ over NER, MR, WR, and TR regions, respectively (Table 2). Previous studies showed an annual growth rate of 3.0–6.8% of SO₂ emissions over China during the periods before 1995 [12,13,26,35,84]. The overall increase in SO₂ over China during 1980–1997 (Figure 6) could be mainly attributed to the rapid growth of fossil fuel and coal consumption [35] due to increasing energy demands, especially in Eastern China (Figure 6). In addition, the percentage of vehicular emissions increased due to population and GDP growth until 1995 [26], which also enhanced SO₂ emissions, as shown below.

During 1997–2001, MERRA-2 SO₂ showed significant negative trends over eastern China (Figure 6). During this period, MERRA-2 SO₂ concentrations decreased by 25%, 28%, 11%, and 18% in YRD, BTH, PRD, and SCB, respectively. Over the rural regions, SO₂ concentrations also showed significant decreasing trends in NER (−0.1 µg m⁻³ yr⁻¹), WR (−0.01 µg m⁻³ yr⁻¹), and TR (−0.03 µg m⁻³ yr⁻¹). However, SO₂ concentrations presented a slight increasing trend in the TR at a rate of 0.03 µg m⁻³ yr⁻¹. The decreasing SO₂ concentrations over a large part of China during 1997–2001 could be mainly attributed to the Asian economic crisis, which started in 1997 with a decrease in GDP [85] (see also the following section), and then a decrease in energy consumption. Similar decreasing trends in industrial emissions and SO₂ concentrations were also observed over southern Europe, mainly in Greece, after the 2008 crisis and the subsequent economic recession [86]. In contrast, after 2001, the SO₂ concentrations increased significantly in all regions of China (Figure 6). Previous studies also noted that the SO₂ emissions increased by approximately 5–7% annually over China during 2001–2006 [87,88]. Our trend analysis showed that MERRA-2 SO₂ concentrations increased rapidly over the four selected urban regions at rates of 2.42, 1.68, 1.52 and 1.23 µg m⁻³ yr⁻¹ in YRD, BTH, PRD, and SCB, respectively. Moreover, SO₂ concentrations also showed significantly increasing trends over the selected rural regions at rates of 0.21, 0.16, 0.04, and 0.014 µg m⁻³ yr⁻¹ in NER, MR, WR, and TR, respectively.

The relative changes showed that the average SO₂ concentration trends were approximately 16 times higher over the four urban regions than in the corresponding rural regions during 2001–2010. This dramatic increase in SO₂ concentrations over the urban regions of China is mainly associated with rapid urbanization, increase in China’s economy, and large fuel consumption to meet increasing energy demands [26].

In contrast, the MERRA-2 SO₂ concentrations decreased slightly over China during the last decade (Figure 5), which is also shown for the period 2007–2021 (Figure 6). Several studies have shown that SO₂ emissions began to decline significantly after 2006, thereby improving air quality [26,28]. To better understand the changes in SO₂ emissions, CAMS and OMI data were also used, which were in general agreement with the MERRA-2 SO₂ variations, indicating downward trends in most parts of Eastern China (Figure 6). More specifically, the total column of SO₂ from CAMS data showed significant negative trends over the examined urban area, with rates of −0.15, −0.05, and −0.06 mg m⁻² yr⁻¹ in YRD, PRD, and SCB, respectively. SO₂ columns from the OMI sensors also showed significant negative trends over large areas close to the urban agglomerations of China with rates of −0.024, −0.03, and −0.02 DU yr⁻¹ in YRD, BTH, and SCB, respectively. Similar to a previous study of the impact of the Chinese Five-Year Plans (FYPs) on anthropogenic aerosols over China [28,89], our study suggests decreasing trends in SO₂ concentrations over most parts of eastern China in response to FYPs. Further details regarding the impact of FYPs on SO₂ concentrations are provided in Section 3.4. Over the four selected rural regions of China, MERRA-2, CAMS, and OMI showed insignificant trends during 2007–2021, except for NER, where CAMS showed a significant decreasing trend at a rate of 0.04 mg m⁻² yr⁻¹. However, there were some differences in the trends between MERRA-2, CAMS, and OMI data in some parts of China during the years 2007–2021, which were mainly attributed to the different databases, retrieval algorithms, and reanalysis data. In [17], the authors also noted differences in trend directions of SO₂ concentrations over India using MERRA-2, CAMS, and OMI data. Such differences could be attributed to missing satellite data and different assimilation algorithms used in the reanalysis data [17,90,91].

3.3. The Impact of Driving Factors on SO₂

To understand the influences of TP, GDP, and UP on SO₂ concentrations in each region, five years of relative changes in SO₂, UP, TP, and GDP in the eight ROIs during 1990–2020 were investigated (Figure 5). Relative changes in MERRA-2 SO₂ during 1990–1995 showed an obvious increase in SO₂ concentrations in BTH, PRD, MR, and TR, while SO₂ concentrations remained approximately stable over YRD and SCB and decreased in NER and WR (Figure 7). During the same period, UPs largely increased over the four selected urban regions compared to the rural regions (Figure 7). The TP data suggested an increase from 1990 to 1995 over all regions, except for YRD. GDP increased by approximately 50% and 40% in PRD and BTH, respectively, in the selected period and increased slightly to moderately over the other regions. Overall, the relationship between the changes in SO₂ and human activity factors could not be easily determined. For example, the most populated region of China, YRD, exhibited a negligible change in SO₂ concentrations during 1990–1995, corresponding to an approximately 25% relative increase in urban percentage, an approximately 10% increase in GDP, and a slight decrease in total population.

During the years 1995–2000, MERRA-2 SO₂ concentrations showed a notable decrease over all the selected regions except for TR, which was highly linked with an average decrease of 7–8% in Chinese GDP. In contrast, the UP and TP data showed increases in the selected regions during 1995–2000. Furthermore, during 2000–2005 and 2005–2010, SO₂ concentrations increased over all ROIs, which was in good agreement with the increase in UP and GDP in the eight ROIs during the same periods. However, TP data showed an average 5% decrease in relative changes in BTH and PRD during 2000–2005, and a 7% decrease in relative changes in the SCB during 2005–2010. During the last two periods, 2010–2015 and 2015–2020, the MERRA-2 SO₂ concentrations exhibited negligible changes or even slight decreases. However, these changes in SO₂ were not in good agreement with those of human factors because UP, TP, and GDP showed increasing tendencies over the different ROIs during 2010–2015 and 2015–2020, thus highlighting the effect of the clean-air policies followed in China during the last decade.

To better understand the changes in SO₂ concentrations with respect to human activity factors, a regression model was applied to China and the eight ROIs; the results are presented in Figure 8. A significant (p < 0.05) positive correlation between UP and SO₂ was observed in Central and Eastern China (Figure 8), which are regions highly dominated by urban areas. However, nonsignificant positive correlations between UP and SO₂ were found in Northern, Western, and Southern China, which are mostly rural areas. Our results showed that the UP in the four selected urban regions (YRD, BTH, PRD, and SCB) was positively correlated (p < 0.05) with SO₂, suggesting that the UP could be a major factor controlling SO₂ variations. Previous studies have suggested that the increase in PM₂.₅ concentrations over the urban regions of China is directly linked to the increase and expansion of urban areas [28,92]. Another factor that could affect SO₂ variation in China is the total population. Our statistical results suggest a significant positive correlation between TP and SO₂ over some urban areas in China with large populations. Therefore, TP in the three selected urban regions (YRD, BTH, and PRD) was positively (p < 0.05) correlated with SO₂, indicating that TP could be another factor controlling SO₂ variations in the selected regions. In the SCB, TP was not significantly correlated with SO₂, which could be attributed to the decrease in TP during a period of high increase in SO₂ emissions, i.e., from 2005 to 2010. Our results showed nonsignificant (p < 0.05) positive correlations between SO₂ concentrations and TP in rural regions (NER, MR, WR, and TR). This leads to the conclusion that the UP and TP may not play a significant role in SO₂ variation in the rural regions of China. In addition, GDP was positively correlated with SO₂ over vast areas of China. However, a significant correlation was found only in the MR. This finding was consistent with previous studies that suggested positive correlations between GDP, internal migration, and air pollutants in China [92,93].

The influence of the human factors was further tested using the GTWR model. The three variables (UP, TP, and GDP) were assimilated into the GTWR model under the conditions of Akaike Information Criterion (AIC) minimization and non-collinearity. Figure 8e shows the determination coefficients (R²) for provincial areas obtained from the GTWR model. The spatial distribution of the R² values revealed the combined effects of UP, TP, and GDP on SO₂ concentrations. The R² varied from 0.36 to 0.75 across China, with the highest values in the urban/populated provinces in eastern China, compared to the rural and remote regions in NW, NE, central China, and the TR (Figure 8e). This reveals that the UP, TP, and GDP seem to have a stronger influence on SO₂ variations in urban regions than in rural ones in China.

The contribution of each factor (UP, TP, and GDP) on SO₂ concentration changes was also investigated. Our analysis showed that the UP is a dominant human activity-driving factor that contributes to SO₂ concentration changes of about 23.3%, 22.2%, 25.4%, and 19.3% in YRD, BTH, PRD, and SCB, respectively (Figure 8f). Previous studies suggested that the urban class was a dominant factor in PM₂.₅ changes in east and south China [28,61]. Over YRD, BTH, PRD, and SCB explains about 19.1%, 17.3%, 20.2%, and 7.9% of the changes of SO₂ concentrations, respectively (Figure 8f). Furthermore, over the urban regions, the GDP exhibited low and insignificant contributions (3–5%) to SO₂ changes, while over the rural regions, the UP and TP contributed less (0.5–6%) to SO₂ variations (Figure 8f). GDP data contributed about 13% to SO₂ changes over the MR region. Previous studies have suggested that human activity factors, including urbanization, have nonsignificantly contributed to PM₂.₅ changes in the rural regions of China [28,61]. On the other hand, meteorological variables are believed to have significant influence on PM₂.₅ decadal changes [76,94], thus weakening the role of the population density, expansion of urban areas, and GDP [61,95]. It is worth mentioning that despite the continuous increase in population, GDP, and urban areas that potentially enhance air pollution, the SO₂ concentrations displayed a declining trend after the year 2010 over large regions in China. The decrease in SO₂ levels could be attributed to the effects of Chinese policies to control air pollution and other clean-air mitigation strategies, which are addressed below.

3.4. The Role of Chinese Government Policies

A series of air pollution prevention and control policies has been implemented in China over the last two decades. The 11th (2006–2010), 12th (2011–2015) and 13th (2016–2020) FYPs have been applied [28,89]. These plans were intended to improve the air quality in Chinese megacities and surrounding areas by reducing anthropogenic emissions.

To better evaluate the temporal changes in SO₂ in response to Chinese emissions control policies, a trend analysis was applied over the selected urban and rural regions using SO₂ concentration data from MERRA-2, CAMS, and OMI during the periods 2006–2010, 2011–2015, and 2016–2021 (

Table 3. Trends based on MERRA-2 SO₂ (µg m⁻³ yr⁻¹), CAMS SO2 (mg m⁻² yr⁻¹), and OMI SO₂ (DU yr⁻¹) in the eight selected regions during different periods. Bold text refers to significant trends (p < 0.05).

		Urban Regions				Rural Regions
Regions/Periods		YRD	BTH	PRD	SCB	NER	MR	WR	TR
		Trend	Trend	Trend	Trend	Trend	Trend	Trend	Trend
11th FYP (2006–2010)	MERRA-2	2.05	1.06	1.16	1.02	0.23	0.025	0.025	0.015
	CAMS	−0.45	−0.56	−0.16	−0.02	−0.001	−0.06	−0.02	−0.01
	OMI	−0.03	−0.08	−0.004	−0.017	−0.003	−0.023	−0.004	0.001
12th FYP (2011–2015)	MERRA-2	0.224	−0.37	−0.3	−0.28	−0.02	−0.04	−0.04	−0.006
	CAMS	−0.24	0.17	−0.15	−0.26	0.08	0.001	−0.05	−0.001
	OMI	−0.02	−0.05	0.002	−0.03	0.008	0.002	−0.002	0.001
13th FYP (2016–2021)	MERRA-2	−0.95	−0.53	−0.65	−0.56	−0.07	0.045	−0.017	−0.015
	CAMS	−0.45	−0.26	−0.014	−0.24	−0.04	−0.04	−0.004	−0.008
	OMI	−0.03	−0.02	−0.007	−0.01	0.01	−0.01	0.002	−0.007
(2013–2021)	MERRA-2	−0.189	−0.32	−0.18	−0.18	−0.07	0.03	−0.009	−0.005
	CAMS	0.01	−0.02	0.04	−0.13	−0.08	0.004	−0.02	0.003
	OMI	−0.023	−0.014	0.003	−0.017	−0.001	−0.005	−0.002	−0.003

). In general, the 11th FYP was planned to reduce the SO₂ concentrations by 10% during 2006–2010 [28,89]. Throughout the 11th FYP, MERRA-2 SO₂ showed nonsignificant positive trends over both the examined urban and rural regions, except for NER where the trend (0.23 µg m⁻³ yr⁻¹) was statistically significant. A previous study also showed that the AOD and PM₂.₅ concentrations from MERRA-2 data did not sufficiently decrease in response to the 11th FYP [28,96].

In comparison, CAMS and OMI SO₂ showed negative trends in the four selected urban areas with significant declining trends of −0.45 and −0.65 mg m⁻². yr⁻¹ in YRD and BTH, respectively, and −0.08 DU yr⁻¹ in BTH, thus coinciding with the implementation of the plans. In response to the 11th FYP, CAMS and OMI SO₂ levels revealed a reduction range of 2–45% over the four examined urban regions during 2006–2010. SO₂ concentrations began to decrease in 2006, which was mainly attributed to the application of flue gas desulfurization in Chinese power plants [97,98]. In addition, the rural regions exhibited negative trends in CAMS and OMI SO₂ data, with significant declining rates of −0.06 mg m⁻² yr⁻¹ and −0.023 DU yr⁻¹ for CAMS and OMI in MR, respectively (Figure 9a). These reductions in rural regions could be mainly associated with a general reduction in urban/industrial SO₂ emissions from the main sources in East China.

During the 12th FYP, the ratio of energy consumption to per-capita GDP was planned to decrease by 16% and SO₂ concentrations were to be reduced by 8–10%. During the 12th FYP implementation period (2011–2015), MERRA-2 SO₂ concentrations began to exhibit negative trends over urban regions, except for the YRD. In the four rural regions, negative MERRA-2 SO₂ trends were also observed (Figure 9b). CAMS SO₂ also showed decreasing trends in YRD and PRD, whereas BTH exhibited an increasing trend. Furthermore, over the urban regions of YRD and BTH, OMI SO₂ concentrations showed a decreasing trend. Nonetheless, SO₂ increased slightly in the PRD during the same period. In 2011–2015 (12th FYP), CAMS and OMI SO₂ concentrations showed significant decreasing trends in SCB with rates of −0.26 mg m⁻² yr⁻¹ and −0.03 DU yr⁻¹, respectively. SO₂ reductions were also noted over the rural regions of China using CAMS and OMI data, which could be mainly attributed to SO₂ reduction from sources in Chinese urban regions.

A 15% reduction in SO₂ by 2020 compared to 2015 levels was planned through the 13th FYP. During 2016–2021, the three datasets showed consistent negative SO₂ trends over the urban regions (Figure 9c), with significant decreasing trends in SCB for MERRA-2 (−0.95 µg m⁻³ yr⁻¹) and OMI (−0.03 DU yr⁻¹). Overall, MERRA-2, CAMS, and OMI SO₂ data revealed a reduction range of 3–30% from 2016 to 2021 over the four urban regions. In the four rural regions, negative trends were also observed during the same period in the three datasets, except for in NER and WR, as OMI showed increasing trends of 0.01 and 0.002 DU yr⁻¹, respectively.

Because there were no ground measurements in 2012, 2013 was used as the primary year to evaluate the achievement of the 12th and 13th FYPs. Furthermore, the 12th and 13th FYPs largely focused on reducing SO₂ levels via desulfurization in the industrial and traffic sectors [28,96,99]. Thus, a trend examination was applied to SO₂ from the CAMS, OMI, MERRA-2, and ground stations in the years 2013–2021 (Figure 9). The results indicate that SO₂ concentrations from MERRA-2 declined over eastern China and in the four urban regions. However, the distribution of the trends was more variable for the CAMS data, which showed increasing trends in the YRD and PRD regions and a significant decreasing trend in the SCB at a rate of −0.13 mg m⁻² yr⁻¹. In contrast, OMI also exhibited a large spatial variation in trend values, with significant decreasing trends in the three urban regions, except for the PRD (Figure 9g). Similar significant decreasing trends in SO₂ concentrations were also observed at the ground monitoring stations during 2013–2021, which were mostly detected in Eastern China but also in the northeast, central, and northwest regions, as well as in the TR and southern rural areas. Therefore, the implementation of FYPs since 2007 has led to a significant drop in SO₂ concentrations during the last decade in China, not only in urban/industrialized regions but also in rural and remote areas. Consistent with several previous studies [26,28,95], the FYPs reduced SO₂ concentrations (among other air pollutants) and improved air quality.

4. Conclusions

In this study, we assessed the long-term trends in SO₂ concentrations over China from 1980 to 2021 based on reanalysis data (MERRA-2). We also examined the SO₂ total column using data from CAMS reanalysis and OMI satellite observations from 2007 to 2021. Eight representative regions, including four urban (PRD, BTH, YRD, and SCB) and four rural (NER, MR, WR, and TR), were selected for analysis. The relationships between SO₂ multi-year variations and influencing factors, including GDP, population, and urbanization, were explored using linear regression and geographical temporally weighted regression models to evaluate the anthropogenic impacts on SO₂ concentrations. Finally, we examined the influence of control policies on SO₂ trends after 2006. The main findings of this study are summarized as follows:

• The multi-year average SO₂ concentrations from MERRA-2 data from 1980 to 2021 over the four rural regions are approximately 11 times lower than those over the four urban regions.
• The SO₂ concentration increased significantly over the selected regions from 1980 to 2021, with major changes occurring in between. It increased during 1980–1997 and 2001–2010, but dropped during 1997–2001 and 2010–2021. The relative change showed that the average MERRA-2 SO₂ concentration trends over the four urban regions were approximately 16 times higher than those in the four rural regions.
• The results revealed that the driving factors associated with human activities, including urbanization, GDP, and population, played significant roles in multi-decadal SO₂ variations and trends in the main urban areas of China.
• The SO₂ concentration significantly decreased in most regions of China after 2010, which was attributed to the control policies of China (12th and 13th FYP).
• The OMI SO₂ data showed significant downward trends in the last decade over most regions of China, exhibiting better agreement with SO₂ variations from ground-based stations than with SO₂ data from the MERRA-2 and CAMS reanalysis.

The assessment of the long-term dynamics of SO₂ concentrations revealed that variations in SO₂ concentrations were highly influenced by human activities in the urban areas of China. Nonetheless, with continued anthropogenic activities associated with urban expansion, SO₂ levels showed a declining trend after 2007, which is credited to the control policies of China. Our project improves our understanding of the long-term changes in SO₂ concentrations in urban regions of China compared to rural regions, as well as the effects of mitigation strategies and air pollution control plans.

5. Limitations and Future Prospects

This study used different SO₂ datasets, including observation data from OMI, and reanalysis data from MERRA-2 and CAMS. For the OMI satellite data, the primary sources of uncertainty were mainly associated with haze and cloud condensation. In addition, some gaps were presented, which may cause some uncertainties in SO₂ retrievals from satellite. On the other hand, MERRA-2 and CAMS assimilate different sources of AOD in Earth system modeling. The differences across the datasets are probably attributable to the algorithms used. Different algorithms may produce some uncertainties. Regardless of the limitation, this study would be of major interest, especially over the Chinese regions, and it comprises important information that could be of use to climate change, air quality, human activities and health in China [100,101]. Therefore, our findings highlight the importance of future studies on the impact of SO₂ on regional climate change in China.