# A Combined Lidar-Polarimeter Inversion Approach for Aerosol Remote Sensing Over Ocean

Feng Xu, Lan Gao, Jens Redemann, Connor Flynn, W. Reed Espinosa, Arlindo da Silva, Snorre Stamnes, Sharon Burton, Xu Liu, Richard Ferrare, et al.

► To cite this version:

Feng Xu, Lan Gao, Jens Redemann, Connor Flynn, W. Reed Espinosa, et al.. A Combined Lidar-Polarimeter Inversion Approach for Aerosol Remote Sensing Over Ocean. Frontiers in Remote Sensing, 2021, 2, 10.3389/frsen.2021.620871. [https://hal.science/hal-03320478](https://hal.science/hal-03320478)

HAL Id: hal-03320478

https://hal.science/hal-03320478

Submitted on 16 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# A Combined Lidar-Polarimeter Inversion Approach for Aerosol Remote Sensing Over Ocean

Feng Xu¹*, Lan Gao¹, Jens Redemann¹, Connor J. Flynn¹, W. Reed Espinosa², Arlindo M. da Silva², Snorre Stamnes³, Sharon P. Burton³, Xu Liu³, Richard Ferrare³, Brian Cairns⁴ and Oleg Dubovik⁵

¹School of Meteorology, The University of Oklahoma, Norman, OK, United States, ²NASA Goddard Space Flight Center, Greenbelt, MD, United States, ³NASA Langley Research Center, Hampton, VA, United States, ⁴NASA Goddard Institute for Space Studies, New York, NY, United States, ⁵Université de Lille, CNRS, UMR 8518 - LOA - Laboratoire d'Optique Atmosphérique, Lille, France

An optimization algorithm is developed to retrieve the vertical profiles of aerosol concentration, refractive index and size distribution, spherical particle fraction, as well as a set of ocean surface reflection properties. The retrieval uses a combined set of lidar and polarimeter measurements. Our inversion includes using 1) a hybrid radiative transfer (RT) model that combines the computational strengths of the Markov-chain and adding-doubling approaches in modeling polarized RT in vertically inhomogeneous and homogeneous media, respectively; 2) a bio-optical model that represents the water-leaving radiance as a function of chlorophyll-a concentration for open ocean; 3) the constraints regarding the smooth variations of several aerosol properties along altitude; and 4) an optimization scheme. We tested the retrieval using 50 sets of coincident lidar and polarimetric data acquired by NASA Langley airborne HSRL-2 and GISS RSP respectively during the ORACLES field campaign. The retrieved vertical profiles of aerosol single scattering albedo (SSA) and size distribution are compared to the reference data measured by University of Hawaii's HiGEAR instrumentation suite. At the vertical resolution of 315 m, the mean absolute difference (MAD) between retrieved and HiGEAR derived aerosol SSA is 0.028. And the MADs between retrieved and HiGEAR effective radius of aerosol size distribution are 0.012 and 0.377 micron for fine and coarse aerosols, respectively. The retrieved aerosol optical depth (AOD) above aircraft are compared to NASA Ames 4-STAR measurement. The MADs are found to be 0.010, 0.006, and 0.004 for AOD at 355, 532 and 1,064 nm, respectively.

**Keywords:** aerosols, remote sensing, lidar-polarimeter, combined retrieval, radiative transfer

## INTRODUCTION

Aerosols are tiny particles suspended in the atmosphere. They stem from a variety of natural and anthropogenic sources and processes, and influence the distribution of radiation by either directly absorbing, scattering and emitting radiation or indirectly by providing cloud condensation nuclei. In addition, the concentration and speciation of aerosols in the boundary layer impact air quality and consequently public health. Characterization of spatial and temporal variations of a few key aerosol characteristics such as abundance, size, absorption, and composition is critical for quantifying their

## OPEN ACCESS

**Edited by:**

Jose Antonio Sobrino,
University of Valencia, Spain

**Reviewed by:**

Pengwang Zhai,
University of Maryland, United States
Xiaoguang Xu,
University of Maryland, United States

*Correspondence:

Feng Xu
fengxu@ou.edu

**Specialty section:**

This article was submitted to
Satellite Missions,
a section of the journal
Frontiers in Remote Sensing

Received: 24 October 2020

Accepted: 07 January 2021

Published: 21 April 2021

**Citation:**

Xu F, Gao L, Redemann J, Flynn CJ, Espinosa WR, da Silva AM, Stamnes S, Burton SP, Liu X, Ferrare R, Cairns B and Dubovik O (2021) A Combined Lidar-Polarimeter Inversion Approach for Aerosol Remote Sensing Over Ocean.
Front. Remote Sens. 2:620871
doi: 10.3389/frsen.2021.620871
global and regional impact on climate and air quality. However, the short life cycle of aerosols and high-degree of variability in their composition and distribution has been a challenge for their accurate remote sensing. According to the most recent IPCC report (IPCC, 2013), aerosols remain the largest radiative forcing uncertainty.

To better anchor aerosol direct and indirect impacts on the Earth's energy balance, the Global Climate Observing System (GCOS, 2011) sets out observational requirements for a few key aerosol characteristics such as aerosol optical depth (AOD, denoted by "τ"), single scattering albedo (SSA), extinction profiles and layer height. For example, the accuracy of AOD at 550 nm is required to be better than the max of (0.1τ, 0.03) at the horizontal resolution of 5-10 km and a temporal resolution of 4 h. The accuracy of SSA is required to be better than 0.03. Aerosol extinction profiles are required to be better than 10% at a horizontal resolution of 200-500 km, a vertical resolution <1 km near the tropopause and of ~2 km in the middle stratosphere. The NASA Aerosol Clouds, and Convection and Precipitation (ACCP) mission that is currently undergoing formulation (https://vac.gsfc.nasa.gov/accp/home.htm) is also setting up criteria on a set of geophysical variables (documented as Milestone SATM, September 16, 2019). For example, the total and planetary boundary layer (PBL) effective total fine mode AODs are desired to be better than 0.05τ<sub>total/fine</sub> + 0.02. In both column and PBL effective sense, the accuracy of SSA is desired to be better than 0.03; aerosol extinction-to-backscatter ratios are required to be better than 25%; the real part of aerosol refractive index is desired to be better than 0.025; and the nonspherical AOD fractions are desired to be better than 10%. Aerosol extinction profile is desired to be better than the max of (0.02 km<sup>-1</sup>, 20%-40%) at the vertical resolution from 30 to 250 m, depending on specific mission objective.

Capitalizing on the rapid advancement of lidar and polarimetric remote sensing techniques and the relevant remote sensing theory over the last few decades, the combination of these two types of measurement is highly promising for meeting the requirements described above. As an active remote sensing technique, lidar sends out laser pulses into atmosphere and detects the backscattered radiation. Accounting for speed of light, the time delay in the received signals determines the distance from the backscatter location, and the range-resolved profiles of attenuated backscatter (e.g., CALIOP (Winker et al., 2007)) or backscatter, extinction and depolarization (e.g., the measurement by the High Spectral Resolution Lidar (HSRL, (Hair et al., 2001; Hair et al., 2008; Rogers et al., 2009) at single or multiple wavelengths resolves aerosol types and microphysical properties, and consequently vertically resolved radiative forcing. The determination of vertically resolved aerosol radiative forcing has direct implications for the thermodynamic structure of atmosphere and consequent effects on monsoonal circulations and the global hydrological cycle. Polarimeters, in contrast, use the Sun as the radiation source and measure radiometric and polarimetric signals of sunlight reflection from top-of-atmosphere (TOA) in multi-spectral/hyper-spectral and multi-angular dimensions. It enables a moderate-to-high horizontal resolution aerosol/cloud optical and microphysical products

during daytime. The global coverage of aerosol products enables its wide applications such as aerosol-cloud interaction and air quality studies. Like a radiometer, a polarimeter can provide large swath aerosol/cloud context for TOA flux assessment. Both theoretical studies and real data analysis have proved polarimeter's capability of constraining column effective aerosol abundance, absorption and microphysical properties together with surface reflection properties (see for example [Mishchenko et al., 2011; Hasekamp et al., 2011; Chowdhary et al., 2012; Knobelspiesse et al., 2012; Dubovik et al., 2019; Remer et al., 2019; among others]). Compared to lidar, however, polarimeter carries much less information about the vertical distributions of aerosols.

It is apparent that lidar and polarimeter have highly complementary strengths: while lidar resolves details of aerosol vertical profiles and types, polarimeter provides constraints on anchoring column-effective aerosol abundance, absorption and microphysical properties across a large swath. Consequently, a combination of these two types of measurements enhances the observational information about aerosol properties. Informing aerosol retrieval by these two types of measurements, we target a comprehensive determination of aerosol size, type, abundance and absorption as well as their vertical variations.

To fully utilize the information contained in lidar and polarimetric measurements, reliable forward and inversion models are prerequisite. Some earlier efforts have been made toward combining lidar and polarimetric/radiometric measurements to improve the retrieval of aerosol properties. For example, Knobelspiesse et al. (2011) investigated the Research Scanning Polarimeter (RSP, Cairns et al., 1999) and HSRL (Hair et al., 2008) data collected during the field campaign - Arctic Research of the Composition of the Troposphere from Aircraft and Satellites (ARCTAS). It is found that the combined use of the two types of measurement in retrieval reduces the likelihood of unsuccessful retrievals caused by significantly biased initial estimate of the aerosol number concentration in optimal estimate inversion approach. In the case of ground-based observations, the LIRIC and GARRLiC algorithms have been developed by Chaikovsky et al. (2016) and by Lopatin et al. (2013), respectively. The latter is currently a branch of the Generalized Retrieval of Aerosol and Surface Properties (GRASP) algorithm (Dubovik et al., 2011; Dubovik et al., 2014). These two algorithms use joint data from a multi-wavelength lidar and an AERONET radiometer and derive two vertical profiles of fine and coarse aerosol components as well as extra parameters on the column-integrated properties of aerosols. Particularly, LIRIC derives vertical profiles using the retrieval results from AERONET but not changing them. GARRLiC/GRASP, on the other hand, fits simultaneously lidar and photometer data and derives both enhanced vertical and columnar aerosol properties. As ongoing efforts, Liu et al., (2016) and Burton et al., (2019) have been developing an optical estimate approach to retrieve vertically resolved particle concentrations, effective radius, and absorption by using both RSP and HSRL observations acquired during DISCOVER-AQ, TCAP and ORACLES field missions. Espinosa et al., (2019) used the GRASP algorithm to simulate the retrieval of attenuated
backscatter lidar measurements in 532 and 1,064 nm in combination with a radiometer or polarimeter with 10 angles spaced over ±57° with 6 wavelengths spanning 0.44–2.2 µm. Compared to stand alone polarimetric or radiometric retrieval, an enhanced accuracy in fine mode aerosol SSA was observed from the combined retrieval. In essence, the combined retrieval capitalizes on the benefit of backscatter observations that contain sensitivity to the profiling properties of aerosols and the benefits of polarimetric observations that contain information about aerosol loading and absorption.

In this paper, we describe a research algorithm that allows a combined use of lidar and polarimetric observations for retrieving vertical profiles of speciated aerosol abundance and aerosol intrinsic properties such as size distribution, refractive index, and spherical particle fraction. The algorithm builds upon our earlier work using only polarimetric observations to retrieve a column effective set of aerosol properties over ocean (Xu et al., 2016). The paper is organized as follows. Following an algorithm overview in *General Structure of the Algorithm*, we introduce our 1D RT model for polarized RT in a coupled atmosphere-ocean system in *Aerosol Scattering and Ocean Surface Reflection. Forward Model for Light Propagation in a Coupled Atmosphere-Ocean System* formulates the inversion and error estimate approach. In *Coupled Retrieval of Aerosols and Ocean Properties*, we test the inversion algorithm using 50 sets of coincident lidar and polarimetric observations made by the HSRL-2 and RSP instruments aboard the NASA ER-2 during the NASA's ORACLES field campaign (Redemann et al., 2021). These observations cover a moderately high bio-mass burning and sea salt aerosol loading with column total AOD$_{532-nm}$ between ~0.2 and ~0.3 and an altitude dependent SSA. Comparisons of the retrieved AOD, SSA, and aerosol size distribution are made to the reference data measured by 4-STAR and HiGEAR which are aboard the NASA P-3 aircraft. A summary is given in *Retrieval Implementation and Validation*.

## GENERAL STRUCTURE OF THE ALGORITHM

The combined lidar and polarimetric retrieval of aerosol loading and properties consists of three major modules: 1) A static light scattering dataset for intrinsic aerosol properties including refractive index and multiple rectangular size bins. The dataset is pre-calculated and called during the inversion to provide on-the-fly single scattering properties for spherical and nonspherical particle properties. These properties include aerosol optical depth, phase matrix and SSA for radiative transfer modeling; and extinction, backscatter and depolarization to fit lidar measurement; 2) A hybrid radiative transfer model that combines the computational strengths of the Markov-chain and adding-doubling approaches in modeling polarized RT in vertically inhomogeneous and homogeneous media, respectively; and 3) Optimization model. The inversion is assisted by allowing an imposition of smoothness constraints on the vertical variations of speciated aerosol refractive index, size parameters and spherical particle fraction. Following retrieval success and the

output of the optimized solution, ancillary aerosol quantities including AOD, SSA, particulate lidar ratio, particulate extinction and backscatter coefficients as well as particulate depolarization are derived. Error estimate is provided for all retrieval outcome.

With these algorithm modules, we outlined in **Figure 1** the algorithm flow of combined lidar-polarimeter inversion. Input into the retrieval consists of 1) lidar signals including extinction, backscatter, and/or depolarization or a subset of them; and 2) polarimetric signals including multi-angle and multi-spectral measurements of radiance and polarimetric components. As retrieval parameters, the state vector contains the column volume concentrations of multiple aerosol species and vertical profiles of 1) speciated aerosol refractive index with its spectral dependence constrained by a linear model in logarithmic space; 2) speciated log-normal volume weighted size distribution constrained by median size parameter and standard deviation; 3) species dependent spherical particle volume fraction; and 4) species column concentrations in vertical layers. Also included in the retrieval state vector are chlorophyll-a (Chl-a) concentration, an empirical spectral term that adjusts the water-leaving radiances to mitigate bio-optical modeling error, and the wind speed.

Our retrieval starts with an initialized state vector $\mathbf{x}_0$. Speciated aerosol abundance, optical properties and microphysical properties are used to calculate vertically resolved aerosol single scattering quantities including 1) AOD, SSA and phase matrix for radiative transfer modeling and fit polarimetric measurement of radiance and degree of linear polarization (DOLP); and 2) aerosol extinction coefficient, scattering coefficient and depolarization ratio for fitting relevant lidar measurements. To guide the direction of search of solution through iterations, we run the forward model multiple times to calculate the Jacobian matrix that contains the derivatives of radiative and lidar quantities with respect to the aerosol and ocean surface properties. Moreover, a smoothness matrix is constructed to improve monotonic convergence of iteration toward the vicinity of globally optimized solution. A linear system is then established that includes the difference between model and observation, the Jacobian matrix, the smoothness matrix and the a priori. By solving the linear system, a stepwise increment of solution is determined to adjust the solution. In an iterative way, retrieval converges to a final solution.

## AEROSOL SCATTERING AND OCEAN SURFACE REFLECTION

Dependent on the information content offered by specific combinations of lidar and polarimeter, our algorithm is designed to be flexible in the aspect of allowing multiple aerosol species to be retrieved. In this section we will introduce the calculation of aerosol single scattering quantities as well as ocean surface reflection.

Adopting a log-normal volumetric size distribution $[dv_l^{(s)}(r)/d\ln r]$ for all ($N_s$) aerosol species, we have the subcolumn volumetric concentrations in $l$th layer,
<figure><img src="image_3.png" /><figcaption>FIGURE 1 | Algorithm flowcharts for retrieving aerosol and surface properties from a combined use of lidar and polarimeter.</figcaption></figure>

$$ \frac{dV_l(r)}{d \ln r} = \sum_{s=1}^{N_s} C_{v,l}^{(s)} \frac{dv_l^{(s)}(r)}{d \ln r}, \qquad (1) $$

where the size distribution of the sth species $dv_l^{(s)}(r)/d \ln r$ is parameterized by a geometric volume mean radius $r_m^{(s)}$ and a logarithmic standard deviation $\sigma_{m,l}^{(s)}$:

$$ \frac{dv_l^{(s)}(r)}{d \ln r} = \frac{1}{\sqrt{2\pi}\sigma_{m,l}^{(s)}} \exp\left[-\frac{(\ln r - \ln r_{m,l}^{(s)})^2}{2(\sigma_{m,l}^{(s)})^2}\right]. \qquad (2) $$

In Eq. 1, $C_{v,l}^{(s)}$ is the subcolumn volume concentration of the sth aerosol species in the lth layer and is related to its volume fraction ($f_{v,l}^{(s)}$) and the total column concentration ($C_{v,\text{tot}}^{(s)}$) by $C_{v,l}^{(s)} = C_{v,\text{tot}}^{(s)} f_{v,l}^{(s)}$.

If we discretize the specified aerosol size distribution in the lth layer into $N_{\text{bin}}$ size bins, the volumetric fraction of the ith bin is equal to $f_{i,l}^{(s)} d \ln r_i$ ($1 \le i \le N_{\text{bin}}$), where $f_{i,l}^{(s)}$ is calculated by

$$ f_{i,l}^{(s)} = \frac{dv_l^{(s)}(r_i)}{d \ln r_i}, \qquad (3) $$

where $r_i$ stands for the central radius of ith size bin and $dv_l^{(s)}(r_i)/d \ln r_i$ can be evaluated from Eq. 2.

## Aerosol Single Scattering Properties for Radiative Transfer Calculation

In this section, we are going to introduce the calculation of aerosol single scattering properties including AOD, SSA and phase matrix elements for radiative transfer modeling.

Accounting for the contribution by both spherical and spheroid-based nonspherical particles, the volumetric

extinction coefficient $K_{\text{ext}, a,l}^{(s)}$ and scattering coefficient $K_{\text{sca}, a,l}^{(s)}$ (in the unit of $m^{-1}$) are calculated by

$$ K_{\text{ext},a,l}^{(s)} = K_{\text{ext},a,l}^{(s),\text{sphere}} + K_{\text{ext},a,l}^{(s),\text{spheroid}} \qquad (4) $$

for sth aerosol species. In the above equation, the volumetric extinction component of spherical particles is evaluated by

$$ \begin{aligned} K_{\text{ext},a,l}^{(s),\text{sphere}} &= C_{n_0,a,l}^{(s),\text{sphere}} \int_{r_{\text{min}}}^{r_{\text{max}}} \frac{dn_l^{(s)}(r)}{d \ln r} c_{\text{ext},a,l}^{(s),\text{sphere}}(m_{a,l}^{(s)}, r) d \ln r \\ &= C_{v_0,a,l}^{(s),\text{sphere}} f_{v,l}^{(s),\text{sphere}} \int_{r_{\text{min}}}^{r_{\text{max}}} \frac{dv_l^{(s)}(r)}{d \ln r} \frac{1}{v(r)} c_{\text{ext},a,l}^{(s),\text{sphere}}(m_{a,l}^{(s)}, r) d \ln r, \end{aligned} \qquad (5) $$

where $C_{n_0,a,l}^{(s),\text{sphere}}$ is the aerosol number concentration ($m^{-3}$), $dn_l^{(s)}(r)/d \ln r$ is the number weighted aerosol size distribution in association with $dv_l^{(s)}(r)/d \ln r$, $r_{\text{min}}$ and $r_{\text{max}}$ are the lower and upper bounds of the size distribution respectively, $c_{\text{ext},a,l}^{(s),\text{sphere}}$ is extinction cross section as a function of particle refractive index $m_{a,l}^{(s)}$ and radius $r$, $f_{v,l}^{(s),\text{sphere}}$ is the volume fraction of spherical particles, $v(r)$ is the volume of a spherical particle of radius $r$, and $C_{v_0,a,l}^{(s)}$ is the aerosol volumetric concentration (in the unit of $m^3/m^3$). It is related to the subcolumn aerosol concentration ($C_{v,a,l}^{(s)}$ in the unit of $m^3/m^2$) in lth layer by

$$ C_{v,a,l}^{(s)} = C_{v_0,a,l}^{(s)} \Delta H_l = C_{v,\text{tot}}^{(s)} f_{v,l}^{(s)}, \qquad (6) $$

where

$$ f_{v,l}^{(s)} = \frac{C_{v,l}^{(s)} \Delta H_l}{C_{v,\text{tot}}^{(s)} \sum_{l=1}^{N_{\text{Layer}}} C_{v,l}^{(s)} \Delta H_l}. \qquad (7) $$
Calculation of volumetric extinction coefficient of spheroidal particles in Eq. 4 should account for both size and aspect ratio distributions of spheroids (Dubovik et al., 2011). Therefore,

$$ 
\begin{aligned}
K_{\text{ext},a,l}^{(s),\text{spheroid}} &= C_{v_0,a,l}^{(s),\text{spheroid}} \int_{r_{\min}}^{r_{\max}} \frac{dn_l^{(s)}(r)}{d \ln r} c_{\text{ext},a,l}^{(s),\text{spheroid}}(m_{a,l}^{(s)}, r) d \ln r \\
&= C_{v_0,a,l}^{(s)} (1 - f_{v,l}^{(s),\text{sphere}}) \int_{\eta_{\min}}^{\eta_{\max}} \frac{dn_l^{(s)}(\eta)}{d \ln \eta} \frac{dv_l^{(s)}(r)}{d \ln r} \frac{1}{v(r)} c_{\text{ext},a,l}^{(s),\text{spheroid}}(m_{a,l}^{(s)}, r, \eta) d \ln \eta d \ln r
\end{aligned}
\quad (8) $$

where r becomes the volume equivalent spherical particle radius, the aspect ratio $\eta < 1$ for oblate spheroids and $\eta > 1$ for prolate spheroids, and $\eta_{\min}$ and $\eta_{\max}$ are the lower and upper bounds of the aspect ratio distribution respectively. Note that we assume spheroids to have fixed aspect ratio distribution (Legrand et al., 2014) and $dn_l^{(s)}(\eta)/d \ln \eta$ is independent of aerosol size. Both $dn_l^{(s)}(\eta)/d \ln \eta$ and $dv_l^{(s)}(r)/d \ln r$ are normalized to unity in the presence of non-zero lower bound or non-infinity upper bound in their integral over r and $\eta$, respectively.

For numerical calculation, the aerosol size distribution is descritized into $N_{\text{bin}}$ rectangular bins so that for spherical particles

$$ K_{\text{ext},a,l}^{(s),\text{sphere}} = C_{v_0,a,l}^{(s)} f_{v,l}^{(s),\text{sphere}} \sum_{i=1}^{N_{\text{bin}}} \frac{dv_l^{(s)}(r_i)}{d \ln r} K_{\text{ext},a,l,i}^{(s),\text{sphere}}(m_{a,l}^{(s)}, r_i), \quad (9) $$

where the component of sth aerosol species in its ith size bin is

$$ K_{\text{ext},a,l,i}^{(s),\text{sphere}}(m_{a,l}^{(s)}, r_i) = \int_{\ln r_i - \Delta \ln r_i / 2}^{\ln r_i + \Delta \ln r_i / 2} \frac{c_{\text{ext},a,l}^{(s),\text{sphere}}(m_{a,l}^{(s)}, r)}{v(r)} d\ln r. \quad (10) $$

For spheroids, we further need to break the integral of volumetric extinction into $N_{\text{AR}}$ bins of aspect ratio so that:

$$ K_{\text{ext},a,l}^{(s),\text{spheroid}} = C_{v_0,a,l}^{(s)} (1 - f_{v,l}^{(s),\text{sphere}}) \sum_{i=1}^{N_{\text{bin}}} K_{\text{ext},a,l,i}^{(s),\text{spheroid}} (m_{a,l}^{(s)}, r_i) \frac{dv_l^{(s)}(r_i)}{d \ln r}, \quad (11) $$

where

$$ K_{\text{ext},a,l,i}^{(s),\text{spheroid}}(m_{a,l}^{(s)}, r_i) = \sum_{q=1}^{N_{\text{AR}}} \frac{dn_l^{(s)}(\eta_q)}{d \ln \eta} \int_{\ln \eta_q - \Delta \ln \eta_q / 2}^{\ln \eta_q + \Delta \ln \eta_q / 2} \int_{\ln r_i - \Delta \ln r_i / 2}^{\ln r_i + \Delta \ln r_i / 2} \frac{c_{\text{ext},a,l}^{(s),\text{spheroid}}(m_{a,l}^{(s)}, r, \eta)}{v(r)} d\ln r d\ln \eta. \quad (12) $$

For simplicity, the layer resolved volumetric extinction coefficient can be rewritten in the following form,

$$ K_{\text{ext},a,l}^{(s)} = \sum_{i=1}^{N_{\text{bin}}} f_{i,l}^{(s)} \left[ f_{v,l}^{(s),\text{sphere}} K_{\text{ext},a,l,i}^{(s),\text{sphere}} + (1 - f_{v,l}^{(s),\text{sphere}}) K_{\text{ext},a,l,i}^{(s),\text{spheroid}} \right], \quad (13) $$

where $K_{\text{ext},a,l,i}^{(s),\text{sphere}}$ and $K_{\text{ext},a,l,i}^{(s),\text{spheroid}}$ are given in Eqs 10 and 12, respectively.

Replacing the subscript “ext” by “sca” in above derivations (Eqs. 4–13) gives relevant equations for the volumetric scattering coefficient:

$$ K_{\text{sca},a,l}^{(s)} = \sum_{i=1}^{N_{\text{bin}}} f_{i,l}^{(s)} \left[ f_{v,l}^{(s),\text{sphere}} K_{\text{sca},a,l,i}^{(s),\text{sphere}} + (1 - f_{v,l}^{(s),\text{sphere}}) K_{\text{sca},a,l,i}^{(s),\text{spheroid}} \right]. \quad (14) $$

In a similar way, the phase matrix element $P_{a,j,k,l}$ of sth aerosol species in the lth layer is calculated from

$$ P_{a,j,k,l}^{(s)} = \frac{\sum_{i=1}^{N_{\text{bin}}} f_{i,l}^{(s)} \left[ f_{v,l}^{(s),\text{sphere}} K_{a,\text{sca},l,i}^{(s),\text{sphere}} P_{a,j,k,l,i}^{(s),\text{sphere}} + (1-f_{v,l}^{(s),\text{sphere}}) K_{a,\text{sca},l,i}^{(s),\text{spheroid}} P_{a,j,k,l,i}^{(s),\text{spheroid}} \right]}{\sum_{i=1}^{N_{\text{bin}}} f_{i,l}^{(s)} \left[ f_{v,l}^{(s),\text{sphere}} K_{a,\text{sca},l,i}^{(s),\text{sphere}} + (1-f_{v,l}^{(s),\text{sphere}}) K_{a,\text{sca},l,i}^{(s),\text{spheroid}} \right]}, \quad (15) $$

In the numerator of the above equation, we have

$$ K_{\text{sca},a,l,i}^{(s),\text{sphere}} P_{a,j,k,l,i}^{(s),\text{sphere}} = \int_{\ln r_i - \Delta \ln r_i / 2}^{\ln r_i + \Delta \ln r_i / 2} \frac{c_{\text{sca},a,l}^{(s),\text{sphere}}(m_{a,l}^{(s)}, r)}{v(r)} P_{a,j,k,l,i}^{(s),\text{sphere}}(m_{a,l}^{(s)}, r) d\ln r, \quad (16) $$

$$ K_{\text{sca},a,l,i}^{(s),\text{spheroid}} P_{a,j,k,l,i}^{(s),\text{spheroid}} = \sum_{q=1}^{N_{\text{AR}}} \frac{dn_l^{(s)}(\eta_q)}{d \ln \eta} \int_{\ln \eta_q - \Delta \ln \eta_q / 2}^{\ln \eta_q + \Delta \ln \eta_q / 2} \int_{\ln r_i - \Delta \ln r_i / 2}^{\ln r_i + \Delta \ln r_i / 2} \frac{c_{\text{sca},a,l}^{(s),\text{spheroid}}(m_{a,l}^{(s)}, r, \eta)}{v(r)} P_{a,j,k,l,q}^{(s),\text{spheroid}}(m_{a,l}^{(s)}, r, \eta) d\ln r d\ln \eta. \quad (17) $$

Summing over the contributions from all aerosol species gives the total volumetric extinction and scattering coefficients in lth layer,

$$ K_{\text{ext}[\text{sca}],a,l} = \sum_{s=1}^{N_s} K_{\text{ext}[\text{sca}],a,l}^{(s)}, \quad (18) $$

where $K_{\text{ext},a,l}^{(s)}$ and $K_{\text{sca},a,l}^{(s)}$ are given by Eqs 13 and 14 respectively.

The speciated aerosol SSA in lth layer is calculated as the ratio of speciated scattering and extinction coefficients by,

$$ \omega_{0,a,l}^{(s)} = \frac{K_{\text{sca},a,l}^{(s)}}{K_{\text{ext},a,l}^{(s)}}, \quad (19) $$

Contributed by all aerosol species, the overall SSA contributed by all aerosol species in lth layer is

$$ \omega_{0,a,l} = \frac{\sum_{s=1}^{N_s} f_{v,l}^{(s)} K_{\text{sca},a,l}^{(s)}}{\sum_{s=1}^{N_s} f_{v,l}^{(s)} K_{\text{ext},a,l}^{(s)}}, \quad (20) $$

where $f_{v,l}^{(s)}$ stands for the volume fraction of sth aerosol species in lth layer, as given by Eq. 7.

The overall aerosol phase matrix elements are evaluated by,

$$ P_{a,j,k,l} = \frac{\sum_{s=1}^{N_s} f_{v,l}^{(s)} K_{\text{sca},a,l}^{(s)} P_{a,j,k,l}^{(s)}}{\sum_{s=1}^{N_s} f_{v,l}^{(s)} K_{\text{sca},a,l}^{(s)}}, \quad (21) $$
where $P_{a,j,k,l}^{(s)}$ is given by Eq. 15. The layer aerosol optical depth (AOD, $\Delta\tau_{a,l}$) is contributed by all aerosol species. Therefore, we have

$$ \Delta\tau_{a,l} = \sum_{s=1}^{N_s} \Delta\tau_{a,l}^{(s)}. \qquad (22) $$

In the above equation, the AOD for sth aerosol species is

$$ \Delta\tau_{a,l}^{(s)} = C_{v, a, l}^{(s)} K_{\text{ext}, a, l}^{(s)}, \qquad (23) $$

where $C_{v,a,l}^{(s)}$ is the subcolumn volumetric concentration of sth aerosol species, evaluated by Eq. 6, and $K_{\text{ext},a,l}^{(s)}$ is column effective extinction coefficient (in the unit of m⁻¹) of sth aerosol species in the lth layer, evaluated by Eq. 13.

## Aerosol Single Scattering Properties for Lidar Equation Use

In this section, we are going to introduce the calculation of aerosol single scattering properties including extinction, backscatter and depolarization for use in the lidar equation.

The total extinction coefficient for lth layer ($\alpha_{a,l}$, in the unit of m⁻¹) is calculated by summing over the contribution from all aerosol species, namely

$$ \alpha_{a,l} = \sum_{s=1}^{N_s} \alpha_{a,l}^{(s)}, \qquad (24) $$

where the specified aerosol extinction $\alpha_{a,l}^{(s)}$ in a layer of geometric thickness $\Delta H_l$ is related to layer resolved aerosol optical depth by

$$ \alpha_{a,l}^{(s)} = \frac{\Delta\tau_{a,l}^{(s)}}{\Delta H_l} = \frac{C_{v,a,l}^{(s)} K_{\text{ext},a,l}^{(s)}}{\Delta H_l}, \qquad (25) $$

where $C_{v,a,l}^{(s)}$ and $K_{\text{ext},a,l}^{(s)}$ are given by Eqs 6 and 7 respectively.

Calculation of the total backscatter $\beta_{a,l}$ (in the unit of m⁻¹Sr⁻¹) also needs to account for the contributions of all aerosol species as well, namely

$$ \beta_{a,l} = \sum_{s=1}^{N_s} \beta_{a,l}^{(s)}, \qquad (26) $$

where the specified aerosol backscatter is evaluated by

$$ \beta_{a,l}^{(s)} = \beta_{\text{sca},a,l}^{(s)} \frac{P_{a,1,1,l}^{(s)} (180^{\circ})}{4\pi}, \qquad (27) $$

where $P_{a,1,1,l}^{(s)}(180^{\circ})$ is given by Eq. 15 and the aerosol scattering coefficient ($\beta_{\text{sca},a,l}^{(s)}$, in the unit of m⁻¹) in lth layer is evaluated by,

$$ \beta_{\text{sca},a,l}^{(s)} = \frac{C_{v,a,l}^{(s)} K_{\text{sca},a,l}^{(s)}}{\Delta H_l}, \qquad (28) $$

where $K_{\text{sca},a,l}^{(s)}$ is given by Eq. 16.

Assuming random orientation of nonspherical particles, the aerosol linear depolarization ratio of lth layer is calculated as

$$ \delta_{a,l}(180^{\circ}) = \frac{1 - P_{a,2,2,l}(180^{\circ})/P_{a,1,1,l}(180^{\circ})}{1 + P_{a,2,2,l}(180^{\circ})/P_{a,1,1,l}(180^{\circ})}, \qquad (29) $$

where the layer resolved phase matrix elements $P_{a,i,j,l}$ are evaluated via Eq. 21.

## Development of a Lookup Table for Aerosol Single Scattering Properties

Calculation of aerosol single scattering properties (AOD, SSA, p) for RT modeling and ($\alpha_a$, $\beta_a$, $\delta_a$) for fitting lidar measurements needs to use the aerosol refractive index as the input. We adopt the following two parameter constrained model for refractive index as a function of wavelength (Dubovik and King, 2000),

$$ \begin{aligned} m_{r, l}(\lambda) &= m_{r0, l}(\lambda/\lambda_0)^{-\kappa_{r0, l}} \\ m_{i, l}(\lambda) &= m_{i0, l}(\lambda/\lambda_0)^{-\kappa_{i0, l}} \end{aligned} \qquad (30) $$

where $m_{r0}$ and $m_{i0}$ are the real and imaginary parts of the refractive index at the reference wavelength $\lambda_0$, and $\kappa_{r0}$ and $\kappa_{i0}$ are used to characterize their spectral dependence.

During the optimization process, the spectrally dependent refractive index and size distribution of aerosols and the spherical particle fraction of aerosols are updated dynamically. To avoid inefficient on-the-fly light scattering computations of aerosol extinction [scattering] coefficients $K_{\text{ext}[\text{sca}], a,l}^{\text{sphere}[\text{spheroid}]}$, and aerosol phase matrix $P_{a,j,k,i}^{\text{sphere}[\text{spheroid}]}$ are precalculated for a set of size bins and refractive indices and saved in a lookup table. The table is then linearly interpolated to determine the scattering properties associated with dynamically updated refractive index during the retrieval. The bin resolved properties are weighted by size distribution parameters to determine $K_{\text{ext}[\text{sca}], a,i}^{\text{sphere}[\text{spheroid}]}$, $P_{a,j,k,i}^{\text{sphere}[\text{spheroid}]}$ and then ($\alpha_a$, $\beta_a$, $\delta_a$) for the whole size distribution (see Aerosol Single Scattering Properties for Radiative Transfer Calculation and Aerosol Single Scattering Properties for Lidar Equation Use).

In the generation of the lookup table, the real part of refractive index is uniformly distributed in linear space from 1.3 to 1.7 with 22 grids. The imaginary part of refractive index is uniformly distributed in logarithmic space from 1 × 10⁻⁸ to 3 × 10⁻⁴ with five grids and from 5 × 10⁻⁴ to 0.5 with 15 grids. The lower and upper bounds for aerosol size are set as $r_{\text{min}} = 5$ nm and $r_{\text{max}} = 50$ µm, respectively. A total of 50 size bins are equally spaced in logarithmic space to calculate aerosol light scattering properties. Mie theory is used to calculate spherical particle scattering (van de Hulst, 1981). For spheroids, the lower and upper bounds for aerosol axis ratio to be 0.33 and 3.0 respectively and a total of 25 bins are set for precalculating spheroidal scattering properties by use of T-matrix method in combination with geometric optics (Dubovik et al., 2006).

## Surface Reflection

The use of bio-optical model to simulate ocean water-leaving radiance effectively reduces the parameter space and provides extra constraints in the retrieval of ocean bulk properties. It has been more and more adopted over the last decade and justified by many research groups (see for example [Chowdhary et al., 2012; Xu et al., 2016; Gao et al., 2018; [van de Hulst, 1981]).
Stamnes et al., 2018; among others)). Our ocean surface reflection model incorporates a polarized specular reflection term, a Lambertian term for depolarizing ocean foam reflection and an empirical Lambertian correction term “Δa” to account for the errors of the single-parameter based bio-optical model (i.e., departures from the predetermined functional relationships to Chl-a). The overall bidirectional ocean surface reflection matrix R<sub>surf</sub> is described by (Xu et al., 2016),

$$ \pi \mathbf{R}_{\text{surf},\lambda} = f_{\text{Foam}} a_{\text{Foam},\lambda} \mathbf{D}_0 + (1 - f_{\text{Foam}}) \mathbf{R}_{\text{W},\lambda} + (1 - f_{\text{Foam}}) \mathbf{R}_{\text{WL},\lambda}^{\text{Bio}} \\ + (1 - f_{\text{Foam}}) \Delta a_{\text{WL},\lambda} \mathbf{D}_0, \qquad (31) $$

where D<sub>0</sub> is a zero matrix except D<sub>0,11</sub> = 1; a<sub>Foam</sub> is foam albedo; f<sub>Foam</sub> is foam coverage fraction related to wind speed V<sub>wind</sub> by f<sub>Foam</sub> = 2.95 × 10<sup>-6</sup>×V<sub>wind</sub><sup>3.52</sup> (Koepke, 1984); R<sub>W,λ</sub> is the polarized specular reflection from water surface (Cox and Munk, 1954a; Cox and Munk, 1954b), and R<sub>WL,λ</sub><sup>Bio</sup> is the reflection matrix of the ocean-interface system with a correction of Raman scattering (Lee et al., 2013; Xu et al., 2016), and “Δa<sub>WL,λ</sub>” is the Lambertian correction term. Note that R<sub>WL,λ</sub><sup>Bio</sup> is a physically based term depending on the Chl-a concentration [Chl_a]. The last two terms of Eq. 31 constitute our water-leaving radiance model. With and without assuming Δa<sub>WL</sub> to be 0 the simplified and the empirically adjusted bio-optical models are formulated, respectively. By setting the Sun at zenith and viewing angle to be nadir, we can evaluate the normalized water-leaving radiance from R<sub>WL</sub><sup>Bio</sup> and Δa<sub>WL</sub>,

$$ \begin{aligned} \text{nLw}_{\lambda} = & \frac{F_{0,\lambda}}{\pi} \left(\frac{d_0}{d}\right)^2 [R_{\text{WL},\lambda,11}^{\text{Bio}} (\theta_v = 0^\circ; \theta_0 \\ & = 0^\circ; [\text{Chl_a}]) + \Delta a_{\text{WL},\lambda}], \end{aligned} \qquad (32) $$

where d<sub>0</sub> is the Earth-Sun distance at which F<sub>0</sub> is reported, and d is the Earth-Sun distance at the time of measurement. Note that nLw, R<sub>W</sub>, R<sub>surf</sub>, R<sub>WL</sub><sup>Bio</sup>, a<sub>Foam</sub>, Δa<sub>WL</sub>, and F<sub>0</sub> in Eqs 31 and 32 are all spectrally dependent.

# FORWARD MODEL FOR LIGHT PROPAGATION IN A COUPLED ATMOSPHERE-OCEAN SYSTEM

## Polarized Radiative Transfer

As illustrated in Figure 2, a coupled atmosphere-surface model is established for 1D RT modeling, which consists of an aerosol/air-molecule atmospheric layer and ocean water at the bottom of the atmosphere. The whole atmosphere is subdivided into N<sub>Layer</sub> layers, each bounded by the altitudes h<sub>l</sub> and h<sub>l+1</sub> (h<sub>l</sub> < h<sub>l+1</sub>) and N<sub>Layer</sub> depends on the vertical resolution of lidar measurement or the desired resolution for retrieval products.

As the input to radiative transfer modeling, vertically resolved AOD, SSA, and phase matrices of aerosols in all

<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | The coupled atmosphere-ocean system for RT modeling. The Sun illuminates the top of atmosphere with zenith angle θ<sub>0</sub> and azimuthal angle φ<sub>0</sub> (φ<sub>0</sub> is assumed to be zero here). The sensor views the atmosphere at viewing angle θ<sub>v</sub> and azimuthal angle φ<sub>v</sub>. The Markov chain model is used for computing polarized RT in the atmosphere. The ocean is assumed optically homogeneous so that the doubling method can be applied to achieve modeling efficiency. Coupling of local radiative fields from ocean and atmospheric layer is completed by using an “adding” strategy.</figcaption></figure>

layers are calculated from specified aerosol refractive index, size distributions and spherical particle fraction. They are mixed with molecule scattering to get bulk scattering properties in all layers. Consequently, the optical thickness (Δτ<sub>l</sub>), SSA (ω<sub>0,l</sub>) and phase matrix (P<sub>l</sub>) of lth layer contributed by both aerosols and Rayleigh-scattering molecules are evaluated by

$$ \Delta\tau_l = \Delta\tau_{a,l} + \Delta\tau_{R,l}, \qquad (33) $$

$$ \omega_{0,l} = \frac{\Delta\tau_{R,l} + \omega_{0,a,l}\Delta\tau_{a,l}}{\Delta\tau_{R,l} + \Delta\tau_{a,l}}, \qquad (34) $$

$$ \mathbf{P}_{l}(\Theta) = \frac{\Delta\tau_{R,l}\mathbf{P}_{R}(\Theta) + \omega_{0,a,l}\Delta\tau_{a,l}\mathbf{P}_{a,l}(\Theta)}{\Delta\tau_{R,l} + \omega_{0,a,l}\Delta\tau_{a,l}}, \qquad (35) $$

where as a function of scattering angle Θ, P<sub>R</sub> and P<sub>a,l</sub> are the Rayleigh and aerosol phase matrices, respectively, the aerosol SSA of lth layer (ω<sub>0,a,l</sub>) is given by Eq. 20 and the AOD of lth layer (Δτ<sub>a,l</sub>) is given by Eq. 22.

Knowing the layer optical depth, SSA and phase matrices for all atmospheric layers, an efficient calculation of reflection and transmission matrices of the vertically inhomogeneous atmosphere is implemented by use of the Markov chain method (Esposito and House, 1978; Esposito 1979; Xu et al., 2010; Xu et al., 2011; Xu et al., 2012). As an approximation, the ocean water constituents are assumed to be homogeneously mixed. Then the doubling method (see for example [Stokes, 1862; van de Hulst, 1963; Hansen, 1971; de Haan et al., 1987;
Evans and Stephens, 1991; among others]) is applied for a fast calculation of the reflection matrices of the ocean bulk. Finally, the interaction between atmosphere and ocean through multiple reflection and transmission are coupled via an “adding” approach (Xu et al., 2016). The output of atmosphere and ocean reflection is the reflection matrix for the coupled atmosphere-ocean system, from which the radiance and polarimetric components are derived to compare to observations. The incorporation of Markov chain model and doubling methods into a hybrid RT modeling scheme capitalizes on the strength of these two methods in modeling RT in vertically inhomogeneous and homogeneous media, respectively. Moreover, it provides an efficient Jacobian calculation and an optimization-based retrieval by saving the local RT fields (reflection and transmission matrices) for each layer from forward modeling, which are then reused for subsequent Jacobian evaluations as long as the aerosol or ocean parameters remain unperturbed. In other words, only the RT in the layer with perturbed parameters requires recalculation.

## Lidar Equation

The propagation of a light pulse from a single-channel lidar in the atmosphere is governed by lidar equation (e.g., (Winker et al., 1996; Vaughan, 2004; Weitkamp, 2005; Burton et al., 2018)):

$$ \frac{P_p(z)}{g_p} = \frac{1}{z^2} [\beta_m(z) + \beta_a(z)] \exp\left\{-2 \int_0^z [\alpha_m(z') + \alpha_a(z')] dz'\right\}, \quad (36) $$

where $P_p(z)$ is the calibrated signal at the detector as a function of the range (distance, z) between the lidar and the target atmospheric volume, $g_p$ is a gain factor which takes into account detector and receiver optical efficiencies, $\beta_a(z)$ and $\beta_m(z)$ are backscatter coefficients of aerosols and air molecules at range z, and $\alpha_a$ and $\alpha_m$ are the aerosol and molecular extinction coefficients respectively. The exponential term on the right-hand side of the above equation represents two-way attenuation of the lidar signal between the detector and the atmospheric volume.

The multiplication of the exponential term on the right hand side of Eq. 36 and the term “[$\beta_a(z) + \beta_m(z)$]” gives attenuated backscatter. The molecular scattering and extinction coefficients $\beta_m(r)$ and $\alpha_m(r)$ can be accurately determined from measured atmospheric density or from climatology. What remain to be resolved are the aerosol extinction and backscatter coefficients. However, with these two unknowns coupled in the attenuated backscatter lidar system via Eq. 36, they cannot be separated from one measured variable on the left-hand-side of the equation. By adopting the HSRL measurement technique, however, the received backscatter is split into two channels, one with a very narrow bandwidth optical absorption filter in the molecular channel to suppress the aerosol backscatter, whereas another is a combined channel that detects the intensity contributed by both aerosol and molecules. In this case, the lidar equation for the combined signal remains the same as the above equation while for the molecular channel with narrow bandwidth optical absorption filter the lidar equation becomes [e.g., Burton et al., 2015],

$$ \frac{P_p(z)}{g_p} = \frac{1}{z^2} \beta_m(z) \exp\left\{-2 \int_0^z [\alpha_m(z') + \alpha_a(z')] dz'\right\}. \quad (37) $$

High-spectral-resolution lidars such as NASA Langley's HSRL-1 and HSRL-2 send a linearly polarized beam into an atmosphere. The returned lidar signals are measured with linear polarization analyzers oriented both parallel (||) and perpendicular (⊥) to the emitted beam. By taking the ratio of light intensities of these two components, the depolarization ratio of the returned transmitted signal is derived as:

$$ \delta = \frac{I_{\perp}}{I_{\parallel}}. \quad (38) $$

Due to the anisotropy of the air molecules, the depolarization ratio of the pure molecular atmosphere is nonzero but can be pre-calculated from climatology. Moreover, the backscatter of a linearly polarized laser beam from spherical particles is totally linearly polarized so that ($\delta = 0$). However, the presence of nonspherical aerosols causes the measured depolarization ratio to deviate from zero which provides extra constraints on retrieving non-sphericity of aerosols.

The lidar equation Eq. 37 is single scattering based. It can break into a discretized form, namely multiplication of local backscatter at range z with the sum of the contributions of extinction in all layers above z. The calculation of aerosol component of layer effective extinction ($\alpha_a$), backscatter ($\beta_a$) and depolarization ($\delta_a$) have been specified in Aerosol Single Scattering Properties for Lidar Equation Use. In addition, to combine lidar observations with polarimetric measurements: one can structure the retrieval to fit either direct lidar measurements such as total extinction, backscatter and depolarization, or its aerosol products such as particulate extinction, backscatter and depolarization ratio. The latter is adopted in the current retrieval demonstration in Retrieval Implementation and Validation.

## COUPLED RETRIEVAL OF AEROSOLS AND OCEAN PROPERTIES

### State Vector

Our algorithm is developed to be flexible and allow for multiple ($N_s$) aerosol species to be retrieved from a combined use of lidar and polarimetric measurements depending on the information content of 1) the specific configuration of lidar measurement capabilities in terms of providing the whole set or a subset of extinction, backscatter/ attenuated backscatter, and depolarization ratio measurements, their vertical resolutions and number of wavelengths; and 2) the number of radiometric and polarimetric view angles and spectral channels provided by a polarimeter. The state vector $\mathbf{x}$ that includes retrieval parameters for $N_s$ aerosol species can be expressed as,

$$ \mathbf{x} = [\mathbf{x}^{(1)}; \mathbf{x}^{(2)}; \dots; \mathbf{x}^{(N_s)}; \mathbf{x}^{(\text{ocean})}]_{N \times 1}, \quad (39) $$
<table><caption>**TABLE 1** | Parameters in coupled aerosol and surface property retrieval, their initial estimates, and order of difference and Lagrange multipliers for imposing smoothness constraints.</caption><thead><tr><th>Aerosol parameters (layer resolved)</th><th></th><th>Range</th><th>Initial guess</th><th>Order of finite<br>difference for vertical/<br>spectral<br>smoothness<br>constraints (m<sub>h</sub>)</th><th>Lagrange regularization<br>factor<br>(γ<sub>h</sub>)</th></tr></thead><tbody><tr><td rowspan="4">Real part of refractive index - parameter<br>1 (m<sub>r0</sub>)</td><td>Species 1</td><td>[1.33,1.60]</td><td>1.45</td><td>m<sub>h</sub> = 1</td><td>γ<sub>h</sub> = 10</td></tr><tr><td>Species 2</td><td>[1.33,1.60]</td><td>1.35</td><td>m<sub>h</sub> = 1</td><td>γ<sub>h</sub> = 10</td></tr><tr><td>Species 1</td><td>[5 × 10<sup>-4</sup>, 5 × 10<sup>-1</sup>]</td><td>1 × 10<sup>-1</sup></td><td>m<sub>h</sub> = 1</td><td>γ<sub>h</sub> = 5</td></tr><tr><td>Species 2</td><td>[5 × 10<sup>-8</sup>, 5 × 10<sup>-2</sup>]</td><td>1 × 10<sup>-5</sup></td><td>m<sub>h</sub> = 1</td><td>γ<sub>h</sub> = 5</td></tr><tr><td rowspan="4">Real part of refractive index - parameter<br>2 (κ<sub>r0</sub>)</td><td>Species 1</td><td>[1 × 10<sup>-5</sup>, 1 × 10<sup>3</sup>]</td><td>1 × 10<sup>-2</sup></td><td>m<sub>h</sub> = 1</td><td>γ<sub>h</sub> = 10</td></tr><tr><td>Species 2</td><td>[1 × 10<sup>-5</sup>, 1 × 10<sup>3</sup>]</td><td>1 × 10<sup>-2</sup></td><td>m<sub>h</sub> = 1</td><td>γ<sub>h</sub> = 10</td></tr><tr><td>Species 1</td><td>[1 × 10<sup>-5</sup>, 1 × 10<sup>3</sup>]</td><td>5 × 10<sup>-1</sup></td><td>m<sub>h</sub> = 1</td><td>γ<sub>h</sub> = 5</td></tr><tr><td>Species 2</td><td>[1 × 10<sup>-5</sup>, 1 × 10<sup>3</sup>]</td><td>5 × 10<sup>-1</sup></td><td>m<sub>h</sub> = 1</td><td>γ<sub>h</sub> = 5</td></tr><tr><td rowspan="2">Imaginary part of refractive index -<br>parameter 2 (κ<sub>i0</sub>)</td><td>Species 1</td><td>[1 × 10<sup>-2</sup>, 5 × 10<sup>-1</sup>]</td><td>1 × 10<sup>-1</sup></td><td>m<sub>h</sub> = 2</td><td>γ<sub>h</sub> = 5 × 10<sup>-1</sup></td></tr><tr><td>Species 2</td><td>[1 × 10<sup>-1</sup>, 5 × 10<sup>1</sup>]</td><td>1 × 10<sup>-1</sup></td><td>m<sub>h</sub> = 2</td><td>γ<sub>h</sub> = 5 × 10<sup>-1</sup></td></tr><tr><td rowspan="2">Imaginary part of refractive index -<br>parameter 2 (κ<sub>i0</sub>)</td><td>Species 1</td><td>[0.9,1.1]</td><td>UHSAS in-situ<br>measurement<sup>a,b</sup></td><td>m<sub>h</sub> = 2</td><td>γ<sub>h</sub> = 5 × 10<sup>-1</sup></td></tr><tr><td>Species 2</td><td>[0.9,1.1]xσ<sub>APS</sub><sup>c,b</sup></td><td>APS in-situ measurement<sup>c,b</sup></td><td>m<sub>h</sub> = 2</td><td>γ<sub>h</sub> = 5 × 10<sup>-1</sup></td></tr><tr><td rowspan="2">Volumetric median aerosol radius (r<sub>m</sub>)</td><td>Species 1</td><td>[1 × 10<sup>-5</sup>, 1 × 10<sup>3</sup>]</td><td>Associated with AOD = 0.005/<br>N<sub>Layer</sub></td><td>m<sub>h</sub> = 3</td><td>γ<sub>h</sub> = 5 × 10<sup>-3</sup></td></tr><tr><td>Species 2</td><td>[1 × 10<sup>-5</sup>, 1 × 10<sup>3</sup>]</td><td>Associated with AOD = 0.005/<br>N<sub>Layer</sub></td><td>m<sub>h</sub> = 3</td><td>γ<sub>h</sub> = 5 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">Standard deviation of size distribution (σ<sub>m</sub>)</td><td>Species 1</td><td>[0.9,1.1]</td><td>UHSAS in-situ<br>measurement<sup>a,b</sup></td><td>m<sub>h</sub> = 2</td><td>γ<sub>h</sub> = 5 × 10<sup>-1</sup></td></tr><tr><td>Species 2</td><td>[0.9,1.1]xσ<sub>APS</sub><sup>c,b</sup></td><td>APS in-situ measurement<sup>c,b</sup></td><td>m<sub>h</sub> = 2</td><td>γ<sub>h</sub> = 5 × 10<sup>-1</sup></td></tr><tr><td rowspan="2">Aerosol subcolumn volume concentration<br>(C<sub>v</sub>, µm<sup>3</sup>/µm<sup>2</sup>)</td><td>Species 1</td><td>[1 × 10<sup>-8</sup>, 5]</td><td>Associated with AOD = 0.005/<br>N<sub>Layer</sub></td><td>m<sub>h</sub> = 3</td><td>γ<sub>h</sub> = 5 × 10<sup>-3</sup></td></tr><tr><td>Species 2</td><td>[1 × 10<sup>-8</sup>, 5]</td><td>Associated with AOD = 0.005/<br>N<sub>Layer</sub></td><td>m<sub>h</sub> = 3</td><td>γ<sub>h</sub> = 5 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">Spherical particle fraction (f<sub>v</sub><sup>sphere</sup>)</td><td>Species 1</td><td>[1 × 10<sup>-3</sup>, 1 × 10<sup>1</sup>]</td><td>9 × 10<sup>-1</sup></td><td>m<sub>h</sub> = 2</td><td>γ<sub>h</sub> = 1 × 10<sup>-1</sup></td></tr><tr><td>Species 2</td><td>[1 × 10<sup>-3</sup>, 1 × 10<sup>1</sup>]</td><td>9 × 10<sup>-1</sup></td><td>m<sub>h</sub> = 2</td><td>γ<sub>h</sub> = 1 × 10<sup>-1</sup></td></tr><tr><td>Surface parameters (ocean)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Adjustment term (Δa<sub>WL</sub>(λ), mw/cm<sup>2</sup>-sr-µm)</td><td></td><td>[1 × 10<sup>-5</sup>, 5 × 10<sup>1</sup>]</td><td>1 × 10<sup>-5</sup></td><td>—</td><td>—</td></tr><tr><td>Chlorophyll a concentration ([Chl_a], mg/m<sup>3</sup>)</td><td></td><td>[1 × 10<sup>-2</sup>, 6 × 10<sup>1</sup>]</td><td>3 × 10<sup>-1</sup></td><td>—</td><td>—</td></tr><tr><td>Surface wind speed (V<sub>wind</sub>, m/sec)</td><td></td><td>[1 × 10<sup>1</sup>, 3 × 10<sup>1</sup>]</td><td>5</td><td>—</td><td>—</td></tr></tbody></table>

where N is the total number of aerosol and ocean properties to be retrieved, the subvector x<sup>(s)</sup> contains aerosol properties for sth aerosol species. For a total of N<sub>aer</sub> = 8 aerosol parameters per species, we have,

$$ \mathbf{x}^{(s)} = [C_{v}^{(s)}; \mathbf{m}_{r_0}^{(s)}; \mathbf{k}_{r_0}^{(s)}; \mathbf{m}_{i_0}^{(s)}; \mathbf{k}_{i_0}^{(s)}; \mathbf{r}_{m}^{(s)}; \mathbf{\sigma}_{m}^{(s)}; f_{v}^{(s),sphere}]_{N_{\text{aer}} \times 1} \quad (40) $$

where for sth aerosol species C<sub>v</sub> contains layer specific column aerosol concentration (m<sup>3</sup>/m<sup>2</sup>), (m<sub>r0</sub>, κ<sub>r0</sub>) and (m<sub>i0</sub>, κ<sub>i0</sub>) are layer resolved model parameters that constrain spectrally dependent real and imaginary parts of aerosol refractive index, respectively (via Eq. 30); r<sub>m</sub> and σ<sub>m</sub> contain layer resolved geometric volume median radius and standard deviation of volumetric log-normal distribution; f<sub>v</sub><sup>sphere</sup> contain the layer resolved volume fraction of spherical aerosols. As the last component of the state vector in Eq. 39, x<sup>(ocean)</sup> contains parameters that contribute to ocean surface reflection. These parameters include a simplified bio-optical model parameter chlorophyll-a concentration, wind speed and the empirical terms (Lambertian) that mitigate the bio-optical modeling errors by adjusting the spectral water-leaving radiance (see Development of a

Lookup Table for Aerosol Single Scattering Properties). Retrieval starts with an initial estimate of these parameters, as specified in Table 1.

## Constraints

Imposition of various types of constraints is important for stabilizing the retrievals and improving accuracy. The combined lidar-polarimeter inversion approach described here uses three types of constraints: 1) direct observation signals provided by lidar and polarimetric remote sensing measurements (see formulation in *Observational Constraints* (i = 1)); 2) a priori constraint for the retrieval parameters (see formulation in *A Priori Constraints* (i = 2)); and 3) smoothness constraint over the variations of certain aerosol properties along altitude (see formulation in *Smoothness Constraints* (i = 3)). Incorporation of all these three types of constraints into our combined inversion approach is described in *Construction of Overall Equation System*.

Imposition of the above three types of constraints (M = 3) can be formulated in the framework of statistical inversion of multi-source constraints. Our

<footnote><sup>a</sup>UHSAS, Ultra-High Sensitivity Aerosol spectrometer; σ<sub>UHSAS</sub> is derived from fitting in-situ size measurement by UHSAS using the log-normal size distribution (Eq. 103).</footnote>

<footnote><sup>c</sup>APS, Aerodynamic Particle Sizer (APS-TSI3320); σ<sub>APS</sub> is derived from fitting in-situ size measurement by APS using the log-normal size distribution (Eq. 103).</footnote>

<footnote><sup>b</sup>σ<sub>UHSSA</sub> and σ<sub>APS</sub> are adopted when measured lidar extinction is larger than its uncertainty; otherwise the initial guess of σ<sub>m</sub> is set to 0.2 and its range is relaxed to [0.05, 0.60].</footnote>
objective is to solve the following system of equations (Dubovik, 2004):

$$ \mathbf{f}_i^* = \mathbf{f}_i(\mathbf{x}) + \Delta\mathbf{f}_i^*, \quad 1 \le i \le M \qquad (41) $$

where $\mathbf{f}_i^*$ denotes the ith type of constraint, $\Delta\mathbf{f}_i^*$ is the error with this type of constraint, and $\mathbf{x} = \mathbf{x}_{\text{state}}$ as defined in Eq. 39. Formally, the statistical independence of different sources of constraints means that the covariance matrix of joint constraint $\mathbf{f}^* = [\mathbf{f}_1^*; \mathbf{f}_2^*; \cdots; \mathbf{f}_M^*]$ has the following structure

$$ \mathbf{C}_{f^*} = \begin{bmatrix} \mathbf{C}_1 & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{C}_2 & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{C}_M \end{bmatrix}, \qquad (42) $$

where $\mathbf{C}_{f^*}$ indicates the covariance matrix of the ith constraint ($\mathbf{f}_i^*$). Following the expressions of $\mathbf{f}_i^*$ and $\mathbf{C}_{f_i^*}$, the probability distribution function (PDF) of joint data ($1 \le i \le M$) can be derived by multiplying PDFs of data from all M sources of constraints, namely,

$$ P(\mathbf{f}(\mathbf{x})|\mathbf{f}^*) = \prod_{i=1}^{M} P(\mathbf{f}_i(\mathbf{x})|\mathbf{f}_i^*) \sim \exp\left\{ -\frac{1}{2} \sum_{i=1}^{M} [\mathbf{f}_i(\mathbf{x}) - \mathbf{f}_i^*]^T (\mathbf{C}_i)^{-1} [\mathbf{f}_i(\mathbf{x}) - \mathbf{f}_i^*] \right\}. \quad (43) $$

Further introducing the weight matrix (W), the objective cost function to be minimized has a quadratic form, namely,

$$ \Psi_{\text{total}}(\mathbf{x}) = \sum_{i=1}^{M} \gamma_i \Psi_i(\mathbf{x}), \qquad (44) $$

where

$$ \Psi_i(\mathbf{x}) = \frac{1}{2} [\mathbf{f}_i(\mathbf{x}) - \mathbf{f}_i^*]^T \mathbf{W}_i^{-1} [\mathbf{f}_i(\mathbf{x}) - \mathbf{f}_i^*], \qquad (45) $$

with

$$ \mathbf{W}_i = \frac{1}{\varepsilon_i^2} \mathbf{C}_i, \qquad (46) $$

$$ \gamma_i = \frac{\varepsilon_i^2}{\varepsilon_i^2 + \varepsilon_i^2}, \qquad (47) $$

In the above equations, $\varepsilon_i^2$ is the first diagonal element of $\mathbf{C}_i$ (i.e., $\varepsilon_i^2 = \{C_i\}_{11}$) and the Lagrange factor $\gamma_i$ weights the contribution of each type of constraint with respective to the first one ($\gamma_1 = 1$).

Minimization of $\Psi_{\text{total}}(\mathbf{x})$ in Eq. 44 means its gradient with respect to the solution $\mathbf{x}$ approaches zero, such that

$$ \nabla \Psi_{\text{total}}(\mathbf{x}) = \sum_{i=1}^{M} \gamma_i \nabla \Psi_i(\mathbf{x}) = 0, \qquad (48) $$

which can be ensured by enforcing the gradient of all components approach zero, namely

$$ \nabla \Psi_i(\mathbf{x}) = \mathbf{K}_i^T \mathbf{W}_i^{-1} (\mathbf{f}_i(\mathbf{x}) - \mathbf{f}_i^*) = 0, \qquad (49) $$

where $\mathbf{K}_i$ is the Jacobian matrix containing the derivatives of the ith type of constraint with respect to the retrieval parameters. To approach the solution iteratively, we replace $\mathbf{x}$ by $\mathbf{x} - \Delta\mathbf{x}$ in Eq. 49 and substitute

$$ \mathbf{f}_i(\mathbf{x} - \Delta\mathbf{x}) = \mathbf{f}_i(\mathbf{x}) - \mathbf{K}_i \Delta\mathbf{x}, \qquad (50) $$

into it. This results in

$$ \nabla \Psi_i(\mathbf{x} - \Delta\mathbf{x}) = \mathbf{K}_i^T \mathbf{W}_i^{-1} (\mathbf{f}_i(\mathbf{x}) - \mathbf{K}_i \Delta\mathbf{x} - \mathbf{f}_i^*) = 0, \qquad (51) $$

or equivalently,

$$ (\mathbf{K}_i^T \mathbf{W}_i^{-1} \mathbf{K}_i) \Delta \mathbf{x} = \mathbf{K}_i^T \mathbf{W}_i^{-1} (\mathbf{f}_i(\mathbf{x}) - \mathbf{f}_i^*) = \nabla \Psi_i(\mathbf{x}). \qquad (52) $$

The Jacobian matrix $\mathbf{K}_i$ in Eqs 49–52 consists of the derivative of the mth observational or a priori data with respect to the nth unknown,

$$ K_{i,(m,n)} = \frac{\partial f_{i,m}}{\partial x_n}. \qquad (53) $$

More explicit evaluation of $\mathbf{f}_i(\mathbf{x})$, $\mathbf{f}_i^*$ and $\mathbf{K}$ and $\mathbf{W}$ matrices for all types of constraints is discussed in the following subsections.

### Observational Constraints (i = 1)

Without the loss of generality, we assume a polarimeter measures multi-angular radiance (L) across multiple spectral channels. Neglecting the contribution from circular polarization, the polarimetric components Q/L and U/L can be combined into DOLP via DOLP = (Q² + U²)¹/²/L. Lidar observations include extinction (α), backscatter or attenuated backscatter (β) and depolarization ratio (δ). Arranging all types of observations into a single column vector we have

$$ \mathbf{f}_{i=1}^{*}(\mathbf{x}) = (\mathbf{f}_{1,L}^{*}; \mathbf{f}_{1,DOLP}^{*}; \mathbf{f}_{1,\alpha}^{*}; \mathbf{f}_{1,\beta}^{*}; \mathbf{f}_{1,\delta}^{*}) = [y_1; y_2; \cdots; y_{N_f}], \qquad (54) $$

where the total number of signals $N_f = N_L + N_{DOLP} + N_\alpha + N_\beta + N_\delta$, where $N_\alpha$ is the number of “x” type of signals.

The aerosol and ocean properties in the state vector are adjusted so that the model prediction of these types of signals $\mathbf{f}_1(\mathbf{x}) = [\mathbf{f}_{1,L}; \mathbf{f}_{1,DOLP}; \mathbf{f}_{1,\alpha}; \mathbf{f}_{1,\beta}; \mathbf{f}_{1,\delta}]$ fit the observational constraints $\mathbf{f}_1^*(\mathbf{x})$. The calculation of $\mathbf{f}_1(\mathbf{x})$ is introduced in Forward Model for Light Propagation in a Coupled Atmosphere-Ocean System. The Jacobian matrix $\mathbf{K}_1$ consists of first order partial derivatives with respect to all elements in the state vector $\mathbf{x}$, namely,

$$ \mathbf{K}_1 = [ \mathbf{K}_L \quad \mathbf{K}_{DOLP} \quad \mathbf{K}_\alpha \quad \mathbf{K}_\beta \quad \mathbf{K}_\delta ]^T, \qquad (55) $$

where $\mathbf{K}_L$, $\mathbf{K}_{DOLP}$, $\mathbf{K}_\alpha$, $\mathbf{K}_\beta$, and $\mathbf{K}_\delta$ are the Jacobian matrices containing derivatives of L, DOLP, α, β and δ with respect to aerosol parameters $x_j$ ($1 \le j \le N$), respectively. Accounting for all N parameters in the state vector (see Eq. 39), $\mathbf{K}_L$ can be expanded into the following matrix form:
$$ \mathbf{K}_{\mathrm{L}} = \begin{bmatrix} \frac{\partial L_1}{\partial x_1} & \frac{\partial L_1}{\partial x_2} & \cdots & \frac{\partial L_1}{\partial x_N} \\ \frac{\partial L_2}{\partial x_1} & \frac{\partial L_2}{\partial x_2} & \cdots & \frac{\partial L_2}{\partial x_N} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial L_{N_L}}{\partial x_1} & \frac{\partial L_{N_L}}{\partial x_2} & \cdots & \frac{\partial L_{N_L}}{\partial x_N} \end{bmatrix}, \quad (56) $$

where $N_L$ depends on the number of views angles ($N_{view}$) and the number of channels ($N_{λ, L}$) of a polarimeter. Assuming the same number of view angles for all spectral bands, then we have $N_L = N_{λ,L} \times N_{view}$. Replacing “L” in the above equation by DOLP gives the Jacobian matrix $K_{DOLP}$ and in this case the total number of DOLP signals depends on the number of polarimetric bands ($N_{λ,DOLP}$) so that $N_{DOLP} = N_{λ,DOLP} \times N_{view}$. Calculation of the derivatives in Eq. 56 can proceed numerically by use of finite difference method, or analytically (and more efficiently) if the RT model is linearized.

The Jacobian matrix for extinction coefficient equals to

$$ \mathbf{K}_{\alpha} = \begin{bmatrix} \mathbf{K}_{\alpha,l=1} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{K}_{\alpha,l=2} & \cdots & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{K}_{\alpha,l=N_{\mathrm{Layer}}} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0}_{\mathrm{Ocean}} \end{bmatrix}, \quad (57) $$

where $N_{Layer}$ is the number of layers in which lidar signals are resolved. Assuming extinction $\alpha$ is provided at $M_\alpha$ wavelengths, the layer resolved Jacobian matrix ($K_{\alpha,l}$) for extinction coefficient equals to,

$$ \mathbf{K}_{\alpha,l} = \begin{bmatrix} \frac{\partial \alpha_l(\lambda_1)}{\partial x_{l,1}} & \frac{\partial \alpha_l(\lambda_1)}{\partial x_{l,2}} & \cdots & \frac{\partial \alpha_l(\lambda_1)}{\partial x_{l,N}} \\ \frac{\partial \alpha_l(\lambda_2)}{\partial x_{l,1}} & \frac{\partial \alpha_l(\lambda_2)}{\partial x_{l,2}} & \cdots & \frac{\partial \alpha_l(\lambda_2)}{\partial x_{l,N}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \alpha_l(\lambda_{M_\alpha})}{\partial x_{l,1}} & \frac{\partial \alpha_l(\lambda_{M_\alpha})}{\partial x_{l,2}} & \cdots & \frac{\partial \alpha_l(\lambda_{M_\alpha})}{\partial x_{l,N}} \end{bmatrix}, \quad (58) $$

The Jacobian matrix for backscatter coefficient and attenuated backscatter coefficients is evaluated by,

$$ \mathbf{K}_{\beta} = \begin{bmatrix} \mathbf{K}_{\beta,l=1} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{K}_{\beta,l=2} & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{K}_{\beta,l=N_{\mathrm{Layer}}} \end{bmatrix}, \quad (59) $$

and

$$ \mathbf{K}_{\beta} = \begin{bmatrix} \mathbf{K}_{\beta,l=1} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{K}_{\beta,l=1} & \mathbf{K}_{\beta,l=2} & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{K}_{\beta,l=1} & \mathbf{K}_{\beta,l=2} & \cdots & \mathbf{K}_{\beta,l=N_{\mathrm{Layer}}} \end{bmatrix}, \quad (60) $$

respectively. Assuming $\beta$ is provided at $M_\beta$ wavelengths, the layer resolved Jacobian matrix for backscatter coefficient ($K_{\beta,l}$) can be expanded into the following specific form,

$$ \mathbf{K}_{\beta,l} = \begin{bmatrix} \frac{\partial \beta_l(\lambda_1)}{\partial x_{l,1}} & \frac{\partial \beta_l(\lambda_1)}{\partial x_{l,2}} & \cdots & \frac{\partial \beta_l(\lambda_1)}{\partial x_{l,N}} \\ \frac{\partial \beta_l(\lambda_2)}{\partial x_{l,1}} & \frac{\partial \beta_l(\lambda_2)}{\partial x_{l,2}} & \cdots & \frac{\partial \beta_l(\lambda_2)}{\partial x_{l,N}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \beta_l(\lambda_{M_\beta})}{\partial x_{l,1}} & \frac{\partial \beta_l(\lambda_{M_\beta})}{\partial x_{l,2}} & \cdots & \frac{\partial \beta_l(\lambda_{M_\beta})}{\partial x_{l,N}} \end{bmatrix}, \quad (61) $$

where $\beta = \beta_{l,a}$ if we use the derived product of particulate backscatter for the $l$th layer from lidar measurement, and $\beta = \beta_{l,att}$ if we used direct measured attenuated backscatter signals, where $\beta_{l,att}$ is the product of the backscatter (contributed by both molecules and aerosols) and the two-way transmission of the atmospheric volume between the lidar and the backscatter volume in question (see Eq. 36).

The Jacobian matrix for depolarization ratio ($\delta$) equals to

$$ \mathbf{K}_{\delta} = \begin{bmatrix} \mathbf{K}_{\delta,l=1} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{K}_{\delta,l=2} & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{K}_{\delta,l=N_{\mathrm{Layer}}} \end{bmatrix}. \quad (62) $$

Assuming extinction $\delta$ is provided at $M_\delta$ wavelengths, the layer resolved Jacobian matrix for depolarization ratio ($K_{\delta,l}$) equals to

$$ \mathbf{K}_{\delta,l} = \begin{bmatrix} \frac{\partial \delta_l(\lambda_1)}{\partial x_{l,1}} & \frac{\partial \delta_l(\lambda_1)}{\partial x_{l,2}} & \cdots & \frac{\partial \delta_l(\lambda_1)}{\partial x_{l,N}} \\ \frac{\partial \delta_l(\lambda_2)}{\partial x_{l,1}} & \frac{\partial \delta_l(\lambda_2)}{\partial x_{l,2}} & \cdots & \frac{\partial \delta_l(\lambda_2)}{\partial x_{l,N}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \delta_l(\lambda_{M_\delta})}{\partial x_{l,1}} & \frac{\partial \delta_l(\lambda_{M_\delta})}{\partial x_{l,2}} & \cdots & \frac{\partial \delta_l(\lambda_{M_\delta})}{\partial x_{l,N}} \end{bmatrix}. \quad (63) $$

where $\delta_l$ is contributed by both molecules and aerosols in $l$th layer and is equal to particulate depolarization ratio ($\delta = \delta_{l,a}$) if the derived aerosol depolarization ratio is used directly.

The overall covariance matrix incorporates the submatrices from all types of signals, namely,

$$ \mathbf{W}_1 = \begin{bmatrix} \mathbf{W}_L & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{W}_{\mathrm{DOLP}} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{W}_\alpha & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{W}_\beta & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{W}_\delta \end{bmatrix}, \quad (64) $$

If absolute measurement uncertainty is used and correlation among signals are neglected, then $W_x$ is a diagonal matrix. Taking $j$th signal ($x_j$) with relative uncertainty $\sigma_j$ as an example, we have $W_x, \mathbf{jj} = \sigma_j^2$. If relative uncertainty ($\varepsilon_j$) is used, then we have $W_x, \mathbf{jj} = \varepsilon_j^2 x_j^2$. Under certain circumstances signals may differ in magnitude. To balance their contributions to the cost function which reflects the quality of fitting, observations are often transferred to logarithmic space for fitting. In this case, we have $W_x, \mathbf{jj} = \varepsilon_j^2$.
## A Priori Constraints (i = 2)

Starting from Eq. 41, the a priori constraint can be explicitly expressed as (Dubovik, 2004),

$$ \mathbf{f}_{i=2}^{*} = \mathbf{x}^{a \text{ priori}} = \mathbf{x} + \Delta \mathbf{x}^{a \text{ priori}}. \quad (65) $$

Moreover, $\mathbf{K}_{i=2} = I$ (identity matrix) and $\mathbf{W}_{i=2} = (1/\sigma_{a_{i1}}^2)\mathbf{C}_a^*$ in Eqs 49–52. In explicit form, $\mathbf{W}_{i=2}$ is constructed from estimated range of each parameter relative to the first one (Dubovik and King, 2000),

$$ \mathbf{W}_{i=2} = \mathbf{W}_a = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & \frac{(x_{2,\text{max}} - x_{2,\text{min}})^2}{(x_{1,\text{max}} - x_{1,\text{min}})^2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{(x_{N,\text{max}} - x_{N,\text{min}})^2}{(x_{1,\text{max}} - x_{1,\text{min}})^2} \end{bmatrix}. \quad (66) $$

where $[x_{i,min}, x_{i,max}]$ are the lower and upper bounds of ith retrieval parameters, respectively. For retrieving combined RSP and HSRL-2 observations, their estimates are specified in Table 1.

## Smoothness Constraints (i = 3)

The third type of constraint reflects the smooth variation of a certain type of parameter as a function of altitude. we have

$$ \mathbf{f}_3^* = \mathbf{0}^* = \mathbf{S}_{i=3, m} \mathbf{x} + \Delta_g^*, \quad (67) $$

where $\mathbf{S}_{i=3,m}$ is the differentiation matrix of mth order for third type of constraint and the delta term $\Delta_g^*$ accounts for the uncertainties associated with the assumption of smooth variation of an aerosol parameter to the mth order. The error of such an assumption is accounted via the cost function component given in Eq. 75 later.

When a retrieval parameter varies smoothly as a function of some variable (e.g., altitude) $z$, it is assumed to be locally approximated by a smooth function $g(z)$, such as a constant, a line, a parabola, etc. With a polynomial form, the mth derivative approaches zero (Dubovik et al., 2011), such that,

$$ \left\{ S_{i,m}(z)z = \frac{d^m g(z)}{dz^m} = 0 \Rightarrow \begin{cases} g_{m=1}(z) = \text{const} \\ g_{m=2}(z) = Az + B \\ g_{m=3}(z) = Az^2 + Bz + C \end{cases} \right. . \quad (68) $$

For a discretized grid of the variable $z$, the explicit form of $S_{i,m}z$, is,

$$ S_{i,m}(z_i)z_j = \begin{cases} \frac{dg}{dz} \approx \frac{\Delta_1^1 g}{\Delta_1(z)} = \frac{g(z_{j+1}) - g(z_j)}{\Delta_1(z_j)}, & \text{for } m=1 \\ \frac{d^m g}{dz^m} \approx \frac{\Delta_m^m g}{\Delta_m(z)} = \frac{\Delta_{m-1}^{m-1} g(z_{j+1}) / \Delta_{m-1}(z_{j+1}) - \Delta_{m-1}^{m-1} g(z_j) / \Delta_{m-1}(z_j)}{[\Delta_{m-1}(z_j) + \Delta_{m-1}(z_{j+1})]/2}, & \text{for } m \ge 2 \end{cases} \quad (69) $$

Taking the orders of difference $m = 1$ and 2 as examples, we have

$$ \left\{ \begin{aligned} \Delta_{m=1}(z_j) &= z_{j+1} - z_j \\ \Delta_{m=2}(z_j) &= [\Delta_1(z_j) + \Delta_1(z_{j+1})]/2 \end{aligned} \right. . \quad (70) $$

Application of the above equation to L discretized grids $r_j$ (namely, $1 \le j \le L$) leads to

$$ \mathbf{f}_i(\mathbf{x}) = \mathbf{S}_{i,m}\mathbf{x}, \quad (71) $$

so that by invoking Eq. 53,

$$ \mathbf{K}_i = \mathbf{S}_{i,m}, \quad (72) $$

where the matrix $\mathbf{S}_{i,m}$ is evaluated by,

$$ \mathbf{S}_{i,m=1} = \begin{bmatrix} \frac{1}{\Delta_1(1)} & -\frac{1}{\Delta_1(1)} & 0 & \cdots & 0 \\ 0 & \frac{1}{\Delta_1(2)} & -\frac{1}{\Delta_1(2)} & \cdots & 0 \\ 0 & 0 & \cdots & \ddots & 0 \\ 0 & 0 & \cdots & \frac{1}{\Delta_1(L-1)} & \frac{1}{\Delta_1(L-1)} \end{bmatrix}, \quad (73) $$

and

$$ \mathbf{S}_{i,m=2} = \begin{bmatrix} \frac{2}{\Delta_1(1)[\Delta_1(1) + \Delta_1(2)]} & -\frac{2}{\Delta_1(1)[\Delta_1(2)]} & \frac{2}{\Delta_1(2)[\Delta_1(1) + \Delta_1(2)]} & 0 & 0 & \cdots & 0 \\ \frac{2}{\Delta_1(2)[\Delta_1(2) + \Delta_1(3)]} & -\frac{2}{\Delta_1(2)[\Delta_1(3)]} & \frac{2}{\Delta_1(3)[\Delta_1(2) + \Delta_1(3)]} & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots & & \vdots \\ \frac{2}{\Delta_1(L-1)[\Delta_1(L-1) + \Delta_1(L)]} & -\frac{2}{\Delta_1(L-1)[\Delta_1(L)]} & \frac{2}{\Delta_1(L)[\Delta_1(L-1) + \Delta_1(L)]} & \cdots & \cdots & & \frac{2}{\Delta_1(L)[\Delta_1(L)]} \end{bmatrix}. \quad (74) $$

The same principle applied to higher orders of difference ($m > 2$) ensures a smooth curve with $d^m g(z)/dz^m = 0$. Substitution of Eqs 71 and 72 and $\mathbf{f}_i^* = \mathbf{0}^*$ into Eq. 45 gives,

$$ \Psi_i(\mathbf{x}) = \frac{1}{2} \mathbf{x}^\mathrm{T} \mathbf{S}_{i,m}^\mathrm{T} \mathbf{W}_{i,m}^{-1} \mathbf{S}_{i,m} \mathbf{x}, \quad (75) $$

where the weighting matrix $\mathbf{W}$ has the following diagonal terms,

$$ \{\mathbf{W}_{i,m}\}_{jj} = \frac{1}{\Delta_m(z_j)}, \quad (76) $$

and $\Delta_m(z_j)$ is specified in Eq. 70 for $m = 1$ and 2, and can be generalized to an arbitrary higher order.

Substitution of Eqs 71 and 72 and $\mathbf{f}_i^* = \mathbf{0}^*$ into Eq. 52 gives

$$ (\mathbf{S}_{i,m}^\mathrm{T} \mathbf{W}_{i,m}^{-1} \mathbf{S}_{i,m}) \Delta \mathbf{x} = \mathbf{S}_{i,m}^\mathrm{T} \mathbf{W}_{i,m}^{-1} (\mathbf{S}_{i,m} \mathbf{x}). \quad (77) $$

Further defining

$$ \mathbf{\Omega}_i = \mathbf{S}_{i,m}^\mathrm{T} \mathbf{W}_{i,m}^{-1} \mathbf{S}_{i,m}, \quad (78) $$

we have

$$ \mathbf{\Omega}_i \Delta \mathbf{x} = \mathbf{\Omega}_i \mathbf{x}. \quad (79) $$

## Construction of Overall Equation System

Accounting for the above three types of constraints, the solution to the non-linear minimization of Eq. 44 can be approached iteratively. At iteration $q$, the solution is updated as,

$$ \mathbf{x}_{q+1} = \mathbf{x}_q - \Delta \mathbf{x}_q, \quad (80) $$
where $\Delta\mathbf{x}_q$ is obtained by accounting for all three types of constraints derived in *Observational Constraints* (i = 1), *A Priori Constraints* (i = 2), and *Smoothness Constraints* (i = 3),

$$ \mathbf{A}_q \times \Delta\mathbf{x}_q = \nabla\Psi_{\text{total}}(\mathbf{x}_q). \quad (81) $$

As the explicit form, we have

$$ \mathbf{A}_q = \mathbf{K}_{1,q}^{\mathrm{T}} \mathbf{W}_1^{-1} \mathbf{K}_{1,q} + \gamma \mathbf{\Omega} + \gamma_a \mathbf{W}_a^{-1}, \quad (82) $$

$$ \nabla\Psi_{\text{total}}(\mathbf{x}_q) = \mathbf{K}_{1,q}^{\mathrm{T}} \mathbf{W}_1^{-1} [\mathbf{f}_1(\mathbf{x}_q) - \mathbf{f}_1^*] + \gamma \mathbf{\Omega} \mathbf{x}_q + \gamma_a \mathbf{W}_a^{-1} (\mathbf{x}_q - \mathbf{x}^a_{\text{priori}}), \quad (83) $$

where “$\gamma\mathbf{\Omega}$” incorporates vertical and spectral smoothness constraints, namely,

$$ \gamma\mathbf{\Omega} = \begin{bmatrix} \gamma_h \mathbf{\Omega}_h & 0 \\ 0 & \mathbf{0}_{\text{ocean}} \end{bmatrix}, \quad (84) $$

where the smoothness matrix $\gamma_h\mathbf{\Omega}_h$ acts on the vertical variations of aerosol parameters; and the multipliers $\gamma_h$ control the strength of the constraints. The explicit forms of $\gamma_h\mathbf{\Omega}_h$ are given by.

$$ \gamma_h \mathbf{\Omega}_h = \begin{pmatrix} \gamma_{h,n=1} \mathbf{\Omega}_{n=1} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \gamma_{h,n=2} \mathbf{\Omega}_{n=2} & \cdots & \mathbf{0} \\ \cdots & \cdots & \ddots & \cdots \\ \mathbf{0} & \mathbf{0} & \cdots & \gamma_{h,N_{\text{aer}}} \mathbf{\Omega}_{N_{\text{aer}}} \end{pmatrix}, \quad (85) $$

where the diagonal terms correspond to the constraints imposed on the $N_{\text{aer}}$ aerosol parameters.

## Convergence Criteria

Ideally, a retrieval is deemed successful when the minimization of the cost function is achieved in such a way that ensures

$$ \Psi_{\text{total}} \le (N_f + N_c + N_{a*} - N_a) \varepsilon_f^2, \quad (86) $$

where $N_f$ is the number of observations, $N_c$ is the number of constraints imposed on retrieval (including vertical smoothness constraints), $N_a$ is the number of retrieval parameters, and $N_{a*}$ is the number of a priori estimates of parameters; $\varepsilon_f^2$ is the expected variance due to measurement errors. In practice, forward RT modeling error and other unmodeled effects can impede the cost function being satisfied by the form of Eq. 86. Taking such a practical situation into account, the retrieval is also terminated when the relative difference of fitting residuals with solutions from two successive iterations drops below a user-specified threshold value, $\varepsilon_c^2$, namely,

$$ \frac{|\Psi_{\text{total}}(\mathbf{x}_{q+1}) - \Psi_{\text{total}}(\mathbf{x}_q)|}{\Psi_{\text{total}}(\mathbf{x}_q)} \le \varepsilon_c^2. \quad (87) $$

## Determination of Lagrange Multipliers

The Lagrange multipliers reflect the strength of smoothness constraints for a given parameter to be retrieved. Each multiplier is defined as,

$$ \gamma_i = \varepsilon_{i=1}^2 / \varepsilon_i^2, \quad (88) $$

where $\varepsilon_i^2$ are the first diagonal elements of the covariance matrices corresponding to ith type of constraints. To estimate $\varepsilon_i^2$ for a given parameter to be retrieved ($f=x$) which is a function of $z$, the most unsmooth known solution $f^{\text{us}}(z)$ over the target domain is used as suggested by Dubovik and King (2000) namely,

$$ \varepsilon^2 = \int_{z_{\min}}^{z_{\max}} \left( \frac{d^m[f^{\text{us}}(z)]}{d^m z} \right)^2 dz, \quad (89) $$

where $z_{\min}$ and $z_{\max}$ specify the lower and upper bound of $z$. For $\gamma_{i=3}$ is evaluated using most unsmooth known solution of parameter $x$ over the target domain. In the current study, the vertical variations of extinction and backscatter from lidar measurements are referred to initialize the Lagrange multipliers for profiled aerosol properties.

Indeed, the Lagrange multipliers for practical implementation has to be modified to account for possible redundancy of the measured and a priori data and to reduce the weight of a priori constraints during an optimization process (Xu et al., 2019):

$$ \gamma_i^{\text{Final}} = \frac{N_f}{N_i} \frac{\tilde{\varepsilon}_f^2}{\varepsilon_f^2} \gamma_i, \quad (90) $$

where

$$ \tilde{\varepsilon}_f^2(\mathbf{x}_q) \approx \frac{\Psi_{\text{total}}(\mathbf{x}_q)}{(N_f + N_c + N_{a*} - N_a)\varepsilon_f^2}. \quad (91) $$

For retrieving a combined set of RSP and HSRL-2 observations, first estimates of Lagrange parameters $\gamma_i$ are specified in Table 1 and then scaled by the ratio of the cost functions associated with the updated and the previous iterative solutions during optimization. Typically, retrieval takes 8-10 iterations to meet the convergence criteria specified *Convergence Criteria*.

## Retrieval Error Estimate

Random errors ($\Delta\mathbf{y}_{\text{rand}}^{\text{obs}}$) and systematic errors ($\Delta\mathbf{y}_{\text{syst}}$) are two major error types. The former is contributed by measurements while the latter is contributed by both measurement ($\Delta\mathbf{y}_{\text{syst}}^{\text{obs}}$) and forward model ($\Delta\mathbf{y}_{\text{model}}$). Following Dubovik's error estimate methodology [2004], we also consider a priori smoothness constraints regarding the vertical distribution of aerosol parameters, and a priori estimate of the solution in retrieval error estimate. The overall covariance matrix for error estimate for combined lidar-polarimeter inversion is expressed as,

$$ \mathbf{C}_{\Delta x, \text{syst}} = \Delta \mathbf{x}_{\text{syst}} (\Delta \mathbf{x}_{\text{syst}})^{\mathrm{T}} + \mathbf{C}_{\Delta x, \text{rand}}, \quad (92) $$

where $\Delta\mathbf{x}_{\text{syst}}$ is approximated by

$$ \Delta\mathbf{x}_{\text{syst}} = \mathbf{A}^{-1} \nabla\Psi, \quad (93) $$

and covariance matrix of the retrieval solution contributed by the random error is evaluated by

$$ \mathbf{C}_{\Delta x, \text{rand}} = \mathbf{A}^{-1} \varepsilon_{\text{rand}}^2, \quad (94) $$

where $\varepsilon_{\text{rand}}$ is the random noise, $\mathbf{A}$ is computed at the retrieved solution $\mathbf{x}$ via Eq. 82 and $\nabla\Psi$ is given by
$$ \nabla \Psi = \mathbf{J}^T \mathbf{W}_f^{-1} \Delta \mathbf{y}_{\text{syst}} + \gamma \Omega \mathbf{x}^{\text{true}} + \gamma_a \mathbf{W}_a^{-1} (\mathbf{x}^{\text{true}} - \mathbf{x}^*), \quad (95) $$

where the total systematic error is expressed as

$$ (\Delta \mathbf{y}_{\text{syst}})^2 = (\Delta \mathbf{y}_{\text{syst}}^{\text{obs}})^2 + (\Delta \mathbf{y}_{\text{model}})^2. \quad (96) $$

By performing the above error analysis, it is assumed that the modeling error and the systematic errors of an instrument in Eq. 28 are well characterized. Indeed, their quantification is rather difficult due to the complexity of error analysis for a sensor that consists of multiple optical units as well as the under-characterization of multiple sources of modeling errors.

For practicality, $\Delta\mathbf{x}_{\text{syst}}$ is calculated by assuming

$$ \Delta \mathbf{y}_{\text{syst}} = \mathbf{f}_1(\mathbf{x}^{\text{retrieved}}) - \mathbf{f}_1^*. \quad (97) $$

so that Eq. 93 becomes

$$ \Delta \mathbf{x}_{\text{syst}} = \mathbf{A}^{-1} \left[ \begin{aligned} & \mathbf{K}_1^T \mathbf{W}_1^{-1} [\mathbf{f}_1(\mathbf{x}^{\text{retrieved}}) - \mathbf{f}_1^*] + \gamma \Omega \mathbf{x}^{\text{true}} \\ & + \gamma_a \mathbf{W}_a^{-1} (\mathbf{x}^{\text{true}} - \mathbf{x}^{\text{a priori}}) \end{aligned} \right]. \quad (98) $$

As the retrieved solution is the closest estimate of the “true” solution $\mathbf{x}^{\text{true}}$, we assume $\mathbf{x}^{\text{true}} = \mathbf{x}^{\text{retrieved}}$. When a priori $\mathbf{x}^{\text{a priori}}$ is unavailable, we further assume $\mathbf{x}^{\text{a priori}} = \mathbf{x}^{\text{true}}$. Under these assumptions, the last term on the right-hand-side of Eq. 95 disappears (namely $\gamma_a \mathbf{W}_a^{-1} (\mathbf{x}^{\text{true}} - \mathbf{x}^{\text{a priori}}) = 0$). Note that Eq. 94 treats measurement errors as random. Indeed, these errors are already contained in the observation $\mathbf{f}_1$ used by Eq. 98. Therefore, by implementing Eqs 92, 94 and 98, it is possible that errors are double-counted in the case that they bias the solution in the same direction as the modeling errors, resulting in a conservative error estimate.

To estimate errors for functions of the retrieved parameters (namely $y = f(\mathbf{x})$, and $y$ can be vertical layer resolved SSA and AOD, vertical layer resolved extinction and scattering coefficients, accumulated AOD from the TOA, effective radius, etc.), the chain rule is applied so that we have the following matrix form,

$$ \Delta y = \sqrt{\mathbf{K}_y^T \mathbf{C}_{\Delta x} \mathbf{K}_y}, \quad (99) $$

where $\mathbf{K}_y$ denotes the Jacobian array containing derivatives of $y$ with respect to all retrieved parameters $x$ that are relevant to the calculation of $y$.

Summing over the contribution by both systematic and random errors, the overall error is estimated by taking the diagonal term of the retrieval error covariance matrix, namely

$$ \Delta x_j = \sqrt{(\mathbf{C}_{\Delta x})_{jj}}, \quad (100) $$

where $\mathbf{C}_{\Delta x}$ is given in Eq. 92.

Our error estimates involve the use of Jacobians and assume 1) the retrieved solution is representative of the solution space; and 2) retrieval errors are linearly proportional to the measurement errors. These assumptions can be problematic in situations where the model and/or observation errors are large. Moreover, the observation covariance matrix assumes no correlations among errors of observational signals so that its off-diagonal terms are

zero in Eq. 57. In practice, however, observation errors can highly correlated between neighboring spectral bands and/or angular measurements. In addition, the a priori covariance matrix also assumes no error dependency between state vector elements so that its off-diagonal terms in Eq. 66 are zero. These are potential error sources of our error estimate model.

# RETRIEVAL IMPLEMENTATION AND VALIDATION

## Observations From HSRL-2 and RSP

Following the algorithm formulation in the previous sections, the retrieval approach is applied to combined sets of polarimetric and lidar observations acquired by RSP and HSRL-2 on September 12, 2016 during the ORACLES field campaign. The polarimeter signals include radiance and DOLP and the lidar signals include particulate extinction coefficient, backscatter, and depolarization ratio.

During the ORACLES mission, RSP and HSRL-2 were aboard NASA’s ER-2 aircraft flying at an altitude of 20 km. NASA’s P-3 Orion aircraft flew underneath the ER-2 at up to 6 km in altitude. It carried a set of instrumentation, including the Hawaii Group for Environmental Aerosol Research (HiGEAR, (McNaughton et al., 2009)) aerosol in situ sampling instruments and NASA Ames Spectrometer for Sky-Scanning, Sun-Tracking Atmospheric Research (4-STAR, (Dunagan et al., 2013)) that combines airborne Sun tracking and sky scanning with grating spectroscopy to improve knowledge of atmospheric constituents. HiGEAR’s suite of in situ sampling instruments measure absorption and scattering coefficients under dry conditions, which can be used to derive aerosol SSA. 4-STAR measures accumulated AOD.

We selected the scenes that met the following criteria: 1) collocated measurements by RSP and HSRL-2 onboard NASA’s ER-2 high altitude aircraft are available; 2) collocated AOD, SSA, and aerosol size distribution data for retrieval validation provided by 4-STAR and HiGEAR are available; and 3) the measurements were made at cloud-free condition, as identified from HSRL-2’s backscatter measurements. Figure 3 provides a curtain of HSRL-2 aerosol backscatter at 532 nm along the ER-2 flight on September 12, 2016. This is the typical aerosol scene from the ORACLES campaign over the southeast Atlantic Ocean. A layer of elevated smoke aerosol was originated from biomass burning and transported to the ocean. After mixing with marine sea salt aerosol (dark red in the lidar curtain) from the top of the boundary layer, the mixture complexity provides us an ideal scene to test our algorithm. As indicated by the white line in the upper right part of the figure, the lidar and polarimeter data meeting these criteria were acquired from 14:25 to 14:34 UTC. During this timeframe, the solar zenith angle is around 52°, and the ER-2 aircraft flew a distance about 100 km from (-18.73°N, 8.72°E) to (-19.39°N, 9.39°E). A total of 50 sets of observations with collocated measurements by HSRL-2 and RSP data were collected to test the combined retrieval. Moreover, 4-STAR measurements of accumulated AODs at different altitudes under 6 km (as traversed by the P-3) and SSA derived from
<figure><img src="image_5.png" /></figure>

HiGEAR products of absorption and scattering coefficients, and number weighted aerosol size distribution are available in this time frame, as indicated in **Figure 4**. They are used for retrieval validation.

The RSP instrument measures the radiance and polarimetric components Q and U in nine wavelength bands from UV (410 nm) to short-wave infrared (2,250 nm) (Cairns et al., 1999). If RSP's measurements are collocated with HSRL-2's (lat, lon) within 1 km horizontal distance, then its angular data are used for combined retrieval. For the 50 sets of observations under test, a total of 152 view angles from -39° to 25° around nadir in the wavelengths 410, 470, 550, 670, and 865 nm are used. Note that due to the yaw of the ER-2 aircraft and high altitude, RSP's view angular range is reduced to [-39°, 25°] around nadir. These measurements are within ±35° around the principal plane. The relative uncertainty is 2% for intensity and the absolute uncertainty 0.002 for DOLP. HSRL-2 is the second generation NASA Langley airborne High Spectral Resolution Lidar (Hair et al., 2008). It uses the HSRL technique (Shipley et al., 1983) to independently measure aerosol backscatter and extinction at 355 and 532 nm. While HSRL-2 products of particulate extinction, backscatter, depolarization ratio are reported at 15 m vertical resolution, the true resolution is 315 m for extinction (α_a) and 15 m for backscatter (β_a) and depolarization (δ_a) in the ORACLES campaign (Burton et al., 2018). The lidar signals that we used in this study include aerosol extinction, backscatter, and depolarization ratio. In the HSRL-2 product, the aerosol extinction at 355 and 532 nm wavelengths are under the name "355_ext" and "532_ext", respectively. To preserve the independence of extinction data used for retrieval, we access α_a at 315 m resolution. The aerosol backscatter data at 355, 532, and 1,064 nm wavelength are under the name "355_bsc_Sa", "532_bsc_Sa", and "1,064_bsc_Sa", respectively in the HSRL-2 product. Their reported values are averaged at the same resolution as extinction for lidar ratio ("S_a") calculation. The depolarization ratios at the three lidar wavelengths are under

the name "355_aer Dep", "532_aer Dep", and "1,064_aer Dep" for the three wavelengths in the product. In computing the error covariance matrix, an uncertainty of 17 Mm⁻¹, 5%, and 10% are assumed for α_a, β_a, and δ_a, respectively. Within the 315 m layer, optical homogeneity is assumed. Such a layer is further divided into sublayers to ensure the radiative transfer modeling accuracy.

As the assumption, two aerosol species are retrieved from the combined use of HSRL-2 and RSP observations. Each species is characterized by an independent and layer resolved set of refractive index, log-normal volumetric size distribution, and spherical particle fraction. These quantities are initialized with values in **Table 1** before retrieval starts. As informed by the range of the imaginary part of the refractive index in **Table 1**, weakly to strongly absorbing fine mode aerosols, and non-absorbing to moderately absorbing coarse mode aerosols are assumed in retrieval. Moreover, the initial guess of the standard deviation of fine and coarse aerosol size distribution are informed by UHSAS measurements.

## Retrieval Validation Against 4-STAR and HiGEAR Measurements

As the first check, a retrieval is performed using RSP and HSRL-2 collocated measurement at 14:25:10 UTC, namely the first case among the total 50 sets of collocated observations. **Figure 5** displays the whole set of retrieved aerosol geophysical parameters as a function of altitude. These parameters include aerosol refractive index, volume concentration, geometric volume mean radius and logarithmic standard deviation of aerosol size distributions, and spherical particle fraction. These quantities are displayed in different panels. From the ranges of the retrieved real and imaginary parts of refractive index, the two aerosol species pre-assumed in our retrieval can be distinguished as smoke and sea salt particles.

Following the transportation of biomass burning smoke particles for a long distance (e.g., >800 km) from the
<figure><img src="image_6.png" /></figure>

**FIGURE 4** | Collocation of HSRL-2, 4-STAR, and the HiGEAR suite of instrumentation in taking measurements. The column aerosol optical depth accumulated from the position of the ER-2 aircraft to the location of P-3 aircraft calculated from HSRL-2 is used to validate the retrieval and to compare to the 4-STAR reference data. The aerosol SSA derived from extinction and scattering coefficients measured by HiGEAR suite of instrumentation is used to validate the retrieved SSA. HiGEAR also provides fine and coarse mode aerosol size distributions measured by UHSAS and APS instruments at dry condition, which are compared to the retrieved size distributions at ambient condition. Approximately at the altitude from 3 to 6 km, P-3 followed the trajectory of ER-2.

<figure><img src="image_7.png" /></figure>

**FIGURE 5** | Demonstration of aerosol geophysical property retrieval using a combined set of lidar (HSRL-2) and polarimeter (RSP) observations acquired at 14:25:10UTC on September 12, 2016 during the ORACLES field campaign. Upper left panel: real part of aerosol refractive index; Upper middle panel: imaginary part of aerosol refractive index; Upper middle right panel: volume concentration (unit: m³/m³); Bottom left panel: geometric volume mean radius (unit: µm) as a function of altitude; Bottom middle panel: logarithmic standard deviation of aerosol size distribution; Bottom right panel: spherical particle fraction. The horizontal bars in all panels indicate the retrieval errors. To estimate the retrieval uncertainties, we use the relative uncertainty 2% for RSP's intensity measurement and the absolute uncertainty 0.002 for DOLP measurement, and use the uncertainties 17 Mm⁻¹, 5%, and 10% for HSRL-2's particulate extinction, backscatter and depolarization, respectively.
continent, the volumetric concentration plot shows the co-existence of smoke and sea salt particles in two regions, namely from 0 to 1 km and from 2.5 to 5.2 km. Moreover, the sea salt particles present a larger size than that of smoke particles. The spherical particle fraction for smoke and sea salt particles are ~0.8 and >0.9, respectively, at the altitude close to the sea surface (with an retrieval uncertainty between 0.1 and 0.2). Increased nonspherical particle fraction for the sea salt particles is observed at a higher altitude from 2.5 to 5.2 km. This could be caused by crystallization as temperature, and relative humidity decreases at higher altitude. We have also noticed that the geometric mean radius of smoke particles is slightly lower in the boundary layer than in the free troposphere. In the context of the ORACLES campaign, the smoke particles in and above the boundary layer can be different in terms of their ages, transport, and microphysical properties. The smoke particles may be entrained into the PBL via large-scale re-circulation (Redemann et al., 2021) instead of directly sinking into PBL through subsidence as previously thought. Thus, the smoke aerosol properties in the PBL may be affected by their interactions with clouds or new particle formation processes. This is a likely cause of observing smaller smoke particles in the PBL than at high altitudes from the case study in Figure 5.

Calculated from the speciated aerosol abundance and intrinsic quantities displayed in Figures 5 and 6A displays the retrieved profiles of accumulated AOD (left panel) and SSA (right panel) at 532 nm (black curves), with comparison to the collocated HiGEAR measurement (red dot) at the altitude around 3 km. The retrieved AOD in the left panel of Figure 6A is accumulated AOD from TOA to different altitudes with 315 m resolution. Fairly consistent agreement between retrieved values and reference data derived from 4-STAR and HiGEAR measurements can be observed. The difference is within the estimated retrieval error for SSA and accumulated AOD. In Figure 6B, the retrieved fine mode (left panel) and coarse mode (right) aerosol size distributions (black curves) are compared to the in-situ HiGEAR measurement averaged within a vertical bin of 315 m (red curves). To compare to the number weighted aerosol size distribution measured by UHSAS and APS, the retrieved volumetric size distribution Eq. 2 needs a conversion. Specifically, the retrieved $r_m^{(s)}$ and $\sigma_m^{(s)}$ are converted to the geometric number mean radius $r_g^{(s)}$ and the standard deviation for sth aerosol species $\sigma_g^{(s)}$ by use of the following relationship (Grainger, 2017):

$$r_g^{(s)} = r_m^{(s)} e^{-3[\sigma_m^{(s)}]^2}, \qquad (101)$$

$$\sigma_g^{(s)} = \sigma_m^{(s)}. \qquad (102)$$

On such a basis, the retrieved number weighted size distribution is calculated as

$$\frac{dn^{(s)}(r)}{d\ln r} = \frac{1}{\sqrt{2\pi}\sigma_g^{(s)}} \exp\left[ -\frac{(\ln r - \ln r_g^{(s)})^2}{2(\sigma_g^{(s)})^2} \right]. \qquad (103)$$

While UHSAS reports number concentrations of aerosols at their geometric diameters, APS reports the concentration at aerodynamic diameters of aerosols. Therefore, the APS

aerodynamic particle diameter ($d_{ae}$) needs to be converted to the geometric diameter $d$. This is done via $d = d_{ae} \times (\chi/\rho_p)^{0.5}$, where $\rho_p$ is aerosol density and $\chi$ is the dynamic shape parameter. Assuming all aerosols measured by APS are sea salt particles, we adopt $\chi = 1.08$ and $\rho_p = 1.8 \text{ g/cm}^{-3}$ (Froyd et al., 2019). Then the conversion is determined as $d = d_{ae}/1.29$. By performing such a conversion, we neglect the impact of highly irregular particles on $\chi$ and assume equivalent spherical diameter for “d” which may cause bias to the estimate of $\chi$ (Reid et al., 2003). In addition, UHSAS and APS measure the number concentrations of fine and coarse mode particles in the ranges of geometric radius [0.03, 0.5] and aerodynamic radius [0.22, 8.1] $\mu$m, respectively. Apparently, the two instruments have an overlapped measurement range of geometric radius from 0.2 to 0.5 $\mu$m. Nevertheless, this range is accounted in determining $r_g^{(s)}$ and $\sigma_g^{(s)}$ for both fine and coarse aerosols. With these preconditions and conversions, Figure 6B compared the retrieved particle size distribution to HiGEAR’s in-situ measurement. While the difference is basically within the retrieval error estimated using the methodology described in Retrieval Error Estimate, the coarse mode aerosol size distribution is subjected to larger retrieval errors than the fine mode aerosol size distribution.

Along with the results demonstrated in Figures 5, 6, and 7A demonstrates the retrieval fit to HSRL-2’s measurements of vertical extinction coefficient, scattering coefficients, and depolarization ratio of aerosols. The vertical mean of absolute fitting error of particulate extinction coefficient is found to be 7.6 and 8.9 $Mm^{-1}$ for 355 and 532 nm, respectively, for altitudes below 5.2 km, and the vertical mean of absolute fitting error of particulate backscatter is 0.16, 0.07 and 0.14 $Mm^{-1}Sr^{-1}$ for 355, 532 and 1,064 nm, respectively. Noticeably, the fit to the particulate backscatter of 1,064 nm at ~1.4 km altitude has a relatively larger error. Indeed, to derive vertically resolved backscatter at 1,064 nm one needs to know the extinction at the same wavelength. Since HSRL-2 does not measure extinction at 1,064 nm directly, it is estimated from an assumed relationship with the measured lidar ratio at 532 nm. Though provided as a best guess, such an estimate may cause extra uncertainty to the 1,064 nm backscatter, which is beyond our 5% error estimate for particulate backscatter at all wavelengths. The comparison of the depolarization ratio fit to the measurements was demonstrated in the bottom panels of Figure 7A. In the altitude range from 2.5 to 5.2 km, the fitting at the 532 nm wavelength is noticeable to be beyond the measurement error. Considering that the depolarization is highly relevant to aerosol non-sphericities, such a deviation may be caused by the modeling errors with the pre-assumed aspect ratio distribution of spheroidal particles (Legrand et al., 2014). In Figure 7B, we compare the fit of aerosol lidar ratio ($S_a = \alpha_a/\beta_a$) to the measurements. Plotted as the horizontal error bars, the measurement uncertainty with aerosol lidar ratio ($\Delta S_a$) is derived from the uncertainties of extinction ($\sigma_{\alpha_a} = 17 \text{ Mm}^{-1}$) and backscatter coefficients ($\epsilon_{\beta_a} = 5\%$) via error propagation $\Delta S_a = \alpha_a/\beta_a \sqrt{(\sigma_{\alpha_a}/\alpha_a)^2 + (\epsilon_{\beta_a})^2}$. In the altitude ranges from 0 to 1 km and from 2.5 to 5.2 km where most aerosols are distributed, the fitting errors are observed to be less than 20
<figure><img src="image_8.png" /></figure>

**FIGURE 6** | (A) Retrieved profiles of AOD (**left panel**) and SSA (**right panel**) at 532 nm (black curves), with comparison to the collocated HiGEAR measurement (red dot); (B): retrieved fine mode (**left panel**) and coarse mode (**right panel**) aerosol size distributions (black curves), with comparison to the *in-situ* HiGEAR measurements via UHSAS and APS (red curves). The retrieved AOD, SSA and size distribution are derived from the speciated aerosol abundance and intrinsic quantities displayed in **Figure 5**. Moreover, the retrieved AOD in the left panel of **Figure 6A** is accumulated AOD from TOA to different altitudes with 315 m resolution. The horizontal bars in all panels indicate the retrieval errors. Collocation of HiGEAR and RSP/HSRL-2 measurements was demonstrated in **Figure 4**. To estimate the retrieval uncertainties, the measurement uncertainties are same as those used in **Figure 5**.

Sr, or less than 20% in relative sense viewing the measured lidar ratio is larger than 100 in 355 and 532 nm. Such a difference is within the measurement uncertainties of the lidar ratio indicated by the horizontal error bars.

The retrieved BRF and DOLP are compared to the measured data by RSP in **Figure 8**. The spectral mean of fitting error for BRF and DOLP is less than ~2% and 0.003 respectively in the four bands centered at 410, 470, 670, and 865 nm. By contrast, the fitting error at 550 nm is relatively larger (~4.5% for BRF and ~0.007 for DOLP), which could be caused by the approximation in estimating the ozone absorption to be ~0.03 in optical depth at 550 nm for radiance calculation, and/or due to the residue error of instrument calibration. Moreover, the fitting errors have certain angular dependence: relatively larger fitting errors can be observed when the polarimeter's view angle gets more and more off nadir. This could be caused by the inhomogeneity in the spatial distribution of aerosols in the measurement zone, as seen from the lidar curtain in **Figure 3**.

Following the case study demonstrated through **Figures 5–8**, **Figure 9** gives a systematic comparison of AOD and SSA at 532 nm from all 50 sets of retrievals to the measurements by 4-STAR and HiGEAR instruments. Starting from the maximum flight altitude (6 km) of P-3 aircraft, which carries 4-STAR and

HiGEAR, the accumulated AOD measured by 4-STAR is sampled at a resolution of 315 m. Plotted in red dots, comparison of retrieval to 4-STAR reference data shows the MAD to be around 0.006, which is within the AOD retrieval errors for all cases. The determinant of regression is about 0.99. As another measure of the quality of retrieved resolution in fitting observation, we also compare retrieved AOD to HSRL-2 AOD calculated from its extinction measurements. As indicated with the blue squares, most data fall very close to the 1:1 line - indicating a high quality of the agreement. Also, at the vertical resolution of 315 m, the right panel of **Figure 9** shows the MAD of SSA retrievals is 0.028. Overall, the retrieved SSAs are slightly larger than HiGEAR reference data derived from its measurement of absorption and scattering coefficients. This could be due to the retrieved aerosol properties from RSP and HSRL-2's observations are under ambient conditions, while the *in-situ* measurements of absorption and scattering coefficients are under dry conditions. However, seeking an exact explanation for such a difference is hard since it is within the estimated retrieval errors of SSA. Indeed, gaining an SSA retrieval accuracy of 0.028 is largely due to the constraints provided by RSP's multiangle polarimetry, by HSRL's direct measurement of extinction, as well as by the smoothness constraints imposed on the variation of refractive
<figure><img src="image_9.png" /></figure>

**FIGURE 7** | Comparison of HSRL-2 measurements with model fit using the retrieved aerosol properties displayed in **Figure 6**. Consistent with **Figures 5** and **6**, the observations were acquired at 14:25:10 UTC on September 12, 2016 during ORACLES campaign. **(A)** Upper panels: extinction coefficients at 355 nm (**left**) and 532 nm (**right**) as a function of altitude; Middle panels: backscatter coefficient at 355 nm (**left**), 532 nm (**middle**), and 1,064 nm (**right**) as a function of altitude. Bottom panels: depolarization ratio at 355 nm (**left**), 532 nm (**middle**), and 1,064 nm (**right**) as a function of altitude. **(B)** Lidar ratio at 355 nm (**left**) and 532 nm (**right**). The horizontal bars in all panels indicate the measurement uncertainties. The measurement uncertainties are same as those used in **Figure 5**.

index and size distribution. The benefits of using polarimetry to constrain the retrieval of aerosol SSA has been demonstrated by multiple airborne polarimeters such as RSP (Wu et al., 2015), AirMSPI (Xu et al., 2017), AirSPEX (Hasekamp et al., 2019) and AirHARP (Puthukkudy et al., 2020). Generally, retrieval error is within ~0.02–0.04 when the absolute error with DOLP measurement is less than ~0.005. However, in this work, we are retrieving vertically dependent aerosol properties. Therefore, the parameter space is significantly expanded, and the polarimeter’s capability may spread increasingly thin as the number of vertical layers for retrieval to resolve increases. To mitigate this issue, extra a priori constraints regarding smooth variations of some aerosol properties are imposed, and they play a

significant role in retrieval accuracy enhancement. Similar to the comparison of AOD in the left panel of **Figure 9**, the left and right panels of **Figure 10** give a comparison of the accumulated AODs against the product of 4-STAR at 355 and 1,064 nm, respectively. The MADs of AOD at these two wavelengths are found to be 0.010 and 0.004, respectively. As a retrieval quality check, the comparison was also made to AOD calculated from HSRL-2's extinction measurement in 355 nm. Compared to the agreement of retrieval AOD and HSRL-2 AOD in 532 nm, the agreement of the two AODs in 355 nm slightly reduces. For AOD at 1,064 nm, the difference in the retrieval and 4-STAR AODs falls beyond the range of retrieval error bars at some points. To fill the gap, we may need to further account for the uncertainties in 4-STAR AOD
<figure><img src="image_10.png" /><figcaption>FIGURE 8 | (A) Comparison of RSP radiance measurements (in BRF unit) with model fit using the retrieved aerosol properties displayed in Figure 5. Consistent with Figures 5 and 6, the observations were acquired at 14:25:10UTC on September 12, 2016 during ORACLES campaign. Comparison is made for five wavelengths from near UV (410 nm) to near-infrared (865 nm) in different panels. (B) Comparison of RSP DOLP measurements with model fits using the retrieved aerosol properties displayed in Figure 5. The vertical bars in all panels indicate the measurement uncertainties.</figcaption></figure>

measurements (e.g., ~0.01 for above-cloud AOD sampled at 500 nm and 1,020 nm according to LeBlanc et al., (2020)).

Figures 11 and 12 compare the effective radius and effective variance of fine and coarse mode aerosols against the reference data derived from in-situ measurements by use of the UHSAS and APS instruments aboard P-3 aircraft. We filtered out the results above 5.2 km where aerosols loading is very low so that the retrieved aerosol microphysical properties are subjected to relatively larger retrieval errors. Moreover, the UHSAS/APS measured the size-resolved number concentration are fit with the number weighted log-normal size distribution expressed in Eq. 103. Following the derivation of $r_g^{(s)}$ and $\sigma_g^{(s)}$ from fitting in-situ measurements, the effective radius and variance is calculated via

$$r_{\text{eff}}^{(s)} = r_g^{(s)} e^{\frac{5}{2} [\sigma_g^{(s)}]^2}, \quad (104)$$

$$v_{\text{eff}}^{(s)} = e^{[\sigma_g^{(s)}]^2} - 1. \quad (105)$$

Comparison of fine and coarse mode aerosol effective radii and variance in the left panel of Figures 11 and 12 indicates the retrieved values are larger than their counterparts as measured by UHSAS and APS. The MADs of effective radii are around 0.012 and 0.377 µm for fine and coarse aerosols, respectively. Such a difference is basically within the estimated range of retrieval errors. MADs of effective variance are around 0.132 and 0.257 for fine and coarse aerosols, respectively. Such a deviation is beyond the range of estimated retrieval errors. However, the in-situ aerosol size measurement by HiGEAR suite of instruments has 5% (Uin 2016) and 10% (Pfeifer et al., 2016) uncertainty for fine and coarse mode aerosols, respectively. Moreover, HiGEAR measurement was carried out under dry conditions while the retrieved aerosol properties from RSP and HSRL-2's observations
<figure><img src="image_11.png" /></figure>

**FIGURE 9** | Systematic comparison of accumulated partial-column AOD and layer-effective SSA from retrieval using 50 combined sets of HSRL-2 and RSP observations to the reference measurement by 4-STAR and HiGEAR. Left panel: comparison of accumulated AOD at 532 nm against 4-STAR and HSRL-2 measurements; Right panel: comparison of SSA at 532 nm against HiGEAR reference data. Collocation of these measurements was demonstrated in **Figure 4**. For SSA comparison, less than 50 retrieval data points are plotted due to the screening of data for altitude above 5.2 km (too low aerosol loading). The vertical bars in all panels indicate the retrieval errors.

<figure><img src="image_12.png" /></figure>

**FIGURE 10** | Same as the left panel of **Figure 9** but the comparison against 4-STAR and HSRL-2 reference data for AOD at 355 nm, and against 4-STAR reference data for AOD at 1,064 nm. Left panel: comparison of AOD at 355 nm; Right panel: comparison of AOD at 1,064 nm. The vertical bars in all panels indicate the retrieval errors.

<figure><img src="image_13.png" /></figure>

**FIGURE 11** | Comparison of effective radius (**left panel**) and variance (**right panel**) of fine mode aerosol size distribution. The vertical bars in all panels indicate the retrieval errors.
<figure><img src="image_14.png" /><figcaption>FIGURE 12 | Same as Figure 11 but for coarse mode aerosol size distribution.</figcaption></figure>

were under ambient conditions. These are potential sources that can cause the difference in the estimated geometric number mean radius and standard deviation from HiGEAR measurements, which further impacts the comparison of aerosol effective radius and variance.

## SUMMARY

An optimization approach is developed to combine the use of lidar and polarimetric observations to retrieve the vertical profiles of speciated aerosol abundance and properties over the ocean. These aerosol properties include aerosol refractive index, size distribution, and spherical particle fractions. Acquired during 2016 ORACLES field campaign, the HSRL-2 measurements of extinction, backscatter and depolarization in combination with RSP measurements of radiance and DOLP are used to test the retrieval. Comparison to the accumulated AODs from TOA to different altitudes measured by 4-STAR aboard the P-3 at different altitudes (up to 6 km) shows an MAD of 0.010, 0.006, and 0.004 for 355, 532 and 1,064 nm, respectively. At the vertical resolution of 315 m, validation of SSA against the in-situ measurement by HiGEAR shows an accuracy of 0.028 for 532 nm. Comparing the retrieved effective radii of fine and coarse mode aerosols to the reference data from in-situ measurements by UHSAS and APS, we found the MADs to be around 0.012 and 0.377 µm for fine and coarse mode aerosols, respectively.

While most comparisons indicate that the differences are within the estimated retrieval errors or within measurement uncertainties of relevant quantities, it is noteworthy that we assume 1) a plane-parallel atmosphere-surface system for computing polarized radiative transfer and 2) sphere and spheroid model for computing aerosol single-scattering properties. These assumptions are subjected to errors. One error source that is particularly relevant to lidar-polarimeter combined retrieval under the context of high spatial resolution is the difference in the atmosphere viewed along nadir direction by lidar and the atmosphere viewed along oblique angles of a polarimeter. In a strict sense, combined retrieval should

only use nadir polarimetric measurements, or all angular measurements but with a 3D forward RT model that accurately accounts for the spatial inhomogeneity with the atmosphere. However, the former option will reduce the benefits of multi-angular benefits in constraining while the latter option will reduce the modeling efficiency. As a trade-off between retrieval accuracy and efficiency, we use all angular measurements and retain the simplicity of assuming a plane-parallel atmosphere-surface system to perform the retrieval. Then the modeling errors are factored into the error estimates (see Retrieval Error Estimate). Nevertheless, seeking an effective way to further mitigate modeling errors while retaining the modeling efficiency is a worthwhile topic to pursue in the next stage.

## DATA AVAILABILITY STATEMENT

Publicly available datasets were analyzed in this study. This data can be found here: https://espoarchive.nasa.gov/archive/browse/oracles.

## AUTHOR CONTRIBUTIONS

FX formulated the inversion approach and developed the prototype codes for retrieval use. LG interfaced the retrieval with OU Supercomputing Center for Education and Research, implemented the retrieval test, and provided inputs on interpreting the results. JR, CF, WE, AS, SS, SB, XL, RF, BC and OD provided insightful discussions on the algorithm development and proper use of HSRL-2 and RSP data for retrieval test, and made editorial changes.

## ACKNOWLEDGMENTS

The authors acknowledge support from National Aeronautics and Space Administration (NASA) for conducting the study presented in this paper. O. Dubovik acknowledges the support of the Labex CaPPA project, which is funded by the French National Research Agency under contract ANR-11-LABX-
0005-01. F. Xu thanks Dr. Yongxiang Hu and Dr. Xiaomei Lu from NASA Langley Research Center for providing helpful discussions on lidar remote sensing algorithm development. The authors are grateful to the reviewers for their constructive

comments on the early version of this manuscript. A majority of the computation for this project was performed at University of Oklahoma's Supercomputing Center for Education and Research (OSCR).

## REFERENCES

Burton, S. P., Hair, J. W., Kahnert, M., Ferrare, R. A., Hostetler, C. A., Cook, A. L., et al. (2015). Observations of the spectral dependence of linear particle depolarization ratio of aerosols using NASA Langley airborne high spectral resolution lidar. *Atmos. Chem. Phys.* 15, 13453-13473. doi:10.5194/acp-15-13453-2015

Burton, S. P., Hostetler, C. A., Cook, A. L., Hair, J. W., Seaman, S. T., Scola, S., et al. (2018). Calibration of a high spectral resolution lidar using a Michelson interferometer, with data examples from ORACLES. *Appl. Optic.* 57, 6061-6075. doi:10.1364/AO.57.006061

Burton, S. P., Stamnes, S., Liu, X., Dawson, K., Chemyakin, E., Ferrare, R., et al. (2019). High spectral resolution lidar for aerosol characterization and combined lidar + polarimeter retrieval. Available at: https://www.giss.nasa.gov/staff/mmishchenko/APOLO-2/Abstracts/Burton.pdf.

Cairns, B., Russell, E. E., and Travis, L. D. (1999). The research scanning polarimeter: calibration and ground-based measurements. *PROC SPIE*. 3754, 186. doi:10.1117/12.366329

Chaikovsky, A., Dubovik, O., Holben, B., Bril, A., Goloub, P., Tanré, D., et al. (2016). Lidar-Radiometer Inversion Code (LIRIC) for the retrieval of vertical aerosol properties from combined lidar/radiometer data: development and distribution in EARLINET. *Atmos. Meas. Tech.* 9, 1181-1205. doi:10.5194/amt-9-1181-2016

Chowdhary, J., Cairns, B., Waquet, F., Knobelspiesse, K., Ottaviani, M., Redemann, J., et al. (2012). Sensitivity of multiangle, multispectral polarimetric remote sensing over open oceans to water-leaving radiance: analyses of RSP data acquired during the MILAGRO campaign. *Remote Sens. Environ.* 118, 284-308. doi:10.1016/j.rse.2011.11.003

Cox, C., and Munk, W. (1954a). Measurement of the roughness of the sea surface from photographs of the sun's glitter. *J. Opt. Soc. Am.* 44, 838-850.

Cox, C., and Munk, W. (1954b). Statistics of sea surface derived from sun glitter. *J. Mar. Res.* 13, 198-227.

de Haan, J. F., Bosma, P. B., and Hovenier, J. W. (1987). The adding method for multiple scattering calculations of polarized light. *Astron. Astrophys.* 181, 371-391.

Dubovik, O., Herman, M., Holdak, A., Lapyonok, T., Tanré, D., Deuzé, J. L., et al. (2011). Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations. *Atmos. Meas. Tech.* 4, 975-1018. doi:10.5194/amt-4-975-2011

Dubovik, O., and King, M. D. (2000). A flexible inversion algorithm for retrieval of aerosol optical properties from Sun and sky radiance measurements. *J. Geophys. Res.* 105, 673-696. doi:10.1029/2000JD0900282

Dubovik, O., Lapyonok, T., Litvinov, P., Herman, M., Fuertes, D., Ducos, F., et al. (2014). GRASP: a versatile algorithm for characterizing the atmosphere. SPIE Newsroom 2014.

Dubovik, O., Li, Z., Mishchenko, M. I., Tanré, D., Karol, Y., Bojkov, B., et al. (2019). Polarimetric remote sensing of atmospheric aerosols: instruments, methodologies, results, and perspectives. *J. Quant. Spectrosc. Radiat. Transfer.* 224, 474-511. doi:10.1016/j.qsrst.2018.11.024

Dubovik, O. (2004). "Optimization of numerical inversion in photopolarimetric remote sensing," in *Photopolarimetry in remote sensing*. Editors G. Videen, Y. Yatskiv, and M. I. Mishchenko (Dordrecht, The Netherlands: Kluwer Academic Publishers), 65-106.

Dubovik, O., Sinyuk, A., Lapyonok, T., Holben, B. N., Mishchenko, M. I., Yang, P., et al. (2006). Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust. *J. Geophys. Res.* 111, D11208. doi:10.1029/2005JD006619

Dunagan, S. E., Johnson, R., Zavaleta, J., Russell, P. B., Schmid, B., Flynn, C., et al. (2013). Spectrometer for sky-scanning sun-tracking atmospheric research (4STAR): instrument technology. *Rem. Sens.* 5, 3872-3895. doi:10.3390/rs5083872

Espinosa, W. R., Castellanos, P., Sayer, A., Colarco, P., Kemppien, O., Nowotnick, E., et al. (2019). An exploration of joint LIDAR and multiangle polarimeter aerosol retrieval capabilities using the GRASP algorithm and OSSE data derived from the GEOS model. Available at: https://www.giss.nasa.gov/staff/mmishchenko/APOLO-2/Abstracts/Espinosa.pdf.

Esposito, L. W. (1979). An "adding" algorithm for the Markov chain formalism for radiation transfer. *Astrophys. J.* 233, 661-663.

Esposito, L. W., and House, L. L. (1978). Radiative transfer calculated from a Markov chain formalism. *Astrophys. J.* 219, 1058-1067.

Evans, K. F., and Stephens, G. L. (1991). A new polarized atmospheric radiative transfer model. *J. Quant. Spectrosc. Radiat. Transfer.* 46, 413-423. doi:10.1016/0022-4073(91)90043-P

Froyd, K. D., Murphy, D. M., Brock, C. A., Campuzano-Jost, P., Dibb, J. E., Jimenez, J.-L., et al. (2019). A new method to quantify mineral dust and other aerosol species from aircraft platforms using single-particle mass spectrometry. *Atmos. Meas. Tech.* 12, 6209-6239. doi:10.5194/amt-12-6209-2019

Gao, M., Zhai, P-W., Franz, B., Hu, Y., Knobelspiesse, K., Werdell, P. J., et al. (2018). Retrieval of aerosol properties and water leaving reflectance from multiangular polarimetric measurements over coastal waters. *Optic Express.* 26, 8968-8989. doi:10.1364/OE.26.008968

GCOS (2011). *Systematic observation requirements for satellite-based data products for climate; 2011 update, supplemental details to the satellite-based component of the "implementation plan for the global observing system for climate in support of the UNFCCC (2010 update)*. Geneva, Switzerland: WMO.

Grainger, R. G. (2017). Some useful formulae for aerosol size distributions and optical properties. Available at: http://eodg.alm.ox.ac.uk/user/grainger/research/aerosols.pdf.

Hair, J. W., Caldwell, L. M., Krueger, D. A., and She, C.-Y. (2001). High-spectral-resolution lidar with iodine-vapor filters: measurement of atmospheric-state and aerosol profiles. *Appl. Optic.* 40, 5280-5294. doi:10.1364/AO.40.005280

Hair, J. W., Hostetler, C. A., Cook, A. L., Harper, D. B., Ferrare, R. A., Mack, T. L., et al. (2008). Airborne high spectral resolution Lidar for profiling aerosol optical properties. *Appl. Optic.* 47, 6734-6752. doi:10.1364/AO.47.006734

Hansen, J. E. (1971). Multiple scattering of polarized light in planetary atmospheres. Part II. Sunlight reflected by terrestrial water clouds. *J. Atmos. Sci.* 28, 1400-1426. doi:10.1175/1520-0469(1971)028<1400:MSOPLI>2.0.CO;2

Hasekamp, O., Litvinov, P., and Butz, A. (2011). Aerosol properties over the ocean from PARASOL multi-angle photopolarimetric measurements. *J. Geophys. Res.* 116, D14204. doi:10.1029/2010JD015469

Hasekamp, O. P., Fu, G., Rusli, S. P., Wu, L., Noia, A. D., aan de Brugh, J., et al. (2019). Aerosol measurements by SPEXone on the NASA PACE mission: expected retrieval capabilities. *J. Quant. Spectrosc. Radiat. Transfer.* 227, 170-184. doi:10.1016/j.qsrst.2019.02.006

Knobelspiesse, K., Cairns, B., Mishchenko, M., Chowdhary, J., Tsigaridis, K., van Diedenhoven, B., et al. (2012). Analysis of fine-mode aerosol retrieval capabilities by different passive remote sensing instrument designs. *Optic Express.* 20, 21457-21484. doi:10.1364/OE.20.021457

Knobelspiesse, K., Cairns, B., Ottaviani, M., Ferrare, R., Hair, J., Hostetler, C., et al. (2011). Combined retrievals of boreal forest fire aerosol properties with a polarimeter and lidar. *Atmos. Chem. Phys.* 11, 7045-7067. doi:10.5194/acp-11-7045-2011

Koeppke, P. (1984). Effective reflectance of oceanic whitecaps. *Appl. Optic.* 23, 1816-1824.

LeBlanc, S. E., Redemann, J., Flynn, C., Pistone, K., Kacenelenbogen, M., Segal-Rosenheimer, M., et al. (2020). Above-cloud aerosol optical depth from airborne observations in the southeast Atlantic. *Atmos. Chem. Phys.* 20, 1565-1590. doi:10.5194/acp-20-1565-2020

Lee, Z. P., Hu, C., Shang, S. L., Du, K. P., Lewis, M., Arnone, R., et al. (2013). Penetration of UV-Visible solar light in the global oceans: insights from ocean color remote sensing. *J. Geophys. Res.* 118, 4241-4255. doi:10.1002/jgrc.20308
Legrand, M., Dubovik, O., Lapyonok, T., and Derimian, Y. (2014). Accounting for particle non-sphericity in modeling mineral dust radiative properties in the thermal infrared. *J. Quant. Spectrosc. Radiat. Transfer.* 149, 219–240. doi:10.1016/j.jqsrt.2014.07.014

Liu, X, Stamnes, S., Ferrare, R. A., Hostetler, C., Burton, S., Chemyakin, E., et al. (2016). Development of a combined lidar-polarimeter inversion algorithm for retrieving aerosol. Abstract A54E-06, AGU.

Lopatin, A., Dubovik, O., Chaikovsky, A., Goloub, P., Lapyonok, T., Tanré, D., et al. (2013). Enhancement of aerosol characterization using synergy of lidar and sun-photometer coincident observations: the GARRLiC algorithm. *Atmos. Meas. Tech.* 6, 2065–2088. doi:10.5194/amt-6-2065-2013

McNaughton, C. S., Clarke, A. D., Kapustin, V., Shinozuka, Y., Howell, S. G., Anderson, B. E., et al. (2009). Observations of heterogeneous reactions between Asian pollution and mineral dust over the Eastern North Pacific during INTEX-B. *Atmos. Chem. Phys.* 9, 8283–8308. doi:10.5194/acp-9-8283-2009

Mishchenko, M. I., Yatskiv, Y. S., Rosenbush, V. K., and Videen, G. (2011). *Polarimetric detection, characterization, and remote sensing*. Cham, Switzerland: Springer.

Pfeifer, S., Müller, T., Weinhold, K., Zikova, N., Martins dos Santos, S., Marinoni, A., et al. (2016). Intercomparison of 15 aerodynamic particle size spectrometers (APS 3321): uncertainties in particle sizing and number size distribution. *Atmos. Meas. Tech.* 9, 1545–1551. doi:10.5194/amt-9-1545-2016

Puthukkudy, A., Martins, J. V., Remer, L. A., Xu, X., Dubovik, O., Litvinov, P., et al. (2020). Retrieval of aerosol properties from airborne hyper-angular rainbow polarimeter (AirHARP) observations during ACEPOL 2017. *Atmos. Meas. Tech.* 13, 5207–5236. doi:10.5194/amt-13-5207-2020

Redemann, J., Wood, R., Zuidema, P., Doherty, S. J., Luna, B., LeBlanc, S. E., et al. (2021). An overview of the ORACLES (observations of aerosols above clouds and their interactions) project: aerosol-cloud-radiation interactions in the Southeast Atlantic basin. *Atmos. Chem. Phys.* 21, 1507–1563. doi:10.5194/acp-2020-449

Reid, J. S., Jonsson, H. H., Maring, H. B., Smirnov, A., Savoie, D. L., Cliff, S. S., et al. (2003). Comparison of size and morphological measurements of coarse mode dust particles from Africa. *J. Geophys. Res.* 108 (D19), 8593. doi:10.1029/2002JD002485

Remer, L. A., Knobelspiesse, K. D., Zhai, P.-W., Xu, F., Kalashnikova, O., Chowdhary, J., et al. (2019). Retrieving aerosol characteristics from the PACE mission, Part 2: multi-angle and polarimetry. *Front. Environ. Sci.* 7, 94. doi:10.3389/fenvs.2019.00094

Rogers, R. R., Hair, J. W., Hostetler, C. A., Ferrare, R. A., Obland, M. D., Cook, A. L., et al. (2009). NASA LaRC airborne high spectral resolution lidar aerosol measurements during MILAGRO: observations and validation. *Atmos. Chem. Phys.* 9, 4811–4826. doi:10.5194/acp-9-4811-2009

Shipley, S. T., Tracy, D. H., Eloranta, E. W., Trauger, J. T., Sroga, J. T., Roesler, F. L., et al. (1983). High spectral resolution lidar to measure optical-scattering properties of atmospheric aerosols. I. theory and Instrumentation. *Appl. Optic.* 22, 3716–3724. doi:10.1364/AO.22.003716

Stamnes, S., Hostetler, C., Ferrare, R., Burton, S., Liu, X., Hair, J., et al. (2018). Simultaneous polarimeter retrievals of microphysical aerosol and ocean color parameters from the “MAPP” algorithm with comparison to high-spectral-resolution lidar aerosol and ocean products. *Appl. Optic.* 57, 2394–2413. doi:10.1364/AO.57.002394

Stokes, G. G. (1862). On the intensity of the light reflected from or transmitted through a pile of plates. *Proc. Roy. Soc. Lond.* 11, 545–556.

IPCC (2013). *Climate change 2013: the physical science basis. Contribution of working group I to the fifth assessment report of the intergovernmental panel on climate change*. New York, NY: Cambridge.

Uin, J. (2016). Ultra-high sensitivity aerosol spectrometer (UHSAS) instrument handbook. ARM Tech. Rep. DOE/SC-ARM-TR-163, 17.

van de Hulst, H. C. (1963). A new look at multiple scattering. Goddard Institute for Space Studies. Tech. Rep.NASA TM-I03044.

van de Hulst, H. C. (1981). *Light scattering by small particles*. New York, NY: Dover.

Vaughan, M. (2004). Algorithm for retrieving lidar ratios at 1064 nm fromspace-based lidar backscatter data. *Proc. SPIE.* 5240, 104–115. doi:10.1117/12.510770

Weitkamp, C. (2005). *LIDAR: range-resolved optical remote sensing of the atmosphere*. New York, NY: Springer.

Winker, D. M., Couch, R. H., and McCormick, P. (1996). An overview of LITE: NASA’s lidar in-space technology experiment. *Proc. IEEE.* 84, 164–180.

Winker, D. M., Hunt, W. H., and McGill, M. J. (2007). Initial performance assessment of CALIOP. *Geophys. Res. Lett.* 34, L19803. doi:10.1029/2007GL030135

Wu, L., Hasekamp, O., van Diedenhoven, B., and Cairns, B. (2015). Aerosol retrieval from multiangle, multispectral photopolarimetric measurements: importance of spectral range and angular resolution. *Atmos. Meas. Tech.* 8, 2625–2638. doi:10.5194/amt-8-2625-2015

Xu, F., Davis, A. B., and Diner, D. J. (2016). Markov chain formalism for generalized radiative transfer in a plane-parallel medium, accounting for polarization. *J. Quant. Spectrosc. Radiat. Transfer.* 184, 14–26. doi:10.1016/j.jqsrt.2016.06.004

Xu, F., Davis, A. B., Suniti, S. V., Martonchik, J. V., and Diner, D. J. (2012). Linearization of Markov chain formalism for vector radiative transfer in a plane-parallel atmosphere/surface system. *Appl. Optic.* 51, 3491–3507. doi:10.1364/AO.51.0003491

Xu, F., Davis, A. B., West, R. A., and Esposito, L. W. (2010). Markov chain formalism for polarized light transfer in plane-parallel atmospheres, with numerical comparison to the Monte Carlo method. *Optic Express.* 19, 946–967. doi:10.1364/OE.19.000946

Xu, F., Davis, A. B., West, R. A., Martonchik, J. V., and Diner, D. J. (2011). Markov chain formalism for vector radiative transfer in a plane-parallel atmosphere overlying a polarizing surface. *Opt. Lett.* 36, 2083–2085. doi:10.1364/OL.36.002083

Xu, F., Diner, D. J., Dubovik, O., and Schechner, Y. (2019). A correlated multi-pixel inversion approach for aerosol remote sensing. *Rem. Sens.* 11, 746. doi:10.3390/rs11070746

Xu, F., van Harten, G., Diner, D. J., Kalashnikova, O. V., Seidel, F. C., Bruegge, C. J., et al. (2017). Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager. *J. Geophys. Res. Atmos.* 122, 7004–7026. doi:10.1002/2017JD026776

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration, Dr. Xu, Dr. Gao, Dr. Redemann, Dr. Flynn, and Dr. Dubovik.

At least a portion of this work is authored by W. Reed Espinosa, Arlindo M. da Silva, Snorre Stamnes, Sharon P. Burton, Xu Liu, Richard Ferrare and Brian Cairns on behalf of the U.S. Government and, as regards Dr. Espinosa, Dr. da Silva, Dr. Stamnes, Dr. Burton, Dr. Liu, Dr. Farrare, Dr. Cairns, U.S. copyright protection does not attach to separable portions of a Work authored solely by U.S. Government employees as part of their official duties. The U.S. Government is the owner of foreign copyrights in such separable portions of the Work and is a joint owner (with any non-U.S. Government author) of U.S. and foreign copyrights that may be asserted in inseparable portions the Work. The U.S. Government retains the right to use, reproduce, distribute, create derivative works, perform and display portions of the Work authored solely or co-authored by a U.S. Government employee. Non-U.S. copyrights also apply. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.