Long-Term Analysis of Aerosol Optical Depth over the Huaihai Economic Region ( HER ) : Possible Causes and Implications

The quality of air is increasingly affecting regional climate and human activity. Collection 6 aerosol products retrieved from the Moderate-Resolution Imaging Spectroradiometer (MODIS) onboard the Terra satellite were validated based on CE-318 sun photometric data to analyze their applicability in the Huaihai economic region (HER) at the Xuzhou and Shouxian sites. The spatio-temporal variations of aerosol optical depth (AOD) were also analyzed over HER from 2000 to 2016, with analyses of the correlation with potential driving factors, including meteorology, vegetation and human factors. HER is an economic cooperation organization with multiple industrial structures, containing coal resource-based cities, a national transportation hub and agricultural and high-altitude areas, which shows regional differences in AOD. The results suggest that MODIS Terra AOD products show good agreement with ground observations, with correlation coefficients of above 0.84 in HER, and the main pollutants for high AOD values are fine particles (the mean Ångström exponent was 1.16). The average annual change in AOD varied with a weak growth trend over the past 17 years, while a transition in 2012 made the tendency change from upward to downward due to the extensive cooperation of cities in the joint prevention and control of the deterioration of the ecological environment. The largest monthly mean AOD value appeared in June, which resulted from significant agricultural residue burning. The spatial distribution of multi-year average AOD occurred with a banded high-value center, extending from the north-west to the south-east. The high aerosol loadings were located in resource-based cities, and industrially developed and south-eastern coastal areas, whereas the regions with relatively low AOD in HER were distributed in the southern agricultural and northern high-altitude areas. The AOD value in the western, northern, and eastern coastal areas of HER showed a significant increasing trend, while no area exhibited a decreasing trend. The average wind speed has the largest negative correlation with the AOD value in terms of the natural driving factors, and GDP (gross domestic product) was more positively correlated with AOD with respect to the human factors, in comparison with population density.


Introduction
Aerosol is a multi-phase system composed of solid or liquid particles and gas carriers suspended in the atmosphere.Aerosol, an important component of atmospheric pollutants, not only affects the radiation and energy budget balance of the Earth-atmosphere system directly by absorbing and scattering solar radiation, but also affects the regional and global atmospheric environment, human health, and visibility indirectly or semi-directly [1][2][3][4].Ground-based observations and satellite remote sensing are two kinds of methods used to detect aerosol at present.The former can provide accurate aerosol optical depth (AOD) for one point in space, but it cannot be extended widely in space due to the limitations of observational conditions and instruments, as well as the spatio-temporal variability of aerosol distributions.The latter has overcome the spatial inadequacy of ground-based observations, and it is possible to obtain large-scale, high-spatial resolution, and multi-channel remote-sensing data to provide strong data support for the researches of the properties of global aerosol.However, the retrieved AOD can also be affected by multiple factors such as surface cover, topography, circulation and human activities, and AOD may have no value in winter in some areas, due to the weather conditions and cloud cover [5,6], and AOD also does not always agree well with surface PM 2.5 (with particle diameters <2.5 µm), the important parameter for health effects.Researchers have developed various techniques based on different sensors to retrieve AOD data over land and sea, such as the Moderate Resolution Imaging Spectroradiometer (MODIS) on the Terra and Aqua satellites.This technique can provide global data about aerosol properties and has been widely used [7,8].The applicability of MODIS products varies greatly in different geographical spaces, and the observations from the ground can serve as a data source for the verification of observations from satellite remote sensing [9].
The validation of MODIS aerosol products and regional suitability analysis provide the foundation for a variety of investigations.Thus, researchers have begun to evaluate MODIS aerosol products by using ground data derived from the Aerosol Robotic Network (AERONET) [10], achieving certain results.Numerous studies [10][11][12][13][14] have also found that AOD retrieved from MODIS generally agrees well with AERONET monitoring data over different regions, suggesting that MODIS aerosol data can be used to analyze the transport and distribution of aerosol.Based on MODIS aerosol data, it is of great value to analyze the transport and distribution of aerosol, and the spatial and temporal patterns in aerosol loadings are revealed at the regional scale [15,16].Acharya and Sreekesh [17] used the MODIS_L2 dataset to analyze the dynamic changes of seasonal AOD in India from 2001 to 2009.Luo et al. [18] analyzed the characteristics of the optical thickness of remote-sensing aerosol on the Chinese mainland during the first decade of the 21st century and found that the average AOD spatial distribution center in China typically exhibits two low peaks and two high peaks.In some studies, a long-term trend analysis of AOD has been conducted and variations in different typical regions across China have been identified, both spatially and temporally [16,19].Bao et al. [20] showed that AOD concentration in central and northern China was high in spring and summer.In the Yangtze River basin and South China, the aerosol distribution is related to industrial emissions and agricultural biomass burning.Atmospheric aerosol concentration in East Asia is also affected by the monsoon and weather in addition to the local emission sources [21,22].
The research areas of aerosol in China are currently concentrated in hot spots such as the North China Plain, Central China, the Yangtze River delta, the Pearl River delta, and the Sichuan basin [13,[23][24][25][26].Many studies have analyzed the influence on AOD in these economically developed areas from the aspects of topography, vegetation, and human factors.However, few studies have reported on the long-term changes in aerosol characteristics across the Huaihai economic region (HER), a transition zone between the south and north of China.The economic development of HER is accelerating and the effect of human factors on aerosols needs to be emphasized, while the impact of meteorological factors on AOD should also be considered.The AOD characteristics of Xuzhou, the central city in HER, were analyzed based on ground and satellite observations [6,27,28].The results revealed that seasonal variation and monthly mean fluctuations of Xuzhou AOD were significant, and the spatial distribution of AOD had a wide range of regional features.A long-term trend analysis of more than 10 years on the spatio-temporal variation of aerosol in HER by satellite remote sensing is still limited at present.Data from the latest MODIS Collection 6 version during the period 2000-2016 were collected and the spatio-temporal characteristics of AOD in HER were analyzed via this long-term series data.This investigation verified the MODIS satellite aerosol data based on ground measurements and conducted a long-term trend analysis of aerosol loadings at the regional scale in order to explore the annual, monthly, and seasonal variations of AOD and the spatial distribution features of the aerosol optical properties over HER.The potential driving forces of the evolution of aerosol on a pixel scale were also considered by combining meteorological conditions, vegetation, and human factors, among which the human factors were selected in view of the level of economic development in regional cities and the availability of data on influence factors.This may provide a reference for the study of aerosol and also provide a theoretical basis for the government to formulate policies to control air pollution in HER.

Study Area
HER, which includes 20 medium-sized cities in four provinces (Jiangsu, Shandong, Henan, and Anhui), is located in the heart of the east coast of China, with an area of 178,100 km 2 .Its climate is semi-arid and semi-humid.The spatial scope is 32.5-36.5 • N and 114-121 • E (see Figure 1).There are large coal mines and thermal power plants in Zaozhuang, Northern Anhui, and Xuzhou, as well as a large number of rubber tire, machinery manufacturing, textile, and chemical enterprises, with large amounts of coal-fired flue gas and exhaust emissions.
variation of aerosol in HER by satellite remote sensing is still limited at present.Data from the latest MODIS Collection 6 version during the period 2000-2016 were collected and the spatio-temporal characteristics of AOD in HER were analyzed via this long-term series data.
This investigation verified the MODIS satellite aerosol data based on ground measurements and conducted a long-term trend analysis of aerosol loadings at the regional scale in order to explore the annual, monthly, and seasonal variations of AOD and the spatial distribution features of the aerosol optical properties over HER.The potential driving forces of the evolution of aerosol on a pixel scale were also considered by combining meteorological conditions, vegetation, and human factors, among which the human factors were selected in view of the level of economic development in regional cities and the availability of data on influence factors.This may provide a reference for the study of aerosol and also provide a theoretical basis for the government to formulate policies to control air pollution in HER.

Study Area
HER, which includes 20 medium-sized cities in four provinces (Jiangsu, Shandong, Henan, and Anhui), is located in the heart of the east coast of China, with an area of 178,100 km 2 .Its climate is semi-arid and semi-humid.The spatial scope is 32.5-36.5°N and 114-121° E (see Figure 1).There are large coal mines and thermal power plants in Zaozhuang, Northern Anhui, and Xuzhou, as well as a large number of rubber tire, machinery manufacturing, textile, and chemical enterprises, with large amounts of coal-fired flue gas and exhaust emissions.As the central city of HER, Xuzhou is a typical coal-industrial city, and there are many industrial parks, mainly concentrated in the northern part of the city.These are located in the windward region during the heating period.However, the urban areas are mostly dominated by plains with relatively flat terrain, and there are several hills around to form a small local basin, which are not conducive to the dissipation of local pollutants.HER is also one of the earliest organizations of the regional economic cooperation and plays an important role in the Chinese regional economic pattern.However, it has struggled to survive under the pressure of the economically booming areas of the Yangtze River Delta to the south and Bohai Rim Economic Region to the north.This has caused the region to be an economic fracture zone and development has been difficult to maintain for a long period of time.Recently, the growth rate of the regional economy has become higher than the national average with the coordinated development of HER, and the total economic volume has continued to rise in proportion to the whole country.At the same time, local air pollution has become increasingly prominent.

Data Sources
We obtained MODIS/Terra Collection 6 Level 2 aerosol datasets (MOD04_L2) from March 2000 to February 2017 from National Aeronautics and Space Administration (NASA) Level 1 and the Atmosphere Archive and Distribution System Distributed Active Archive Center (LAADS DAAC) [29].The aerosol ground-based data observed from the CE-318 sun photometer at the atmospheric remote-sensing observatory in the China University of Mining and Technology (34.22 • N, 117.14 • E) were also chosen.AERONET provides three levels of AOD: Level 1.0, Level 1.5, and Level 2.0 [30].The Level 1.5 observation data, the highest level in Xuzhou site, were selected for verification of MODIS aerosols [11,31,32], because AERONET Level 2.0 AOD real-time data from July 2013 to December 2016 are not available.The Level 2.0 data were used in Shouxian site.The observation error of AERONET AOD was between 0.01 and 0.02 [33], and it was considered to be the true value for validation and evaluation of the precision of the satellite product.
Eight factors that could form AOD spatial variability were selected from both natural and anthropogenic aspects, including five main meteorological factors obtained from surface observations during 2000 to 2015: average wind speed (AWS), sunshine duration (SD), precipitation (P), average relative humidity (ARH), and average temperature (AT).These data were provided by the Lake-Watershed Science Data Center, National Earth System Science Data Sharing Infrastructure, National Science and Technology Infrastructure of China.AWS, ARH and AT are monthly mean data.The sunshine duration and precipitation were cumulative values.The spatial resolution of meteorological factors was the same as that of AOD by interpolation.The normalized difference vegetation index (NDVI) monthly synthetic product (MODND1M) with 500-m spatial resolution was provided by International Scientific and Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences [34].Two annual socioeconomic indices for cities-the gross domestic product (GDP) and population density (PD)-were taken from the statistical yearbook of 2000-2015 in Jiangsu, Shandong, Anhui, and Henan provinces (the statistical data for 2016 was temporarily lacking).The calculation of the Pearson correlation coefficient [35] and the test of significance between driving factors were made with the SPSS19.0software package [36].Because the data magnitude of anthropic factors is different from that of natural factors, the correlation between anthropic factors was calculated separately.

Methods
The AOD at 550 nm, used in this study, is one of the most widely used and reliable aerosol optical parameters [7,37].While AERONET ground-based observation data had no corresponding wavelength to match it exactly, the AERONET AOD data in the 550 nm band was obtained by band interpolation between 440 nm and 870 nm according to the Ångström equation [38], based on the assumption of a uniform spectral dependence between these two wavelengths [39].
where τ a (λ) represents AOD in the band of wavelength λ; β represents Ångström turbidity; and α is Angstrom wavelength that describes the particle size.An Ångström exponent less than 1 represents dominance of coarse mode aerosols, while an Ångström exponent greater than 1 represents the dominance of fine mode aerosols [40].
We extracted the 550-nm spectral data from a daily 10-km aerosol dataset.The statistics of monthly, seasonal (March to May for spring, June to August for summer, September to November for autumn, December to the following February for winter), and annual mean series AOD for all pixels were obtained by projecting, mosaicking, vector cropping, and averaging operations with ArcGIS 10.3 (http://www.esri.com/arcgis)(see Section 1 in the Supplementary Material).The continuous spatial distribution of each meteorological factor was obtained by using the cokriging interpolation method based on the measured data from the pertinent meteorological stations [41].The multi-year NDVI value was calculated by the mean value method and the NDVI time-series dataset with the same AOD resolution was resampled.The ArcGIS platform was used to load the statistical data into the vector layers of each city and convert the raster data into the same type of projection and grid size as AOD.
We calculated the trend rate using unary linear regression analysis based on AOD pixel values to determine the overall trend of AOD variation from 2000 to 2016 [42], recorded as slope.If slope is greater than zero, it indicates that the AOD value in the studied period is increasing.A larger slope indicates a greater intensity of change; otherwise, the AOD decreases.The F-test was used to determine the significance of the spatial variation trend of AOD because the variance ratio of regression follows the F distribution with degrees of freedom of (1, n−2) (n = 17 in this study).According to the different thresholds at the significant levels of 0.01 and 0.05, we separated the trend rates of the slope into non-significant, significant, and extremely significant levels.Meanwhile, because slope values can be divided into negative (slope < 0) and positive (slope > 0) values, the results can be classified as one of the following six levels: extremely significant decrease, significant decrease, no significant decrease, no significant increase, significant increase, and extremely significant increase.
Natural factors and human activities that have spatial variability on the smaller pixel scale will present different effects on AOD in different regions.Therefore, analyzing the spatial and temporal configuration of the AOD evolution and the driving factors can help to grasp the specific influence scope and degree of different driving forces.Pearson product-moment correlation analysis was used to identify the relationship between AOD and potential impact factors [32], and the correlation coefficient is denoted as R. AOD is closely related to the driving factor when R is large.The correlation coefficient was tested by means of the chartable method according to the known degree of freedom.Then, the correlation was classified as significant negative correlation (0.01-negative, 0.05-negative), significant positive correlation (0.01-positive, 0.05-positive), no significant negative correlation (no-negative), and no significant positive correlation (no-positive), combined with the critical value of ρ = 0 at the significant levels of 0.01 and 0.05 [43].

Comparison and Validation
MOD04_L2 aerosol product value represents the instantaneous AOD in the space of 10 km × 10 km, while AERONET AOD is the instantaneous observation at 15-min intervals.The transit time of the Terra satellite in the Chinese area is around 10:30 AM locally, and the aerosol observations of MODIS and AERONET are different both in temporal and spatial scale.The mean MODIS AOD in the buffer zone with a radius of 10 km centered on the ground observation sites were used to compare with the mean AOD values that were observed at the ground sites within ±0.5 h of satellite overpass times.
Finally, a total of 110 pairs of matched samples were utilized for analysis and application in the Xuzhou site, but just 26 valid values could match with the ground data observed in the Shouxian site (as shown in Table 1).One of the main reasons for this limited data matching is that only a very short period of observation data can be obtained in this site.The upper and lower limits of the expectation error range of aerosol product are taken as ∆τ = ±0.05± 0.15 • τ AERONET [37].Figure 2 presents the linear fitting equation and coefficient of determination (R 2 ) of MODIS AOD and AERONET data, from which we can see that MODIS C6 aerosol products are highly consistent with ground-based measurements, with a small root-mean-square error (RMSE) and mean-absolute-error (MAE) [44,45].A total of 61.8% of the samples fall into the NASA expectation error range in the Xuzhou site, and the correlation coefficient reaches 0.87.This result is similar to the validation for the Xuzhou area of Li et al. [27].The correlation coefficient is 0.84 in the Shouxian site, lower than the verification result (R 2 = 0.92) obtained by Deng et al. [11], but the slope (~0.97) and intercept (~0.02) are closer to the ideal state and 76.9% of dots meet the allowed error requirement.Therefore, the data derived from MODIS can be applied to investigate the spatial and temporal distribution of aerosols in HER.However, the intercept of the line of best fit is greater than zero, which indicates that the C6 version of aerosol products is overvalued in the study area.This may be due to the underestimation of the surface reflectivity and the error of assumptions of the aerosol model [37].compare with the mean AOD values that were observed at the ground sites within ±0.5 h of satellite overpass times.Finally, a total of 110 pairs of matched samples were utilized for analysis and application in the Xuzhou site, but just 26 valid values could match with the ground data observed in the Shouxian site (as shown in Table 1).One of the main reasons for this limited data matching is that only a very short period of observation data can be obtained in this site.The upper and lower limits of the expectation error range of aerosol product are taken as ∆τ = ±0.05± 0.15 • [37].Figure 2 presents the linear fitting equation and coefficient of determination (R 2 ) of MODIS AOD and AERONET data, from which we can see that MODIS C6 aerosol products are highly consistent with ground-based measurements, with a small root-mean-square error (RMSE) and mean-absolute-error (MAE) [44,45].A total of 61.8% of the samples fall into the NASA expectation error range in the Xuzhou site, and the correlation coefficient reaches 0.87.This result is similar to the validation for the Xuzhou area of Li et al. [27].The correlation coefficient is 0.84 in the Shouxian site, lower than the verification result (R 2 = 0.92) obtained by Deng et al. [11], but the slope (~0.97) and intercept (~0.02) are closer to the ideal state and 76.9% of dots meet the allowed error requirement.Therefore, the data derived from MODIS can be applied to investigate the spatial and temporal distribution of aerosols in HER.However, the intercept of the line of best fit is greater than zero, which indicates that the C6 version of aerosol products is overvalued in the study area.This may be due to the underestimation of the surface reflectivity and the error of assumptions of the aerosol model [37].As shown in Figure 3, the comparison results indicate that MODIS_550 nm AOD basically showed the same change trends of matched CE-318 sun photometer monthly average AOD in different bands.The MODIS_550 nm AOD value was slightly higher than that of the CE-318 sun photometer at 550 nm; however, they both show the similar characteristic "peak valley" type, and peaks appeared in March, July, and September while the three low values occurred in January, May, and August.The temporal variation analysis of AOD is discussed in Sections 3.2 and 3.3.As shown in Figure 3, the comparison results indicate that MODIS_550 nm AOD basically showed the same change trends of matched CE-318 sun photometer monthly average AOD in different bands.The MODIS_550 nm AOD value was slightly higher than that of the CE-318 sun photometer at 550 nm; however, they both show the similar characteristic "peak valley" type, and peaks appeared in March, July, and September while the three low values occurred in January, May, and August.The temporal variation analysis of AOD is discussed in Sections 3.2 and 3.3 A scatter plot between AOD and the Ångström exponent of the matched data in HER is shown in Figure 4.It can be seen that the maximum value of the Ångström exponent may also increase as AOD increases.We divided the Ångström exponent with 0.3 as the interval, and the mean value of the corresponding AOD within each interval from 0.2 to 1.7 were 0.377, 0.294, 0.428, 0.604, 0.657, respectively, and the average Ångström exponent was 1.16, indicating that the major factor causing the high AOD value in the study area is the increase of fine particles.The AOD is mainly concentrated in about 0.35 when α is less than one, and the value of α greater than one accounts for 74.76% of the total, showing that the specific proportion of fine particles in this area is larger.

Annual Variation and Monthly Trend
Figure 5a shows the temporal variations of MODIS AOD in HER, both monthly and yearly, from January 2000 through December 2016, indicating the trend of change over time.The vertical blue bars at the top denote mean AOD over the entire period 2000-2016 for each year, from which we can see the average AOD value of each year, a total of 17, ranges from 0.591 to 0.836, and there are two peak periods, 2007 and 2011.The monthly mean AOD tended to increase gradually from January, peaking in June, and then decreased until December, with the lowest AOD of approximately 0.44.The period May-July dominated the maximum aerosol loadings, and June was the highest AOD month with an AOD value of approximately 0.97 over the entire study period.The A scatter plot between AOD and the Ångström exponent of the matched data in HER is shown in Figure 4.It can be seen that the maximum value of the Ångström exponent may also increase as AOD increases.We divided the Ångström exponent with 0.3 as the interval, and the mean value of the corresponding AOD within each interval from 0.2 to 1.7 were 0.377, 0.294, 0.428, 0.604, 0.657, respectively, and the average Ångström exponent was 1.16, indicating that the major factor causing the high AOD value in the study area is the increase of fine particles.The AOD is mainly concentrated in about 0.35 when α is less than one, and the value of α greater than one accounts for 74.76% of the total, showing that the specific proportion of fine particles in this area is larger.A scatter plot between AOD and the Ångström exponent of the matched data in HER is shown in Figure 4.It can be seen that the maximum value of the Ångström exponent may also increase as AOD increases.We divided the Ångström exponent with 0.3 as the interval, and the mean value of the corresponding AOD within each interval from 0.2 to 1.7 were 0.377, 0.294, 0.428, 0.604, 0.657, respectively, and the average Ångström exponent was 1.16, indicating that the major factor causing the high AOD value in the study area is the increase of fine particles.The AOD is mainly concentrated in about 0.35 when α is less than one, and the value of α greater than one accounts for 74.76% of the total, showing that the specific proportion of fine particles in this area is larger.

Annual Variation and Monthly Trend
Figure 5a shows the temporal variations of MODIS AOD in HER, both monthly and yearly, from January 2000 through December 2016, indicating the trend of change over time.The vertical blue bars at the top denote mean AOD over the entire period 2000-2016 for each year, from which we can see the average AOD value of each year, a total of 17, ranges from 0.591 to 0.836, and there are two peak periods, 2007 and 2011.The monthly mean AOD tended to increase gradually from January, peaking in June, and then decreased until December, with the lowest AOD of approximately 0.44.The period May-July dominated the maximum aerosol loadings, and June was the highest AOD month with an AOD value of approximately 0.97 over the entire study period.The

Annual Variation and Monthly Trend
Figure 5a shows the temporal variations of MODIS AOD in HER, both monthly and yearly, from January 2000 through December 2016, indicating the trend of change over time.The vertical blue bars at the top denote mean AOD over the entire period 2000-2016 for each year, from which we can see the average AOD value of each year, a total of 17, ranges from 0.591 to 0.836, and there are two peak periods, 2007 and 2011.The monthly mean AOD tended to increase gradually from January, peaking in June, and then decreased until December, with the lowest AOD of approximately 0.44.The period May-July dominated the maximum aerosol loadings, and June was the highest AOD month with an AOD value of approximately 0.97 over the entire study period.The high AOD is closely related to the study area of burning rape and wheat straw, and the emissions of crop fires include aerosols, specifically carbon monoxide, carbon dioxide, nitrates, and black carbon [46].The daily average mass concentrations of PM 2.5 and PM 10 (with particle diameters <10 µm) turned out to be high during the period of agricultural residue burning [47], which was correspond with the result found by Wang et al. [48] that there was small peak seen in PM 2.5 from May to July for harvesting wheat, and east China were shown to have higher regional PM 2.5 concentrations compared with other regions in China.It's worth noting that there were small but visible peaks in winter during the years 2013-2015, because of the severe haze in eastern China that began in 2013 as well as the next two years, and the high value of AOD appeared in the winter of these years [49,50].The government realized the seriousness of air pollution and took a series of measures (e.g., shutting factories down and burning gas instead of coal) to alleviate the pressure of air pollution on people's lives.
As illustrated in Figure 5b, the frequent alternation between high and low monthly AOD values indicates noticeable fluctuations called "sawtooth cycles", which suggest a marked seasonal variation in AOD and have been found in previous studies [16,19].The valley bottom values were approximately 0.3-0.4 and peaks in the sawtooth cycle were 2-3 times higher than those low values.Periodic fluctuations were present within the monthly timescale, indicating inter-monthly AOD variation characteristics.
A linear trend analysis was conducted based on the time series of monthly-averaged AOD to reflect the overall trend of the entire study area.As shown in Figure 5b, there was a weak upward trend (+0.006/year) from 2000 to 2016 over the entire region.Meanwhile, from the annual average AOD in Figure 5c we can clearly find that AOD increased before 2012 (the time corresponding to the maximum mean value point is the turning point and AOD reached the highest value in 2011) and decreased after 2012.In the first stage, 2000-2011, there was an upward trend, with a yearly AOD slope of +0.018, and in the second, 2012-2016, there was an apparent downward tendency with a larger magnitude of 0.042 per year.The possible causes for this phenomenon are as follows.There are nine resource cities and six old industrial bases in HER, coupled with convenient traffic conditions, which means the coal-mining industry occupies an important position in the regional industrial structure [51].The coal gangues, discharged from the coal mine, were commonly stacked nearby, and the stock increased constantly.The weathered material of coal gangue entered the atmosphere under the action of wind force after weathering, rain erosion, and spontaneous combustion, significantly polluting the atmospheric environment.The national government has attached great importance to promoting regional coordinated development, and has made unremitting efforts to promote industrial adjustment and resource-based urban transformation since the 18th National Congress of the Communist party of China in November 2012.Simultaneously, cities in the research area have carried out extensive cooperation in the joint prevention and control of the deterioration of the ecological environment, and achieved certain results.Thus, an inflection point occurred during the ascent of the AOD value and the AOD decreased in the second stage.In the study of eastern and central China, Zhao et al. [52] found that AOD increased before 2006 and decreased after 2011, which was similar to the turning point (2012) found in this article.The reason for the increased AOD is that emission increases due to rapid economic development and for the decrease is the effective emission reduction in primary aerosols, SO 2 , and NO X .

Seasonal Characteristics
The seasonal means of AOD in the period 2000-2016 are shown in Figure 6.AOD in HER exhibited seasonal variability with the maximum occurring during summer, followed by spring, winter, and autumn, when AOD was slightly lower than in winter.It can be noted that seasonal variability showed a similarity over the two periods, with the exception of an anomaly in summer, which indicated that AOD in the 2000-2011 period was higher than in the 2012-2016 period.
The AOD value was also high in spring, with the long-term average annual value reaching 0.80.The northern and southern parts of the study area are affected by sand and dust, and the climate is dry with a large wind speed, which is also conducive to the impacting of local sand weather [53].The overall rate of increase in AOD was lowest in spring over the 17-year period, and the rate of decline was also minimal between 2012 and 2016 (see Table 2).AOD was highest in summer with a mean value of 0.86 (Figure 6), possibly due to the effect of higher air temperatures on the precursor of the secondary organic aerosols and photochemical interactions, which increase aerosol concentrations in the atmosphere [54].During the peak time of the Asian summer monsoon season, more aerosols can also be transported from southern China to the north by the dispersal functions of the monsoon circulation [55].The summer AOD varied from 0.62 to 1.17 with the highest value appearing in 2011 throughout this 17-year period.Meanwhile, the lowest summer AOD value occurred in 2000, and the annual average value increased by 0.56%.The AOD is also affected by the peak in straw burning, which leads to an increase in particulate concentration in the air and aggravates the air pollution further.At the same time, the anthropogenic aerosol emissions reach their highest level [11].However, there was a general downward trend for AOD in summer over the period 2012-2016, declining by an annual rate of 6.17%.It can be seen that the number of fire points in Anhui, Henan, Jiangsu and Shandong has been significantly reduced in the agricultural harvest season during 2014-2016 (see Tables S1 and S2 in the Supplementary Material) as a result of daily monitoring for fire points of straw burning by the Ministry of Environmental Protection of the People's Republic of China [56].The reason is that the government had vigorously promoted and popularized the mechanized chopping of wheat straw and the act of returning to the field, achieving initial success.This is consistent with the findings of Lijuan et al. [57] and Xin et al. [58].Compared with other seasons, the AOD in autumn did not show significant annual fluctuation.However, AOD decreases significantly compared with summer, which may be related to the increase in precipitation after August and the increase in atmospheric aerosol scavenging.In addition, the atmospheric temperature drops rapidly after autumn, in which the secondary aerosol transformation is far less than in summer [54].In winter, the atmospheric environment is beneficial to the spread of aerosols under the effect of the northern winter monsoon.The monsoon, combined Compared with other seasons, the AOD in autumn did not show significant annual fluctuation.However, AOD decreases significantly compared with summer, which may be related to the increase in precipitation after August and the increase in atmospheric aerosol scavenging.In addition, the atmospheric temperature drops rapidly after autumn, in which the secondary aerosol transformation is far less than in summer [54].In winter, the atmospheric environment is beneficial to the spread of aerosols under the effect of the northern winter monsoon.The monsoon, combined with short sunshine duration, results in less formation of secondary aerosols and lower average relative humidity, impeding the absorption and growth of aerosol particles.Thus, the winter AOD value is low compared with other seasons as a whole.

Spatial Variation of Aerosol Optical Depth (AOD)
Figure 7 illustrates the geographic distribution of the multi-year average AOD in HER from 2000 to 2016.The highest AOD occurred over Kaifeng, Heze, Jining, Xuzhou, Huai'an, and Yancheng, extending from north-west to south-east.This distribution may be represented as a banded high-value center.The 17-year averaged AOD in the north-western and south-eastern coastal areas reached 0.80-0.93,which was mainly related to the major coal bases in Jining, Xuzhou, and other north-western areas in HER.The air pollution caused by dust particles and harmful gases from coal mining and coal-producing processes made an important contribution to the increase in AOD.A large number of coal-fired boilers and coal-fired power plants in Yancheng will also emit coal flue gases to the atmosphere, which raises the AOD value in the area.In the south, the AOD loadings of cities in Anhui are relatively lower, but the overall average was still at a high level, with a multi-year AOD of 0.65-0.70 and higher in the northern part of the study area.Anhui is a large agriculture province where crops, nursery stock, flowers, and Chinese herbal medicine are planted in a large area, so vegetation coverage is better with weak industrial activity that decreases the value of AOD.Open-field burning of crop residue entering into the atmosphere also makes AOD higher in the harvest season [11], as explained in Section 3.2.Some cities with higher terrain in Shandong (Figure 8) have a low average AOD value, with values less than 0.50, coinciding with the conclusion from other studies that areas with higher elevation have lower AOD [16,26].Figure 9 indicates the seasonal average distribution of AOD in HER between 2000 and 2016, which is similar to the spatial pattern of the multi-year averaged AOD.The distribution of high AOD extends from north-west to south-east, but the seasonal mean AOD loadings differ greatly in different seasons.The AOD value was higher in spring and summer but lower in autumn and winter.The AOD loadings were mainly high in the eastern and western areas in spring.In summer, the research area was basically covered by large AOD values, and the AOD loading in eastern Jiangsu and western Shandong exceeded 0.90.The spatial distribution of AOD was relatively consistent in autumn and winter, which shows that the AOD value decreased significantly in these two seasons, compared with summer, with AOD in the entire area mainly lower than 0.60.Figure 9 indicates the seasonal average distribution of AOD in HER between 2000 and 2016, which is similar to the spatial pattern of the multi-year averaged AOD.The distribution of high AOD extends from north-west to south-east, but the seasonal mean AOD loadings differ greatly in different seasons.The AOD value was higher in spring and summer but lower in autumn and winter.The AOD loadings were mainly high in the eastern and western areas in spring.In summer, the research area was basically covered by large AOD values, and the AOD loading in eastern Jiangsu and western Shandong exceeded 0.90.The spatial distribution of AOD was relatively consistent in autumn and winter, which shows that the AOD value decreased significantly in these two seasons, compared with summer, with AOD in the entire area mainly lower than 0.60.As shown in Figure 10, the regional change and spatial difference of AOD have been remarkable during the past 17 years.The AOD presented a growth trend (the average slope is +0.006•year −1 , corresponding to the overall long-term pattern of change, as presented in Figure 5b) in most areas of HER, and the area with increased AOD accounts for 95.96% of the total area (no significant increase areas account for 59.04%).This growth trend was related to the gradual increase in support policies from the state to the local level in HER, along with the process of regional cooperation and the acceleration of economic development.The AOD values in the northern region, parts of the western region, and eastern coastal areas of the study area showed significant growth trends at a 95% confident level, accounting for 36.48%,whereas the AOD in central and southern parts changed with weak trends, covering 63.08% of the total area.In recent years, the government's support for the east coast areas and Shandong has increased, and perhaps the growth of the industries has increased the amount of particulate matters from industrial production to the atmosphere.However, there is almost no noticeable decreasing trend in this research area, where As shown in Figure 10, the regional change and spatial difference of AOD have been remarkable during the past 17 years.The AOD presented a growth trend (the average slope is +0.006 year −1 , corresponding to the overall long-term pattern of change, as presented in Figure 5b) in most areas of HER, and the area with increased AOD accounts for 95.96% of the total area (no significant increase areas account for 59.04%).This growth trend was related to the gradual increase in support policies from the state to the local level in HER, along with the process of regional cooperation and the acceleration of economic development.The AOD values in the northern region, parts of the western region, and eastern coastal areas of the study area showed significant growth trends at a 95% confident level, accounting for 36.48%,whereas the AOD in central and southern parts changed with weak trends, covering 63.08% of the total area.In recent years, the government's support for the east coast areas and Shandong has increased, and perhaps the growth of the industries has increased the amount of particulate matters from industrial production to the atmosphere.However, there is almost no noticeable decreasing trend in this research area, where industrial and anthropogenic activities are strengthening.

The Correlation between AOD and Natural Factors
The Pearson correlation coefficients of natural factors in Table S3 (see Supplementary Material) show that precipitation has a strong correlation with average temperature, average relative humidity and NDVI [35,59].Average temperature is also highly correlated with precipitation and NDVI.So precipitation and average temperature are no longer the natural factors of choice.
Figure 11(A)a identifies the spatial relationship between the multi-year averaged AOD and the average wind speed over HER.The 0.01 and 0.05 levels of significant negative correlation accounted for 67.2% and 8.1%, respectively, and the maximum absolute value of R reached 0.89.This phenomenon is particularly prominent in Anhui, Henan, Shandong cities, eastern Xuzhou, and the western region of Lianyungang.The correlation between AOD and wind speed reveals notable spatial heterogeneity.It shows a significant negative correlation in the western part of the study area, which gradually decreased from the western part to the middle, and turns out to be a positive correlation in the eastern part.The AOD in western Yancheng is positively related to the average wind speed.The speed of the wind affects the velocity of the dispersion of sulfate and soot emitted by many coal-fired boilers and thermal power plants [60].The monsoon from the south-east coast accelerates the expansion of the AOD large-value area, and the high humidity of the coastal air is conducive to the growth of aerosol particles [61,62].
The correlation between sunshine duration and AOD has obvious spatial differentiation (Figure 11(A)b and Figure 11(B)b).There is a significant negative relationship from the south-west of the study area in Fuyang and Kaifeng towards the north-east, transitioning to the significant positive relationship in Jining.The correlation between sunshine and AOD was also significant in the local areas of Tai'an, Rizhao, and Bengbu, but the other regions did not show a strong spatial relationship.The sunshine duration in the southwest of the study area is negatively correlated with AOD.This is caused by the increase in human activity in the region, especially in industrial emissions (in 2015, considering the average concentration of PM2.5 in the provinces, Henan was ranked 1st and Anhui was ranked 9th), which reduces the amount of sunlight and, subsequently, increases the aerosol concentration.The previous study has found that the increasing anthropogenic secondary aerosol precursor experienced better "gas particle transformation" in sufficient sunshine and higher temperature [54].The AOD in Jining and Tai'an increased with the increase in sunshine duration, and whether it is also for this reason, yet to be further investigated.

The Correlation between AOD and Natural Factors
The Pearson correlation coefficients of natural factors in Table S3 (see Supplementary Material) show that precipitation has a strong correlation with average temperature, average relative humidity and NDVI [35,59].Average temperature is also highly correlated with precipitation and NDVI.So precipitation and average temperature are no longer the natural factors of choice.
Figure 11(A)a identifies the spatial relationship between the multi-year averaged AOD and the average wind speed over HER.The 0.01 and 0.05 levels of significant negative correlation accounted for 67.2% and 8.1%, respectively, and the maximum absolute value of R reached 0.89.This phenomenon is particularly prominent in Anhui, Henan, Shandong cities, eastern Xuzhou, and the western region of Lianyungang.The correlation between AOD and wind speed reveals notable spatial heterogeneity.It shows a significant negative correlation in the western part of the study area, which gradually decreased from the western part to the middle, and turns out to be a positive correlation in the eastern part.The AOD in western Yancheng is positively related to the average wind speed.The speed of the wind affects the velocity of the dispersion of sulfate and soot emitted by many coal-fired boilers and thermal power plants [60].The monsoon from the south-east coast accelerates the expansion of the AOD large-value area, and the high humidity of the coastal air is conducive to the growth of aerosol particles [61,62].
The correlation between sunshine duration and AOD has obvious spatial differentiation (Figure 11(A)b and Figure 11(B)b).There is a significant negative relationship from the south-west of the study area in Fuyang and Kaifeng towards the north-east, transitioning to the significant positive relationship in Jining.The correlation between sunshine and AOD was also significant in the local areas of Tai'an, Rizhao, and Bengbu, but the other regions did not show a strong spatial relationship.The sunshine duration in the southwest of the study area is negatively correlated with AOD.This is caused by the increase in human activity in the region, especially in industrial emissions (in 2015, considering the average concentration of PM 2.5 in the provinces, Henan was ranked 1st and Anhui was ranked 9th), which reduces the amount of sunlight and, subsequently, increases the aerosol concentration.The previous study has found that the increasing anthropogenic secondary aerosol precursor experienced better "gas particle transformation" in sufficient sunshine and higher temperature [54].The AOD in Jining and Tai'an increased with the increase in sunshine duration, and whether it is also for this reason, yet to be further investigated.The average relative humidity and AOD showed a significant negative correlation (Figure 11(A)c and Figure 11(B)c), and the ratio of these areas accounts for 33.9% of the total.AOD was increased with increasing relative humidity on account of hygroscopic growth of aerosol particles according to Yoon and Kim [63] and Liu et al. [3].However, AOD decreases with the increase in average relative humidity in some part of Suqian, Huai'an, and Yancheng, probably due to the scavenging of aerosol particles by precipitation, which leads to the decrease in large particles and AOD.AOD in the inland regions of Kaifeng, Zhoukou, Shangqiu, and Heze were particularly significantly negatively correlated with relative humidity.Although the relative humidity in these areas is probably smaller than in the coast, the negative correlation of AOD and relative humidity also reflects wet removal of the aerosol.
There was a weak relationship between the spatially averaged NDVI and AOD in general (Figure 11(A)d and Figure 11(B)d), but there exists a significant negative correlation area, distributing from east to west, with striped sporadic distribution in the middle and south of HER.The average correlation coefficient in the northern region of Lianyungang and Suqian is 0.40.Previous research shows that the good vegetation cover and the dust retention-effect of plants can help reduce pollution and inhibit the increase of aerosol [12], but in this article, the effect is not obvious.

Correlation between AOD and Anthropogenic Factors
The regions in which GDP is positively correlated with AOD are discretely distributed throughout the study area and concentrated in the western, northern, and eastern coastal areas, as illustrated in Figure 11(A)e and Figure 11(B)e.The average correlation coefficient of Tai'an, Jining, and Heze is 0.55, and the maximum positive coefficient of Jining reaches 0.71.The GDP of Zhoukou and Fuyang are associated with AOD, and the average correlation coefficient is 0.43.Human and industrial activities in Xuzhou and Yancheng are stronger than in Bozhou (a major agricultural region), and GDP is positively correlated with the AOD of these two cities in Jiangsu province.However, in the south, the cities in Anhui, Suqian and Huai'an have a low correlation between GDP and AOD.This shows that AOD is not greatly affected by GDP in regions where economic development is slow.The previous result found by Xu et al. [12] shows that the faster the process of urbanization, the higher the economic level and the greater the loading of AOD.
As Figure 11(A)f shows, the significant positive correlation between AOD and population density, accounting for 54.9% of the total area, is much greater than the negative correlation in HER.The average correlation coefficient of Tai'an, Jining and Heze is 0.23, indicating that the population distribution and activity of the region have little impact on AOD.In contrast, the average correlation coefficient of Zhoukou and Fuyang was 0.65, and the impact is much greater than that of GDP on AOD.It is noteworthy that areas with significant negative correlation between AOD and population density appeared in Yancheng, the city that has a high population migration rate with a declining fertile population and an accelerating aging population, year by year.Yancheng possesses a vast territory, so the population density has remained relatively stable since it reached its lowest in 2005.However, the AOD loadings of Yancheng have increased in recent years, followed by a gradual decrease, with the development of the urban economy.As such, the Yancheng population density and AOD loading showed a significant negative correlation.

Conclusions
This study carried out a validation of the applicability of MODIS product and identified the spatio-temporal variability and trend of AOD over the Huaihai economic region in China, using Terra MOD04_L2 aerosol data obtained from 2000 to 2016, in relation to natural and human driving factors.The findings can be summarised as follows: AOD provided by the MOD04 product was highly correlated with AERONET AOD data (R 2 = 0.87, 0.84) in HER, and most of the matched samples are in the expected error range, with small RMSE and MAE.The average Ångström exponent was 1.16 and fine aerosol particles dominated in the atmosphere becoming the main pollutants for high AOD values.AOD varied with a weak growth trend, whereas analyses of AOD variations over two periods, an uptrend in 2000-2011 and a downtrend in 2012-2016, revealed a remarkable transition in the long-term overall tendency under the intervention of government's environmental protection policies and actions.The highest AOD values were mainly in June which resulted from serious straw burning, and the lowest appeared in December, suggesting an obvious monthly and seasonal variation cycle.The AOD of all four seasons exhibited a slight upward trend, with the highest growth rate in winter.The seasonally averaged AOD reached its maximum in summer, followed by spring, winter and autumn.High AOD values appeared in the eastern and western regions of HER in spring, while the research area was basically covered by high AOD value in summer and AOD decreased significantly in autumn and winter.The spatial distribution of AOD featured a "banded high-value center" extending from the north-west to the south-east.AOD is relatively low in the northern part of HER because of high terrain and shows a markedly increasing trend in the south-western, northern, and eastern coastal areas.However, only a handful of areas in the central and southern regions show a decreasing trend.
The average wind speed has the largest negative correlation with the AOD value among these natural factors, while NDVI is weakly related to AOD through relevant analysis on the pixel scale.Aerosol variation was positively and closely correlated with GDP in human elements, and the influence of population density was relatively weak, but it still had a great impact on AOD loadings.In addition, we selected six natural factors and two man-made factors to analyze the driving forces in this research, but the emission of coal and fossil fuel combustion in the research area will also have certain effect on aerosol loadings.Therefore, we may do further research on the impact factors of AOD combined with the 3-km spatial resolution MODIS products.

Figure 1 .
Figure 1.Geographical map and Aerosol Robotic Network (AERONET) sites (the red stars) used in this study.The shaded parts are the Huaihai Economic Region (HER) and it consists of 20 prefectural cities from 4 provinces: 5 from Jiangsu (Xuzhou, Lianyungang, Yancheng, Huai'an and Suqian), 7 from Shandong (Zaozhaung, Heze, Jining, Taian, Laiwu, Linyi and Rizhao), 3 from Henan (Kaifeng, Zhoukou and Shangqiu), and 5 from Anhui (Huaibei, Bozhou, Fuyang, Bengbu and Suzhou).The black and green solid triangles symbolizes thermal power plants and heating plants, respectively, from the Thermal Power Generation Business Directory (2016, version) of four provinces.The pink solid dots represent the local pillar industries including machinery, rubber tires, textiles and chemical enterprises.The blue parts are the Beijing-Tianjin-Hebei Region and the Yangtze River Delta, respectively.

Figure 2 .
Figure 2. Linear fit between Terra aerosol optical depth (AOD) and AERONET AOD ground observation at 550 nm at the (a) Xuzhou and (b) Shouxian sites.The error bars represent the standard deviations.

Figure 2 .
Figure 2. Linear fit between Terra aerosol optical depth (AOD) and AERONET AOD ground observation at 550 nm at the (a) Xuzhou and (b) Shouxian sites.The error bars represent the standard deviations.

Figure 3 .
Figure 3. MODIS AOD (at 550 nm) and AERONET AOD (at 440, 675, 870 and 1020 nm) for Xuzhou from 2013 to 2016 by month.The uncertainty is the standard deviation.

Figure 4 .
Figure 4. Relationship between the Ångström exponent and AOD at the wavelength of 550 nm (the solid black line is a cubic polynomial fitting curve of the scatter plot).

Figure 3 .
Figure 3. MODIS AOD (at 550 nm) and AERONET AOD (at 440, 675, 870 and 1020 nm) for Xuzhou from 2013 to 2016 by month.The uncertainty is the standard deviation.

Figure 3 .
Figure 3. MODIS AOD (at 550 nm) and AERONET AOD (at 440, 675, 870 and 1020 nm) for Xuzhou from 2013 to 2016 by month.The uncertainty is the standard deviation.

Figure 4 .
Figure 4. Relationship between the Ångström exponent and AOD at the wavelength of 550 nm (the solid black line is a cubic polynomial fitting curve of the scatter plot).

Figure 4 .
Figure 4. Relationship between the Ångström exponent and AOD at the wavelength of 550 nm (the solid black line is a cubic polynomial fitting curve of the scatter plot).

Figure 5 .Figure 5 .
Figure 5. (a) Domain-averaged AOD distribution with respect to month and year (filled with blue color).The horizontal blue bars on the right (on the top) demonstrate mean AOD over the period 2000-2016 for each month (year).(b) Time series of the territory-averaged monthly mean AOD at 550 nm in HER from March 2000 to February 2017.The trend line (dashed) calculated by using unary linear regression analysis is also shown.(c) Time series of the regional-averaged yearly mean

Figure 6 .
Figure 6.Seasonal distribution of multi-year mean AOD in HER (Huaihai Economic Region) between 2000 and 2016.The error bars represent the standard deviation of the mean.

Figure 6 .
Figure 6.Seasonal distribution of multi-year mean AOD in HER (Huaihai Economic Region) between 2000 and 2016.The error bars represent the standard deviation of the mean.

Figure 7 .
Figure 7. Average annual distribution map for AOD in HER during 2000-2016.

Figure 7 .
Figure 7. Average annual distribution map for AOD in HER during 2000-2016.

Figure 7 .
Figure 7. Average annual distribution map for AOD in HER during 2000-2016.

Figure 8 .
Figure 8.The topographic map of HER.

Figure 8 .
Figure 8.The topographic map of HER.

Figure 10 .
Figure 10.Significance level for AOD long-term change trend of spatial distribution over HER between 2000 and 2016.

Figure 10 .
Figure 10.Significance level for AOD long-term change trend of spatial distribution over HER between 2000 and 2016.

Figure 11 .
Figure 11.(A) Relationships between AOD and each driving force and (B) correlation coefficients in HER: AWS, average wind speed; SD, sunshine duration; ARH, average relative humidity; NDVI, normalized difference vegetation index; GDP, gross domestic product; PD, population density.The time range for the data used here is 2000-2015.

Figure 11 .
Figure 11.(A) Relationships between AOD and each driving force and (B) correlation coefficients in HER: AWS, average wind speed; SD, sunshine duration; ARH, average relative humidity; NDVI, normalized difference vegetation index; GDP, gross domestic product; PD, population density.The time range for the data used here is 2000-2015.

Table 1 .
The details of the sunphotometer observation sites.

Table 1 .
The details of the sunphotometer observation sites.

Table 2 .
The AOD growth rate in different year periods of HER (Huaihai Economic Region) (%).