Evaluation and comparison of MAIAC, DT and DB aerosol products over China

. A new Multiangle Implementation of the Atmospheric Correction (MAIAC) algorithm has been applied in the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor and has recently provided globally high spatial resolution 10 Aerosol Optical Depth (AOD) products at 1 km. Moreover, several improvements have been modified in the classical Dark Target (DT) and Deep Blue (DB) aerosol retrieval algorithms in MODIS collection 6.1 products. Thus, validation and comparison of the MAIAC, DT and DB algorithms are urgent in China. In this paper, we present a comprehensive assessment and comparison of AOD products at a 550 nm wavelength based on three aerosol retrieval algorithms in the MODIS sensor using ground-truth measurements from Aerosol Robotic Network (AERONET) sites over China from 2000 to 2017. In general, 15 MAIAC products achieved better accuracy than DT and DB products in the overall validation and accuracy improvement of DB products after the QA filter, demonstrating the highest values among the three products. In addition, the DT algorithms had higher aerosol retrievals in cropland, forest, and ocean land types than the other two products, and the MAIAC algorithms were more accurate in grassland, built-up, unoccupied land and mixed land types among the three products. In the geometry dependency analysis, the solar zenith angle, scattering angle and relative azimuth angle, excluding the view zenith angle, 20 significantly affected the performance of the three aerosol retrieval algorithms. The three products showed different accuracies with varying regions and seasons. Similar spatial patterns were found for the three products, but the MAIAC retrievals were smaller in the North China Plain and higher in Yunan Province compared with the DT and DB retrievals before the QA filter. After the


Introduction
Aerosols are a multi-compartment system consisting of suspended solid and liquid particles in the atmosphere, which play an important role in radiative forcing (Rajeev et al., 2001), regional climate (Qian et al., 1999;Feng et al., 2019a) and urban air pollution (Dominici et al., 2014).The aerosol optical depth (AOD) is the key aerosol optical parameter, defined as the vertical integration of the aerosol extinction coefficient from the ground to the top of the atmosphere (TOA).Ground measurements 2 Data Description Three aerosol products, e.g., MAIAC, DT and DB are stored in Hierarchical Data Format (*.hdf) files, and we obtain corresponding *.hdf files in the China region from 2000 to 2017 from the NASA Earthdata Search website (https://search.earthdata.nasa.gov/search/).Due to snow, cloud and desert ground surface types will increase the retrieval uncertainty, and three product provides a quality assurance (QA) flag to indicate the retrieval quality.Ground measurement aerosol data obtained from the AERONET website (https://aeronet.gsfc.nasa.gov/)were used to validate the accuracy of three satellite aerosol products.Additionally, land cover data from the Geographical Information Monitoring Cloud Platform (GIMCP, http://www.dsac.cn/)were utilized to analyse the land cover dependency for three satellite aerosol products.

DT products
The DT algorithm retrieves AOD parameters based on the assumption that the surface reflectance in two visible bands, e.g., 470 nm and 644 nm, presents a good linear relationship with the surface reflectance in the shortwave infrared (SWIR) band, e.g., 2119 nm, in dark, dense vegetated area, and the measurement in SWIR band is transparent with the aerosol particle (Kaufman et al., 1997, Levy et al., 2013).The surface and aerosol information can then be decoupled from the TOA spectral reflectance.Compared with the DT algorithm in collection 6.0, the DT algorithm in collection 6.1 mainly revises the surface characterization over the land surface when the urban percentage is larger than 20% (Gupta et al., 2016).
The DT algorithm produces two resolution aerosol products in collection 6.0 and 6.1, e.g., 3 km×3 km and 10 km×10 km.The two resolution products share the same retrieval protocol except the use of different retrieval boxes.For example, the 10 km product organizes 20×20 group pixels with the three aforementioned band measurements at 500 resolution into the retrieval box, whereas the 3 km product combines three band measurements in the 6×6 pixel group into a retrieval box (Remer et al., 2013).The comparison between the 10 km product and 3 km product from collection 6.0 on the global scale (Remer et al., 2013) and the China region (He et al., 2017) shows that the accuracy of the 10 km product is superior to one of the 3 km product, even though the 3 km product provides finer resolution aerosol retrievals.In this study, we consider the 10 km product of the newest collection 6.1 version from the Terra satellite.In DT products, QA=3 indicates high confidence data, and QA=1 indicates marginal confidence data (Levy, et al., 2013).In this paper, the scientific data set (SDS), named the "Image_Optical_Depth_Land_And_Ocean" QA filter and "Optical_Depth_Land_And_Ocean" with a QA filter (QA>1 for ocean and QA=3 for land) are extracted to compare the accuracy with and without the QA filter.

DB products
The DB algorithm retrieves the AOD parameter under the hypothesis that the surface reflectance in the deep blue band, e.g., 412 nm, is much smaller than in longer bands over bright surfaces, such as urban and desert regions (Hsu et al., 2004).First, the DB algorithm retrieves 1 km aerosol properties using the global surface reflectance database in visible bands, e.g., 412 nm, 470 nm and 650 nm, and then aggregates 1 km pixels into a 10 km scale.In collection 6.0, the surface reflectance database is improved using knowledge of the normalized difference vegetation index, scattering angle and season (Hsu et al., 2013).The ability to retrieve aerosol data over a bright surface for the DB algorithm greatly expands the coverage of aerosol retrieval.
The general principles for collection of the 6.1 DB products are still the same as those in the collection 6.0 version.The major improvements for collection 6.1 DB products are in the radiometric calibration, heavy smoke detection, artefact reduction over heterogeneous terrain, surface model in elevated terrain and regional/seasonal aerosol optical models (Hsu, 2017).
The same to the DT products, SDS named "Deep_Blue_Aerosol_Optical_Depth_550_Land" without the QA filter and "Deep_Blue_Aerosol_Optical_Depth_550_Land_Best_Estimate" with the QA filter (QA=2, 3 for land) in collection 6.1 from Terra satellite were selected for our study to validate the accuracy improvement by the QA filter.The solar zenith angle in "Solar_Zenith" SDS datasets, view zenith angle in "Sensor_Zenith" SDS datasets, solar azimuth angle in "Solar_Azimuth" SDS datasets, sensor azimuth angle in "Sensor_Azimuth" SDS datasets and scattering angle in "Scattering_Angle" SDS datasets were also assessed to determine the geometry dependence for DT and DB products.

MAIAC products
The MAIAC algorithm relies on the assumption that the surface reflectance changes slowly over time and shows high variability over space, whereas the aerosol loading changes very fast over time and varies only on a limited space scale.The main procedure of MAIAC is as follows: first, MAIAC resamples MODIS L1B measurements into a fixed 1 km grid, and then it adopts 4~16 day time series of resampled MODIS measurement to retrieve the surface Ross-Thick Li-Sparse (RTLS) bidirectional reflectance distribution function (Lucht et al., 2000) using the measurements in SWIR band.Subsequently, the linear spectral regression coefficient (SRC) between 470 nm and 2119 nm for each 1 km grid is retrieved instead of using the empirical regression coefficient in the DT algorithm.Finally, the AOD parameter at 470 nm can be retrieved by searching the minimum spectral residual between the theoretical TOA reflectance of the look-up table and the measurements in the red and SWIR bands.The AOD is originally retrieved at 470 nm, and the AOD parameter at 550 nm is computed using the AOD parameter at 470 nm based on spectral properties, expressed by the regional aerosol model from the MAIAC look-up table .The detailed MAIAC algorithm has been described by Lyapustin et al., 2011.Data used in this study were from the "Optical_Depth_055" and "AOD_QA" SDS datasets, and data were collected from the Terra satellite.The datatype of the "AOD_QA" SDS datasets is a 16-bit unsigned integer, and the best retrieved quality can be selected if 8~11 bytes of "AOD_QA" SDS dataset bits is "0000", which indicates the retrieval pixel and its adjacent pixel is clear (Lyapustin et al., 2018).The solar zenith angle in the "cosSZA" SDS datasets, view zenith angle in the "cosVZA" SDS datasets, relative azimuth angle in the "RelAZ" SDS datasets and scattering angle in the "Scattering_Angle" SDS datasets were also selected to analyse the view geometry dependence for MAIAC products.

AERONET data
AERONET is a global ground-based aerosol monitoring network that provides continuous optical and microphysical properties of aerosols at a 15-min sampling rate.The total uncertainty for the AERONET AOD parameter under cloud-free conditions is lower than ±0.01 for a wavelength longer than 440 nm and less than ±0.02 for shorter wavelengths (Holben et al., 1998).
Some studies were also conducted to examine the properties of these high-quality measurements in China (Liu et al., 2011;Xia et al., 2007;Li et al., 2007;Che et al., 20018, 2014;Bi et al., 2014).These high-accuracy datasets support various satellite AOD product evaluation research in the China region (Tao et al., 2015;Tian et al., 2018;He et al., 2017;Sogacheva et al., 2018).AERONET provides three quality levels of data, e.g., level 1.0, level 1.5 and level 2.0, in version 3. Here, we only 5 selected quality-assured level 2.0 data as ground-truth data to validate the satellite data.Figure 1 shows the locations of the selected 50 AERONET sites across China in this study.Table 1 reports the site name, longitude, and latitude of the selected sites.However, the AERONET site does not record aerosol measurements at 550 nm, and thus we interpolated the AOD parameter at 550 nm using the Ångström exponent in the two neighbouring bands at 500 nm and 675 nm (Ångström et al., 1929;Eck et al., 1999)  where  500 and  675 are the AOD parameter at 500 nm and 675 nm, respectively,  550 is the interpolated AOD parameter at 550 nm,  500−675 is the corresponding Ångström exponent and ln(*) is the logarithmic operator.

Land cover data
One key difficulty in the aerosol retrieval algorithm is to decouple surface and atmosphere information in the satellite apparent reflectance.Land cover information greatly affects atmosphere properties (Xu et al., 2018;Feng et al., 2019b).Understanding the uncertainties in a satellite aerosol retrieval algorithm for different land cover types is necessary.GIMCP land cover data with 30 m resolution in the years 2000, 2005, 2008, 2010, and 2013 were used in this study.The first level of GIMCP land cover data includes cropland, forest, grassland, water, built-up and unoccupied land.Among them, unoccupied land includes desert, gobi, saline-alkali soil, swampland, bare land and bare rock gravel, which mainly includes bright surfaces.The high spatial resolution and abundant land cover types support our studies.Figure 1 shows the first level land cover type across the China mainland in 2013.

The selected spatiotemporal window
There is only a small amount of matchup data between the satellite data and ground data when using the direct matching method, e.g., use only one pixel where the AERONET sites is located and ground measurement at the exact satellite overpass time, due to large amounts of missing data in AERONET or satellite data and the time delay between the satellite overpass time and AERONET sampling time.Therefore, under the assumption that aerosol information is homogeneous in a limited spatial and temporal lag, a suitable spatiotemporal window is often adopted to increase the matchup data number.Thus, satellite measurements in the spatial window around the AERONET sites are averaged, and ground measurements in the temporal window centred on the satellite overpass time are averaged.
For 10 km DT and DB products, the selected spatial window is often 50 km×50 km, and the temporal window is ±30 min (Ichoku et al., 2002, He et al., 2017, Tao et al., 2015).For MAIAC products, Matins et al. described five different spatial windows, e.g., 3 km, 15 km, 25 km, 75 km and 125 km, and four temporal windows, e.g., 30 min, 60 min, 90 min and 120 min, to validate the MAIAC product over South America (Matins et al., 2017).The results showed that 25 km×25 km and ±60 min are reasonable for the Terra satellite.For comparison with 10 km DT and DB products, we selected 30 km×30 km as the spatial window closest to the best spatial window for the MAIAC product and employed the best temporal window ±30 min of 10 km product because we also noticed that the validation accuracy is very close for ±30 min and ±60 min temporal window in the results of Matins et al., although the matchup data number of the ±60 min temporal window is more than one of the ±30 min temporal window (Matins et al., 2017).

Land cover types for the AERONET sites
The first level of GIMCP land cover data were used to label the AERONET site group.Due to the selected 30 km×30 km spatial window in section 3.1, we followed Matins et al.'s work and labelled the AERONET sites if the proportion of one land cover type exceeded 50% in the spatial window around the AERONET site.If there was no dominant land cover type, we defined the land cover type for this AERONET site as a mixed group (Matins et al., 2017).Except for the defined first level type in the GIMCP land cover data, we found some coastal AERONET sites in which the dominant region was ocean, so we defined the land cover type for these sites as the ocean group.
Table 2 shows the land cover types for each AERONET site in 2013.There were no land cover type changes for most sites except Hangzhou_City, Muztagh_Ata and NAM_CO site.For the Hangzhou_City site, the land cover type changed from cropland to mixed group from 2005 to 2008, potentially due to the process of urbanization.For the Muztagh_Ata site, the land cover type changed from unoccupied land to grass land from 2008 to 2010, and the land cover type for the NAM_CO site varied from grassland to the mixed group between 2008 and 2010.We labelled each matchup dataset for the three sites using the land cover type in the nearest year to the AERONET sampling time.

Statistical approach
The expected error (EE) envelope is often used to validate satellite retrieval uncertainties.More than 66% of retrievals falling within the expected error lines indicates good accuracy.For the DT algorithm, the EE envelope is generally defined as ±(0.05+0.15×AOD)over land, and over 66% of retrievals meet the defined expected error limits at the global scale (Levy et al., 2010(Levy et al., , 2013;;Remer et al., 2005).In the global scale validation for MAIAC product, over 66% of retrievals satisfy the ±(0.05+0.1×AOD)EE limits, demonstrating that the accuracy of MAIAC is relatively higher than the DT algorithm over land (Lyapustin et al., 2018).In the regional validation of South America for the MAIAC product, the EE envelope is defined as ±(0.05+0.05×AOD),and the fraction of retrievals within these EE limits are also over 66% (Matins et al., 2017).In our study, to compare DT and DB products, we adopted ±(0.05+0.15×AOD)as the EE envelope and calculated the proportion within the EE envelope (Within_EE) using equation ( 2): In addition to the EE envelope, we also adopted coefficient of determination (R 2 ) and Pearson correlation coefficient (R) to study the correlation between the satellite AOD and AERONET AOD.The root mean square error (RMSE) was also utilized to analyse the dispersion degree of accuracy of the satellite AOD.The mean bias (Bias) was used to describe the bias of the satellite AOD.These statistical indicators were calculated using equation ( 3)-( 6), respectively.
The AOD  and AOD  are the satellite AOD retrievals and AERONET data, respectively.The AOD ̅̅̅̅̅̅  and AOD ̅̅̅̅̅̅  are the corresponding mean values. is the matchup data number.
In order to compare the spatiotemporal completeness of three products, daily spatial completeness and the temporal completeness are defined by equation ( 7)-( 8).
All the statistical indicators are calculated for three products before and after QA filter to indicate the accuracy improvements and the reduction of spatiotemporal completeness by QA flag.

Overall accuracy comparison
Figure 2 shows the overall evaluation for MAIAC, DT and DB products before and after the QA filter.In total, MAIAC products have more matchup data than DT and DB products, which indicates the completeness of the MAIAC product may be higher than the DT and DB products.Before the QA filter, the statistic showed that 69.84% of retrievals fall within the EE envelope, indicating a good accuracy for MAIAC products in China.Compared with DT and DB products, only 53.64% and 55.66% of retrievals were determined for DT and DB products.Based on the R statistical result, the results for the three products were all greater than 0.9, indicating that the three products are all well correlated with the ground-truth AERONET data.By contrast, the  2 statistical result for MAIAC products, e.g., 0.847, was superior to those for the DT and DB products, e.g., less than 0.8.From the  statistical result, no significant bias was observed for the overall MAIAC product.However, according to the corresponding bias boxplot in different AOD bins, a slight overestimation was observed when the AODs were less than 0.5, and a slight underestimation when the AODs were between 0.5 and 1.The DT and DB products appeared to be less overestimated based on the  result.From the corresponding bias boxplot, the mean bias result for each different AOD bin was also almost greater than zero.After the QA filter, the correlation for the MAIAC product slightly improved, but the Within_EE result was slightly reduced, and the RMSE and  results increased.From the corresponding bias boxplot subfigure, the positive mean biases when the AODs were less than 0.2 increased compared with corresponding results before the QA filter, and the negative biases when the AODs were between 0.5 and 1 were reduced.These phenomena resulted in the reduced overall accuracy.The reason for the changes in these statistical indicators will be explained in section 4.2.For the DT and DB products, the overall accuracies were all improved after the QA filter.The improvement of the DB product was more obvious than one for the DT product.The Within_EE result was improved from 57.66% to 63.32%, and the mean biases in the bias boxplot showed no obvious overestimation trend after the QA filter.However, the DT product was still overestimated after the QA filter, and only a little improvement was achieved in the Within_EE result.3 shows the corresponding statistical results.At the low, moderate and moderate-high levels, all statistical indicators showed that the MAIAC product had better accuracy than the DT and DB products before the QA filter.At the high level, the DT product achieved the highest correlation with the ground-truth data and low RMSE results, 10 but the positive bias result for the DT product revealed that the overestimation phenomenon was more serious than for the other two products.After the QA filter, the accuracy of the DB product was higher compared with the other two products at the low level because the positive bias phenomenon became more severe for the MAIAC product at this level.At the moderate level, the MAIAC product demonstrated the best correlation and lowest RMSE results with a slightly higher positive bias than the DB product.At the moderate-high level, the MAIAC product remained the best quality product among the three products.
At the high level, the DT product achieved the best correlation and lowest RMSE with the highest positive bias.

Land cover type dependency analysis
Figure 3 shows a scatterplot figure of the MAIAC products in different land cover types before and after the QA filter.In total, MAIAC retrievals in cropland, built-up, grassland, ocean types were more accurate than forest, unoccupied land and mixed types according to the Within_EE results.After the QA filter, except for grassland, the accuracies all improved, and the improvement effect in ocean type was more obvious.
The high aerosol loading, e.g., AODs > 1, mostly emerged in cropland (Figure 3 a-i and a-ii) and built-up (Figure 3 d-i and dii) types due to biomass burning in the dry season and multiple human activities in the built-up area (Zhang et al., 2010;Wang et al., 2018).MAIAC retrieved AODs with a very high accuracy for the two land cover types.The R and R 2 results were over 0.93 and 0.84, respectively, and the Winthin_EE results showed that more than 74% of retrievals fell within the EE limits.In comparison, retrievals in cropland showed little bias, in contrast to a small positive bias in the built-up area, and RMSE results in the built-up area were smaller than those in the cropland area.This high retrieval accuracy in cropland and built-up regions can support relative studies on biomass burning and anthropogenic emissions.In evergreen forest areas (Figure 3 b-i and b-ii), the retrievals showed a good correlation with ground measurements, with 5  _ = 0.874,   = 0.904.However, the  2 results without and with the QA filter were all lower than 0.8, and only approximately 45% of retrievals fell within the EE envelope.The result is opposite to the conclusion that the MAIAC algorithm improves the dark target retrieval accuracy better than the DT algorithm in Lyapustin et al., 2011(Lyapustin et al., 2011).To eliminate the influence of retrieval accuracy in the specific site, Figure 4 shows a scatterplot figure of the forest AERONET site, ignoring the sites with matchup numbers less than 10.We can observe good performance in the Chiayi, Qiandaohu and 10

AOD
Xinglong sites, and the corresponding Within_EE results were all higher than 70%.The relatively inferior performance sites were Banqiao, Taipei_CWB.After the QA filter, the accuracies were improved to 76.19% and 61.79%, respectively.The site with the worst performance was only the Lulin site, where the MAIAC retrievals were systemically higher than the ground measurements, and less than 4% of retrievals fell within the EE limits.The percentage of forest type in the 30 km×30 km spatial window around the Lulin site always exceeded 80% in 2000, 2005, 2008, 2010 and 2013.This high proportion of forest type eliminates the influence of other mixed land cover type.The Lulin site is located in Taiwan peninsula, and thus the improper aerosol type in the MAIAC algorithm and cloud cover may explain the overestimation in the Lulin site (Lyapustin et al., 2018).For the grassland type (Figure 3 c-i and c-ii), over 83.68% of MAIAC retrievals fell into the EE lines before the QA filter, and the  2 =0.750, =0.875,RMSE=0.085, and Bias=-0.018results all showed good accuracy in the grassland type.However, after the QA filter, the accuracy was dramatically decreased with Within_EE=46.02%, 2 =0.687, =0.868,Bias = 0.051 and RMSE=0.114,representing the main reason for some of the decreased overall statistical results shown in Figure 2 for the MAIAC product after the QA filter.It is noteworthy that some values were underestimated when the AODs were less than 0.5, and these values were discarded after the QA filter.However, some overestimated values emerged when the AODs were very small.To identify the reason, we also performed a statistical validation for each grassland type site in Figure 5, excluding the site with a matchup number less than 10.Before the QA filter, the underestimated values were mainly in the NAM_CO and QOMS_CAS sites.These two sites are located in the Tibetan plateau.The MAIAC algorithm filled the AOD retrievals using climatology values, e.g., 0.014, in high altitude regions, e.g., over 4.2 km, and the QA for climatology values was 0111 (Lyapustin et al., 2018).After the QA filter, the climatology values were thrown away in the NAM_CO site.For the QOMS_CAS site, nearly 2.13% of pixels still had altitudes less than 4.2 km in the spatial window.MAIAC retrievals in these pixels were overestimated compared with the ground measurements.After the QA filter, the Within_EE results decreased from 92.26% to 38.53%.A severe underestimation phenomenon was found in Lanzhou_City site, in contrast to the positive bias in its closest SACOL site.The small matchup number for the Lanzhou_City site might be the reason for the underestimation phenomenon.A great improvement was found in the Muztagh_Ata site after the QA filter.MAIAC had good accuracy in unoccupied land cover type (Figure 3 e-i and e-ii), with Within_EE results of 67.44% and 71.43% before and after the QA filter, and R and R 2 results over 0.9 and 0.8, respectively.Figure 3 f-i and f-ii indicate that MAIAC also achieved better performance in the mixed land cover area, with Within_EE=66.80%and R=0.882.In the ocean area (Figure 3 f-i and f-ii), MAIAC algorithm retrievals seemed to be overestimated when AODs were small, and the R=0.796 result was a little worse than one of the other land types.After the QA filter, the overestimated values were discarded, and the accuracy was greatly improved from R=0.796, Within_EE=67.96% to R=0.921, and Within_EE=78.22%.
In comparison to DT and DB products, Table 4 shows the validation of the statistical results for the MAIAC, DT and DB products with different land type covers.In cropland area, the accuracy of the DT product was evidently better than that of the MAIAC and DB products according to the R 2 , R, and RMSE results.However, the values seemed to be overestimated compared with the MAIAC product, and the Within_EE result was a little smaller compared with the MAIAC product.In the forest area, the DT algorithm also achieved optimal accuracy compared with the MAIAC and DB products.However, only 5 56.23% of the retrievals met the EE limits, which was less than the DB product.In the grassland type region, the accuracies for the three products were all decreased after the QA filter, and we consider the validations of the three products to have all been influenced by the NAM_CO and QOMS_CAS sites.Compared with DT and DB products, the MAIAC product obtained the best retrieval accuracy.Owing to the overestimation phenomenon in the QOMS_CAS sites after the QA filter, the Within_EE result dramatically dropped from 83.68% to 46.02%.In Built-up, unoccupied land and mixed regions, the MAIAC 10 product performed better than the DB product, and the DB product was more accurate than the DT product.In the ocean region, the DT product was clearly more accurate than the DB and MAIAC products.-------" means no matchup pairs or that the matchup pairs number less than 10.

Cropland
Table 5 shows the validation accuracy for three products after the QA filter in four seasons.In cropland, the retrieval accuracies in autumn for the three products were better than in other seasons.For forest land types, three products showed a higher correlation in autumn than the other seasons, but the Within_EE values demonstrated the best results in winter, and the 5 corresponding results for DB products were clearly higher than for the other two products.In terms of grassland type, MAIAC and DB products were more accurate in summer and spring, respectively.In the built-up region, all products showed a high correlation in all seasons, but DT products were seriously overestimated.In unoccupied land, matchup pairs for MAIAC and DB products were more focused in spring, and MAIAC products performed better than DB products.A high correlation was also found for the three product in mixed and ocean regions in all seasons, but more MAIAC retrievals met the EE envelope 10 line.-------" means no matchup pairs or that the matchup pairs number less than 10.

Cropland
The Ångström exponent (AE) is a key parameter to describe aerosol particle size, and in general, local aerosol sources play a dominant role in aerosol regimes.We followed the work of Martins et al., 2017 to discover aerosol particle sizes in different land covers (Martins et al., 2017), Figure 6 shows a scatterplot of the AE (440 nm-675 nm) parameter versus AOD for different 5 land cover types.Our results were similar to those of Martins et al., 2017.The aerosol types in China are mainly fine-mode aerosol particles (AE>1).Some coarse-mode particles (AE<0.5)are mainly found in some regions with sparse vegetation, e.g.,

View of the geometry dependency analysis
To determine how the view geometry influence the accuracy for three retrieval algorithms, we analysed view geometry dependency using the following four angles: solar zenith angle (SZA), view zenith angle (VZA), scattering angle (SA) and relative azimuth angle (RAA) (Superczynski et al., 2017).We separated each kind angle into 10 bins and statistically analysed the AOD bias distribution in each bin.The results are displayed in Figure 8.
In terms of the solar zenith angle, the three retrieval algorithms all showed a strong dependency with different characteristics.
A slight downtrend along with SZA was found in the MAIAC algorithm, and the MAIAC retrievals seemed slightly overestimated when SZA was less than 40° and underestimated when it as larger.The mean biases only fluctuated between -0.05 and 0.05.For the DT algorithm, the mean bias first arose when the SZAs were small, and the mean bias reached the maximum at SZA≈25°.Then, the mean biases decreased as SZA increased.The mean biases were close to zero when SZA reached the maximum value.With regard to the DB algorithm, the mean bias first slowly decreased when the SZAs were less than 35° and then rapidly rose as SZA increased.After the QA filter, the whole mean bias line shifted down.
The MAIAC and DB algorithm showed no dependency on the view zenith angle, and the corresponding mean bias lines did not fluctuate much along with VZA.Compared with the results obtained before and after the QA filter, the mean bias line for the MAIAC algorithm slightly increased, and the mean bias line for the DB algorithm moves down to a relatively great degree.
VZA slightly affected the DT performance with a little downtrend.After the QA filter, the mean bias line slightly declined.
The scattering angle also greatly impacted the performance of the three retrieval algorithms.MAIAC retrievals seemed to be underestimated when the SAs were less than 100° and slightly overestimated when they were between 100° and 155°.When the SAs were larger than 155°, the retrievals tended to be underestimated.After the QA filter, the corresponding retrievals at large SAs tended to be overestimated.For the DT and DB retrievals, a significant uptrend was observed for the mean bias along with SAs.Small positive biases were found when the SAs were very small, and large positive biases occurred when the SAs were very large.After the QA filter, the significant uptrend was alleviated for DB retrievals, but a large negative bias was found when SA approached 180°.We consider the scarce matchup number of DB retrievals to be the main reason for the negative bias.
For the MAIAC algorithm, positive biases occurred as RAA approached the extremes of 0°, 180°, and negative bias emerged as RAA neared 90°, where the matchup numbers were very limited in the three angle intervals.In the rest angle intervals, MAIAC showed no dependency on RAA.After the QA filter, a downtrend of the mean bias was apparent along with RAA during backscattering (RAA<90°), and an uptrend of the mean bias was observed during forward-scattering (RAA>90°).For the DT algorithm, the positive mean bias decreased as RAA increased upon backscattering and first increased and then decreased upon forward-scattering.After the QA filter, the downward trend tended to be alleviated upon backscattering.For the DB algorithm, upon backscattering, the positive mean bias first decreased from very high to zero and then increased to become somewhat high.Upon forward-scattering, the positive mean biases were all larger than 0.05.After the QA filter and upon backscattering, no dependency on RAA was observed for the DB algorithm, but the highest mean bias was lower than zero.Upon forward-scattering, an obvious linear downtrend from positive to negative bias was observed as RAA increased.

Analysis of the spatiotemporal retrieval accuracy
To investigate retrieval accuracy of the three algorithms at different regions and different times, Figure 9 shows the R, RMSE, Bias and Within_EE results for each AERONET site, ignoring the sites with fewer than 10 matchup numbers, which might cause unreliable statistical results.Three products presented different retrieval accuracies in different regions.In the BTH region (marked by the black box in Figure 1), three products showed a good correlation with the ground measurements, e.g., R>0.9.There were, however, more retrievals for MAIAC and DT products falling within the EE limits than the DB product.Based on the Bias results, the DT and DB products seemed to be overestimated compared with the MAIAC product.The DT product was more positively biased compared with the MAIAC product.In the YRD region, the within_EE results showed that more MAIAC retrievals met the EE limits than DT and DB products.A good correlation of the three products was also found in this region.However, the DT product was overestimated, and DB was underestimated.In the PRD region, the MAIAC retrievals were obviously more accurate than the DT and DB retrievals.The Within_EE results for the MAIAC retrievals in this region were all greater than 70%.The Within_EE results for the DT retrievals were relatively low for some sites before the QA filter.After the QA filter, the Within_EE results were greatly promoted.DB retrievals in this region demonstrated the worst performance with low Within_EE results, a bad correlation and a negative bias.In addition, the MAIAC product was also the most accurate product in the NW area.The Within_EE and R results overall were higher than for the DT and DB products.Additionally, the RMSE results for the MAIAC product in this region were also relatively lower than those for the BTH and YRD region.The Within_EE results for the MAIAC product for most sites in the west of Taiwan were higher than 66% after the QA filter, demonstrating a high accuracy compared with DT and DB products.However, according to the east site, e.g., Lulin, the MAIAC retrievals seemed to be overestimated with low Within_EE and R results.Additionally, DB retrievals in the Lulin site seemed to be unbiased with high Within_EE (over 70%) and relatively high R (over 0.8) results.In the Tibet area, three algorithms all failed to retrieve AODs according to the statistical results due to the high latitude and snow cover.
Figure 10 presents the monthly validation results for the three products.We overlooked the specific QOMS_CAS site for this purpose due to its poor performance after the QA filter, which would affect the overall accuracy.Three products showed a good correlation with the ground measurements for all months with R>0.85.The AOD deviation for the DT and MAIAC products was higher in spring and summer than autumn and winter, consistent with the results of He et al., 2017(He et al., 2017).The RMSE results for the DB products were generally higher than the DT and MAIAC products before the QA filter.
After the QA filter, the RMSE results decreased with no obvious seasonal variability law.The DT product seemed to be systematically overestimated, and the positive biases were extremely high in spring and summer.The MAIAC product was positively biased from June to October with a Bias<0.1.The DB product was positively biased in all seasons before the QA filter, but the Bias results from June to October were significantly reduced after the QA filter.Before the QA filter, the Within_EE results for the MAIAC product were higher than the DT and DB products in all months.However, less than 60% of the MAIAC retrievals fell within the EE limits in summer.After the QA filter, the Within_EE results for the DB product from June to September were superior to those of the MAIAC and DT products.The R 2 results for the MAIAC products were stable for all months, and most R 2 results were over 0.8.The DB product had a lower R 2 in the cold season from November to February, and in April and May, the R 2 results for the DT product were generally lower than those in the other months.After the QA filter, the DB product achieved higher R 2 results from April to September.According to the matchup number results, the MAIAC product had more matchup numbers than the DT and DB products.However, all products had fewer matchup 5 numbers in summer due to the increased cloud cover in the rainy season.In summary, the MAIAC product was more accurate than the DT and DB products expect for in summer season.In contrast to the positive bias of MAIAC retrievals in summer, the DB product after the QA filter could achieve unbiased results with higher Within_EE and R 2 than the MAIAC product.To investigate the annual change in retrieval accuracy for three products to ascertain whether MODIS instrument maintain its performance due to it exceeding its designed lifetime.However, according to Table 1, the time durations of each AERONET site were significantly different.Thus, the matchup observation pair during each year were from different sites.This phenomenon may result from incomparable validation results for each year.However, if only considering the sites with the same monitoring time, most sites will be discarded, and fewer matchup numbers will cause unreliable corresponding statistical results.Thus, we still adopted all site measurements.We ignored the results for the years 2000, 2001, 2002 and 2003 due to fewer matchup numbers in these years.According to Figure 11, three products showed a high correlation with ground measurements according to the R and R 2 results, except in 2009.The reason for the sharp decline in the R and R 2 results in 2009 was mainly that some sites, e.g., Qiandaohu, SACOL, Kaiping, Shouxian, and Zhangye, did not have matchup pairs in this year, and matchup pairs containing bad retrieval satellite pixels around the Lanzhou_City and NAM_CO sites appeared in 2009.Thus, the high correlation revealed MODIS instrument results consistent the AOD retrievals from 2000 to 2017.The MAIAC AOD deviation was generally small, with most RMSE results being lower than 0.2 and larger than 0.15.The RMSE results for the DB product were generally larger than 0.2 before the QA filter.After the QA filter, the RMSE results varied in a large range from 0.15 to 0.25.Based on the Bias result, there was a significant uptrend for the three products over the year.
The MAIAC Bias results were generally smaller than the DT and DB products, and most Bias results for the MAIAC product were within ±0.05, with a negative bias before 2010 and a positive bias after 2010.To eliminate the influence of the contribution of some specific sites in the specific year, Figure 12 plots bias time series for five AERONET sites with a monitoring time covering most study years and ignoring data with matchup numbers less than 20.The bias uptrend seemed to appear in three products for all selected sites except the EPA-NCU site for DT products.Thus the significant uptrend of Bias results is not caused by the significantly different time durations of AERONET sites.For the Within_EE results, MAIAC also showed better accuracy than DT and DB products, and for Within_EE results, a slight declining trend was observed over the year.The matchup numbers for the three products revealed an increasing trend due to the establishment of greater numbers of AERONET sites in the China region over time.

Analysis on spatial pattern variation difference
To compare the difference in spatial variations for the three products, we upscaled the MAIAC product to match the grid of the DT and DB products; thus, 1 km pixels falling within the 10 km grid were averaged.Such a protocol can aid in investigating differences in different regions between the three products.
Figure 13 presents multiyear averaged and difference results between MAIAC, DT and DB products, with aerosol loading presenting a noteworthy assembly characteristic.Higher AOD values were concentrated in the North China Plain and Szechwan Basin where the land cover types were mainly cropland-oriented, as shown in Figure 1.Before the QA filter, compared with the DT and DB observations, the MAIAC AODs were smaller in the North China Plain and larger in Yunnan Province and East Taiwan.After the QA filter, the DB AODs became smaller in the North China Plain and southeast region.
Compared with the DB AODs, the MAIAC AODs became slightly higher in the North China Plain (difference over 0.1) and obviously higher in Southeast China (difference over 0.3).Recall the statistical result presented in Figure 9, in which the DT and DB products were overestimated in the BTH region, the DB product was underestimated in the YRD region, and the MAIAC product seemed to be overestimated in East Taiwan.These findings indicate that MAIAC retrievals are more accurate than DT and DB in the North China Plain and southeast region, and DB retrievals are more accurate than MAIAC in East Taiwan.However, due to the lack of the AERONET site in Yunnan Province, we could not evaluate the accuracy of the three products in Yunnan Province.The difference before and after the QA filter for the MAIAC product was very small, except for some individual pixels in the Tibet region.In addition, there was an obvious boundary in the 30° latitude according to the AOD spatial variation difference subfigure between MAIAC and the other two products.This boundary was caused by the different regional aerosol model used above and below 30° latitude (Lyapustin et al., 2018).Figure 14 show the seasonal comparison results among three products before and after the QA filter.The AOD spatial variation for the three products showed apparent seasonal characteristics.The AODs in the North China Plain in summer were higher 5 than in other seasons, and the AODs in the Tarim Basin in spring were higher than in other seasons.Based on the AOD spatial variation difference map, the difference between MAIAC and DT in the North China Plain evolved gradually from negative in spring to positive in winter.The negative difference between MAIAC and DB in the North China Plain was higher in summer and winter than in spring and autumn.The positive difference in Yunnan Province between MAIAC and DT was slightly lower than that between MAIAC and DB.After the QA filter, AODs in South China for the DB product were extremely 10 low compared with those for the MAIAC product.

Analysis of spatiotemporal completeness
Based the upscale MAIAC 10 km data in section 4.4, the spatial completeness in equation ( 7) and temporal completeness in equation ( 8) for three products are showed in Figure 15 and Figure 16.According to Figure 15, the spatial completeness of the MAIAC product was higher than the DT and DB products before and after the QA filter.The spatial completeness of the DT product was smallest due to its retrieval failure on a bright surface.The spatial completeness for all the products showed an obvious periodical trend change.Table 6 shows the statistics for the spatial completeness of the three products in different seasons.Before the QA filter, the averaged spatial completeness of MAIAC (46.87%) was higher than DT (16.66%) and DB (34.80%).After the QA filter, the reduced proportion of MAIAC ( 17    Plateau, the temporal completeness of the MAIAC product in this region was very high (over 80%).After the QA filter, the temporal completeness rapidly decreased in this region.In the other region, the declining proportions of temporal completeness 5 for MAIAC were mostly lower than 10%, except for Yunnan Province (nearly 15%), Hainan Province (nearly 20%) and East Taiwan (nearly 20%).Compared with the MAIAC and DB products, DT retrievals were very scarce in Tarim Basin due to failure on the bright desert surface.DT retrievals were more concentrated on the North China Plain and Yunan Province.After the QA filter, a dramatically reduced proportional area of temporal completeness (nearly 30%) for DT products was observed in the cropland region in Northeast China.The severely reduced proportional area (nearly 40%) for the DB product after the 10 QA filter was mainly focused on unoccupied land, e.g., gobi, saline-alkali soil, etc., at the top of the Tibet Plateau.Compared with the MAIAC product, before the QA filter, the DB product showed more retrievals in the Tarim Basin, North China Plain and Southeast China, and fewer retrievals in Yunnan Province and Northeast China.After the QA filter, the temporal completeness of the MAIAC product was better than the DB product in all regions.

Conclusion
In this study, we present the first comprehensive validation and comparison of three MODIS aerosol retrieval algorithms (i.e., 5 MAIAC, DT and DB) across China in terms of overall accuracy, land cover dependency, viewed geometry dependency, spatiotemporal retrieval accuracy, spatial distribution difference and spatiotemporal completeness.These validation results may guide users to utilize the three products appropriately.The main results and conclusion are presented below.

−
In terms of overall accuracy, the MAIAC product is more accurate than the DT and DB products.The DT and DB products are positively biased before the QA filter, and the positive bias for the DB product is alleviated by the QA filter.10 − DT retrievals in cropland, forest and ocean seem to be more accurate but with a positive bias than retrievals by the MAIAC and DB algorithms.The MAIAC algorithm performs better in grassland, built-up and mixed areas than the DT and DB algorithms.

Figure 1 :
Figure 1: Locations of the selected AERONET sites around China displayed on the land cover map from 2013.BTH: Beijing-Tianjin-Heibei; YRD: Yangtze River Delta; PRD: Pearl River Delta; NW: northwestern China.

Figure 2 :
Figure 2: Overall accuracy evaluation of MAIAC, DT and DB AOD versus AERONET AOD at 550 nm before and after the QA filter.The black line, red line and dashed line in the scatterplot are the 1:1 reference line, regression line and expected error (EE=±(0.05+0.15×AOD))line, respectively.The AOD bias boxplot uses the 25% and 75% percentiles with 50 bins.The red points in the boxplot are the mean bias for 50 bins.5

Figure 3 .
Figure 3. Evaluation of the MAIAC accuracy for different land cover types before and after the QA filter.The black line, red line and dashed line in the scatterplot are the 1:1 reference line, regression line and expected error (EE=±(0.05+0.15×AOD))line, respectively.

Figure 4 .
Figure 4. Evaluation of the MAIAC accuracy in the forest area for each AERONET site before and after the QA filter.The black line, red line and dashed line in the scatterplot are the 1:1 reference line, regression line and expected error (EE=±(0.05+0.15×AOD))line, respectively.

Figure 5 .
Figure 5. Evaluation of the MAIAC accuracy in the grassland area for each AERONET site before and after the QA filter.The black line, red line and dashed line in the scatterplot are the 1:1 reference line, regression line and expected error (EE=±(0.05+0.15×AOD))line, respectively.

Figure 7 .
Figure 7. Scatterplot of the AOD bias from matchup data versus the AERONET Ångström exponent (440 nm-675 nm) before and after the QA filter.

Figure 8 .
Figure 8. Dependency of the AOD bias on the solar zenith angle, view zenith angle, scattering angle and relative azimuth angle for the (a) MAIAC product, (b) DT product and (c) DB product before and after the QA filter.The dark blue bar is the histogram bin, red points in the shadowed area are the mean bias of the corresponding bin, and the top and bottom blue line are the 75% and 25% percentiles of the AOD biases in the corresponding bin, respectively.

Figure 9 .
Figure 9. Evaluation results for MAIAC, DT, and DB after and before the QA filter in each AERONET site.The subscript QA denotes the corresponding results after the QA filter.

Figure 10 .
Figure 10.Validation of MAIAC, DT and DB in different months before and after the QA filter.

Figure 11 .Figure 12 .
Figure 11.Validation of MAIAC, DT and DB in different years before and after the QA filter from 2004 to 2017.5

Figure 13 .
Figure 13.Averaged AOD distributions throughout the year for MAIAC, DT, and DB before the QA filter and their differences after the QA filter from 2000 to 2017.The subscript QA denotes the corresponding results after the QA filter.

Figure 14 .
Figure 14.Seasonal averaged AOD distributions for MAIAC, DT, and DB and their differences before and after the QA filter from 2000 to 2017.The subscript QA denotes the corresponding results after the QA filter.
.18%) exceeded DB (15.30%) and DT (8.66%) because many climatology values in the Tibet Plateau were discarded.Comparison of the spatial completeness in four seasons revealed a higher spatial completeness for the three products in autumn than the other three seasons due to the reduced cloudiness in the dry autumn season.The spatial completeness in winter was smallest due to the influence of the surface snow cover and large deciduous trees.Compared with MAIAC and DB products, the spatial completeness of the DT product in winter was minimal due to the bright surface in winter.

Figure 15 .
Figure 15.Daily spatial completeness for MAIAC, DT and DB from 2000 to 2017 before and after the QA filter.

Figure 16
Figure 16 presents the temporal completeness in China for the three products.Due to the climatology values in the Tibet

Figure 16 .
Figure 16.Spatial distributions of temporal completeness for MAIAC, DT, and DB before and after the QA filter and their differences from 2000 to 2017.The subscript QA denotes the corresponding results after the QA filter.