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

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.


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 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τ total/fine + 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 −1 , 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-ofatmosphere (TOA) in multi-spectral/hyper-spectral and multiangular 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 multiwavelength 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 onthe-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 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 polarizartion (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 (s) l (r)/d ln r] for all (N s ) aerosol species, we have the subcolumn volumetric concentrations in lth layer, where the size distribution of the sth species dv (s) l (r)/d ln r is parameterized by a geometric volume mean radius r (s) m and a logarithmic standard deviation σ (s) m,l : In Eq. 1, C (s) v,l is the subcolumn volume concentration of the sth aerosol species in the lth layer and is related to its volume fraction (f (s) v,l ) and the total column concentration If we discretize the speciated aerosol size distribution in the lth layer into N bin size bins, the volumetric fraction of the ith bin is where r i stands for the central radius of ith size bin and dv (s) l (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 where C (s),sphere n 0 ,a,l is the aerosol number concentration (m −3 ), dn (s) l (r)/d ln r is the number weighted aerosol size distribution in association with dv (s) l (r)/d ln r, r min and r max are the lower and upper bounds of the size distribution respectively, c (s),sphere ext,a,l is extinction cross section as a function of particle refractive index m (s) a,l and radius r, f (s),sphere v,l is the volume fraction of spherical particles, v(r) is the volume of a spherical particle of radius r, and C (s) v0,a,l is the aerosol volumetric concentration (in the unit of m 3 /m 3 ). It is related to the subcolumn aerosol concentration (C (s) v,a,l in the unit of m 3 /m 2 ) in lth layer by where 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, where r becomes the volume equivalent spherical particle radius, the aspect ratio η < 1 for oblate spheroids and η > 1 for prolate spheroids, and η min and η 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 (s) l (η)/d ln η is independent of aerosol size. Both dn (s) l (η)/d ln η and dv (s) l (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 η, respectively.
For numerical calculation, the aerosol size distribution is descritized into N bin rectangular bins so that for spherical particles where the component of sth aerosol species in its ith size bin is For spheroids, we further need to break the integral of volumetric extinction into N AR bins of aspect ratio so that: For simplicity, the layer resolved volumetric extinction coefficient can be rewritten in the following form, ext,a,l,i where K (s),sphere ext,a,l,i and K (s),spheroid ext,a,l,i are given in Eqs 10 and 12, respectively.
Replacing the subscript "ext" by "sca" in above derivations  gives relevant equations for the volumetric scattering coefficient: sca,a,l,i In a similar way, the phase matrix element P a,j,k,l of sth aerosol species in the lth layer is calculated from In the numerator of the above equation, we have Summing over the contributions from all aerosol species gives the total volumetric extinction and scattering coefficients in lth layer, where K (s) ext, a,l and K (s) sca,a,l 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, Contributed by all aerosol species, the overall SSA contributed by all aerosol species in lth layer is where f (s) v,l 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 Frontiers in Remote Sensing | www.frontiersin.org April 2021 | Volume 2 | Article 620871 5 where P (s) a, j,k,l is given by Eq. 15. The layer aerosol optical depth (AOD, Δτ a,l ) is contributed by all aerosol species. Therefore, we have Δτ a,l Ns s 1 In the above equation, the AOD for sth aerosol species is where C (s) v,a,l is the subcolumn volumetric concentration of sth aerosol species, evaluated by Eq. 6, and K (s) ext,a,l is column effective extinction coefficient (in the unit of m −1 ) 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 (α a,l , in the unit of m −1 ) is calculated by summing over the contribution from all aerosol species, namely where the speciated aerosol extinction α (s) a,l in a layer of geometric thickness ΔH l is related to layer resolved aerosol optical depth by where C (s) v,a,l and K (s) ext,a,l are given by Eqs 6 and 7 respectively. Calculation of the total backscatter β a,l (in the unit of m −1 Sr −1 ) also needs to account for the contributions of all aerosol species as well, namely where the speciated aerosol backscatter is evaluated by where P (s) a,1,1,l (180 + ) is given by Eq. 15 and the aerosol scattering coefficient (β (s) sca,a,l , in the unit of m −1 ) in lth layer is evaluated by, where K (s) sca,a,l is given by Eq. 16. Assuming random orientation of nonspherical particles, the aerosol linear depolarization ratio of lth layer is calculated as δ a,l (180 + ) 1 − P a,2,2,l (180 + ) P a,1,1,l (180 + ) 1 + P a,2,2,l (180 + ) P a,1,1,l (180 + ) , 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 (α a , β a , δ 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), where m r0 and m i0 are the real and imaginary parts of the refractive index at the reference wavelength λ 0 , and κ r0 and κ 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 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 −8 to 3 × 10 −4 with five grids and from 5 × 10 −4 to 0.5 with 15 grids. The lower and upper bounds for aerosol size are set as r min 5 nm and r 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 Xu et al., 2016;Gao et al., 2018; Frontiers in Remote Sensing | www.frontiersin.org April 2021 | Volume 2 | Article 620871 6 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 surf is described by (Xu et al., 2016), where D 0 is a zero matrix except D 0,11 1; a Foam is foam albedo; f Foam is foam coverage fraction related to wind speed V wind by f Foam 2.95 × 10 −6 ×V wind 3.52 (Koepke, 1984); R W,λ is the polarized specular reflection from water surface (Cox and Munk, 1954a;Cox and Munk, 1954b), and R Bio WL,λ 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 WL,λ " is the Lambertian correction term. Note that R Bio WL,λ 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 WL to be 0 the simplified and the empirically adjusted biooptical 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 Bio WL and Δa WL , where d 0 is the Earth-Sun distance at which F 0 is reported, and d is the Earth-Sun distance at the time of measurement. Note that nLw, R W , R surf , R Bio WL , a Foam , Δa WL , and F 0 in Eqs 31 and 32 are all spectrally dependent.

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 Layer layers, each bounded by the altitudes h l and h l+1 (h l < h l+1 ) and N Layer 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 layers are calculated from speciated 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 (Δτ l ), SSA (ω 0,l ) and phase matrix (P l ) of lth layer contributed by both aerosols and Rayleigh-scattering molecules are evaluated by where as a function of scattering angle Θ, P R and P a,l are the Rayleigh and aerosol phase matrices, respectively, the aerosol SSA of lth layer (ω 0,a,l ) is given by Eq. 20 and the AOD of lth layer (Δτ a,l ) 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; FIGURE 2 | The coupled atmosphere-ocean system for RT modeling. The Sun illuminates the top of atmosphere with zenith angle θ 0 and azimuthal angle ϕ 0 (ϕ 0 is assumed to be zero here). The sensor views the atmosphere at viewing angle θ v and azimuthal angle ϕ v . 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.
Frontiers in Remote Sensing | www.frontiersin.org April 2021 | Volume 2 | Article 620871 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]): 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, β a (z) and β m (z) are backscatter coefficients of aerosols and air molecules at range z, and α a and α 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 "[β a (z) + β m (z)]" gives attenuated backscatter. The molecular scattering and extinction coefficients β m (r) and α 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], 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: Due to the anisotropy of the air molecules, the depolarization ratio of the pure molecular atmosphere is nonzero but can be precalculated from climatology. Moreover, the backscatter of a linearly polarized laser beam from spherical particles is totally linearly polarized so that (δ 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 (α a ), backscatter (β a ) and depolarization (δ 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.

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 x that includes retrieval parameters for N s aerosol species can be expressed as, where N is the total number of aerosol and ocean properties to be retrieved, the subvector x (s) contains aerosol properties for sth aerosol species. For a total of N aer 8 aerosol parameters per species, we have, where for sth aerosol species C v contains layer specific column aerosol concentration (m 3 /m 2 ), (m r 0 , κ r 0 ) and (m i 0 , κ i 0 ) are layer resolved model parameters that constrain spectrally dependent real and imaginary parts of aerosol refractive index, respectively (via Eq. 30); r m and σ m contain layer resolved geometric volume median radius and standard deviation of volumetric log-normal distribution; f v sphere contain the layer resolved volume fraction of spherical aerosols. As the last component of the state vector in Eq. 39, x (ocean) 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 -a UHSAS, Ultra-High Sensitivity Aerosol spectrometer; σ UHSAS is derived from fitting in-situ size measurement by UHSAS using the log-normal size distribution (Eq. 103). c APS, Aerodynamic Particle Sizer (APS-TSI3320); σ APS is derived from fitting in-situ size measurement by APS using the log-normal size distribution (Eq. 103). b σ UHSSA and σ APS are adopted when measured lidar extinction is larger than its uncertainty; otherwise the initial guess of σ m is set to 0.2 and its range is relaxed to [0.05, 0.60].
objective is to solve the following system of equations (Dubovik, 2004): where f p i denotes the i th type of constraint, Δf p i is the error with this type of constraint, and x x state as defined in Eq. 39. Formally, the statistical independence of different sources of constraints means that the covariance matrix of joint constraint fp [f p 1 ; f p 2 ; /; f p M ] has the following structure where C i indicates the covariance matrix of the ith constraint (f p i ). Following the expressions of f p i and C f p i , the probability distribution function (PDF) of joint data (1 ≤ i ≤ M) can be derived by multiplying PDFs of data from all M sources of constraints, namely, Further introducing the weight matrix (W), the objective cost function to be minimized has a quadratic form, namely, where In the above equations, ε 2 i is the first diagonal element of C i (i.e., ε 2 i {C i } 11 ) and the Lagrange factor c i weights the contribution of each type of constraint with respective to the first one (c 1 1).
Minimization of Ψ total (x) in Eq. 44 means its gradient with respect to the solution x approaches zero, such that which can be ensured by enforcing the gradient of all components approach zero, namely where 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 x by x − Δx in Eq. 49 and substitute into it. This results in or equivalently, The Jacobian matrix K i in Eqs 49-52 consists of the derivative of the mth observational or a priori data with respect to the nth unknown, More explicit evaluation of f i (x), f p i and K and 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 2 + U 2 ) 1/2 /L. Lidar observations include extinction (α), backscatter or attenuated backscatter (β) and depolarization ratio (δ). Arranging all types of observations into a single column vector we have where the total number of signals N f N L + N DOLP + N α + N β + N δ , where N x 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

introduced in Forward Model for Light Propagation in a Coupled
Atmosphere-Ocean System. The Jacobian matrix K 1 consists of first order partial derivatives with respect to all elements in the state vector x, namely, where K L , K DOLP , K α , K β , and K δ are the Jacobian matrices containing derivatives of L, DOLP, α, β and δ with respect to aerosol parameters x j (1 ≤ j ≤ N), respectively. Accounting for all N parameters in the state vector (see Eq. 39), K L can be expanded into the following matrix form: Frontiers in Remote Sensing | www.frontiersin.org April 2021 | Volume 2 | Article 620871 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 × 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 × 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 where N Layer is the number of layers in which lidar signals are resolved. Assuming extinction α is provided at M α wavelengths, the layer resolved Jacobian matrix (K α,l ) for extinction coefficient equals to, The Jacobian matrix for backscatter coefficient and attenuated backscatter coefficients is evaluated by, and respectively. Assuming β is provided at M β wavelengths, the layer resolved Jacobian matrix for backscatter coefficient (K β,l ) can be expanded into the following specific form, where β β l,a if we use the derived product of particulate backscatter for the lth layer from lidar measurement, and β β l,att if we used direct measured attenuated backscatter signals, where β 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 (δ) equals to Assuming extinction δ is provided at M δ wavelengths, the layer resolved Jacobian matrix for depolarization ratio (K δ,l ) equals to where δ l is contributed by both molecules and aerosols in lth layer and is equal to particulate depolarization ratio (δ δ l,a ) if the derived aerosol depolarization ratio is used directly. The overall covariance matrix incorporates the submatrices from all types of signals, namely, If absolute measurement uncertainty is used and correlation among signals are neglected, then W x is a diagonal matrix. Taking jth signal (x j ) with relative uncertainty σ j as an example, we have W x, jj σ 2 j . If relative uncertainty (ε j ) is used, then we have W x, jj ε 2 j x 2 j . 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, jj ε 2 j .

A Priori Constraints (i 2)
Starting from Eq. 41, the a priori constraint can be explicitly expressed as (Dubovik, 2004), Moreover, K i 2 I (identity matrix) and W i 2 (1/ε 2 a p 11 )C ap in Eqs 49-52. In explicit form, W i 2 is constructed from estimated range of each parameter relative to the first one (Dubovik and King, 2000), 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 where S i 3,m is the differentiation matrix of mth order for third type of constraint and the delta term Δ p 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, For a discretized grid of the variable z, the explicit form of S i,m z, is, Taking the orders of difference m 1 and 2 as examples, we have Application of the above equation to L discretized grids r j (namely, 1 ≤ j ≤ L) leads to so that by invoking Eq. 53, where the matrix S i,m is evaluated by, and 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 f p i 0p into Eq. 45 gives, where the weighting matrix W has the following diagonal terms, and Δ 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 f p i 0 p into Eq. 52 gives Further defining we have

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, where Δ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), As the explicit form, we have where "cΩ" incorporates vertical and spectral smoothness constraints, namely, where the smoothness matrix c h Ω h acts on the vertical variations of aerosol parameters; and the multipliers γ h control the strength of the constraints. The explicit forms of c h Ω h are given by.
where the diagonal terms correspond to the constraints imposed on the N 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 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; ε 2 f 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, ε 2 c , namely,

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, where ε 2 i are the first diagonal elements of the covariance matrices corresponding to ith type of constraints. To estimate ε 2 i for a given parameter to be retrieved (f x) which is a function of z, the most unsmooth known solution f us (z) over the target domain is used as suggested by Dubovik and King (2000) namely, where z min and z max specify the lower and upper bound of z. For c 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 : For retrieving a combined set of RSP and HSRL-2 observations, first estimates of Lagrange parameters c 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 (Δy obs rand ) and systematic errors (Δy syst ) are two major error types. The former is contributed by measurements while the latter is contributed by both measurement (Δy obs syst ) and forward model (Δy 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, where Δx syst is approximated by and covariance matrix of the retrieval solution contributed by the random error is evaluated by where ε rand is the random noise, A is computed at the retrieved solution x via Eq. 82 and ∇Ψ is given by where the total systematic error is expressed as Δy syst 2 Δy obs syst 2 + Δy model 2 .
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 undercharacterization of multiple sources of modeling errors. For practicality, Δx syst is calculated by assuming so that Eq. 93 becomes As the retrieved solution is the closest estimate of the "true" solution x true , we assume x true x retrieved . When a priori x a priori is unavailable, we further assume x a priori x true . Under these assumptions, the last term on the right-hand-side of Eq. 95 disappears (namely c a W −1 a (x true − x a priori ) 0). Note that Eq. 94 treats measurement errors as random. Indeed, these errors are already contained in the observation 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(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, where 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 where C Δ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 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°t o 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 . 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 −1 , 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 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. ; 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 −1 , 5%, and 10% for HSRL-2's particulate extinction, backscatter and depolarization, respectively.
Frontiers in Remote Sensing | www.frontiersin.org April 2021 | Volume 2 | Article 620871 16 continent, the volumetric concentration plot shows the coexistence 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 (s) m and σ (s) m are converted to the geometric number mean radius r (s) g and the standard deviation for sth aerosol species σ (s) g by use of the following relationship (Grainger, 2017): On such a basis, the retrieved number weighted size distribution is calculated as 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 ×(χ/ρ p ) 0.5 , where ρ p is aerosol density and χ is the dynamic shape parameter. Assuming all aerosols measured by APS are sea salt particles, we adopt χ 1.08 and ρ p 1.8 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 χ and assume equivalent spherical diameter for "d" which may cause bias to the estimate of χ (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] μm, respectively. Apparently, the two instruments have an overlapped measurement range of geometric radius from 0.2 to 0.5 μm. Nevertheless, this range is accounted in determining r (s) g and σ (s) g for both fine and coarse aerosols. With these preconditions and conversions, Figure 6B compared the retrieved particle size distribution to HiGEAR's insitu 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 α a /β a ) to the measurements. Plotted as the horizontal error bars, the measurement uncertainty with aerosol lidar ratio (ΔS a ) is derived from the uncertainties of extinction (σ αa 17 Mm −1 ) and backscatter coefficients (ε β a 5%) via error propagation ΔS a α a /β a (σ α a /α a ) 2 + (ε β 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 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 , 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.
Frontiers in Remote Sensing | www.frontiersin.org April 2021 | Volume 2 | Article 620871 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 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 (s) g and σ (s) g from fitting insitu measurements, the effective radius and variance is calculated via r (s) eff r (s) g e 5 2 σ (s) 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 insitu 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 Frontiers in Remote Sensing | www.frontiersin.org April 2021 | Volume 2 | Article 620871 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 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 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.