Evaluating the assumptions of surface reflectance and aerosol type selection within the MODIS aerosol retrieval over land: The problem of dust type selection

Abstract. Aerosol Optical Depth (AOD) and Angstrom exponent (AE) values derived with the MODIS retrieval algorithm over land (Collection 5) are compared with ground based sun photometer measurements at eleven sites spanning the globe. Although, in general, total AOD compares well at these sites (R2 values generally over 0.8), there are cases (from 2 to 67% of the measurements depending on the site) where MODIS clearly retrieves the wrong spectral dependence, and hence, an unrealistic AE value. Some of these poor AE retrievals are due to the aerosol signal being too small (total AOD


202
T. Mielonen et al.: MODIS aerosol type selection measurements can provide global aerosol information. However, because of the huge variety of aerosol and surface conditions, passive satellite algorithms are underdetermined in their information content. Therefore, to retrieve realistic aerosol properties, the algorithms require reasonable assumptions and constraints on aerosol properties and the contributions of the underlying surface to the observed radiance field. The quality of the assumptions is critical, especially over land where the surface reflectance can dominate the reflectance measured at the top of the atmosphere. Due to the stronger signal from the surface, the error in the derived aerosol optical depth (AOD) is typically 10 times larger than the error in the retrieved surface reflectance (Kaufman et al., 1997). Consequently, developing satellite retrieval algorithms is hard.
One of the most widely used passive satellite instruments in aerosol remote sensing is the Moderate Resolution Imaging Spectroradiometer (MODIS, Salomonson et al., 1989). MODIS instruments are aboard the Terra and Aqua satellites. Having a swath width of 2300 km, Terra MODIS and Aqua MODIS cover the Earth's surface every 1 to 2 days. Beginning in early 2007, a consistent algorithm (Levy et al., 2007a) was used to process the entire MODIS data record, providing the so-called Collection 5 (hereafter C5) data set. Levy et al. (2007a, and the on-line ATBD) have provided a detailed documentation of the algorithm and the changes from previous versions.
The C5 over land algorithm differs from prior versions in its retrieval process logic, as well as in its assumptions about aerosol model optical properties and surface reflectance. The algorithm makes a simultaneous inversion of MODIS-observed reflectance in three channels (0.47, 0.66 and 2.13 µm), to retrieve total aerosol loading, proportion of fine/coarse aerosol particles, and surface reflectance. The products are three semi-independent parameters, the total AOD (or τ ) at 0.55 µm, the nondust or fine model weighting (FMW), and the surface reflectance at 2.13 µm. FMW is the fractional contribution of fine mode aerosol to the total AOD at 0.55 µm.
C5 uses five aerosol models (Levy et al., 2007b), four of which are used in the standard algorithm: dust, weakly absorbing, moderately absorbing and strongly absorbing. Each of these four models has been determined by analysis of AERONET almucantar inversions (e.g. Dubovik et al., 2002), and are comprised of two lognormal modes in their size distribution. The dust model is non-spherical and dominated by coarse mode, whereas the other three are spherical and dominated by the fine mode. The three fine models are different mainly in their single scattering albedo (SSA), such that the weakly, moderately and strongly absorbing models have approximate SSA at 550 nm of 0.95, 0.91 and 0.86, respectively. The dust model has a SSA of 0.95 . The value of the FMW represents the "mixture" of one of the fine-dominated models (hereafter fine model) and the dust model, where the fine model is prescribed by geog-raphy and season. Thus, the FMW is not a weight between fine and coarse modes, but represents the proportion of nondust aerosol model. A FMW value of zero represents the case where the dust model provides the best solution. The retrieved spectral AOD represents the optical contributions of the dust/non-dust mixture, from which theÅngström exponent (AE) can be derived.
Validation with data from the Aerosol Robotic Network (AERONET, Holben et al., 1998) shows that retrieval of AOD has improved compared to C4 (Levy et al., 2007a;Mi et al., 2007;Papadimas et al., 2008;Jethva et al., 2007). Moreover, MODIS C5 measurements have been extensively compared with AERONET measurements globally Levy et al., 2007a;Jethva et al., 2010;Jethva et al., 2007;Mi et al., 2007;Kaskaoutis et al., 2007;Papadimas et al., 2008;Misra et al, 2008;Drury et al., 2008;Remer et al., 2008, Mishchenko et al., 2010. Several studies, for example Jethva et al. (2010) and Oo et al. (2008) have concluded that the surface assumptions and aerosol models in the MODIS retrieval are not representative for all locations. As for AE, there have been no comprehensive studies comparing MODIS and AERONET values of AE. Levy et al. (2007a) did not find high levels of correlation from their small sample of comparisons. Mishchenko et al. (2007), Liu and Mishchenko (2008) and Mishchenko et al. (2009) compared MODIS to retrievals from the Multi-angle Imaging SpectroRadiometer (MISR, Kahn et al., 2001) and found no correlation between the two instruments, further indicating that the information observed by passive satellite techniques over land is weak. Although Kahn et al. (2009) showed that some of the poor comparisons were resulting from data misinterpretation, retrieval of global AE is still considered a difficult problem.
Why are the AE discrepancies so large, even when total AOD is comparable? Why does MODIS retrieval report 100% dust, even in regions, where no significant dust loading should be expected? This study addresses the question why the MODIS algorithm wrongly selects the dust model in such cases. Previous studies have suggested, and  have stated, that the MODIS observations do not contain enough information to robustly separate dust from non-dust. Uncertainties in the surface reflectance, as well as in aerosol absorption are too large. However, from the collocations with sun photometer data, we can begin to quantitatively evaluate these systematic errors. In some sites the misidentification occurs because the surface reflectance assumptions are consistently wrong, whereas at other sites, the error is dominated by the initial assumption of aerosol model type.
In this paper, we study some of the sites where MODIS is consistently picking the wrong aerosol model, which in turn leads to errors in total AOD and especially in AE. We determine the dominant cause of aerosol type mis-identification, and determine whether the retrieval can be tuned in a physically meaningful way at these sites.

MODIS data
MODIS C5 data from both satellites are used in the comparisons with the ground-based measurements. The expected error over land in the MODIS AOD (or τ ) is ±0.05 ± 0.15τ . We used level 2 AOD data which has a spatial resolution of 10 × 10 km 2 at nadir. Only AOD data with best quality (Quality flag = 3) were considered in the analysis . For the comparisons with sun photometer measurements we used Optical Depth Land And Ocean product. In the more detailed analysis of the MODIS retrieval algorithm we used Corrected Optical Depth Land product. These data were downloaded from LAADS Web (http://ladsweb.nascom.nasa.gov/ data/search.html).

Sun photometer data
Two different types of sun photometer instruments were used in this study. The Precision Filter Radiometer (PFR; Wehrli, 2000;Kim et al., 2008) is a passive instrument which measures direct solar irradiance in four narrow spectral bands (368, 412, 500, and 862 nm). The bandwidth of the instrument is 5 nm and the full field of view angle is 2.5 degrees. Derived products from the measurements are AOD for the four wavelengths and AE. The time resolution of the products is 1 min. The absolute uncertainty in AOD of the PFR instrument is between 0.01 and 0.02 (Carlund et al., 2003). For the comparison with MODIS, the AOD measured by the PFR at 500 nm was interpolated to 550 nm using theÅngström power law, τ 550 = τ 500 (550/500) −α , where τ 550 is the AOD at 550 nm, τ 500 is the AOD at 500 nm and α is the AE based on the linear fit taken over all four PFR wavelengths. Cloud screening of the PFR measurements is based on the method presented by Smirnov et al. (2000). The Aerosol Robotic Network (AERONET; Holben et al., 1998) uses Cimel sun photometers which measure AOD at 340, 380, 440, 500, 675, 870 and 1020 nm. Measurements are provided every 15 min during daytime. AERONET also provides the angular distribution of sky radiances at four wavelengths (440, 670, 870 and 1020 nm), and aerosol properties such as aerosol size distribution, complex refractive index, and single scattering albedo (SSA), once every hour in clear sky conditions. The direct-sun AOD from AERONET are accurate to within ±0.01 for wavelengths larger than 400 nm and ±0.02 for shorter wavelengths (Holben et al., 1998;Eck et al., 1999). Holben et al. (1998) have described AERONET measurements in more detail. For the comparison with MODIS, the AOD measured by AERONET at 500 nm was interpolated to 550 nm by using theÅngström power law. We used cloud screened and quality assured level 2.0 data and level 2.0 inversion products in the analysis.
We compared the AERONET (Cimel) and PFR AOD measurements carried out at Sodankylä between March and September, 2007. For 1790 coexistent measurements with less than 3 minute time difference, the comparison showed that the ground based instruments agree extremely well at all comparable wavelengths. At visual and UV wavelengths Table 2. Statistics of the concurrent AERONET/PFR and MODIS AOD measurements. Subscript "all" refers to all coexistent measurements, "over" refers to measurements with MODIS AE over 1, and "under" refers to measurements with MODIS AE under 1. n is the number of measurements, R 2 is the correlation coefficient, mean(gb-M) is the average difference between the ground-based and MODIS measurements, RMS diff is the Root-Mean-Square Difference, and Pixels is the percentage of pixels within the expected accuracy of MODIS. Coarse AERONET and Coarse MODIS are the percentages of coarse particle measurements at each site based on AERONET and MODIS data, respectively. (368, 412 and 500 nm) the correlation coefficients (R 2 ) values were 0.99. For the longer wavelength band (862 nm) the correlation is almost as good, with R 2 = 0.97. For the all the wavelengths the bias between the instruments was less than 0.01. This and previous analyses (e.g. Kim et al., 2008) show that PFR and AERONET measured AOD are of similar quality.

Collocation of MODIS and sun photometer data
The ground-based measurements were used in the analysis only if the measurements were made within an one hour window centered at the MODIS overpass. The averaging of the satellite and surface measurements was based on the method described by Ichoku et al. (2002). If MODIS had more than 5 aerosol retrievals inside a grid of 5 × 5 pixels (50 × 50 km 2 ), centered at an AERONET/PFR station, these retrievals were averaged and compared to ground-based measurements. In addition, ground-based measurements were used in the analysis only if there were more than 6 PFR or 3 AERONET measurements within the one hour window. This was done to ensure that one instrument did not measure in the beginning of the time window and the other one in the end. Comprehensive global comparisons between AERONET and MODIS are presented by Levy et al. (2010), here we compare MODIS and sun photometer data for the specific sites, where there are many cases where MODIS is not identifying dust properly. Table 2 summarizes the differences be-tween MODIS and AERONET/PFR AOD measurements at the sites studied. All parameters are presented for all measurements, for the measurements with MODIS AE over 1, and for the measurements with MODIS AE under 1. The measurements were stratified by the MODIS AE in order to study how errors in AE are tied to errors in AOD. In addition, Table 2 includes the percentage of coarse particle measurements at each site based on AERONET and MODIS data (amount of measurements with AE under 1). The AERONET values show that in Kanpur 49% of the coexistent measurements has coarse aerosols while in Xianghe the fraction is 30%. The other sites do not have significant amounts of coarse aerosols. However, AE data from MODIS suggests that the coarse fraction is larger than 15% at all the sites, except in Ispra where the fraction is only 2%. In Kanpur the MODIS coarse fraction is as high as 92% while in Xianghe it is 84%.
At Alta Floresta, GSFC, Ispra, Mexico City, Mongu, and Tõravere the R 2 values are the lowest for the cases with AE under 1. At Jokioinen, Rome, Sodankylä, and Xianghe the values of R 2 for the measurements with AE under 1 are slightly lower or in the same range as in the other classes. At Kanpur the R 2 values are equal for the AE > 1 and AE < 1 classes. The averaged difference between the instruments (ground based -MODIS, henceforth mean(gb-M)) and the Root-Mean-Square of the difference (RMS diff) are very similar for all categories. In addition, Table 2   accuracy. At all sites the percentages are larger or in the same range as the results presented by Remer et al. (2008) and Jethva et al. (2007). Both reported that more than 70% of the measurements at 550 nm over land were within expected accuracy. Figure 1 shows examples of comparisons of the AODs from MODIS and ground-based instruments at 550 nm wavelength stratified by MODIS AE. For some sites (e.g. Jokioinen, Mongu), dust is known to be rare (see Table 2). However, it is common at those sites that MODIS retrieves AE < 1. The question is why? We found that: (1) Retrievals with AE < 1 are connected with AOD < 0.5. (2) MODIS retrieves AE < 1 at all studied sites, thus the flawed dust model selection is not confined to specific aerosol types or surface types. For the sites where also AERONET indicated the presence of coarse aerosol particles (e.g. Kanpur and Xianghe, Table 2) the MODIS AE is under 1 also for larger AOD values and the R 2 values for the two size classes (under and over 1) are equal or in the same range. Figure 2 shows the same data as Figure 1 but stratified by AERONET AE. Evidently, the number of measurements with AE < 1 is significantly lower than with the MODIS data at all sites. In addition, most of the measurements with AERONET AE < 1 are associated with very low AOD values.

SSA studies
We used single scattering albedo (SSA) data from AERONET level 2.0 inversion products to study whether absorption could explain the differences between MODIS and AERONET AOD values. To avoid measurement uncertainties associated with low AOD values, only measurements with AOD at 550 nm over 0.3 (both MODIS and AERONET) were considered. The AERONET SSA value was thought to be representative for the MODIS retrieval if the time difference between the measurements was less than 10 h. Figure 3 shows two examples of comparisons of SSA spectra for MODIS a-priori aerosol models with SSA values from the AERONET inversion categorized by MODIS AE. At Tõravere, the a-priori SSA spectra in the MODIS algorithm are in agreement with the AERONET data for both size classes. The classification of the coarse aerosol particles coincides with a strong dust plume advected from the Sahara on the 22 March 2007 (based on HYSPLIT backtrajectories (Draxler and Rolph, 2010;Rolph, 2010), thus for these cases MODIS selects the aerosol models correctly. On the other hand, the SSA spectra measured by AERONET in Mexico City are lower than those for the a-priori model, at each wavelength. Similar behavior was also detected at Ispra and Mongu. On the basis of AERONET AE, none of these sites were influenced by coarse particles during these measurements. This can also be seen from the spectral behavior of the AERONET SSA. At Tõravere, the averaged SSA of the cases with MODIS AE < 1 looks like the dust model, and the average of the cases with MODIS AE > 1 looks like the fine model. At Mexico City, on the other hand, the spectral dependence of the AERONET SSA does not look like either of the MODIS assumed models. However, both cases have the same wavelength dependence, indicating more absorbing fine mode aerosols than assumed in the model. The only difference between the types is that on average the measurements where MODIS detects AE < 1 are associated with smaller AERONET SSA values. This raises the question whether the absorption characteristics of the fine mode aerosols could affect the selection of the aerosol model combination in the MODIS aerosol retrieval. However, before we  Fig. 3. MODIS aerosol model SSAs compared with AERONET SSA products at two sites: Tõravere and Mexico City. SSA spectra for a-priori dust and fine models (moderately absorbing) used in MODIS retrieval are shown in red and blue, respectively. The magenta line represents the averaged SSA spectrum from AERONET measurements for observations whose MODIS AE is over 1 (fine model). The green color is an average for observations whose MODIS AE is smaller than 1 (dust model). At Tõravere the dust measurements are from a single dust episode (22nd March 2007). Only measurements with AOD(550 nm) values larger than 0.3 (MODIS and AERONET) are considered. 28 Fig. 3. MODIS aerosol model SSAs compared with AERONET SSA products at two sites: Tõravere and Mexico City. SSA spectra for a-priori dust and fine models (moderately absorbing) used in MODIS retrieval are shown in red and blue, respectively. The magenta line represents the averaged SSA spectrum from AERONET measurements for observations whose MODIS AE is over 1 (fine model). The green color is an average for observations whose MODIS AE is smaller than 1 (dust model). At Tõravere the dust measurements are from a single dust episode (22 March 2007 can answer this question we have to study how the aerosol model combination is affected by surface parameterizations, see below.

Surface reflectance studies
In the next step, we investigated how the C5 surface reflectance estimates influence the retrieval of AE. We used individual pixels and divided the data into two classes based on MODIS AE as before. Then we calculated the monthly mean surface reflectances at 660 nm from all measurements and normalized the surface reflectance values with these mean values to remove seasonal variations. This was done to remove the effects of natural variation in surface reflectance properties due to seasonal cycle of vegetation. Figure 4 presents the normalized surface reflectances for measurements with MODIS AE over 1, and under 1, as a function of AERONET AE. The surface reflectance data was binned into 0.1 wide AE classes. This example from Mongu is also representative for other sites which had a sufficient number of measurements from both classes (Alta Floresta, Kanpur, Xianghe). At these sites, measurements with MODIS AE over 1 had larger normalized surface reflectance values than the measurements with MODIS AE under 1. This was also seen in the 2.13 µm surface reflectance data. Surface reflectance should not depend on the AE, as the AERONET AE data shows in Fig. 4 where the normalized surface reflectance for the MODIS measurements with AE > 1 is almost constant as a function of AERONET AE. However, in the MODIS retrieval, larger particles (dust model) are associated with lower surface reflectance values than small particles (fine model). Based on these results, it seems that both aerosol absorption and surface parameterizations could have effects on the aerosol model combination in the MODIS retrieval. In order to separate their effects, we have to study the retrieval algorithm in more detail.

Analysis of the MODIS retrieval algorithm
To determine why MODIS tends to wrongly identify dust in certain scenes, we need to focus on the details of the MODIS algorithm. Specifically, we try to determine which algorithm assumptions are violated, and why the retrieval behaves as it does.
The MODIS aerosol retrieval assumes that the measured reflectance (ρ * ) at the top of the atmosphere (TOA) is composed of atmospheric "path radiance" (ρ a (θ,θ 0 ,φ)) and transmission of surface reflectance through the atmosphere:  Here θ is the view zenith angle, θ 0 is the solar zenith angle and φ is the azimuth angle of the scattered radiation from the solar beam. F d (θ 0 ) is the normalized downward flux for zero surface reflectance, equivalent to the total downward transmission, T (θ ) is the upward total transmission into the satellite field of view, and s is the atmospheric backscattering ratio (Kaufman et al., 1997).
Since values of ρ a (θ,θ 0 ,φ), F d (θ 0 ), T (θ), and s are unknown properties of the atmosphere, and the surface reflectance is also unknown, Eq. (1) cannot be solved directly. Instead, the MODIS retrieval algorithm is designed to solve Eq. (1), by referring to the lookup tables (LUT) containing simulated values of ρ a (θ,θ 0 ,φ), F d (θ 0 ), T (θ), and s. Therefore, Eq. (1) can be inverted, solving for the surface reflectance, e.g., Starting with the 2.13 µm wavelength, the MODIS algorithm determines the surface reflectance that solves the equation, repeating for each of seven aerosol loadings (corresponding to AOD at 0.55 µm of 0.0, 0.25, 0.5, 1.0, 2.0, 3.0, and 5.0), and for the fine and dust models, separately. .
By these equations, the surface reflectance in the two visible channels (470 and 660 nm) can be estimated from the reflectance at 2.13 µm. Equation 1 is then applied to estimate TOA reflectance (ρ * ) in the two visible wavelengths. Note that the result is the estimated TOA reflectance in three channels (470 nm, 660 nm and 2.1 µm), corresponding to seven aerosol loadings, for the fine and coarse models, separately. These estimated TOA values, can be compared with the observed TOA reflectance. By interpolating between entries of aerosol loading, the algorithm finds the exact loading which satisfies equality between the estimated and observed TOA reflectance at 470 nm. Now, the algorithm can mix the fine and coarse models. By iterating the model mixing (FMW) through steps of 0.1, while simultaneously recalculating surface reflectance that satisfies the surface reflectance parameterization, the algorithm arrives at different TOA estimates for different mixtures. Equality is required at 470 nm and 2.1 µm, however there will be differences between estimated and observed values at 660 nm. The particular loading/mixture that minimizes the 660 nm difference is considered as the best solution, thus retrieving both the AOD and FMW. Since spectral AOD is consistently represented by the LUT, it is easy to report AOD at 470 and 660 nm, and derive AE from these values with the Angström power law.
Having AE of 0.6, the dust model has a flatter spectral dependence than any of the fine models (AE = 1.6). Thus, for a given loading index (AOD at 550 nm), the coarse (dust) model contributes more to the TOA reflectance at 660 nm than does the fine model. Since the algorithm has to find the best match at 660 nm, this means that for a given value of surface reflectance at 660 nm, the algorithm must either adjust the total aerosol loading or the relative contribution of the coarse model. If the surface target is "brighter" than expected via the surface parameterization, MODIS may observe larger TOA reflectance than modeled by Eq. ( 1). In this case, MODIS will compensate by picking a dust only aerosol mixture. As a result, MODIS will indicate dust aerosol even when other observations do not indicate dust in the column.

MODIS standalone algorithm
In order to trace the performance of the MODIS retrieval algorithm, we obtained a "standalone" version of the code. By standalone, we mean a version of the C5 process code that includes all the physics and mathematics of the operational aerosol retrieval, without the subroutines related to cloud screening, pixel selection, or any of the other MODIS processing overhead (making HDF files, etc.). For inputs, the standalone code requires only values of the observation geometry (angles, location, month) and measured spectral TOA reflectance. Its outputs include only AOD, AE and FMW, plus values of the 660 nm matching error. Therefore, the standalone version enables us to study the mechanics of the retrieval, and test how changes in assumptions can affect the retrieved parameters. Specifically, the standalone algorithm allows us to easily modify the surface reflectance parameterization and determine how changes affect retrievals of total AOD and AE.

Evaluation of the surface assumptions
The surface parameterization of Levy et al. (2007a), given by Eqs. (3-9), includes a dependence on the NDVI SWIR , (Karnieli et al., 2001). This vegetation index is similar to the standard Normalized Difference Vegetation Index (NDVI; red and NIR; Tucker, 1979) except the use of the two longer wavelengths (2.13 and 1.24 µm) helps to reduce aerosol contamination. Basically, the surface parameterization assumes that the slope 660/2130 increases as a function of surfacegreenness.
To evaluate the quality of the assumed surface properties in the algorithm we used a MODIS based surface albedo climatology created by Moody et al. (2005). The global climatology contains spatially complete white-sky albedo (bihemispherical reflectance) data from the years 2000-2004. Moody et al. have also made a five-year aggregate of high quality MOD43B3 data which was filled using an interpolation technique (accuracy of 3-8%; Moody et al., 2008). This aggregate is assumed to be representative of an average year.
In addition, the data set contains surface albedo statistics (average and standard deviation) sorted by ecosystems. These are available at various resolutions of each 16-day time period's snow-free filled albedo maps for each separate year as well as for the five-year climatology. The statistics are available at various resolutions between 0.5 × 0.5 and 90 × 360 (degrees of latitude × degrees of longitude). For our comparisons we selected the albedo statistics for the five-year climatology with the 60 × 120 resolution. This data provided us averaged albedo values at three wavelengths (0.66, 1.24 and 2.13 µm) categorized by ecotypes at 9 different areas. Combined, these 9 areas cover the whole globe. With the 16-day time resolution we get 23 surface albedo values representing the yearly cycle for each ecotype in each of the 9 areas. The overall absolute accuracy of the MODIS surface albedo data is within 0.05 .
Surface albedo is the ratio of the radiant flux reflected from a unit surface area into the whole hemisphere to the incident radiant flux of hemispherical angular extent (Schaepman-Strub et al., 2006), thus it is not view angle dependent. The reflectance in the MODIS retrieval, however, depends on the solar and viewing geometry. To minimize the discrepancies resulting from the comparison of slightly different quantities we decided to compare ratios instead of absolute values. The selected ratios were the slope NDVI SWIR 660/2130 and NDVI SWIR . These parameters were selected because the slope NDVI SWIR 660/2130 is the most important parameter governing the wavelength dependence of the surface reflectance in the MODIS retrieval and it depends on NDVI SWIR as shown in Eq. 9. The slope 660/2130 also depends on the scattering angle (as can be seen from Eq. 5), however, the effect of the scattering angle is not as large as that of the NDVI SWIR . In addition to the slope 660/2130 , the retrieval algorithm uses a ratio of 470 nm and 660 nm surface reflectance (slope 470/660 ). However, this second slope does not have a strong effect on the retrieved fine model weighting because Rayleigh optical depth dominates the signal at the shorter wavelengths.

Analysis with the standalone algorithm
We carried out a number of sensitivity studies with the stand-alone version of the MODIS aerosol retrieval algorithm. We modified the surface reflectance parameters (slope 660/2130 and slope 470/660 ) to check how they affect the aerosol model combination (i.e. FMW), AE, and overall AOD. The FMW parameter was most sensitive to the changes in the slope 660/2130 values. Because the slope 660/2130 depends mainly on the slope NDVI SWIR 660/2130 , we decided to concentrate on it. Figure 5  from retrievals where the slope NDVI SWIR 660/2130 has been multiplied by 0.9 and 1.1, respectively. Modification of the slope does not affect the surface reflectance at 2.13 µm, whereas an increase of 10% in the slope increases surface reflectance at 660 nm and 470 nm by 10%. The errorbars at 660 nm show how much combining of the fine and dust models affects the results. The top of the bar indicates FMW of zero (dust model) while the bottom of the bar is for FMW of one (fine model). In this case, the original retrieval resulted in the best fit with the unmixed dust model because the measured TOA reflectance at 660 nm was larger than the calculated one. When the slope NDVI SWIR 660/2130 was 10% larger than the original, the calculated TOA reflectance at 660 nm became larger and closer to the measured value, thus changing the selected aerosol model from dust to fine. On the other hand, if the slope NDVI SWIR 660/2130 was 10% smaller than the original, the calculated TOA reflectance at 660 nm moved further away from the measured value without affecting the aerosol type selection.
As illustrated in Figure 5, the combination of aerosol models (and therefore AE) is sensitive to the value of the slope NDVI SWIR 660/2130 , therefore it is extremely important that the slope value represents the surface under investigation very accurately. We made a sensitivity study which showed that in order to have less than 20% error in the retrieved AE value, the error in the slope NDVI SWIR 660/2130 has to be less than 10% when AOD is between 0.2 and 0.4. For larger AODs the error in the AE is smaller but for smaller AODs it can be significantly higher. The error of 20% in AE means that AE of 0.6 could vary between 0.48 and 0.72, and AE of 1.5 could vary between 1.2 and 1.8. This accuracy would still enable the separation between fine and coarse aerosol particles.
To study the quality of the slope-NDVI relationship (slope NDVI SWIR 660/2130 in Eq. 9) in the algorithm, we compared it with the MODIS surface albedo climatology. Figure 6 presents the dependence of the slope NDVI SWIR 660/2130 on NDVI SWIR based on global MODIS surface albedo data. The different ecotypes are separated by colors and symbols. In the surface albedo climatology, the globe was divided into 9 areas (60 × 120 degrees each) and averaged surface albedo values for all ecotypes in these areas were calculated. Each ecotype from each of the 9 areas has 23 points representing the yearly cycle of the vegetation. The solid black line in each of these figures shows how the relationship of the parameters is taken into account in the MODIS aerosol retrieval. The albedo data in Figure 6 show that the relationship between NDVI SWIR and slope NDVI SWIR 660/2130 has a large variation between different ecotypes. However, for all ecotypes the slope NDVI SWIR 660/2130 decreases as NDVI SWIR increases which is opposite to the relationship in the retrieval. This behavior was confirmed with linear regression models which all had negative slope values for each ecotype. In addition, the retrieval assumes smaller slope NDVI SWIR 660/2130 values than the surface albedo data shows. To test how the AE and AOD results are affected if we assume the slope-NDVI relationship to behave more like the albedo data suggest, we applied, as a first try, an inverse relationship that fits the data slightly better (shown in Fig. 6 with dashed lines). As one can see, the inverse relationship fits well to number of ecotypes, e.g. savanna and shrubland. However, it is not able to represent, for example, urban, wetland and needleleaf ecotypes. In an attempt to evaluate the MODIS AE quality we used AERONET AE data. The instruments were in agreement on the particle size if both AEs were either larger than 1.4 or smaller than 0.9. In other cases the instruments disagreed. The measurements with AEs between 0.9 and 1.4 were left out of the analysis to avoid the difficult comparison of mixed aerosols. In this study, we compared single measurements, whose time difference was less than 10 min. The number of usable measurement (Table 3) varies for the modified retrievals due to the "AE buffer zone" (from 0.9 to 1.4) and the exclusion of retrievals with AOD less than 0.3 or negative surface reflectance at 2.13 µm. The comparison results are presented in Fig. 7 Table 3, which present the AE and AOD agreement parameters for 6 sites. The parameter AE agreement refers to the fraction of AE agreement while mean(gb-M) is the averaged AOD difference between the instruments and R 2 is the correlation coefficient for AOD. Retrievals done with the original slope-NDVI relationship as used in the C5 MODIS aerosol retrieval are called "Original" while the retrievals done with the inverse of that relationship are called "Inverse". Both of these main retrievals have subretrievals where the fine aerosol model was assigned to be either strongly absorbing (strong, SSA = 0.86) or weakly absorbing (weak, SSA = 0.95) regardless of the fine model originally assigned for the location. Table 3 and Fig. 7 basically show how the selection of the slope-NDVI relationship and fine aerosol model affect the AE and AOD agreement at each site. Figure 7 and Table 3 show that at all studied sites (Alta Floresta, Kanpur, Mexico City, Mongu, Rome, Xianghe) the modification of the slope-NDVI relationship improved the AE agreement. For example, in Alta Floresta the original AE agreement fraction is 0.83 and with the inverse relationship the AE agreement fraction increases to 0.92. When the inverse relationship is combined with the weakly absorbing model, the AE agreement increases to 0.94. Unfortunately, the AOD correspondence at the studied sites does not improve as clearly. The R 2 values between MODIS and AERONET AOD (R2 in Fig. 7) only improve at Alta Floresta when using the inverse slope-NDVI relationship. In Kanpur and Mexico City the use of the inverse relationship improves the correlation but even better correlation is achieved by merely using a different fine model. However, the averaged AOD difference (mean(gb-M)) decreases at all sites, except in Mongu and Xianghe. In Mongu and Xianghe all three parameters give mixed results. The best AOD correspondence is achieved with the original slope-NDVI relationship while the best AE agreement is achieved with the inverse relationship. In Alta Floresta, AOD and AE correspondence can be improved at the same time by using the inverse of the slope and weakly absorbing model. The results for Rome are similar as for Alta Floresta, although the R 2 value for the Inverse(weak) case is not the best possible. In Kanpur and Mexico the two AOD correspondence parameters (R 2 and mean(gb-M)) give mixed results. By using the inverted slope-NDVI relationship and strongly absorbing model in Kanpur we get the best values for AE agreement and mean AOD difference, but R 2 indicates the best AOD correspondence with the original slope-NDVI relationship and absorbing model. Jethva et al. (2010) also concluded that the MODIS retrievals over India could be improved by using a more absorbing fine model and increased surface reflectance in the visible channels. These results show that the absorption of the fine model does not affect the AE agreement, however, it has a significant effect on the AOD correspondence at some sites, like Alta Floresta and Xianghe. The results in Mexico City and Rome are poor in all the cases because they are urban sites with bright surfaces. In addition, Mexico City is an elevated site which also increases the uncertainty in the MODIS retrieval . Therefore, we studied how much the slope NDVI SWIR 660/2130 should be increased at these two sites in order to improve the AE and AOD retrievals. We did the retrievals several times and for each new run we increased the slope by 10%. Then we compared these new retrieval values with AERONET measurements and found Table 3. Comparison of the AE and AOD agreement parameters from the original retrieval (Original), and the retrieval with the inverse NDVI-slope relationship (Inverse) at 6 sites. In addition, retrievals where aerosol models were assigned to be either strongly absorbing (strong, SSA 550 nm = 0.86) or weakly absorbing (weak, SSA 550 nm = 0.95) are presented for both classes. AE agreement refers to the fraction of AE agreement while mean(gb-M) is the averaged AOD difference between the instruments (AERONET-MODIS). R 2 is the correlation coefficient for AOD. These values are also presented in Fig. 6. The best values for each of the parameters are shown in bold. In addition, the number of measurements used in the comparisons are presented. 660/2130 should be as much as 1.6 times larger than the original to introduce significant improvements into the AE (AE agreement fraction 0.80) and AOD (R 2 = 0.56 and mean(gb-M) = −0.03) agreements. In Rome, the best AE (AE agreement fraction 0.95) and AOD agreement (R 2 = 0.57 and mean(gb-M) = −0.01) is achieved with 30% larger slope NDVI SWIR 660/2130 . Similar results were presented by de Almeida et al. (2007) who analyzed the effect of the surface reflectance ratio between the 660 nm and 2.13 µm on AOD retrievals over Mexico City. They found that by increasing the 660/2130 ratio from approximately 0.55 to 0.73 the agreement between the MODIS and sun photometer AOD improved significantly.  Fig. 7. Comparison of the AE and AOD agreement parameters from the original retrieval (Original), and the retrieval with the inverse NDVI-slope relationship (Inverse) at 6 sites. In addition, retrievals where aerosol models were assigned to be either a strongly absorbing (strong, SSA 550nm =0.86) or weakly absorbing (weak, SSA 550nm =0.95) are presented for both classes.

and in
Colors of the bars refer to these different retrievals. AE agreement (top panel) refers to the fraction of AE agreement while mean(gb-M) (middle panel) is the averaged AOD difference between the instruments and R2 (bottom panel) is the correlation coefficient for AOD. These values are also presented in Table 3. 32 Fig. 7. Comparison of the AE and AOD agreement parameters from the original retrieval (Original), and the retrieval with the inverse NDVI-slope relationship (Inverse) at 6 sites. In addition, retrievals where aerosol models were assigned to be either a strongly absorbing (strong, SSA 550 nm = 0.86) or weakly absorbing (weak, SSA 550 nm = 0.95) are presented for both classes. Colors of the bars refer to these different retrievals. AE agreement (top panel) refers to the fraction of AE agreement while mean(gb-M) (middle panel) is the averaged AOD difference between the instruments and R2 (bottom panel) is the correlation coefficient for AOD. These values are also presented in Table 3. correlation between MODIS and ground based AOD at the sites where no significant amount of dust is expected. This happened especially when AOD values were small. According to Levy et al. (2010) MODIS does not provide quantitative information about aerosol size (AE) over land. However, it is an important parameter if one tries to estimate the anthropogenic component of aerosols.
Our findings are in agreement with the recent studies presented by Levy et al. (2010) and Jethva et al. (2010). In contrast to their studies, we also analyzed single MODIS retrievals and therefore were able to study the combining of the fine and dust models in more detail. When employing the spatio-technique of Ichoku et al. (2002), it is harder to see what happens in the combining of the aerosol models due to the averaging of at least 5 measurements. Single measurements, on the other hand, show clearly that the retrieval algorithm usually selects either fine model weighting of 1 or 0. Because AE depends only on the combination of the models, the AE retrieval fails if the combining fails. Our study suggests that the surface reflectance assumptions, especially the slope NDVI SWIR 660/2130 in the MODIS algorithm is the major reason for the inaccurate AE values and the flawed aerosol model combining in the algorithm. Comparison of the surface parameters shows that the slope-NDVI relationship used in the MODIS retrieval algorithm is not supported by the MODIS albedo climatology: 1. the albedo data show that the slope-NDVI relationships depend strongly on ecotypes. The possibility to take this into account in future Collections should be studied further.
2. the slope NDVI SWIR 660/2130 values based on the climatology were larger than the ones in the algorithm.
3. the albedo data also show that the slope NDVI SWIR 660/2130 decreases as the NDVI increases for all ecotypes which is opposite to the relationship in the algorithm. 4. changing the slope-NDVI relationship in the algorithm to agree better with the albedo data (inverse relationship) increased the correspondence of MODIS and AERONET AE at the studied sites. AOD correlation did not improve as clearly, but the mean difference decreased at all sites, except in Mongu and Xianghe.
5. the uncertainty of the fine model's SSA does not have a significant effect on the AE agreement between MODIS and AERONET, however, it affects the AOD correspondence.
Although, the modification of the slope NDVI SWIR 660/2130 improved the AE retrieval at all studied sites and AOD retrieval at some sites, the effect on global scale is not yet known. This will be studied further in the future.