Retrievals of aerosol optical depth over the western North Atlantic Ocean during ACTIVATE

Aerosol optical depth was retrieved from two airborne remote sensing instruments, the Research Scanning Polarimeter (RSP) and Second Generation High Spectral Resolution Lidar (HSRL-2), during the National Aeronautics and Space Administration (NASA) Aerosol Cloud meTeorology Interactions oVer the western ATlantic Experiment (ACTIVATE). The field campaign offers a unique opportunity to evaluate an extensive 3-year dataset under a wide range of meteorological conditions from two instruments on the same platform. However, a long-standing issue in atmospheric field studies is that there is a lack of reference datasets for properly validating field measurements and estimating their uncertainties. Here we address this issue by using the triple collocation method, in which a third collocated satellite dataset from the Moderate Resolution Imaging Spectroradiometer (MODIS) is introduced for comparison. HSRL-2 is found to provide a more accurate retrieval than RSP over the study region. The error standard deviation of HSRL-2 with respect to the ground truth is 0.027. Moreover, this approach enables us to develop a simple, yet efficient, quality control criterion for RSP data. The physical reasons for the differences in two retrievals are determined to be cloud contamination, aerosols near the surface, multiple aerosol layers, absorbing aerosols, non-spherical aerosols, and simplified retrieval assumptions. These results demonstrate the pathway for optimal aerosol retrievals by combining information from both lidars and polarimeters for future airborne and satellite missions.

1 Introduction

Aerosol particles constitute an important component of the Earth's atmosphere by altering its radiative energy balance. Aerosols impact our climate system directly by scattering and absorbing solar and terrestrial radiation and indirectly by serving as cloud condensation nuclei (CCN) and ice nuclei (IN) (Lohmann and Feichter, 2005; Pöschl, 2005). Changing CCN and IN concentrations can influence the distribution of liquid water, ice, and mixed-phase clouds and the frequency of precipitation. The complex interaction of aerosols with clouds, radiation, and meteorology makes it difficult to probe the feedback and response of our climate system under a perturbation of anthropogenic or natural aerosols (Bellouin et al., 2020; Chen et al., 2014). In the most recent Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC), aerosols still represent the largest uncertainty in the global radiative forcing (IPCC, 2023).

To disentangle aerosol–cloud interactions, it is imperative to collect adequate observations for robust statistics across a wide range of cloud regimes. The recent Aerosol Cloud meTeorology Interactions oVer the western ATlantic Experiment (ACTIVATE) epitomizes a coordinated effort to respond to this need (Sorooshian et al., 2023). Being one of the National Aeronautics and Space Administration (NASA) Earth Venture Suborbital-3 (EVS-3) missions, ACTIVATE has three objectives:

• first, to quantify the underlying relationship between the aerosol number concentration, CCN concentration, and cloud drop number concentration;

• second, to develop a better process-level representation of cloud properties and aerosol–cloud interactions in a hierarchy of numerical models;

• finally, to evaluate state-of-the-art remote sensing instruments that are built for retrieving the aerosol and cloud properties (Sorooshian et al., 2019).

Aerosol optical depth (AOD) is a vertically integrated quantity of aerosol extinction coefficient, representing the column aerosol loading. It is one of the most widely used parameters to monitor the long-term evolution of aerosols at both global and regional scales (Andreae, 2009; Seinfeld et al., 2016). It is important to evaluate AOD retrievals and their associated uncertainties. An uncertainty of AOD retrievals of 0.01 could lead to an uncertainty of aerosol radiative forcing of ∼ 0.5 W m⁻² (Chylek et al., 2003; Hansen et al., 1995; Mishchenko et al., 2004). Two approaches are regularly used to quantify uncertainties of AOD retrievals: propagated (prognostic) uncertainty and truth-based (diagnostic) uncertainty (Gao et al., 2022; Sayer et al., 2020). The former is a theoretical estimate which is obtained by considering measurement uncertainties, forward model uncertainties, and a priori assumptions via minimizing a suitable cost function. The latter is to compare the retrieval with a reference dataset. The two approaches are complementary to each other as the reference dataset can assess the theoretical estimate using appropriate statistical methods (Gao et al., 2022).

This work is focused on the second approach as two remote sensing instruments, the Research Scanning Polarimeter (RSP) and Second Generation High Spectral Resolution Lidar (HSRL-2), retrieved AOD over the western North Atlantic Ocean during ACTIVATE. Before proceeding, however, there is a catch, that is, the following question: how do we objectively assess these two datasets? Ideally, the error of a dataset should be quantified by comparing it with the “ground truth” (Caires and Sterl, 2003). Practically, the ground truth for AOD measurements is replaced by measurements from some ground-based instruments such as sun photometers in the AErosol RObotic NETwork (AERONET; Holben et al., 1998). But ground-based observations are unavailable over the study region, primarily an open ocean. Furthermore, all instruments and retrievals are subject to errors consisting of systematic bias and random noises (Stoffelen, 1998).

We propose to use a statistical method called triple collocation (TC) which partially circumvents the aforementioned issues. A third dataset is introduced, and the error characteristics of each dataset are determined with respect to the unknown ground truth. The most accurate dataset is then used to help design a simple, yet efficient, data flag to improve the quality of the other datasets. The physical processes that lead to the difference between RSP and HSRL-2 retrievals are discussed along with three case studies.

2 Data

2.1 Overview of the ACTIVATE field campaign

The ACTIVATE field campaign started in February 2020 and concluded in June 2022 with six deployments. There were winter and summer deployments each year, roughly covering November to March and May to September, respectively (Sorooshian et al., 2023). While it is recommended to explicitly label the year and season of each deployment, for simplicity we refer to the three winter and three summer deployments collectively as the “winter” and “summer” deployments.

Two aircraft, the King Air (high-flying) and Falcon (low-flying) from the NASA Langley Research Center (LaRC), were used to carry out spatially coordinated flights over the western North Atlantic Ocean. The total number of flights for the King Air and Falcon are 174 and 168, respectively, with 162 being joint flights. The remote sensing instruments in this study were flown on the King Air which logged a total of 592 flight hours (Sorooshian et al., 2023).

Most of the 3–4 h flights were based at NASA LaRC, with a notable exception being 17 flights based in Bermuda during June 2022 (Fig. 1). Due to air traffic considerations associated with military restricted areas, the aircraft usually had to transit through one of three waypoints, OXANA, ZIBUT, and ZIZZI.

Figure 1 The ACTIVATE study region and NASA Langley King Air flight tracks between 2020 and 2022. Summer and winter deployments are indicated in red and blue, respectively. The two campaign bases, NASA Langley Research Center (LaRC) and Bermuda, and three commonly used waypoints, OXANA, ZIBUT, and ZIZZI, are annotated with arrows.

2.2 Research Scanning Polarimeter

The Research Scanning Polarimeter (RSP) is a remote sensing instrument that uses three pairs of optical assemblies to simultaneously measure linear polarization of the intensity at four polarization azimuths and nine spectral channels centered at 410, 470, 555, 670, 865, 960, 1590, 1880, and 2264 nm for each scene (Cairns et al., 1999, 2003). The first three components of the Stokes vector (I, Q, and U) are then derived from the measurement. RSP operates by continuously scanning the field of view (14 mrad) over ± 55° from viewing zenith (∼ 140 views) along the flight track, which samples at 0.8° interval (Alexandrov et al., 2012; Cairns et al., 2003). Each RSP scan consists of measurements from all different viewing angles. That is, each target at the cloud top or ground surface is swept from different scans at different viewing angles. These actual scans are assembled into virtual scans for each target at nadir (Alexandrov et al., 2012).

The aerosol properties are retrieved using the RSP Microphysical Aerosol Properties from Polarimetry (MAPP) algorithm version 1.48, consisting of a vector radiative transfer model and a Mie code (Stamnes et al., 2018; Schlosser et al., 2022). The retrieval algorithm is geared for field campaigns which desire both speed and flexibility. Speed is important for field campaigns because data are needed as soon as possible for preliminary data analysis and flight planning. To speed up computations, the Mie code assumes particles are spherical and homogeneous. Flexibility is also important because the traditional pre-computed Mie lookup table approach requires constant aircraft altitudes and aerosol locations, which limits the choice of flight patterns. To remain flexible, both the vector radiative transfer and Mie computations are done on the fly.

The retrieval algorithm uses an optimal-estimation method (Rodgers, 2000) to minimize a cost function χ² that describes the difference between the observations and a forward model for aerosol parameter retrievals and includes a term that takes into account of a priori knowledge of aerosol properties. The cost function χ² is formulated and described in detail by Stamnes et al. (2018) and includes a metric describing the goodness of fit between the forward model and RSP measurements, called the normalized cost function χ’. Here we only consider retrievals successful when χ’ is less than 0.15. As regards the a priori best guess, we adopt an elementary and conservative approach by taking the mean of the permissible range of each aerosol parameter, which is also utilized as the first guess.

The current version of algorithm assumes that the aerosol size distribution is bimodal: one fine mode and one coarse mode (Stamnes et al., 2018). Aerosol properties such as AOD, refractive index, and effective radius are retrieved separately for the two modes. The coarse mode is assumed to be composed of non-absorbing sea salt particles with a refractive index equivalent to that of water. The coarse-mode aerosol top height is assumed to be 1 km, while the fine-mode aerosol top height is retrieved (Schlosser et al., 2022).

The theoretical ±1σ accuracy of the MAPP algorithm for AOD is ± 0.02 (Stamnes et al., 2018). The practical accuracy of RSP has been estimated in several field campaigns through comparison with sun photometers, but most of them took place over land. Knobelspiesse et al. (2011) found that the mean error of RSP relative to the Ames Airborne Tracking Sunphotometer (AATS-14) is 0.057 at 532 nm for a forest fire event during the Arctic Research of the Composition of the Troposphere from Aircraft and Satellites (ARCTAS) field campaign. Di Noia et al. (2017) compared RSP aerosol retrievals against AERONET during the POlarimeter DEfinition EXperiment (PODEX) and Studies of Emissions and Atmospheric Composition, Clouds and Climate Coupling by Regional Surveys (SEAC⁴RS) field campaigns in 2012 and 2013, respectively. The bias and root mean squared error of these retrievals at 500 nm are 0.02 and 0.04, respectively. Fu et al. (2020) compared measurements from four multi-angle polarimeters, including RSP, against AERONET during the Aerosol Characterization from Polarimeter and Lidar (ACEPOL) field campaign over the western part of the United States in 2017. The standard deviations between RSP and AERONET are 0.035, 0.027, and 0.014 at 380, 440, and 675 nm, respectively.

In this study, we use the RSP level 2 aerosol product which includes column-averaged aerosol optical and microphysical parameters and ocean color parameters. The RSP data are averaged horizontally over 10 scans from the level 1C product, with a temporal resolution of ∼ 9.3 s and an along-track spatial resolution of ∼ 1 km.

2.3 Second Generation High Spectral Resolution Lidar

The Second Generation High Spectral Resolution Lidar (HSRL-2) is an airborne multi-wavelength instrument which uses the HSRL technique to separate the aerosol and molecular backscatter signals at 355 and 532 nm and the standard backscatter lidar technique at 1064 nm. HSRL-2 can vertically resolve the backscatter coefficient β and depolarization δ at 355, 532, and 1064 nm and extinction coefficient α at 355 and 532 nm (Burton et al., 2018; Hair et al., 2008; Ferrare et al., 2023). Because of the incorporation of three backscatter and two extinction measurements, the instrument is also referred to as a “3β+2α” lidar (Burton et al., 2016). The field of view of HSRL-2 is 1 mrad (Schlosser et al., 2022). AOD is derived from the aerosol backscattering coefficient using the difference in molecular return signals (Hair et al., 2008).

The performance of HSRL-2 measurements has been compared with various instruments. For example, Shinozuka et al. (2013) showed that the root mean squared difference of 532 nm AOD between the Spectrometer for Sky-Scanning, Sun-Tracking Atmospheric Research (4STAR) and HSRL-2 is 0.01 during the 2012 Two-Column Aerosol Project (TCAP) campaign. Sawamura et al. (2017) showed outstanding agreement between the AERONET and HSRL-2 from over 300 profiles collected during the 2013 Deriving Information on Surface Conditions from COlumn and VERtically Resolved Observations Relevant to Air Quality (DISCOVER-AQ) campaign. Stamnes et al. (2018) found that the root mean squared deviations between RSP and HSRL-2 AOD at 532 nm are ∼ 0.07 during TCAP and ∼ 0.04 during the 2014 Ship-Aircraft Bio-Optical Research (SABOR) campaign.

In this study, we use the 10 s aerosol profile product of which spatial resolution is ∼ 1 km. We also use the following derived products: mixing layer height (MLH) and aerosol identifier. MLH is estimated from the sharp gradients of HSRL-2 aerosol backscatter profiles at 532 nm (Scarino et al., 2014). A Haar wavelet transform is first applied to the vertical profiles to identify the transition zone above the top of the aerosol mixed layer, and then MLH is defined as the altitude of the maximum value after the wavelet transform. The information of aerosol identifiers (or types) is based on the HSRL aerosol classification scheme (Burton et al., 2012) which identifies aerosols using the lidar ratio, depolarization, backscatter color ratio, and spectral depolarization ratio from HSRL-2 measurements. Eight types of aerosols are classified, including dust, polluted marine, fresh smoke, smoke, urban, marine, a dusty mix, and ice. The algorithm has been shown to improve a similar scheme for the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) instrument (Burton et al., 2013; Ferrare et al., 2023).

To expedite comparison with gridded global reanalysis data, selected instantaneous meteorological fields from the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2; Gelaro et al., 2017), are interpolated to the HSRL-2 curtains along the flight tracks. These 3-hourly MERRA-2 variables have a horizontal resolution of 0.5° × 0.625° (latitude × longitude) on 42 pressure levels from 1000 to 0.1 hPa.

Figure 2 Distribution of RSP (blue) and HSRL-2 (red) AOD for the whole field campaign period (2020–2022). Each bin is 0.025. Mean values (dashed colored lines) and selected statistics (mean, median, SD, and IQR) are also shown.

2.4 Other ACTIVATE datasets

The King Air navigational data for RSP and HSRL-2 are provided by an Applanix 610 system at 2 Hz, including flight date, time, longitude, latitude, altitude, and other aircraft motion variables (Sorooshian et al., 2023). Two cameras were mounted on the King Air: one facing nadir and the other facing forward (Sorooshian et al., 2023). The images collected from the airborne cameras are used to identify the presence of clouds above and below the aircraft. Recently, a cloud detection neural network (CDNN) algorithm was developed to speed up the cloud identification process by training and testing the camera imageries via convolutional neural network, which resulted in above-aircraft and below-cloud cloud mask products at 1 Hz (Nied et al., 2023). The accuracy of the forward-facing CDNN is 96 %, calculated through validation against human-labeled testing data. The Falcon relative humidity data are provided by a diode laser hygrometer (DLH) with an uncertainty of 0.1 ppmv (Diskin et al., 2002).

Table 1 Monthly RSP and HSRL-2 AOD statistics. SD and IQR denote standard deviation and interquartile range, respectively.

2.5 Moderate Resolution Imaging Spectroradiometer

The Moderate Resolution Imaging Spectroradiometer (MODIS) is an imaging instrument with 36 spectral channels ranging from 0.41 to 14.5 µm. The instrument collects data at three nadir spatial resolutions (250 m, 500 m, and 1 km) (Remer et al., 2005). As part of the NASA Earth Observing System, two almost identical MODIS sensors were launched on the Terra and Aqua satellites in December 1999 and May 2002, respectively. Both satellites observe the Earth from a polar-orbiting, sun-synchronous orbit at an altitude of ∼ 700 km. The equatorial overpass times for the descending mode of Terra and the ascending mode of Aqua are 10:30 and 13:30 local solar time, respectively. The swath of each scan is 2330 km × 10 km (across track × along track at nadir) (Levy et al., 2013). Seven bands are used for aerosol retrievals.

In this study, we use the following level 2 aerosol products from the MODIS Collection 6.1: MOD04_3K from Terra and MYD04_3K from Aqua. The horizontal resolution is 3 km × 3 km at nadir. The Dark Target (DT) algorithm is used to retrieve the 3 km aerosol product (Levy et al., 2013; Remer et al., 2013). The algorithm is complicated (see, for example, Remer et al., 2005). Generally, each grid box consists of $\mathrm{6}\times \mathrm{6}=\mathrm{36}$ pixels (i.e., each pixel is at 500 m resolution at nadir). The algorithm checks the number of good pixels in each grid box using multiple criteria before doing the actual inversion. Therefore, the existence of cloudy pixels does not stop the retrieval as long as there are enough good pixels within the same grid box (Remer et al., 2013). Note also that the cloud pixels would not be included during the inversion.

To compare with the RSP and HSRL-2 AOD at 532 nm, the MODIS AOD at 550 nm is converted to 532 nm using the Ångström exponent from the MODIS AOD at 470 and 550 nm. Compared to AERONET, the expected error for the DT algorithm is $\pm (\mathrm{0.03}\pm \mathrm{0.05}\times \text{AOD}$) (Remer et al., 2013).

2.6 Cloud-Aerosol Lidar with Orthogonal Polarization

The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) is the primary instrument that was launched on the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) satellite in April 2006 (Winker et al., 2009). We use the CALIOP level 3 stratospheric aerosol profile products which are based on version 4 of CALIOP level 1B and level 2 datasets (Kar et al., 2019). The horizontal resolution of the monthly stratospheric AOD data is 5° × 20° (latitude × longitude). To maximize the length of available data, we stitch the datasets from versions 1.00 and 1.01. The version 1.01 begins in July 2020.

3 Evaluation methods

3.1 Traditional metrics

To summarize the distribution of an AOD dataset, it is common to use mean and standard deviation (SD) to measure its central tendency and dispersion, which are known to be affected by outliers. The median and interquartile range (IQR= the third quartile minus the first quartile) are also used as they are more robust against outliers.

To evaluate AOD between two collocated datasets X and Y, each having n elements, we use the following metrics:

$$\begin{array}{}\text{(1)}& {\displaystyle }r& {\displaystyle }={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(X_i-\stackrel{\mathrm{‾}}{X})(Y_i-\stackrel{\mathrm{‾}}{Y})}{\sqrt{{\sum }_{i=\mathrm{1}}^n(X_i-\stackrel{\mathrm{‾}}{X}{)}^{\mathrm{2}}}\sqrt{{\sum }_{i=\mathrm{1}}^n(Y_i-\stackrel{\mathrm{‾}}{Y}{)}^{\mathrm{2}}}}},\text{(2)}& {\displaystyle }\text{MB}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(X_i-Y_i)=\stackrel{\mathrm{‾}}{X}-\stackrel{\mathrm{‾}}{Y},\text{(3)}& {\displaystyle }\text{RMSD}& {\displaystyle }=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(X_i-Y_i{)}^{\mathrm{2}}},\end{array}$$

where r, MB, and RMSD are Pearson's correlation coefficient, the mean bias, and the root mean squared deviation, respectively, and the overbar denotes the mean value of a dataset.

Figure 3 Two-dimensional histograms of RSP and HSRL-2 AOD. (a) All data. (b) Winter deployments only. (c) Summer deployments only. Bin size is 0.025 for both instruments. In each plot, the linear fit (dashed black), 1 : 1 line (thick dashed gray), ±0.1 intervals (thin dashed gray), and selected error metrics are shown.

3.2 Triple collocation

Triple collocation (TC) is a statistical method that estimates error statistics of a target variable from three independent and collocated datasets (Stoffelen, 1998). The method does not require any one of the triplet items to be the ground truth. It has been used for characterizing errors in a wide range of geophysical variables such as sea surface temperature (O'Carroll et al., 2008), soil moisture (Draper et al., 2013), and ocean wave height (Caires and Sterl, 2003; Janssen et al., 2007).

In this study, we take advantage of three available AOD retrievals from RSP, HSRL-2, and MODIS. To start with, each dataset is considered as an estimate of the (unknown) ground truth. Therefore, a functional relationship (or an error model) is assumed to exist between each dataset and the ground truth to account for the potential deviation from the truth. The linear relationship is one of the simplest error models. Some linear error models are discussed in detail in Zwieback et al. (2012). We adopt an affine (also linear) model for the AOD datasets:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\tau }}_i=a_i+b_i\mathrm{\Theta }+{\mathit{ϵ}}_i,\end{array}$$

where τ_i ($i\in \mathit{\{}\text{RSP},\text{HSRL-2},\text{MODIS}\mathit{\}}$) represents the AOD from dataset i; Θ is the true AOD; and a_i, b_i, and ϵ_i are the additive bias, multiplicative bias, and random error from dataset i, respectively. To simplify the notation, we assign RSP=1, HSRL-2=2, and MODIS=3.

Together with the error model being linear, the following assumptions are made (Gruber et al., 2016). First, the random error has zero mean (E[ϵ_i]=0). Second, the random error is uncorrelated with the true AOD $\left(\text{Cov}\right[{\mathit{ϵ}}_i,\mathrm{\Theta }]=\mathrm{0})$. Third, the random error of different products is uncorrelated with each other $\left(\text{Cov}\right[{\mathit{ϵ}}_i,{\mathit{ϵ}}_j]=\mathrm{0},i\ne j)$. Lastly, the additive and multiplicative biases are time invariant.

The equations of variance and covariance between the collocated datasets can then be obtained (McColl et al., 2014; Stoffelen, 1998) by

$$\begin{array}{ll}{\displaystyle }\text{Cov}[{\mathit{\tau }}_i,{\mathit{\tau }}_j]& {\displaystyle }=b_i{b}_j\text{Var}\left[\mathrm{\Theta }\right]+b_i\text{Cov}[\mathrm{\Theta },{\mathit{ϵ}}_j]+b_j\text{Cov}[\mathrm{\Theta },{\mathit{ϵ}}_i]\\ \text{(5)}& {\displaystyle }& {\displaystyle }+\text{Cov}[{\mathit{ϵ}}_i,{\mathit{ϵ}}_j],\text{(6)}& {\displaystyle }\text{Var}\left[{\mathit{\tau }}_i\right]& {\displaystyle }=b_i^{\mathrm{2}}\text{Var}\left[\mathrm{\Theta }\right]+\mathrm{2}{b}_i\text{Cov}[\mathrm{\Theta },{\mathit{ϵ}}_i]+\text{Var}\left[{\mathit{ϵ}}_i\right],\end{array}$$

where $\text{Var}\left[\mathrm{\Theta }\right]={\mathit{\sigma }}_{\mathrm{\Theta }}^{\mathrm{2}}$ is the truth variance and $\text{Var}\left[{\mathit{ϵ}}_i\right]={\mathit{\sigma }}_{{\mathit{ϵ}}_i}^{\mathrm{2}}$ is the error variance. Note that $\text{Cov}[{\mathit{\tau }}_i,{\mathit{\tau }}_j]=\text{Cov}[{\mathit{\tau }}_j,{\mathit{\tau }}_i]$.

The error variances are then solved with the help of the above assumptions (McColl et al., 2014):

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{\sigma }}_{{\mathit{ϵ}}_{\text{RSP}}}^{\mathrm{2}}& {\displaystyle }=\text{Var}\left[{\mathit{\tau }}_{\mathrm{1}}\right]-{\displaystyle \frac{\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{2}}]\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{3}}]}{\text{Cov}[{\mathit{\tau }}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{3}}]}},\text{(8)}& {\displaystyle }{\mathit{\sigma }}_{{\mathit{ϵ}}_{\text{HSRL-2}}}^{\mathrm{2}}& {\displaystyle }=\text{Var}\left[{\mathit{\tau }}_{\mathrm{2}}\right]-{\displaystyle \frac{\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{2}}]\text{Cov}[{\mathit{\tau }}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{3}}]}{\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{3}}]}},\text{(9)}& {\displaystyle }{\mathit{\sigma }}_{{\mathit{ϵ}}_{\text{MODIS}}}^{\mathrm{2}}& {\displaystyle }=\text{Var}\left[{\mathit{\tau }}_{\mathrm{3}}\right]-{\displaystyle \frac{\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{3}}]\text{Cov}[{\mathit{\tau }}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{3}}]}{\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{2}}]}}.\end{array}$$

McColl et al. (2014) extended the theory by solving the correlation coefficient of each dataset with respect to the ground truth, resulting in

$$\begin{array}{ll}\text{(10)}& {\displaystyle }{r}_{\text{RSP}}& {\displaystyle }=\pm \sqrt{{\displaystyle \frac{\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{2}}]\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{3}}]}{\text{Var}\left[{\mathit{\tau }}_{\mathrm{1}}\right]\text{Cov}[{\mathit{\tau }}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{3}}]}}},{\displaystyle }{r}_{\text{HSRL-2}}& {\displaystyle }=\pm \text{sgn}\left(\text{Cov}\right[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{3}}\left]\text{Cov}\right[{\mathit{\tau }}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{3}}\left]\right)\\ \text{(11)}& {\displaystyle }& {\displaystyle }\times \sqrt{{\displaystyle \frac{\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{2}}]\text{Cov}[{\mathit{\tau }}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{3}}]}{\text{Var}\left[{\mathit{\tau }}_{\mathrm{2}}\right]\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{3}}]}}},{\displaystyle }{r}_{\text{MODIS}}& {\displaystyle }=\pm \text{sgn}\left(\text{Cov}\right[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{2}}\left]\text{Cov}\right[{\mathit{\tau }}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{3}}\left]\right)\\ \text{(12)}& {\displaystyle }& {\displaystyle }\times \sqrt{{\displaystyle \frac{\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{3}}]\text{Cov}[{\mathit{\tau }}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{3}}]}{\text{Var}\left[{\mathit{\tau }}_{\mathrm{3}}\right]\text{Cov}[{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{2}}]}}},\end{array}$$

where sgn is the signum function; there is a sign ambiguity for each r_i, but practically all the AOD datasets are expected to be positively correlated with the unknown true value. If the ground truth is known, the correlation coefficient can be easily computed from Eq. (1). However, if the ground truth is unknown, the TC method enables the computation of the correlation coefficient from Eqs. (10)–(12).

Table 2 Mean squared deviation (MSD) decomposition of RSP and HSRL-2 AOD.

Table 3 Triple collocation analysis of RSP, HSRL-2, and MODIS AOD.

3.3 Data alignment

To ensure the robustness of triple collocation analysis, we attempt to minimize the representation error due to the different support scales between the AOD datasets. The 10 s HSRL-2 data are coarsened to 30 s resolution by horizontally averaging three consecutive data points, which results in a 3 km product. Then, all RSP data points within the same 30 s window are averaged. To minimize cloud impact, all HSRL-2 and RSP data points must be valid within the same window. The data alignment and aggregation result in 6988 pairs of data from 134 flights.

4 Results and discussion

4.1 Comparison of RSP and HSRL-2

The design of ACTIVATE to take place over the northwest Atlantic across different seasons was partly to characterize aerosol–cloud interactions across a wide range of aerosol conditions (Sorooshian et al., 2019). The overall distribution of RSP and HSRL-2 AOD is shown in Fig. 2, which demonstrates that a broad range of values was observed. The mean values of RSP and HSRL-2 AOD from 6988 pairs for the whole field campaign period are 0.173 and 0.120, respectively. RSP has a higher standard deviation than HSRL-2 (0.109 vs. 0.074). Both the mean and standard deviation are sensitive to outliers, so the median and interquartile range are also computed. The mean and median differences between RSP and HSRL-2 are about the same (0.053 vs. 0.047) and the interquartile range of RSP is larger than HSRL-2.

Table 1 shows the monthly breakdown of the AOD statistics. More than 80 % of data points come from March, May, and June. The largest mean AOD happens in May, August, and September for both instruments, and the median AOD corroborates the corresponding maxima. The same months also exhibit the largest variations as confirmed from the values of both standard deviation and interquartile range. The mean AOD is smallest in January and November, although the number of samples in winter is somewhat inadequate. The seasonal difference is largely consistent with Corral et al. (2021), who used coastal monitoring sites and a reanalysis product to study the regional AOD properties. During summer, there are more African dust events (Aldhaif et al., 2020) and more biomass burning events originating from North America (Mardi et al., 2021).

Figure 4 Two-dimensional histograms of the triple-collocated AOD datasets. Bin sizes are 0.02 for both instruments. (a) RSP vs HSRL-2. (b) MODIS vs HSRL-2. (c) MODIS vs RSP. In each plot, the linear fit (dashed black), 1:1 line (thick dashed gray), ±0.1 intervals (thin dashed gray), and selected performance metrics are shown.

Figure 5 Histograms of the MODIS AOD dataset. (a) Quality flag. (b) Cloud fraction. (c) Local coefficient of variation. Red color represents all the triple-collocated MODIS data points. Blue color represents all MODIS data points within a 25 km radius of RSP/HSRL-2.

The discrepancy between RSP and HSRL-2 AOD is shown in Fig. 3. Around 80 % of data points are within ± 0.1 intervals (Fig. 3a). The correlation coefficient r, RMSD, and MB are 0.727, 0.092, and 0.053, respectively. Most of the extreme outliers occur during the summer deployments (Fig. 3b–c).

A complementary approach of analyzing the deviation between RSP and HSRL-2 is to examine their mean squared deviation (MSD=RMSD²) decomposition. If two datasets x and y are identical to each other, all data points fall onto the 1 : 1 line (y=x, thick dashed gray lines in Fig. 3), and MSD equals zero. If there are discrepancies between two datasets (i.e., MSD>0), a linear fit $y=a+bx$ (dashed black lines in Fig. 3), where a is intercept and b is slope, quantifies the deviations from 1 : 1 line with a≠0 and/or b≠1. Here we use the MSD decomposition proposed by Gauch et al. (2003),

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\text{MSD}& =\text{SB}+\text{NU}+\text{LC}\\ & =(\stackrel{\mathrm{‾}}{y}-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+(\mathrm{1}-b{)}^{\mathrm{2}}{\mathit{\sigma }}_x^{\mathrm{2}}+(\mathrm{1}-r^{\mathrm{2}}){\mathit{\sigma }}_y^{\mathrm{2}},\end{array}\end{array}$$

where SB is squared bias, NU is non-unity slope, LC is lack of correlation, overlines denote the mean values, σ is standard deviation, r is the correlation coefficient between x and y, and b is the slope of the regression line of y on x.

Such partitioning has a straightforward geometric meaning for all three components. SB represents the translation of the linear fit from 1 : 1 line when there is a systematic bias, NU represents the rotation of linear fit when the slope b becomes different from 1, and LC represents the scattering of data points. The MSD decomposition in Table 2 shows that data scattering is the main contribution to the MSD for both winter and summer deployments. Systematic data bias is more severe in winter than in summer although the magnitude of bias is smaller in winter.

A question naturally arises from the above analysis as to which dataset is more accurate. Traditional metrics cannot answer this question because there is no available reference dataset over the western North Atlantic Ocean. To overcome this issue, we are motivated to carry out the triple collocation analysis.

Figure 6 Sensitivity of the performance of RSP filter criterion candidates. (a) Maximum normalized cost function (left column). (b) Maximum coarse-mode AOD (middle column). (c) Maximum fine-mode AOD (right column). In each plot, we calculate a performance metric as a function of criterion candidate threshold. As a result, the number of data points decreases when the criterion gets more stringent (smaller value). Intersecting dashed red lines indicate 50 % of the data points (y axis) meet the criterion threshold (x axis). First to fourth rows show, in order, the number of valid data points, mean bias, root mean squared deviation, and correlation coefficient.

Figure 7 Same as Fig. 3 but for HSRL-2 and filtered RSP data.

Table 4 Same as Table 2 but for HSRL-2 and filtered RSP AOD.

4.2 Triple collocation analysis

We introduce MODIS AOD as the third dataset. For each pair of RSP and HSRL-2 data, we collocate the closest MODIS data point within ± 60 min and 25 km, which is deemed suitable for mesoscale aerosol variabilities (Anderson et al., 2003). The collocation process results in 2344 triplets, roughly 34 % of the original data points. Zwieback et al. (2012) suggested that at least 500 triplets are needed for robust error estimates.

Table 3 summarizes the error metrics of the three instruments. HSRL-2 has the lowest error standard deviation (${\mathit{\sigma }}_{{\mathit{ϵ}}_{\text{HSRL-2}}}=\mathrm{0.0273}$) and highest correlation coefficient (r_HSRL-2=0.926) with respect to the true underlying values. After triple collocation, data scattering not only prevails between RSP and HSRL-2 (Fig. 4a) but also occurs between RSP and MODIS (Fig. 4c). HSRL-2 also maintains relatively higher correlation coefficients with both RSP and MODIS (Fig. 4a–b). The results suggest that HSRL-2 is the most accurate dataset among the three. Note that MODIS AOD corresponds to the whole atmospheric column, whereas RSP and HSRL-2 AOD correspond to the column below the aircraft. Therefore, MODIS AOD should be slightly larger than RSP and HSRL-2 AOD due to the contribution above the aircraft. The mean bias shows that MODIS is larger than HSRL-2 (Fig. 4b), but RSP is larger than MODIS (Fig. 4c). This suggests that it is more likely that RSP has a high bias over HSRL-2, not that HSRL-2 has a low bias over RSP (Fig. 3).

The result of the triple collocation analysis can be affected by a multitude of factors. Several factors related to MODIS data quality are discussed as follows.

First, MODIS provides a data quality flag which ranges from 0 to 3, with 3 being the best quality. Levy et al. (2013) recommend using data with non-zero quality flags. Figure 5a shows that all triple-collocated data points are at least 1 (red), which is much higher than all MODIS data points within a 25 km radius of RSP/HSRL-2 data points (blue). A related quantity is MODIS cloud fraction, which indicates possible cloud contamination. Figure 5b shows that the cloud fraction of over 55 % of triple-collocated data points is less than 0.1. The influence of excluding inferior MODIS collocated points is summarized in Table 3. While the error metrics of MODIS improve, there is little change for RSP and HSRL-2 which does not justify giving up 30 %–40 % of collocated data.

Second, since the Terra and Aqua satellites together provide only two overpasses during daytime because of their sun-synchronous orbits, we combine both datasets to maximize the number of files for collocation. This combination approach was adopted and tested in Antuña Marrero et al. (2018) by comparing 7-year Aqua, Terra, and Aqua + Terra data separately with measurements from an AERONET site. Generally, RMSD, MB, and r do not differ much in all cases (Table 5 in Antuña Marrero et al., 2018), which gives us confidence to combine Terra and Aqua data. We also performed sensitivity tests on the data quality of the third dataset and on the potential bias between two satellites using a synthetic data approach, which is included in the Supplement Sect. S1. The major conclusion is that the robustness of the error metrics is more sensitive to the number of data points rather than the data quality.

Third, collocation with the nearest MODIS data point may introduce bias when there is a spatial gradient of aerosols. Although it is generally difficult to evaluate the influence of the spatial gradient of aerosols, we adapt the local coefficient of variation (LCOV) from Anderson et al. (2003), which is the ratio of standard deviation and mean of all MODIS data points within a certain radius R of collocated RSP/HSRL-2 data points. In this case, R=25 km. Figure 5c shows that around 50 % and 90 % of LCOV are less than 0.1 and 0.2, respectively. An LCOV of 0.1 means that for a normally distributed random variable, 68 % of data points are within 10 % of the mean (Anderson et al., 2003). Similarly, 68 % of data points are within 20 % of the mean when the LCOV is 0.2. The mean HSRL-2 AOD is ∼ 0.1; therefore, most of the local AOD variation is estimated to be within 0.02, which is smaller than the expected error of the Dark Target algorithm. Interestingly, the current result is similar to the LCOV distributions obtained in Anderson et al. (2003) at similar spatial scale. Table 3 shows that including MODIS data points at high LCOV environs does not affect the assessment on the error metrics of RSP and HSRL-2.

Figure 8 Cloud impact on the bias between RSP and HSRL-2. (a) Bias-grouped by above-aircraft cloud mask fraction. (b) Bias-grouped by below-aircraft cloud mask fraction. (c) Cloud cover grouped by below-aircraft cloud mask fraction. Numbers at the top indicate the number of each data group. The box and whiskers indicates the variability of each group (box, first and third quartiles; black line in box, median; red line in box, mean; whiskers, a distance of 1.5 × IQR beyond the first and third quartiles; red circles, outliers).

4.3 RSP filter

The triple collocation analysis indicates that the most accurate AOD dataset over the western North Atlantic Ocean is HSRL-2 and paves the way for improving the data quality of RSP: the aim is to design a simple RSP filter to reduce disagreement between the two instruments. Three guiding principles are set. First, the filter should effectively get rid of outliers as suggested from the result of MSD decomposition (Table 2). Second, at the same time, the filter should keep as many data points as possible. Third, the criterion should be independent of HSRL-2 properties.

Three filter candidates are sought, including RSP normalized cost function χ^′, coarse-mode AOD τ_c, and fine-mode AOD τ_f. The normalized cost function χ^′ quantifies the performance of the vector radiative transfer model in the RSP-MAPP algorithm (Stamnes et al., 2018). Although it is tempting to use the total RSP AOD to filter out outliers, it is anticipated to be only effective for summer data because on average AOD is higher and fluctuates more during summer deployments (Table 1). Instead, τ_c and τ_f are chosen as they are part of retrieved parameters. For a successful retrieval, other than χ^′ being lower than 0.15, the algorithm also requires that $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\le {\mathit{\tau }}_{\mathrm{c}}\le \mathrm{0.3}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\le {\mathit{\tau }}_{\mathrm{f}}\le \mathrm{0.7}$. Because τ_c and τ_f are rarely larger than 0.2 and 0.5, respectively, lowering either may help eliminate outliers.

The filter candidates are tested by iteratively decreasing the maximum value allowed in the dataset, and three performance metrics (r, MB, and RMSD) are evaluated (Fig. 6). For a balanced candidate, it is expected that when the criterion gets more stringent, r should increase, and MB and RMSD should decrease at the same time. It is clear that the fine-mode AOD does not qualify because r generally decreases, while the maximum allowed fine-mode AOD is reduced (Fig. 6c). In other words, constraining the fine-mode AOD will lead to more data outliers, which goes against one of our principles. The remaining cost function and coarse-mode AOD are balanced candidates. A relatively fair way to compare the two balanced candidates is to focus on their performance metrics when 50 % of data points are left (intersection of dashed red lines in Fig. 6a–b). That means the ranges of cost function and coarse-mode AOD are 0 to 0.049 and 0 to 0.014, respectively. Compared to the coarse-mode AOD, the cost function performs better on the whole. Therefore, the cost function is regarded as the best simple filtering criterion. It has two additional benefits. First, it is globally applicable. Second, it also allows offline post-processing. For simplicity, the cost function threshold is set to be 0.05 thereupon.

Figure 9 Same as Figure 8 but for aerosol impact. (a) Average HSRL-2 aerosol column total thickness grouped by HSRL-2 aerosol type. (b) Bias-grouped by the number of HSRL-2 aerosol types. (c) Bias-grouped by RSP fine-mode single-scattering albedo (SSA) at 532 nm. (d) Bias-grouped by the HSRL-2 column-maximum 532 nm depolarization ratio δ₅₃₂. (e) The maximum HSRL-2 δ₅₃₂ at the lowest 500 m grouped by the collocated MERRA-2 relative humidity (RH). (f) Bias-grouped by the difference of aircraft altitude and HSRL-2 (effective) aerosol top height (ATH). (g) Bias-grouped by HSRL-2 mixed layer height (MLH) during summer deployments. (h) Bias-grouped by HSRL-2 MLH during winter deployments. Numbers on the right indicate the number of each data group.

To evaluate the cost function filter, we perform the comparison between RSP and HSRL-2 again (Fig. 7). All metrics generally improve. Compared to 80 % in Fig. 3a, more than 89.8 % of data points are now within ± 0.1 intervals after filtering. By analyzing MSD decomposition before and after filtering (Tables 2 and 4), it is found that the filter is effective in reducing MSD during summer deployments by cutting down around 80 % of the systematic bias.

4.4 Why do RSP and HSRL-2 disagree with each other?

Stamnes et al. (2018) suggested some possible sources of the discrepancies between two instruments: (i) aerosols above the aircraft, (ii) clouds above the aircraft, (iii) aerosols near the aircraft, (iv) three-dimensional effects due to clouds surrounding the aircraft, (v) aerosols close to the surface (within 100 m), (vi) multiple aerosol layers from different types of aerosol, (vii) presence of absorbing aerosols, (viii) presence of non-spherical aerosols, and (ix) intrinsic retrieval and measurement uncertainties. These factors are examined using the available field campaign data.

To streamline the analysis, the aforementioned factors can be broadly organized into two groups: cloud and aerosol impacts. We use the above-aircraft and below-aircraft cloud masks to detect the presence of cloud that leads to the contamination of aerosol retrievals. Cloud mask is a binary variable being either zero (absence) or one (presence). Because the temporal resolution of cloud mask is 1 Hz, we compute the fraction of cloud mask being 1 (i.e., cloud is present) during each 30 s window of a 3 km HSRL-2 data point. There is only 1 year of below-aircraft cloud mask data for 2021 summer and winter deployments, but it includes other cloud properties such as cloud fraction. The cloud fraction is computed as the number of cloud pixels divided by the total number of pixels on a given camera image which does not have sun glint.

Figure 10 Schematic of altitudes associated with the high-flying aircraft King Air and its onboard instruments RSP and HSRL-2.

Above the aircraft, over 65 % (4546 out of 6988) of the cloud mask fraction are greater than zero (Fig. 8a). It is interesting that once the fraction is above zero, the mean and median AOD biases are almost invariant, which suggests that the cloud influence on aerosol retrievals is largely independent of its cloud mask fraction. Most of the outliers occur when the cloud mask is on for over 80 % of time. Below the aircraft, around 34 % (724 out of 2156) of the cloud mask fraction is greater than zero (Fig. 8b). The mean and median biases are again similar across a wide range of cloud mask fraction groups. It is not necessary that the presence of cloud is related to high cloud cover fraction, but higher cloud mask fraction is associated with higher below-aircraft cloud cover (Fig. 8c). This implies that there is not a linear relationship between the cloud cover and bias. Note that RSP normalized cost function χ^′ is sensitive to below-aircraft cloud mask fraction but not above-aircraft cloud mask fraction (not shown).

For aerosol impact, first it is helpful to be aware of the diversity of aerosols encountered during the field campaign because all eight aerosol types are detected using the HSRL aerosol classification scheme (Fig. 9a). Pure dust and ice are scarcely identified, but all other types occur around 50 %–60 % of time. The total bin thickness for each aerosol type in each column is estimated by the product of the number of bin associated with that aerosol type and the thickness of each bin (∼ 15 m). The bin thickness among different aerosol types varies a lot. The whiskers of urban aerosol and the dusty mix are more than 2500 m. The median of other aerosols is usually below 200 m. The vertical aerosol variability can be very complicated, so we computed the number of aerosol types in each column, which can be viewed as the minimum number of aerosol layers in each column. Over 73 % of time there are more than two aerosol types (Fig. 9b), which can potentially invalidate one of the RSP aerosol model assumptions. Note that the zero type occurs because the classification scheme occasionally fails to distinguish and identify it as an ambiguous type which is not counted. The median bias of all groups with different number of aerosol type is similar, but the bias gets more dispersed when the number of aerosol type increases.

Figure 11 Annual cycle of CALIPSO stratospheric AOD centered at 37.5° N and 70° W. The 2006–2021 climatology is drawn in black, and all available ACTIVATE years are drawn in other colors.

The presence of absorbing aerosols may be indicated by the presence of four HSRL-2 aerosol types: smoke, fresh smoke, pure dust, and the dusty mix (Fig. 9a). The dusty mix is a mixture of pure dust with other types of aerosols (Burton et al., 2012). Another way is to examine the RSP fine-mode single-scattering albedo (SSA) at 532 nm, which is defined as the ratio of scattering coefficient and extinction coefficient. The larger the absorption efficiency, the more SSA departs from 1. An SSA value from 0.90 to almost 1 represents an aerosol transitioning from moderately to weakly absorbing (Stamnes et al., 2018). It is clear that the minimum mean bias occurs when SSA is close to 1, and the mean bias gradually increases when SSA decreases (Fig. 9c). Around 47 % of SSA is below 0.90, which gives a similar result compared to HSRL-2 absorbing aerosols (the dusty mix accounts for around 53 % of time).

Figure 12 Time series of selected measurements during ACTIVATE RF67 on 19 May 2021. (a) HSRL-2 (red), RSP total (blue), RSP coarse-mode (orange) AOD, and RSP normalized cost function (gray). (b)–(c) King Air (red) and Falcon (black) altitude, RSP aerosol top height (purple), and HSRL-2 mixing layer height (gray). All altitude values are in kilometers. (b) HSRL-2 532 nm aerosol scattering coefficient β_a,532 in Mm⁻¹ sr⁻¹ (color scale). (c) HSRL-2 532 nm depolarization ratio δ₅₃₂ (color scale).

The presence of non-spherical aerosols may be indicated by the presence of three HSRL-2 aerosol types: ice, pure dust, and the dusty mix (Fig. 9a). Generally, the linear depolarization ratio is zero for spherical particles but markedly departs from zero for non-spherical particles (Mishchenko and Hovenier, 1995). Therefore, all these aerosol types are characterized by a high degree of aerosol depolarization, although there is not a clear-cut threshold for these three aerosol types. Schlosser et al. (2022) deemed 0.13 a suitable threshold for the maritime study region of ACTIVATE. Burton et al. (2013) found that the first quartiles of depolarization for ice, pure dust, and the dusty mix are 0.23, 0.31, and 0.13, respectively, from measurements collected on 109 flights during 2006 and 2012. The HSRL-2 column-maximum depolarization is more than 0.1 for over 84 % of measurements (Fig. 9d). However, unless the column-maximum depolarization is less than 0.05, all other data groups show a large bias range.

As part of marine aerosols, sea salt is another potential source for contributing to high depolarization in the marine boundary layer due to its shape. Pure dry sea salt is non-spherical, but its shape is sensitive to ambient relative humidity (RH). When RH increases, the sea salt particle remains solid until it reaches the deliquescence RH at which the solid particle absorbs water and forms aqueous solution. The aqueous-phase sea salt is nearly spherical and therefore has low depolarization. On the contrary, when RH decreases, the aqueous sea salt solution can return to its crystalline form at the efflorescence RH. The deliquescence and efflorescence RHs are usually different. Therefore, the shape of a sea salt particle depends not only on the magnitude but also the history of ambient RH (Haarig et al., 2017).

Figure 13 Time series of selected measurements during ACTIVATE RF143 on 22 March 2022. (a) HSRL-2 (red), RSP total (blue), RSP coarse-mode (orange), MODIS (green) AOD, and RSP normalized cost function (gray). (b)–(c) King Air (red) and Falcon (black) altitude, RSP aerosol top height (purple), and HSRL-2 mixing layer height (gray). All altitudes are in kilometers. (b) HSRL-2 532 nm aerosol scattering coefficient β_a,532 in Mm⁻¹ sr⁻¹ (color scale). (c) HSRL-2 532 nm depolarization ratio δ₅₃₂ (color scale). (d) Lowest 500 m MERRA-2 (blue) and Falcon (black) relative humidity in percent.

Ferrare et al. (2023) reported that sea salt may cause high depolarization at low RH (i.e., more likely to be non-spherical) within the boundary layer for around one-third of flights during the first 2 years of ACTIVATE based on HSRL-2 and dropsonde measurements. Here we extend the result by collocating the maximum HSRL-2 depolarization at the lowest 500 m with the interpolated MERRA-2 RH for the whole campaign period (Fig. 9e). Seethala et al. (2021) showed that despite having a rather coarse spatial resolution, MERRA-2 adequately reconstructs the observed lower tropospheric thermodynamic structure during ACTIVATE, especially the relative humidity and specific humidity fields. Median depolarization gradually increases when RH decreases, suggesting that sea salt may change its shape over a wide range of RH due to both deliquescence and efflorescence. Note that low RH conditions (<60 %) only occur for around 25 % of time, which means non-spherical aerosols do not often play a big role in the RSP retrieval errors.

The high depolarization associated with sea salt provides evidence of the existence of aerosols close to the surface because sea surface is the major source of sea salt in marine boundary layer. Conversely, it is more difficult to detect aerosols near and above the aircraft. The near-range limit for all ACTIVATE HSRL-2 extinction profiles is 1500 m so that the transmitter-to-receiver overlap function is close to 1 (Burton et al., 2012; Hair et al., 2008). To detect the possibility of aerosols within the near range of aircraft, we compute an effective HSRL-2 aerosol top height h_t, which is analogous to the RSP MAPP aerosol top height (Stamnes et al., 2018):

$$\begin{array}{}\text{(14)}& {\displaystyle }\underset{\mathrm{0}}{\overset{h_t}{\int }}\mathit{\alpha }\mathrm{d}z=\mathrm{0.95}\underset{\mathrm{0}}{\overset{h}{\int }}\mathit{\alpha }\mathrm{d}z,\end{array}$$

where h_t is aerosol top height, h is aircraft altitude, and α is extinction coefficient.

The aerosol top height difference Δh_t can then be defined as the difference between the aircraft altitude and effective aerosol top height (Fig. 10). Obviously, the difference is at least 1500 m due to the near-range limit (Fig. 9f). A small difference indicates that the aerosol profile is top heavy, but the result for 1500–2000 m does not show larger bias compared to other groups. Note that the dispersion of bias is much smaller when the effective top height is very low.

Aerosols above the aircraft can be estimated from the CALIPSO stratospheric AOD (Fig. 11). Climatologically, the mean stratospheric AOD in the study region is around 0.011 and swings between 0.008 and 0.013 throughout a year. During ACTIVATE, the 2020 winter deployment experienced the record-high AOD in February and March, and the 2021 summer deployment experienced the record-high AOD in May for all CALIPSO years (2006–2021). Major possible sources of stratospheric AOD perturbation in midlatitudes include volcanic eruption and deep convection (Kremser et al., 2016).

We speculate that the large stratospheric AOD perturbations are more likely due to volcanic eruption, not deep convection. First, overshooting convection climatologically occurs more often over land than ocean in midlatitudes (Liu et al., 2020). Second, most of the overshooting convection in the 30–40° latitude band occurs during MAM and JJA (Liu et al., 2020), which is opposite to the CALIPSO climatology, i.e., maxima during SON and DJF but minima during MAM and JJA. The peak AOD in 2020 may be due to the stratospheric eruption at Raikoke (48° N, 153° E) in June 2019 (Kloss et al., 2021). Its estimated plume top height is about 16–17 km from the Himawari-8 satellite. Subsequently the stratospheric AOD increased up to 0.027 in October 2019 (not shown) and gradually decreased in the ensuing months. The spike in May 2021 may be due to the eruption at La Soufrière (13° N, 61° W) in April. Although there are large stratospheric AOD perturbations during ACTIVATE, their magnitude (∼ 0.01–0.02) cannot explain all the biases between the datasets in this study. Therefore, the influence of aerosols above the aircraft should be relatively small compared to other factors.

Figure 14 Time series of selected measurements during ACTIVATE RF170 on 10 June 2022. (a) HSRL-2 (red), RSP total (blue), RSP coarse-mode (orange), MODIS (green) AOD, and RSP normalized cost function (gray). (b)–(c) King Air (red) and Falcon (black) altitude, RSP aerosol top height (purple), and HSRL-2 mixing layer height (gray). All altitudes are in kilometers. (b) HSRL-2 532 nm aerosol scattering coefficient β_a,532 in Mm⁻¹ sr⁻¹ (color scale). (c) HSRL-2 532 nm depolarization ratio δ₅₃₂ (color scale). (d) Lowest 500 m MERRA-2 (blue) and Falcon (black) relative humidity in percent.

Other than the single homogeneous aerosol layer assumption, the assumption of coarse-mode sea salt aerosol being located between the surface and 1 km may also lead to retrieval uncertainties. During summer deployments, the mean HSRL-2 mixed layer height (MLH) is around 534 m. The bias is much larger when the RSP coarse-mode aerosol layer height (assumed to be 1 km) is much higher than the HSRL-2 MLH (Fig. 9g). During winter deployments, the mean HSRL-2 MLH (around 764 m) is much closer to the RSP coarse-mode aerosol layer height, but the large bias prevails when there is big mismatch between two heights (Fig. 9h). Note that on average the coarse-mode RSP AOD accounts around 14 % of total RSP AOD so the overall influence is likely to be tiny.

It is beneficial to consolidate the result by revisiting all factors. (i) There is no direct evidence for the presence of above-aircraft aerosols, but elevated stratospheric AOD is seen at least in 2020 and 2021, though the overall impact should be small. (iii) There is no direct evidence for aerosols near aircraft. (ii) and (iv)–(viii) There is evidence for the presence of above-aircraft clouds, clouds surrounding the aircraft, sea salt near the sea surface, more than two aerosol layers, absorbing aerosols, and non-spherical aerosols. (ix) There is evidence for simplified and unwarranted assumptions, but the impact may be small.

4.5 Case studies

Three case studies are arranged to complement the previous analysis. All three cases are statistical surveys in which Falcon collects data at different levels of marine boundary layer. Full details of the Falcon flight strategy and rationale for statistical survey flights are provided in Dadashazar et al. (2022) and Sorooshian et al. (2023).

The first case study is research flight (RF) 67 on 19 May 2021 (King Air: 12:31–15:55 UTC; Falcon: 12:27–15:49 UTC). The two aircraft took off from NASA LaRC, headed southeast via OXANA, and returned near 33° N. RF67 provides a baseline for good agreement between RSP and HSRL-2 AOD. Except near the turning point and towards the end of the trip, the whole survey is mostly cloud free and under light wind.

Generally, the aerosol loading is low, and the difference between RSP and HSRL-2 AOD is small (Fig. 12a). The RSP cost function is small most of the time. The resemblance between RSP and HSRL-2 can be attributed to the relative simple aerosol vertical distribution. Most of the RSP AOD comes from the fine mode (Fig. 12a). The RSP aerosol top height stays at around 3 km (Fig. 12b). The HSRL-2 532 nm backscattering coefficient signal mostly concentrates under the mixing layer at 1 km. Also, in the mixing layer, the 532 nm depolarization ratio is low, indicating a small influence from sea salt (Fig. 12c).

The second case study is RF143 on 22 March 2022 (King Air: 17:12–21:00 UTC; Falcon: 17:36–21:12 UTC). The afternoon flight constitutes the return trip from Bermuda to LaRC. While the morning flight (RF142) from LaRC to Bermuda observes a pronounced aerosol gradient, the afternoon homebound flight (RF143) provides a more complicated aerosol structure.

Overall, the aerosol loading is above the monthly average (Fig. 13a). The RSP aerosol top height is higher compared to RF67 with a coupled aerosol layer above the mixing layer (Fig. 13b). The aerosol layer becomes decoupled around 18:30 UTC. Initially, RSP and HSRL-2 AOD agree with each other, which can also be confirmed by the collocated MODIS data. Starting from 19:00 UTC, enhanced aerosol scattering is observed in the boundary layer for around 45 min (Fig. 13b). At the same time, RSP AOD fluctuates a lot with multiple spikes from its coarse mode (Fig. 13a). The rapidly changing RSP AOD is most likely due to both aerosol layers in the mixing layer and near 4 km aloft. Both layers do have a high depolarization ratio which indicates non-spherical aerosols (Fig. 13c). There are not enough MERRA-2 relativity humidity data collocated with HSRL-2 at the lowest 500 m (Fig. 13d), so we make use of in situ relative humidity observations from Falcon near 19:20 UTC (highlighted in gray shade). The high depolarization ratio coincident with relatively low relative humidity suggests that it might be related to sea salt. The HSRL aerosol classification scheme indicates that the decoupled aerosol layer consists of ice and the dusty mix, which partly explains the observed high depolarization ratio.

The final case study is RF170 on 10 June 2022 (King Air: 11:57–15:35 UTC; Falcon: 12:20–15:37 UTC). By that time, the flight base had moved to Bermuda. The two aircraft took off and headed southwest for a statistical survey.

RSP and HSRL-2 AOD at times diverge from each other (Fig. 14a). The occurrence of extreme differences is mostly associated with spikes in RSP coarse-mode AOD. Most of the backscattering signal is below 2 km (Fig. 14b), and the mixing layer is very shallow (around 500 m). A remarkably high depolarization ratio is seen above the mixing layer near 13:45 UTC (Fig. 14c). The high MERRA-2 and Falcon near-surface relative humidity indicates that it is not associated with sea salt (Fig. 14d). The HSRL aerosol classification scheme identifies that most of the high depolarization aerosols between 13:00 and 14:30 UTC consist of a dusty mix, suggesting some seasonal influence from African dusts.

In short, the RSP and HSRL-2 show good agreement in relatively clean environments in the first case study, but as demonstrated by the second and third cases, the RSP retrieval degrades under a high polarization ratio, which indicates the influence of non-spherical aerosols.

5 Conclusions

We analyzed 3 years of AOD data over the western North Atlantic Ocean and attempted to estimate and understand their uncertainties, using two retrieval products collected during ACTIVATE and one satellite retrieval product. The key findings are as follows:

1. RSP and HSRL-2 AOD have large seasonal variations but also exhibit considerable deviations between the two retrievals. The absence of a reference dataset over the western North Atlantic Ocean, and arguably all global marine areas, prompts us to use triple collocation to determine which dataset is more accurate.

2. Triple collocation analysis indicates that HSRL-2 is the most accurate dataset over the study region. The correlation coefficient and error standard deviation of AOD with respect to the expected ground truth are 0.93 and 0.027, respectively.

3. The RSP retrieval quality can be improved considerably by applying a simple, yet efficient, filtering criterion: a more stringent normalized cost function of 0.05. There are two additional benefits. First, it can be easily applied globally. Second, it is suitable for offline post-processing.

4. The reasons of disagreement between RSP and HSRL-2 are statistically diagnosed using all ACTIVATE datasets. There is a basket of potential factors, including cloud contamination, aerosols near the surface, multiple aerosol layers, absorbing aerosols, non-spherical aerosols, and simplified retrieval assumptions. The statistical analysis is accompanied with case studies of both subpar and high-quality RSP retrievals. By combining information from lidar, polarimeter, and other pertinent datasets, it is possible to single out the most probable causes of poor retrievals, such that RSP AOD retrievals can be optimally improved.

Motivated by the last point, future efforts are needed along two lines: multi-angle polarimeter and HSRL (or HSRL-2) data from other field campaigns should be analyzed to test the robustness of the conclusions here, and a more efficient retrieval of HSRL-2 MLH data should be developed potentially through two-dimensional image processing. Earlier recommendations such as combining HSRL-2 elevated aerosol layer altitudes with RSP retrievals should be reconsidered (Knobelspiesse et al., 2011). The current HSRL-2 MLH retrieval algorithm detects the edge of aerosol gradient in each one-dimensional column. Instead, reading the whole two-dimensional flight curtain at once and using a two-dimensional wavelet transform will potentially speed up the computation and improve the detection accuracy. The above results will provide a better constraint on aerosol properties and help prepare concurrent deployments of lidars and polarimeters for future spaceborne missions.