The Influence of Underlying Land Cover on the Accuracy of MODIS C6.1 Aerosol Products—A Case Study over the Yangtze River Delta Region of China

Due to the significant spatial variation of the performance of Moderate Resolution Imaging Spectrometer (MODIS) aerosol optical depth (AOD) retrievals, validation is very important for applications of MODIS AOD products at regional scales. This study presents a comparative analysis of Collection 6.1 MODIS AOD retrievals and ground measurements from five local sites and one Aerosol Robotic Network (AERONET) site in the Yangtze River Delta (YRD) region, which significantly complements the previous validation that utilized limited AERONET measurements. Generally, MODIS AOD retrievals showed a reasonable agreement with collocated ground measurements (R2 > 0.7), with 66% of Dark Target (DT) 10 km retrievals, 56% of Deep Blue (DB) 10 km retrievals, and 69% of DT 3 km retrievals falling within the expected error (EE = ±(0.05 + 0.2 × AOD)). Nevertheless, it was found that the DT AOD retrievals tended to be overestimated over urbanized and lakeside sites, while the DB AOD retrievals tended to be underestimated over all ground sites except for lakeside sites. Such patterns appeared to be linked with the systematic biases of the single-scattering albedo estimation in the AOD retrieval algorithms. Another significant finding of this study is that the uncertainties of the MODIS AOD retrievals were highly correlated with the land cover proportions of urbanized features and water (LCP_UW) in the surrounding region, especially for the DT products. An empirical correction method based on these correlations could substantially reduce the uncertainties of DT AOD products over high LCP_UW areas. The results not only highlight the significant impacts of both urban and water areas on the MODIS AOD retrieval algorithms but also create new possibilities to correct such impacts once the universal correlations between LCP_UW and the uncertainty measures are established.


Introduction
Atmospheric aerosols, which are defined as fine solid or liquid particles suspended in the atmosphere, play an important role in the Earth-atmosphere system through scattering and absorbing solar radiation (known as the aerosol direct effect) and modifying the microphysical properties and lifetimes of clouds (known as the aerosol indirect effect) [1]. Aerosol is also closely connected to many significant issues, including climate change, air quality, and public health [2][3][4]. However, owing to the wide spatial and temporal variability, complicated physical and chemical compositions, and short life cycle, atmospheric aerosol continues to be a large source of uncertainty in the climate system [5]. Although satellite It is technically challenging to understand the impacts of land cover on the accuracy of satellite-derived AOD because the uncertainties of the final AOD retrievals may not only be attributed to the land surface changes but also depend on the selection of aerosol type [40]. The MODIS DT/DB AOD algorithms simplify the aerosol conditions into five types (continental, weak absorption fine, moderate absorption fine, strong absorption fine, and dust coarse), where the uncertainties may also change between different assumed aerosol types [44]. An added difficulty is related to the fact that aerosol properties (scattering, particle size distributions, etc.) could also demonstrate remarkable differences even when the simplified aerosol type is the same [45], prohibiting quantitative estimates of the correlations between AOD uncertainties and land surface features.
As one of the most important economic centers in China, the Yangtze River Delta (YRD) has suffered from deteriorated air pollution issues and the rapid growth of urbanization and industrialization in recent years [46]. Investigating the performance of satellite-derived aerosol products over the YRD is conducive to understanding the spatiotemporal status of the air quality in the region. However, although several AERONET sites have been established in the YRD region over the past years, most of the sites provide only a small amount of quality-assured observations for a few days, prohibiting their uses in comprehensive validation efforts. To overcome this limitation, five ground sites over the Zhejiang Province region have been constructed by the local meteorological administrations to routinely monitor the air quality. Continuous aerosol measurements from these sites could not only provide air quality data for the surrounding areas but also support the assessments of MODIS AOD products in this region. Moreover, the spatial heterogeneity of the aerosol conditions could be considered limited between different sites, as they are located within a relatively small area. Therefore, the changes in the uncertainty levels of MODIS AOD retrievals among these ground sites appear to be linked to only the land cover dynamics, given the similar aerosol properties within a relatively small area.
Motivated by the lack of effective validations of MODIS AOD products in the YRD and the limited knowledge of the land-cover-associated AOD errors, the current study is designed to (1) evaluate the performance of Terra/Aqua MODIS DB10, DT3, and DT10 over the YRD region; (2) determine the correlations between the uncertainty measures of MODIS AOD and the percentages of urban and water areas in the surrounding region; and (3) discuss other potential factors that may impact the accuracy levels of MODIS AOD products.

In Situ AOD Data from Six Ground Sites
In situ measurements collected from six ground sites were used to validate MODIS AOD products in this study. These sites are mainly located in the core regions of the YRD (see locations in Figure 1), and the sites are named Hangzhou (HZ), Linan (LA), Tonglu (TL), Chunan (CA), Jiande (JD), and Taihu (TH). A Cimel-318 (CE-318) sun photometer was mounted at each location to measure the aerosol optical properties. Typically, the total uncertainty of AOD derived from CE-318 measurements is within ±0.01 for wavelengths longer than 440 nm [47]. Five of the six sites (i.e., HZ, LA, TL, CA, and JD) located in Zhejiang Province were maintained by the local meteorological administrations, and the measurements were obtained from 2013 to 2016. Level 1.0 AOD data (without cloud screening) over the five sites were calculated by ASTPwin software (Cimel Ltd. Co., Paris, France). Level 1.5 AOD data (cloud-screened and quality-controlled) were derived following the methodology discussed in detail by Smirnov et al. (2000) [48]. The number of days with valid AOD measurements in the five local sites ranged from 665 to 726 (Shown in Table 1). An AERONET site (TH) located in the neighboring Jiangsu Province was also included in the study. To guarantee the data quality, we chose the Version 3 database from the AERONET website and used the level 2.0 (quality-assured) inversion data in the validation. The number of days with valid AOD measurements during 2013-2016 in TH site was 316, just nearly a half in comparison with the other five sites. All ground sites presented a similar AOD frequency distribution (shown in Figure 2), and the average AOD over the six sites ranged from 0.51 to 0.67 (Table 1). Detailed statistical information on the datasets of the six sites is listed in Table 1 Taihu (TH), and they are all marked using two black boxes at sizes of 9 km × 9 km and 30 km × 30 km centered over the ground site, corresponding to sampling windows of 3 × 3 pixels for MODIS AOD products at spatial resolutions of 3 and 10 km.  Taihu (TH), and they are all marked using two black boxes at sizes of 9 km × 9 km and 30 km × 30 km centered over the ground site, corresponding to sampling windows of 3 × 3 pixels for MODIS AOD products at spatial resolutions of 3 and 10 km.

MODIS AOD Data
MODIS AOD products that were concurrent with the in situ measurements were downloaded from the Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/search/, accessed on 10 June 2021). Three types of MODIS AOD products (at 0.55 μm) are available for MODIS instruments aboard both Terra and Aqua, including DT10 and DB10 products  The land cover maps over the six ground sites during 2010 were obtained from the GlobeLand30 product (http://www.globallandcover.com, accessed on 21 October 2020), which has a spatial resolution of 30 m and an overall reported accuracy above 80% [49]. Five land cover types surround the ground sites: artificial surface (or urbanized area), forest, grassland, water, and cultivated land. Specifically, the land cover of four of the stations (LA, TL, CA, and JD) is dominated by vegetated areas, accounting for 50% of both the 9 km × 9 km and 30 km × 30 km boxes. In contrast, the HZ site is located in an urbanized area, and the TH site is close to water, and artificial surface features and water account for the most coverage in the selected windows centered over these stations. Indeed, the land cover proportions of urbanized features and water (LCP_UW) represent 10~80% of the total coverage within the windows, and such a large dynamic range makes it possible to analyze the impacts of LCP_UW on the accuracy of MODIS AOD retrievals.

MODIS AOD Data
MODIS AOD products that were concurrent with the in situ measurements were downloaded from the Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/search/, accessed on 10 June 2021). Three types of MODIS AOD products (at 0.55 µm) are available for MODIS instruments aboard both Terra and Aqua, including DT10 and DB10 products (MOD04_2L/MYD04_2L) and DT3 products (MOD04_3K/MYD04_3K). The MODIS DT algorithm takes advantage of the linear relationship of the spectral reflectance at the 0.47, 0.66, and 2.12 µm bands to retrieve the AOD over most vegetated and dark soiled surfaces [50]. Data with two resolutions (3 and 10 km) are provided by the DT algorithm. The DT3 AOD product differs from the 10 km product in only the size of the window during the generation processes. Validations with data collected from the global AERONET sites showed that the uncertainties of the 3 and 10 km DT AOD products are within ±(0.05 + 0.20τ) and ±(0.05 + 0.15τ) [7,30], respectively. Meanwhile, due to the adoption of a hybrid approach, the C6.1 MODIS DB algorithm could retrieve AOD over entire land surfaces including urbanized and vegetated regions [6]. The uncertainty of the DB AOD product is estimated to be within ±(0.03 + 0.20τ) [29]. Indeed, we selected only the C6.1 MODIS AOD products with the highest quality (quality assurance flag = 3) for validation in this study.

Selection of the Match-Ups between Satellite and Ground Measurements
Allowing for that both spatial and temporal matching strategy would have a significant influence on the accuracy assessment of satellite-derived AOD products, we adopted rigorous match-up criteria to minimize the impact of sampling mode on the accuracy assessment results. (1) Spatially, sampling windows of 3 × 3 pixels centered on the ground site were first extracted from the MODIS AOD product, and the windows with less than 2 pixels of valid AOD retrievals were excluded. A heterogeneity test was also applied to assure that no significant aerosol gradient occurred within each box, and we selected only the windows with associated standard deviations of less than 0.01. Then, the mean values of the sampling windows were calculated and used to compare with the ground-measured AOD. (2) Temporally, the difference in observation time between satellite and ground data was limited within ±30 min, assuming that the changes in aerosol conditions should be negligible within such a short period.

Accuracy Measures
The accuracy measures used to evaluate the MODIS AOD retrievals include the square of the correlation coefficient (R 2 ), root mean square error (RMSE), mean relative error (MRE), and expected error (EE). They are defined as follows: where AOD M and AOD G represent the collocated MODIS and ground-derived AOD retrievals and AOD M and AOD G represent the means of the MODIS and in situ AOD retrievals, respectively. The EE defined in this study is the target accuracy of the MODIS 3 km AOD product [30], which is the least restrictive among the three MODIS AOD products. The percentages of MODIS AOD retrievals that fall above, within, and below the EE envelope are calculated by Equations (7)- (9).

Correlations between MODIS AOD Accuracies and Land Cover Proportions
For each ground site, we calculated the coverage percentages for different cover types within two different boxes (i.e., 9 km × 9 km and 30 km × 30 km) centered over the site ( Figure 3). The sizes of the two selected boxes correspond to 3 × 3-pixel windows in 3 and 10 km MODIS AOD products. To quantitatively estimate how the surface features could impact the performance of the MODIS AOD algorithms, we analyzed the relationships between the percentage of urban and water areas within a region and the accuracy measures of the MODIS AOD retrievals (in terms of R 2 , within EE, and MRE) over the sampling windows centered over the ground sites. The proportion of urban cover has been applied in some DT-style inversion algorithms as a correction factor to improve the accuracy of Remote Sens. 2022, 14, 938 7 of 18 MODIS AOD retrievals over urban areas [38,51]. However, the presence of water could also serve as an important factor to interrupt the assumed correlations between visible and SWIR bands, increasing the bias of MODIS AOD retrievals [24]. Therefore, the coverage percentages of both water and urban areas were considered in this study.
measures of the MODIS AOD retrievals (in terms of R 2 , within EE, and MRE) over the sampling windows centered over the ground sites. The proportion of urban cover has been applied in some DT-style inversion algorithms as a correction factor to improve the accuracy of MODIS AOD retrievals over urban areas [38,51]. However, the presence of water could also serve as an important factor to interrupt the assumed correlations between visible and SWIR bands, increasing the bias of MODIS AOD retrievals [24]. Therefore, the coverage percentages of both water and urban areas were considered in this study.

Overall Accuracy of MODIS C6.1 AOD Product
The comparisons between three types of MODIS Terra AOD retrievals and concurrent ground site measurements are illustrated in Figures 4-6. The performances of the satellite-derived AOD products varied significantly among the different ground sites and algorithms. Generally, when the land cover was dominated by vegetation, pronounced agreements could be found between the satellite and in situ data. Numerically, the correlation coefficients (R 2 ) for the four vegetated sites were 0.77~0.81, 0.8~0.86, and 0.62~0.72 for the DT10, DT3, and DB10 products, respectively. Moreover, the fractions of the AOD retrievals that fell within the EE were 71~83%, 63~87%, and 49~60% for the DT10, DT3, and DB10 products, respectively. Slightly lower accuracies were found for urbanized and lakeside sites (HZ and TH). The R 2 values for these two sites were 0.39~0.71 for DT10, 0.28~0.59 for DT3, and 0.62~0.76 for DB10, and the fractions of the AOD retrievals within the EE were 17~50% for DT10, 24~27% for DT3, and 61~62% for DB10. On the other hand, while the DT algorithm (including 3 and 10 km products) outperformed the DB algorithm over the vegetated region, the opposite was found to be true for urbanized and lakeside sites. Such comparisons were also conducted for the MODIS Aqua AOD products ( Figures  S1-S3), where the patterns appeared to be very similar to those of the MODIS Terra products.

Overall Accuracy of MODIS C6.1 AOD Product
The comparisons between three types of MODIS Terra AOD retrievals and concurrent ground site measurements are illustrated in Figures 4-6. The performances of the satellite-derived AOD products varied significantly among the different ground sites and algorithms. Generally, when the land cover was dominated by vegetation, pronounced agreements could be found between the satellite and in situ data. Numerically, the correlation coefficients (R 2 ) for the four vegetated sites were 0.77~0.81, 0.8~0.86, and 0.62~0.72 for the DT10, DT3, and DB10 products, respectively. Moreover, the fractions of the AOD retrievals that fell within the EE were 71~83%, 63~87%, and 49~60% for the DT10, DT3, and DB10 products, respectively. Slightly lower accuracies were found for urbanized and lakeside sites (HZ and TH). The R 2 values for these two sites were 0.39~0.71 for DT10, 0.28~0.59 for DT3, and 0.62~0.76 for DB10, and the fractions of the AOD retrievals within the EE were 17~50% for DT10, 24~27% for DT3, and 61~62% for DB10. On the other hand, while the DT algorithm (including 3 and 10 km products) outperformed the DB algorithm over the vegetated region, the opposite was found to be true for urbanized and lakeside sites. Such comparisons were also conducted for the MODIS Aqua AOD products ( Figures S1-S3), where the patterns appeared to be very similar to those of the MODIS Terra products.

Correlations between Accuracy Levels and the Proportions of Water and Urban Areas
The accuracy measures (including R 2 , MRE, and percentage of data that fell within the EE (PWEE)) of the different MODIS AOD products are plotted against LCP_UW in Figure 7, and the corresponding correlations are listed in Table 2. Significant relationships between LCP_UW and R 2 were only found in all DT10 AOD retrievals, indicating that the impact of LCP_UW on the correlations between MODIS AOD retrievals and ground measurements is limited. In contrast, the MRE of the Terra and Aqua DT AOD retrievals presents statistically significant relationships (i.e., R 2 > 0.8, p < 0.05) with LCP_UW for both the 3 and 10 km resolution products. Pronounced correlations were also revealed between LCP_UW and PWEE, especially for products from MODIS Terra, where the R 2 values reached 0.85 and 0.95 for the 10 and 3 km resolution data, respectively. None of the accuracy measures of the MODIS DB10 products were correlated with LCP_UW. The same correlation analysis was also conducted for land cover proportion of water (LCP_W) and land cover proportion of urbanized features (LCP_U) (see Tables S1 and S2). Apparently, the accuracy levels of the MODIS AOD products were less correlated with LCP_W/LCP_U than they were with LCP_UW, further highlighting the impacts of both water and urban coverage on the accuracy of the AOD retrieval algorithms.

Correlations between Accuracy Levels and the Proportions of Water and Urban Areas
The accuracy measures (including R 2 , MRE, and percentage of data that fell within the EE (PWEE)) of the different MODIS AOD products are plotted against LCP_UW in Figure 7, and the corresponding correlations are listed in Table 2. Significant relationships between LCP_UW and R 2 were only found in all DT10 AOD retrievals, indicating that the impact of LCP_UW on the correlations between MODIS AOD retrievals and ground measurements is limited. In contrast, the MRE of the Terra and Aqua DT AOD retrievals presents statistically significant relationships (i.e., R 2 > 0.8, p < 0.05) with LCP_UW for both the 3 and 10 km resolution products. Pronounced correlations were also revealed between LCP_UW and PWEE, especially for products from MODIS Terra, where the R 2 values reached 0.85 and 0.95 for the 10 and 3 km resolution data, respectively. None of the accuracy measures of the MODIS DB10 products were correlated with LCP_UW. The same correlation analysis was also conducted for land cover proportion of water (LCP_W) and land cover proportion of urbanized features (LCP_U) (see Tables S1 and S2). Apparently, the accuracy levels of the MODIS AOD products were less correlated with LCP_W/LCP_U than they were with LCP_UW, further highlighting the impacts of both water and urban coverage on the accuracy of the AOD retrieval algorithms.  The strong correlations between LCP_UW and MRE found in the DT AOD retrievals also make it possible to predict the MRE and thus correct the potential uncertainties of DT AOD based on the underlying land cover information. The corrections can be conducted as: where and are the MODIS AOD and corrected MODIS AOD; is derived from the regression coefficients listed in the Table 3. The regression coefficients were calculated by linear regression between MRE and LCP_UW with a least square method.
A comparison of the bias of DT AOD retrievals before and after the correction over all the ground sites is shown in Figure 8. It is obvious that the bias of AOD products has  The strong correlations between LCP_UW and MRE found in the DT AOD retrievals also make it possible to predict the MRE and thus correct the potential uncertainties of DT AOD based on the underlying land cover information. The corrections can be conducted as: where AOD m and AOD corrected m are the MODIS AOD and corrected MODIS AOD; MRE predict is derived from the regression coefficients listed in the Table 3. The regression coefficients were calculated by linear regression between MRE and LCP_UW with a least square method. A comparison of the bias of DT AOD retrievals before and after the correction over all the ground sites is shown in Figure 8. It is obvious that the bias of AOD products has been substantially reduced through the correction with Equation (10). Especially for the lakeside and urbanized sites (TH and HZ), the overestimations have been decreased from the range of 0.1~0.3 to the range of −0.05 to 0.05 when the correction method was applied. The biases of DT AOD retrievals over vegetated sites (LA, TL, CA and JD) are not changed obviously before and after correction. This is principally because original DT AOD retrievals tend to have relatively small bias over vegetated areas and the performance of the correction scheme is limited.
The biases of DT AOD retrievals over vegetated sites (LA, TL, CA and JD) are not c obviously before and after correction. This is principally because original DT AOD als tend to have relatively small bias over vegetated areas and the performance correction scheme is limited.
To further demonstrate the availability of the proposed correction method, w pared the validation results of DT AOD retrievals at TH site before and after cor during year 2012 (shown in Figure 9). The performance of DT AOD retrievals was icantly improved when the correction method was applied. For instance, the MRE been decreased from the range of 27~39% to the range of −14 to −6%, the RMSEs fo AOD have been decreased from 0.19~0.20 to 0.13~0.14, and the PWEEs have b creased from 29~43% to 67~73%. The comparison results demonstrate the effective the LCP_UW-based MRE prediction model in improving the accuracy levels of D retrievals over high LCP_UW areas.   To further demonstrate the availability of the proposed correction method, we compared the validation results of DT AOD retrievals at TH site before and after correction during year 2012 (shown in Figure 9). The performance of DT AOD retrievals was significantly improved when the correction method was applied. For instance, the MREs have been decreased from the range of 27~39% to the range of −14 to −6%, the RMSEs for DT10 AOD have been decreased from 0.19~0.20 to 0.13~0.14, and the PWEEs have been increased from 29~43% to 67~73%. The comparison results demonstrate the effectiveness of the LCP_UW-based MRE prediction model in improving the accuracy levels of DT AOD retrievals over high LCP_UW areas.
Remote Sens. 2022, 14, x FOR PEER REVIEW 12 of Figure 9. Comparison of MODIS DT10 and DT3 with collocated ground measurements at TH s before correction and after correction during 2012. The abbreviations BC and AC in the figure re resent before correction and after correction, respectively.

Comparison of the Current Accuracy Levels with Previous Studies and Products
Considerable overestimations were found in the urbanized YRD site for the DT AO retrievals, which is similar to the results of previous validations [22,51,52]. A total of 44.9 of the DT10 AOD retrievals estimated using MODIS Terra and 42.5% estimated usin Aqua were above the EE, while these fractions were less than 20% for vegetated sites. Th overestimation was significant for the DT3 AOD product at the HZ site, with 75.7% (Terr and 75.8% (Aqua) of the AOD estimates falling above the EE. These results also agree well with the results of pioneer studies, which recognized urban cover as a key factor th impacts the assumed relationships between the spectral reflectance of visible and SWI bands, leading to reduced accuracies of MODIS DT AOD retrievals [38]. The same reaso could also explain the lower uncertainties of the 10 km DT data than those of the 3 k data, where the coverage of artificial surface within the 3 × 3 pixel region is lower in th former (41.9%) than that in the latter (62.2%).
Former studies demonstrated that the PWEE values (EE: ±(0.05 + 20%)) of MOD Collection 6 (C6) DT10 AOD products were 63% and 69% when validated using Asia and global AERONET datasets [7,23], respectively. As EE was generally determined u der the normal distribution assumption when the EE envelope contained about 68% (a proximately one standard deviation) of the match-ups [53], the C6 DT10 AOD's unce tainty could be indicated by the EE defined in the paper. The PWEE values of C6.1 DT1 AOD were generally greater than 70% in the vegetated sites of YRD and less than 52% the urbanized and lakeside sites of YRD, indicating that the uncertainties were overes mated in vegetated regions and underestimated in urbanized and lakeside regions und the current EE definition. In addition, the MODIS AOD C6.1 data used in the current stud

Comparison of the Current Accuracy Levels with Previous Studies and Products
Considerable overestimations were found in the urbanized YRD site for the DT AOD retrievals, which is similar to the results of previous validations [22,51,52]. A total of 44.9% of the DT10 AOD retrievals estimated using MODIS Terra and 42.5% estimated using Aqua were above the EE, while these fractions were less than 20% for vegetated sites. The overestimation was significant for the DT3 AOD product at the HZ site, with 75.7% (Terra) and 75.8% (Aqua) of the AOD estimates falling above the EE. These results also agreed well with the results of pioneer studies, which recognized urban cover as a key factor that impacts the assumed relationships between the spectral reflectance of visible and SWIR bands, leading to reduced accuracies of MODIS DT AOD retrievals [38]. The same reason could also explain the lower uncertainties of the 10 km DT data than those of the 3 km data, where the coverage of artificial surface within the 3 × 3 pixel region is lower in the former (41.9%) than that in the latter (62.2%).
Former studies demonstrated that the PWEE values (EE: ±(0.05 + 20%)) of MODIS Collection 6 (C6) DT10 AOD products were 63% and 69% when validated using Asian and global AERONET datasets [7,23], respectively. As EE was generally determined under the normal distribution assumption when the EE envelope contained about 68% (approximately one standard deviation) of the match-ups [53], the C6 DT10 AOD's uncertainty could be indicated by the EE defined in the paper. The PWEE values of C6.1 DT10 AOD were generally greater than 70% in the vegetated sites of YRD and less than 52% in the urbanized and lakeside sites of YRD, indicating that the uncertainties were overestimated in vegetated regions and underestimated in urbanized and lakeside regions under the current EE definition. In addition, the MODIS AOD C6.1 data used in the current study demonstrated noticeable improvements in the YRD when compared with the previous version of Collection 5 (C5) products. Specifically, validation by He et al. (2010) showed that the slopes of the linear fits between C5 DT AOD and concurrent AERONET measurements in this region ranged between 0.5 and 0.74, and C5 DT algorithm tended to overestimate the AOD in low aerosol loadings and underestimate the AOD in high aerosol loading conditions. Our results showed that the slopes of the C6.1 DT retrievals increased to~0.9 for the 10 km product and~1.0 for the 3 km AOD product, indicating that the systematic bias was effectively removed through the release of the new product.
As for the performance of the MODIS DB AOD products, when compared with large-scale validations, the PWEE values of the current regional study were generallỹ 56% for the observed sites, while global validations demonstrated higher PWEE values (EE: ±(0.05 + 20%)) of >72% [29]. This indicates that the uncertainties of DB AOD are overestimated globally and underestimated over the YRD region under the current EE definition. The PWEE of the DB products remained relatively stable given the rapid changes in land cover types among the six sites. Although the DB products showed slightly worse performance over vegetated sites than the DT products, the overestimation of the DB AOD product was significantly reduced over urbanized and lakeside sites, which is primarily due to the hybrid approach adopted in the DB algorithm [6,29]. The better retrievals by the DT algorithm over vegetated areas are likely linked with the use of more stable and effective spectral reflectance relationships between visible and SWIR bands than those utilized in the DB algorithm [6,7,30]. Similar to the validations with AERONET observations in China [22], India [44], and Africa [54], the DB algorithm tends to underestimate AOD over the YRD region. These problematic retrievals could either result from the overestimation of surface reflectance or the improper selection of aerosol models [40].
Another significant difference between the DT and DB AOD products is the number of match-ups. The total number of match-ups between the DT AOD retrievals and in situ measurements (DT10: n_Terra = 1896, n_Aqua = 1847; DT3: n_Terra = 1303, n_Aqua = 1165) was significantly lower than that for the DB retrievals (n_Terra = 2340, n_Aqua = 2324). Such remarkable improvements in data coverage tend to be associated with the updated C6.1 DB algorithm resulting in a rapid increase in the number of valid retrievals under certain interrupted conditions, such as minimal haze or thin clouds [6].

Impacts of Aerosol Types
To determine if the aerosol conditions were identical between the six ground sites, we compared the aerosol types over the ground sites for the DT3, DT10, and DB10 AOD products. Indeed, for any given time, the assumed aerosol types of the six ground sites were identical, indicating that the impact of aerosol type on the accuracy of MODIS AOD retrievals should be spatially uniform, minimizing the impacts of the changes in aerosol type on the relationships shown in Figure 7 and Table 2. Indeed, the aerosol type in the DB10 product remained stable (i.e., nonabsorbing urban/industrial model) for all ground sites throughout the study period. In contrast, seasonality of the aerosol type occurred for each site in the DT3 and DT10 products, with the moderate absorption fine (MAF) model dominating in the spring (March to May) and the weak absorption fine (WAF) model dominating in other seasons. To examine whether these seasonal changes in aerosol conditions could impact the relationships between the LCP_UW and the uncertainty levels of AOD products, we plotted the 10 km/3 km DT AOD against the sun photometer measurements under different aerosol types ( Figure 10). However, the points from the two aerosol models (i.e., WAF and MAF) exhibited certain degrees of differences in statistical accuracy measures (i.e., R 2 , RMSE, slope) and could not be distinguished from each other, regardless of the AOD algorithm and resolution, suggesting that the seasonal changes in the aerosol model have limited impacts on the correlations between LCP_UW and the MODIS DT AOD uncertainties. To further investigate the impact of aerosol type on the accuracies of the MODIS AOD retrievals, we compared the magnitudes of single-scattering albedo (SSA) between MODIS observations and AERONET measurements. SSA is a ratio between scattering and total extinction, ranging from 0 (purely absorbing aerosols) to 1 (purely scattering aerosols). Figure 11 shows the climatological seasonal SSA (675 nm) values derived from AER-ONET measurements at the TH site from 2005 to 2016. With mean SSA values of ~0.96 in the summer and ~0.95 in other seasons, the study region appeared to be dominated by nonabsorbing or weakly absorbing aerosols. Obviously, the SSA tends to be underestimated in the DT algorithm (especially in spring) and overestimated in the DB algorithm (see Figure 11), leading to the overestimation of DT AOD and underestimation of DB AOD [40]. Such comparison could also partially explain the systematic overestimation or underestimation of MODIS AOD retrievals, as shown in Figures 4-6. To further investigate the impact of aerosol type on the accuracies of the MODIS AOD retrievals, we compared the magnitudes of single-scattering albedo (SSA) between MODIS observations and AERONET measurements. SSA is a ratio between scattering and total extinction, ranging from 0 (purely absorbing aerosols) to 1 (purely scattering aerosols). Figure 11 shows the climatological seasonal SSA (675 nm) values derived from AERONET measurements at the TH site from 2005 to 2016. With mean SSA values of~0.96 in the summer and~0.95 in other seasons, the study region appeared to be dominated by nonabsorbing or weakly absorbing aerosols. Obviously, the SSA tends to be underestimated in the DT algorithm (especially in spring) and overestimated in the DB algorithm (see Figure 11), leading to the overestimation of DT AOD and underestimation of DB AOD [40]. Such comparison could also partially explain the systematic overestimation or underestimation of MODIS AOD retrievals, as shown in Figures 4-6.

Conclusions
In the present study, we provide a comprehensive comparison between the C6.1 MODIS AOD retrievals and ground-based measurements from six neighboring sites (HZ, LA, TL, CA, JD, and TH) in the YRD region from 2013 to 2016. Overall, the accuracies of the MODIS DT10, DT3, and DB10 products in the study region were comparable to those of the global validation results [29,30]. A significant difference was found in the performance of MODIS AOD retrievals over ground sites with different land cover types. While the AOD products derived from the DT algorithm (both 3 and 10 km resolutions) agreed well with the ground measurements over vegetated sites (LA, TL, CA, and JD), overestimations were found over urbanized and lakeside sites (HZ and TH). In contrast, underestimations occurred with the DB AOD algorithm for all ground sites, and this method demonstrated similar performances in both vegetated and urbanized sites. The overestimation and underestimation patterns of the DT and DB algorithms were due to the systematic errors of the SSA calculation during the AOD retrieval process.
The most striking finding of this study is that the LCP_UW within a region appeared to be correlated with the uncertainty measures of the MODIS AOD products. For example, the MRE of the Terra DT10 and DT3 products demonstrated statistically significant relationships (i.e., R 2 > 0.8, p < 0.05) with LCP_UW. These correlations would not have been found without two factors: (1) the relatively small distances between different ground sites, minimizing the impacts of changes in aerosol properties on the correlations, and (2) the significant differences in the land cover types surrounding the selected sites, with the LCP_UW ranging between 10 and 80%, making it possible to establish relationships between the accuracy measures.
This study demonstrated that the accuracies of satellite AOD estimates are not only associated with the coverage of urban areas but also modulated by the presence of water. Furthermore, we also showed that the performance of DT AOD retrievals over lakeside and urbanized areas could be significantly improved based on the relationship between MRE and LCP_UW. Future efforts are required to establish robust correlations between LCP_UW and uncertainty measures of AOD products at both regional and continental scales to correct the potential impacts from both urban and water areas.
Supplementary Materials: The following supporting information can be downloaded at: www.mdpi.com/xxx/s1, Figure S1: Comparison of Aqua MODIS DT 10 km AOD with collocated ground measurements over the six ground sites: (a) HZ, (b) LA, (c) TL, (d) CA, (e) JD, and (f) TH; Figure 11. Boxplot of seasonal single-scattering albedo (SSA) derived from AERONET measurements at the TH site. The SSA values estimated by the DT and DB algorithms are also marked in the graph. Spring (March, April, and May); summer (June, July, and August); autumn (September, October, and November); winter (December, January, and February).

Conclusions
In the present study, we provide a comprehensive comparison between the C6.1 MODIS AOD retrievals and ground-based measurements from six neighboring sites (HZ, LA, TL, CA, JD, and TH) in the YRD region from 2013 to 2016. Overall, the accuracies of the MODIS DT10, DT3, and DB10 products in the study region were comparable to those of the global validation results [29,30]. A significant difference was found in the performance of MODIS AOD retrievals over ground sites with different land cover types. While the AOD products derived from the DT algorithm (both 3 and 10 km resolutions) agreed well with the ground measurements over vegetated sites (LA, TL, CA, and JD), overestimations were found over urbanized and lakeside sites (HZ and TH). In contrast, underestimations occurred with the DB AOD algorithm for all ground sites, and this method demonstrated similar performances in both vegetated and urbanized sites. The overestimation and underestimation patterns of the DT and DB algorithms were due to the systematic errors of the SSA calculation during the AOD retrieval process.
The most striking finding of this study is that the LCP_UW within a region appeared to be correlated with the uncertainty measures of the MODIS AOD products. For example, the MRE of the Terra DT10 and DT3 products demonstrated statistically significant relationships (i.e., R 2 > 0.8, p < 0.05) with LCP_UW. These correlations would not have been found without two factors: (1) the relatively small distances between different ground sites, minimizing the impacts of changes in aerosol properties on the correlations, and (2) the significant differences in the land cover types surrounding the selected sites, with the LCP_UW ranging between 10 and 80%, making it possible to establish relationships between the accuracy measures.
This study demonstrated that the accuracies of satellite AOD estimates are not only associated with the coverage of urban areas but also modulated by the presence of water. Furthermore, we also showed that the performance of DT AOD retrievals over lakeside and urbanized areas could be significantly improved based on the relationship between MRE and LCP_UW. Future efforts are required to establish robust correlations between LCP_UW and uncertainty measures of AOD products at both regional and continental scales to correct the potential impacts from both urban and water areas.