Two decades observing smoke above clouds in the south-eastern Atlantic Ocean: Deep Blue algorithm updates and validation with ORACLES field campaign data

Abstract

This study presents and evaluates an updated algorithm for quantification of absorbing aerosols above clouds (AACs) from passive satellite measurements. The focus is biomass burning in the south-eastern Atlantic Ocean during the 2016 and 2017 ObseRvations of Aerosols above CLouds and their intEractionS (ORACLES) field campaign deployments. The algorithm retrieves the above-cloud aerosol optical depth (AOD) and underlying liquid cloud optical depth and is applied to measurements from the Sea-viewing Wide Field-of-view Sensor (SeaWiFS), Moderate Resolution Imaging Spectroradiometer (MODIS), and Visible Infrared Imaging Radiometer Suite (VIIRS) from 1997 to 2017. Airborne NASA Ames Spectrometers for Sky-Scanning, Sun-Tracking Atmospheric Research (4STAR) and NASA Langley High Spectral Resolution Lidar 2 (HSRL2) data collected during ORACLES provide important validation for spectral AOD for MODIS and VIIRS; as the SeaWiFS mission ended in 2010, it cannot be evaluated directly. The 4STAR and HSRL2 comparisons are complementary and reveal performance generally in line with uncertainty estimates provided by the optimal estimation retrieval framework used. At present the two MODIS-based data records seem the most reliable, although there are differences between the deployments, which may indicate that the available data are not yet sufficient to provide a robust regional validation. Spatiotemporal patterns in the data sets are similar, and the time series are very strongly correlated with each other (correlation coefficients from 0.95 to 0.99). Offsets between the satellite data sets are thought to be chiefly due to differences in absolute calibration between the sensors. The available validation data for this type of algorithm are limited to a small number of field campaigns, and it is strongly recommended that such airborne measurements continue to be made, both over the southern Atlantic Ocean and elsewhere.

1 Introduction

Spaceborne monitoring of absorbing aerosols above clouds (AACs), typically smoke or mineral dust aerosols above liquid-phase clouds, has been a topic of increasing research interest in recent years. provide a review of the field, and a more recent list of approaches to their quantification. These AACs are important for multiple reasons. Their direct radiative effects can be very different from those above cloud-free surfaces , and they can have indirect and semi-direct effects on cloud formation, life cycle, and precipitation . Their presence can lead to biases in retrieval of cloud optical depth (COD) and cloud effective radius (CER) if they are not accounted for, as they alter the brightness and spectral shape of the top-of-atmosphere (TOA) signal observed by passive sensors in a systematic way . Additionally, they are largely missing from satellite aerosol optical depth (AOD) data sets derived from passive spaceborne imaging radiometers, which typically process only cloud-free scenes. Global aerosol and cloud fields tend to show similar regional and seasonal variations year after year, and AACs frequently occur downwind of some important aerosol source regions. These include, for example, smoke outflow from south-eastern Asia or southern Africa, as well as dust from the Sahara, Arabian Peninsula, and deserts in north-eastern Asia e.g.. This interannual repeatability means that AOD data sets can have a persistent coverage gap in these regions, which biases estimates of the total atmospheric aerosol burden and hinders aerosol transport analyses.

Semi-quantitative AAC observations from space began with the Total Ozone Monitoring Spectrometer (TOMS) sensor series, which used an ultraviolet aerosol index (UVAI) to take advantage of the spectral darkening of AACs . The large footprint size of TOMS (24–62 $\mathrm{km}$ at nadir, dependent on sensor), however, was a limiting factor to quantitative applications. Similar observations are available from the Global Ozone Monitoring Instrument (GOME) sensor series. While simple to calculate, UVAI is only a semi-quantitative measure of AOD as it depends in a non-linear way on aerosol, cloud, and surface properties as well as solar/view geometry . Quantitative analysis benefited from the 2006 launch of the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP), which is able to provide vertical profiles of aerosol and cloud backscatter and depolarisation , and opened up a new era of quantitative spaceborne AAC research e.g.. More recently, this was supplemented by analyses based on the Cloud Aerosol Transport System (CATS) lidar on the international space station from 2015 to 2017 . While these sensors still have some limitations, the particular features of AACs provide constraints which can obviate some of the assumptions required for these standard backscatter lidar aerosol retrieval algorithms , improving the quantification of AOD and lidar ratio for these cases.

Over the past decade or so, novel algorithmic techniques have been developed for spaceborne quantification of AACs from other sensors. used the improved spatial, spectral, and radiometric capabilities of the Ozone Monitoring Instrument (OMI) over TOMS/GOME to use UVAI to make a more quantitative assessment of the AOD from AACs. This approach was subsequently refined to improve regional assumptions by , enabling global application. also applied a conceptually similar approach to Moderate Resolution Imaging Spectroradiometer (MODIS) measurements. used Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) data to estimate the radiative effect of smoke AACs above the south-eastern Atlantic. Here, AOD and COD were not retrieved, but rather the total shortwave radiative effect was estimated by considering separately those parts of the spectrum measured by SCIAMACHY strongly and weakly influenced by AACs and inferring the aerosol-induced contribution. developed an extension of the MODIS cloud optical properties retrieval algorithm for the south-eastern Atlantic, with a goal to remove the biases in retrieved COD and CER resulting from the neglect of AACs in the standard MODIS cloud data set. developed a similar technique but focused on filling AAC-related gaps in the Deep Blue (DB) aerosol data set. This was demonstrated with MODIS data but was in principle also applicable to the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) and Visible Infrared Imaging Radiometer Suite (VIIRS) sensors to which DB AOD retrieval algorithms have also been applied e.g.. The Polarisation and Directionality of Earth's Reflectance (POLDER) instrument's multidirectional and polarimetric measurement capabilities provide greater information content for aerosols and clouds compared to single-view passive sensors. As a result, several POLDER-based techniques have also been used to quantify AACs .

Figure 1

VIIRS true-colour image from 4 September 2017 showing smoke generated in central/southern Africa transported above marine stratocumulus clouds in the south-eastern Atlantic Ocean. Red dots indicate active fire detections. Region shown corresponds approximately to 36–2${}^{\circ }$ S, 3${}^{\circ }$ W–38${}^{\circ }$ E. Image obtained from NASA Worldview, https://worldview.earthdata.nasa.gov (last access: 27 June 2019).

[Figure omitted. See PDF]

Figure 2

Long-term (2002–2015) mean MODIS Aqua (a) daytime cloud fraction, (b) clear-sky total column AOD at 550 $\mathrm{nm}$, and (c) total fire counts for the month of September for central and southern Africa and surrounding regions. Cloud and aerosol data are at 1${}^{\circ }$ horizontal resolution, while fire counts are at 0.5${}^{\circ }$ horizontal resolution. The green box (25${}^{\circ }$ S–0${}^{\circ }$ N, 15${}^{\circ }$ W–15${}^{\circ }$ E) denotes the approximate region of focus for the ORACLES campaign flights.

[Figure omitted. See PDF]

Much of this research has focussed on African biomass burning. From approximately June to October, agricultural fires move south from central Africa, generating large volumes of smoke which is transported into the south-eastern Atlantic Ocean where it passes over persistent large-scale marine stratocumulus cloud decks . Figure  shows a case from 4 September 2017 where smoke (appearing greyish-brown) from widespread fires is seen blanketing much of Angola and northern Namibia and covering part of a marine stratocumulus cloud deck which has formed along the coastline. Taking a larger perspective, Fig.  shows the long-term (2002–2015) average daytime cloud fraction from the MODIS Collection 6.1 cloud mask;, clear-sky total column AOD at 550 $\mathrm{nm}$ from the MODIS Collection 6.1 Deep Blue/Dark Target merged product;, and cloud-corrected overpass-corrected MODIS Collection 5 fire counts for the month of September. Intense burning across the continent causes large-scale AOD features over land, which are transported both over the stratocumulus deck to the west and in a so-called “river of smoke” to the south-east into the Indian Ocean . Discontinuities in the AOD field in this composite are due in part to sampling (due to the coastal discontinuity in cloud cover), as well as land–ocean algorithm differences. Cloud fraction over portions of the Atlantic approaches 100 %, meaning few clear-sky AOD retrievals are possible; cloud cover over the southern Indian Ocean is lower.

These features make this region a natural laboratory for AAC studies, and several field campaigns have been carried out to better understand aerosol–cloud–precipitation–radiation interactions in this region. Of most interest to the present analysis are the Southern African Regional Science Initiative (SAFARI) year 2000 campaign and the ObseRvations of Aerosols above CLouds and their intEractionS (ORACLES) campaign , which had deployments in the 2016–2018 burning seasons. These campaigns included suites of airborne instrumentation for characterisation of AACs, which have also provided invaluable data for the validation of AAC retrieval algorithms. Indeed, SAFARI-2000 data were used by in the evaluation of the demonstration AAC retrieval algorithm further developed here. Additional field campaigns with different foci related to the southern African aerosol–cloud system have been carried out during the same period as ORACLES ; these include Layered Atlantic Smoke Interactions with Clouds (LASIC), CLoud-Aerosol-Radiation Interactions and Forcing (CLARIFY), and AErosol RadiatiOn and CLOuds in Southern Africa (AeroClo-SA). Deployments and flights generally took place within the area outlined in green in Fig. . The measurements from ORACLES are most directly suited to the evaluation of AAC retrieval algorithms, so they are used here.

The purpose of this study is to describe updates to the initial AAC retrieval algorithm presented by , in preparation for its implementation in the DB aerosol data product suite, and use data collected during the 2016 and 2017 ORACLES deployments to further evaluate the algorithm. The study is organised as follows. Section  describes relevant features of the SeaWiFS, MODIS, and VIIRS satellite sensors; provides a summary of the retrieval algorithm introduced by ; and describes recent updates. Section  details the airborne data obtained during ORACLES and uses these observations to evaluate the updated AAC retrieval algorithm. Finally, the updated algorithm has been applied to process SeaWiFS, MODIS, and VIIRS observations across the large domain (40${}^{\circ }$ S–10${}^{\circ }$ N, 30${}^{\circ }$ W–60${}^{\circ }$ E) shown in Fig.  from the start of the satellite missions to the end of 2017. Section  presents an initial look at this time series and compares results from the different platforms. These AAC retrievals are available upon request to the authors. A separate multi-algorithm comparison exercise is planned for a follow-up study; the purpose here is to introduce, evaluate, and examine the updated algorithm which will eventually be included within the DB data sets. The 2018 ORACLES deployment, and evaluation of the COD retrievals, will likewise be considered in a future study.

2 Satellite AAC retrieval algorithm summary and updates

2.1 Relevant sensor characteristics

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 , MODIS , and VIIRS 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 $\mathrm{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 $\mathrm{nm}$ for MODIS, 490 $\mathrm{nm}$ for SeaWiFS and VIIRS), green (550 $\mathrm{nm}$), red (650 $\mathrm{nm}$ for MODIS, 670 $\mathrm{nm}$ for SeaWiFS and VIIRS), and near-infrared (nIR, 865 $\mathrm{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 $\mathrm{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 $\mathit{\rho }$ for band $i$ is defined as

1$${\mathit{\rho }}_i={\displaystyle \frac{\mathit{\pi }{D}_{\odot }^{\mathrm{2}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}{L}_{\mathit{\lambda }}\left(\mathit{\lambda }\right){\mathrm{\Phi }}_i\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\mathit{\mu }}_{\mathrm{0}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}{E}_{\mathit{\lambda }}\left(\mathit{\lambda }\right){\mathrm{\Phi }}_i\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},$$ where $L_{\mathit{\lambda }}$ is the spectral radiance passing into the satellite field of view at TOA; $E_{\mathit{\lambda }}$ the downwelling solar spectral irradiance at TOA, perpendicular to the Sun and at 1 astronomical unit (AU); and ${\mathrm{\Phi }}_i$ the sensor RSR for band $i$, all functions of wavelength $\mathit{\lambda }$. The factor $D_{\odot }$ is the Earth–Sun distance in AU (variable throughout the year) and ${\mathit{\mu }}_{\mathrm{0}}$ the cosine of the solar zenith angle, which affect the total solar radiation received. Note that $L_{\mathit{\lambda }}$ and so ${\mathit{\rho }}_i$ depend on solar/observation geometry (and of course surface and atmospheric state), omitted here for simplicity of notation.

Figure 3

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

[Figure omitted. See PDF]

For MODIS, nominal horizontal pixel sizes vary from 0.25 to 1 $\mathrm{km}$ (dependent on band); here, the finer-resolution bands are aggregated and coregistered to 1 $\mathrm{km}$. For VIIRS, the nominal pixel size for the relevant bands is 0.74 $\mathrm{km}$. For SeaWiFS, pixel sizes are 1.1 $\mathrm{km}$ but on board resampling performed for Global Area Coverage (GAC) mode subsamples these to provide an effective horizontal resolution of $\sim \mathrm{4.4}$ $\mathrm{km}$ at nadir. As GAC data are a subsampling rather than an average, it is most appropriate to consider these as 1.1 $\mathrm{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 and smallest for VIIRS . Swath widths are 1502 $\mathrm{km}$ for SeaWiFS (GAC mode), 2330 $\mathrm{km}$ for MODIS, and 3040 $\mathrm{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  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. ) 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 $\sim \mathrm{25}$${}^{\circ }$; 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 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 $\mathrm{nm}$; in this study mentions of AOD without a specific wavelength indicated refer to 550 $\mathrm{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 . 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 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 . 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 ; 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 $\mathrm{3}\times \mathrm{3}$, $\mathrm{10}\times \mathrm{10}$, and $\mathrm{6}\times \mathrm{6}$ L1b pixels for SeaWiFS, MODIS, and VIIRS, respectively, giving L2 at-nadir horizontal pixel sizes of 13.5, 10, and 6 $\mathrm{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 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 $\mathrm{nm}$).

retrieved AOD and COD at 550 $\mathrm{nm}$ (the state vector, $\mathit{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 . 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 . OE propagates measurement and forward model uncertainties to provide an estimate ${\mathbf{S}}_x$ of the uncertainty on the retrieved state $\mathit{x}$,

2$${\mathbf{S}}_x={\left({\mathbf{K}}^T{\mathbf{S}}_y\mathbf{K}\right)}^{-\mathrm{1}},$$ where the covariance matrix ${\mathbf{S}}_y$ represents the uncertainty on the measurements (both radiometric and terms arising from forward model assumptions), and $\mathbf{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 . 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: 3$$J\left(\mathit{x}\right)={\left({\mathit{y}}_{\mathrm{m}}\left(\mathit{x}\right)-{\mathit{y}}_{\mathrm{LUT}}\left(\mathit{x}\right)\right)}^T{\mathbf{S}}_y^{-\mathrm{1}}\left({\mathit{y}}_{\mathrm{m}}\left(\mathit{x}\right)-{\mathit{y}}_{\mathrm{LUT}}\left(\mathit{x}\right)\right).$$ This is the sum of square residuals between measured (${\mathit{y}}_{\mathrm{m}}$) and modelled LUT (${\mathit{y}}_{\mathrm{LUT}}$) reflectances relative to their expected magnitudes given in ${\mathbf{S}}_y$, for a given point $\mathit{x}$ in state space (i.e. combination of AOD and COD). Assuming ${\mathbf{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 . 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 $\mathbf{K}$ is calculated numerically. The first guess at $\mathit{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 ofand references therein. The cloud is assumed to consist of a homogeneous and fully overcast layer with a top altitude of 1.5 $\mathrm{km}$ above surface level, geometric thickness of 0.3 $\mathrm{km}$, and be composed of liquid water droplets with size conforming to a gamma distribution with an effective radius of 12 $\mathrm{\mu }$m and effective variance of 0.1. The underlying surface is assumed to be Lambertian (see also Sect.  and ). The aerosol is assumed to lie in a homogeneous layer with a top height 1 $\mathrm{km}$ above the cloud top and be 0.5 $\mathrm{km}$ thick.

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 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 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 . An exception is periodic additional contributions from mineral dust in the northern part from December to February . Over the main ORACLES domain (green box in Fig. ), and during the peak burning season, however, other AAC sources are expected to be negligible.

The Mongu smoke aerosol model was described by . 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 4$${\displaystyle \frac{\mathrm{d}V\left(r\right)}{\mathrm{d}\ln \left(r\right)}}={\displaystyle \frac{C_{\mathrm{v}}}{\sqrt{\mathrm{2}\mathit{\pi }}\mathit{\sigma }}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\ln \left(r\right)-\ln \left(r_{\mathrm{v}}\right)}{\mathit{\sigma }}\right)}^{\mathrm{2}}},$$ for particles of size $r$, given total mode particle volume $C_{\mathrm{v}}$, mode (which is also median and geometric mean) radius $r_{\mathrm{v}}$, and width $\mathit{\sigma }$. found that the modal radius (in $\mathrm{\mu }$m) of the fine mode was dependent on the fine-mode AOD at 550 $\mathrm{nm}$ (${\mathit{\tau }}_{\mathrm{f}}$) as follows: 5$$r_{\mathrm{v},\mathrm{f}}=\mathrm{0.161}+\mathrm{0.013}\ln \left(\mathrm{0.63}{\mathit{\tau }}_{\mathrm{f}}\right).$$

The spread $\mathit{\sigma }$ (dimensionless) of the fine mode was found to have a weak dependence on (${\mathit{\tau }}_{\mathrm{f}}$), 6$${\mathit{\sigma }}_{\mathrm{f}}=\mathrm{0.469}+\mathrm{0.023}\ln \left(\mathrm{0.074}{\mathit{\tau }}_{\mathrm{f}}\right);$$ i.e. for higher smoke loadings the fine-mode particles were larger on average and had a broader distribution. In contrast, $r_{\mathrm{v},\mathrm{c}}$ and ${\mathit{\sigma }}_{\mathrm{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 $\mathrm{\mu }$m respectively. assumed a representative fine-mode fraction (FMF) of AOD at 550 $\mathrm{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 . The assumed aerosol refractive index is shown for AERONET wavelengths in Table . 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. ) is weakly dependent on the AOD and in general varies from $\sim \mathrm{0.86}$ in the blue band to $\sim \mathrm{0.82}$ in the nIR band.

Table 1

Spectral complex refractive index for the smoke aerosol optical model used in this study, following .

Wavelength	Fine mode	Coarse mode
440 $\mathrm{nm}$	1.51–0.024$i$	1.45–0.0035$i$
675 $\mathrm{nm}$	1.52–0.022$i$	1.45–0.0015$i$
870 $\mathrm{nm}$	1.52–0.021$i$	1.45–0.0015$i$
1020 $\mathrm{nm}$	1.52–0.021$i$	1.45–0.0015$i$

Table  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 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  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. ) 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 2

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

Uncertainty in factor	Typical importance to total retrieval uncertainty budget
Sensor absolute calibration	Systematic ($\sim \mathrm{2}$ %) and always a factor
Aerosol SSA	Around 1.5 % at TOA; possible regional and/or seasonal biases from static assumption
Cloud effective radius	Around 1.5 % at TOA; assumption contributes about 20 % of total AOD uncertainty
Cloud structure	Possible 3-D effects and masking errors near edges; hard to quantify
Surface albedo	Negligible ($\ll \mathrm{1}$ % at TOA) provided cloud is opaque (COD $\mathit{⪆}\mathrm{4}$)
Aerosol vertical layering	Potentially significant if multiple layers with distinct SSA/aerosol type present
Aerosol FMF	Generally $<\mathrm{1}$ % at TOA unless aerosol type (smoke/dust) incorrect
Cloud altitude	Generally negligible ($\ll \mathrm{1}$ % at TOA)
Trace gas absorption correction	Generally negligible ($\ll \mathrm{1}$ % at TOA)
Cloud effective variance	Generally negligible ($\ll \mathrm{1}$ % at TOA)
Sensor random error (noise)	Generally negligible ($\ll \mathrm{1}$ % at TOA)
LUT interpolation	Generally negligible ($\ll \mathrm{1}$ % at TOA)
Surface pressure error	Generally negligible ($\ll \mathrm{1}$ % at TOA)

2.3 Algorithm updates since Sayer et al. (2016)

The same overall approach and RT forward model described in Sect.  is used in the present study, with updates described below. These are intended to improve upon approximations made in , 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.1 Pixel selection and aggregation

For the two test case scenes in , the AAC retrieval algorithm was applied to MODIS data at full (nominal 1 $\mathrm{km}$) resolution. Here, to prepare for integration into the main DB/SOAR AOD data sets, the retrieval is instead performed at the equivalent pixel aggregations for the L2 products for each sensor (see Sect. ). This is achieved by taking the median spectral TOA reflectance for suitable L1b pixels within each L2 cell. Use of medians rather than means decreases sensitivity to cloud masking errors and 3-D RT effects which are not accounted for by the forward model . A cell is only processed if the proportion of suitable pixels within the cell is greater than 75 % (i.e. at least 75/100 for MODIS, 48/64 for VIIRS, or 7/9 for SeaWiFS), as the forward model is less appropriate for broken clouds. A suitable pixel is defined as one which is thought to represent a liquid-phase cloud (with or without an overlying absorbing aerosol layer).

For MODIS Terra and Aqua, the standard cloud mask product is used, and cloud phase is taken from the standard MODIS cloud optical properties data sets . Both of these are available within the C6.1 MOD06_L2 (for Terra) and MOD06_L2 (for Aqua) data files. For VIIRS, the equivalent cloud mask product (VNPCLDMK) is used from the current version 1. No VIIRS cloud-phase product is available at the time of writing, so water clouds were identified empirically by assuming that any cloudy pixel with a brightness temperature (BT) in the VIIRS 12 $\mathrm{\mu }$m band below 270 $\mathrm{K}$ corresponded to an ice or mixed-phase cloud and discarding it. While empirical, this threshold seems appropriate in this case based on manual examination of the data, as the vast majority of AAC cases in this region correspond to marine stratocumulus clouds with warmer BTs. For both MODIS and VIIRS, only pixels identified as “probably cloudy” or “confidently cloudy” are considered.

SeaWiFS has no equivalent cloud mask product and lacks the tIR bands useful for determining cloud phase. The historical background for this is that SeaWiFS data were mostly intended and used for monitoring of ocean colour over water, as well as land vegetation indices, for which a clear-sky conservative mask (i.e. few missed clouds) was necessary. Because of this, SeaWiFS cloud masking in those data products is simple and aims to identify and remove not only clouds but aerosol-laden scenes, as well as pixels close to those scenes e.g.. Here, the focus is different, as the desire is to retain optically thick clouds which are likely to be liquid water, and so tests and thresholds are modified, although follow similar principles to the above references. As a result, a separate cloud mask has been developed, drawing from that developed for the DB/SOAR SeaWiFS aerosol products . Specifically, land and water surfaces have different TOA reflectance brightness tests, such that a pixel is marked as cloudy and suitable if

7$${\mathit{\rho }}_{\mathrm{412}}\ge {\displaystyle \frac{\mathrm{0.1}{\mathit{\mu }}_{\mathrm{0}}}{\mathit{\pi }}}$$ over land or 8$${\mathit{\rho }}_{\mathrm{865}}\ge {\displaystyle \frac{\mathrm{0.07}{\mathit{\mu }}_{\mathrm{0}}}{\mathit{\pi }}}$$ over water. The factor of ${\mathit{\mu }}_{\mathrm{0}}$ in the numerator accounts for the fact that reflectance approaches infinity as the Sun approaches the horizon (Eq. ), while with this normalisation the reflectance of an optically thick cloud is less dependent on solar zenith angle. The specific bands chosen for land and water are those at which the surface reflectance tends to be smallest, offering the best discrimination between cloudy and cloud-free scenes, as well as thresholds robust to the presence of AACs.

These tests have been found to be fairly effective at identifying optically thick clouds, over the range of solar zenith angles encountered in the study region (typically from 10 to 60${}^{\circ }$ with an average around 30${}^{\circ }$). Alternate cloud masking strategies may be needed for SeaWiFS for other regions. As there is also no cloud-phase mask for SeaWiFS, additional tests are implemented to identify optically thin clouds (such as cirrus, but also residual optically thin liquid phase). This is based on spatial variability at the 412 $\mathrm{nm}$ band (where clouds tend to show greater spatial variability than cloud-free scenes). This considers $\mathrm{3}\times \mathrm{3}$ groupings of L1b pixels, marking the central pixel cloudy if the test is passed, and again has separate tests for land and water pixels. Over land, the pixel is marked as cloudy but unsuitable if the ratio between the maximum and minimum reflectance at 412 $\mathrm{nm}$ is greater than 1.5 but the absolute brightness test is not passed. Over water, it is marked as cloudy but unsuitable if the standard deviation of reflectance over water pixels within the $\mathrm{3}\times \mathrm{3}$ pixel box is greater than $\mathrm{0.01}{\mathit{\mu }}_{\mathrm{0}}/\mathit{\pi }$ but the absolute brightness test is not passed.

Only detected clouds passing the brightness tests are processed with the AAC retrieval algorithm for SeaWiFS (provided the cell they are in meets the 75 % suitability threshold described above). The TOA reflectance thresholds remove optically thin cirrus clouds from consideration, and output QA filtering (described below) removes others. The spatial variability tests are intended to provide a summary view of the true (i.e. total suitable plus unsuitable) cloud fraction, for comparison with the other sensors. However, this limitation does mean that the total cloud fraction and suitable pixels on which AAC retrieval is attempted may differ between SeaWiFS and the MODIS/VIIRS applications of the algorithm.

2.3.2 Surface elevation

In , the two MODIS test cases examined were predominantly over water, for which the assumption of 1 standard atmosphere pressure is reasonable. This is not necessarily the case over land; Fig. S1 in the Supplement shows that much of the study region is above 1 $\mathrm{km}$ in altitude. Not accounting for this has the potential for regional biases in the algorithm results, as atmospheric pressure determines the total Rayleigh scattering and its interaction with atmospheric multiple scattering and absorption. This could be particularly evident across land–ocean boundaries, e.g. off the coasts of Namibia and Angola where the stratocumulus deck is often encountered. As a result, an additional dimension has been added to the retrieval LUT to account for elevation-dependent changes in surface pressure. Surface elevation ($z$) provided within the L1b files for each pixel is converted to surface pressure ($p$) according to the relationship

9$$p=p{}_{\mathrm{0}}{e}^{\left(\frac{-z}{H}\right)},$$ where the reference pressure $p_{\mathrm{0}}$ is taken as 1013.25 $\mathrm{mbar}$ and the atmospheric scale height $H$ is assumed to be 7.4 $\mathrm{km}$ (sensitivity to this number is small). The LUT contains nodes at 1013.25, 700, and 400 $\mathrm{mbar}$ surface pressure, sufficient to cover the range of elevations encountered here with minimal (generally < 0.5 %) interpolation error in TOA reflectance, and is (as in other dimensions) interpolated linearly.

2.3.3 Ancillary meteorological data

As in the routinely produced DB/SOAR AOD data sets, ancillary meteorological data are needed to correct the TOA reflectance for absorption by trace gases (for the bands considered here, chiefly ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ and ${\mathrm{O}}_{\mathrm{3}}$) and provide a near-surface (10 $\mathrm{m}$) wind speed to calculate Sun glint reflectance over water (see Sect. ). For MODIS and VIIRS, these are obtained from the NASA Goddard Earth Observing System Model, version 5 GEOS-5,, Forward Processing for Instrument Teams (FP-IT) data stream, available from http://gmao.gsfc.nasa.gov/products (last access: 27 June 2019), which is also used in VIIRS DB/SOAR data processing . This begins in 2000, so it is unavailable for the initial years of the SeaWiFS mission; as a result, the Modern-Era Retrospective analysis for Research and Applications, version 2 MERRA2,, available from https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2 (last access: 27 June 2019), is used for the full SeaWiFS record instead. MERRA2 is built on an earlier version of the GEOS-5 model; for the quantities used here (column ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$, column ${\mathrm{O}}_{\mathrm{3}}$, and 10 $\mathrm{m}$ wind speed) the differences between FP-IT and MERRA2 are generally small and the differences introduce negligible additional uncertainty. These fields are at 0.5${}^{\circ }$ latitude by 0.625${}^{\circ }$ longitude resolution, with timesteps of 1 h (MERRA2) and 3 h (FP-IT), and the data are interpolated linearly in space and time to each L1b pixel.

Trace gas absorption correction follows the method and coefficients of , as in VIIRS DB/SOAR data, for MODIS and VIIRS. For SeaWiFS, coefficients from the SeaDAS software (which follows the same basic approach) are used. The purpose of this correction is to simplify retrieval LUT generation by removing the need to include variations in these gas absorbers within the LUTs. The assumption is made that trace gas absorption can be decoupled from other (Rayleigh, aerosol, cloud, surface, and their interaction) scattering and absorption. Then, the TOA reflectances are brightened by dividing by the estimated transmittance as a result of these absorbers, giving the TOA reflectance which would be observed in the absence of these species. For ${\mathrm{O}}_{\mathrm{3}}$ this is reasonable because the bulk of the absorption occurs in the stratosphere and is separated from the bulk of the atmospheric signal; in addition, ozone varies fairly smoothly in space and time. For ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ this is less valid as water vapour varies on finer spatiotemporal scales and is more heterogeneous in its vertical distribution through the atmosphere. Here, as in , the assumption is made that half the water vapour lies below the cloud (and is not seen) and half above. For the bands used in the AAC retrieval, ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ absorption is fairly weak and, except for the nIR band, ${\mathrm{O}}_{\mathrm{3}}$ is the dominant absorber. Total atmospheric gaseous transmittance varies depending on band, solar/view geometry, and atmospheric constituents but generally ranges from $\sim \mathrm{0.99}$ (for the blue bands) to $\sim \mathrm{0.8}$ (for a low Sun and oblique view in the green bands, with a high ozone concentration). Hence, while large errors in the ${\mathrm{O}}_{\mathrm{3}}$ absorption correction are thought to be unlikely, a larger potential error of order 50 % in ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ absorption causes an error of only 1 % or less in TOA reflectance at these bands.

Although ${\mathrm{NO}}_{\mathrm{2}}$ absorbs in the blue part of the spectrum, no absorption correction is applied. This is (as with many other AOD retrieval algorithms) in part due to no availability of this parameter in standard reanalysis data streams ingested for satellite data processing and in part because for the present application it is expected to be a second-order effect. Although potentially significant for fields such as ocean colour analysis near source regions , ${\mathrm{NO}}_{\mathrm{2}}$'s short lifetime means it often has a low abundance away from sources and outside the boundary layer. Since the majority of scenes here are far from potential strong ${\mathrm{NO}}_{\mathrm{2}}$ sources (e.g. industry), and AAC cases are typically around the top of the boundary layer (i.e. the neglected absorption would be below cloud), this is expected to be a second-order contribution to the total uncertainty in TOA reflectance in the blue band.

2.3.4 Land surface reflectance

When a cloud is opaque, the TOA reflectance across the visible part of the spectrum is largely insensitive to the underlying surface albedo. Hence, the demonstration algorithm in made the simplifying assumption of a spectrally neutral surface albedo of 0.05 in all bands. However, when the cloud is optically thin, there is a surface contribution to the TOA signal and assumptions about surface albedo become more important. While the QA tests described below filter out low-COD scenes, if the underlying surface reflectance is brighter than assumed, it is possible that a low-COD cloud could be erroneously retrieved as a combination of higher COD with a higher AOD and pass the QA tests under some circumstances. As a result, in the present study, the surface albedo assumption over land is updated using a climatology derived from MODIS data.

For this purpose the gap-filled snow-free albedo product MCD43GF, is used as a basis. MCD43GF is derived using the MODIS bidirectional reflectance distribution function (BRDF) Terra and Aqua combined product MCD43A1, and applying additional filtering and spatial/temporal constraints to provide BRDF model parameters at 30 arcsec resolution, once every 8 days from 2003 to 2015. Note that the inputs used for the currently available version of MCD43GF derive from the MODIS Collection 5 processing.

While MCD43GF provides full BRDF model parameters, for computational simplicity these are used to calculate white-sky (Lambertian) albedo for use in the retrieval forward model. This approximation is justifiable because under a cloudy sky it is likely that most of the light field below the cloud will be diffuse. For example, for a COD of 1 and vertical incidence only 37 % ($\sim e^{-\mathrm{1}}$) of photons entering the top of the cloud will be directly transmitted without being scattered at least once or absorbed. Above- and below-cloud aerosol and Rayleigh scattering and absorption will further decrease this proportion.

Additionally, to decrease the storage overhead and enable processing outside the 2003–2015 time frame, the source MCD43GF products are spatiotemporally aggregated to provide a database for a representative year (retaining the 8 d time steps) at 0.05${}^{\circ }$ resolution. The spatial aggregation is done first, taking the source MCD43GF products and recording the median albedo within each 0.05${}^{\circ }$ grid cell. After the spatial aggregation, for each grid cell, spectral band, and 8-day period (out of 46 in a year), the median albedo from up to 13 years is taken as representative of that location and time of year. This collapses the interannual variation to provide, for each point, the annual cycle of surface albedo, which is used in the AAC retrieval.

As a measure of the uncertainty introduced by the spatial coarsening, Fig. S2 shows the mean of the spatial standard deviation of albedo within each grid cell. For all bands except 865 $\mathrm{nm}$, this is generally small ($\mathit{<}\mathrm{0.02}$) – even at 865 $\mathrm{nm}$, generally $\mathit{<}\mathrm{0.04}$. This indicates reasonable homogeneity of surface brightness on these scales. The exceptions tend to be salt pans, e.g. the Makgadikgadi Pans in Botswana. Figure S3 shows the mean (across all 8-day periods) temporal standard deviation (across the 13 years) of surface albedo, i.e. a measure of the interannual variability at each location. Spatial patterns are broadly in line with S2, although the magnitudes tend to be slightly higher. This is expected as interannual variability in weather patterns influences vegetation growth and harvest, which influences the surface albedo, especially at 865 $\mathrm{nm}$, which is strongly linked to vegetation cover . assessed the gap-filling procedure in MCD43GF by randomly removing input data and comparing the gap-filled result with that withheld data. For white-sky albedo, this gave root mean square errors (RMSEs) of 0.020 and 0.027 for red and nIR bands respectively. These are similar to or smaller than the quadrature sum of the spatial and temporal aggregation variabilities shown in Figs. S2 and S3. Many of the areas with higher spatiotemporal variability are also associated with lower cloud cover (e.g. Fig. ), meaning they are areas less likely to have an AAC retrieval in the first place, although it is possible that these regions do not represent real cases of smaller variability, but rather more cloudiness means less source data available as input to MCD43GF. did not show results for blue or green bands, but based on the results here it is likely they would be similar or smaller. It is therefore reasonable to assume that the method applied here to generate a climatological database for the AAC retrieval does not significantly degrade the utility of the MODIS albedo product for this application. The resulting annual cycles of surface albedo are shown for four sample locations, representing different surface types, in Fig. S4. In all cases the annual cycle tends to be larger than the interannual variability, which is encouraging as the year-to-year changes are neglected in the present approach.

In the AAC retrieval, surface albedo is assigned at full L1b resolution using the nearest 0.05${}^{\circ }$ grid cell in the climatology from the 8 d period of the year into which the granule falls. Analogously to the aggregation of TOA reflectance, the cell median surface albedo is also used during the retrieval process. The same database is used for all three sensors, as the source BRDF products are at present only available for MODIS, although an equivalent VIIRS data processing suite is in development. Differences in band centres and widths thus have the potential to introduce additional error for SeaWiFS and VIIRS retrievals using this MODIS-derived database, although these are expected to be smaller than 0.02. This is a second-order contribution to the total forward model error in terms of TOA reflectance, especially for optically thick clouds.

2.3.5 Water surface reflectance

Analogously to the over-land surface reflectance treatment, the assumption of a spectrally neutral albedo of 0.05 from is also updated over water surfaces. The reflectance is instead modelled as a combination of a wind-roughened surface using the wind-isotropic model of , with the ancillary data described in Sect.  as input, added to a reflectance of 0.05, 0.04, 0.03, and 0.03 to represent ocean colour and whitecap contributions for the blue, green, red, and nIR bands respectively. Real deviations from this are expected to be of the order $\pm \mathrm{0.01}$, which is again a second-order contribution to the total forward model error in terms of TOA reflectance under optically thin clouds and becomes negligible for opaque clouds.

2.3.6 Retrieval QA

As in , 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 , 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. . Pixels are only retained if the following criteria are all met.

• The retrieval cost (Eq. ) 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 $\ge$ 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 , 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 $\mathrm{3}\times \mathrm{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. ) for these retrievals are often (but not always) large (Fig. i). 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 4

AAC 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 $\mathrm{nm}$. The green box in (a) shows the region used for comparison with airborne data by . Pixels without valid retrievals are shaded in grey.

[Figure omitted. See PDF]

The granule in Fig.  shows one of the test cases compared by against airborne data in SAFARI-2000; in that case, the airborne measurements gave an estimate of the above-cloud AOD at 550 $\mathrm{nm}$ of 0.49, with a spatiotemporal standard deviation of 0.04, within the area outlined by a green box in Fig. a. 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 . 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. noted that the good agreement (AOD within 0.02) for this case may be fortuitous as the estimated uncertainty on the retrieved AOD (Eq. ) is $\pm \mathrm{0.18}$, which is somewhat larger.

3 Validation with airborne data during ORACLES

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

Figure 5

Flight 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.

[Figure omitted. See PDF]

The NASA Ames Spectrometers for Sky-Scanning, Sun-Tracking Atmospheric Research (4STAR) instrument is an aircraft-mountable hyperspectral Sun photometer and sky radiometer . It is a successor to the multichannel Ames Airborne Tracking Sunphotometer (AATS) instruments , which were used by in validation of the initial version of the Deep Blue AAC retrieval algorithm (and cf. Fig. ). 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${}^{\circ }$ , with a radiometric deviation of less than 1 % in this span, compared to 3.7${}^{\circ }$ for AATS . The smaller FOV reduces uncertainties due to scattered light in the direct-beam signal . As operated during ORACLES, 4STAR has 1556 overlapping and continuous bands ranging from 350 to 1700 $\mathrm{nm}$, compared to 6 or 14 distinct non-overlapping spectral bands on the AATS instruments .

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. . More information about the 2016 deployment, as well as 4STAR-derived aerosol data, can be found in . 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 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 $\mathrm{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 and transmission-based cloud property retrievals 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;. HSRL2 provides vertical profiles of atmospheric backscatter, depolarisation, and extinction; an advance of the 2016 deployment was the addition of a 355 $\mathrm{nm}$ channel alongside the 532 and 1064 $\mathrm{nm}$ channels. Note that the 1064 $\mathrm{nm}$ channel lacks HSRL capability and is backscatter only, so above-cloud AOD is only provided at 355 and 532 $\mathrm{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 .

During the 2016 ORACLES deployment , HSRL2 flew on the NASA ER-2 high-altitude aircraft (Fig. ), also based out of Walvis Bay, Namibia. The ER-2 typically flew at an altitude around 20 $\mathrm{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 $\mathrm{km}$ (flight altitudes were generally below 8 $\mathrm{km}$; when not spiraling, a reasonable number of legs were between 5 and 6 $\mathrm{km}$, and few were above 6 $\mathrm{km}$). The current data versions used here are R7 for 2016 and R1 for 2017.

Profiles are measured at 15 $\mathrm{m}$ vertical resolution and 10 $\mathrm{s}$ temporal resolution; the data contain a flag to identify profiles containing AAC cases. The spectral AOD was determined as described in by , 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 and . In brief, there is a random component (which is quantified within the data and typically negligible, $\ll \mathrm{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. evaluated 532 $\mathrm{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. found a similar level of agreement against AERONET from HSRL-2 at both 355 and 532 $\mathrm{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.  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 $\pm \mathrm{3}$ $\mathrm{h}$ is used as a cut-off for a matchup to be valid. This is longer than the $\pm \mathrm{0.5}$–1 $\mathrm{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. ), 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 $\pm \mathrm{1}$ $\mathrm{h}$ (not shown). compared an OMI-based algorithm with HSRL2 measurements from ORACLES 2016; they found (their Fig. 4) that imposing a time difference threshold of $\pm \mathrm{1}$–2 $\mathrm{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 .

Table 3

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

ORACLES	Count
data set	MODIS Terra	MODIS Aqua	VIIRS
4STAR, 2016	532 (20)	432 (15)	835 (17)
4STAR, 2017	190 (12)	285 (15)	561 (18)
HSRL2, 2016	1918 (13)	1896 (14)	4441 (14)
HSRL2, 2017	1066 (16)	156 (10)	484 (12)

The AOD is evaluated at the satellite wavelengths used (i.e. bands centred near 470/490, 550, 650/670, and 865 $\mathrm{nm}$, dependent on sensor), as well as 500 nm, as the latter is (along with 550 $\mathrm{nm}$) a commonly used reference wavelength in aerosol analyses. For the HSRL2 data the available AOD at 355 and 532 $\mathrm{nm}$ are interpolated to 470/490 and 500 $\mathrm{nm}$, and extrapolated to 550 $\mathrm{nm}$, using the Ångström exponent (AE, denoted $\mathit{\alpha }$) where 10$$\mathit{\alpha }=-{\displaystyle \frac{\log \frac{{\mathit{\tau }}_{{\mathit{\lambda }}_{\mathrm{1}}}}{{\mathit{\tau }}_{{\mathit{\lambda }}_{\mathrm{2}}}}}{\log \frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}}}}$$ over the wavelength range ${\mathit{\lambda }}_{\mathrm{1}}$–${\mathit{\lambda }}_{\mathrm{2}}$ (here 355–532 $\mathrm{nm}$). For 4STAR, up to 12 AOD measurements are available across the relevant wavelength range. Therefore, following 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: 11$$\log \left({\mathit{\tau }}_{\mathit{\lambda }}\right)=a_{\mathrm{0}}+a_{\mathrm{1}}\log \left(\mathit{\lambda }\right)+a_{\mathrm{2}}\log (\mathit{\lambda }{)}^{\mathrm{2}}.$$

Coefficients $a_{\mathrm{0}}$, $a_{\mathrm{1}}$, and $a_{\mathrm{2}}$ 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 \left(\mathit{\tau }\right)$ and $\log \left(\mathit{\lambda }\right)$ is not linear but curved, depending on fine-mode particle size . 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. () 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. () (and is typically around $\pm \mathrm{0.01}$). For HSRL2, the uncertainty is taken as $\pm \mathrm{0.02}$ at 470/490 and 500 nm and $\pm \mathrm{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. ) 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., are as follows.

• 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,.

• 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.

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

• 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.

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

• 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 ${\mathit{\sigma }}_{\mathrm{ret}}$ (square root of the relevant element of ${\mathbf{S}}_x$ in Eq. ) and aforementioned airborne uncertainty/variability ${\mathit{\sigma }}_{\mathrm{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, ${\mathit{\tau }}_{\mathrm{ret}}$ and ${\mathit{\tau }}_{\mathrm{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 ($\sim \mathrm{68}$ %) of matchups should be in agreement within this bound and two standard deviations ($\sim \mathrm{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 6

Scatter plots and summary statistics for the comparison between 550 $\mathrm{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 $\mathrm{1}:\mathrm{1}$ line is dotted grey.

[Figure omitted. See PDF]

Figure  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 $\mathrm{1}:\mathrm{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. b) tends to be somewhat smaller than the uncertainty on the individual matchups (vertical bars in Fig. a). 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 7

Summary 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.

[Figure omitted. See PDF]

In Fig. , summary statistics equivalent to those presented in Fig. , 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. ), 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.  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. ) 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. ) 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 ) are reasonable.

Figure 8

Histograms of normalised retrieval error (i.e. actual error divided by expected difference ED) for AOD at 550 $\mathrm{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 $\pm \mathrm{1}$ and $\pm \mathrm{3}$ standard deviations, respectively.

[Figure omitted. See PDF]

Figure 9

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

[Figure omitted. See PDF]

As noted, theoretically the ED should indicate the one standard deviation ($\mathrm{1}\mathit{\sigma }$, $\sim \mathrm{68}$th 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 $+\mathrm{1}$ means that the retrieved AOD was 1 $\times$ ED higher than the airborne AOD for a particular matchup, for example. This distribution is assessed for the 550 $\mathrm{nm}$ data in Fig. . 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  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 $\mathrm{1}\mathit{\sigma }$ errors (i.e. 68th percentile of absolute retrieval error in each bin) and the bottom row the same for $\mathrm{2}\mathit{\sigma }$ (i.e. 95th percentile) errors. For a perfectly characterised retrieval system, these points should lie on the $\mathrm{1}:\mathrm{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 . 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 10

Quantile–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. . The $\mathrm{1}:\mathrm{1}$ line is dotted grey.

[Figure omitted. See PDF]

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. , 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. . 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. ). 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.  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.  shows that the tails of the uncertainty distribution ($\mathrm{2}\mathit{\sigma }$ errors) tend to be quantitatively better estimated than the ($\mathrm{1}\mathit{\sigma }$ 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 11

AE (470–870 $\mathrm{nm}$) assumed in the retrieval as a function of the AOD at 550 $\mathrm{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.

[Figure omitted. See PDF]

The AOD dependence of the size distribution in the aerosol optical model assumed in the retrieval (Sect. ) 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 $\mathrm{nm}$ (Eq. ), is shown in Fig. . The decline from values near 2 in low-AOD conditions to $\sim \mathrm{1.7}$ in high-AOD conditions is a result of the AOD dependence of Eqs. () and ().

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 $\mathrm{nm}$. The second is four Maritime Aerosol Network MAN, 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 $\pm \mathrm{0.02}$ . 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 $\mathrm{nm}$ respectively), which is then subtracted from the total MAN AOD. The remaining AOD is assumed to be smoke, from which the 550 $\mathrm{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 $\pm \mathrm{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 $\mathrm{nm}$ below 0.1, as AE calculation is highly uncertain when the AOD is low e.g., 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 $\mathrm{nm}$) would be expected to result in AE values lower by $\sim \mathrm{0.5}$ compared to the 440–870 $\mathrm{nm}$ spectral range , which would be less directly comparable.

Figure  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 $\sim \mathrm{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 $\mathrm{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, observed a general tendency for increasing AE with altitude for AAC cases (ranging from $\sim \mathrm{1.5}$ at 0.5 $\mathrm{km}$ to $\sim \mathrm{2.0}$ near 4 $\mathrm{km}$). The data collected were most dense from around 1–2 $\mathrm{km}$, over which the AE was fairly flat around 1.7; altitudes below 1 $\mathrm{km}$ or above 3 $\mathrm{km}$ were comparatively poorly sampled and so possibly less representative. This structure is also a secondary contribution to the retrieval uncertainty (Table ). 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 12

Relevant 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 , are in black diamonds. Airborne in situ measurements from SAFARI-2000 reported by are shown for fresh and aged smoke in light and dark green boxes, respectively; airborne remotely sensed measurements from SAFARI-2000 reported by are shown in purple.

[Figure omitted. See PDF]

Similarly to the AE case, the SSA assumed in the AAC retrieval also varies with AOD due to the changing particle size (Eqs.  and ). This AOD dependence is fairly small: a range of $\sim \mathrm{0.02}$ at 470 $\mathrm{nm}$ and $\sim \mathrm{0.04}$ at 870 $\mathrm{nm}$. This is shown in Fig. , 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. 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 .

The first are AERONET inversions from nine sites; four of these sites tend to sample near-source burning (and see Fig. ) 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., and the aerosol loading (particularly in the south-eastern part of the region) can contain additional sulfate contributions from power plants . 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 ; these are pre-filtered to remove poor-quality retrievals, as well as retrievals with an AOD at 440 $\mathrm{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 $\pm \mathrm{0.03}$ . The bulk of this is due to calibration uncertainty and is therefore systematic within a given (roughly year-long) deployment ; 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 . These estimated the SSA at 529 $\mathrm{nm}$ from nephelometer measurements of aerosol scattering and particle soot absorption photometer (PSAP) measurements of absorption. 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 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 $\mathrm{nm}$ for two SAFARI-2000 flights from Table 2 of 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 $\pm \mathrm{0.04}$. The Otavi case (but not the Ascension Island flight) also included nephelometer/PSAP measurements which gave a very similar SSA; noted that refractive index assumptions used in the PCASP calculation were informed by the PSAP data. More recent work by revealed an error in the corrections applied in PSAP data processing which result in a typical underestimate of reported absorption by $\sim \mathrm{13}$ %; for the cases, this translates to an overestimate of SSA by $\sim \mathrm{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 shown here apply corrections accounting for this error.

SAFARI-2000 results for flight a786 near Mongu reported by , 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 $\sim \mathrm{0.02}$ for this case; the best-fit curve is plotted in the figure. also presented results for that flight, and the 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.and references therein there is an increase in SSA from August to October (by $\sim \mathrm{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 $\sim \mathrm{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 $\sim \mathrm{0.06}$ in the blue spectral region and a larger range of up to $\sim \mathrm{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 , 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 $\mathrm{nm}$). As with the results in Fig. , 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 $\mathrm{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. . However, 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 is retained for the present work.

3.5.3 Vertical structure

As noted in Sect. , the algorithm assumes a cloud 0.3 $\mathrm{km}$ thick with a top altitude of 1.5 $\mathrm{km}$ above surface level, with an overlying aerosol layer 0.5 $\mathrm{km}$ thick with a top height 1 $\mathrm{km}$ above the cloud top. Both cloud and aerosol are assumed to be vertically homogeneous. Initial analyses of aerosol/cloud altitudes in this region were generally performed by CALIOP data; provide a summary of some of these results. However, a combination of sensor design and algorithmic limitations means that CALIOP products have been shown to overestimate the bottom of optically thick aerosol layers, leading both to an underestimate of the above-cloud AOD and an overestimate of the gap between cloud top and the overlying aerosol layer bottom . This is due to attenuation of the laser signal through the layer, meaning that the returns from the bottom portion can be too weak for the layer detection algorithm to work correctly. Improvements in the recent CALIOP version 4 data mean that this issue has been slightly ameliorated, although not fully bypassed .

also analysed aerosol/cloud vertical structure in this region using both CALIPSO and CATS data during the 2015–2016 burning seasons (July–October). As the CATS lidar used a wavelength of 1064 $\mathrm{nm}$, the aerosol signal is generally somewhat weaker than at the 532 $\mathrm{nm}$ used by CALIOP, lessening the impact of the attenuation issue on layer detection. Overall, they found CATS reported liquid cloud top heights around 1–1.5 $\mathrm{km}$, typical separations between cloud top and aerosol layer base height of 0.25–0.5 $\mathrm{km}$, and aerosol layer top heights around 3.5–4.5 $\mathrm{km}$ (for a geometric thickness of 1.5–3.25 $\mathrm{km}$). The ranges quoted here arise from longitudinal gradients: as the layer moved west from the coast of Africa into the Atlantic, they found decreases in aerosol top height and increases in cloud top height. Thus, nearer the coast the separation between aerosol and cloud was larger. Meridionally, found higher cloud tops and aerosol layer bases nearer the Equator than toward the southern end of the study region, although the separation between the two layers was relatively constant. One limitation was that to decrease solar noise only nighttime CATS data were used, but as these are large-scale features, and CALIPSO day–night differences were not large (aside from known detection sensitivity issues), it is plausible that these results also hold for daytime measurements. did not examine cloud geometric thickness, although the AAC retrieval algorithm presented here is insensitive to that for opaque clouds.

In light of this, the assumptions made in the AAC retrieval algorithm presented here seem reasonable, although refinements might consider a longitudinal variation of vertical structure and expanding the geometric thickness of the aerosol layer. found that the algorithm was less sensitive to this assumption than other error sources such as SSA assumptions. use a CALIPSO-based climatology of aerosol height data in their OMI data set (but not cloud height or aerosol geometric thickness), which is helpful as retrievals using OMI's UV wavelengths are more sensitive to vertical structure assumptions than the visible/nIR bands used here.

4 A 20-year record from SeaWiFS, MODIS, and VIIRS

4.1 Time series

Figure 13

Monthly time series of various satellite data sets over the green box (25${}^{\circ }$ S–0${}^{\circ }$ N, 15${}^{\circ }$ W–15${}^{\circ }$ E) in Fig. . 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 $\mathrm{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 $\mathrm{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${}^{\circ }$ 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.

[Figure omitted. See PDF]

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. . These are constructed by averaging daily data over the green box (25${}^{\circ }$ S–0${}^{\circ }$ N, 15${}^{\circ }$ W–15${}^{\circ }$ E) in Fig. , 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 . Also shown are the latest version 1.8.9.1 OMI UVAI data, described by . 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. 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 . 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 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${}^{\circ }$ E (as no fires occur over ocean). The fire data are also not presently available for the full MODIS records.

Table 4

Correlation coefficients between monthly mean AAC 550 $\mathrm{nm}$ AOD and other time series shown in Fig. .

Pairing	SeaWiFS	MODIS Terra	MODIS Aqua	VIIRS
MS UVAI	0.84	0.91	0.93	0.91
OMI UVAI	0.92	0.92	0.93	0.89
Fraction MS UVAI $>$ 0.75	0.84	0.86	0.88	0.89
Fraction OMI UVAI $>$ 0.75	0.93	0.92	0.94	0.92
MODIS Terra total column AOD	0.90	0.88	0.91	0.89
MODIS Aqua total column AOD	0.91	0.90	0.91	0.91
MODIS Terra below-cloud AOD	0.13	0.069	0.14	0.23
MODIS Aqua below-cloud AOD	0.062	0.034	0.098	0.25
MODIS Terra fire counts	0.91	0.85	0.86	0.78
MODIS Aqua fire counts	0.90	0.88	0.88	0.84

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 . 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 and ascribed to diurnal variations in fire activity. Table  shows the correlation between each of the time series in Fig.  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. . 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 $\mathrm{nm}$ AOD estimated from MODIS Terra and Aqua are $\mathrm{0.15}\pm \mathrm{0.05}$ and $\mathrm{0.13}\pm \mathrm{0.04}$ respectively, which is only slightly larger than ship-based measurements of AOD in maritime environments without significant continental influence .

4.2 Spatial patterns and offsets

Table 5

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

Pairing	Correlation	Median	rms
		offset	difference
SeaWiFS/MODIS Terra	0.95	0.077	0.097
SeaWiFS/MODIS Aqua	0.96	0.052	0.078
MODIS Terra/Aqua	0.99	$-\mathrm{0.027}$	0.032
MODIS Terra/VIIRS	0.96	0.005	0.029
MODIS Aqua/VIIRS	0.97	0.040	0.052

Figure 14

Histograms of AAC 550 $\mathrm{nm}$ AOD retrieved over the green box (25${}^{\circ }$ S–0${}^{\circ }$ N, 15${}^{\circ }$ W–15${}^{\circ }$ E) in Fig. , 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.

[Figure omitted. See PDF]

Figure  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.  (aside from the SeaWiFS mission which ended in 2010 and so cannot be directly validated with ORACLES data). Table  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  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  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 and using OMI and CALIOP data, respectively, and is consistent with known transport patterns . 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 . Recent work by 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 . In brief, this method assumes that the calibration at the 865 $\mathrm{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 $\mathrm{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 $\mathrm{nm}$ bands respectively, i.e. tilting reflectance downwards at shorter wavelengths compared to 865 $\mathrm{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. found spectral biases in the SeaWiFS atmospheric correction, which is used in the vicarious calibration, and 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 15

Multiannual monthly mean maps of AAC 550 $\mathrm{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.

[Figure omitted. See PDF]

The VIIRS data used here were cross calibrated against MODIS Aqua as described in , 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 $\mathrm{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 $\sim \mathrm{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 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 , so it is also important that these are updated as better spectroscopic measurements or atmospheric reanalyses become available.

5 Conclusions

The ORACLES field campaign and others have provided a wealth of valuable information for the evaluation and refinement of AAC retrieval algorithms for smoke in the south-eastern Atlantic Ocean. This study has detailed updates to an AAC retrieval algorithm and then evaluated it largely using data collected during the 2016 and 2017 ORACLES deployments. This builds on the initial algorithm presented and evaluated with SAFARI-2000 field campaign data by , providing the largest-scale validation possible to date, and can further be supplemented by future analyses of ORACLES and CLARIFY data as these become available. One of the key drives behind the development of this algorithm was to extend coverage of Deep Blue aerosol data products to include AAC cases and thereby fill in some systematic gaps in these global data sets. The algorithm was developed with this in mind, explaining the choice of spatial resolution as well as the spectral range of bands used (470–870 $\mathrm{nm}$). The validation and time series results reveal a reasonable degree of consistency in the resulting data sets, although with some offsets which are likely due to small systematic calibration differences. Calibration assessment and correction (for both absolute calibration and on-orbit degradation) remain a challenge to creating consistent multi-sensor data sets (for AOD and other quantities), and small AOD offsets can persist despite similarities in revealed seasonal and interannual variability. This points to the need for continued traceable calibration against a common reference source, with quantified uncertainties, for satellite measurements in the solar spectrum. Ideally this might be achieved on orbit, as has been done using hyperspectral data for the thermal infrared e.g., in order to enable consistent cross calibration of multiple instruments against a high-quality reference. Consistency of calibration is expected to be the primary factor driving offsets between the retrieved time series.

Overall, the validation and comparison exercise has revealed that the AAC algorithm presented here performs roughly within expectations, based on sensitivity analyses and the uncertainty estimates provided via the optimal estimation retrieval technique. Specific areas for potential refinement have been identified, chiefly sensor calibration and potential adjustments to assumed aerosol and cloud optical properties and structure. As one example, the optical model used in this region could account for intra-annual variations in SSA following , or ancillary data sources could be used as a switch between possible optical models as in . SeaWiFS would also benefit from a cloud mask of the same calibre as that available to MODIS and VIIRS. These refinements would be expected to improve the consistency between the different sensors to which the algorithm has been applied and reduce some sources of systematic uncertainty at certain times and locations. Moving forward to a global application would require the development of equivalent appropriate assumptions globally, which can be done by leveraging climatologies of vertical structure from lidar (spaceborne or ground-based), as well as representative aerosol optical properties from AERONET and potentially global model-based climatologies e.g.. Similar spectral bands are also present on the new generation of geostationary sensors launched in recent years, providing a rapid revisit of some important AAC systems . It would be advantageous to apply Deep Blue and this AAC algorithm to those sensors, improving knowledge of the diurnal cycle of AACs, which is not fully sampled by the instruments on Sun-synchronous platforms used here. In the meantime, the AAC data set generated in this work is available for interested researchers.

Unfortunately, the available validation data for these algorithms remain highly sparse. The results here suggest that the available ORACLES flights, while a significant important milestone and far ahead of the characterisation of other AAC systems, may not yet represent sufficient sampling to provide a robust regional validation. Validation and assessment of derived product uncertainties is inherently a statistical exercise, both due to the nature of an uncertainty estimate and the difficulty in simultaneously characterising each factor affecting the TOA satellite signal to a degree sufficient to confidently ascribe the reason(s) for errors on individual retrievals. The available validation from field campaigns in other regions is much more limited. Aircraft observations are a powerful tool to provide data-rich, thorough characterisation of sampled air masses. Such campaigns could be supplemented by instrumentation carried upon unmanned aerial vehicles (UAVs), a technology which has advanced greatly in recent years for a recent review see. Frequent launches of instrumented UAVs, combined with geostationary satellites, could provide an important temporal sampling component to further refine understanding of the processes influencing the evolution of these systems. It is critical, both for answering science questions about the role of AACs in the Earth system and for a robust quantitative validation of spaceborne AAC data sets, that such observations continue to be made.

Data availability

The Deep Blue aerosol-above-cloud retrieval data shown in this work are available upon request to the authors, in advance of their eventual incorporation to the main Deep Blue data products. ORACLES airborne data are registered as . More information about Deep Blue is available from https://deepblue.gsfc.nasa.gov (last access: 27 June 2019). ORACLES field campaign data and further information are available from https://espo.nasa.gov/oracles (last access: 27 June 2019). ORACLES investigators and their team members are thanked for their considerable efforts in planning, executing, and archiving the results of the ORACLES field campaign deployments. SeaWiFS data are available from https://oceancolor.gsfc.nasa.gov (last access: 27 June 2019). MODIS and VIIRS atmospheres data are available from https://ladsweb.modaps.eosdis.nasa.gov/ (last access: 27 June 2019). NASA GMAO are thanked for GEOS-5 FP-IT and MERRA2 reanalysis data, available from https://gmao.gsfc.nasa.gov (last access: 27 June 2019). Crystal Schaaf (UMB) is thanked for useful discussions about the MODIS albedo products; the MCD43GF data used in this study are available from https://www.umb.edu/spectralmass (last access: 27 June 2019). AERONET and MAN data are available from https://aeronet.gsfc.nasa.gov (last access: 27 June 2019); Alexander Smirnov (SSAI) and the AERONET/MAN investigators are thanked for the creation and stewardship of these data records, as well as useful discussions concerning data from the region of interest. MODIS fire and OMI data are available from https://earthdata.nasa.gov/about/daacs/daac-ges-disc (last access: 27 June 2019) and MS UVAI data from http://www.icare.univ-lille1.fr (last access: 27 June 2019).

The supplement related to this article is available online at: https://doi.org/10.5194/amt-12-3595-2019-supplement.

Author contributions

AMS developed the algorithm, performed the analysis, and wrote the manuscript. NCH is the PI of the Deep Blue aerosol project; NCH, JL, and WVK are the other key developers for the DB and SOAR clear-sky aerosol retrieval algorithms and provided guidance. SCT provided scientific suggestions and general guidance through the project and algorithm development. SB, MAF, and RAF were involved in HSRL2 data collection, processing, and analysis. MK, SL, KP, JR, MSR, and YS were involved in 4STAR data collection, processing, and analysis. All authors contributed to reviewing and editing the manuscript.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “New observations and related modelling studies of the aerosol–cloud–climate system in the Southeast Atlantic and southern Africa regions (ACP/AMT inter-journal SI)”. It is not associated with a conference.

Acknowledgements

This research was funded by NASA's Radiation Science Program, managed by Hal Maring. The VIIRS Atmospheres SIPS at the University of Wisconsin are thanked for assistance in obtaining and processing VIIRS data. Richard Frey (University of Wisconsin), Steve Platnick (NASA GSFC), and Kerry Meyer (NASA GSFC) are thanked for access to and advice about the VIIRS cloud mask data product. Tom Eck (USRA), Kate Szpek (UK Met Office), and Connor Flynn (PNNL) are thanked for useful discussions about remote and in situ measurements of aerosol SSA in the study region. We are grateful to three anonymous reviewers for comments that improved the clarity and readability of the paper.

Financial support

This research has been supported by NASA (Radiation Science Program).

Review statement

This paper was edited by J. M. Haywood and reviewed by three anonymous referees.