Investigation of Spatiotemporal Variation and Drivers of Aerosol Optical Depth in China from 2010 to 2020

Abstract

China has experienced rapid economic growth and serious control of aerosol emissions in the past decade. Thus, the spatiotemporal variations and driving factors of aerosol optical depth (AOD) are urgently needed to evaluate the effectiveness of aerosol control activities. The innovation of this study is a detailed spatial and temporal analysis of aerosol pollution in eight major regions of China from 2010 to 2020 using the MERRA-2 AOD reanalysis product and the driving mechanism based on the Granger causality test, sensitivity, and contribution analysis. The results show that the spatial distribution of AOD varied across the areas. Divided by the Hu Line, the AOD values of the Eastern areas were significantly higher than those of the Western areas. The temporal trend in the last eleven years was dominated by a continuous decline and moderate fluctuations at both annual and seasonal scales. The relationship between socioeconomic factors and AOD drivers was more significant in economically developed regions, suggesting that China pays more attention to haze control while developing its economy. The driving relationship between AOD and temperature was weak, while wind speed and relative humidity were more influential. For vegetation factors, Granger effects were mainly observed in the Northeast, Beijing-Tianjin-Hebei, Guangdong, Central China, and Southwest regions. In the Guangdong and Southwest regions, vegetation and economic factors were the more influential drivers. This study provides a scientific basis for the detection of aerosol changes, driving mechanisms and pollution management in China.

1. Introduction

Aerosol refers to a variety of small solid and liquid particles suspended in the atmosphere with sizes ranging from 0.1–100 nm [1,2,3]. By absorbing and scattering solar radiation, aerosols affect the radiation budget between Earth and the atmosphere and consequently influence terrestrial carbon transfer [4,5,6]. Moreover, aerosols, as cloud condensation nuclei, directly influence the formation and microphysical characteristics of clouds, such as cloud albedo, which can further induce changes in weather and climate [2,4,7,8]. Heavy aerosol loadings can cause air pollution, reduce visibility, and increase the risk of respiratory diseases in humans [2,9,10,11,12,13].

The aerosol optical depth (AOD) is computed as the integral of the solar extinction coefficient in the vertical direction based on the weakening of incident solar light by aerosols. Therefore, AOD is a fundamental property indicating the concentration of aerosols and a key physical quantity characterizing the degree of atmospheric turbidity [14,15,16]. Recent developments in aerosol-related research have shown an increasing need for AOD data with high spatial and high temporal resolutions [17,18]. Traditional methods to retrieve AOD include ground measurements and satellite remote sensing. The former is limited by the sparse distribution of aerosol observation networks such as AERONET [19,20], which makes it difficult to obtain aerosol information on a large scale. The latter is limited by the availability and quality of the acquired satellite data [21,22], such as the influence of clouds. In contrast, reanalysis AOD products are particularly advantageous for providing spatially and temporally consistent AOD information derived from data assimilation systems that integrate global climate models and multisource AOD measurements [17]. Previous studies have successfully used reanalysis AOD products to investigate the spatiotemporal changes in AOD [23,24,25], such as the Modern-Era Retrospective Analysis for Research and Applications Aerosol (MERRA) reanalysis product [25], yet the accuracy of the reanalysis product needs to be carefully evaluated.

As a direct indicator of air pollution, the dynamics and drivers of AOD have been the focus of research [26,27]. Numerous studies have attempted to investigate the spatial and temporal changes in AOD, as well as the driving forces, in hazy areas using satellite-derived or reanalysis AOD products [28,29]. Existing studies have shown that human activities such as industrial production, transportation, straw burning, and energy consumption are important factors influencing AOD [30]. In addition, air quality is not only related to anthropogenic aerosol emissions but is also influenced by meteorological conditions. While pollution sources are relatively stable, meteorological factors, including temperature (TEM), relative humidity (RHU), and wind speed (WIN), have important effects on the evolution of aerosols [31]. In addition, vegetation can effectively absorb harmful particles in the atmosphere, such as sulfur dioxide. The inverse relationship between vegetation cover and dust emissions is undeniable, with varying levels of influence in different regions [32,33]. China is located across a wide latitude and under different environmental conditions, which leads to each parameter having different effect levels under different conditions [14,30,34].

The methods commonly used for driving force analysis mainly include correlation or regression analysis [6,35,36], the Granger causality test, and sensitivity analysis [37,38]. However, most regression-based analyses actually reveal a statistical relationship rather than a causal relationship between AOD and related factors [39] since they cannot determine the direction of causality between two related variables [40,41]. Granger causality tests reveal causal relationships between variables [39] and identify the direction of causality between two related variables [40,41], but they do not provide a sufficiently deep quantitative description to measure the strength of causality. Sensitivity analysis quantitatively describes the relationship between the variables and the drivers [42], but it does not consider a causal relationship. Therefore, this paper combines the Granger causality analysis and sensitivity analysis methods together to investigate AOD driving mechanisms to enhance our understanding of AOD dynamics.

China has experienced rapid economic growth and serious air pollution in the past few decades [2,28]. Since the 18th National Congress of the Communist Party of China (CPC) convened in 2012, China has made the response to climate change a higher priority in state governance and has adopted a series of policies, measures, and actions to tackle climate change, and the mitigation of air pollution is one of the specific goals [43]. To verify the effectiveness of air pollution control and provide guidance for further actions, this paper aims to outline a detailed investigation of the spatiotemporal variations and driving mechanisms of the AOD in mainland China from 2010 to 2020. The innovations of this study include (1) using the MERRA-2 AOD reanalysis product with a 3 h temporal resolution to describe highly dynamic variations in AOD and a detailed evaluation of the trend, robustness, and variability of the changes in AOD, and (2) detailed analysis of the socioeconomic, meteorological, and vegetation drivers combining the Granger causality test [44], sensitivity, and contribution rates. This study provides an up-to-date analysis of the long-term trends in aerosols and their driving factors in major areas in China, which will provide a reference for related studies on environmental changes and for the ongoing clean air acts in China.

2. Materials and Methods

2.1. Study Area

To adapt to the differences in aerosols across space, regional analysis was performed to identify long-term trends in representative areas in China. The China Atmospheric Environment Bulletin [45] identifies eight key regions that are close in economic level and form industrial zones of a certain scale. At the same time, these eight regions have different degrees of aerosol pollution and are the key areas for the Chinese government to combat haze. This includes Beijing-Tianjin-Hebei, Fenway Plain, Yangtze River Delta, Guangdong Province, Northeast China, Central China, Southwest China, and Northwest China, which we will refer to as BTH, FWP, YRD, GD, NE, CC, SW, and NW hereafter, as shown in Figure 1 and

Table 1. Administrative districts of eight major areas.

Region	Administrative District
Beijing-Tianjin-Hebei (BTH)	Beijing city, Tianjin city, Hebei Province
Fenway Plain (FWP)	Shanxi Province, Shaanxi Province and Henan Provinces
Yangtze River Delta (YRD)	Shanghai city, Jiangsu Province, Zhejiang Province, Anhui Province
Guangdong province (GD)	Guangdong Province
Northeast (NE)	Heilongjiang Province, Jilin Province, Liaoning Province
Central China (CC)	Hunan Province, Hubei Province, Jiangxi Province
Southwest (SW)	Yunnan Province, Sichuan Province, Guizhou Province, Chongqing city
Northwest (NW)	Shaanxi Province, Gansu Province, Ningxia Province, Xinjiang Province

. In addition, this paper divided mainland China into two parts using the Hu Line, which divided the country into two parts with significantly different population densities and socioeconomic development [46]. Changes in AOD in different areas were analyzed using zonal analysis.

2.2. Data Description and Processing

2.2.1. MERRA-2 AOD Data

The MERRA-2 reanalysis product (available at https://disc.gsfc.nasa.gov, accessed on 15 May 2022) is derived from a data assimilation system that integrates the GEOS-5 climate model and a variety of satellite- and ground-observed AOD data, such as the satellite-derived AOD data from AVHRR, MISR, and MODIS sensors and ground measurements from AERONET. It provides 21 categories of data from 1980 to the present, including water vapor, radiation, AOD, etc. [25]. MERRA-2 AOD data with a 0.5° × 0.625° spatial resolution and a 3 h temporal resolution from 2010 to 2020 in mainland China were used in this study. Daily averages of MERRA-2 AOD data at 3 h intervals were calculated for subsequent analysis.

2.2.2. AREONET AOD Data

AERONET (available at http://aeronet.gsfc.nasa.gov, accessed on 18 May 2022)) is a global ground-based network that provides long-term measurements of aerosol optical properties at wide-range wavelengths from 340 to 1640 nm with a high temporal resolution [47]. AERONET uses the CE-318 sun photometer to measure and provide data on aerosol optical, microphysical, and radiometric properties for more than 25 years [7,17,26]. AERONET aerosol data have three different levels, including level 1.0 unprocessed raw data, level 1.5 cloud filtered and quality controlled data, and level 2.0 cloud screened and quality assured data, corresponding to different sites [28]. A total of 87 AERONET sites in China (Figure 1), measuring instantaneous AOD at a temporal resolution of less than one hour, were used in this study. Since very few sites can provide complete records of AOD data for 2010–2020, we included all AERONET data of levels 1.5 and 2.0 at the moments near the time of MERRA-2 AOD data at 3 h intervals for validation.

To maintain consistency with the MERRA-2 AOD data at 550 nm, we first interpolated the AERONET AOD measurements at 500 nm and 675 nm to 550 nm using:

$$AO{D}_{550}=\beta \times {550}^{-\alpha }$$

(1)

$$\left\{\begin{array}{c}\alpha =-\frac{\ln (\frac{AO{D}_{500}}{AO{D}_{675}})}{\ln (\frac{500}{675})}\\ \beta =\frac{AO{D}_{500}}{{500}^{-\alpha }}=\frac{AO{D}_{675}}{{675}^{-\alpha }}\end{array}\right.$$

(2)

where $\alpha$ is the Angstrom exponent, estimated from the slope of the spectral AOD plots in logarithmic scales, $\beta$ is the turbidity coefficient, AOD₅₀₀ means AOD at 550 nm, and AOD₆₇₅ means AOD at 675 nm.

2.2.3. Socioeconomic Data

As both aerosol emissions and controls are related to socioeconomic development, this paper selected two socioeconomic indicators, gross domestic product per capita (GDPPC) and industrial output per square kilometer (IOPSK), to investigate the feedback between AOD changes and socioeconomic development. Data were collected from the China Statistical Yearbook (http://www.stats.gov.cn/, accessed on 20 May 2022) from 2010 to 2020, which provides academic statistics issued by the National Bureau of Statistics. However, the China Statistical Yearbook has not publicized provincial industrial output data since 2020, thus the IOPSK data from 2010–2019 were included in our analysis.

2.2.4. Meteorological Data

This study collected daily meteorological data from the 699 primary meteorological stations from the China Meteorological Administration website (https://cmdp.ncc-cma.net, accessed on 22 May 2022). Wind speed (WIN) affects the velocity of airflow, temperature (TEM) indicates the intensity of particle motion, and relative humidity (RHU) indicates the percentage of water vapor pressure and saturated water vapor pressure in the air, with variables that have a significant impact on AOD evolution. Therefore, we chose WIN, TEM, and RHU as the weather driving factor, as claimed in previous studies [30]. All daily meteorological data were temporally aggregated to different time scales to adapt to the analysis of the AOD spatiotemporal variations. Figure 2 shows the distribution of the meteorological sites and the seasonal mean values of the three factors. In our analysis, spring refers to the months from March to May, summer from June to August, autumn from September to November, and winter from December to February. All the meteorological data were spatially interpolated to raster images using the Kriging interpolation method [48,49,50] and then resampled to the raster cells of MERRA-2 AOD images for comparative analysis.

2.2.5. Vegetation Continuous Fields Data

The Terra MODIS Vegetation Continuous Fields (VCF) product (MOD44B) is a subpixel-level representation of surface vegetation cover estimates globally. It is designed to continuously represent Earth’s terrestrial surface as a proportion of basic vegetation traits and provides a gradation of three surface cover components: percent tree cover, percent nontree cover, and percent bare. Considering the large-scale reforestation activities in China, this study used the yearly percent tree cover data of MOD44B with 250 m spatial resolution in China from 2010 to 2020 to evaluate the effect of increasing trees on AOD changes.

2.3. Methodology

Based on MERRA2-AOD data, this paper investigates the spatial and temporal evolution of AOD in China from 2010 to 2020 at the pixel and regional scales and analyzes the driving mechanisms of meteorological, economic, and vegetation factors on AOD. The overall technical process is shown in Figure 3.

2.3.1. Mann–Kendall (M–K) Test

The Mann–Kendall (M–K) test was used to analyze the temporal trend in the seasonal and annual means of AOD at each pixel. The M–K test is a nonparametric statistical test that detects long-term trends as well as sudden changes in a meteorological element [51]. The mathematical formulations of the M–K test are as follows.

First, assume a set of time series data m₁, …, m_n in which each variable is independent and randomly distributed. If k, j ≤ n, and k ≠ j, the statistics of the test are calculated as follows:

$$S={\displaystyle \sum }_{k=1}^{n-1}{\displaystyle \sum }_{j=k+1}^{n}Sgn\left(m_j-m_k\right)$$

(3)

Among them,

$$Sgn\left(m_j-m_k\right)=\left\{\begin{array}{c}1,\left(m_j-m_k\right)>0\\ 0,\left(m_j-m_k\right)=0\\ -1,\left(m_j-m_k\right)<0\end{array}\right.$$

(4)

S follows a normal distribution with a mean of 0 and a variance of Var(S) = n(n − 1)(2n + 5)/18. The original M–K statistic, designated by Z, was computed as follows:

$$Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{Var\left(S\right)}},S>0\\ 0,S=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}},S<0\end{array}\right.$$

(5)

where Z > 0 indicates a growing trend and Z < 0 indicates a decreasing trend. When |Z| > 1.96 and 2.58, the trend identified by the Z value is statistically significant with confidence levels of 95% and 99%, respectively [52].

2.3.2. Trend Sustainability

Based on the M–K test, a rescaled polar difference (R/S) analysis-based approach was used to calculate the Hurst index to characterize the sustainability of the AOD change trend. Details can be found in Li, et al. [53]. A Hurst index greater than 0.5 means that the future trend is consistent with the past trend, and a Hurst index less than 0.5 means that the future trend is opposite to the past trend.

2.3.3. Coefficient of Variation

The coefficient of variation (CV) is used to quantify the variation in the AOD at annual scales, expressed as:

$${\mathrm{CV}}_{\mathrm{AOD}}=\frac{{\sigma }_{AOD}}{\overline{AOD}}$$

(6)

where CV_AOD represents the coefficient of variation of the AOD from 2010 to 2020, ${\sigma }_{AOD}$ is the standard deviation of annual means of AOD, and $\overline{AOD}$ is the multiyear average of the AOD values from 2010 to 2020. Larger CV_AOD values indicate more significant fluctuations.

2.3.4. Granger Causality Test

The Granger causality test was the first attempt to quantify the causal relation between time series, in which causality and feedback are defined in an explicit and testable manner [40,44]. The Granger causality test method uses different lag orders for two time series and then evaluates the impact of one time series on the other time series. The regression between the two series can be expressed as:

$$y_t={\displaystyle \sum }_{i=1}^l{\gamma }_i{x}_{t-i}+{\displaystyle \sum }_{i=1}^l{\tau }_i{y}_{t-i}+{\epsilon }_t$$

(7)

where $\gamma$ and $\tau$ are regression coefficients; $\epsilon$ is random error; and l is the lag length (=1 in this study). In Equation (7), Y is a function of previous values of the dependent variables (i.e., Y) and independent variables (i.e., X). To test whether X causes Y, we eliminate X by restricting $\gamma$ in Equation (7) to zero and derive the regression function as:

$$y_t={\displaystyle \sum }_{i=1}^l{\tau }_i^{\prime }{y}_{t-i}+{\epsilon }_t^{\prime }$$

(8)

To test whether the restricted estimates are statistically significantly different from the unrestricted estimates [i.e., Equation (7) or Equation (8)], the following F ratio is calculated:

$$F_{X~Y}=\frac{\left(R{S}_R-R{S}_U\right)/H}{R{S}_U/\left(p-R\right)}~F\left(H,p-R\right)$$

(9)

where RS_R and RS_U are the sum of squared errors of restricted and unrestricted versions of Equations (7) and (8), respectively. The term H is the number of coefficients set to zero in the restricted version, R is the number of predictors in the unrestricted version, and p is the number of observations. Taking the confidence level as 0.95, if F_X~Y ≤ F_0.05 (H, p-R), we accept the null hypothesis that the variable eliminated from Equation (7) is not the cause of the dependent variable according to the Granger causality test. Otherwise, we reject the null hypothesis.

The Granger causality between AOD and related factors was investigated on different temporal and spatial scales. The Granger causalities between AOD and socioeconomic factors (i.e., GDPPC and IOPSK) were analyzed by the variations in their regional means and annual means from 2010 to 2020. The Granger causality tests between AOD and meteorological factors were analyzed by the seasonal means at both the pixel and regional scales. For VCF, Granger causality tests were conducted using the annual means at both the pixel and regional scales. This can help reveal the operational scale of the factors driving AOD changes.

2.3.5. Sensitivity and Contribution Rates

The sensitivity rate (SR) of the AOD to the driving factors was defined as [37,38,42]:

$$\mathrm{SR}=\frac{\overline{DF}}{\overline{AOD}}\times \frac{\sum \left(D{F}_i-\overline{DF}\right)\left(AO{D}_i-\overline{AOD}\right)}{\sum {\left(D{F}_i-\overline{DF}\right)}^2}$$

(10)

where DF_i is the average value of the driving factors in year i, including economic, meteorological, and vegetation driving factors. $\overline{DF}$ and $\overline{AOD}$ are the mean values of the driving factors and AOD from 2010 to 2020, respectively. SR > 0 means AOD increases with increasing driving factors, while SR < 0 means AOD decreases with increasing driving factors. A larger |SR| indicates a higher sensitivity of the driving factors to AOD. The contribution rate (CR) of the driving factors to the AOD is defined as [42]:

$$\mathrm{CR}=\frac{L\times Se{n}_{slope}}{\overline{DF}}\times \mathrm{SR}\times 100\%$$

(11)

where Sen_slope represents the rate of change of the driving factor and L is the length of the time series. In this paper, Formulas (10) and (11) are used to calculate the sensitivity rate of AOD to the drivers and the contribution rate of the drivers to the AOD. CR > 0 indicates that the driving factors positively contributes to the AOD, i.e., an increase in driving factors leads to an increase in AOD, while CR < 0 indicates that the driving factors negatively contributes to the AOD, i.e., an increase in the driving factors leads to a decrease in AOD. A larger |CR| indicates a higher contribution.

3. Results

3.1. Data Validation

To validate the accuracy of the MERRA-2 product, the AOD values derived from AERONET measurements and the MERRA-2 product were compared at the site level, and the AERONET data that are closer to the times of MERRA-2 AOD were selected. Figure 4 compares AERONET AOD data of levels 1.5 and 2.0 and MERRA-2 AOD at close moments. A good consistency between MERRA-2 AOD and AERONET measurements was found, with a root mean square error (RMSE) of 0.16 and a correlation coefficient (R²) of 0.73. Overall, the data accuracy is good and consistent with previous studies [14,30]. A number of data points deviated from the 1:1 line, probably due to numerous reasons, such as uncertainty in measurement and models and rapidly changing atmospheres. The large difference in spatial resolutions between the site-level AERONET data and the coarse-resolution MERRA-2 AOD may be a major reason.

3.2. Seasonal and Annual Mean AOD

Figure 5 shows the seasonal and annual average AOD values in China for the period 2010 to 2020. Divided by the Hu Line, the regional mean AOD of the Eastern and Western areas were 0.446 ± 0.162 and 0.188 ± 0.115, respectively, indicating significant differences in the two areas. For the Eastern and Western regions, the variation in AOD was similar across seasons. The AOD is higher in spring and summer, and a high AOD area in the Western region is located in the Taklamakan Desert region, where high wind speed in spring and summer cause frequent windy and sandy weather in the Taklamakan Desert [25,54]. The lowest wind speed and slow air movement in winter resulted in the lowest AOD in the Western region in winter. The overall high AOD in spring in the Eastern region and the high AOD in summer are concentrated in the BTH, which is related to local human activities [2]. The standard deviation (STD) values of AOD followed the order of winter > summer > autumn > spring, which was different from the order of mean AOD values. The mean value reflects the overall level of AOD, while the standard deviation indicates the regional variations in AOD. This indicates that in the Eastern region, there are distinct geographical characteristics of AOD across seasons due to differences in aerosol sources, topography, and climate. AOD values are also higher in summer in the BTH region, which may be due to the high aerosol emissions from local industrial production, straw burning, traffic, and energy consumption. This is the direct cause of the high standard deviation in the summer in the Eastern region. In contrast, the overall AOD in spring is high, and the standard deviation is low. This indicates that the AOD values in winter and summer have large differences in different regions [30].

To account for the regional differences, Figure 6 shows the statistical distribution of seasonal and mean AOD values in different major areas. The annual mean AOD values ranged from 0.271 to 0.587 and followed the order of CC > YRD > GD > FWP > BTH > SW > NE > NW. The areas of the FWP, NW, and YRD showed very small seasonal variations, while the other areas showed larger seasonal variations. The AOD in GD is higher in spring than in other seasons because coarse dust particles transported from the north over long distances in spring adversely affect the air quality in GD [3]. The spatial variability of AOD in the SW and BTH regions is large. Seven out of eight major areas had the highest AOD values in spring, while only the BTH area had the highest AOD values in summer. Such differences might be attributed to different driving forces.

3.3. Temporal Trend of AOD

Figure 7 shows the temporal trend of daily average AOD values in the Eastern and Western areas of China divided by the Hu Line from 2010 to 2020. It clearly shows that the AOD has been decreasing from 2010 to 2020 in both the Eastern and Western areas of China, indicating effective control for aerosol emissions. Such a decreasing trend is more significant in the Eastern area than in the Western area. This is reasonable because the AOD in the Eastern areas was higher due to the greater intensity of human activity and might have experienced a larger reduction after control. Figure 8 shows the monthly average and standard deviations of AOD values across the 11 years from 2010 to 2020 in the Eastern and Western regions of the Hu Line. The highest AOD in the Western region was in May (0.261 ± 0.159), which was lower than the lowest AOD of the Eastern region, which occurred in December (0.359 ± 0.21). Overall, the highest AOD value occurred in April (0.3775 ± 0.2305), and the lowest value occurred in December (0.2024 ± 0.1908).

Figure 9a shows the temporal trends detected by the M–K test method from the annual means and seasonal means from 2010 to 2020. In the Eastern areas, the AOD is generally dominated by anthropogenic aerosols, and the effective control of aerosol emissions led to a significant decrease in these areas. In contrast, aerosols in the Western part of the Hu Line are mainly natural since these areas are less populated and developed. The increasing trends of AOD in the Xinjiang area in four seasons are related to special topography, aerosol transfer from other regions, and changes in local climatic factors [54]. For temporal variations in annual mean AOD values, 60.7% of the area in China, especially 34.4% of the Eastern area of the Hu Line, showed a decreasing trend with a 95% confidence level, confirming the significant reductions in AOD temporally. Across different seasons, the largest reductions in AOD were found in summer, as 54.6% of the area showed a decreasing trend. The AOD in spring was also reduced significantly and was mainly distributed in CC and NE. In autumn, only 29.02% of the area showed a decreasing trend, while a small portion of the area in Xinjiang even showed an increasing trend. The reductions in winter were the smallest, and only 20.53% of the area showed a decreasing trend. Figure 9b shows the spatial distribution of Hurst indices in China from 2010 to 2020. On the annual scale, 88.6% of the areas derived Hurst indices greater than 0.5. This indicates that the general decreasing trend in most areas and the local increasing trend in the Taklamakan Desert will likely continue. On seasonal scales, the current trend in most areas seemed to continue deriving Hurst indices > 0.5, while the exception is that the Taklamakan Desert region in autumn deriving the Hurst index < 0.5, indicating that the increasing trend is not sustainable.

Figure 9c shows the annual and seasonal average CV values of AOD in China from 2010 to 2020. On the annual average scale, the overall CV value was 0.11, which shows higher fluctuations in FWP, CC, BTH, and the border areas of NE and Inner Mongolia. Seasonally, the CV in NE and northern Inner Mongolia are higher in spring and summer, while those in CC, southern FWP, and western YRD are higher in summer, which indicates a large fluctuation of the time series. The possible reason is that straw burning in NE and CC affects the air quality in the surrounding areas. With the gradual improvement of air pollution control measures in China, aerosol pollution has decreased, thus producing large fluctuations in the time series. The fluctuations were most significant in the Taklamakan Desert region in autumn, probably due to the windy and sandy weather in the Taklamakan region. In winter, the highest fluctuation was found in Yunnan Province.

Table 2. Temporal trends in the AOD for the eight major regions detected by the M–K test and Hurst index.

	NE	FWP	CC	BTH	NW	SW	YRD	GD
Spring	−1.713	−3.425	−2.958	−2.335	−0.467	−2.491	−3.270	−1.090 *
Summer	−2.491	−3.425	−2.803	−3.270	−0.467	−3.114	−2.335	−2.335 *
Autumn	−2.335	−1.713	−2.647	−2.491	1.246 *	−2.647	−2.335	−2.491
Winter	−2.024	−2.180	−1.401	−1.713	−2.024	−1.401	−2.024	−2.803
Annual	−2.491	−3.425	−2.647	−3.581	−0.934	−3.114	−3.270	−2.803

Note: The color of each cell represents different trends. White indicates no significant change (|Z| < 1.96), light green indicates a decreasing trend (−2.58 < Z < −1.96), and light blue indicates a significant decreasing trend (Z < −2.58). * means that the Hurst index is less than 0.5.

shows the trends of AOD and sustainability in eight major areas. According to the Z values, seven out of eight major areas showed decreasing trends in AOD with a 95% confidence level. The NW area showed insignificant changes, possibly due to a lower intensity of human activities and a low level of aerosol emissions in the past decade. For each major area, the decreasing trends were different across different seasons. The most significant decreasing trends were found for the FWP area in spring and summer, CC area in spring, summer, and winter, BTH area in summer, SW area in summer and autumn, YRD area in spring, and GD area in winter. Combining the seasonal mean AOD in eight major areas in Figure 5, the seasons that had the most significant decreasing trends coincided with the hazy seasons in each major area. This indicates that aerosol control measures might have been taken effectively for the hazy seasons in each major area in China. According to the Hurst indices, the AOD trend in the GD region in spring and summer and the NW region in autumn will not be sustained. On the annual scale, the Hurst indices in all regions are greater than 0.5, which indicates that aerosol pollution will continue to decrease in major regions of China.

3.4. Analysis of Driving Factors

3.4.1. Socioeconomic Factors

Figure 10 shows the temporal changes in the annual values of AOD, IOPSK, and GDPPC in eight major areas from 2010 to 2020. Generally, the AOD showed decreasing trends, while the IOPSK and GDPPC values showed increasing trends in all eight major areas from 2010 to 2020. This indicates that China has effectively controlled aerosol emissions while maintaining upward economic growth. A positive relationship between the AOD and socioeconomic factors (i.e., IOPSK and GDPPC) was found in most areas, in which an area with a high level of economic development normally means a high level of AOD. For example, the developed YRD and GD areas were hazier, while the less developed NE, NW, and SW areas were less hazy. However, CC is an exception, which was less developed but had the highest AOD values from 2010 to 2020 among the eight areas. This also necessitates further analysis of the factors driving the spatiotemporal changes in AOD.

Table 3. Driving force analysis performed on a regional scale between AOD and socioeconomic factors of major regions.

	GDPPC			IOPSK
	GCT	SR	CR/%	GCT	SR	CR/%
NE	/	−0.083	−2.0	/	−0.202	−3.573
NW	/	0.012	0.979	/	−0.030	−3.394
SW	/	−0.269	−26.398	√	−0.223	−20.999
YRD	/	−0.431	−36.907	/	−0.439	−22.623
FWP	/	−0.591	−24.408	/	−0.325	−28.175
CC	/	−0.358	−34.187	/	−0.329	−22.270
BTH	√	−0.378	−26.156	/	−0.503	−33.358
GD	/	−0.436	−31.280	/	−0.316	−28.924

Note: X is GDPPC/IOPSK, Y is AOD. √ indicates that causality is established (p ≤ 0.05) and / indicates that causality is not established (p > 0.05).

shows the results of the driving mechanism between AOD and socioeconomic factors. Among the eight regions, the Granger causality between AOD and GDPPC was confirmed in BTH, while the causality between AOD and IOPSK was confirmed in SW. According to the results of SR, AOD responded most significantly to GDPPC in the FWP region and to IOPSK in the BTH. According to the CR results, the contribution of GDPPC to AOD was higher in the YRD, CC, and GD regions. The contribution of IOPSK to AOD was higher in the BTH, FWP, and GD regions. Overall, economic development contributes to AOD in most of the regions, which indicates that local economic development benefits haze control activities. From 2010 to 2020, while China’s industrial structure changed, the use of fossil fuels decreased, and new energy sources were gradually exploited, which are also important for mitigating aerosol pollution.

3.4.2. Meteorological Factors

Figure 11a shows the results of the Granger causality test performed on a pixel basis between AOD and meteorological factors. The blank areas represent the pixels that failed the Granger causality test, while the colored areas confirmed a Granger causality between AOD and related factors. Granger causality has been found between AOD and WIN/RHU in most areas of northern China, including the NW, NE, BTH, FWP, and Taklamakan Desert areas. In the SW, YRD, CC, and GD areas, AOD was related to the TEM, WIN, and RHU. Particularly for the CC area, the RHU and WIN were major causes for AOD change, which might indicate that WIN and RHU can influence AOD at a pixel scale.

Figure 11b shows the SR of AOD to meteorological drivers. The overall SR of AOD to TEM is small, approximately −0.25 in approximately 26.4% of the regions, mainly in CC, GD, and the Eastern part of the SW region. Approximately 30.8% of the regions in China have a positive SR of AOD to RHU, mainly in the BTH, GD, and Qinghai-Tibet Plateau regions. This indicates that AOD increases with increasing RHU in these regions. The sensitivity of AOD to WIN was significantly different on both sides of the Hu Line. This phenomenon indicates that the windy and sandy weather to the west of the Hu Line promotes aerosol pollution. On the east side of the Hu Line, the WIN accelerates airflow and mitigates aerosol pollution. Figure 11c shows the CR of TEM, RHU, and WIN on AOD. TEM contributed negatively to AOD in the Southern region, while TEM contributed positively to AOD in the Qinghai-Tibet Plateau and part of NE. In the SW, CC, and Southern parts of the NE, RHU reduced AOD pollution. WIN contributed more to AOD in the NE, Qinghai-Tibet Plateau, and Taklamakan Desert regions, where WIN caused more windy and dusty weather and intensified aerosol pollution. In the SW, CC, and coastal regions, WIN had a significant mitigation effect on aerosol pollution.

Table 4. Driving force analysis performed on a regional scale between AOD and meteorological factors of major regions.

	CGT			SR			CR/%
	TEM	RHU	WIN	TEM	RHU	WIN	TEM	RHU	WIN
NE	/	/	/	0.096	0.259	0.361	7.04	−1.7932	6.273
NW	/	/	√	0.177	−1.372	1.376	6.214	−4.723	7.421
SW	/	/	√	−0.324	−2.681	1.209	−5.434	−12.221	15.593
YRD	/	√	/	−0.09	−0.692	0.795	−2.580	−4.759	−5.659
FWP	/	√	/	0.092	0.228	−0.307	1.439	0.364	−3.832
CC	√	√	√	−0.371	−1.677	−1.297	−8.378	−10.936	−13.079
BTH	/	/	√	0.285	1.078	−0.492	4.662	−3.802	0.949
GD	√	√	/	−0.481	0.532	−0.618	−5.484	3.126	−2.474

Note: √ indicates that causality is established (p ≤ 0.05), / indicates that causality is not established (p > 0.05).

shows the results of the driving mechanism performed on a regional scale between AOD and meteorological factors. The GCT results clearly show that different areas have different causality relations. No causality was found in the NE area. In the NW and SW areas, WIN caused changes in AOD. The changes in AOD in the GD area were only attributed to TEM, and AOD changes in the YRD and FWP areas were only attributed to RHU. In the BTH area, WIN caused changes in AOD. In the CC area, TEM, RHU, and WIN were all causes attributed to changes in AOD. Overall, the SR of AOD to TEM is low. For TEM, the SR and CR values were negative in the SW, YRD, CC, and GD but were positive in the NE, NW, FWP, and BTH regions. AOD showed a more significant negative sensitivity to RHU in the NW, SW, and CC regions, indicating that RHU inhibits AOD. In the NW and SW regions, AOD showed high positive sensitivity to WIN, and the contribution of WIN to AOD was high in these two regions, which indicated that WIN was the main driver. In the CC region, AOD showed a high negative sensitivity to WIN, and the contribution of WIN to AOD in this region was high at −13.079, indicating that WIN accelerated the airflow and alleviated aerosol pollution.

3.4.3. Vegetation Continuous Fields Factor

Figure 12 shows the results of the driving mechanism performed on a pixel scale between AOD and VCF. The distribution of the Granger effect area is more fragmented, mainly in the NE, BTH, GD, CC, and SW. Overall, the sensitivity of AOD to VCF is negative in most areas, which indicates that VCF suppresses aerosol pollution, and the sensitivity of AOD to VCF was higher in southern areas, such as GD, CC, and Southern YRD, and weaker in areas west of the Hu Line, probably because the overall level of AOD was lower in Western areas, and thus, the sensitivity of AOD to VCF was not significant. The contribution of VCF to AOD was negative in most areas of China, and VCF had a more significant inhibitory effect on AOD in parts of the Qinghai-Tibet Plateau, SW, GD, and NE, with CR values exceeding 10%.

Table 5. Driving force analysis performed on a regional scale between AOD and VCF in major regions in China.

	GCT	SR	CR/%
NE	√	−0.946	−10.88
NW	/	0.134	2.31
SW	√	−1.523	−20.93
YRD	/	1.080	−4.76
FWP	/	0.226	1.63
CC	√	−1.047	−11.1
BTH	√	0.171	−0.74
GD	√	−0.999	−24.74

Note: √ indicates that causality is established (p ≤ 0.05), / indicates that causality is not established (p > 0.05).

shows the results of the driving mechanism performed on a regional scale between AOD and VCF. It shows that there is Granger causality between VCF and AOD in SW, BTH, CC, NE, and GD. In the areas where Granger causality exists, the results of CR and SR indicate that VCF suppresses aerosol pollution. The sensitivity of AOD to VCF is high in the SW, GD, and CC areas with SR of −1.523, −0.999, and −1.047, respectively. In the YRD, AOD had a significant positive sensitivity to VCF with an SR of 1.08. The sensitivity of VCF to VCF was high in the SW, GD, CC, and NE regions with significant contributions of AOD with CR values of −20.93%, −24.74%, −11.1%, and −10.88%, respectively.

4. Discussion

This study investigated the spatiotemporal variations and driving factors of AOD in mainland China from 2010 to 2020 using the MERRA-2 reanalysis product.

4.1. Environmental Impact of AOD

In addition to examining the causal relationship between drivers and AOD, the Granger causality test can also explain the feedback of the environment on AOD. Therefore, this section discusses the effects of AOD on WIN, RHU, TEM, and VCF, as shown in Figure 13. The variation in TEM due to AOD is small and only concentrated in the coastal areas and the Northern parts of SW. The Granger effects of AOD on RHU and WIN are similar and mainly concentrated in the Northern parts of China, CC and the Western part of SW. The impact of VCF on AOD was mainly concentrated in parts of the FWP, SW, and NW. However, whether AOD exerts a promoting or inhibiting effect on vegetation growth is related to the level of AOD itself [5]. Therefore, whether the decrease in AOD leads to an increase or decrease in VCF needs to be further investigated.

4.2. Spatial and Temporal Evolutionary Characteristics of AOD

The spatial distribution of AOD varied across the areas. Divided by the Hu Line, the AOD values of the Eastern areas were significantly higher than those of the Western areas. The high AOD values of the Western areas were mainly distributed in the Taklamakan Desert, where sparse vegetation, high dust, perennial drought, and low rainfall were major causes of the more dusty weather [55]. Both the linear trend analysis on a regional scale and the M–K trend analysis on pixel and regional scales showed a significant decreasing trend in most areas of mainland China, especially in all the Eastern areas of the Hu Line. In autumn, Hurst < 0.5 in the Taklamakan Desert region was mainly due to the unpredictable local wind and sandy weather, so the AOD trend was not sustainable. For most of China, especially inland cities with stable pollution sources and a Hurst index > 0.5, AOD will continue to decrease in the coming period.

4.3. Driving Mechanism of AOD by Region

Based on the results of the Granger causality test, a Granger causality relationship between the changes in AOD and socioeconomic factors was found for the BTH and SW areas. This indicates that economic growth benefited aerosol control activities. Such a causal relationship was not found in other major areas because the influence of economic growth on aerosol control is confounded with many other factors, such as government-initiated policies.

For meteorological factors, the changes in AOD in different major areas were attributed to different reasons, which might be related to local geographic factors such as climate and topography. In the Western area of the Hu Line, the AOD in the Taklamakan Desert was caused by RHU and WIN, which was consistent with previous findings by Fan, et al. [54] that the coupled effects of meteorological and topographic factors in the Tarim Basin explain the source and evolution of aerosols in the area. For the three meteorological factors, WIN caused changes in AOD in the SW, NW, and BTH areas, while RHU caused changes in AOD in the YRD and FWP areas, and TEM caused changes in AOD in the GD area. In the CC area, TEM, RHU, and WIN were all causes of the changes in AOD.

Vegetation can prevent wind and sand and regulate climate, which helps to mitigate aerosol pollution [6,56,57,58]. This paper reveals the geographic characteristics of the effect of vegetation on AOD at the pixel scale and regional scale. Overall, VCF plays a positive role in mitigating aerosol pollution, which is consistent with previous conclusions [30,59]. The contribution of VCF to AOD is more significant in NE, NW, and parts of the Qinghai-Tibet Plateau, which indicates that China has achieved initial results in controlling windy and sandy weather through afforestation. In addition, the continued implementation of the policy of returning farmland to forest has both increased the proportion of forestland and reduced atmospheric pollution activities such as straw burning, which has alleviated aerosol pollution to some extent [30].

Although meteorological factors were found to cause decreases in AOD in all major regions, it is worth noting that the overall decreasing trend of AOD is related to aerosol control activities. Existing studies have confirmed that human activities are the leading cause of air pollution [3,31,60]. With rapid economic growth in the past decade, the AOD in mainland China has shown a significant decreasing trend, which was a combined result of intensified aerosol emissions and enhanced aerosol controls by human activities. Since 2012, the governments of China have taken a series of measures to reduce the emission of aerosol gases. These control activities, such as industrial structure adjustment, limiting the number of oil-burning motor vehicles, environmental protection, vegetation restoration, and manual spraying operations to reduce dust in the air, were taken at different spatial and temporal scales and were difficult to quantify. However, these environmental protection activities can also affect the radiation budget between the biosphere and atmosphere and even influence the circulation of atmospheric currents [61,62]. Changes in meteorological factors, such as RHU and WIN, are consequences of these environmental protection activities. This study only elucidates the driving relationship between AOD and environmental factors, while a further study of the effects of environmental protection activities on climate is expected to enhance our understanding of the feedback between human activities and climate.

4.4. Limitation and Further Efforts

Based on MERRA2-AOD data, this paper analyzes the spatial and temporal evolution and driving mechanisms of AOD in China from 2010 to 2020. At the same time, the evolution and driving mechanisms of AOD in China’s major economic regions are studied at the regional scale. However, the study at the regional scale makes us pay insufficient attention to the evolution and driving mechanisms of AOD in key cities. In the Eastern region, human activities are an important factor influencing AOD. Therefore, the study of AOD in large cities is of great importance. The outbreak of COVID-19 has brought the world a series of changes. The Chinese government has taken control measures, which have greatly affected the way of life of people and socioeconomic activities. The spatial and temporal evolution of AOD during COVID-19 should be discussed separately; otherwise, it will affect our analysis of AOD for 2010–2020. Therefore, studying the spatial and temporal evolution and driving mechanisms of AOD at the urban scale and considering the impact of COVID is the next step of our research.

5. Conclusions

Based on MERRA-2-AOD reanalysis data, this paper investigated the spatial and temporal variations in AOD in China from 2010 to 2020 and the factors driving these changes. The main conclusions are as follows:

(1) The spatial distribution of AOD varied across the areas. Divided by the Hu Line, the AOD values of the Eastern areas were significantly higher than those of the Western areas. The high AOD values of the Western areas were mainly distributed in the Taklamakan Desert, which followed the order of spring > summer > autumn > winter. The annual mean AOD values of the eight major regions ranged from 0.271 to 0.587 in the order of CC > YRD > GD > FWP > BTH > SW > NE > NW.

(2) AOD in China was dominated by continuous decline and moderate fluctuations (0.1 < CV < 0.15) on the annual average scale and continuous decline and moderate to high fluctuations on the seasonal scale. From 2010 to 2020, AOD in Eastern and Western China showed a sustainable decline in general, while the decline was more pronounced in the Eastern region of the Hu Line than in the Western region.

(3) The Granger causality between AOD and GDPPC was confirmed in the BTH region, while the causality between AOD and IOPSK was confirmed in the SW region.

(4) The driving mechanism of meteorological factors on AOD had obvious geographic characteristics. Overall, the driving relationship between AOD and TEM was weak. The sensitivity of AOD to RHU was positive in BTH, the Qinghai-Tibet Plateau, the Sichuan Basin, and GD and negative in the rest of the regions. Divided by the Hu Line, the sensitivity of AOD to WIN was positive in the Western areas and was negative in most Eastern areas. The actual contribution of TEM and RHU to AOD was small overall.

(5) Overall, VCF had a certain ameliorating effect on aerosol pollution, which indicated that the measures of afforestation to control aerosol loadings in China had achieved initial results. This paper focuses on the spatial and temporal evolution of AOD in China from 2010 to 2020 and analyzes the driving mechanisms of economic, vegetation, and meteorological factors on AOD. However, the factors affecting AOD are complex, and other influencing factors are not thoroughly investigated in this paper. In addition, the layout and influencing factors of AOD during COVID-19 may change, and this paper does not investigate AOD in depth during this period, which will be the focus of our next study.