2.1 Relevant sensor characteristics

Sayer et al. (2016) developed the AAC retrieval algorithm with a goal of implementation being as similar as feasible across the different sensors, relying on only those bands common to the three instrument types. SeaWiFS (McClain et al., 2004), MODIS (Barnes et al., 1998), and VIIRS (Cao et al., 2013) are all passive broad-swath imaging radiometers. SeaWiFS operated from late 1997 to December 2010; MODIS provides data on the Terra platform from late February 2000, MODIS on the Aqua platform from July 2002, and VIIRS on the Suomi National Polar-orbiting Partnership (SNPP) from March 2012. Both MODIS sensors and SNPP VIIRS are still operational. SNPP VIIRS was followed with an additional VIIRS sensor launched in late 2017 (not considered in this study), and more are scheduled for the future.

SeaWiFS measured reflected solar radiation at the top of atmosphere (TOA) in eight bands with centres from 412 to 865 nm; MODIS and VIIRS have additional solar bands, as well as thermal infrared (tIR) channels. The AAC retrieval relies on common bands centred near blue (470 nm for MODIS, 490 nm for SeaWiFS and VIIRS), green (550 nm), red (650 nm for MODIS, 670 nm for SeaWiFS and VIIRS), and near-infrared (nIR, 865 nm) wavelengths. These calibrated and geolocated measurements are referred to as level 1b (L1b) data. Note that in this study these approximate wavelengths and/or band colour names (e.g. green for 550 nm) are sometimes referred to for simplicity, although all radiative transfer (RT) calculations use full sensor relative spectral response (RSR) functions. Specifically, the TOA reflectance ρ for band i is defined as

where L_λ is the spectral radiance passing into the satellite field of view at TOA; E_λ the downwelling solar spectral irradiance at TOA, perpendicular to the Sun and at 1 astronomical unit (AU); and Φ_i the sensor RSR for band i, all functions of wavelength λ. The factor D_⊙ is the Earth–Sun distance in AU (variable throughout the year) and μ₀ the cosine of the solar zenith angle, which affect the total solar radiation received. Note that L_λ and so ρ_i depend on solar/observation geometry (and of course surface and atmospheric state), omitted here for simplicity of notation.

Figure 3Fraction of days a given point (longitude) in the retrieval domain (Fig. 2) is observed at least once by the individual SeaWiFS, MODIS, and VIIRS sensors, as a function of latitude.

Download

For MODIS, nominal horizontal pixel sizes vary from 0.25 to 1 km (dependent on band); here, the finer-resolution bands are aggregated and coregistered to 1 km. For VIIRS, the nominal pixel size for the relevant bands is 0.74 km. For SeaWiFS, pixel sizes are 1.1 km but on board resampling performed for Global Area Coverage (GAC) mode subsamples these to provide an effective horizontal resolution of ∼4.4 km at nadir. As GAC data are a subsampling rather than an average, it is most appropriate to consider these as 1.1 km pixels with gaps between them (as opposed to the continuous swaths of MODIS and VIIRS). All quoted pixel sizes are for nadir viewing geometries, at which pixels are approximately square. Away from nadir the pixels enlarge, begin to overlap, and become more distorted in shape due to the scan geometry and Earth's curvature. This distortion is largest for MODIS (Sayer et al., 2015a) and smallest for VIIRS (Wolfe et al., 2013). Swath widths are 1502 km for SeaWiFS (GAC mode), 2330 km for MODIS, and 3040 km for VIIRS (meaning VIIRS has no inter-orbit gaps). Around the Equator SeaWiFS also tilted to decrease the fraction of the swath affected by Sun glint in each hemisphere; this tilt led to several scan lines near the Equator being missing. Figure 3 shows the resultant fraction of days when each of the sensors sampled a given location within the region, as a function of latitude. For SeaWiFS, coverage over the core of the stratocumulus deck (Fig. 2) was obtained on about 60 % of days, potentially leading to larger sampling biases than the other sensors. For MODIS this value is closer to 85 % at the Equator and becomes 100 % poleward of ∼25^∘; for VIIRS, the whole region is imaged at least once per day. All the sensors are on platforms in Sun-synchronous polar orbits; MODIS Terra has a daytime equatorial crossing time of 10:30 (local solar time) while MODIS Aqua and VIIRS have an overpass time of approximately 13:30 local solar time. SeaWiFS crossed near local noon at launch, although it drifted in the later years of the mission (and ended around 13:30–14:00 in 2010). It is possible that these differences in overpass time will lead to differences in the retrieval results; for AAC cases, however, downwind travel takes places over periods of several days, and so it is unlikely that, far from sources, the AOD will have changed significantly between satellite overpasses.

The DB aerosol retrieval algorithm (Hsu et al., 2013) has also been applied to all these sensors to retrieve AOD for cloud-free scenes over land. The main data product from DB is the AOD at 550 nm; in this study mentions of AOD without a specific wavelength indicated refer to 550 nm. For the SeaWiFS and VIIRS applications of DB (but not MODIS, at present), a Satellite Aerosol Retrieval Algorithm (SOAR) is applied over water surfaces to provide a near-global picture (Sayer et al., 2012, 2018a, b). This combination of DB and SOAR is often colloquially referred to as the Deep Blue data product suite, even though DB and SOAR are separate algorithms which use different bands and assumptions due to the differing characteristics of the aerosol retrieval problem over land and water surfaces.

This study uses the latest L1b data versions. For SeaWiFS this is obtained from the SeaWiFS Data Analysis System (SeaDAS) software package version 7.5 (available at https://seadas.gsfc.nasa.gov/, last access: 27 June 2019). SeaDAS applies vicarious calibration coefficients obtained as described in Franz et al. (2007) to SeaWiFS TOA reflectances. For MODIS and VIIRS the current L1b data versions are Collection 6.1 (C6.1) and version 2 respectively. The main difference between C6.1 and the previous Collection 6 (C6) L1b data is the development and implementation of a fix for crosstalk in MODIS Terra's tIR bands due to sensor degradation, which hindered the cloud mask in C6 in recent years (Moeller et al., 2017). For the VIIRS application of SOAR, and for the AAC retrievals discussed here, VIIRS TOA reflectances are cross calibrated to bring them into radiometric consistency with MODIS Aqua using the method and coefficients of Sayer et al. (2017); the residual uncertainty on this correction is approximately 0.5 %–1 % for the bands used here (not counting any error on the MODIS Aqua calibration itself).

DB/SOAR aerosol retrieval processing uses these L1b data at full resolution but provides output level 2 (L2, geophysical data) products at coarser resolution. These L2 aggregations are 3×3, 10×10, and 6×6 L1b pixels for SeaWiFS, MODIS, and VIIRS, respectively, giving L2 at-nadir horizontal pixel sizes of 13.5, 10, and 6 km respectively. To distinguish between native L1b pixels and the coarser L2 resolution, these latter sizes are often known as L2 “cells” rather than “pixels”. The bulk of retrieval uncertainty (for both total column AOD and AAC cases) is not due to radiometric noise but rather algorithmic assumptions; the coarsening has therefore been historically mostly to aid in pixel selection and post-retrieval quality filtering via analysis of L2 cell statistics (discussed later) and decrease the computational and data storage burden. This corresponds to one and two cells per scan line for MODIS and VIIRS, respectively; SeaWiFS imaged only one along-track pixel per scan. These output resolutions are also adopted here, due to the motivation for incorporation into the main DB data products.

2.2 Summary of the Sayer et al. (2016) AAC retrieval algorithm

The physical principle behind the demonstration AAC retrieval algorithm presented in Sayer et al. (2016) is that, in the presence of light-absorbing aerosols above a liquid-phase cloud, increases in COD brighten the TOA signal (as clouds tend to be bright and white) while AACs darken the signal as AOD increases. This darkening is more pronounced at the shorter wavelengths due to the tendency for increased absorption AOD (AAOD) at shorter wavelengths. For smoke aerosols this is due to the rapid increase in AOD with decreasing wavelength, while for dust it arises from the low single-scattering albedo (SSA; strong absorption) at blue and green wavelengths but SSA close to 1 at red and nIR wavelengths. Quantitative information about AACs can be extracted from the magnitude and spectral shape of TOA reflectance across this wavelength range (470–870 nm).

Sayer et al. (2016) retrieved AOD and COD at 550 nm (the state vector, x) simultaneously by a weighted multispectral least-squares fit of TOA reflectances in the four (blue, green, red, nIR) bands to modelled TOA reflectances stored in a lookup table (LUT). The RT calculations used to create the LUT were performed with a tool based on the Vectorised Linear Discrete Ordinates (VLIDORT) code (Spurr, 2006). The same VLIDORT-based tool is used in the present study. The solution is found by iterative minimisation of the squared residuals (differences between measured and LUT TOA reflectances) using the optimal estimation (OE) technique (Rodgers, 2000). OE propagates measurement and forward model uncertainties to provide an estimate S_x of the uncertainty on the retrieved state x,

where the covariance matrix S_y represents the uncertainty on the measurements (both radiometric and terms arising from forward model assumptions), and K (also known as the weighting function or Jacobian matrix) is the gradient of observations with respect to state measurements at the solution. More detail is given in Sect. 3.1 of Sayer et al. (2016). OE also provides a metric describing the level of agreement between measurement and modelled reflectances at the retrieval solution relative to the expected level of disagreement (retrieval cost, J), which is what is minimised iteratively in the retrieval:

This is the sum of square residuals between measured (y_m) and modelled LUT (y_LUT) reflectances relative to their expected magnitudes given in S_y, for a given point x in state space (i.e. combination of AOD and COD). Assuming S_y has realistic values and the measurements are informative on the state variables, J is expected to take values around the number of degrees of freedom (here two, as four measurements are being used to retrieve two parameters). These metrics are useful for quality assessment (QA) of the retrieval output; they are only quantitatively robust if the underlying forward model is appropriate, the input uncertainties well-quantified, and the forward model approximately linear near the solution (Povey and Grainger, 2015). OE can optionally also account for a priori information on the state vector, but that has not been included in the present implementation. LUTs are interpolated linearly during the retrieval, and K is calculated numerically. The first guess at x is taken as the LUT node point with the lowest cost, and convergence is typically obtained within three to four iterations.

Based on typical features of aerosol–cloud systems, instrument capabilities, sensitivity analyses, and retrieval simulations, the RT forward model is set up as follows (cf. Sect. 3.1 of Sayer et al., 2016, and references therein). The cloud is assumed to consist of a homogeneous and fully overcast layer with a top altitude of 1.5 km above surface level, geometric thickness of 0.3 km, and be composed of liquid water droplets with size conforming to a gamma distribution with an effective radius of 12 µm and effective variance of 0.1. The underlying surface is assumed to be Lambertian (see also Sect. 2.3.4 and 2.3.5). The aerosol is assumed to lie in a homogeneous layer with a top height 1 km above the cloud top and be 0.5 km thick.

Sayer et al. (2016) considered six different optical models for AACs, corresponding to four different types of smoke aerosols from different source regions, and dust aerosols with two different SSAs. These optical models were based on results from Aerosol Robotic Network (AERONET) almucantar scan inversions (Dubovik and King, 2000) representative of various source regions and aerosol types. Other sources of AACs such as volcanic ash were not included due to their comparative rarity and the fact that they have less repeatable and well-defined optical properties. Retrieval simulations in Sayer et al. (2016) indicated that the information content of the measurements was not always sufficient for the retrieval to select the correct aerosol type out of the six using the cost function alone. Therefore here the optical model representing strongly absorbing smoke derived from AERONET inversions in Mongu (Zambia) is used in all cases; based on previous studies this is expected to be reasonably representative of the smoke AACs encountered in the study region (Piketh et al., 1999; Eck et al., 2003, 2013; Swap et al., 2003; Reid et al., 2005; Queface et al., 2011). An exception is periodic additional contributions from mineral dust in the northern part from December to February (Pandithurai et al., 2001; Ben-Ami et al., 2009). Over the main ORACLES domain (green box in Fig. 2), and during the peak burning season, however, other AAC sources are expected to be negligible.

The Mongu smoke aerosol model was described by Sayer et al. (2014a). It is a bimodal log-normal optical model, such that fine (subscripted f throughout) and coarse (subscripted c throughout) mode volume size distributions are each of the form

for particles of size r, given total mode particle volume C_v, mode (which is also median and geometric mean) radius r_v, and width σ. Sayer et al. (2014a) found that the modal radius (in µm) of the fine mode was dependent on the fine-mode AOD at 550 nm (τ_f) as follows:

The spread σ (dimensionless) of the fine mode was found to have a weak dependence on (τ_f),

i.e. for higher smoke loadings the fine-mode particles were larger on average and had a broader distribution. In contrast, r_v,c and σ_c were found to be AOD-independent (across the small range of coarse-mode AOD observed) and take typical values of 3.34 and 0.67 µm respectively. Sayer et al. (2016) assumed a representative fine-mode fraction (FMF) of AOD at 550 nm, i.e. ${\mathit{\tau }}_{\mathrm{f}}/({\mathit{\tau }}_{\mathrm{f}}+{\mathit{\tau }}_{\mathrm{c}})$, of 0.9 for these smoke AACs, based on typical values from Sayer et al. (2014a). The assumed aerosol refractive index is shown for AERONET wavelengths in Table 1. These values are interpolated in log–log space to the satellite band centres for calculation of aerosol phase matrix elements and SSA. The resulting SSA (discussed further in Sect. 3.5.2) is weakly dependent on the AOD and in general varies from ∼0.86 in the blue band to ∼0.82 in the nIR band.

Table 1Spectral complex refractive index for the smoke aerosol optical model used in this study, following Sayer et al. (2014a).

Download Print Version | Download XLSX

Table 2 provides a brief summary of the main factors contributing to the retrieval error budget (in terms of effect on TOA reflectance) and when they are important. This is arranged in rough order of severity, based on results in Sayer et al. (2016) and discussion here, and with a focus on application to the ORACLES study region. Due to non-linearity of the retrieval system it is not trivial to map these into uncertainty on retrieved AOD/COD, as it is quite context-dependent; e.g. a large error on surface albedo would be important for AOD/COD retrieval for an optically thin cloud but is negligible for an opaque cloud with a COD of 10. As such Table 2 cannot be too specific as this would be misleading. The ability to assess sensitivities and provide uncertainty estimates on a case-by-case basis (Eq. 2) is an advantage of the OE retrieval framework applied here. The three leading factors added in quadrature provide the 3 % uncertainty on TOA reflectance assumed in the retrieval.

Table 2Magnitude of the main contributions to the total uncertainty on the TOA signal.

Download Print Version | Download XLSX

2.3 Algorithm updates since Sayer et al. (2016)

The same overall approach and RT forward model described in Sect. 2.2 is used in the present study, with updates described below. These are intended to improve upon approximations made in Sayer et al. (2016), in particular for retrievals for clouds above land surfaces, and prepare the AAC retrieval for integration with the standard DB and SOAR data products.

2.3.6 Retrieval QA

As in Sayer et al. (2016), QA metrics are used to filter the retrievals to remove scenes where the retrieval was not able to find a good fit between measured and modelled reflectances, or where unphysical spatial structure suggests that the forward model may have been inappropriate. These tests are similar to those described in Sayer et al. (2016), with updates based on examination of the larger data volume processed for this study. An example showing the overall QA flag and results for individual tests is given in Fig. 4. Pixels are only retained if the following criteria are all met.

• The retrieval cost (Eq. 3) is less than 5, indicating that the forward model is able to match the spectral TOA reflectance well. In practice cost function values tend to cluster in the range 0–2 or else be much higher than 5, so the results are only weakly sensitive to the value of this threshold.

• The COD ≥ 2, as for optically thin clouds the retrieval solution is often ambiguous and more sensitive to errors in surface reflectance assumptions. These factors do not always lead to a high value of the cost function. This is a slight relaxation of the COD threshold of 4 used in Sayer et al. (2016), due to the improved surface reflectance models used in this work. It can increase the potential data volume by 50 % or more in some cases, although some of these retrievals are subsequently removed by other QA tests.

• The retrieval has two or more (out of a possible eight) neighbours. Cases with zero or one neighbours are often found in conditions of broken cloudiness (e.g. cloud fragments in the middle of open-celled stratocumulus), which again may mean the forward model is not appropriate but does not always result in a high retrieval cost.

• The absolute difference between the retrieved AAC AOD and the median of AOD retrieved in the 3×3 retrieval box around it is smaller than 0.2. This removes spikes of high or low AOD which can result from isolated thin clouds, cloud mask errors, or poor surface assumptions. In practice, these retrievals are often around the edges of cloud fields. The physical basis behind this is that the AOD fields are expected to be spatially smooth on the scales of several retrievals. Note that the OE-provided uncertainty estimates (Eq. 2) for these retrievals are often (but not always) large (Fig. 4i). Sayer et al. (2016) implemented a test based instead on the estimated retrieval relative uncertainty, which had similar results for high-AOD artefacts but was less effective at identifying low-AOD outliers. This test might be less appropriate in other regions of the world where spatial variability in the aerosol field is higher.

Figure 4AAC retrieval for a MODIS Terra granule during SAFARI-2000 from 13 September 2000. Panels show (a) a true-colour image; (b) the overall QA flag; (c–f) results of individual QA tests, as described in the text; (g, h) the retrieved AOD before and after applying QA tests; and (i) the estimated uncertainty on retrieved AOD at 550 nm. The green box in (a) shows the region used for comparison with airborne data by Sayer et al. (2016). Pixels without valid retrievals are shaded in grey.

Download

The granule in Fig. 4 shows one of the test cases compared by Sayer et al. (2016) against airborne data in SAFARI-2000; in that case, the airborne measurements gave an estimate of the above-cloud AOD at 550 nm of 0.49, with a spatiotemporal standard deviation of 0.04, within the area outlined by a green box in Fig. 4a. The current version of the algorithm retrieves mean (median) AOD across this box of 0.48 (0.51), in very good agreement and close to the results of the demonstration algorithm shown in Sayer et al. (2016). The small difference from those prior results for this example is expected as the only relevant differences are the updates to the MODIS L1b data version (Collection 6 to 6.1) and aggregation/QA tests. Sayer et al. (2016) noted that the good agreement (AOD within 0.02) for this case may be fortuitous as the estimated uncertainty on the retrieved AOD (Eq. 2) is ±0.18, which is somewhat larger.

3.1 NASA Ames Spectrometers for Sky-Scanning, Sun-Tracking Atmospheric Research

Figure 5Flight tracks for the 2016 and 2017 ORACLES deployment, coloured by date. From top to bottom, panels indicate the 2016 P3, 2016 ER-2, and 2017 P3 aircraft flight tracks.

Download

The NASA Ames Spectrometers for Sky-Scanning, Sun-Tracking Atmospheric Research (4STAR) instrument is an aircraft-mountable hyperspectral Sun photometer and sky radiometer (Dunagan et al., 2013). It is a successor to the multichannel Ames Airborne Tracking Sunphotometer (AATS) instruments (Redemann et al., 2003; Schmid et al., 2003), which were used by Sayer et al. (2016) in validation of the initial version of the Deep Blue AAC retrieval algorithm (and cf. Fig. 4). 4STAR combines the Sun-tracking ability of AATS with a sky-scanning ability similar to that of ground-based AERONET Sun/sky photometers. Its full-width field of view (FOV) when measuring direct solar beam irradiance is 2.4^∘ (LeBlanc et al., 2019), with a radiometric deviation of less than 1 % in this span, compared to 3.7^∘ for AATS (Segal-Rozenhaimer et al., 2013). The smaller FOV reduces uncertainties due to scattered light in the direct-beam signal (Segal-Rozenhaimer et al., 2013; Smirnov et al., 2018). As operated during ORACLES, 4STAR has 1556 overlapping and continuous bands ranging from 350 to 1700 nm, compared to 6 or 14 distinct non-overlapping spectral bands on the AATS instruments (Dunagan et al., 2013).

The instrument was mounted on the NASA P3 aircraft for both the 2016 (based out of Walvis Bay, Namibia) and 2017 (based out of São Tomé) ORACLES deployments. Flight tracks where scientific data were collected are shown in Fig. 5. More information about the 2016 deployment, as well as 4STAR-derived aerosol data, can be found in LeBlanc et al. (2019). Flights included spiral profiles through smoke layers above clouds (as well as ramps and level legs), to enable the airborne instrumentation to measure atmospheric properties at different points within the smoke layers. The data set includes a flag (described in LeBlanc et al., 2019) to indicate measurements when the aircraft was flying above the top of a cloud and below a smoke layer. These flagged data points comprise a fairly small fraction of the total data set but allow estimates of the above AOD suitable for validation of the satellite retrievals.

The 4STAR data product used here is the spectral AOD from direct-Sun measurements. Exact measurement characteristics change between deployments, but in general the data are provided at around two dozen wavelengths (outside of strong gas absorption features) with a temporal resolution of 1 s. The uncertainty in spectral AOD is estimated on a point-by-point basis, and is largely driven by uncertainties on radiometric calibration and trace gas absorption, but is typically of order 0.005–0.02 (decreasing with increasing wavelength), much smaller than the expected uncertainty on the satellite retrieval. The current data versions used here are R3 for 2016 and R1 for 2017. Note that 4STAR measurements can also provide AERONET-like aerosol inversions (Pistone et al., 2019) and transmission-based cloud property retrievals (LeBlanc et al., 2015) which will be considered in a separate study.

3.2 NASA Langley High Spectral Resolution Lidar

The second airborne instrument used is the NASA Langley High Spectral Resolution Lidar version 2 (HSRL2; Hair et al., 2008; Burton et al., 2018). HSRL2 provides vertical profiles of atmospheric backscatter, depolarisation, and extinction; an advance of the 2016 deployment was the addition of a 355 nm channel (Burton et al., 2018) alongside the 532 and 1064 nm channels. Note that the 1064 nm channel lacks HSRL capability and is backscatter only, so above-cloud AOD is only provided at 355 and 532 nm. An advantage of the HSRL technique is that it is able to determine both backscatter and extinction, removing the need for a lidar ratio assumption, which can be a large source of uncertainty in backscatter lidar AOD retrievals such as from CALIOP (Omar et al., 2009).

During the 2016 ORACLES deployment (Burton et al., 2018), HSRL2 flew on the NASA ER-2 high-altitude aircraft (Fig. 5), also based out of Walvis Bay, Namibia. The ER-2 typically flew at an altitude around 20 km, well above the clouds and the bulk of the aerosols. As the ER-2 was flying at high altitude, a larger proportion of the flight provides data suitable for validating the AAC algorithm compared to 4STAR, for which only data collected immediately above a cloud top are relevant. During the 2017 ORACLES deployment, the HSRL2 instrument flew on the P3 with 4STAR at lower altitudes; this means that an unknown amount of aerosol above the plane will have been missed in the 2017 deployment. This should be borne in mind when examining the 2017 matchup statistics, along with the fact that in 2017 HSRL2 and 4STAR coverage is mutually exclusive. To decrease the contribution from this unobserved aerosol, 2017 HSRL2 data are only used here when the P3 was flying above 5 km (flight altitudes were generally below 8 km; when not spiraling, a reasonable number of legs were between 5 and 6 km, and few were above 6 km). The current data versions used here are R7 for 2016 and R1 for 2017.

Profiles are measured at 15 m vertical resolution and 10 s temporal resolution; the data contain a flag to identify profiles containing AAC cases. The spectral AOD was determined as described in by Hair et al. (2008), from the difference of the molecular channel signals at the top of the profile and at the top of the cloud. Assessments of the uncertainties of AOD determined from HSRL2 data are provided by Hair et al. (2008) and Burton et al. (2018). In brief, there is a random component (which is quantified within the data and typically negligible, ≪0.01) and a larger, locally systematic component. This systematic component is expected to be dominated by uncertainties in the molecular profile used in the retrieval and is difficult to quantify. As such, the uncertainty in HSRL AOD data is typically estimated by comparing against other simultaneous observations. Rogers et al. (2009) evaluated 532 nm AOD from an older version of the HSRL instrument in Mexico City and found a rms difference of 0.008 against AATS data and 0.056 against AERONET; in the latter case, some of the disagreement with AERONET was thought to result from small amounts of aerosol above the plane's flight altitude, and one outlying point (of only 10 total) contributed disproportionately to the higher rms difference. Sawamura et al. (2017) found a similar level of agreement against AERONET from HSRL-2 at both 355 and 532 nm from two field campaigns over urban areas in California and Texas.

3.3 Validation approach

Only retrievals passing the QA tests described in Sect. 2.3.6 are considered. As the airborne data have a higher spatial and temporal resolution than the satellite retrievals, the satellite data are validated by checking for and aggregating the 4STAR and HSRL2 data inside each individual retrieval footprint. Although this leads to a large number of matchups, it is important to bear in mind that the resulting matched data are not independent, due to the large autocorrelation in the underlying aerosol field, and retrieval errors are similarly likely to be autocorrelated on these length scales. The airborne data are available only for a limited spatial domain over a short time period within each year. This is a different picture from total column AOD validation using ground-based AERONET sites, which are composed of individual dispersed sites as opposed to flight tracks. For this reason, as well as statistics for all matchups individually, granule-average statistics (i.e. statistics calculated using averages of all matchups from individual granules) are also presented. These should exhibit reduced autocorrelation compared to the all-matchups data. Note that these are calculated averaging all matched retrievals and airborne data within individual granules, not simply averaging all retrievals within the granules.

The satellite overpasses and flight tracks were mostly not simultaneous, and a time difference threshold of ±3 h is used as a cut-off for a matchup to be valid. This is longer than the ±0.5–1 h often used for comparison against AERONET sites and is adopted as the temporal variability of these large-scale smoke plumes is expected to be somewhat limited. Part of the rationale for a shorter time threshold in AERONET validation analyses is the potential for an incoming cloud field to remove or modify the aerosol during the time between measurement and overpass; as the AAC retrieval is concerned with those aerosols above (and less likely to be modified by) clouds, that rationale is less relevant here. Using a stricter time threshold in this analysis essentially has the effect of removing individual flight legs from consideration; due to the limited number of flights available (Fig. 5), it is difficult to disentangle contributions from true temporal variability from those due to individual flight characteristics (i.e. sampling differences) in the changes in comparative statistics, although the overall picture does not significantly change with a threshold of ±1 h (not shown). Jethva et al. (2018) compared an OMI-based algorithm with HSRL2 measurements from ORACLES 2016; they found (their Fig. 4) that imposing a time difference threshold of ±1–2 h improved some comparison statistics, compared to no time difference threshold. That appears to be driven in part by the removal of some high-AOD events when either time difference threshold was imposed. This suggests some aerosol motion over the course of a day but less over the course of several hours, so it is not inconsistent with the results here. The total number of matchups (and individual granules containing matchups) is shown in Table 3.

Table 3Number of individual retrieval matchups (and number of contributing granules, in parentheses) for each satellite sensor and ORACLES data set.

Download Print Version | Download XLSX

The AOD is evaluated at the satellite wavelengths used (i.e. bands centred near 470/490, 550, 650/670, and 865 nm, dependent on sensor), as well as 500 nm, as the latter is (along with 550 nm) a commonly used reference wavelength in aerosol analyses. For the HSRL2 data the available AOD at 355 and 532 nm are interpolated to 470/490 and 500 nm, and extrapolated to 550 nm, using the Ångström exponent (AE, denoted α) where

over the wavelength range λ₁–λ₂ (here 355–532 nm). For 4STAR, up to 12 AOD measurements are available across the relevant wavelength range. Therefore, following Eck et al. (1999) a least-squares fit of all available AODs to a quadratic polynomial is performed and used to calculate the AOD at each wavelength of interest:

Coefficients a₀, a₁, and a₂ are calculated on a point-by-point basis. This quadratic formulation is more robust to calibration problems in individual channels and accounts for the fact that in fine-mode dominant aerosol conditions the relationship between log (τ) and log (λ) is not linear but curved, depending on fine-mode particle size (Eck et al., 1999; Schuster et al., 2006). The longer wavelengths are not considered for the HSRL2 comparison to avoid the potentially larger extrapolation errors due to this spectral curvature; likewise, the availability of only two wavelengths means that Eq. (11) cannot be applied for HSRL2.

For 4STAR, the uncertainty on an individual matchup is taken as the median of the uncertainties on the spectral AOD used for the fitting in Eq. (11) (and is typically around ±0.01). For HSRL2, the uncertainty is taken as ±0.02 at 470/490 and 500 nm and ±0.03 at 550 nm, to allow for a small extrapolation error. In both cases, the standard deviation of measurements within each satellite footprint is added to this in quadrature to account for potential spatiotemporal heterogeneity. This latter term is typically 0.01 or smaller, and the total uncertainty is likewise typically much smaller than the estimated uncertainty on the satellite retrievals.

Due to the different flight locations (Fig. 5) and potential for different systematic uncertainties in the airborne data between deployments, results are reported separately for 2016 and 2017. The main metrics used here to evaluate the AAC retrievals, which are as often used in AOD validation exercises, including DB (e.g. Sayer et al., 2018b, 2019), are as follows.

1. The correlation coefficient (R), as a measure of how well the satellite data track the variability of the airborne data. Spearman's rank correlation coefficient is used rather than the more common Pearson linear correlation coefficient. The reasons for this include the facts that the relationship between airborne and satellite AOD may not be linear, and also that Spearman's correlation is less sensitive to extreme outliers (either sampling-related or retrieval problems) which may be unrepresentative of the behaviour of the data set. While Pearson correlation has historically been the more frequently used one in aerosol data analyses, other fields are increasingly appreciating the use of Spearman correlation for situations where this is better supported by the nature of the data (e.g. the medical literature, Schober et al., 2018).

2. The median bias between the data sets, defined satellite–airborne, as a measure of the general offset. Again, medians are more robust to outliers which can skew the means.

3. The median relative bias between the data sets, defined (as above) relative to the airborne data.

4. The root mean square error (RMSE), which is a commonly reported metric, although is dependent upon the typical level of AOD as well as the presence of outliers.

5. The mean absolute error (MAE), similar to RMSE but less weighted by outliers.

6. The fraction (f) of points matching within the total expected level of difference (ED). The ED is taken as the quadrature sum of the expected retrieval uncertainty σ_ret (square root of the relevant element of S_x in Eq. 2) and aforementioned airborne uncertainty/variability σ_air under the assumption that these two are independent, i.e. $\mathrm{ED}=\sqrt{{\mathit{\sigma }}_{\mathrm{ret}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{air}}^{\mathrm{2}}}$. For satellite-retrieved and airborne AOD, τ_ret and τ_air respectively, the relevant inequality assessed is therefore the fraction satisfying $|{\mathit{\tau }}_{\mathrm{ret}}-{\mathit{\tau }}_{\mathrm{air}}|\le \mathrm{ED}$. If these uncertainties are appropriate, then one standard deviation (∼68 %) of matchups should be in agreement within this bound and two standard deviations (∼95 %) within twice this bound. Again, however, the spatiotemporal autocorrelation in the observations and limited sample size mean that this is unlikely to be true for this particular set of data. Still, the metric provides a general guideline on how quantitatively similar the estimated uncertainties are to the actual retrieval errors. This statistic is not presented for the granule-average comparison, because it is not meaningful for that case.

3.4 Validation results

Figure 6Scatter plots and summary statistics for the comparison between 550 nm AOD from MODIS Aqua AAC retrievals and 4STAR data, during ORACLES 2016. Statistics are as defined in the text. Panel (a) shows the comparison for all individual matchups; horizontal and vertical error bars indicate the estimated uncertainties on the 4STAR and satellite retrievals, respectively. Panel (b) shows median 4STAR and MODIS data from matchups obtained within each granule, and horizontal and vertical error bars show the standard deviation of matched 4STAR and satellite AOD within each granule, respectively. The 1:1 line is dotted grey.

Download

Figure 6 shows one example of instantaneous and granule-averaged results, for the case of MODIS Aqua and 4STAR data in 2016. Here, 15 granules contributed a total of 432 matchups. The bulk of the points in both cases cluster around the 1:1 line. For the instantaneous matchups, there are some outliers, which tend to be retrieved with a large estimated uncertainty; in general, the estimated uncertainty on the satellite retrievals is, as expected, larger than that due to uncertainty and variability in the airborne data. A lot of the scatter is decreased when going to granule-averaged statistics, such that correlation increases and MAE and RMSE decrease. The absolute bias does not change much. Interestingly, the variability on the granule-averaged satellite data (vertical bars in Fig. 6b) tends to be somewhat smaller than the uncertainty on the individual matchups (vertical bars in Fig. 6a). This is likely due to high autocorrelation in the retrieval uncertainties (i.e. an error source on a given retrieval is likely to be very similar to the errors on retrievals adjacent to it), which is a result of the flight-track sampling of airborne data. This also indicates that, as with many other AOD retrieval algorithms, the bulk of the error is not true random noise but rather locally systematic due to the context (i.e. geometry, atmospheric, and surface conditions) of the retrieval. Similar patterns (not shown) are observed for the other satellite sensors and airborne deployments.

Figure 7Summary line plots of spectral AOD validation statistics. Columns show (left to right) comparisons for 4STAR 2016, 4STAR 2017, HSRL2 2016, and HSRL2 2017. Rows show (top to bottom) rank correlation, median (satellite–airborne) bias, median relative bias, RMSE, MAE, and fraction f agreeing within the ED. In all panels, solid lines denote statistics for all matchups and dashed for granule-average comparisons. Data for MODIS Terra, MODIS Aqua, and VIIRS are shown in red triangles, green diamonds, and blue squares respectively.

Download

In Fig. 7, summary statistics equivalent to those presented in Fig. 6, but for all wavelengths and satellite/airborne comparisons assessed, are presented. Several statistics (e.g. correlation, f) show limited spectral dependence. Others (e.g. RMSE, MAE, and in some cases the bias) shrink with increasing wavelength, which is expected due to the rapid decrease in AOD of smoke with increasing wavelength. Results for the granule-averaged comparison are often similar to those from the instantaneous comparison (sometimes slightly better, sometimes slightly worse); the same tendencies are seen between satellite sensors and across wavelengths. This also points to the bulk of the errors in the retrieval being contextual rather than truly random (aside from a few individual outlying pixels). The HSRL2 comparison shows a smaller difference between instantaneous and granule-averaged comparison statistics than 4STAR, perhaps due to the generally larger number of matchups, but a smaller number of contributing granules for HSRL2.

Interestingly, the different ORACLES comparison data sets reveal some different patterns. For example, the 2016 data (both 4STAR and HSRL2) indicate near-zero (MODIS) and negative (VIIRS) bias tendencies in the satellite data, while for the 2017 data the biases tend to be more positive. In this sense, the different deployments do not paint identical pictures about the retrieval error characteristics. Recalling the facts that in 2016 4STAR and HSRL2 were on separate aircraft but flying in similar locations and at similar times, while in 2017 they were on the same aircraft and flying over a different region (Fig. 5), this suggests that data sets from single deployments may not be providing sufficient sampling of the aerosol–cloud system to fully characterise satellite retrieval uncertainties. The differences might be partially coincidental due to the particular cases sampled on the flights or may reflect more persistent differences in the locations of the two deployments; it is difficult to disentangle these two possibilities with the available data. It is therefore cautioned that the validation results presented here may not have sufficient sampling to be comprehensive, and further field campaigns in this region (and others) would be desirable to obtain a fuller validation of AAC retrievals. Note that the 2018 ORACLES flight tracks followed a similar pattern to those in 2017; the 2018 data are not publicly available at the time of writing and an initial release is expected later in 2019.

The similarity between MAE and RMSE lines in Fig. 7 indicates that there are few extreme outliers, as RMSE is sensitive to outliers while MAE is more robust. This is encouraging and provides further evidence that the QA tests (Sect. 2.3.6) are reasonably successful at removing cases where the forward model is inappropriate. The fraction f of matchups agreeing within ED is similar to the theoretical value of 68 %, indicating that on average the estimated uncertainties provided by the OE technique (Eq. 2) and uncertainty characterisation of the airborne data are reasonable. The spectral stability of f (as well as the other statistics) is further evidence that the uncertainty characterisation and retrieval assumptions (Table 2) are reasonable.

Figure 8Histograms of normalised retrieval error (i.e. actual error divided by expected difference ED) for AOD at 550 nm. Panels show (left–right) data for MODIS Terra, MODIS Aqua, and VIIRS matchups. In all cases matchups from 4STAR 2016, 4STAR 2017, HSRL2 2016, and HSRL2 2017 are shown in red, green, blue, and purple respectively. The black line shows the theoretical Gaussian distribution with mean 0 and variance 1, and dotted and dashed lines indicate ±1 and ±3 standard deviations, respectively.

Download

Figure 9Comparison between magnitudes of expected difference (ED) and actual absolute retrieval errors. The top row shows ED (i.e. 1σ uncertainty) against the 68th percentile (i.e. 1σ) retrieval error, binned as a function of ED. The bottom row shows 2 × ED (i.e. 2σ uncertainty) against the 95th percentile (i.e. 2σ) retrieval error, for the same bins. Panels show (left–right) data for MODIS Terra, MODIS Aqua, and VIIRS matchups. Colours are as in Fig. 8. The 1:1 line is dotted grey.

Download

As noted, theoretically the ED should indicate the one standard deviation (1σ, ∼68th percentile) expectation of disagreement between satellite and airborne data. Collectively, the distribution of normalised retrieval error $({\mathit{\tau }}_{\mathrm{ret}}-{\mathit{\tau }}_{\mathrm{air}})/\mathrm{ED}$ should approximate a Gaussian distribution with mean 0 and variance 1. A normalised error of +1 means that the retrieved AOD was 1 × ED higher than the airborne AOD for a particular matchup, for example. This distribution is assessed for the 550 nm data in Fig. 8. The distributions appear reasonable, although they tend to peak too strongly near a normalised error of 0 and (particularly for VIIRS) have more negative outliers than expected. Differences between the statistics for the different ORACLES deployments are again also visible. Figure 9 examines this another way, comparing the actual and expected retrieval errors as a function of ED (in 10 equally populated bins, in each case). Here, the top row compares actual vs. expected 1σ errors (i.e. 68th percentile of absolute retrieval error in each bin) and the bottom row the same for 2σ (i.e. 95th percentile) errors. For a perfectly characterised retrieval system, these points should lie on the 1:1 line. They share a common tendency for underestimating the retrieval error when the ED is low and overestimating when it is high, with the crossover point being an ED around 0.15–0.2. This latter point (i.e. if a large ED is estimated, it tends to be too large) was also found in the retrieval simulations performed in Sayer et al. (2016). This may be due to non-linearity in the retrieval system in these conditions, in which case the validity of the OE uncertainty estimate is expected to break down.

Figure 10Quantile–quantile (QQ) plots comparing distributions of AODs from co-located satellite and airborne measurements, from 5th to 95th percentiles of the matched data. Panels show (left–right) data for MODIS Terra, MODIS Aqua, and VIIRS matchups. Colours are as in Fig. 8. The 1:1 line is dotted grey.

Download

The opposite case (i.e. if a very small ED is estimated, it tends to be too small) most commonly occurs when the satellite-retrieved AAC AOD is near zero but the airborne data report an AOD around 0.1–0.15. This suggests that the error budget is missing some component which can be important in fairly low-AOD conditions, perhaps related to calibration uncertainty, the cloud model, or some correlation between forward model error at different wavelengths. From Fig. 8, these large negative outliers tend to occur more frequently in VIIRS than in MODIS. This is further supported by quantile–quantile (QQ) plots of the matched data, shown in Fig. 10. The QQ plots reveal that for MODIS Terra/Aqua the distributions of satellite and airborne AOD are fairly similar (although satellite AOD are often slightly higher). In contrast for VIIRS it is common for the retrieval to report near-zero AOD a disproportionately high fraction of the time. The reasons for this are not yet known; it is plausible that they are related to limitations of the current cloud mask used (Sect. 2.3.1). VIIRS also has a broader swath than MODIS, although retrieval errors as a function of viewing and scattering angles were examined for all sensors and no patterns could be found with the available sampling (not shown). Since the sensitivity to AOD comes largely from the magnitude of spectral darkening across the visible wavelength range, it is also possible that a small calibration of forward model bias is responsible. Overall, these results indicate that the MODIS-derived AAC record is presently likely to be more reliable than the VIIRS-derived AAC record. Note that in Fig. 10 the lines belonging to data for the same year are more similar to each other than the lines for the same instruments (i.e. 4STAR or HSRL2) for different years, further implying that apparent differences in performance are likely related to the specific scenes observed each year.

Nevertheless, the bottom row of Fig. 9 shows that the tails of the uncertainty distribution (2σ errors) tend to be quantitatively better estimated than the (1σ errors). This indicates that the current uncertainty estimates do have some quantitative value for identifying retrievals with larger errors. The combination of occasional large positive and negative outliers and the fact that the ED is somewhat linked to the retrieved AOD (low-AOD cases tend to have a low ED, high-AOD cases a higher ED) suggests that for calculating daily level 3 aggregate data, medians may be a better option than either simple means or error-weighted means. This is because the AOD fields tend to be fairly spatially coherent, while either a simple or weighted mean may bias the aggregate.

3.5 Evaluation of retrieval assumptions

3.5.1 Spectral dependence of AOD

Figure 11AE (470–870 nm) assumed in the retrieval as a function of the AOD at 550 nm (black), together with airborne 4STAR data from the 2016 (red) and 2017 (blue) ORACLES deployments, as well as the estimated smoke component (see text) of MAN cruises in the region (green). For the 4STAR and MAN data, points and lines indicate bin medians and the central 68 % of data, respectively.

Download

The AOD dependence of the size distribution in the aerosol optical model assumed in the retrieval (Sect. 2.2) results in the wavelength dependence of AOD being a function of aerosol loading. This dependence, illustrated as the AE over the wavelength range 470–870 nm (Eq. 10), is shown in Fig. 11. The decline from values near 2 in low-AOD conditions to ∼1.7 in high-AOD conditions is a result of the AOD dependence of Eqs. (5) and (6).

These data are compared with two other sources; the first is the AE calculated over the same wavelength range (from all available spectral AODs) from the 2016 and 2017 4STAR deployments. These data are then divided into 25 evenly populated bins as a function of the AOD at 550 nm. The second is four Maritime Aerosol Network (MAN, Smirnov et al., 2009) cruises (Saint Helena 2016 and Research Vessel Meteor 2013, 2015, and 2016) within the region in recent years. These cruises included measurements near the coast, in the open ocean, and under varying levels of smoke influence. MAN AOD data consist of measurements with a handheld Microtops Sun photometer, with a typical level of uncertainty around ±0.02 (Knobelspiesse et al., 2004). As they measure total column AOD, they include the contribution from maritime aerosol as well as smoke layers. Here, the maritime contribution is taken as the 20th percentile of all AODs measured in these cruises (0.168, 0.148, 0.112, and 0.093 at 440, 500, 675, and 870 nm respectively), which is then subtracted from the total MAN AOD. The remaining AOD is assumed to be smoke, from which the 550 nm AOD and AE are calculated and the same binning exercise carried out (using only five bins due to the small data volume). Using the 20th percentile as a threshold is somewhat arbitrary, although the resulting mid-visible AOD is similar to other estimates, and reasonable changes to this threshold result in changes of order ±0.1 to bin-median AE. Failing to remove a maritime component would give a misleading AOD dependence, as low-AOD conditions would be dominated by the background (low-AE) sea spray aerosol while high-AOD conditions would be more dominated by smoke (but still influenced by the maritime contribution).

All these data exclude points with AOD at 550 nm below 0.1, as AE calculation is highly uncertain when the AOD is low (e.g. Wagner and Silva, 2008), leaving 34 397, 14 988, and 220 4STAR 2016, 2017, and MAN data points respectively. The variability within these binned data arises from both noise in the AE calculation and real spatiotemporal variability. Note that some of the low-AOD 4STAR data from 2017 are thought to have sampled a small amount of dust, which would explain the anomalously low AE in the lowest bins. HSRL2 data are not used here because its wavelength range (355–532 nm) would be expected to result in AE values lower by ∼0.5 compared to the 440–870 nm spectral range (Eck et al., 1999; Schuster et al., 2006), which would be less directly comparable.

Figure 11 shows that the optical model assumed in the retrieval lies on the upper end of the 4STAR and MAN observations. There is a closer match with MAN values, although these are more uncertain than 4STAR due to the subtraction of the estimated maritime contribution. The offset from bin-median 4STAR AE values over much of the AOD range is ∼0.2. The practical implications of an AE overestimate of 0.2 translate to an approximate 3 % overestimate and 10 % underestimate of AOD at the most extreme wavelengths of 470 and 870 nm, respectively, somewhat smaller than the total estimate of retrieval uncertainty in most cases. Therefore the spectral dependence of AOD in the aerosol optical model assumed in the retrieval seems reasonable. Using 4STAR data from the 2016 deployment, LeBlanc et al. (2019) observed a general tendency for increasing AE with altitude for AAC cases (ranging from ∼1.5 at 0.5 km to ∼2.0 near 4 km). The data collected were most dense from around 1–2 km, over which the AE was fairly flat around 1.7; altitudes below 1 km or above 3 km were comparatively poorly sampled and so possibly less representative. This structure is also a secondary contribution to the retrieval uncertainty (Table 2). The uncertainty and variation within the 4STAR and MAN data are insufficient to determine whether the small decrease in AE with increasing AOD in the assumed optical model is reproduced by these direct-Sun measurements. The 2017 4STAR data do show this decline, although as this draws from a small number of flights it may not be representative.

3.5.2 Aerosol SSA

Figure 12Relevant spectral SSA data collected in (a) August, (b) September, and (c) October. The grey shaded region indicates the range of the assumed SSA for mid-visible AOD between 0.1 and 1. Mean and standard deviation of AERONET inversions at sites which tend to sample near-source and transported smoke are shown in red and blue triangles respectively. Monthly mean and standard deviation of surface measurements during the LASIC field campaign, reported by Zuidema et al. (2018), are in black diamonds. Airborne in situ measurements from SAFARI-2000 reported by Haywood et al. (2003) are shown for fresh and aged smoke in light and dark green boxes, respectively; airborne remotely sensed measurements from SAFARI-2000 reported by Bergstrom et al. (2003) are shown in purple.

Download

Similarly to the AE case, the SSA assumed in the AAC retrieval also varies with AOD due to the changing particle size (Eqs. 5 and 6). This AOD dependence is fairly small: a range of ∼0.02 at 470 nm and ∼0.04 at 870 nm. This is shown in Fig. 12, together with SSA estimated from a range of other sources within this region. These are split by month of year, as previous work (e.g. Eck et al., 2013; Zuidema et al., 2018) has reported a gradual increase in SSA through the burning season. SSA variations are linked to differences in fuel types, ageing, and a transition between flaming and smoldering combustion through the burning season resulting in part from meteorological patterns (Zheng et al., 2018).

The first are AERONET inversions from nine sites; four of these sites tend to sample near-source burning (and see Fig. 2) while the other five tend to sample transported smoke. Further background information and locations for these sites can be found at https://aeronet.gsfc.nasa.gov (last access: 27 June 2019). These classifications are somewhat fuzzy; sites can sample smoke from a variety of ages, due to meteorological patterns which lead to recirculation of air masses over the continent (e.g. Swap et al., 1996; Tyson, 1997), and the aerosol loading (particularly in the south-eastern part of the region) can contain additional sulfate contributions from power plants (Piketh et al., 1999; Eck et al., 2013). Still, the classification and set of sites provide an indication of the regional variation of SSA encountered.

The AERONET data shown are monthly means and standard deviations of version 3 level 2.0 inversions (Giles et al., 2019); these are pre-filtered to remove poor-quality retrievals, as well as retrievals with an AOD at 440 nm < 0.4, for which the SSA is quantitatively less reliable. Data were further filtered to remove points with AE < 1.2, to restrict to smoke-dominated cases, although this had a negligible effect on the results. The uncertainty on the level 2.0 SSA is expected to be ±0.03 (Dubovik et al., 2000). The bulk of this is due to calibration uncertainty and is therefore systematic within a given (roughly year-long) deployment (Dubovik et al., 2000; Eck et al., 2013); most of these sites are multi-year records, such that these uncertainties may partially cancel out. Comparisons with the previous version 2 AERONET data (not shown) reveal quantitatively similar climatological results for this region.

Also shown are the monthly mean and standard deviation of surface-based estimates made at Ascension Island from the LASIC field campaign in 2016, reported by Zuidema et al. (2018). These estimated the SSA at 529 nm from nephelometer measurements of aerosol scattering and particle soot absorption photometer (PSAP) measurements of absorption. Zuidema et al. (2018) cautioned that these near-surface measurements (sampling air masses from the boundary layer) may not always be representative of the total column and noted that Leahy et al. (2007) found that airborne in situ observations from SAFARI-2000 tended to report lower SSA than total column estimates.

Next, SSA at 450, 550, and 700 nm for two SAFARI-2000 flights from Table 2 of Haywood et al. (2003) are shown. These correspond to flights sampling fresh smoke (the Otavi plume, flight a790, 13 September) and aged smoke at Ascension Island (flight a794, 19 September). These were computed from size distributions measured by a passive cavity aerosol spectrometer probe (PCASP) with an assumed refractive index, resulting in an estimated uncertainty of ±0.04. The Otavi case (but not the Ascension Island flight) also included nephelometer/PSAP measurements which gave a very similar SSA; Haywood et al. (2003) noted that refractive index assumptions used in the PCASP calculation were informed by the PSAP data. More recent work by Ogren (2010) revealed an error in the corrections applied in PSAP data processing which result in a typical underestimate of reported absorption by ∼13 %; for the Haywood et al. (2003) cases, this translates to an overestimate of SSA by ∼0.02, within the notional uncertainty but systematic. It is plausible that this influences the PCASP results as they were informed by PSAP data. Newer PSAP analyses such as Zuidema et al. (2018) shown here apply corrections accounting for this error.

SAFARI-2000 results for flight a786 near Mongu reported by Bergstrom et al. (2003), calculated from Solar Spectral Flux Radiometer (SSFR) data, are also shown. The SSFR technique uses measurements of flux above and below the smoke layer, constrained by meteorological data and AATS AOD. The SSA uncertainty on individual wavelengths was ∼0.02 for this case; the best-fit curve is plotted in the figure. Haywood et al. (2003) also presented results for that flight, and the Bergstrom et al. (2003) data lie in between the PCASP and nephelometer/PSAP results for that case (not shown); it is not clear whether data from the same parts of this flight were used for those cases, and the age of the smoke sampled is uncertain.

Consistent with previous studies (e.g. Reid et al., 2005; Eck et al., 2013, and references therein) there is an increase in SSA from August to October (by ∼0.05), and data sampled near source regions tends to be more strongly absorbing than data from aged air masses. The AERONET site at Windpoort is the most strongly absorbing of those placed into the predominantly transported sites, perhaps because it is closer to the source regions than the others. AERONET monthly mean results from Mongu and Mongu Inn (the latter being a replacement site for the former, in the same location) are offset by up to 0.04. The range of interannual variability in SSA at Mongu (not shown) is typically ∼0.05, although some of that would be expected to average out over the time period available (10 years for Mongu, 4 for Mongu Inn). It is therefore possible that the two sites have slightly different error characteristics. The LASIC data are somewhat lower than the others, possibly due to the aforementioned surface sampling.

Within a given month, the AERONET and SAFARI-2000 results tend to span a range of ∼0.06 in the blue spectral region and a larger range of up to ∼0.1 at red and nIR wavelengths, with variability decreasing through the burning season. Recent SSA observations from the ORACLES 2016 deployment using remote-sensing and in situ instrumentation, presented by Pistone et al. (2019), were found to agree within the range bounded by the previous SAFARI-2000 values, though overall the ORACLES SSAs tend towards the lower range (instrument medians between 0.85 and 0.88 at 500 nm). As with the results in Fig. 12, the different techniques did show some variation, due to a combination of real variability and retrieval/measurement uncertainty; SSA is a difficult quantity to measure. The AAC retrieval's SSA assumption (which is centred near 0.875, 0.87, 0.86, and 0.85 at 470, 550, 650, and 870 nm respectively) is in the middle of this range, although the large variability in the reference data suggests large spatiotemporal variability in aerosol optical properties. This implies the potential for spatial/temporal structure in the AAC retrieval error which may explain some of the differences in bias between 2016 (mostly September) and 2017 (mostly August) flights in Sect. 3.4. However, Sayer et al. (2016) found that the AAC algorithm was less sensitive to errors in SSA assumptions in cases of strongly absorbing aerosols (such as here) than weakly absorbing aerosols. This is because the sensitivity (darkening) of TOA reflectance to changes in AOD is stronger for more strongly absorbing aerosols.

A summary of these results is that SSA is variable in space and time in this region, but the retrieval assumptions are broadly in-family and likely a reasonable approximation for typical smoke conditions both near source and downwind. Analysis of in situ and remotely sensed SSA data from CLARIFY is ongoing (Kate Szpek, personal communication, 2018). When complete, the airborne results from both campaigns can be used to inform possible updates to the optical model, such as considering a seasonally and/or longitudinally dependent SSA. This must be balanced against the danger of over-tuning results to a limited data set, which is why the original AERONET-based optical model (Sayer et al., 2014a, 2016) is retained for the present work.

4.1 Time series

Figure 13Monthly time series of various satellite data sets over the green box (25^∘ S–0^∘ N, 15^∘ W–15^∘ E) in Fig. 2. Panel (a) shows the monthly mean UVAI from the multi-sensor (MS) and OMI data sets, and (b) shows the fraction of days in each data set where the box-averaged UVAI is over 0.75. Panel (c) shows the mean AAC 550 nm AOD using the algorithm presented in this work, applied to SeaWiFS, MODIS Terra/Aqua, and VIIRS measurements. Panel (d) is a time series of monthly mean total column (cloud-free) over-water 550 nm AOD from the MODIS Terra/Aqua DT data sets. Panel (e) is the difference between total column and above-cloud AOD (i.e. d–c), estimated for MODIS Terra and Aqua. Panel (f) shows monthly total corrected fire counts from MODIS Terra and Aqua (box shifted 25^∘ E from the others). Throughout, MS data are shown in magenta, OMI in purple, SeaWiFS in teal, MODIS Terra in red, MODIS Aqua in green, and VIIRS in blue. Months with fewer than three contributing days are excluded.

Download

This section briefly examines spatiotemporal patterns in the 20-year record obtained by applying the AAC retrieval algorithm presented here to the four satellite sensors. A broader study comparing results against other satellite AAC AOD/COD data products is planned for the future. First, monthly time series of the retrievals and other relevant satellite data sets are shown in Fig. 13. These are constructed by averaging daily data over the green box (25^∘ S–0^∘ N, 15^∘ W–15^∘ E) in Fig. 2, which corresponds to the core of the stratocumulus cloud deck and the main flight region for ORACLES deployments, and then computing monthly averages from these.

Figure 13a shows time series of UVAI from two data records: the multisensor (MS) UVAI data set version 1.7 combines TOMS, GOME, SCIAMACHY, and OMI observations, dividing each of the fairly coarse resolution sensor pixels into several subpixels to produce a long-term (starting 1978) data set with consistent spatial resolution (Tilstra et al., 2012, 2013). Also shown are the latest version 1.8.9.1 OMI UVAI data, described by Torres et al. (2018). This latest OMI data version updates the UVAI calculation to decrease variations associated with changes in solar/sensor geometry and particle shape which influence apparent seasonality. Figure 13b provides the fraction of days within each month where the box-average UVAI was above 0.75, a subjective but reasonable threshold (e.g. Tilstra et al., 2013; Torres et al., 2018) for the presence of non-negligible levels of absorbing aerosols. Figure 13c and d present the AAC AOD from the algorithm presented here and the total column (i.e. from cloud-free scenes) AOD from the MODIS Dark Target (DT) Collection 6.1 over-water algorithm (Levy et al., 2013). Figure 13e provides an estimate of the below-cloud AOD from MODIS, by subtracting the AAC data in Fig. 13c from the total AOD in Fig. 13d. Strong caution is required in this as the DT and AAC algorithms are independent and have different error characteristics, although it provides a crude proxy for the relative partitioning of aerosol above and below cloud level. Finally, Fig. 13f shows fire counts from the cloud-corrected overpass-corrected MODIS current Collection 5 (Giglio et al., 2003, 2006) data set (MOD14CM1/MYD14CM1). The fire counts represent total detections rather than an average, and the longitude range is shifted to cover the source region 10–40^∘ E (as no fires occur over ocean). The fire data are also not presently available for the full MODIS records.

Table 4Correlation coefficients between monthly mean AAC 550 nm AOD and other time series shown in Fig. 13.

Download Print Version | Download XLSX

In all of these time series (aside from the below-cloud AOD estimates in Fig. 13e), the annual cycle of fires and associated emissions, strongest from June–September, is evident. Interannual variability is comparatively limited but generally consistent between data sets. A smaller secondary peak from December to February is also seen, likely due to a combination of Sahelian fires and dust transport (Pandithurai et al., 2001; Ben-Ami et al., 2009). MODIS Terra fire counts are around a factor of 5 lower than those observed by Aqua; this pattern was observed in multiple global source regions by Giglio et al. (2006) and ascribed to diurnal variations in fire activity. Table 4 shows the correlation between each of the time series in Fig. 13 and the four AAC AOD data sets. Correlation coefficients are high (0.78–0.94) and show small variability between the four AAC AOD records generated here. Due to the small number of points in the time series, the differences in correlation coefficients between sensors are not statistically significant and it is not possible to state robustly which quantity is most strongly correlated with the retrieved AAC AOD. These results indicate that these quantities may provide a useful proxy for variations in aerosols transported above clouds, if AAC retrievals are not available. However the strength of the relationships might not hold for other regions where aerosol and cloud properties covary differently.

The exception is the estimated below-cloud AOD, which is only very weakly correlated with the above-cloud AOD. This might imply that very little of the smoke is transported within the marine boundary layer, which is generally consistent with the discussion in Sect. 3.5.3. However, as mentioned previously, due to large uncertainties caution should be used in interpreting these data. The mean and standard deviation of below-cloud 550 nm AOD estimated from MODIS Terra and Aqua are 0.15±0.05 and 0.13±0.04 respectively, which is only slightly larger than ship-based measurements of AOD in maritime environments without significant continental influence (Smirnov et al., 2009, 2011).

4.2 Spatial patterns and offsets

Table 5Comparative statistics for monthly mean AAC 550 nm AOD between the four AAC data sets generated in this work. Offsets are defined subtracting the second indicated sensor from the first.

Download Print Version | Download XLSX

Figure 14Histograms of AAC 550 nm AOD retrieved over the green box (25^∘ S–0^∘ N, 15^∘ W–15^∘ E) in Fig. 2, aggregated from level 2 retrievals over the full satellite records processed during the ORACLES campaign season. Panels show data for August, September, and October. Throughout, SeaWiFS data are shown in teal, MODIS Terra in red, MODIS Aqua in green, and VIIRS in blue.

Download

Figure 13 also shows offsets between the AAC retrievals, with SeaWiFS being the highest and VIIRS the lowest. This is consistent with the validation results in Sect. 3 (aside from the SeaWiFS mission which ended in 2010 and so cannot be directly validated with ORACLES data). Table 5 quantifies the consistency between these time series, revealing very high correlation coefficients (0.95–0.99) despite these offsets. Pearson correlation coefficients are calculated for this instance, as the data sets are notionally inferring the same quantity using a similar technique and should be subject to the same causes for outliers (e.g. extreme events in a given month). Figure 14 shows histograms of the AOD retrieved by all four sensors from August to October. In all cases a significant fraction of the data retrieve near-zero AOD and have a secondary roughly log-normal distribution of non-zero AOD. All show the decline in AOD (shift to the left) from August to October; the SeaWiFS and VIIRS histograms are shifted to the right and left respectively, compared to the others. Figure 15 shows that these offsets are found across the broader spatial domain, with the four sensors reporting consistent spatial and temporal patterns of AOD. As well as the main smoke plume in the ORACLES domain, a secondary river of smoke outflow into the southern Indian Ocean is seen peaking in September. This feature was also observed by Jethva et al. (2018) and Kar et al. (2018) using OMI and CALIOP data, respectively, and is consistent with known transport patterns (Swap et al., 2003). AOD magnitudes are more different over land, although due to lower cloud cover the data volume is significantly lower and so sampling differences may dominate.

While pixel selection and differences in sensor resolution likely also contribute, the shifts in histogram shape may be plausibly ascribed to uncertainties in the absolute calibration of the sensors. The aerosol signal is small compared to that of the underlying cloud, and a spectral bias in calibration or the retrieval forward model could lead to a systematic bias in the retrieved AOD. This is also an issue with clear-sky AOD retrieval algorithms; e.g. despite identical algorithms there are known systematic offsets in AOD between MODIS Terra and Aqua (Levy et al., 2013; Sayer et al., 2015b). Recent work by Chang et al. (2017) used stable ground sites and identified relative offsets of up to around 2 % in the calibration of these two sensors.

As noted previously, SeaWiFS was calibrated vicariously against ground-based data as described by Franz et al. (2007). In brief, this method assumes that the calibration at the 865 nm band is correct and then adjusts the gain of the other bands such that water-leaving radiance retrievals at this site are unbiased. While an effective method for the ocean colour applications which were the main focus of SeaWiFS, this technique has two main disadvantages for others: first, the untested assumption that the 865 nm band is unbiased; and, second, that the process propagates errors in the ocean colour retrieval atmospheric correction (e.g. aerosol and trace gas scattering/absorption assumptions) into the derived vicarious gain. The latest coefficients used here apply scaling factors of 0.982, 0.9948, and 0.9648 to the SeaWiFS 490, 550, and 670 nm bands respectively, i.e. tilting reflectance downwards at shorter wavelengths compared to 865 nm, which is the direction which would increase the retrieved AAC AOD. If either of the limitations described above are important, this could explain the rightward shift of the SeaWiFS histograms and positive offset seen in the data. Recent work supports these possibilities. Kahn et al. (2016) found spectral biases in the SeaWiFS atmospheric correction, which is used in the vicarious calibration, and Voss and Flora (2017) illustrated that simplifications in current water-leaving radiance processing in the reference data used for the vicarious calibration process lead to small spectral biases.

Figure 15Multiannual monthly mean maps of AAC 550 nm AOD constructed from the data records processed at the present time. Columns show (left–right) data for August, September, and October. Rows show (top–bottom) data for SeaWiFS (1997–2010), MODIS Terra (2000–2017), MODIS Aqua (2002–2017), and VIIRS (2012–2017). Grid cells with fewer than 5 years contributing are shown in grey.

Download

The VIIRS data used here were cross calibrated against MODIS Aqua as described in Sayer et al. (2017), consistent with the main VIIRS Deep Blue data processing. This applied corrections of 0.992, 0.956, 0.941, and an average of 0.963 (with some small temporal dependence) to the 490, 550, 670, and 865 nm bands respectively. This scaling would be expected to decrease the retrieved AAC AOD compared to the uncorrected case. The uncertainties on these corrections were estimated to be ∼0.5 %–1 %, and similar results were found using analyses of cloudy scenes (Kerry Meyer, personal communication, 2018). It is plausible that there is a residual spectral bias in the derived calibration which is leading to biased above-cloud AOD, although Sayer et al. (2017) did find that applying this cross calibration improved clear-sky AOD retrievals, but there was a residual spectral dependence to the AOD bias.

It is difficult to say from the available validation data which data set is closest to the truth. However it seems reasonable to assume that adopting a consistent calibration method for the sensors – whether against a satellite or ground target reference – may improve the consistency of the time series generated. Trace gas absorption corrections can manifest in a similar way to calibration issues, as they are systematic adjustments to bands. Differences between spectroscopic databases or correction parameterisations can also lead to offsets in retrievals (Patadia et al., 2018), so it is also important that these are updated as better spectroscopic measurements or atmospheric reanalyses become available.