Interactive comment on “ An overland aerosol optical depth data set for data assimilation by filtering , correction , and aggregation of MODIS Collection 5 optical depth retrievals ”

optical depth data set for data assimilation by filtering, correction, and aggregation of MODIS Collection 5 optical depth retrievals” by E. J. Hyer et al. Anonymous Referee #1 Received and published: 11 November 2010 1. This is a very ambitious paper that develops corrections to the MODIS AOD product relative to AERONET over land, with the aim of providing better correlation with the sun photometer observations. This approach offers a good, practical way to produce an “assimilation quality” global satellite AOD product. This study is very carefully done, and demonstrates good understanding of the products used. It represents an important contribution to the literature.


Introduction
The Moderate Resolution Imaging Spectroradiometer (MODIS) instruments on board the Terra and Aqua platforms provide nearly daily global coverage of key atmospheric and land surface parameters.Retrievals of aerosol optical depth (AOD) by MODIS are the most commonly used of any satellite AOD product.Although MODIS is technically a research instrument, its use in operational applications is increasingly widespread.Remote sensing and modeling technology has progressed to a point where operational aerosol data assimilation methods are being implemented for forecasting applications (e.g., Zhang et al., 2008;Hollingsworth et al., 2008).
The use, aggregation, and statistical reduction of remotely sensed aerosol data are application specific.Crucial to the application of AOD data in a data assimilation system is a thorough understanding of the data's error characteristics.Not only is a global estimate of measurement uncertainty needed, but also a point by point uncertainty determination gained through a thorough understanding of retrieval processes that can cause both random and systematic errors.AOD data with consistent biases or poorly characterized uncertainties can degrade analyses and forecasts if used in a data assimilation system (Reid et al., 2009).
For over ocean data, Zhang and Reid (2006) identified MODIS MOD04/MYD04 retrieval bias from sources including lower boundary condition (e.g., white-capping and glint), cloud contamination, and microphysics.From these analyses, a series of quality control/assurance procedures and empirical corrections were applied to devise a data assimilation grade Level 3 product.These data are now assimilated operationally at the Fleet Numerical Meteorology and Published by Copernicus Publications on behalf of the European Geosciences Union.
Oceanography Center (FNMOC) through the Navy Variational Data Assimilation System for AOD (NAVDAS-AOD) into the Navy Aerosol Analysis and Prediction System (NAAPS) (Zhang et al., 2008).
Assimilation of over-land AOD data, however, is a very different problem, both in terms of the retrieval of AOD from satellite observations and assimilation into an atmospheric model.Even without systematic bias, the precision of retrieved AOD over land is lower than over the ocean.The land surface has strong aerosol sources, and observations in proximity to these sources will have stronger spatial and temporal gradients.These gradients reduce the accuracy and precision of AOD retrievals, and interfere with pair-wise comparison to point validation datasets.But the principal challenge of aerosol retrieval over the land surface is the surface itself.Whereas the ocean surface is dark and well characterized, the land surface is brighter and more heterogeneous with strong temporal variability.
The strongest aerosol signal in visible wavelengths comes from scattering of incoming solar radiation.Inverting to obtain the signal of aerosol effects from reflected light requires simultaneous retrieval of the surface reflectance.As surface brightness increases, the relative contribution of the atmospheric radiance decreases, leading to a differential signalto-noise limitation of the atmospheric component.Thus, retrieval over a brighter surface requires a tighter constraint on surface reflectance to achieve a comparably precise retrieval of AOD.The complexity and variability of the lower boundary condition make over-land retrieval of AOD with passive optical observations very challenging.
The MODIS Collection 5 over-land aerosol retrieval (Levy et al., 2007b;Remer et al., 2006) corrects many of the shortcomings of earlier versions, and now represents the state of the art for global retrievals of AOD.It is made publically available in near real time through the NASA Land-Atmosphere Near real-time Capability for EOS (LANCE) (http://lance.nasa.gov).Based on a preliminary analysis, Collection 5 reached the minimum efficacy requirements to be considered for inclusion into the US Navy's operational aerosol modeling.Preliminary studies without major correction have shown that the inclusion of MODIS over-land data can improve model scores in some regions (Reid et al., 2009).
This study expands the Navy's over ocean data assimilation quality dataset (Zhang and Reid, 2006) to include a data assimilation quality Level 3 over land data product based on data collection 5 for ultimate use in NAVDAS-AOD (Zhang et al., 2008).This study is concerned with identifying the various sources of random and systematic error in retrieval of AOD, developing filters and corrections to improve data quality when possible, and characterizing residual uncertainty for use in aerosol data assimilation systems.
The center of the current study is an evaluation of the MODIS Collection 5 data product using 4 years of Aerosol Robotic Network (AERONET, Holben et al., 1998) Sun photometer data and available ancillary data sources.Random and systematic errors in the MODIS Collection 5 AOD retrieval are examined in detail, with a detailed consideration of both global and regionally specific sources of error.Corrections are developed and evaluated for errors caused by the parameterization of the surface boundary condition and by regional microphysical bias.Based on this analysis, a data assimilation ready level 3 product and error model for operational use is presented as the outcome of this study.

Data and methods
The methods and rationales for this analysis are similar to the over-ocean study of Zhang and Reid (2006), investigating lower boundary condition, microphysical and cloud bias.However, due to differences between the ocean and land problems two important differences dictate modifications to the protocol.First, the MODIS over-land aerosol retrieval yields very little information about aerosol size, even compared with the limited coarse-mode/fine-mode fractionation derived by the over-ocean retrieval (Kahn et al., 2009).Second, surface conditions affecting AOD retrieval will co-vary geographically with aerosol properties.For instance, high aerosol loads over brighter surfaces (arid or barren ecosystems) will tend to be made up of coarser particles (dust).These differences require some simplifications to be made to the analysis, and the resulting QA/QC system.
Similar to older versions, the Collection 5 algorithm uses the MODIS 2.12 µm reflectance to estimate surface reflectance at 0.47 and 0.66 µm, and retrieves aerosol properties by inverting the observed reflectance in the visible wavelengths.The Collection 5 retrieval uses a linear model of the relationship between infrared and visible surface reflectances, informed by two empirical relationships estimated from a database of atmospherically corrected MODIS surface reflectance data from cloud-free areas with τ M < 0.2 (Levy et al., 2007b).The first correction relates the 0.66 µm/2.12 µm regression line to a reflectancebased Normalized Difference Vegetation Index of vegetation condition, NDVI SWIR , calculated as NDVI SWIR = (ρ 1.24 µm − ρ 2.12 µm )/(ρ 1.240 µm + ρ 2.12 µm ), chosen for its relatively greater sensitivity to vegetation condition and lower sensitivity to atmospheric conditions (including aerosol particles).The second correction modifies the slope and intercept estimated from the first relationship by as much as ±50% according to scattering angle, based on the observations of Remer et al. (2001).Note that because the information used in the MODIS AOD retrieval all comes from the same viewing geometry, the phenomenon corrected by the scattering angle dependence is the relative anisotropy between the visible and near-infrared reflectances.
MODIS Collection 5 products have a ground footprint spatial resolution of 10 × 10 km at nadir, increasing to more than 20 × 40 km at the edge of the swath.Datasets are divided into 5-min granules covering approximately 2330 km across the satellite ground track and 2030 km along the ground track.Data used in this study are the 0.55 µm "Corrected Optical Depth -Land" from the Level 2 product (hereafter τ M ).One year of data from one sensor is approximately 105 000 Level 2 granules.
The MODIS Level 2 aerosol product includes numerous pieces of auxiliary information about the retrieval conditions.In this study, the following were used: (1) the MODIS Mandatory quality assurance (QA) flag, which assigns each retrieval an estimated quality of "Bad", "Marginal", "Good" or "Very Good," (2) the cloud fraction information, indicating the fraction of pixels within the retrieval footprint with MODIS-detected cloud, (3) viewing angle and (4) the scattering angle for each retrieval.
For the purposes of data assimilation MODIS products from Terra and Aqua behave similarly, and in this paper joint statistics are sometimes presented.When potentially significant differences between Terra and Aqua are diagnosed, separate statics are presented and discussed.Appendix D presents discussion and analysis specifically targeting differences between Terra and Aqua MODIS AOD.
Note that the MODIS retrieval algorithm permits negative values of retrieved AOD, and in certain locations this condition can be quite common (e.g.forested area in South America, see Sect.4.3).Our statistical analysis of the MODIS Level 2 data uses these negative retrievals.The aggregated product produced after correction of the Level 2 data does not include negative AODs; if the AOD is negative after aggregation, it is truncated to 0.

AERONET sun photometer data
The basis for our evaluation of MODIS retrieved AOD is direct measurements of AOD from the Aerosol Robotic Network (AERONET) (Holben et al., 1998) of Sun photometer instruments.This study uses AERONET Level 2 quality controlled data collected throughout the entire network for 2005-2008.Quality control procedures are as described in Smirnov et al. (2000).Depending on the exact instrument used, there are a variety of wavelengths collected ranging from 0.34-1.6µm.AOD at the MODIS retrieval wavelength of 0.55 µm (hereafter τ A ) was derived from the quadratic interpolation method of O'Neill et al. (2003).Typically AERONET level 2 AOD data have uncertainties of < ±0.015 (Eck et al., 1999;Schmid et al., 2003) spectral deconvolution algorithm and are used in this analysis.
MOD04 retrievals were matched to AERONET measurements of AOD with tolerances of ±30 min and ±30 km.Errors resulting from transport or localization conditions were found to be small at the broad statistical level (see Appendix A).Every possible match between instantaneous AERONET measurements and MODIS retrievals was included in the matched dataset.Because many AERONET stations have a high sampling rate (15 min or less), individual MODIS retrievals are often paired with multiple AERONET retrievals.For the final step of estimating instrument error variance for data assimilation, as described in Sect.7, AERONET data from each site were aggregated into six-hour bins.
The AERONET Level 2.0 algorithm also includes steps to remove cloud contamination, creating a potential sampling bias in our analysis of cloud effects in the MODIS retrieval.This is discussed further in Sect.3.2.

MODIS albedo and snow products
The lower boundary condition is an important component of the aerosol retrieval, and the complexity and variability of the land surface pose significant challenges for accurate retrievals.This study is concerned with two potential issues: (1) Surface albedo characterization and in particular how it relates to the IR-visible regression, and (2) potential for snow contamination.Surface properties were characterized for this study using the MODIS MCD43 albedo product (Schaaf et al., 2002).This product is produced as a Level 3 16-day composite product based on the MOD09 atmospherically corrected surface reflectance product (Vermote et al., 1997).The 16-day compositing algorithm is designed to systematically eliminate the influence of clouds and aerosols on the observations.The algorithm uses observations from multiple viewing geometries to estimate the bidirectional reflectance function (BRDF), by fitting observations to a RossThick-LiSparse model of surface albedo and BRDF (Lucht et al., 2000).The model parameters can be used to calculate black-and white-sky albedoes for the seven primary MODIS bands, as well as integrated visible, nearinfrared, and shortwave albedoes.This study used the blacksky hemispheric albedoes at 0.47 µm, 0.66 µm, and 2.12 µm from the MODIS albedo product (product MCD43C3), calculated using the mean solar zenith angle for each location and time.The QA information included with the MODIS albedo product was used to exclude all albedo data with quality other than "very good." For snow, the MODIS albedo product incorporates tests for snow contamination in each 500-m MODIS pixel that www.atmos-meas-tech.net/4/379/2011/Atmos.Meas. Tech., 4, 379-408, 2011 goes into the MCD43 albedo product.The MCD43 product calculates albedo using only snow-free pixels, and includes a diagnostic variable indicating the fraction of 500 m pixels in each 0.05 • cell where snow was indicated.This information was used to diagnose residual snow-related bias in the MODIS AOD.
Because the footprint of MODIS aerosol retrievals has a highly variable size, shape, and orientation, no attempt was made to match MODIS AOD retrievals precisely to the Level 3 snow/albedo data.Instead, the center of each AOD retrieval was matched to the nearest 0.05 • grid cell from the MCD43C3 climate modeling grid (CMG) albedo product.For the snow analysis, a spatial and temporal window around the matched observation was used (see details in Sect.3.3).
Since both the AOD and snow/albedo products are derived from the same MODIS observations, there is the possibility that the interaction between AOD and albedo is reciprocal.Appendix B to this paper includes an analysis to determine whether MODIS albedo data were contaminated by aerosol.Results in Appendix B show that persistent high AOD often resulted in failed albedo retrievals, but that AOD contamination in albedo retrievals with "very good" QA flags was negligible.

Statistical evaluation of MODIS AOD
To evaluate the utility of AOD information for use in a data assimilation grade product, standards of accuracy must be established.In terms of interaction with an aerosol transport model, both random and systematic errors will affect analysis and forecast outcomes.Producing a dataset for data assimilation applications involves three goals: (a) elimination of outliers, (b) elimination of systematic bias, (c) quantitative characterization of residual errors.Because of the physics of the retrieval problem, the relative strengths of the land surface, microphysics and cloud biases are AOD dependent.Nominally aerosol regimes can be broken down into: <0.2, low AOD near signal to noise thresholds; 0.2-0.6,sufficient aerosol signal in the single scattering regime; 0.6-1.4multiple scattering regime; and >1.4 extreme AOD events.These categories are used as necessary in the statistical analysis.
A simple intuitive metric for analysis of the retrieved AOD values is the tolerance which bounds AOD errors.For this study, a target accuracy is defined by: and sets of retrievals are evaluated according to the fraction that fall within, above, and below this target accuracy.This is a weaker constraint than the proposed 0.05 ± 0.15 × AOD used in some studies (e.g., Levy et al., 2007bLevy et al., , 2010)); however, those studies apply spatial and temporal averaging to the MODIS and AERONET AOD data before calculation of errors.This metric provides a simple accounting for the frequency of positive and negative outliers, but provides little information about the magnitude of errors.Where appropriate, cumulative distribution functions (CDF) of errors are used for a more quantitative description.
Positive and negative errors in retrieved AOD are often asymmetrical, which violates a key assumption of ordinary least-squares regression.Linear fits of MODIS and AERONET AOD consistently have intercepts significantly different from zero.These non-zero intercepts are driven by errors in retrieved AOD in relatively clean conditions, due to problems with the surface boundary condition, cloud contamination, or other problems.Because the error budget of retrieved AOD varies as a function of AOD, the nonzero intercepts often artificially skew the regression slope away from the relationship indicated by the high-AOD data.Therefore, to derive a linear representation of the relationship between AOD values, the slope forced through zero is calculated using the following equation: Slope calculations in this analysis are made using only moderate to high AOD (0.2 < τ A < 1.4).Where data volume is sufficient, a separate slope calculation is made for extreme AOD (τ A > 1.4).
Finally, this analysis makes extensive use of the Root Mean Squared Error (RMSE), calculated as: RMSE is sensitive to both systematic and random errors and provides an estimate of the expected error in AOD in the absence of information about specific sources of errors.In nonbackground aerosol conditions, RMSE has a strong relationship with AOD (see Sect. 3.1).For this analysis, the background observations (τ AERONET < 0.2) are used to calculate a "noise floor" for RMSE, and then a linear relationship of the form is estimated for higher AOD.Prognostically, "estimated RMSE" (ε) is calculated for a MODIS observation of a given optical depth as the higher of the regressed value or the RMSE "noise floor": When data volumes are sufficient, this approach is extended to use a separate linear relationship for very high AOD values.

Aggregation and textural filtering
For assimilation into an aerosol model, AOD observations are aggregated into a gridded Level 3 product for assimilation.This is done for two reasons: first, to avoid artifacts caused by assimilation of subgrid features the model cannot accurately resolve, and second, to reduce random error through averaging.The method of Zhang and Reid (2006) is used to aggregate the MODIS AOD to a resolution of 1 • by 6 h.Three textural checks are used to avoid assimilation of sub-grid features likely to create anomalous features in the 1 • model, and also to reduce residual cloud contamination.The first is a "buddy check" in the aggregation to exclude AOD retrievals with no adjacent retrievals.The second sets a minimum number of retrievals per grid cell.Zhang and Reid (2006) imposed a minimum of 5 retrievals per 1 • grid cell for the overocean data.Results from a compliance test to evaluate data loss vs. quality improvement indicated that requiring more than 3 retrievals per grid cell resulted in substantial data loss with minimal improvements to error statistics (3 AOD retrievals = nominally 300 km 2 = 2.5% of the area of a 1 • cell at the equator, in the worst case).The final textural filter uses a local variance check to eliminate all grid cells with mean τ M > 0.2 and coefficient of variation of AOD within the grid cell (standard deviation/mean) greater than 0.5.The effects of these textural filters are described and discussed in Sect.7.3.

Results 1: global error terms for MODIS Collection 5 aerosol product
Figure 1a shows the global extent of coverage of the MODIS Collection 5 AOD retrieval over land, with the AERONET stations used in this study marked with "+" and overlaid on the annual average AOD for the years 2001-2009.Also overlaid are the boundaries defining regions for analyses that are presented in subsequent sections.This analysis uses over 3.8 million data points from more than 290 sites.Figure 1b gives the subsequent global scatter/density plot.Note that both axes of this plot actually begin at −0.05; as discussed in Sect.2.1, negative AOD values can be returned by the MODIS retrieval, and are not truncated in our analysis of the Level 2 data.Each vertex of the black solid line represents 50 000 AOD retrievals, binned in order of τ A .The black dotted curves represent the 25th and 75th percentile τ M of each τ A bin.The green solid lines indicate the tolerances from Eq. ( 1).The slope of τ M vs. τ A is near unity, although a small cloud of high biased AOD retrievals are visible for AOD < 0.2 and a clear upward break is discernable for AOD greater than 1.4.Note that representativeness error should result in negative bias at high values of τ A , so the positive bias in τ M at extreme τ A is likely larger than it appears from this figure (also see Appendix A).
Global error statistics for MODIS retrieved AOD from this data segregated by MODIS QA flag values are given in Table 1a-c.Table 1a gives the distribution of retrieved τ M for each subset, regression slopes calculated using Eq. ( 2), and regression coefficients r 2 .fractions relative to the criterion of Eq. (1).Table 1c gives the "noise floor" RMSE estimated using only retrievals with τ A < 0.2, and the diagnostic and prognostic regressions of RMS error vs. AOD.As expected, higher QA flags have better performance with compliance and r 2 values steadily increasing with data quality.Key to data assimilation is a quantitative estimate of the uncertainty of the data.RMSE is plotted as a function of diagnostic AERONET (RMSE vs. τ A ) and prognostic MODIS AOD (RMSE vs. τ M ) in Fig. 2a and b, respectively.This is followed by a CDF of MODIS bias in Fig. 2c is seen to be highly linear with AOD, except for a minimum at around 0.05 below which RMSE actually increases with decreasing AOD.This linear relationship of Fig. 2b is the basis of a prognostic estimate of error for use in data assimilation.For the region of low AOD where RMSE is non-linear, a "noise floor" uncertainty is assigned using the RMSE for all cases where τ A < 0.2.For a given value of τ M , the RMSE is estimated as the largest of the estimates from linear regression and the "noise floor" RMSE.The gray lines in Fig. 2b illustrate the RMSE estimated by this approach for τ M with "very good" QA.
The baseline RMSE reaches a high of 0.15 for the lowest QA value ("Bad"), and improves to as low as 0.08 and 0.07 for 'Very Good' Terra and Aqua respectively.RMS error increases with τ A with slopes ranging from 0.22 to 0.27, and generally decreasing with stricter QA.Retrievals with very high optical depths (AOD > 1.4) appear to have different error characteristics, with larger errors in most cases.
Figure 2c shows cumulative distribution functions (CDFs) of fractional error, defined as (τ M − τ A )/τ A , for QA = "very good" only.At each range, the median error indicates a small negative bias in the global data.Note that for "Very Good" data, negative errors are more common at low AODs than positive errors, in contrast with other QA flags (Table 1b).When all retrievals are included, median biases are closer to zero over all ranges (data not shown).This indicates that the "Very Good" QA flag removes more positive errors than negative errors.The surplus of negative errors at low AOD will be shown to be partly related to problems with the surface boundary condition in the retrieval.
Each of the statistical categories in Table 1 indicates that τ M retrieved from MODIS-Aqua are slightly higher than from MODIS-Terra.This can also be seen in the relationships shown in Fig. 2.This difference is generally very small, but extremely persistent.Note that higher AOD from MODIS-Aqua does not always mean poorer error statistics; at low optical depth, the "Very Good" retrievals from both sensors have a negative bias, and MODIS-Aqua has better compliance in these cases.
The compliance statistics in Table 1c indicate that the algorithm performs slightly better in the 0.2-0.6AOD range compared with higher AOD ranges.In the range of 0.2-0.6 aerosol signal is good, and radiative transfer is in a near-linear single scattering regime.At higher optical depths, multiple scattering becomes increasingly important and leads to amplifications of aerosol microphysical biases.Scale effects also become important at high AOD, because the spatial extent of high-AOD features is frequently smaller or of the same order of magnitude as the comparison here (30 km radius around each AERONET site).
The summary global statistics presented here are in general agreement with the findings of other comparisons (Levy et al., 2007a(Levy et al., ,b, 2010;;Remer et al., 2005), and also demonstrate the importance and utility of an AOD specific error metric such as RMSE that is not emphasized in previous analyses.The Data Quality Statement for the MODIS Collection 5 Aerosol product (ftp://ladsweb.nascom.nasa.gov/allData/5/MOD04 L2/README) states that only data flagged "Very Good" should be used for scientific analysis, and the results here corroborate that these data provide systematically better agreement with AERONET.Unless explicitly stated, the remaining analysis in this manuscript uses only data flagged "Very Good."

Factors affecting global signal-to-noise ratio of retrieved AOD
For an atmospheric modeling application, systematic errors are often more damaging than random errors.In the case of a global model, the model resolution is typically coarser than the AOD observations, permitting extensive averaging to reduce random noise.High levels of random noise can also obscure systematic bias.The dominant signal of aerosols for optical retrieval is the scattering of direct sunlight back to the sensor.The MODIS aerosol retrieval is based on inversion of a full radiative transfer simulation that includes aerosol effects on all pathways except multiple scattering (Levy et al., 2007a), but the primary determinant of the signal strength for retrieval of AOD is the ratio of scattered direct sunlight to light reflected from the surface.This is a reason why retrieval of AOD over ocean is more precise than over land, and also a reason why this MODIS algorithm is unable to retrieve aerosol properties over bright surfaces.
Figure 3 shows how the compliance of the MODIS retrieved AOD varies as a function of surface albedo at 0.47, 0.66, and 2.12 µm.The gray bars at the top and bottom of the graph indicate the fraction of retrievals above and below the target accuracy.Thus, the area between the shaded regions is indicative of the compliant fraction of the data.The solid lines (each vertex is an average of 20 000 retrievals) indicate the mean AOD bias, and the dotted lines indicate the 25th and 75th percentiles of each bin.The graphs in Fig. 3 illustrate two effects.First, the increase in MODIS-AERONET bias with increasing albedo at 0.47 µm and 0.66 µm clearly points to mis-parameterization of the lower boundary condition in the retrieval.This is examined in detail in Sect. 4. Second, the increased spread of the errors at higher albedo illustrates the decrease in differential signal to noise ratio.
Figure 4a illustrates the differential signal-to-noise consideration for AOD retrievals through the optical path length of the atmosphere.Compliance statistics and mean AOD bias are shown as a function of the sensor scan angle (angle from nadir), which determines the atmospheric path length Table 1c.Prognostic and diagnostic regression of RMS error in MODIS AOD as a function of AOD, stratified by MODIS QA value.A single RMS error estimate for low-AOD conditions (τ A < 0.2) is also shown.).This means that comparisons of MODIS aerosol retrievals with narrowswath instruments such as the Multi-Angle Scanning Radiometer (MISR) or the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) will overestimate random error against the whole MODIS product.This discrepancy between nadir and scan edge is caused by the relative contribution of surface reflected light to the total radiance at the sensor.Thus, it diminishes at increasing optical depth: for τ A > 0.6, compliance fractions are 60% and 61% at nadir and scan edge, respectively.At even higher τ A , spatial mismatch between MODIS and AERONET is a larger factor; for τ A > 1.0, compliance is better at nadir.A related phenomenon is shown in Fig. 4b, which depicts the bias and compliance statistics as a function of the scattering angle.In the global dataset the retrieval has almost no systematic bias associated with scattering angle, indicating that the model used to account for anisotropy in the surface reflectance (Levy et al., 2007b) appears to be sufficient.But unavoidably, as scattering angle increases, shadows over vegetated surfaces diminish, the surface brightness increases, and the precision of the retrieval declines.The interaction between solar geometry and scattering angle means that scattering angle distributions are not stationary with latitude, which may also influence this result.At scattering angles smaller than 100 • , 84% of retrievals are compliant.At very high scattering angles, where the sun is almost directly behind the sensor, there is a sharp spike in retrieved AOD: compliance is only 41% for scattering angles above 170 • .This is caused by the "hot spot" of vegetation reflectance (Vermote and Roy, 2002), and retrievals with scattering angles over 170 • (0.5% of our matched retrieval data set) should be avoided because of this problem.

Influence of MODIS-detected clouds on retrieved AOD
Undetected clouds, subpixel or otherwise, can cause a positive bias in the retrieval.Conversely, cloud shadows can result in a negative bias.Overall, however, we expect a predominantly positive bias from cloud effects.The MODIS retrieval includes auxiliary information on the fraction of pixels within the retrieval footprint with MODIS-detected cloud.Only 16% of successful AOD retrievals have MODISdetected cloud cover within the retrieval footprint.This fraction decreases with increasingly strict QA, from 26% of "Bad" retrievals to 10% of "Very Good" retrievals.Retrievals with MODIS-detected cloud have a slight positive bias relative to the complete dataset (Fig. 4c).The difference in mean τ M for retrievals with indicated clouds versus no detected clouds is +0.04, while the corresponding difference in τ A is less than 0.01.While this elevated AOD may indicate undiagnosed subpixel clouds affecting the reflectances used in the retrieval (Zhang and Reid, 2006), or may be an artifact of three-dimensional scattering not included in the retrieval model (Varnai and Marshak, 2009), other studies (Koren et al., 2007) contend that this may be the result of elevated aerosol particle concentrations in the vicinity of clouds.But, given the scale of such aerosol features, we might expect a larger response in τ A .For over water cases, τ A did not show an increase in AOD with increasing cloud cover nearly as large as τ M (Zhang et al., 2005), suggesting artifact may be the more dominant factor.
One important caveat of this analysis is the cloud-clearing undertaken in the processing of the AERONET data.While not 100%, the cloud-clearing algorithms used in AERONET Level 2 processing are highly effective.Because no AERONET retrieval is made in conditions flagged as cloudy, the matched MODIS-AERONET dataset systematically excludes many observations in partially cloudy or thin cirrus cloud conditions.Thus, the true impact of clouds on MODIS Error in τ M as a function of the ratio of MODIS albedo at 0.65 µm and 2.1 µm.The lines were constructed by sorting the albedo of matched retrievals into bins of 20 000, each of which is shown as a vertex of the solid (mean bias) and dotted (25th and 75th percentiles of bias) lines.The gray bars at the top and bottom of the graph represent the fraction above and below, respectively, the compliance tolerances of Eq. ( 1).All statistics were calculated using only matched data with MODIS QA of "Very Good".
retrieved AOD may be more severe than indicated by the matched MODIS-AERONET comparison.
With the exception of the small positive bias, the error statistics for retrievals with MODIS-detected cloud are no worse than for the entire data set (data not shown).Data assimilation studies have shown that in some regions, such as SE Asia, cloud bias can be so severe that it actually decreases model efficacy when assimilated (Reid et al., 2009).Because of the contextual bias associated with using cloudcleared AERONET data for comparison, the extent of cloud contamination issues can only be revealed by data assimilation studies.Further research is needed to determine the extent and localization of cloud contamination, and how best to avoid cloud contamination without unduly reducing data coverage.

Diagnosis of snow contamination in MODIS AOD
Snow cover disrupts the optical retrieval of AOD because the variation in snow reflectance with wavelength is completely different from that of vegetated surfaces, which disrupts the estimation of surface boundary condition at the core of the AOD retrieval.Also, snow has a very high reflectance in visible wavelengths, reducing the differential signal-to-noise ratio of aerosol retrievals.Snow-contaminated pixels are identified and removed from the Collection 5 MODIS AOD retrieval using an empirical thresholding technique based on observed radiance in near-infrared and thermal wavelengths.This technique was shown to be a great improvement over the snow masking used in the Collection 4 MODIS AOD retrieval (Li et al., 2005).Snow contamination is extremely difficult to completely eliminate because snow can cover any fraction of the retrieval footprint.For this study, the snow contamination checks in the MODIS Collection 5 algorithm (Li et al., 2005) were extended using the snow flag in the MODIS 0.05 • albedo product, which indicates the fraction of 500-m pixels affected by snow during the 16-day compositing period of this product.Using this product makes it possible to extend the test for snow contamination both in time, by considering data from a period before the AOD retrieval, and in space, by looking for snow over a wider area.These spatially and temporally extended checks do not give proof of snow contamination of a given retrieval.However, as will be shown, they permit identification and removal of retrievals where the probability of snow contamination is greater.
For AERONET retrievals in Boreal North America (above 49 • N) with τ A < 0.2, compliance is 70%, and the mean bias in the matched dataset is +0.013, or 14%.Of the noncompliant retrievals, 70% are biased high.Note that the matched dataset includes no retrievals from this region during December through February, and only 4% of retrievals are from November through April.This relates to limitations in both datasets: AERONET sites in northern regions are customarily shut down during winter months to prevent damage by snow, and MODIS retrieves AOD only when the sun is less than 72 • from nadir.
Only 0.2% of matched retrievals from this region are flagged as having snow in the corresponding MODIS 16day albedo product, but these flagged points are 59% compliant and have a mean bias of +0.052, and the non-compliant retrievals are 95% biased high.Thus, the snow flag in the MODIS albedo product is consistent with snow contamination in the MODIS AOD, but the number of points flagged is very small, and their effect on the bulk error statistics is minimal.This indicates that this test for snow is effective but likely incomplete.
If we extend the snow test in time by checking the snow flag in the MODIS albedo product for the 32 days prior to the AOD observation, and in space by looking for snow in a 0.35 × 0.35 • box around the AOD retrieval, we capture nearly 12% of the retrievals in the matched dataset, including nearly all the retrievals in March and April.The retrievals captured in this wider net have 59% compliance, a positive bias of +0.031, and 85% of non-compliant retrievals are biased high.Thus, even this larger set of retrievals is consistent with snow contamination.When these retrievals are removed, the compliance of the remaining retrievals in Boreal NA is 73%, 3% higher than before filtering.The mean bias is reduced to +0.009 (9%).Non-compliant retrievals are still 66% biased high.These results indicate that some snow contamination has been removed, but positive biases persist in the data, whether related to snow contamination or other causes.
Table 2 shows statistics for retrievals excluded using the matched and extended filters based on the MODIS albedo product, for all regions that extend to mid-latitudes or above.

Basic QA filtering for subsequent testing
Based on the results in this section, the remainder of our analysis at the retrieval level (Sects.4-6) will exclude the following MODIS AOD data: 1. Data with MODIS mandatory QA other than "Very Good"; 2. Data with MODIS-indicated cloud; 3. Data with scattering angle above 170 • ; This represents the extent of filtering that can be achieved using only the auxiliary data included in the MODIS aerosol product.This level of filtering will serve as the baseline for our evaluation of the additional filtering and correction steps described in the next sections.

Description of regional biases
Global statistics average out considerable regional variability in the retrieval quality, and are not geographically representative because they are weighted to locations with high AERONET coverage.Even at the regional level, the covariance of particle properties with AOD results in mixed individual site efficacy as a function of AOD.Table 3 shows statistics, calculated using QA filtering as described in the previous section, for 14 different land regions.Statistics for each AERONET site for the 2005-2008 study period are provided in the Supplement associated with this paper, or can be sent by request.Higher compliance regions include northern Eurasia, Southern Europe-Mediterranean, and the Eastern Contiguous United States (E CONUS).Low compliance regions include South America, Africa above the equator, Central America, and Peninsular Southeast Asia.Just as important, the balance of positive and negative non-compliance shows large regional variability.In Northern Africa, 89% of non-compliant retrievals are biased low, while in the Indian subcontinent, 86% of non-compliant retrievals are biased high.Regional slope is also variable, ranging from a low of 0.69 for Africa above the equator to 1.17 for North American Boreal (the slope for the Western continental US is 1.31, but the regression is weak, with r 2 = 0.25 for that region).
Just as the global mean statistics are composed of regions with different behavior, so the slopes calculated for each region are the aggregate of the sites within that region, whose behavior is not always homogeneous.Our analysis is necessarily coarse-grained; in addition to variations between sites, aerosol properties in our selected regions will also have seasonal and interannual variation that we do not explicitly address (seasonal breakdowns for each site are included in the Supplement).This section provides a few key pieces of context from our analysis and from the literature that can aid in interpreting the statistics for each region.Regional biases will be analyzed and quantified in the next section.

North America
North America is split between a Northern Boreal Region, Eastern CONUS and Western CONUS.Owing to its background nature, the northern boreal AOD is typically low, lowest of all of the regions compared (mean τ A = 0.11).Individual site performance varies widely.For very clean background sites -e.g., Yellowknife (62 • N, 114 • W), Opal (80 • N, 86 • W) -AOD rarely exceeds 0.15, yielding no information from regressions.Sporadic high AOD events occur in the spring and summer months at some sites -e.g., Bonanza Creek (65 • N, 148 • W), Bratts Lake (50 • N, 105 • W), Pickle Lake (51 • N, 90 • W) -due to boreal biomass burning activity (Eck et al., 2009), which can drive r 2 values to nearly 0.9.In very rare occasions, AOD can be extremely high (>8), overwhelming the AERONET sensor (O'Neill et al., 2002) and driving aerosol layer reflectivity to the AOD "Semi-infinite" regime .In addition to these real events, the boreal is also susceptible to potential snow bias which some sites -e.g., Kelowna (50 • N, 119 • W) -clearly show as strong early season high biases.Regressions are dominated by burning activity, and a fairly consistent high slope bias is present among most sites.These slopes are event driven, and are dependent on the microphysical properties of smoke particles constantly evolving in size and absorption (Reid et al., 2005a,b).As a consequence, some sites show slope biases higher than 1.5 (e.g., Bratts Lake), and others close to 1.0 (e.g., Bonanza Creek).Terra and Aqua appear to perform similarly, with slope deviations between them at ±0.1.For the sites affected by extreme AOD events, slope biases ranged from 1.4 to 1.9, with Aqua slopes higher than Terra by 0.2 or more in all cases.
MODIS performs best in the Eastern CONUS.With the more uniform mix of sulfate and organic pollution, slopes for most sites are within within ±0.2 of unity.Forty percent of sites have over 80% compliance.RMS error is also very good, one of the lowest regionally.

Central America
Central America has only four sites with a significant number of data points of AOD > 0.2 (e.g., >500 data points): Tuxtla Gutierrez (17 La Parguera is an excellent site for monitoring the transport of African dust over the Caribbean.In the matched data set, only MODIS-Aqua has sufficient range of AOD to establish a robust correlation (slope = 0.94, r 2 = 0.63).Compliance for MODIS-Aqua is also much better (72% for all AOD vs 57% for MODIS-Terra).

South America
South America has four key aerosol regimes: The August-October burning in Rhondonia and Mato Grosso Brazil; February-May northern biomass burning in Columbia and Venezuela; an Argentinean dust regime; and urban super plumes from such major cities as San Paulo and Buenos Aires.
At low AOD, compliance is highly variable within the region: for AOD < 0.2, compliance varies from 14% at Cuiaba-Miranda (15 • S, 56 • W) to 52% at Petrolina-SONDA (9 • S, 41 • W) and Santa Cruz -UTEPSA (18 • S, 63 • W).At low AODs, negative retrieved AOD is common.Overall RMSE is good in South America for moderate AOD.The massive Brazilian biomass burning signal is the primary determinant of validation statistics for this region.Each of the 10 sites in the biomass burning region have among the highest r 2 values over the globe, yet each of the sites has a significant slope bias, and these biases vary from site to site; this leads to the overall poor r 2 and RMSE values for the continent.Also, differences between Terra and Aqua derived AOD are highly localized, with Aqua being generally higher.The heavily impacted Ji-Parana site (11 • S, 62 • W) in central Brazil (average AOD > 1.0 for all data where AOD > 0.2), shows the strongest biases, with a moderate-AOD and extreme-AOD slopes of 1.13 and 1.54 for Terra, and 1.28 and 1.93 for Aqua.For other smoke receptor sites, the slopes for Terra vs. Aqua show differences on the order of 0.1, with Aqua consistently higher.The positive bias in τ M is considerably stronger at high AOD; the high-AOD slope is always higher than the moderate-AOD slope (e.g., Fig. 7p where slopes are presented for data where τ M > 1.4).The regional mix of forest, cerrado, and pasture burning coupled with the rapid evolution of the smoke plume makes smoke aerosol properties highly variable (Ferek et al., 1998).Since the MODIS retrieval cannot resolve these variations in aerosol properties, the MODIS microphysical bias is highly variable.
No AERONET sites are available in Northern South America.Southern South America is prone to occasional dust emissions in an otherwise pristine environment (Gasso and Stein, 2007).Background AODs are typically low (<0.1).In such conditions over a semi-arid region, we are near the noise threshold for MODIS, which shows poor skill in general, often overestimating AODs by a factor of three.Compliance at most sites in this region is low and is driven by biases at low optical depth, whether high biases (48% of retrievals at Trelew (43 • S, 65 • W)with τ M < 0.2 are biased high) or low biases (55% of retrievals at CEILAP-UTN (35 • S, 58 • W) with τ M < 0.2 are biased low).
Lastly, we examined the urban regions of San Paulo (Brazil) (24 • S, 47 • W) and La Paz (Bolivia) (17 • S, 68 • W).The La Paz site is extremely clean: in the matched MODIS-AERONET dataset, less than 3% of τ A were above 0.2.The MODIS AOD has no skill in regression against AERONET at this site, and compliance for τ M > 0.2 is poor (<20%), but overall compliance is good because of the dominance of very low AOD.Sao Paulo performs moderately well in both regression and compliance, but results for low AOD are worse than most other stations (34% compliant for τ M < 0.2).

Europe, the Mediterranean, and northern Eurasia
Like Eastern CONUS, there are many AERONET sites in Europe.As expected, urban sites often show poorer performance relative to more remote areas, but in general the fraction of complaint data points is very high.RMSEs are quite good, both in the noise floor (0.06 for Terra and Aqua) and the diagnostic slope (0.03 + 0.17 τ A for Terra, 0.03 + 0. The three sites in equatorial Africa -ICIPE-Mbita (0 • S, 34 • E), Nairobi (1 • S, 37 • E), and Kibale (1 • N, 30 • E), typically exhibited AOD < 0.2.They are typically compliant for 50-80% of retrievals with estimated slopes ranging from 0.8 to 1.1.Given the small dynamic range of data for these sites, however, r 2 values are <0.5.Central Africa hosts the world's largest biomass burning features in the months of July-October.Despite the size of these features, there are very few AERONET sites in the region and it is difficult to determine if there is the same performance heterogeneity as South America.However, there is some advection into southern Africa.
In addition to episodic smoke events from the north, high sulfate pollution is present year around in the South African Highveld and the Johannesburg regions.For biomass burning, Mongu (15 • S, 23 • W) is the only representative site of the central biomass burning region, with moderate r 2 value (∼0.6 for both Terra and Aqua), nearly 50% RMS errors, and a slope of 0.8.Bias between MODIS and AERONET at southern African sites appears to be correlated to particle absorption (Tom Eck, NASA GSFC, personal communication, 2010).In the polluted Skukuza site (25 • S, 32 • E) in the Highveld r 2 values are good (>0.75), but a low bias is also persistent (slope = 0.88).The urbanized Witswatersrand University site in Johannesburg (26 • S, 28 • E) shows positive slope bias 1.32 (r 2 = 0.30) but still exhibits 70% compliance overall.

Indian sub-continent
Over our the study period, the number of AERONET sites around the Indian Sub-continent was not sufficient to determine spatial performance.The major aerosol features of the Indo-Gangenic Plain do appear to be well represented in MODIS AOD data.At the Kanpur (27 • N, 80 • E) and Gandhi College (26 • N, 84 • E) sites correlations are good (r 2 > 0.7), and slopes are within ±0.1 of unity.Compliance is above 77% for both sites and satellites.The limited data (N ∼ 500) at Nainital (29 • N, 79 • E) and Pantnagar (29 • N, 80 • E), right at the base of the Himalayan mounts also performs adequately, with some slope bias indicated (ranging from 0.95 to 1.39), and r 2 values of ∼0.25-0.75.All other sites, including all sites outside of the Indo-Gangenic plane (rest of India, Pakistan) have very limited data points (<200) and very poor scores.It is noteworthy, however, that studies with handheld Sun photometers have been conducted in western India for Aqua (Misra et al., 2008).They showed interseasonal variability and moderate r 2 values (∼0.5) in MODIS scores, with in general a slope of 0.8 (but see discussion of linear regression in Sect.2.4).

East and Southeast Asia
East and Southeast Asia encompasses several of the most complex aerosol and land surface environments.Land surface varies from extreme arid to dense vegetation to urban, with occasional strong pollution, smoke, or dust events aloft.Values for individual site r 2 range from 0.38 -Taipei, Taiwan (25 • N, 122 • E) to 0.84 -Xianghe, China (40 • N, 117 • E), and slopes range from 0.75 (Taipei, Taiwan) to 1.2 -Beijing (40 • N, 116 • E).Regions susceptible to dust, smoke, and sulfate aerosol particles can have considerable scatter and extraction of exact biases is well outside the scope of this paper.Because of the environmental heterogeneity of Asia, regressions can be reasonable, but compliance fraction low, and visa versa.Readers are encouraged to review the individual site statistics given in the Supplement.SE Asia also shows mixed results.For most sites r 2 values are above 0.7, yet compliance is typically less than 60%.Fine mode dominates here, yet slopes ranged from 0.76 -Mukdahan (17 • N, 105 • E) to 1.52 -Bac Lieu (9 • N, 106 • E).Terra and Aqua are reasonably consistent (slope difference <0.2 for most sites).

Australia
The bulk of Australian AERONET sites are in arid regions and typically have very low AODs (<0.2).Consequently, regression statistics are poor, and calculated slopes are not meaningful.Jabiru (13 • S, 133 • E), the only location in the matched data set impacted by biomass burning smoke, is the only site with enough moderate AOD retrievals in the matched dataset for a meaningful comparison (slope near 1.0, r 2 = 0.24, compliance of 56% for τ M > 0.2).Australia exhibits the same pattern seen in other arid landscapes where many sites have predominantly low bias at τ M < 0.2, but almost 100% positive bias for 0.2 < τ M < 0.6, indicating a www.atmos-meas-tech.net/4/379/2011/Atmos.Meas.Tech., 4, 379-408, 2011 boundary condition bias as well as significant positive errors in retrieved AOD under clean conditions.

Diagnosis and correction of errors related to surface boundary condition and microphysics
The site by site biases described in Sect. 4 imply that biases in the MODIS aerosol products are spatially and temporally correlated, thus violating a common assumption of many data assimilation schemes (Dee and Da Silva, 1998;Daley and Barker, 2001).Because the RMSE includes both the variance and the bias, systematic biases in the product greatly increase the magnitude of RMSE and hence reduce product impact in the model.Further, systematic bias can result in "hot" and "cold" spots in model analyses which in turn propagate in model forecasts.To the fullest extent possible, data must be debiased before assimilation.Regions which cannot be adequately debiased must be masked.Two effects that warrant a diagnostic debiasing involve the lower boundary condition and aerosol microphysical bias.Because the lower boundary condition and aerosol microphysical bias often co-vary, we must be careful to examine them independently.In this section we isolate lower boundary condition biases which are subsequently removed from the signal for evaluation of residual microphysical biases.

Lower boundary condition
Figure 3 illustrates how shortcomings in the parameterization of the lower boundary condition result in systematic biases in the MODIS aerosol products.As discussed in Sect.2.1, the MODIS AOD retrieval uses a simple model to describe the relationship between surface reflectance in the near-infrared (2.12 µm) and visible (0.66 µm) wavelengths of MODIS.This relationship is affected by numerous environmental factors, not all of which are captured in the parameterization.Viewing conditions such as the scattering angle are important, but reflectance is primarily determined by the properties of the surface.Further, Fig. 4b suggests that with the exception of the "hotspot" the treatment of surface anisotropy is sufficient.
Using MODIS albedo data, an empirical correction for the biases shown in Fig. 3 was estimated.Numerous correction schemes were examined, and the best results were obtained with the following form: where A 0.66 µm and A 2.12 µm are the MODIS black-sky albedo values at those wavelengths, and m 1 , m 2 and b are parameters to be fit empirically.
To estimate these parameters and test the correction, matched AOD retrievals where τ A < 0.2 were subdivided into estimation and validation subsets (selected by AERONET site name: "A-K" for estimation, "L-Z" for validation).Subsets were made geographically independent to avoid overfitting to confounding regional variation not related to surface properties.Regressing the estimation subset, we obtained a correction of the following form: τ M,corrected = τ M − 2.66A 0.66 µm + 1.25A 2.12 µm + 0.056 (7) The regression coefficient for this relationship (r 2 = 0.24) is of the same order as the r 2 for τ M vs. τ A for the low-AOD estimation subset.Regression of albedo data at other wavelengths (as well as using the 0.66 µm/2.1 µm albedo ratio) was also attempted, but the correlation was not improved (results not shown).Regression results using only data from MODIS-Terra or MODIS-Aqua were very similar (results not shown).The contribution of the surface reflectance to observed radiance and thus to AOD error decreases with increasing AOD.Therefore, we applied the calculated albedo correction only to retrievals with τ M < 0.6.
Evaluation of this correction with the geographically independent validation subset showed robust reduction of differences between τ M and τ A .For 124 AERONET sites in the validation subset where the mean correction to τ M was at least 0.005, the albedo correction reduced the mean AOD error in 85 sites.For those 124 sites, the mean site AOD error and the estimated AOD correction were correlated with r = −0.65,indicating the importance of the surface properties in the error budget at low AOD.
Figure 5 illustrates the nature and effects of the correction shown in Eq. (7).As an example, Fig. 5b shows the magnitude of the correction estimated from one global data day of MCD43 albedo data from May 2008.Well known high biases in desert regions are visible, especially in such regions as the desert western United States, the Andes Mountains, arid Australia and the desert belt across Africa, and the Taklamakan and Gobi deserts.Low biases are visible over the more forested regions over the globe.To demonstrate impact over all AERONET sites, Fig. 5a shows the overall effect of the correction on matched MODIS-AERONET site statistics from the complete 2005-2008 data set.Contour lines in Fig. 5a show the magnitude of the correction as a function of the 0.66 µm and 2.12 µm albedo values.Symbols indicate specific AERONET stations, color coded to indicate the mean AOD bias for τ A < 0.2 (opposite in sign to the correction.That is, we want cold colors to counteract warm).The dashed line indicates a 2:1 ratio of near-infrared to visible albedo, a standard assumption in older over-land AOD retrievals (Kaufman et al., 1997).
Of 224 sites where the mean absolute AOD correction was greater than 0.005, 70% had reduced absolute bias, 77% had improved compliance, and 79% had reduced RMS error.For most albedo regimes, this correction effectively debiases the albedo error term.Figure 6 presents the same series of plots as Fig. 3 after correction.These graphs show a much weaker relationship between AOD errors and surface albedo properties.
As applied, the correction increases the range of surface conditions over which errors are small; however, problems The light green area shows regions that fall outside the "strong" limits (50% or more of cases) but within the "weak" limits of Fig. 11.The gray area shows area where valid albedo data were available, but albedo fell outside the "weak" constraints of Fig. 11 (50% or more of cases).Unshaded areas had no valid albedo for the dates shown.
persist at high albedo.Exclusion of certain retrievals based on surface properties is necessary even after correction.The graphs in Figs. 3 and 6 show two cutoffs for albedo values, shown as vertical lines.Based on the error regimes, it is clear that the retrieval works best in the darkest backgrounds to the left of the dashed line (0.06 µm, 0.11 µm, 0.25 µm, and 0.50 for 0.47, 0.66, 2.1 and A 0.66 µm /A 2.12 µm , respectively).This "Strict" threshold (dashed lines, to the left) reflects the surface conditions corresponding to the best retrieval performance in the uncorrected product; the more lax "Weak" threshold (dotted lines, to the right) provides significantly higher coverage while still maintaining acceptable error statistics in most regions.The geographic ramifications of these thresholds are shown in Fig. 5c, for the same data as shown in Fig. 5b: blue areas have at least 50% of albedos within the "Strong" thresholds, green areas have at least 50% within the "Weak" thresholds, and gray areas indicate places where albedo data are available, but fewer than 50% fall within the thresholds of Figs. 4 and 6.Note that gray areas do not indicate no usable retrievals; areas in green and gray lose more than 50% of data volume when albedo thresholds are applied.
The albedo filtering and correction described here, as well as the snow filtering described in Sect.3, depend on datasets not available in a timely fashion for operational use.Appendix C to this paper discusses the creation and evaluation of an alternative approach to filtering and correction using an 8-year data record of MODIS snow and albedo data.

Residual microphysical bias
To support bias attribution, we examined each site's AERONET derived fine/coarse partition as well as regional surface albedo properties.The noise floor statistics, included for each site in the Supplement, are often indicative of issues with the lower boundary condition.Figure 7 shows the r 2 versus slope for the regression of all τ M (with albedo correction applied) to τ A pairs for the range 0.2 < τ M < 1.4.Terra and Aqua MODIS are in red and blue, respectively.Slopes are largely determined by the higher AOD values, which should be more sensitive to microphysical bias.Within some regions, the variability in estimated slope can be large, even for sites with strong correlations to AERONET observations.This suggests spatially and temporally correlated bias in AOD at scales finer than our regional analysis can resolve.However, for many regions, the slope statistics for individual sites are clustered around values significantly different from 1.0.It is these broad slope biases that we hope to quantify and, to the extent possible, correct.
Once we make our albedo correction as described in Sect.5.1, we must remove the residual slope bias.However, as shown in Fig. 7, there are significant site by site Error in τ M as a function of albedo, after application of albedo correction.These figures can be directly compared with Fig. 3a-d.Vertical lines indicate thresholds of albedo used for filtering of AOD; the dashed lines represent the "strong" constraint, and the dotted lines the "weak" constraint.Note that the "strong" constraint imposes both a lower and upper limit on the 0.65 µm/2.1 µm albedo ratio.
differences in MODIS efficacy.These slope biases are a result of microphysical bias in the retrieval and are related to aerosol properties such as particle size distribution and single scattering albedo.We compared mean MODIS bias to the AERONET derived fine mode fraction (O'Neill et al., 2003) for Terra and Aqua separately (Fig. 8a and b, respectively).
Figure 8 shows a clear microphysical bias in the retrieval, with MODIS underestimating AOD for dusty environments, and overestimating in environments where fine-mode pollution is dominant.On average, these biases are on the order of ±10-20%.Two regions are highlighted in Fig. 8: Sahelian Africa (highlighted in brown) and South America (in green).These regions represent clearly disparate regimes of aerosol physical and optical properties, and the MODIS algorithm includes assumed optical properties for these regions that reflect these differences, but may not be entirely accurate (Levy et al., 2007a).In these regions, deviations from unity are quite large, and may indicate an additional misparameterization in particle absorption.Note that some other areas where fine-mode aerosols dominate may not have the same positive bias (e.g., southern Africa).Unlike over ocean, there is no reliable size information available from MODIS to help correct data based on fine/coarse partition on a retrieval by retrieval basis.Given the divergence of site slope performance even for adjacent sites, there is no basis to generate high resolution spatial correction maps based on the AERONET data.We must therefore attempt large scale regional corrections.We expect that RMSE will be elevated even after correction, because of the currently irresolvable residual bias associated with variation within regions.For the most part, regional corrections are on the order of ±10%, with the largest corrections in the North American boreal region and in Sahelian Africa.Given the dramatic overestimation of AOD in the South American burning season, an additional 35% bias correction is applied to all retrievals with τ M > 1.4 in that region (see Fig. 7p).
Because of the empirical nature of the albedo correction, the slope bias correction is only applied to τ M > 0.2.

Effects of filtering and correction
Table 4 shows the regional impact of successive stages of filtering and correction on two statistical measures, the overall compliance of the Level 2 retrievals in the matched data set, and the regression coefficient r 2 indicating the correlation between MODIS and AERONET AOD.Additional statistics are given in the Supplement.The application of basic QA filtering brings the compliance fraction globally from 62% to 67% (69%) for MODIS-Terra (Aqua), and strongly improves the correlation between τ M and τ A .The additional filtering and correction steps described in this section provide an even larger improvement to the global compliance and correlation.
In all regions, compliance is improved after filtering and correction.The albedo correction, acting on the low AOD values that occur most frequently in every region, yields a global increase of 8% in compliance, with improvements in some regions as strong as 19% (South America).
The albedo correction improves the correlation between MODIS and AERONET in all regions except India and Australia, where it is slightly reduced.The slope correction is linear in nature, but can affect the correlation within regions because it is applied only where τ M > 0.2.On the global scale, the correlation change with application of the slope correction indicates the effects of reducing systematic biases between regions, resulting in an increase of 0.1 in the global r 2 over basic QA alone.In some regions, the albedo correction and slope corrections act in opposite directions.In Sahelian Africa, application of the albedo correction actually reduces compliance by 3-4%, but subsequently the slope correction brings a 9-12% improvement in compliance, resulting in a net benefit of 5-9% over basic QA alone.In eastern CONUS, the application of the albedo correction improves the r 2 values by 0.04-0.05,but reveals a slope bias in the data (estimated slope = 1.04 before albedo correction, 1.09 after), which actually degrades performance for τ M in the range 0.2-1.4.When this slope bias is then corrected, the resulting AOD values have improved compliance for all ranges of τ M over the basic QA alone.
The empirical corrections developed here to compensate for errors in the retrieval lower boundary condition and aerosol microphysical properties were evaluated using different data from that used to calculate the corrections (as described above).However, since both training and validation data was from the AERONET network, it is possible that retrieval errors that occur predominantly in areas far from AERONET Sunphotometers may go uncorrected or undiagnosed.However, the use of spatially explicit albedo data as well as regional covariance of aerosol properties, offers at least some physical basis for global extrapolation of these corrections.
Other recent publications have achieved refinement of MODIS AOD by the use of neural network approaches, where the best mathematical transformation between inputs and outputs is derived mathematically without an explicit theoretical underpinning.These approaches have achieved good validation statistics using MODIS data together with meteorology as inputs to calculate surface PM 2.5 over the southeast United States (Gupta and Christopher, 2009), and have been used to approve the relationship between MODISestimated AOD and AERONET (Lary et al., 2009).
A black-box approach such as a neural network can be useful when validation data are very dense, and these approaches are also useful for demonstrating the limitations to retrieval accuracy imposed by the basic features of the satellite data.The accuracy achieved with a neural-net approach using MODIS radiances (Radosavljevic et al., 2010) tells us the information content of the satellite radiances relative to the AERONET data, that is, an accuracy limited only by representativeness error between the satellite and Sun photometer observations, and uncertainty in the satellite radiances.
However, AOD retrievals based on neural-net approaches have an important limitation for global aerosol applications.As Radosavljevic et al. (2010) acknowledge, the performance obtained with their approach is strongly tied to the exact range of conditions represented in the training sample.The AERONET network is sparse, and while non-absorbing pollution aerosol conditions are relative well-captured by the network, sampling of absorbing aerosols and especially mixed aerosols is limited, and sampling of high-AOD cases is very limited for all aerosol types.A physically-based retrieval method, while it may have larger errors in some cases due to uncertainty in assumed aerosol optical properties, is a more robust method for global applications, because the uncertainties derived from comparison to limited validation data have a physical basis for extrapolation to regions remote from the validation data set.Also, a physically based approach to AOD correction, based on the physics behind the retrieval process, can be used to improve future generations of aerosol retrievals.For example, the wind speed bias correction to the over ocean aerosol retrievals in the MODIS collection 5 product, identified by Zhang and Reid (Zhang and Reid, 2006), is the basis for planned improvements to the MODIS collection 6 aerosol product (R.Levy and L. Remer, personal communication, 2010).
6 Analysis of Level 3 product

Scenarios used in Level 3 comparison: filtering and correction
For the remaining analyses, results are shown for the following filters and corrections.Details of filters and corrections are in Sect. 5.
-"STRONG": "NEW" + "Strong" Matched Albedo Filter The RAW scenario represents all MODIS AOD without any filters or corrections.The BASE is uncorrected MODIS AOD, filtered with the auxiliary datasets included with the MODIS product as well as the contextual filters described in Sect.2.5.1.The BASE scenario represents the best results obtainable with the MODIS product alone.The NEW and CLIM scenarios represent the results of this study, incorporating filters and corrections based on the ancillary data as well as empirical corrections based on the comparison with AERONET data.The STRONG scenario is included to show the effects of strict filtering of surface properties.
Figures 9 and 10 give an overview of the BASE and NEW scenarios, based on data from December 2007-November 2008.Figure 9 shows the seasonal mean AOD for the BASE scenario, and the relative difference of the NEW scenario.Figure 10 shows the data density of the BASE scenario, and the impact on coverage of the additional filtering in the NEW scenario.Effects of filtering on data volume and data quality, and construction of a prognostic error model for aggregated MODIS AOD, are discussed below.

Effects of specific filters and corrections on AOD error and data volume
Any exclusion of data from the assimilation system must consider the balance of data volume and coverage and data errors.This section addresses the effects of specific filtering steps and presents the evidence for the efficacy of these steps.This analysis uses the error statistics from comparison between gridded MODIS data and AERONET data.The compliance statistics before and after filtering are used to calculate the "effective compliance" of the excluded retrievals, which combines the effect of filtering on errors and data volume in a single metric.Table 5 shows these statistics for a cumulative series of filters and corrections.Note that for generation of the gridded MODIS AOD product, the textural filters are applied last, after all other filtering steps.

Textural filters
The textural filters impose a cost in gridded data volume between 10% (0.2 < τ M < 0.6) and 6% (τ M < 0.2).After filtering, compliance in all AOD ranges is improved.Textural filtering has little effect on the distribution of positive and negative errors and is nearly neutral for τ M < 0.2, but removes substantial positive biases at higher AOD ranges.

Quality assurance filtering
Quality assurance filtering includes both the exclusion of MODIS retrievals with mandatory QA values other than "Very Good" as well as the exclusion of retrievals with scattering angles above 170 • (see Sect. 3.1).These filters, applied after textural filtering, result in exclusion of gridded data ranging from 28% of τ M > 1.4 to 50% of τ M 0.2-0.6.For τ M < 0.2, the excluded data have overall compliance only slightly worse than the unfiltered dataset, but exhibit a strong surplus of positive errors.For higher AOD ranges, compliance is poor and positive errors dominate in the excluded data.
Table 4. Effect of filtering and correction steps on overall compliance and correlation in the matched MODIS-AERONET data set.Slope correction in parentheses for South America was applied only to data with τ M > 1.4 (see Fig. 7p).

Exclusion of partially cloudy retrievals
Because this filter is applied at the retrieval level, it will both reduce the data volume of the gridded dataset, and also modify the AOD values in grids with both cloudy and cloud-free retrievals.Exclusion of partially cloudy retrievals imposes an overall data cost of 7-9% of gridded data in the matched dataset, but the spatial distribution of these data are highly non-uniform.Figure 11 shows the impact of cloud exclusion from the BASE data set in terms of change in seasonal mean gridded AOD for Terra and Aqua.Retrievals excluded by application of this filter have variable compliance, and a large surplus of positive AOD errors.Application of this filter does not change the overall compliant fraction, but shifts the balance of errors from positive to negative.

Albedo filters
QA filtering, albedo filtering has the most dramatic effects on data volume among the suite of filters tested in this study.The "Weak" albedo filter (matched) reduces the data volume by 7% for AOD < 0.2, and as much as 13% for AOD > 1.4."Weak" albedo filtering results in large improvements to compliance of the gridded product.Excluded retrievals have relatively high compliance, but errors are overwhelmingly biased positive.Application of the "weak" albedo filter improves compliance at all AOD levels.The "Strong" albedo filters are extremely aggressive, and result in a reduction of data volume by 27% for AOD < 0.2 and by up to 53% for AOD > 1.4 (compared with the "weak" filters).This filter also has a definite spatial pattern of exclusion, as illustrated by Fig. 5c, which leaves large areas with little or no available AOD data.This extremely selective filtering delivers some additional improvement in compliance, but the quality of the excluded data is generally high.For some highly demanding applications, this strict filtering may be appropriate.

Effect of albedo and regional microphysical corrections
Corrections applied to the data have no effect on data volume, although they can change the distribution of AOD values.The albedo correction increases compliance for τ M < 0.2 from 78% (80%) to 83% (84%) for MODIS-Terra (Aqua).Compliance at higher optical depths is slightly improved as well.
Regional slope corrections improve AOD for all but the highest AOD values (τ M > 1.4).Correlation between gridded MODIS and AERONET AOD is also improved by the regional slope correction, from r 2 = 0.61 (0.60) to r 2 = 0.65 (0.65) for MODIS-Terra (Aqua).

Prognostic error model for gridded product
Uncertainty in the gridded MODIS AOD data is estimated using Eq.(5) (Sect.2.4).For each region, a linear estimate of RMSE as a function of τ M and a "noise floor" RMSE used as a minimum value were calculated.Table 6 shows the parameters of the uncertainty estimation calculation for each region with sufficient data volume, for the RAW, BASE, and NEW scenarios.Note that while Fig. 3 shows a bi-linear estimation to account for the different error characteristics of very high AOD values, the uncertainty estimates for the gridded product use only a single linear estimate, because data volumes are insufficient for a robust calculation at high AOD.Regional estimates of the "noise floor" RMSE are made for all cases having at least 100 points with τ M < 0.2 in the matched aggregated data set.Estimates of the linear relationship are made only for cases with at least 100 points with τ M > 0.2.
The parameters in Table 6 can be used to estimate the uncertainty (ε) of any gridded MODIS AOD. Figure 12 shows maps, based on the same data as Figs.9-10, representing the mean estimated uncertainty for the BASE scenario, and the fractional difference in uncertainty for the NEW scenario.The filtering and correction described in this study reduces the estimated uncertainty in the MODIS AOD over nearly the entire globe.For τ M = 0.2, estimated uncertainty is lower in all regions except for Australia (ε BASE = 0.064, ε NEW = 0.072).For τ M = 1.0, ε is lower in all regions except South America, Southern Europe/Mediterranean, Eurasian Boreal, and Peninsular SE Asia.In each of those cases, the difference between ε BASE and ε NEW is less than 0.05.Globally, one standard deviation of "Very Good" data falls within the (0.05 + 0.2 × τ A ) error thresholds, and we concur with the MODIS science team recommendation that only "Very Good" data be used (Table 1a).However, for most applications a prognostic RMS error model with a noise floor is more appropriate.For global applications using MODIS Level 2 data over land, we recommend the use of the greater of 0.08 or 0.02 + 0.22 × τ M for Terra and 0.07 or 0.01 + 0.26 × τ M for Aqua (Table 1c).

Global noise issues:
The amount of scatter in the MODIS-AERONET comparison was sensitive to observing conditions, especially viewing geometry.Owing to an increase in optical path length and pixel size at larger scan angles, MODIS products have higher compliance at higher scan angles (Fig. 3).Conversely, errors are largest at nadir owing to shorter optical path length and perhaps increases in BRDF gap probabilities.One consequence of this is that comparisons between MODIS and MISR are not necessarily indicative of behavior of the MODIS product at all scan angles.Scattering angle does not appear to affect retrieval error except for scattering angles >170 • owing to the hotspot effect (Fig. 3).

Global cloud bias:
The MODIS retrieval is highly clear sky conservative with very few retrievals with "Very Good" QA and MODIS-detected cloud fraction >15%.Even so, retrievals with non-zero MODISdetected cloud fractions still have perceptible high biases (Fig. 3).Elimination of partially cloudy retrievals reduces mean seasonal AOD by as much as 0.05 over large regions of the tropics (Fig. 11).

Global snow bias:
While Collection 5 has a much improved snow filter, we still find periodic positive perturbations at high latitudes.By using a spatially and temporally extended snow filter based on the snow flag in the MODIS MCD43 albedo product, we can reduce the incidence of positive biased AOD in northern latitudes (Table 2). 5. Regional variability: Despite good overall global statistics, MODIS products have widely varying regional efficacy.
Widely diverging statistics for adjacent AERONET sites suggest spatially and temporally correlated bias.Regions such as the Eastern CONUS and Europe perform best.Sahelian Africa shows the poorest performance.Sites in East Asia have highly mixed efficacy.Urban located sites also tend to have poor efficacy.Regions which experience intermittent smoke events show significant slope bias on a per event basis, which is manifested in highly variable (but uniformly positive) slope bias for sites in South America and boreal North America.This regional variation in retrieved AOD has been seen in other AOD products; comparison of these products has been used to identify regions and seasons where more detailed observations, in situ or otherwise, might be necessary to understand the behavior of the satellite retrievals.
6. Albedo correction: Some bias at low AOD can be corrected based on the use of surface albedo maps.An empirical relationship between surface albedo at 0.66 µm and 2.1 µm and MODIS-AERONET differences in AOD explained more than 20% of the variance in those differences (Eq.7).This correction was shown to be robust using geographically independent data (Sect.5.1), and improves compliance of MODIS AOD by around 8% (Table 4).7. Regional slope correction: A clear, if noisy, relationship between AERONET estimated fine/coarse partitioning and MODIS slope bias was identified (Fig. 8), and used as the basis of a regional correction.This correction Screening and correction for snow, albedo, and regional slope bias, as described in this study, reduces data volume by a further ∼10% while improving the compliant fraction to 80% for both Terra and Aqua (Table 5).
9. Bias between Terra and Aqua: Lastly, using these products, small but perceivable time dependent differences between Terra and Aqua (globally on the order of 0.02) are visible in statistical analysis of the Level 2 products, but are within the uncertainty of the products.These differences can be seen more clearly using the Level 3 products with extensive averaging (see Appendix D).
The regional uniformity of this difference suggests a shortcoming in radiometric calibration of the two instruments, which must be resolved or at least corrected to perform a long-baseline study of global aerosol (Zhang and Reid, 2010).
Taken together, the MODIS Collection 5 aerosol product, with the QA procedures laid out in this study, can be used to produce a product with desirable qualities for data assimilation: low systematic bias and random error and relatively well-characterized residual uncertainties.The next step is testing of this product in the NAVDAS-AOD data assimilation system, which is currently underway.The MODIS aerosol science team is also currently producing an updated version of the MODIS AOD retrieval (Collection 6) (L.Remer, personal communication, 2010), which will include significant changes that hopefully will address some of the biases identified in this study.The lessons from this data evaluation exercise can also be applied to data from other sensors both current and future.

Sampling considerations for matching MODIS 10 km AOD retrievals to AERONET
Representativeness error is an unavoidable consideration for comparison of datasets with different sampling properties.In the case of AERONET and MODIS AOD, the former is a direct-sun measurement representing atmospheric conditions along a line between the sensor and the sun, while the latter is a product made using data covering an area 10 km square on the ground (although much larger away from nadir), but with limited sampling within that footprint due to pixels rejected in the retrieval process.Nominally, one could specify a distance between the nominal center of the MODIS AOD retrieval and the AERONET station that would ensure that the AERONET station fell within the ground footprint of the MODIS retrieval, but this would not ensure a spatial match, and depending on the sun angle, the atmospheric column sampled could still be completely different for the two sensors.
Rather than attempt to capture all of these factors in building a comparison dataset, they are instead rolled into random error which indicates the precision with which individual retrievals can be analyzed in the comparison dataset.Of greater concern are the systematic biases in a comparison dataset caused by the interaction between the scales of the observations and the scales of the phenomena observed.Over land, these biases can potentially be large as fine plumes near sources will often have characteristic dimensions smaller than 10 km.The overall impact of these biases is examined in this section.
Figure A1 shows 4 compliance plots, indicating the change in mean bias and spread between MODIS and AERONET as a function of spatial and temporal separation of the observations.Effects of spatial separation are evaluated for 3 different ranges of τ A .In all cases, a slight decrease in τ M is seen with distance from the AERONET station.This reflects the location of AERONET stations; while systematic attempts are made to place these stations in locations at some remove from local sources, the requirements of accessibility compromise those attempts to a small degree.At moderate AOD values (Fig. A1b), the bias crosses over from positive to negative at a distance of approximately 15 km; at high τ A (Fig. A1c)), the bias is negative at all separations.AERONET retrievals with high AOD (Fig. A1c) are often indicative of near-source features; the strong gradients associated with these features are expected to result in higher representativeness error and negative bias of the large-footprint MODIS AOD compared to the AERONET AOD.
For this study, a maximum distance threshold of 30 km between the nominal center of the MODIS retrieval and the AERONET station is used.This threshold balances increasing random variability with distance, and the need for as many comparison data points as possible.
Figure A1d shows bias and compliance as a function of temporal separation, calculated using only matched retrievals within 10 km of AERONET stations.For a time period of one hour before or after the AERONET observation, only very slight changes in bias and compliance can be observed.For this study, MODIS and AERONET retrievals were paired with a maximum delay of ±30 min.
"high-AOD" case was 0.006 lower than for the "low-AOD" case.For 89% of cases, the absolute difference in the calculated corrections was less than 0.02.The more positive correction corresponded to the low-AOD case in 67% of cases.These results suggest that elevated AOD can result in a small bias to retrieved albedo, but this effect is generally within the range of uncertainty of the MODIS albedo.The interaction with the albedo correction to MODIS AOD calculated in this study is slight, and very rarely larger than the noise level of the MODIS AOD.

Use of climatological data for snow and albedo filtering and albedo correction
The primary purpose of our work with the MODIS AOD product is the improvement of the NAVDAS-AOD operational aerosol data assimilation system.The demands of this system emphasize removal of outliers and quantitative characterization of observation errors.The error correction demonstrated in Sect. 5 expands the range of surface conditions over which MODIS can retrieve albedo within acceptable error limits (compare Figs. 3 and 6).The statistics in Sect.4.2 indicate that with the albedo correction in place, we can use the less stringent ("weak") albedo thresholds.However, this correction is dependent on the MODIS MCD43 albedo product, which is not available in a sufficiently timely fashion for operational use.Thus, the initial operational version of the MODIS AOD data for data assimilation will need an alternative means of surface albedo correction.In addition, the snow detection algorithm described in Sect. 3 above also relies on the MODIS albedo product, and so will also require a climatological substitute method for operational implementation.
In order to develop a filtering and correction method usable for real-time applications, we used the MODIS MCD43C3 data record from 2000-2007.Because of the 8-day production schedule for the 16-day product, this period yields as many as 16 observations for each 0.05 • grid cell.This albedo climatology included data points matching 98% of the matched MODIS-AERONET dataset for 2008.Note that the occurrence of positive non-compliant errors in the retrievals with no climatological albedo data is very high (89%).This is because nearly all of the locations without any successful albedo retrievals during 2000-2007 are in permanently snow-and ice-covered regions (see Appendix B).
Our potential snow climatology was constructed by flagging all locations and dates where 2000-2007 data indicated snow cover at any time.The matched dataset indicated snow contamination in less than 0.1% of 2007 MODIS AOD retrievals (Table 2).The climatological snow filter captures more than 99% of these retrievals, while reducing overall 2007 data volume by 4%.The fraction of retrievals flagged by the climatological filters with real snow contamination is likely to be low, but the excluded data points still have a positive bias, and a poor correlation with AERONET, relative to the entire dataset (see Table C1, and compare Table 2).Note that some of the error statistics for retrievals flagged by the climatological snow filter are better than the overall dataset.This is because many areas susceptible to snow contamination are in darker, denser vegetation, where the AOD retrieval is generally more accurate (in the absence of snow or other disrupting factors).The climatological snow filter captures fewer retrievals than the spatially and temporally extended filter.This could indicate a more effective filter, but further evaluation with the data assimilation system is needed to determine how best to eliminate snow contamination without excessive loss of data.We used the albedo data to calculate a correction based on all valid albedo measurements for each date from 2000-2007.Because missing data are less common with the climatological approach, the multi-year dataset can be used to implement an additional consistency check without excessive data loss.A filter was created to exclude retrievals where fewer than three data points were used to calculate the mean albedo correction, as well as retrievals where the range of corrections exceeded 0.1.These two exclusions together accounted for 7% of data volume in the matched dataset.With this filter applied, the climatological correction was within 0.03 of the matched correction for 98% of cases.
With the quality filter, the climatology also delivers a slightly better correction.When the estimated correction is larger than 0.05, the correction moves τ M closer to τ A in 73% of cases, compared with 69% of cases for the matched albedo correction.
Any locations/dates where the climatology indicated albedo above the threshold in 20% of more of cases were excluded.The combination of the consistency checks applied for application of the albedo correction and the "weak" albedo thresholds reduced data volume by 17%, compared with 16% using the matched data.Data volume reduction using the "strong" albedo thresholds was much greater.
The Level 3 product generated with this climatological approach (CLIM scenario) has similar characteristics to the NEW scenario made using matched snow/albedo data (see   Sect. 7 for descriptions of Level 3 product scenarios).Table C2 shows the prognostic error model derived for this product, and Fig. C1 shows the seasonal mean estimated uncertainty relative to the BASE scenario (compare to Fig. 12).Remaining uncertainty for the CLIM scenario is generally higher than the NEW scenario, but degradation from the BASE scenario, as seen in Australia in Fig. 12, is absent from the CLIM scenario.

Appendix D Bias between terra and aqua MODIS retrieved AOD
The existence of two MODIS sensors, sampling late morning and early afternoon conditions, has led to the application of MODIS data to study diurnal variation in important Earth system processes.Important processes with significant diurnal variability include fires (e.g., Giglio, 2007;Roberts et al., 2009), cloud cover (e.g., Jin et al., 2009), and cloud properties (e.g., Meskhidze et al., 2009).The interaction of variation between aerosol sources and cloud properties is an object of intense scientific interest, as it may shed light on the interactions between aerosol particles and clouds that can strongly affect radiation budgets in polluted regions.
Comparison of Terra and Aqua MODIS to examine diurnal cycles of Earth system processes is dependent on detection and correction of any artifacts of calibration between the two sensors.Because the two sensors never overlap except near the poles, and even there only with large differences in viewing geometry, indirect methods must be used.There is a large literature on these methods, and studies generally conclude that Terra and Aqua MODIS calibration in visible and infrared bands is basically within the precision of the indirect methods used, generally cited as on the order of 2% in topof-atmosphere reflectance (e.g., Wu et al., 2004).However, the AOD retrieval is extremely sensitive to small differences in reflectance under some conditions.
The analyses performed in this study, for virtually all cases where Terra and Aqua results differ, have shown results consistent with τ M from MODIS-Aqua being slightly higher than τ M from MODIS-Terra.For instance, in every region, the fraction of low-τ A retrievals with negative errors below the compliance threshold is higher or equal for MODIS-Terra in every region except Mid-Latitude East Asia (Table 3).The same is true for every QA level (Table 1).In general, these differences could not pass tests of statistical significance, and they are often beyond the least significant digit of our results.However, the accumulation of multiple statistical results gives qualitative support to the idea that these two www.atmos-meas-tech.net/4/379/2011/Atmos.Meas.Tech., 4, 379-408, 2011

Fig. 1 .
Fig. 1.(a) Geographic extent of MODIS AOD retrieval coverage over land (shaded).Colored contours indicate the mean annual AOD retrieved from MODIS at each location, averaged using data from 2001-2008.AERONET stations used in this study are marked with "x".Solid black lines divide the land surface into regions, which are discussed in Sect. 5. (b) MODIS 0.55 µm AOD versus AERONET 0.55 µm AOD (AERONET AOD is interpolated, see Sect.2.2).Symbols are as follows: CONTOURS: each shaded region represents 10% of the matched MODIS-AERONET data points, organized by data density in bins of 0.05 optical depth.The contoured area includes 90% of the matched retrieval, the remaining 10% are shown as POINTS.SOLID LINE: Each vertex represents 50 000 paired retrievals, sorted by τ A .DOTTED LINES: indicate the 25th and 75th percentiles of each bin along the solid line.DASHED LINES: indicate the compliance metric of Eq. (1).The 1:1 line is also shown.

Fig. 3 .
Fig. 3.Error in τ M as a function of surface albedo at (a)0.47 µm, (b) 0.66 µm, (c) 2.12 µm.(d) Error in τ M as a function of the ratio of MODIS albedo at 0.65 µm and 2.1 µm.The lines were constructed by sorting the albedo of matched retrievals into bins of 20 000, each of which is shown as a vertex of the solid (mean bias) and dotted (25th and 75th percentiles of bias) lines.The gray bars at the top and bottom of the graph represent the fraction above and below, respectively, the compliance tolerances of Eq. (1).All statistics were calculated using only matched data with MODIS QA of "Very Good".

Fig. 4 .
Fig. 4. Error in τ M as a function of (a) MODIS view angle and (b) scattering angle.Solid line indicates mean bias; each vertex is the average of 20,000 retrievals.Dotted lines indicate 25th and 75th percentiles of each bin.Gray bars at the top and bottom of the graph represent the fraction above and below, respectively, the compliance tolerances.All statistics were calculated using only matched data with MODIS QA of "Very Good".(c) Error in τ M as a function of MODIS retrieved cloud fraction.The straight dashed lines indicate the mean AOD bias, as well as the non-compliant fraction (low and high) for all retrievals with no indicated cloud (Cloud Fraction = 0).

Fig. 5 .
Fig. 5. Albedo correction and filtering and effects on MODIS AOD.(a) Surface albedo and AOD bias for individual sites.Symbols indicate the mean surface albedo for each AERONET site, with colors indicating the mean bias in AOD for τ A < 0.2 at that site.Contours behind indicate the estimated AOD correction based on the surface albedo (Eq.7).The dashed line indicates the relationship Albedo (0.65 µm) = Albedo(2.1µm)/2.The symbols are placed to indicate the mean albedo for each AERONET site for the matched MODIS-AERONET dataset.The colors of the symbols indicate the mean bias of τ M at each site, on a scale that is the reverse of the contours.Sites marked with "+", are labeled with the name of the site and the mean bias of τ M for the matched data at that site.(b) Example of the estimated albedo correction calculated using Eq.(7), using the 16-day MODIS albedo product for days 177-193 of 2008.(c) Effect on geographic coverage of albedo filtering of MODIS AOD product.The map above illustrates 3 zones based on the MODIS albedo product from days 177-193 of 2008.The dark blue area highlights regions where the surface albedo falls within the "strong" limits shown in Fig. 11.The light green area shows regions that fall outside the "strong" limits (50% or more of cases) but within the "weak" limits of Fig.11.The gray area shows area where valid albedo data were available, but albedo fell outside the "weak" constraints of Fig.11(50% or more of cases).Unshaded areas had no valid albedo for the dates shown.
Fig.6.Error in τ M as a function of albedo, after application of albedo correction.These figures can be directly compared with Fig.3a-d.Vertical lines indicate thresholds of albedo used for filtering of AOD; the dashed lines represent the "strong" constraint, and the dotted lines the "weak" constraint.Note that the "strong" constraint imposes both a lower and upper limit on the 0.65 µm/2.1 µm albedo ratio.
Fig. 7. (A-O)Slope bias and correlation of τ M vs τ A , after albedo correction, using matched retrievals with 0.2 < τ M < 1.4.For each region, statistics for sites within that region are plotted with the number of matched data points.Red (Blue) indicates MODIS-Terra(Aqua).Vertical lines indicate statistics for each satellite calculated for the region as a whole.(P) For South America, slope and correlation statistics are shown calculated using only τ M > 1.4.

Fig. 8 .
Fig. 8. Slope of τ M /τ A (Eq. 2) as a function of mean AERONET fine-mode fraction, calculated using only matched retrievals with τ A > 0.2.Solid line indicates trend of slope for all retrievals.Symbols indicate mean values for individual AERONET sites, solid symbols highlight those sites where τ M and τ A are correlated with r 2 > 0.25.Gray areas indicate non-compliant fraction, as in Figs.6-9.Sites in Northern Africa are highlighted in brown; sites in South America are highlighted in green.Statistics were calculated after application of albedo correction (See Sect.4).

Fig. 9 .
Fig. 9. Mean AOD for 1-• grid cells.The upper map shows the absolute mean AOD for the BASE scenario of basic QA filtering and no corrections applied.The lower map indicates the relative effect on seasonal mean AOD from the filtering and correction applied for the NEW scenario.Data used were from MODIS-Terra for the period June-August 2008 (results for other seasons are included in the Supplement).

Fig. 10 .
Fig. 10.Number of days with data for 1-• grid cells.The upper map shows the absolute number of days with data in each season for the BASE scenario of basic QA filtering and no corrections applied.The lower map indicates the relative effect on data frequency from the filtering and correction applied for the NEW scenario.Data used were from both MODIS instruments for the period June-August 2008.

Fig. 11 .
Fig. 11.Effect of excluding retrievals with MODIS-detected cloud cover on mean MODIS AOD.Maps show the difference between the BASE scenario and a scenario which is identical but includes retrievals with MODIS-detected cloud.Upper map shows data from MODIS-Terra; lower map shows MODIS-Aqua.All maps use data from June-August 2008.

Fig. 12 .
Fig. 12. Mean estimated uncertainty in AOD for 1-• grid cells.The upper map shows the mean uncertainty estimated for the BASE scenario using the prognostic error model.The lower map shows the ratio of the mean uncertainty for the NEW scenario to the BASE scenario.Data used were from MODIS-Terra for the period June-August 2008.

Fig. C1 .
Fig. C1.Estimated uncertainty in CLIM scenario relative to BASE scenario.Compare to Fig. 12. Rows indicate different seasons.Data used were from MODIS-Terra for the period December 2007-November 2008.
Table 1b lists the compliant

Table 1a .
Distribution of AOD values and linear regression results for matched MODIS-AERONET dataset, stratified by MODIS QA value.

Table 2 .
Statistical analysis of filters to remove snow contamination in MODIS AOD.Data shown were calculated using only MODIS AOD with MODIS QA of "Very Good".

Table 3 .
Statistical evaluation of MODIS AOD by region, after application of basic QA filtering.
Conversely, Western CONUS is one of MODIS's greater challenges and shows www.atmos-meas-tech.net/4/379/2011/Atmos.Meas.Tech., 4, 379-408, 2011 one of the highest RMSEs.Previous studies have reported a significant bias (+0.2 or greater) over the arid regions due to shortcomings in the lower boundary condition (Drury et al., 2008).This is manifested in our statistics as extremely poor compliance and strong positive bias for above-background values of τ M , which are often retrieved in conditions where the true AOD is a low background value.Regression fits are generally weak because the AOD in the pristine desert has little range.Heavily urbanized sites such as around Southern California or the San Francisco Bay Area perform nearly as poorly.Like the boreal zone, large biomass burning events can dominate individual regressions, with slopes ranging from 0.6 to more than 1.6.Remote sites with less arid landscapes -e.g., HJ Andrews (44 • N, 122 • W), Missoula (47 • N, 114 • W), and Fresno (37 • N, 120 • W) -have more reasonable performance, with slopes generally within ±0.1 of unity and >70% compliance.
Regression slopes for moderate AOD are typically within +/0.1 of unity.Only a few sites (e.g., Forth Crete -35 • N, 25 • W), have slope biases larger than 1.2.Site specific Terra and Aqua slopes also correlate very well, within ±0.1 of each other for 34/44 sites.Scores are also excellent in the Eurasian Boreal, with greater than 75% compliance for 17/26 sites.African aerosol regimes observed by MODIS can be broken down into the northern dust/smoke impacted Sahelian region, biomass burning impacted equatorial Africa, and the burning and pollution region of southern Africa (No retrievals are made in the pure dust regimes of the Sahara).The Sahelian Africa environment strains the MODIS algorithms, with both high background albedo and variable fine/coarse partition particle size.RMSEs are more than 50% for all τ A < 0.8.Slopes for moderate AOD vary between 0.52 and 1.20, although in general they are in the 0.6-0.8range.AOD is negatively biased in all ranges at all sites except for Izana (28 • N, 16 16 τ A Atmos.Meas.Tech., 4, 379-408, 2011 www.atmos-meas-tech.net/4/379/2011/for Aqua).• W).

Table 5 .
Effects of filtering and correction steps on compliance and data volume of gridded MODIS AOD product.For each step, compliance statistics are given before and after.For filtering steps, the before and after compliance is used to calculated an effective compliance of excluded retrievals.

Table 6 .
Prognostic RMS error model based on comparison of L3 gridded MODIS AOD to AERONET.Model has the form ε EST = MAX(A, B + C τ M ) (Eq. 4)."N/A" indicates regions where data volume was insufficient to calculate model parameters.Error for these regions is estimated using the Global parameters.

Table C1 .
Analysis of snow filter based on MODIS MCD43 climatology.Compare to Table 2.