Exploring systematic offsets between aerosol products from the two MODIS sensors

Long-term measurements of global aerosol loading and optical properties are essential for assessing climate-related questions. Using observations of spectral reflectance and radiance, the dark-target (DT) aerosol retrieval algorithm is applied to Moderate Resolution Imaging Spectroradiometer sensors on both Terra (MODIS-T) and Aqua (MODIS-A) satellites, deriving products (known as MOD04 and MYD04, respectively) of global aerosol optical depth (AOD at 0.55 µm) over both land and ocean, and an Ångström exponent (AE derived from 0.55 and 0.86 µm) over ocean. Here, we analyze the overlapping time series (since mid-2002) of the Collection 6 (C6) aerosol products. Global monthly mean AOD from MOD04 (Terra with morning overpass) is consistently higher than MYD04 (Aqua with afternoon overpass) by ∼ 13 % (∼ 0.02 over land and ∼ 0.015 over ocean), and this offset (MOD04 – MYD04) has seasonal as well as long-term variability. Focusing on 2008 and deriving yearly gridded mean AOD and AE, we find that, over ocean, the MOD04 (morning) AOD is higher and the AE is lower. Over land, there is more variability, but only biomass-burning regions tend to have AOD lower for MOD04. Using simulated aerosol fields from the Goddard Earth Observing System (GEOS-5) Earth system model and sampling separately (in time and space) along each MODIS-observed swath during 2008, the magnitudes of morning versus afternoon offsets of AOD and AE are smaller than those in the C6 products. Since the differences are not easily attributed to either aerosol diurnal cycles or sampling issues, we test additional corrections to the input reflectance data. The first, known as C6+, corrects for long-term changes to each sensor's polarization sensitivity and the response versus the scan angle and to cross-calibration from MODIS-T to MODIS-A. A second convolves the detrending and cross-calibration into scaling factors. Each method was applied upstream of the aerosol retrieval using 2008 data. While both methods reduced the overall AOD offset over land from 0.02 to 0.01, neither significantly reduced the AOD offset over ocean. The overall negative AE offset was reduced. A collection (C6.1) of all MODIS Atmosphere products was released, but we expect that the C6.1 aerosol products will maintain similar overall AOD and AE offsets. We conclude that (a) users should not interpret global differences between Terra and Aqua aerosol products as representing a true diurnal signal in the aerosol. (b) Because the MODIS-A product appears to have an overall smaller bias compared to ground-truth data, it may be more suitable for some applications. However (c), since the AOD offset is only ∼ 0.02 and within the noise level for single retrievals, both MODIS products may be adequate for most applications.

1 Introduction

Measurements of aerosol loading and optical properties are essential for many applications, including quantifying global direct aerosol radiative forcing for climate studies (e.g., Belloiun et al., 2005; Chung et al., 2005; Yu et al., 2006; Kahn, 2011; Boucher et al., 2013), investigating the effect of aerosols on cloud microphysical properties and lifetimes (Nakajima et al., 2001; Lohmann and Feichter, 2005; Koren et al., 2008, 2012; Niu and Li, 2012) and estimating global exposure to air pollution (van Donkelaar et al., 2010; Evans et al., 2012; Boys et al., 2014; Ma et al., 2016). Because aerosols vary significantly by region (Kaufman et al., 2002) and have a lifetime on the order of days (Haywood and Boucher, 2000; Croft et al., 2014), near-daily observations over the entire globe are necessary to characterize the global aerosol system. The Global Climate Observing System (GCOS, 2011, 2016) has designated particular aerosol parameters as essential climate variables (ECVs) for quantifying Earth's climate system and change. To be considered a viable Climate Data Record (CDR), an ECV must be measured globally, with specified accuracies, precisions, spatial and temporal resolution. The ECV also must be measured over the long-term (e.g., multiple decades).

Aerosol optical depth (AOD – a measure of column-integrated aerosol loading) is a designated ECV. To meet requirements as a CDR (Popp et al., 2016; GCOS, 2016), AOD must be measured globally, with a spatial resolution of 10 km or finer and accuracy better than 0.03 or 10 %. In addition, this AOD record must be multi-decadal and drift less than 0.01 decade⁻¹. Polar-orbiting passive satellite sensors are able to provide spatial coverage, frequent sampling and long continuity of data that could be the basis for such a record. In particular, the Moderate Resolution Imaging Spectroradiometer (MODIS) instruments on board the polar-orbiting satellites Terra (since 2000) and Aqua (since 2002) provide state-of-the-art spatial resolution and near-daily retrievals of AOD and other aerosol parameters on a global scale. The length of the aerosol records has prompted studies of the trends in global and regional aerosol loading, and subsequently, estimates of changes in aerosol forcing or radiative effects (Zhang and Reid, 2010; Hsu et al., 2012; Chin et al., 2014; Alfaro-Contreras et al., 2017; Colarco et al., 2014).

Using remote sensing to detect changes or trends in the physical world (e.g., the ambient aerosol), however, requires confidence that the algorithms and retrievals are consistently applied and that the sensors themselves (e.g., calibration, sampling and orbital characteristics) are also consistent. Creation of long-term climate data records often requires the observations from different instruments and platforms to be combined, because a single instrument may not provide sufficient spatial, temporal, or long-term coverage (e.g., the Global Precipitation Climatology Project; Adler et al., 2018). As the community moves towards creating aerosol CDRs that span the lifetime of more than one sensor, we need to pay even more attention to systematic biases and offsets.

In this study, we compare the aerosol climatology from the two identically designed MODIS sensors that fly simultaneously for over 15 years. The specifications of the instruments are essentially identical (sensor characteristics, calibration methods), and the retrieval algorithms are identical. In Sect. 2, we show that there are systematic differences in the derived global aerosol products for Collection 6. Although each sensor shows an insignificant global drift, their differences appear as a small but statistically significant trend. More alarming is that the two data sets are offset from each other, on average by 13 % of their global mean. This is larger than the GCOS requirements for accuracy (GCOS, 2016) and will introduce greater uncertainty than needed in narrowing error bars on current estimates of global aerosol radiative forcing (Boucher et al., 2013). In Sect. 3 we sample model output data to show that the differences between Terra and Aqua aerosol climatology are most likely unphysical, and in Sect. 4 we tested two methods of calibration correction to reduce the problem. Section 5 offers a discussion and conclusion, including suggestions for future calibration efforts from a product-based perspective.

2 MODIS and the Collection 6 aerosol time series

2.1 MODIS

Terra and Aqua were launched in late 1999 and early 2002, and are polar-orbiting sun-synchronous satellites. Terra (Aqua) has a 10:30 (13:30) local equator crossing time and descends (ascends) on the sunlit part of the Earth. From each satellite, MODIS observes top-of-atmosphere (TOA) reflectance (solar origin) and radiance (terrestrial origin) in 36 wavelength bands (0.41 < λ < 14.2 µm): 19 are reflective solar bands (RSBs; λ ≤ 3.9 µm) and the remainder are thermal emission bands (TEBs). Nominal (at nadir view) spatial resolution is 0.25 km for two bands (0.65 and 0.86 µm), 0.5 km for five bands (0.47, 0.55, 1.24, 1.63 and 2.11 µm) and 1 km for the remainder. From orbit ∼ 700 km and ±55^∘ scan angle (0^∘ is nadir view), MODIS observes a ground swath of 2300 km, which provides near-global coverage every day and complete coverage every 2 days. It should be noted that the original mission lifetimes for Terra and Aqua were nominally 5 years.

In terms of sensor specifications, including spectral wavelength characteristics, calibration methods and presentation of data and file formats, the two MODIS instruments (MODIS-Terra or MODIS-T and MODIS-Aqua or MODIS-A) are twins. MODIS data products (from raw data through high-level aggregations), MODIS-T and MODIS-A, are known as MOD and MYD, respectively.

As the flagship sensor aboard two high-profile satellites, MODIS is a complex enterprise. Though the scientific literature is immense, most relevant information can be gleaned across the myriad of NASA websites. The general sensor concept and design are presented at https://modis.gsfc.nasa.gov/about/ (last access: 27 June 2018). Sensor characterization and calibration, up to the processing of geophysically relevant reflectance and radiance data (known as level 1B or L1B) are handled by the MODIS Characterization and Support Team (MCST; https://mcst.gsfc.nasa.gov/, last access: 27 June 2018). Retrievals and derivation of geophysical parameters, known as level 2 (L2) products, are described under Land (https://modis-land.gsfc.nasa.gov/, last access: 27 June 2018), Ocean (https://oceancolor.gsfc.nasa.gov/, last access: 27 June 2018) and Atmosphere (https://modis-atmos.gsfc.nasa.gov/, last access: 27 June 2018) disciplines. The aerosol retrieval follows the Atmosphere discipline and is collectively known as the MxD04_L2 product. Level 3 (L3) daily (MxD08_D3) and monthly aggregations (MxD08_M3) of the aerosol (and other Atmosphere products including clouds) are collectively known as MxD08. All data processing is handled by the MODIS Adaptive Processing System (MODAPS; (https://earthdata.nasa.gov/about/sips/sips-modaps, last access: 27 June 2018) and archival of retrieved products is handled at the Level 1 and Atmosphere Archive and Distribution System (LAADS; https://ladsweb.modaps.eosdis.nasa.gov/, last access: 27 June 2018).

As briefly described on the MODIS design page (https://modis.gsfc.nasa.gov/about/design.php, last access: 27 June 2018), MODIS operates via a Scan Mirror Assembly (SCA), which uses a continuously rotating double side mirror (MS-1 and MS-2). The optical system directs the radiation to four assemblies, one for each of the VIS, NIR, SWIR/MWIR and LWIR spectral regions (to cover the visible, near, shortwave/mediumwave and longwave infrared spectra respectively). For the purpose of the aerosol retrieval, we are primarily concerned with the VIS, NIR and SWIR portions, which include the RSBs. To maintain calibration of the RSBs (originally performed in-lab prior to launch), the system includes a view of space, along with on-board calibration via a Solar Diffuser (SD), a blackbody source (BB), a Spectroradiometric Calibration Assembly (SRCA) and a Solar Diffuser Stability Monitor (SDSM). The redundant system (Xiong and Barnes, 2006) allows the MCST to continually update the calibration coefficients (including gain and offsets), along with satellite and viewing geolocation (latitude/longitude, angles, altitude, etc.). This provides the best possible accuracy and uncertainty estimate for the reflectance and radiance data (L1B), which are expected to be accurate to 2 % for reflectance and 5 % for absolute radiance under typical magnitude conditions. Note that the methodology for MODIS calibration, especially in regards to extending the mission from a nominal 5 years to the current 18+ years, has been continuously evolving. This is reflected within the MCST web page and associated literature.

MODIS products are grouped together in collections, in that a consistent protocol is used to derive L1B data, and then consistent algorithms are used to derive the L2 products and L3 aggregations. The same combination of L1B, L2 and L3 production rules are maintained so that all data in a collection are created the same way. This includes production of new data (collected forward in time) known as “forward processing” and archived data known as “reprocessing”. Under a collection, the entire time series of derived parameters (e.g., AOD) from MODIS-T (2000–present; MOD04 product) and MODIS-A (2002–present; MYD04 product) should be consistent with each other and presumably provide consistent global climatologies of the parameters. The most recent complete collection is known as Collection 6 (C6), and it encompasses time series from both MODIS sensors. Collection 6.1, discussed in Sect. 5, began processing in late October 2017.

Aerosols are ubiquitous in the atmosphere, and there are multiple algorithms for removing the aerosol effect (known as atmospheric correction, or AC) when retrieving properties of land (e.g., Vermote and Kotchenova, 2008; Lyapustin et al., 2011) or ocean (e.g., Ahmad et al., 2010) surfaces. While these AC algorithms report the aerosol information, they are not necessarily focused on providing a global (land + ocean) aerosol product. For global aerosol coverage, NASA uses three separate algorithms to create the MxD04_L2 product. Two of these are considered dark target (DT; https://darktarget.gsfc.nasa.gov, last access: 27 June 2018) because they seek conditions in which the surface appears dark in visible wavelengths. These include retrieval over remote ocean (DT-O; Tanre et al., 1997; Remer et al., 2005; Levy et al., 2013) and retrieval optimized over vegetated or dark-soiled land surfaces (DT-L; Kaufman et al., 1997; Levy et al., 2007a, b, 2010, 2013). There is also the deep-blue algorithm (DB; https://deepblue.gsfc.nasa.gov, last access: 27 June 2018) over land (DB-L) which was developed for brighter surfaces (Hsu et al., 2004) and more recently extended to dark surfaces as well (enhanced DB: Hsu et al., 2013). Here, we focus on the two DT algorithms (DT-L, DT-O) and specifically on the climatology and statistics of the products.

The MODIS DT aerosol retrieval operates primarily by using observations from the seven RSBs with spatial resolutions of 0.5 km or finer. These bands are known as B3, B4, B1, B2, B5, B6 and B7; near λ=0.47 they are 0.55, 0.65, 0.86, 1.24, 1.63 and 2.11 µm, respectively. All are atmospheric window bands with minimal gas absorption. Additional RSBs and TEBs (at 1 km resolution) are used for tasks like cloud masking and snow identification. With this masking, the algorithm discards pixels unsuitable for aerosol retrieval and derives mean spectral reflectance (in the seven bands) that represents cloud-, snow- and ice-free dark-target scenes. Based on prelaunch signal-to-noise tests, global aerosol retrieval is optimized at 10 km (at nadir) spatial resolution, which is the resolution of the MxD04_L2 standard global product. Although there is also a more recent high-resolution 3 km aerosol product (MxD04_3K; Remer et al., 2013), here the MxD04 product (MOD04 or MYD04) refers to the standard (10 km) product, or to a level 3 aggregation of the 10 km products (MOD08 or MYD08).

The DT algorithm (both land and ocean) follows a lookup table (LUT) approach. This means that prior to retrieval, TOA spectral reflectance (in a subset of the seven bands depending on the surface) is simulated using scattering and radiative transfer codes (Wiscombe, 1980; Dubovik et al., 2002; Evans and Stephens, 1991; Ahmad and Fraser, 1982). These LUTs represent realistic combinations of aerosol, molecular and surface reflectance, which during the retrieval are compared with the observations. The solution is the LUT scenario (or multiple scenarios), which minimizes the cost function. From the LUT, one infers the total column loading (the aerosol optical depth or AOD or τ, reported at 0.55 µm), the spectral AOD (at multiple wavelengths), the Ångström exponent, AE or α) and estimates of the relative mixing between fine-sized (e.g., radius < 1 µm) and coarse-sized (radius > 1 µm) aerosol (known as fine-mode weighting or FMW or η). These retrieved aerosol properties, along with diagnostics describing the number of pixels used, the goodness of fit and the confidence in the retrieval product (quality assurance and confidence, known as QAC), are contained as separate quantities within the MxD04 product file. QAC ranges from 0 (no confidence) to 3 (high confidence) in a retrieval.

The details of the retrieval algorithm have evolved over time. However, for each collection, the same retrieval algorithms are applied to both MODIS-T and MODIS-A. For Collection 4 (C4), Remer et al. (2006) compared the two data sets and showed that they derived essentially the same monthly mean AOD over ocean. For C5 data, however, Levy et al. (2010) noted that there were discrepancies between the two data sets and that the MOD04 product appeared to have a high bias compared to ground-based AERONET data in 2003 and a low bias by 2008. There was no apparent overall bias to the MYD04 data. By 2013, the C5 Aqua products continued to show little or no apparent AOD trend over either ocean or land. However, Terra showed a −0.05 decade⁻¹ (−27 %) trend in global mean AOD over land (Lyapustin et al., 2014). Even more striking was that the differences between the two C5 time series was changing. In 2003, Over-land MOD04 showed higher AOD (e.g., offset of +0.02) than MYD04, but by 2013, MOD04 was lower (offset of −0.04). In other words, there was a trend in the offset (−0.06 decade⁻¹) over the period. While not changing in sign, the offset (MOD04-MYD04) over ocean also decreased, from +0.015 to +0.005 (−0.01 decade⁻¹). Such discrepancies between Terra and Aqua, including initial offsets, trends of the offsets and differences between land and ocean trends, were noted in many studies (e.g., Zhang and Reid, 2010; Remer et al., 2008; Yoon et al., 2012). Because the Terra and Aqua satellites have different viewing times over different regions of the world, convolved with the global spatial variability of aerosol distributions and diurnal cycles, we might expect offsets between the two MODIS time series. At the same time, these offsets may vary seasonally due to covarying diurnal cycles of aerosols and clouds (say, heavy dust or smoke being preferentially uncovered in either morning or afternoon). However, systematic trending of the offset is troubling. Considering the requirement that the AOD record should drift by less than 0.01 decade⁻¹ (GCOS, 2016), the differences between the C5 MODIS-T and MODIS-A products were unacceptable for deriving an aerosol CDR. To put this into perspective, a difference in 0.015 AOD is equivalent to ∼ 2–3 W m⁻² offset in estimating global direct aerosol radiative effect (Remer and Kaufman, 2006; Yu et al., 2006).

By the late 2000s, it was increasingly clear that, in addition to aerosols, other C5-derived data records were showing signs of non-physical trends (e.g., Lyapustin et al., 2014; Wang et al., 2012). The redundant on-board calibration protocol appeared to be insufficient for capturing degradation of the MODIS sensors, leading to artificial drifts in observed reflectance and subsequent derived geophysical parameters. To mitigate these drifts, MCST embarked on a new calibration protocol for Collection 6 (C6). In addition to regular observations of the moon and the on-board solar diffuser, MCST began monitoring observations over quasi-stable calibration desert targets, presumed to be nearly invariant (no rain, no changes in vegetation, etc.). Over such invariant Earth View (EV) targets, by compiling statistics of observed reflectance one could monitor long-term drifts in MODIS-observed reflectance. At the same time, the bidirectional reflectance function (BRDF) of such surfaces should be quasi-stable over time, so that, in addition to overall trending, MCST could characterize any trending response vs. scan angle (RVS) (Lyapustin et al., 2014; Toller et al., 2013; Sun et al., 2012; Wu et al., 2011). Corrections would be applied to any MODIS wavelength bands and scan angles that appeared to be drifting by more than 2 % since the first year of each mission. For MODIS-T, nearly all visible bands were drifting, with the shortest wavelengths drifting more rapidly. By the early 2010s, the shortest wavelength (e.g., blue bands) for MODIS-A also required correction.

In synergy with the overhaul of the upstream calibration method, the aerosol retrieval was updated for C6. Levy et al. (2013) introduced changes to the land/sea masking, the upstream cloud mask (e.g., MxD35; Frey et al., 2008) as well as the ancillary data inputs. There were also changes to the aerosol retrieval algorithm: some were made in response to upstream changes but others improved the physical aerosol retrieval. However, while there may be differences in the detrending coefficients applied to each sensor and a particular band may require detrending, the C6 aerosol algorithm is applied independently of upstream processing. The aerosol retrieval for C6, while different from C5, is applied in the same way to both sensors. This results in two sets of level 2 aerosol products (MxD04_L2) which are aggregated into daily (MxD08_D3) and monthly (MxD08_M3) gridded products, following the level 3 protocol for aerosol (L2 → D3 → M3; Levy et al., 2013). Assuming equal-area weighting (e.g., Levy et al., 2009), we further derive global monthly mean AOD, separately over ocean and land. From Fig. 1, we can observe that each pair of C6-derived AOD (at 0.55 µm) time series (land and ocean, separately) tracks more closely than their respective C5-derived versions (Lyapustin et al., 2014).

Yet, large positive offsets (MOD04 – MYD04) remain in C6. From Fig. 1, we see that over land this offset averages about +0.025 (approximately 12 % of the MYD04 global mean). While there is no significant overall trend to this offset, there are short periods of increase and decrease, and its variability appears to increase. A similar pattern was observed by Alfaro-Contreras et al. (2017). Over ocean, the offset averages +0.018, which is also ∼ 12 % of the MYD04 global mean. Although it is not plotted, we note that at 0.86 µm, the offset averages +0.014, which is ∼ 12 % of the global mean at that wavelength. The seasonal variability of the over-ocean offset is regular (maximum during northern summer months) throughout the time series, but there are identifiable periods of both increase and decrease. At 0.55 µm, the offset increased by ∼ 0.005 decade⁻¹ until about 2014 and then dropped suddenly in 2015. By 2017 the difference between sensors was at 2004 levels. While there are remaining (and puzzling) trends to the offsets, the magnitudes of those trends are less than 0.005 decade⁻¹, which suggest that the relative stability of the combined MODIS data records are approaching GCOS specifications for data product drift. From here on, we focus only on the general offsets and not on the trend of the offsets.

Figure 1 Time series of Collection 6, monthly global mean AOD (at 0.55 µm) over land (a) and ocean (b). For each panel, mean AOD (left axis) derived from MOD04 (Terra) is in red, from MYD04 (Aqua) it is in blue and the differences (right axis, MOD04-MYD04) are in black.

2.2 C6 Comparison with AERONET

The offsets to AOD and spectral AOD appear to be pervasive globally and have a significant enough magnitude to be of concern in the creation of climate data records. The accuracy of satellite retrievals is generally assessed through comparison with the ground-based sun photometer aerosol measurements from the AErosol RObotic NEtwork (AERONET; Holben et al., 1998), where the AERONET measurements are considered the “truth” with an uncertainty of ±0.02 (Ichoku et al., 2002; Petrenko et al., 2012). Based on scatterplots, the expected error (EE) of the global satellite product is an envelope that contains approximately two-thirds (or 1σ) of the collocated points. Levy et al. (2013) compared EE at coastal and island sites and estimated DT-O to be ±(0.04 + 0.10τ), where τ is the true AOD (±0.01) as observed by the sun photometer. Note that there are both absolute (±0.04) and relative (10 %) components for describing the EE of retrieved AOD. Compared at inland AERONET sites, Levy et al. (2013) compared EE at coastal and island sites and estimated DT-L to be ±(0.05 + 0.15τ), where both fixed and relative portions are larger than those for DT-O. Because the ocean surface optical properties are well constrained by the models (Cox and Munk, 1954; Koepke, 1984), the DT-O aerosol retrieval has a smaller EE than DT-L.

Table 1 Summary of scatterplots (not shown) of collocated MODIS and AERONET measurements of AOD at 0.55 µm, showing the performance of MYD04 (Aqua) and MOD04 (Terra) products relative to AERONET sun photometer data. Reported variables include N as the number of collocations, and %EE, % > EE and % < EE as percentages of collocations falling within, above and below EE envelopes. Bias is average difference (MODIS–AERONET), and slope, Y-INT, R² are parameters of least-squares linear regression.

While some AERONET sites exhibit a diurnal cycle, Kaufman et al. (2000) show that using AERONET data sampled at the MODIS passing time, the global AOD diurnal cycle is within 2 % of the daily mean AOD. This difference is at the same magnitude or much smaller than the discrepancies between Terra and Aqua retrievals that we discovered in this study, depending on the time span of AERONET data used. Although aerosol regional diurnal cycles may range more widely depending on location and/or season (Smirnov et al., 2002; Zhang et al., 2012), we expect global differences between morning and afternoon to be less than the offsets to MODIS seen here. To test this, we separately compare each MODIS data set to AERONET.

Table 1 summarizes the statistics of scatterplots (and linear regressions), comparing each MODIS C6 data set with appropriate AERONET data from the period 2003–2014. Obviously inappropriate collocation sites are excluded (e.g., Mauna Loa at elevation 3397 m being compared with sea-level retrievals), and data are filtered by QAC recommended by Levy et al. (2013). From these statistics, we see that, while both Terra- and Aqua-retrieved data sets perform similarly (R², RMSE) over different surface types, the overall bias (compared to AERONET) is larger for MOD04 than it is for MYD04. The magnitude of differences in bias (0.027 over land and 0.019 over ocean) are very similar to the overall global offsets we see in the Fig. 1 time series plots. Also, we note that these differences in bias are very similar to a recent MODIS 3 km product comparison (Gupta et al., 2018).

We note that there are more collocations (N) for MOD04/AERONET than MYD04/AERONET. Over land, this is consistent with King et al. (2013), showing that MODIS-T observes smaller cloud fractions over land. However, King et al. (2013) also reports larger cloud fractions observed by MODIS-T over ocean, which is not consistent, except that most AERONET sites are along coastlines.

In coordination with Table 1, Fig. 2 provides the MODIS–AERONET differences as a function of AERONET-measured AOD. The bins of AERONET AOD are set so that there are a nearly equal number of points in each bin. Over both land (panel a) and ocean (panel b), the biases for MOD04 (red dots/shaded envelopes) are larger than that of MYD04 (blue dots/shaded envelopes). Over ocean, both products appear to have positive biases at low AOD over ocean, which is due to DT-O not being allowed to retrieve zero or negative AOD (leading to automatic positive bias). Over land, while both products have a median positive bias, it is larger for MOD04. The difference in the median bias (blue and red dots) is relatively constant across all AOD bins over land but increases with AOD over ocean. Overall, the difference in median bias (each MxD04 collocated with AERONET) is roughly equal to the overall offset in AOD between the two MxD04 time series. The statistics of each MODIS-retrieved AE compared to AERONET are very similar, as shown in panel c, except that AE from MOD04 is lower (by about 0.05) compared to MYD04.

Figure 2 Statistics of MxD04–AERONET difference as a function of AERONET values for AOD (0.55 µm) over land (a) and ocean (b), and for AE over ocean (c). For each panel, data from MOD04 (MYD04) are plotted in red (blue). For each sensor, the dots (and connecting lines) represent the mean of the MxD04–AERONET difference, whereas the shaded area represents the middle ±1σ of the difference. Note that each MxD04 is compared separately to AERONET and that the AERONET data are ordered into bins with an equal number of points.

2.3 Spatial distribution of C6 offsets

It is still possible that offsets between the two MODIS time series are tied to unequal sampling of heavy aerosol events across the globe. To compare offsets as a function of location, we focus on 2008 data. Here, we derive a yearly mean AOD (per grid box) from monthly mean data, assuming valid data in at least 2 months (e.g., L2 → D3 → M3 → Y3). Figure 3 displays the 2008 annual mean for MYD04 AOD (at 0.55 µm) over both land and ocean, as well as the absolute and relative differences between yearly MOD04 and MYD04. Note that, instead of 1^∘ × 1^∘ aggregations like standard MODIS L3 data (MxD08_M3), we derive them at 0.5^∘ × 0.625^∘ resolution (to use in Sect. 3).

Except for the well-known aerosol hotspots (African dust/smoke, Asian pollution/dust, etc.), Fig. 3a shows that most of the globe experiences annual mean low AOD (τ < 0.1). Figure 3b shows the gridded absolute differences (MOD04 – MYD04) and Fig. 3c shows the relative differences. Over most of the globe the absolute differences are 0.015–0.025, showing that the global mean values for 2008 (seen in Fig. 1) arise from a global distribution of offsets of the same small magnitude rather than a residual of widely fluctuating large positive and negative offsets. There are, however, some notable areas with the opposite polarity, primarily regions of known biomass burning (Amazon, southern Africa).

Figure 3 Gridded (0.625^∘ × 0.5^∘) global mean AOD (at 0.55 µm) for 2008, derived from MYD04 (a), the difference between MOD04 and MYD04 (b) and the relative difference (c).

As mentioned in Sect. 2.1, the DT retrieval over ocean reports spectral AOD. Using annual mean AOD at 0.86 and at 0.55 µm, we derive the annual mean Ångström exponent (AE) at each grid. Gridded AE and AE differences are shown in Fig. 4. We observe lower AEs (larger relative particle size) over well-known dust belts as well as over the most remote ocean. Higher AEs (smaller particle size) are observed where there is continental pollution or smoke outflow. Figure 4b shows the differences in mean AE (MOD04 vs. MYD04), showing that, with a few exceptions, MOD04 consistently reports smaller AE (larger particles) than MYD04 by about 0.05.

Figures 3 and 4 demonstrate that the AOD from MODIS-T is consistently higher than MODIS-A and AE is consistently lower, and the uniformity of the offsets is suspicious. However, it is interesting that the southeastern Atlantic downwind of the southern African savanna (Meyer et al., 2013) shows opposite AOD and AE offsets.

Figure 4 Gridded (0.625^∘ × 0.5^∘) global mean AE (at 0.55/0.86 µm) for 2008 derived from MYD04 (a) and the difference between MOD04 and MYD04 (b).

3 Using modeling to study morning versus afternoon offsets

Section 2 identified significant differences between the aerosol products derived from Terra and Aqua. There could be many causes for these discrepancies, from instrument calibration and sampling to physical causes. While MODIS-T and MODIS-A processing are identical, differences could arise from differences in orbits and satellite overpass times. Terra is in descending orbit with a daytime equator crossing (southward) at 10:30 local solar time and Aqua is in ascending orbit (northward) at 13:30 LST. Because of the different headings, although the local overpass time difference is 3 h at the equator, it is closer to 1.5 h in the Northern Hemisphere (NH) midlatitudes and 4.5 h in the SH (Fig. 5). Therefore, it is possible that different aerosol statistics might arise, whether due to diurnal cycles of aerosol or clouds (leading to different sampling). Interestingly, because of the symmetry of the orbits (±N h from local noon at every location), the actual geometrical sampling (e.g., statistics of solar zenith, relative sun/sensor azimuth and resulting scattering and glint angles) of the two sensors is very similar. Indeed, if computing average angles over the entire year, there is on average only a 0.8^∘ difference in the solar zenith angle (MOD04 < MYD04), and a 0.3^∘ difference in the scattering angle (MOD04 > MYD04). This means that, although the aerosol retrieval may have biases as function of angle (e.g., Hyer et al., 2011), the symmetry of Terra and Aqua orbits should not lead to a consistent difference in retrieved AOD.

Figure 5 Gridded average MODIS local observation time (LST) for Aqua (a), Terra (b) and the difference between the two (c).

Cloud types and cloud properties show significant diurnal variation (e.g., Eastman and Warren, 2014). In fact, King et al. (2013) catalogued differences between cloud statistics from the two MODIS sensors, which can be repeated using C6 MODIS data. To explore the differences identified in Sect. 2 and to discount the possibility of diurnal sampling (related to cloud fraction differences) being the root cause of the Terra–Aqua offset, we use results from aerosol simulations performed with the NASA Goddard Earth Observing System, version 5 (GEOS-5) Earth system model (Molod et al., 2015). GEOS-5 is run here in a replay mode, using winds, temperature and pressure fields from the recent Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2) joint aerosol and meteorological reanalysis (Gelaro et al., 2017; Randles et al., 2017). Using MERRA-2 meteorological constraints ensures simulation of real weather events and realistic cloud fields. The model is run globally at a c180 horizontal resolution (∼ 0.5^∘ × 0.625^∘ latitude × longitude resolution) on its cubed-sphere native grid and produces high time resolution (hourly) aerosol output based on the prognostic Goddard Chemistry, Aerosol, Radiation and Transport (GOCART) module (Colarco et al., 2010), which is run online and radiatively coupled. GOCART simulates the mass of various aerosol species (dust, sea salt, sulfate and carbonaceous aerosol) at every model time step. Diurnal cycles in the aerosol distributions arise through prescribed diurnal variability in emissions (e.g., biomass-burning emissions tend to peak in the afternoon) or meteorology. Conversion of the simulated mass to AOD is accomplished through precomputed LUTs that include mass extinction efficiencies as a function of species, particle size and hygroscopicity. The model run used here does not invoke aerosol data assimilation and so is essentially a chemical transport model driven by reanalysis. Unlike satellite products, which have gaps due to swath sampling or decisions by the retrieval algorithm, the model has no gaps in its computed AOD field. There is aerosol under clouds, in glint and over bright surfaces including ice and snow. Also, unlike the satellite, the model is not limited to polar-orbiting overpass times of 10:30 and 13:30 LST. However, the beauty of the model is that we can sample the outputs in any way we want, including in a satellite-like manner (e.g., Schutgens et al., 2016; Colarco et al., 2014).

Consider the partial orbits of Aqua and Terra around 12:00 UTC, specifically, the samplings of each MODIS between 11:30 and 12:30 UTC (Fig. 6). The light colors represent the swath of the MODIS track, and the dark colors represent where AOD was retrieved. Not only do the two tracks cover different parts of the world, each DT product retrieves less than 10 % (due to clouds, glint, bright surfaces, etc.) of the possible opportunities along the swath. Repeating the analysis of Fig. 6 at each hour for the entire year 2008 leads to two aggregations of the model for each satellite: the first represents the full MODIS swaths and the second represents the retrieval of the MODIS-DT products.

Figure 7 is analogous to Figs. 3 and 4, except that instead of aggregating MYD04 versus MOD04 products, we have aggregated MERRA-2 outputs along the full MODIS swaths. Figure 7a and c show global AOD (at 0.55 µm) and AE (0.55 vs. 0.86 µm) at the afternoon (PM) overpass, analogous to global MODIS-A swath sampling in 2008. The aerosol hotspots are obvious, and most of the globe has low AOD. Figure 7b and d show the AOD and AE differences between MERRA-2 as if they were sampled during the morning (AM; like MODIS-T) instead of the afternoon (like MODIS-A) observation time.

From Fig. 7, we can make some generalizations. First, the general patterns of the afternoon AOD (Fig. 7a) are similar to the aggregated MYD04 DT data (Fig. 3a). However, there are no gaps, because there are no DT retrieval decisions (e.g., masking). Most importantly, unlike the MODIS retrieval product (Fig. 6b), there is no overall AM–PM offset to the AOD (Fig. 3b). There are, however, regional differences to the offsets. Morning AOD tends to be lower for the biomass-burning regions over land, which is expected due to diurnal cycle of fire emissions (Boschetti and Roy, 2009). Over the ocean, there is even less variability from zero offset. For over-ocean AE (Fig. 7c), although the general patterns are similar to the MODIS retrieval products (Fig. 4a), the model outputs show lower AE, suggesting that the model has simulates larger particles than the retrievals. The outflow from equatorial Africa is one exception; the model reports much higher AE than the MODIS product, suggesting finer-sized particles. In addition, there are small positive and negative AM–PM differences, with no apparent systematic pattern. Generally, by comparing the model sampled in the morning to the afternoon sample, we see little evidence of global offsets to either AOD or AE.

Figure 6 Map of hourly swath and retrieval aggregation during ±30 min of 12:00 UTC on 28 May 2008. MODIS-T and MODIS-A swaths are in light red and blue shading, whereas retrieved pixels are a darker color. The arrows represent the direction of the satellite orbit across the equator (descending for Terra, ascending for Aqua).

Figure 7 Gridded (0.625^∘ × 0.5^∘) global mean AOD (at 0.55 µm – a, c) and AE (0.55/0.86 – b, d) µm for 2008, derived from sampling of MERRA-2 along the MODIS swaths. (a, b) Derived from PM sampling (like MODIS-A). (c, d) Difference between AM (MODIS-T) and PM (MODIS-A) swaths.

Figure 8 Gridded (0.625^∘ × 0.5^∘) global mean AOD (at 0.55 µm – a, c) and AE (0.55/0.86 – b, d) µm for 2008, derived from sampling of MERRA-2 along the MODIS retrievals. (a, b) Derived from PM sampling (like MYD04). (c, d) Difference between AM (MOD04) and PM (MYD04) sampling.

However, due to clouds, glint and bright surfaces, less than 10 % of the area sampled is actually retrieved by the DT algorithm. Due to differences in cloud fraction between morning and afternoon orbits (e.g., King et al., 2013), there may be systematic differences in the aerosol sampling. For example, while heavy smoke conditions (high AOD, high AE) may be present throughout the day, preferentially cloudy conditions (AM or PM) would affect the sampling of these aerosol events and thus the AM–PM offsets. Figure 8 is analogous to Fig. 7 but represents the model being sampled where and when there is AOD reported in the MYD04 or MOD04 products. The overall AOD and AE patterns (e.g., Fig. 8a and c) are much like those from the entire swath (Fig. 7a and c) but with gaps exactly like the satellite retrievals (Figs. 3a and 4a). Comparing Fig. 8b and d with their counterparts in Fig. 7 shows that, by imposing satellite sampling, the variability of both AOD and AE offsets more resembles the satellite regional distributions. However, the mean offsets to both AOD and AE have not increased due to the imposition of satellite sampling. Sampling alone cannot explain the overall offset seen in the satellite data products.

Except for the Amazon region, where both show negative offsets, there is not much resemblance between the AOD differences shown in Figs. 8c and 3b. Since Fig. 8c represents the expected offsets, the overall positive offset in Fig. 3b has probably masked some of the diurnal cycles expected in that figure. The overwhelming positive offset in Fig. 3b, especially over the oceans where the model shows very little difference, indicates there is a systematic difference in the two retrievals that could only be attributed to instrument calibration.

4 Testing calibration corrections

The C6 MODIS products report persistent systematic offsets in the AOD and AE that cannot be explained by diurnal sampling differences, as was explored in the modeling exercise of Sect. 3. The next possible explanation for the offset is calibration. As explained by Lyapustin et al. (2014), although the MODIS Characterization Support Team (MCST) updated MODIS calibration to account for the severe trending observed in C5 data, there still may be offsets in C6.

4.1 C6+ corrections

The MODIS DT retrieval algorithm is an inversion on multi-spectral data. The reality is that, due to the retrieval being a multi-channel inversion, changing one wavelength at a time leads to nonlinear changes in retrieved AOD and AE. However, the over-ocean-retrieved AOD is most sensitive to changes in the 0.86 µm channel (B2), because of the requirement that measured reflectance must exactly match retrieved reflectance at that wavelength. Changes to the 0.55 µm channel (B3), in turn preferentially impacts the retrieved AE. Over land, retrieval of AOD is most sensitive to changes in the blue (B3 = 0.47 µm) band.

The spectral channels used in the retrieval algorithm are calibrated independently for each sensor and may drift differently over time. Based on monitoring bidirectional-reflectance function (BRDF) over the same pseudo-invariant (remote desert) surfaces as used by MCST, Lyapustin et al. (2014) devised a method for correcting the L1B reflectances. This method, known as C6+, accounts for changing sensitivity to polarization, corrects (as a function of wavelength band) residual trends in both Terra and Aqua, and then applies cross-calibration, assuming Aqua to be the more stable and better-characterized sensor. Thus, C6+ can be applied directly to the C6 L1B data, offering a corrected L1B that can be substituted for the standard L1B and be used to create alternative L2 data. In fact, the C6+ corrected data are already being applied upstream of the MODIS land-surface retrieval products in C6. There are corrections to both MODIS-T and MODIS-A data.

Figure 9 Impact of applying C6+ calibrations to MODIS-T reflectance data on 25 July 2008 at 16:25 UTC. Absolute differences in each wavelength band.

The C6+ calibration involves three steps. The first step is conducting polarization correction (e.g., Meister et al., 2014, 2005; Kwiatkowska et al., 2008). Polarization correction is complicated, because there are both angular (dependence on scan angle and across-track position) and mirror side/optics dependencies (dependence on scanline/detector and along-track position). The corrections may even be of opposite polarities depending on position across-track and along-track. At the same time, the polarization correction is dependent on the scene itself. Rayleigh-scattering (molecular)-dominated scenes (minimal aerosol over dark surfaces) require the largest relative correction. After polarization correction and for each sensor, C6+ assumes that performance was optimal at the beginning of its mission, and that MODIS-A is more stable than MODIS-T overall. Using the quasi-stable desert scenes corrects for the drifts as well as the initial offsets.

Lyapustin et al. (2014) presented formulas (polarization + detrend + cross-calibration) to correct four of the bands (B1, B2, B3, B4) used in the DT retrieval, plus B8 (0.41 µm) used in the deep-blue retrieval (Hsu et al., 2012). More recently this team (Yujie Wang, personal communication, 2016) expanded the correction to include the remaining bands (B5, B6, B7) used in the DT retrieval. The magnitude and polarity of the correction for any MODIS-T or MODIS-A pixel depends on the wavelength band, scan angle, across-track position, along-track detector/scan and the scene itself. Due to the complicated nature of the C6+ correction, and its convolution with the non-linearity of the DT aerosol retrievals over land versus over ocean, we cannot easily perform sensitivity tests. Therefore, we have chosen to use brute force and have applied the C6+ correction upstream of our C6 aerosol retrieval algorithm. We use the same operational processing structure as MODAPS but substitute the corrected L1B data for the archived (LAADS) data.

Figure 10 Impact of applying C6+ calibrations to MODIS-T reflectance data on 25 July 2008 at 16:25 UTC. Relative differences in each wavelength band.

Figures 9 and 10 show absolute difference (${\mathit{\rho }}_{\mathrm{C}\mathrm{6}+}-{\mathit{\rho }}_{\mathrm{C}\mathrm{6}}$) and relative (${\mathit{\rho }}_{\mathrm{C}\mathrm{6}+}-{\mathit{\rho }}_{\mathrm{C}\mathrm{6}})/{\mathit{\rho }}_{\mathrm{C}\mathrm{6}}$ difference to reflectance when applying the C6+ correction for a mostly clear-sky case over the Caribbean Sea and the Gulf of Mexico (25 July 2008). Note that these reflectance units are different to those reported in L1B, because they are from the L2 aerosol product and are normalized for solar zenith angle. Note also that reflectance is only plotted where an aerosol retrieval was attempted, so cloud-masked pixels are not plotted (as per the DT aerosol retrieval). At the same time, we intercept the reflectances before the glint mask is applied during the aerosol retrieval, so one can see the glint patterns in the reflectance.

Figure 11 Impact of applying C6+ calibration corrections to the MxD04 AOD (a, c) and AE (b, d) products. Panels (a, b) show impact on MOD04 (Terra), whereas panels (c, d) show impact on MYD04 (Aqua).

Figure 12 Histograms of the C6+ corrections on the Terra–Aqua AOD differences during 2008 over land (a) and ocean (b). Histograms are derived from comparing gridded AOD.

The C6+ correction provides only a tiny relative change (0.2 %) to B2 (0.86 µm) over this scene. Whether over land, dark ocean or glint (near center of image), the C6+ makes approximately the same overall correction to B2. However, for B3 (0.47 µm), the C6+ correction has scan-angle dependence and leads to reduced (most of the image) or increased reflectance (just to the left of glint). Some of these changes approach 3 %. The relative patterns of decrease and increase shift from wavelength to wavelength, although most reflectance values are reduced using C6+. The high scattering angle dependence of correction for shortwave channels (0.47–0.65 µm in Fig. 9a, b and c) indicates that correction is mainly due to polarization, while for a longer wavelength, the correction is more homogenous over the entire scene, mainly due to detrending of these channels.

These correction characteristics clearly should have implications for the retrieval of global AOD and AE during 2008. We apply the C6+ corrections to both MODIS-T and MODIS-A data sets for the entirety of 2008. Figure 11 shows the differences to the global AOD when C6+ is applied to each sensor. There is very little overall change to the over-ocean AOD in either MODIS data set. Over land, both sensors show a reduced AOD with C6+ but with a larger decrease in MODIS-T. Since plotting gridded C6+ offsets will be indistinguishable from Figs. 3b or 4b, we plot histograms of the offsets (Fig. 12). Over ocean, because of only tiny changes to MODIS-T in the B2 band, the gridded annual mean differences have not budged. Over land, after a decrease in both sensors' AODs, the peak of the offsets decreased by ∼ 0.01 – in the right direction.

Because the C6+ calibration did not have a large impact on the 0.86 µm band, the distribution of the AOD offsets (MOD04-MYD04) did not decrease significantly over ocean. However, due to more significant changes to the 0.55 and 0.65 µm bands (Fig. 9), AE changed in unequal and opposite directions for each sensor (Fig. 11c and d). This significantly changed the AE offsets (Fig. 13). Instead of being consistently negative in C6 (Fig. 4b), there are both regions of positive and negative AE offsets (Fig. 13) after the C6+ correction. Smoke regions (especially in the tropics) show larger AE (smaller particles) for MODIS-T versus MODIS-A, while dust regions and the extratropics show smaller AEs (larger particles) for MODIS-T. This AE offset pattern, while averaging closer to zero, still has significantly more variability than that expected by the model.

4.2 Other corrections

Applying the C6+ calibration correction appears to have reduced the average AOD offset over land. Over ocean, the average AOD offset was unchanged, but the average AE offset decreased. The C6+ was based on obtaining BRDF over pseudo-invariant desert sites. Other investigators have attempted to improve the C6 calibration, comparing measurements over other types of surfaces or scenes, including a site in Antarctica (Dome C) and over deep convective clouds.

For example, Doelling et al. (2015) compared the two MODIS sensors using observations from nearly simultaneous nadir overpasses (NSNOs). Although Terra and Aqua have orbits in opposite directions, their orbits nearly cross each other 14 times per day. Nadir crosses only happen at ∼ 68.3^∘ latitude; however, if small angular and time differences are tolerated, then off-nadir comparisons and corrections can also be performed. The NSNO method, while straightforward, assumes all differences are radiometric and not polarization or angle dependent. Nonetheless, this method provides time-dependent cross-sensor coefficients, effectively tuning MODIS-T to MODIS-A.

Figure 13 Impact of the C6+ corrections on the Terra–Aqua AE (0.55/0.86 µm) differences during 2008. Compared to Fig. 4b, the overall global bias is reduced.

Overall, the mean scaling providing Table III of Doelling et al. (2015) appears to be similar in polarity to the cross-calibration factors provided by Table 3 of Lyapustin et al. (2014). For example, both studies suggest that in B3 (0.47 µm), MODIS-T has a high bias (by ∼ < 1 %) and should be multiplied by ∼ 0.991. Both studies also suggest MODIS-T has a low bias (by ∼ 1.5 %) in B1 (0.65 µm) and should be multiplied by ∼ 1.015. However, they appear to differ in their B2 (0.86 µm) corrections, with Lyapustin et al. (2014) suggesting they should be multiplied by 1.006 but Doelling et al. (2015) suggesting 0.994. Of course, these apparent discrepancies may be canceled out due to different methods: Lyapustin et al. (2014) detrends each sensor independently and then applies cross-calibration, whereas Doelling et al. (2015) convolves the two processes. Nonetheless, both papers suggest that B2 needs correction.

We used Table III of Doelling et al. (2015) to estimate scaling coefficients appropriate for 2008, and we tested them by applying them directly within the aerosol retrieval (rather than the upstream C6+ code). When reading L1B, we applied coefficients (based on nadir NSNO) to each band of MODIS-T. Except for B5 (1.24 µm), which we reduced by approximately 3 %, and B1, which was increased by 1.5 %, most other bands were adjusted by less than 1 % either way. Overall, the results were similar to those when applying C6+. The mean AOD offset over ocean remained, but it was reduced over land (but not to zero). The mean offset for AE was adjusted toward zero, but the spatial patterns of Fig. 4b generally remained.

Figure 14 Differences between the C6.1 and C6 MOD04 (Terra) product for AOD (a) and AE (b). Note that the color scales are identical to those in Figs. 3b and 4b.

5 Collection 6.1

The radiometric calibration for C6 was based on the combination of prelaunch, solar diffuser, moon observations and selected targets on Earth. When C6 processing began in 2012, calibration coefficients were derived in order to smoothly connect the beginning of the mission through to 2012 (reprocessing). As long as the instrument performance did not change too quickly, the C6 methodology could be used for forward processing. However, by early 2016, some of the TEBs for MODIS-T were becoming unusable.

Since the beginning of our analysis of the C6 differences between Terra and Aqua aerosol products, the MODIS team has released an updated collection denoted by C6.1 (https://modis-atmosphere.gsfc.nasa.gov/documentation/collection-61, last access: 27 June 2018). The primary purpose of C6.1 was to correct for the TIR issues which had resulted in a failure of the standard MODIS cloud mask algorithm, affecting all other downstream algorithms using MOD35 as input. This included the DT aerosol retrieval for MOD04. The DT team used the opportunity to make modest improvements to the retrieval algorithm, including (a) to include the corrections for urban surfaces (Gupta et al., 2016), (b) to revise the logic regarding detecting or rejecting ocean pixels using the 1.63 µm band and (c) additional diagnostic changes that affected the output-retrieved AOD or AE. Therefore, if applied to the same C6 L1B input (not corrected with C6+), there are minimal global differences between products of the C6.1 and C6 DT aerosol algorithms. However, with the recent completion (December 2017) of C6.1 reprocessing for MODIS-T (including 2008), we tested whether the updates to upstream L1B and cloud mask (MOD35) would together affect the MOD04 data and therefore help to reduce the MOD04 – MYD04 offsets in the C6 products. While there were also changes to calibration coefficients for the MODIS-T RSB bands, we can confirm that the difference to L1B reflectance is negligible for 2008 data. Therefore, any global differences between C6.1 and C6 aerosol products would be dominated by the upstream pixel selections and not by RSB calibration.

Figure 14 shows the differences (panel a: AOD and panel b: AE) between C6.1 and C6 for MOD04 (Terra) during 2008, showing only small changes in global AOD and AE. On average, AOD increased over ocean by about 0.001 and decreased over land by similar magnitude, which are much less than the desired changes (e.g., Figs. 3b and 4b). The changes from C6 to C6.1 may alleviate some of the AOD offsets over land in 2008 (maybe as much as 10–20 % of the bias in some places) but will exacerbate the bias over ocean by about the same percentage. Likewise changes from C6 to C6.1 will only affect the annual mean biases in AE by 10–20 %, both positively and negatively. The changes introduced by C6.1 are just too small to eliminate the Terra–Aqua differences that were identified and explored in the analysis presented above.

We note that the changes to C6.1 L1B products are temporally dependent, so we might expect larger differences between C6 and C6.1 in later years (especially after 2015). We also have not yet analyzed C6.1 MYD04 data (Aqua C6.1 reprocessing began on 28 December 2017). Thus, there may be slightly different consequences to the aerosol products than are shown here for 2008. However, given the small magnitude seen in 2008 and expected through the entire time series, it is unlikely that the C6.1 changes will provide the fix necessary to bring the Terra and Aqua aerosol products into agreement.

6 Discussion and conclusion

The DT aerosol retrieval has been applied to MODIS-T data since 2000 and MODIS-A data since 2002. Time series of the C6 products (MOD04 and MYD04) are almost in lockstep (Fig. 1). However, compared to MODIS-A (afternoon overpass), the global mean MOD04 (morning overpass) shows consistently higher AOD at 0.55 µm (by ∼ 0.015–0.02 or ∼ 13 %) over both land and ocean. At the same time, there is a 0.005 decade⁻¹ trend to this offset over ocean and increasing seasonal variability over land after 2011.

Focusing on 2008, we studied the AOD offset. Over ocean, the offset appears everywhere, regardless of the overpass time difference (4.5 h in SH midlatitudes, 1.5 h in NH). Over land, there is more variability in the offset, but only known biomass-burning regions display a negative offset (morning AOD is lower than afternoon). Over ocean, we also see that there are consistent offsets in the spectral AOD, as demonstrated by the Ångström exponent (AE), showing that MODIS-T reports lower AE globally by about 0.05.

We used the GEOS-5 replay model output to question the observed global offsets in AOD and AE. When we sample the model along the MODIS swaths (Terra and Aqua separately, then take the differences), most of the globe appears to have no AOD or AE offsets. However, we might expect to observe negative offsets (morning AOD lower than afternoon) in the biomass-burning regions. As we sample the model only for the MxD04 retrievals (MOD04 and MYD04 separately and then take the differences), we see increasing variability in both AOD and AE offsets. This is due to differences in cloudiness between the morning and afternoon, which is convolved into the MODIS data. There is the suggestion of a more generalized offset to AE but not approaching the magnitudes seen from the satellite retrieval products.

Of course, we cannot yet rule out other physical reasons for the offsets. For example, although the retrieval algorithm corrects for gas absorptions (column water vapor, ozone, etc.; Patadia et al., 2018), unaccounted differences between morning and afternoon (for example if 12:00 UTC water vapor was assumed for both 10:30 and 13:30 overpasses) could lead to systematic biases in retrieved AOD. We should consider that the aerosol optical properties themselves (e.g., refractive index, size/shape distribution) could be wrong and also lead to generalized AOD bias compared to AERONET (e.g., Ichoku et al., 2003; Eck et al., 2013). Additionally, if there were differences in optical properties which were not accounted for and due to very late-morning cloud processing (e.g., Eck et al., 2012), one might see offsets between AM and PM, and different offsets between AM and PM versus AERONET. Although the current modeling framework (e.g., our MERRA-2 sampling) does not suggest that it leads to a global offset, this is definitely a topic for further study.

Since the Terra–Aqua bias is so similar to the difference between Terra–AERONET and Aqua–AERONET (Terra–Aqua = Terra–AERONET − Aqua–AERONET), we suspected the MODIS calibration. We tried two alternative calibration efforts, both of which could be applied upstream of the aerosol retrieval. The first, known as C6+ (Lyapustin et al., 2014), included polarization and angular correction for each sensor, a detrend for each sensor, and then cross-calibration to normalize Terra to Aqua. The second (Doelling et al., 2015) does not apply a polarization correction and convolves the detrending and cross-calibration into scaling factors. Each method was applied upstream of the aerosol retrieval using 2008 data. Both methods reduced the overall AOD offset over land from 0.02 to 0.01 but did not significantly affect the overall offset over ocean. The mean negative offset for AE was reduced to zero; however, this led to positive offsets in AE for smoke outflow regions. After calibration, this would translate to MODIS-T interpreting small particles as being larger than MODIS-A, and interpreting large particles as being smaller than MODIS-A. That the two calibration efforts did not remove the offsets entirely, however, does not mean that calibration is not the culprit. Not only is calibration in the bands used for the aerosol retrieval (e.g., B1–B7), but thermal infrared channels and 1.38 µm bands are used for cloud detection and masking. Clearly, more analysis is required.

The MODIS Atmosphere Science Team recently began processing the C6.1 family of products, primarily to address issues related to thermal infrared bands and impacts on the standard cloud mask for MODIS-T. There was no major change to the methodology of the MCST reflective band calibration for MODIS-T. Except for improvements over urban regions, the C6.1 aerosol retrieval is also nearly unchanged. Thus, based on comparing the C6.1 aerosol product with C6 during 2008, we expect there to be no change to the overall offsets to both AOD and AE. However, since 2011 (beginning of C6 processing), additional reflective bands (on both MODIS-T and MODIS-A) have strayed more than 2 %, so there are revisions to overall calibration that may show apparent effects of the two time series in later years (well after 2008) (https://modis-atmosphere.gsfc.nasa.gov/sites/default/files/ModAtmo/C061_L1B_Combined_v10.pdf, last access: 27 June 2018). Thus, although we expect continued overall offsets between Terra and Aqua DT aerosol products, the trend and variability of the offset (e.g., Fig. 1) may change. The Science Team will continue to monitor, compare and attempt alternative calibrations for the Terra and Aqua aerosol products to the end of the satellite missions. At the same time, we will test the aerosol retrieval with new versions of C6+ or other types of additional corrections and determine whether offsets/biases/trends of the aerosol product can be reduced for future collections.

In the meantime, users of the products should not interpret differences between Terra and Aqua aerosol products as representing a true diurnal signal in the aerosol, unless magnitudes of the observed signal greatly exceed the biases described here. However, because collocated comparisons between MODIS aerosol retrievals and AERONET observations show Terra with a larger high bias, the recommendation is to rely more on Aqua retrievals for quantitative long-term climate-related applications. On the other hand, we note that the bias in AOD is only ∼ 0.02, which is the noise level for short-term applications such as air quality forecasting, and thus both Terra and Aqua aerosol products provide adequate quantification for these types of uses.