Multi-Time Scale Analysis of Regional Aerosol Optical Depth Changes in National-Level Urban Agglomerations in China Using Modis Collection 6 . 1 Datasets from 2001 to 2017

With the rapid development of China’s economy and industry, characterizing the spatial and temporal changes of aerosols in China has attracted widespread attention from researchers. The national-level urban agglomerations are the most concentrated areas of China’s economic, population and resource. Studying the spatial and temporal changes of aerosol optical depth (AOD) in these regions has practical guiding significance for effective monitoring of atmospheric particulate pollution. This paper analyzed the spatial and temporal variations of AOD in China’s urban agglomerations during 2001–2017 by using Terra Moderate resolution Imaging Spectroradiometer (MODIS) Collection 6.1 (C6.1) Level 2 aerosol products (MOD04_L2). Five national-level urban agglomerations were chosen: Yangtze River Delta (YRD), Pearl River Delta (PRD), Beijing-Tianjin-Hebei (BTH), Yangtze River Middle-Reach (YRMR) and Cheng-Yu (CY). We analyzed the change patterns of AOD in different urban agglomerations at multi-time scales and built a time series decomposition model to mine the long-term trend, seasonal variation and abnormal change information of AOD time series. The result indicated that averaged AOD values in the five urban agglomerations were basically increased first and then decreased at the annual time scale during 2001–2017. The averaged AOD showed strong seasonal differences and AOD values in spring and summer were typically higher than those in autumn and winter. At the monthly time scale, the AOD typically varied from low in cold months to high in warm months and then decreased during the rainy periods. Time series decompositions revealed that a notable transition around 2007–2008 dominated the long-term overall trend over the five selected urban agglomerations and an initial upward tendency followed by a downward tendency was observed during 2001–2017. This study can be utilized to provide decision-making basis for atmospheric environmental governance and future development of urban agglomerations.


Introduction
Atmospheric aerosols are solid or liquid particles that are suspended in the air.Atmospheric aerosols have great spatial and temporal variability and their distribution characteristics have a significant influence on the formation of cloud and precipitation [1,2].Aerosols can reduce the amount of solar radiation and sunlight reaching the ground by absorbing and scattering solar radiation, thereby responding to the Earth's radiation balance [3].Additionally, aerosol particles can act as cloud condensation nuclei to affect the optical properties of the cloud, the amount of clouds and the lifetime of the cloud, thereby affecting precipitation and water cycles [4,5].Furthermore, aerosols can alter atmospheric chemical processes, which in turn affect the concentration and distribution of greenhouse gas [6,7].Therefore, the monitoring of the spatial and temporal distribution of aerosols has been used as an important basis for accurately assessing the climatic effects of aerosols [8,9] and the impacts of aerosols on air quality [10].
Aerosol optical depth (AOD), which is defined as the total attenuation of absorption and scattering of solar radiation by aerosols along the radiation transmission path, is one of the core parameters widely used for characterizing atmospheric aerosol properties [11] and evaluating their radiation radiative forcing and climatic effects [12].In addition, AOD has been used as a proxy for particulate pollution by the air quality community [13].Satellite-derived AOD products provide a unique tool to monitor aerosol properties at regional and global scales.Moderate resolution Imaging Spectroradiometer (MODIS) AOD is widely used to characterize the spatiotemporal change of atmospheric aerosols because it has been extensively validated and provides continuous observations extending for nearly two decades [14][15][16].Numerous studies have used MODIS AOD to study the long-term spatiotemporal variation of aerosols at regional and global scales [17][18][19][20][21][22][23][24][25].
As the world's largest developing country, China is experiencing rapid economic expansion, industrial development and urbanization, accompanied by increasing air pollution and aerosol particles are the major harmful pollutants that are responsible for the poor air quality of most urban areas in China.Several studies have quantitatively analyzed the spatial and temporal changes of aerosols in China.Xie and Xia (2008) analyzed the spatial and temporal variations of AOD in north China from 1980 to 2001 by using Total Ozone Mapping Spectrometer (TOMS) monthly AOD data [26].Guo et al. (2011) conducted a long-term trend analysis of AOD at one-degree over eight typical regions in China, both spatially and temporally, using TOMS AOD products , along with MODIS AOD data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) [27].Luo et al. (2014) analyzed the trends and seasonal variations of AOD over 10 regions in China from 2000 to 2010 by using the Level 3 MODIS AOD products [28].Che et al. (2015) characterized the aerosol optical properties over a large area of China based on an extensive set of ground-based observations collected from the China Aerosol Remote Sensing Network (CARSNET) during 2002-2013 [29].He et al. (2016) identified the spatio-temporal trends of AOD in China and quantified the plausible factors that influence AOD, using Aqua MODIS Collection 6 aerosol data obtained between 2002 and 2015 [30].Zhao et al. (2017) investigated inter-annual AOD trends during 2001-2015 over Eastern and Central China using observations from Multi-angle Imaging SpectroRadiometer (MISR) and MODIS [31].Sogacheva et al. (2018) analyzed the spatial patterns and inter-annual and seasonal variations of AOD over China based on one-degree Level 3 monthly AOD data retrieved from AlongTrack Scanning Radiometer (ATSR -2, 1995-2002) and Advanced ATSR (AATSR, 2002(AATSR, -2012) ) and Level 3 Terra MODIS merged Dark Target and Deep Blue (DTDB) AOD product (2000-2017) [24,25].
However, few studies have analyzed the distribution and dynamic characteristics of aerosols over national-level urban agglomerations in China, which are regions with the most developed economy and the highest population density.The five national-level urban agglomerations accounted for approximately half of China's GDP in 2010 and the population in these regions has been steadily increasing [25,32].The migration of people from rural areas to these large cities has resulted in fast urbanization and further industrialization and economic growth in these regions [25].Along with the rapid economic growth and expanding urbanization, these regions are also the most polluted areas in China.Therefore, analysis of the spatial and temporal changes of long-term AOD time series can lay a foundation for the study of atmospheric aerosols in these areas and has practical guiding significance for the sustainable development of urban agglomerations and the construction of an environment-friendly society.
In this study, the latest MODIS Collection 6.1 (C6.1) aerosol products are used to analyze the spatiotemporal variations of AOD over five national-level urban agglomerations in China.The analysis is performed at annual, seasonal and monthly time scales and a time series analysis model is built to explore the dynamic characteristics of the monthly AOD time series.The main goal of this study is to provide a reference for understanding the spatial and temporal change patterns of atmospheric aerosols over these highly urbanized areas.study is to provide a reference for understanding the spatial and temporal change patterns of atmospheric aerosols over these highly urbanized areas.

Study areas
Five national-level urban agglomerations in China are selected for this study: the Yangtze River Delta urban agglomeration (YRD), the Pearl River Delta urban agglomeration (PRD), the Beijing-Tianjin-Hebei urban agglomeration (BTH), the Yangtze River Middle-Reach urban agglomeration (YRMR) and the Cheng-Yu urban agglomeration (CY).The spatial pattern of these urban agglomerations was proposed in the National New Urbanization Plan (2014-2020) issued on March 16, 2013 by the Chinese Central Government [32].YRD sits along the central-eastern coastline of China.YRMR is located in central China and includes 3 medium-scale urban agglomerations: the Chang-Zhu-Tan urban agglomeration, the Wuhan metropolitan area and the Poyang Lake city group.YRD and YRMR have a subtropical monsoon climate characterized by hot, humid summers.PRD is located in southeastern China and has a humid subtropical climate with frequent cloudy and rainy weather.BTH is in the northern part of the North China Plain, bordering on the Bohai Sea.It has a monsoon-influenced humid continental climate with severe, dry winters.CY lies in the upstream of Yangtze River Basin in southwest China.It has a subtropical monsoon climate characterized by high humidity and many overcast days.The five selected urban agglomerations have stimulated the development of the national economy and have been the current foci of China's New Urbanization strategy [32,33].Their boundaries are indicated by different colors in Figure 1.

MODIS Collection 6.1 Aerosol Products
The MODIS instruments onboard polar orbiting and sun-synchronous Terra and Aqua satellites have 36 spectral channels and provide nearly daily global coverage at spatial resolutions of 250 m, 500 m and 1000 m [34].The MODIS Adaptive Processing System (MODAPS) generates global AOD products over land at 3 km or 10 km spatial resolutions based on the Dark Target (DT) algorithm and a 10 km spatial resolution based on the Deep Blue (DB) algorithm.It has recently released an improved C6.1 for all MODIS Level 1 (L1) and higher-level Level 2 (L2) & Level 3 (L3) Atmosphere Team products.
In this study, Terra MODIS C6.1 Level 2 aerosol products (MOD04_L2) at 10 km spatial resolution from 1 January 2001 to 31 December 2017 were used.All data were downloaded from the Level 1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC) at NASA's Goddard Space Flight Center in Greenbelt, Maryland, United States.The merged DTDB AOD data were extracted and only land pixels were used.We validated the DTDB AOD product over China with Aerosol Robotic Network (AERONET) data for the sites with an observation period longer than six months.A total of 28 sites were selected and AERONET Version 3 Level 2.0 AOD data were used as reference.Satellite data within a 25 km radius around the AERONET site were averaged to matchup with a one-hour AERONET data segment centered on the satellite overpass time [35].Figure 2 S1 in the Supplement.
The MODIS instruments onboard polar orbiting and sun-synchronous Terra and Aqua satellites have 36 spectral channels and provide nearly daily global coverage at spatial resolutions of 250 m, 500 m and 1000 m [34].The MODIS Adaptive Processing System (MODAPS) generates global AOD products over land at 3 km or 10 km spatial resolutions based on the Dark Target (DT) algorithm and a 10 km spatial resolution based on the Deep Blue (DB) algorithm.It has recently released an improved C6.1 for all MODIS Level 1 (L1) and higher-level Level 2 (L2) & Level 3 (L3) Atmosphere Team products.
In this study, Terra MODIS C6.1 Level 2 aerosol products (MOD04_L2) at 10 km spatial resolution from 1 January 2001 to 31 December 2017 were used.All data were downloaded from the Level 1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC) at NASA's Goddard Space Flight Center in Greenbelt, Maryland, United States.The merged DTDB AOD data were extracted and only land pixels were used.We validated the DTDB AOD product over China with Aerosol Robotic Network (AERONET) data for the sites with an observation period longer than six months.A total of 28 sites were selected and AERONET Version 3 Level 2.0 AOD data were used as reference.Satellite data within a 25 km radius around the AERONET site were averaged to matchup with a one-hour AERONET data segment centered on the satellite overpass time [35].Figure 2

Calculation of mean AOD distributions and regional time series
To analyze the spatial and temporal variation of AOD at annual, seasonal and monthly time scales, the annual, multi-year seasonal and multi-year monthly mean AOD distributions over China are calculated using a direct averaging of all of the valid retrievals.A minimum sampling frequency of 3% for each pixel is applied to exclude poorly sampled pixels.

Calculation of Mean AOD Distributions and Regional Time Series
To analyze the spatial and temporal variation of AOD at annual, seasonal and monthly time scales, the annual, multi-year seasonal and multi-year monthly mean AOD distributions over China are calculated using a direct averaging of all of the valid retrievals.A minimum sampling frequency of 3% for each pixel is applied to exclude poorly sampled pixels.
The regional mean AOD values for the five selected urban agglomerations are aggregated based on their boundaries.Since MODIS AOD values are retrieved only in cloud-free conditions and the derived AOD suffers from non-uniform and incomplete sampling issues, it is expected that a day with more valid pixels will enable a better estimate of the mean than a poorly sampled day.As a result, the area-weighted mean is calculated to reduce the impact of poorly sampled days.That is, each day's contribution is weighted by the number of valid pixels.The generic regional mean AOD (AOD mean ) is: where AOD i is the daily regional mean AOD for the day i, N i is the number of valid pixels on day i.

Time Series Decomposition
In order to find the intrinsic statistical characteristics of AOD variations over the five selected national-level urban agglomerations, time series decomposition is implemented to decompose monthly AOD time series.We build an additive decomposition model to separate monthly AOD time series into three distinct components: where A t is the AOD value at time t, T t is the trend component, S t is the seasonal component and I t represents the irregular component.The trend component is the trend of long-term changes in the AOD time series that gradually increased or decreased with time.It is a trend curve that reflects the general direction of change of the monthly AOD time series.The seasonal component is the fixed variation exhibited by the AOD time series within one year and the main reason for the seasonal changes is the impact of the natural seasonal alternation.Irregular component reflects random changes and it is unpredictable.Irregular changes may occur due to natural disasters, human activities and sudden changes in climate.T t is estimated by using a quadratic trend model: where α 0 , α 1 and α 2 are regression coefficients.Then T t is subtracted from the original data A t to obtain the de-trended series D t .S t is modeled as a function of seasonal indicator variables: where γ i is the coefficient associated with the seasonal indicator variable and δ it is an indicator variable taking the value 1 in the ith month of the year and 0 in all other months.The de-trended series D t is regressed against the seasonal indicators to estimate the coefficients of the indicator variables and then the seasonal component S t can be calculated.The irregular component I t is estimated by subtracting T t and S t from A t .

Annual AOD Changes in Urban Agglomerations
The spatial distributions of annual averaged AOD in China during 2001-2017 are presented in Figure 3.It reveals that the spatial pattern of annual averaged AOD has obvious regional characteristics and is strongly influenced by topography.Since the topography of China is high in the west and low in the east (Figure 1), it disrupts the transmission of high concentrations of aerosol particles from the east to the west.In addition, the eastern part of China has many large industrial cities which are gathering centers for energy consumption and anthropogenic emission.These areas are highly prone to form areas with high AOD values.In YRMR, large industrial cities such as Wuhan, Nanchang and Changsha are areas with high AOD values.This is because the diffusion of pollutants is hindered by the surrounding mountains and the hygroscopic growth of aerosols is enhanced by excessive humidity caused by intensive lakes and rivers.As an important platform for the great development of western China, CY is the strategic support of the Yangtze River Economic Belt.However, the relatively occluded geographical environment of the Sichuan Basin makes it difficult for the aerosol particles generated in CY to be transported outside the basin and aerosol particles aggregate to some extent due to the foggy weather.As a result, the observed AOD is high in this region.In addition, Figures 3 and 4 show that AOD typically increases before 2006 and fluctuates during 2006-2011 and then decreases after 2011 over the selected urban agglomerations.The upward trend observed before 2006 is mainly due to emission increases caused by rapid economic and industrial development.In the 11th five-year plan (2006-2010), China set targets to improve energy efficiency and reduce emissions of sulfur dioxide and primary aerosols.During the 12th five-year plan (2011-2015), the released environmental protection policies extended the previous policies and aimed to improve the air quality further [36].Due to the emission reductions, the trend reversal observed in 2011 might be attributed to the implementation of these environmental protection policies [37,38].In September 2013, China issued the Air Pollution Prevention and Control Action Plan, which aimed to reduce particulate matter concentrations across the country [36].The Action Plan strengthened the targets of air quality improvement and raised stringent air pollution control policies.It seems that the Action Plan has led to a significant drop of AOD since 2014.

Seasonal AOD changes in urban agglomerations
The spatial distributions of seasonal averaged AOD in China during 2001-2017 are presented in Figure 5.It reveals that the averaged AOD values in most parts of the five selected urban agglomerations in spring and summer are significantly higher than autumn and winter.The spatial distribution of areas with high and low AOD values of different urban agglomerations is basically unchanged in the four seasons but the intensity of these areas reflects their respective regional seasonal changes.
Figure 6 compares the seasonal averaged AOD values in five urban agglomerations during 2001-2017.For the five selected regions, AOD is highest in spring, except for BTH where the maximum AOD is observed in summer.In PRD and YRD, AOD values are similar in spring and summer.The high summer AOD over the two regions can be partially caused by the hygroscopic growth of aerosols under high relative humidity conditions.For PRD, YRD and YRMR, AOD is lowest in winter, whereas the minimum AOD is observed in autumn for BTH and CY.In BTH, intensive wintertime warming activities enhance the anthropogenic emission, thus the AOD increases in winter.For CY, the high winter AOD is mainly caused by temperature inversions [39].Since the Sichuan Basin is surrounded by mountains, the strong inversion layer above the atmospheric boundary layer will trap any pollutant that enters the air and therefore haze gradually builds up.

Seasonal AOD Changes in Urban Agglomerations
The spatial distributions of seasonal averaged AOD in China during 2001-2017 are presented in Figure 5.It reveals that the averaged AOD values in most parts of the five selected urban agglomerations in spring and summer are significantly higher than autumn and winter.The spatial distribution of areas with high and low AOD values of different urban agglomerations is basically unchanged in the four seasons but the intensity of these areas reflects their respective regional seasonal changes.
Figure 6 compares the seasonal averaged AOD values in five urban agglomerations during 2001-2017.For the five selected regions, AOD is highest in spring, except for BTH where the maximum AOD is observed in summer.In PRD and YRD, AOD values are similar in spring and summer.The high summer AOD over the two regions can be partially caused by the hygroscopic growth of aerosols under high relative humidity conditions.For PRD, YRD and YRMR, AOD is lowest in winter, whereas the minimum AOD is observed in autumn for BTH and CY.In BTH, intensive wintertime warming activities enhance the anthropogenic emission, thus the AOD increases in winter.For CY, the high winter AOD is mainly caused by temperature inversions [39].Since the Sichuan Basin is surrounded by mountains, the strong inversion layer above the atmospheric boundary layer will trap any pollutant that enters the air and therefore haze gradually builds up.

Monthly AOD changes in urban agglomerations
Figure 7 presents the spatial distributions of monthly averaged AOD in China during 2001-2017.The spatial distribution of areas with high and low AOD values shows similar patterns as the seasonal averaged AOD.However, monthly maps can provide useful information for qualifying the temporal variations of AOD due to different weather conditions.AOD values show an increasing trend during spring and early summer months due to the continuous increase of temperature and humidity.The monthly averaged AOD is prone to yielding large amounts of data gaps in June and July over PRD due to its frequent cloudy and rainy days.In addition, significant amounts of missing

Monthly AOD changes in urban agglomerations
Figure 7 presents the spatial distributions of monthly averaged AOD in China during 2001-2017.The spatial distribution of areas with high and low AOD values shows similar patterns as the seasonal averaged AOD.However, monthly maps can provide useful information for qualifying the temporal variations of AOD due to different weather conditions.AOD values show an increasing trend during spring and early summer months due to the continuous increase of temperature and humidity.The monthly averaged AOD is prone to yielding large amounts of data gaps in June and July over PRD due to its frequent cloudy and rainy days.In addition, significant amounts of missing  The spatial distribution of areas with high and low AOD values shows similar patterns as the seasonal averaged AOD.However, monthly maps can provide useful information for qualifying the temporal variations of AOD due to different weather conditions.AOD values show an increasing trend during spring and early summer months due to the continuous increase of temperature and humidity.The monthly averaged AOD is prone to yielding large amounts of data gaps in June and July over PRD due to its frequent cloudy and rainy days.In addition, significant amounts of missing AOD areas are observed in December and January over CY because of frequent fog events in the Sichuan basin.

Monthly AOD Changes in Urban Agglomerations
Figure 8 compares the monthly averaged AOD values in five urban agglomerations.Typically, AOD varies seasonally from low in cold months to high in warm months and then decreases during the rainy periods.In BTH, YRD and YRMD, the highest AOD is observed in June.This is because the progressive increase of air temperature from January to June will accelerate photochemical reactions and increase aerosol concentrations [40].In addition, the dust aerosol transport from the arid and semi-arid regions also contributes to the increase in aerosol loading [24].Although high relative humidity during summer months strengthens the local aerosol accumulation, East Asian summer monsoon rainfalls wash out aerosol particles and result in low AOD values over these regions.Variations of monthly averaged AOD are much more irregular over PRD and CY than over the other regions, the zigzag pattern may be caused by the sampling issues of satellite measurements due to the frequent foggy, cloudy or rainy days over the two regions.
AOD areas are observed in December and January over CY because of frequent fog events in the Sichuan basin.
Figure 8 compares the monthly averaged AOD values in five urban agglomerations.Typically, AOD varies seasonally from low in cold months to high in warm months and then decreases during the rainy periods.In BTH, YRD and YRMD, the highest AOD is observed in June.This is because the progressive increase of air temperature from January to June will accelerate photochemical reactions and increase aerosol concentrations [40].In addition, the dust aerosol transport from the arid and semi-arid regions also contributes to the increase in aerosol loading [24].Although high relative humidity during summer months strengthens the local aerosol accumulation, East Asian summer monsoon rainfalls wash out aerosol particles and result in low AOD values over these regions.Variations of monthly averaged AOD are much more irregular over PRD and CY than over the other regions, the zigzag pattern may be caused by the sampling issues of satellite measurements due to the frequent foggy, cloudy or rainy days over the two regions.The trend component of monthly averaged AOD time series in the five selected urban agglomerations indicates that AOD increases at first and then decreases between 2001 and 2017 and the transition occurs during 2007-2008.In November 2007, the State Council of China approved the 11 th five-year plan for environmental protection, which aimed to improve energy efficiency and reduce emissions of sulfur dioxide and primary aerosols [41].As a result, the observed trend reversal may be attributed to the implementation of environmental measures.The seasonal component typically has positive values in spring and summer and negative values in autumn and winter.This is in accordance with the fact that high AOD is observed in spring and summer, as is shown in Figure 5. Extreme values of the irregular component are typically observed during rainy months and this may be caused by the sampling issues of satellite measurements.That is, the   The trend component of monthly averaged AOD time series in the five selected urban agglomerations indicates that AOD increases at first and then decreases between 2001 and 2017 and the transition occurs during 2007-2008.In November 2007, the State Council of China approved the 11 th five-year plan for environmental protection, which aimed to improve energy efficiency and reduce emissions of sulfur dioxide and primary aerosols [41].As a result, the observed trend reversal may be attributed to the implementation of environmental measures.The seasonal component typically has positive values in spring and summer and negative values in autumn and winter.This is in accordance with the fact that high AOD is observed in spring and summer, as is shown in Figure 5. Extreme values of the irregular component are typically observed during rainy months and this may be caused by the sampling issues of satellite measurements.That is, the The trend component of monthly averaged AOD time series in the five selected urban agglomerations indicates that AOD increases at first and then decreases between 2001 and 2017 and the transition occurs during 2007-2008.In November 2007, the State Council of China approved the 11th five-year plan for environmental protection, which aimed to improve energy efficiency and reduce emissions of sulfur dioxide and primary aerosols [41].As a result, the observed trend reversal may be attributed to the implementation of environmental measures.The seasonal component typically has positive values in spring and summer and negative values in autumn and winter.This is in accordance with the fact that high AOD is observed in spring and summer, as is shown in Figure 5. Extreme values of the irregular component are typically observed during rainy months and this may be caused by the sampling issues of satellite measurements.That is, the regional mean AOD is biased because a large number of pixels are missed due to frequent cloudy and rainy days.

Discussion
Sogacheva et al. ( 2018) investigated the long-term AOD variations over China during 1995-2017 by combining AOD data retrieved from ATSR-2, AATSR and MODIS instruments.Linear regression was applied to estimate the AOD tendencies for two periods: 1995-2006 and 2011-2017 [25].Although the choice of study regions is different from this paper, it also covers major urban/industrial regions.Their results also revealed an increasing tendency before 2006 and a reverse tendency after 2011 in BTH, YRD, PRD and CY.
To calculate the regional mean AOD at different time scales, we weighted the daily regional mean AOD by the number of valid pixels, bearing in mind that a day with more valid pixels would enable a better estimate of the mean than a poorly sampled day.However, the calculated mean AOD might still be biased when large amounts of data are missing during the selected period.Figure 14 indicates that the sampling frequencies of AOD show significant seasonal changes.Typically, the sampling frequency is lower than 10% in summer over PRD and in winter over CY.This is because AOD data are only available during clear-sky conditions.Clouds, fog and precipitation will lead to a large amount of missing AOD data.Figure 15 presents the seasonal distributions of the standard error of the mean (SEM) derived from AOD data.SEM is used as an indirect measure of the reliability of the estimated mean.It depends on both the variability of AOD and the number of observations.For example, the cloudy fog in the Sichuan Basin in winter makes it difficult to obtain adequate amounts of effective observations and large SEM values indicate that mean AOD values have a low reliability over this region.That is, the derived mean AOD might be unrepresentative and the sampling issue of satellite-derived AOD products should be taken into account.
Remote Sens. 2019, 11 FOR PEER REVIEW 14 regional mean AOD is biased because a large number of pixels are missed due to frequent cloudy and rainy days.

Discussion
Sogacheva et al. ( 2018) investigated the long-term AOD variations over China during 1995-2017 by combining AOD data retrieved from ATSR-2, AATSR and MODIS instruments.Linear regression was applied to estimate the AOD tendencies for two periods: 1995-2006 and 2011-2017 [25].Although the choice of study regions is different from this paper, it also covers major urban/industrial regions.Their results also revealed an increasing tendency before 2006 and a reverse tendency after 2011 in BTH, YRD, PRD and CY.
To calculate the regional mean AOD at different time scales, we weighted the daily regional mean AOD by the number of valid pixels, bearing in mind that a day with more valid pixels would enable a better estimate of the mean than a poorly sampled day.However, the calculated mean AOD might still be biased when large amounts of data are missing during the selected period.Figure 14 indicates that the sampling frequencies of AOD show significant seasonal changes.Typically, the sampling frequency is lower than 10% in summer over PRD and in winter over CY.This is because AOD data are only available during clear-sky conditions.Clouds, fog and precipitation will lead to a large amount of missing AOD data.Figure 15 presents the seasonal distributions of the standard error of the mean (SEM) derived from AOD data.SEM is used as an indirect measure of the reliability of the estimated mean.It depends on both the variability of AOD and the number of observations.For example, the cloudy fog in the Sichuan Basin in winter makes it difficult to obtain adequate amounts of effective observations and large SEM values indicate that mean AOD values have a low reliability over this region.That is, the derived mean AOD might be unrepresentative and the sampling issue of satellite-derived AOD products should be taken into account.A time series decomposition model was used to analyze the intrinsic statistical characteristics of AOD dynamics in this study, whereas there are two issues still deserving to be further discussed.One issue is about the premise of time series decomposition.The time series decomposition model assumes that the seasonal factors are the same for each year, whereas climate change can lead to different seasonal factors in different year, such as changes in clouds and water vapor, wind strength and the amount of sand and dust.This study does not address these seasonal changes and the impact of seasonal factors on aerosols.Another issue is associated with outliers.Although time series decomposition can accurately locate outliers, it does not handle outliers well.Further efforts should be made to eliminate the effects of outliers in the time series analysis model.

Conclusions
Using the latest MODIS Collection 6.1 (C6.1)Level 2 aerosol products, this paper analyzed the change patterns of AOD over five national-level urban agglomerations in China at multi-time scales and built a time series decomposition model to study the long-term dynamics of AOD in the five national-level urban agglomerations in China from 2001 to 2017.The spatial distribution of AOD in China had obvious regional characteristics and high AOD regions were mainly concentrated in east China where most urban agglomerations located.During the study period, changes of annual averaged AOD values in five urban agglomerations were basically increased first and then decreased.Time series decompositions also indicated that long-term AOD trends in the five selected urban agglomerations were basically rising first and then falling and the transition occurred during 2007-2008.The AOD of the five selected urban agglomerations was higher in spring and summer than that in autumn and winter.Although the seasonal variations of AOD in these urban agglomerations were different, the AOD typically varied from low in cold months to A time series decomposition model was used to analyze the intrinsic statistical characteristics of AOD dynamics in this study, whereas there are two issues still deserving to be further discussed.One issue is about the premise of time series decomposition.The time series decomposition model assumes that the seasonal factors are the same for each year, whereas climate change can lead to different seasonal factors in different year, such as changes in clouds and water vapor, wind strength and the amount of sand and dust.This study does not address these seasonal changes and the impact of seasonal factors on aerosols.Another issue is associated with outliers.Although time series decomposition can accurately locate outliers, it does not handle outliers well.Further efforts should be made to eliminate the effects of outliers in the time series analysis model.

Conclusions
Using the latest MODIS Collection 6.1 (C6.1)Level 2 aerosol products, this paper analyzed the change patterns of AOD over five national-level urban agglomerations in China at multi-time scales and built a time series decomposition model to study the long-term dynamics of AOD in the five national-level urban agglomerations in China from 2001 to 2017.The spatial distribution of AOD in China had obvious regional characteristics and high AOD regions were mainly concentrated in east China where most urban agglomerations located.During the study period, changes of annual averaged AOD values in five urban agglomerations were basically increased first and then decreased.Time series decompositions also indicated that long-term AOD trends in the five selected urban agglomerations were basically rising first and then falling and the transition occurred during 2007-2008.The AOD of the five selected urban agglomerations was higher in spring and summer than that in autumn and winter.Although the seasonal variations of AOD in these urban agglomerations were different, the AOD typically varied from low in cold months to high in warm months and then decreased during the rainy periods.The detailed analysis of trends in AODs over different urban agglomerations could

Five
national-level urban agglomerations in China are selected for this study: the Yangtze River Delta urban agglomeration (YRD), the Pearl River Delta urban agglomeration (PRD), the Beijing-Tianjin-Hebei urban agglomeration (BTH), the Yangtze River Middle-Reach urban agglomeration (YRMR) and the Cheng-Yu urban agglomeration (CY).The spatial pattern of these urban agglomerations was proposed in the National New Urbanization Plan (2014-2020) issued on March 16, 2013 by the Chinese Central Government [32].YRD sits along the central-eastern coastline of China.YRMR is located in central China and includes 3 medium-scale urban agglomerations: the Chang-Zhu-Tan urban agglomeration, the Wuhan metropolitan area and the Poyang Lake city group.YRD and YRMR have a subtropical monsoon climate characterized by hot, humid summers.PRD is located in southeastern China and has a humid subtropical climate with frequent cloudy and rainy weather.BTH is in the northern part of the North China Plain, bordering on the Bohai Sea.It has a monsoon-influenced humid continental climate with severe, dry winters.CY lies in the upstream of Yangtze River Basin in southwest China.It has a subtropical monsoon climate characterized by high humidity and many overcast days.The five selected urban agglomerations have stimulated the development of the national economy and have been the current foci of China's New Urbanization strategy [32,33].Their boundaries are indicated by different colors in Figure 1.Remote Sens. 2019, 11 FOR PEER REVIEW 3

Figure 1 .
Figure 1.Locations of five selected urban agglomerations in China.

Figure 1 .
Figure 1.Locations of five selected urban agglomerations in China.
depicts a good agreement between the MODIS and AERONET AOD values.The coefficient of determination (R 2 ) is 0.85 and the root mean squared error (RMSE) is 0.21.The validation results indicate that MODIS C6.1 DTDB AOD is suitable for characterizing long-term AOD variations over China.Statistics of the site-level MODIS/AERONET collocations are listed in Table depicts a good agreement between the MODIS and AERONET AOD values.The coefficient of determination (R 2 ) is 0.85 and the root mean squared error (RMSE) is 0.21.The validation results indicate that MODIS C6.1 DTDB AOD is suitable for characterizing long-term AOD variations over China.Statistics of the site-level MODIS/AERONET collocations are listed in TableS1in the Supplement.

Figure 2 .
Figure 2. Density scatterplot of MODIS C6.1 merged Dark Target and Deep Blue (DTDB) AOD vs. AOD from AERONET stations in China during 2001-2017.A MODIS-AERONET collocation is the spatial average of MODIS AOD within a 25 km radius around the AERONET site and the temporal average of all AERONET observations within ±30 min of satellite overpass.

Figure 2 .
Figure 2. Density scatterplot of MODIS C6.1 merged Dark Target and Deep Blue (DTDB) AOD vs. AOD from AERONET stations in China during 2001-2017.A MODIS-AERONET collocation is the spatial average of MODIS AOD within a 25 km radius around the AERONET site and the temporal average of all AERONET observations within ±30 min of satellite overpass.

Figure 4
indicates that the temporal variations of annual averaged AOD values in five urban agglomerations are quite different.The annual averaged AOD values of YRD vary from 0.49 to 0.76 during 2001-2017.The change of AOD values in YRD shows a growth trend from 2001 to 2007, while a downward trend is observed during 2011-2017.A similar pattern is observed in YRMR, except the initial upward trend ends at 2006 and AOD values in this region are lower than that of YRD between 2007 and 2017.The annual averaged AOD values of CY vary from 0.39 to 0.83.The peak AOD value is observed in 2006 and a continuous downward trend is observed from 2011 onward.The annual averaged AOD of PRD is observed of the order of 0.5-0.7 during 2001-2017.The AOD in PRD gradually decreases after 2012 and the minimum value of 0.35 is observed in 2016.The annual averaged AOD values of BTH are lower than other regions.It is because the western and northern parts of this region are the Xishan and Yanshan mountain ranges, which are low AOD areas.Therefore, regional mean AOD values are biased toward low values despite the AOD is high over the southern part of this region.

Figure 3 .
Figure 3.The annual AOD changes over five urban agglomerations in China.

Figure 3 .
Figure 3.The annual AOD changes over five urban agglomerations in China.

Figure 4 .
Figure 4. Annual averaged AOD values in five urban agglomerations.

Figure 4 .
Figure 4. Annual averaged AOD values in five urban agglomerations.

Figure 5 .
Figure 5.The seasonal AOD changes over five urban agglomerations in China.

Figure 6 .
Figure 6.Seasonal averaged AOD values in five urban agglomerations.

Figure 5 .
Figure 5.The seasonal AOD changes over five urban agglomerations in China.

Figure 5 .
Figure 5.The seasonal AOD changes over five urban agglomerations in China.

Figure 6 .
Figure 6.Seasonal averaged AOD values in five urban agglomerations.

Figure 6 .
Figure 6.Seasonal averaged AOD values in five urban agglomerations.

Figure 7
Figure 7 presents the spatial distributions of monthly averaged AOD in China during 2001-2017.The spatial distribution of areas with high and low AOD values shows similar patterns as the seasonal averaged AOD.However, monthly maps can provide useful information for qualifying the temporal variations of AOD due to different weather conditions.AOD values show an increasing trend during spring and early summer months due to the continuous increase of temperature and humidity.The monthly averaged AOD is prone to yielding large amounts of data gaps in June and July over PRD due to its frequent cloudy and rainy days.In addition, significant amounts of missing AOD areas are observed in December and January over CY because of frequent fog events in the Sichuan basin.

Figure 7 .
Figure 7.The monthly AOD changes over five urban agglomerations in China.

Figure 7 .
Figure 7.The monthly AOD changes over five urban agglomerations in China.

Figure 8 .
Figure 8. Monthly averaged AOD values in five urban agglomerations.
Time series decompositions of monthly AOD time series in the five selected urban agglomeration are shown in Figures 9-13.In these figures, SC represents the seasonal component and IC represents the irregular component.In YRD, the trend component of AOD time series fluctuates between 0.51 and 0.70.AOD is in the rising phase during 2001-2008 and then it decreases.The maximum value of the seasonal component is 0.35 occurred in June and the minimum seasonal component is −0.19 occurred in December.

Figure 9 .
Figure 9.Time series decomposition of monthly AOD time series in the Yangtze River Delta urban agglomeration.

Figure 8 .
Figure 8. Monthly averaged AOD values in five urban agglomerations.

3. 2 .
AOD Changes Using Time Series Decomposition Time series decompositions of monthly AOD time series in the five selected urban agglomeration are shown in Figures 9-13.In these figures, SC represents the seasonal component and IC represents the irregular component.In YRD, the trend component of AOD time series fluctuates between 0.51 and 0.70.AOD is in the rising phase during 2001-2008 and then it decreases.The maximum value of the seasonal component is 0.35 occurred in June and the minimum seasonal component is −0.19 occurred in December.

Figure 7 .
Figure 7.The monthly AOD changes over five urban agglomerations in China.

Figure 8 .
Figure 8. Monthly averaged AOD values in five urban agglomerations.
Time series decompositions of monthly AOD time series in the five selected urban agglomeration are shown in Figures 9-13.In these figures, SC represents the seasonal component and IC represents the irregular component.In YRD, the trend component of AOD time series fluctuates between 0.51 and 0.70.AOD is in the rising phase during 2001-2008 and then it decreases.The maximum value of the seasonal component is 0.35 occurred in June and the minimum seasonal component is −0.19 occurred in December.

Figure 9 .
Figure 9.Time series decomposition of monthly AOD time series in the Yangtze River Delta urban agglomeration.

Figure 9 .
Figure 9.Time series decomposition of monthly AOD time series in the Yangtze River Delta urban agglomeration.In PRD, the trend component of AOD time series fluctuates between 0.41 and 0.67.AOD is in the rising stage between 2001 and 2007 and then a decline phase is observed.The maximum value of the seasonal component is 0.19 occurred in April and the minimum seasonal component is −0.16 occurred in November.

Figure 10 .
Figure 10.Time series decomposition of monthly AOD time series in the Pearl River Delta urban agglomeration.

Figure 11 .
Figure 11.Time series decomposition of monthly AOD time series in the Beijing-Tianjin-Hebei urban agglomeration.

Figure 10 .
Figure 10.Time series decomposition of monthly AOD time series in the Pearl River Delta urban agglomeration.In BTH, the trend component of AOD time series fluctuates between 0.39 and 0.52.AOD is in the rising stage during 2001-2008 and then it decreases.The maximum value of the seasonal component is 0.24 occurred in June and the minimum seasonal component is −0.21 occurred in November.

Figure 10 .
Figure 10.Time series decomposition of monthly AOD time series in the Pearl River Delta urban agglomeration.

Figure 11 .
Figure 11.Time series decomposition of monthly AOD time series in the Beijing-Tianjin-Hebei urban agglomeration.

Figure 11 .
Figure 11.Time series decomposition of monthly AOD time series in the Beijing-Tianjin-Hebei urban agglomeration.In YRMR, the trend component of AOD time series fluctuates between 0.41 and 0.68.AOD is in a rising phase between 2001 and 2007 and then decreases in the period 2008-2017.The maximum value of the seasonal component is 0.21 occurred in June and the minimum seasonal component is −0.17occurred in December.

Figure 12 .
Figure 12.Time series decomposition of monthly AOD time series in the Yangtze River Middle-Reach urban agglomeration.

Figure 13 .
Figure 13.Time series decomposition of monthly AOD time series in the Cheng-Yu urban agglomeration.

Figure 12 .
Figure 12.Time series decomposition of monthly AOD time series in the Yangtze River Middle-Reach urban agglomeration.In CY, the trend component of AOD time series fluctuates between 0.37 and 0.69.AOD is in rising phase during 2001-2007 and then a decline phase is observed.The maximum value of the seasonal component is 0.18 occurred in March and the minimum seasonal component is −0.20 occurred in December.

Figure 12 .
Figure 12.Time series decomposition of monthly AOD time series in the Yangtze River Middle-Reach urban agglomeration.

Figure 13 .
Figure 13.Time series decomposition of monthly AOD time series in the Cheng-Yu urban agglomeration.

Figure 13 .
Figure 13.Time series decomposition of monthly AOD time series in the Cheng-Yu urban agglomeration.

Figure 14 .
Figure 14.Seasonal distributions of AOD sampling frequency over five urban agglomerations in China.

Figure 14 .
Figure 14.Seasonal distributions of AOD sampling frequency over five urban agglomerations in China.

Figure 15 .
Figure 15.Seasonal distributions of the standard error of the mean (SEM) over five urban agglomerations in China.

Figure 15 .
Figure 15.Seasonal distributions of the standard error of the mean (SEM) over five urban agglomerations in China.