Joint retrieval of aerosol and water-leaving radiance from multi-spectral , multi-1 angular and polarimetric measurements over ocean 2

Abstract. An optimization approach has been developed for simultaneous retrieval of aerosol properties and normalized water-leaving radiance (nLw) from multispectral, multiangular, and polarimetric observations over ocean. The main features of the method are (1) use of a simplified bio-optical model to estimate nLw, followed by an empirical refinement within a specified range to improve its accuracy; (2) improved algorithm convergence and stability by applying constraints on the spatial smoothness of aerosol loading and Chlorophyll a (Chl a) concentration across neighboring image patches and spectral constraints on aerosol optical properties and nLw across relevant bands; and (3) enhanced Jacobian calculation by modeling and storing the radiative transfer (RT) in aerosol/Rayleigh mixed layer, pure Rayleigh-scattering layers, and ocean medium separately, then coupling them to calculate the field at the sensor. This approach avoids unnecessary and time-consuming recalculations of RT in unperturbed layers in Jacobian evaluations. The Markov chain method is used to model RT in the aerosol/Rayleigh mixed layer and the doubling method is used for the uniform layers of the atmosphere–ocean system. Our optimization approach has been tested using radiance and polarization measurements acquired by the Airborne Multiangle SpectroPolarimetric Imager (AirMSPI) over the AERONET USC_SeaPRISM ocean site (6 February 2013) and near the AERONET La Jolla site (14 January 2013), which, respectively, reported relatively high and low aerosol loadings. Validation of the results is achieved through comparisons to AERONET aerosol and ocean color products. For comparison, the USC_SeaPRISM retrieval is also performed by use of the Generalized Retrieval of Aerosol and Surface Properties algorithm (Dubovik et al., 2011). Uncertainties of aerosol and nLw retrievals due to random and systematic instrument errors are analyzed by truth-in/truth-out tests with three Chl a concentrations, five aerosol loadings, three different types of aerosols, and nine combinations of solar incidence and viewing geometries.


Introduction
Aerosols exist in the form of airborne suspensions of tiny particles that scatter and absorb sunlight, leading to significant impacts on Earth's energy and water cycles.Quantifying aerosol influences on climate requires accurate determination of their abundances and optical/microphysical properties, which are highly variable spatially and temporally.Aerosol characterization is also crucial for ocean color remote sensing, as the spectral water-leaving radiances account for only 10-15 % of the signal observed at the top of the atmosphere (TOA) and most of the signal arises from atmospheric scattering.Chlorophyll a concentration, colored dissolved organic matter (CDOM), and other ocean optical properties retrieved from spectral water-leaving radiance provides a measure of ocean productivity and health of ocean ecosystem.Small over-or underestimates of the aerosol contribution can bias the determinations of these quantities.
F. Xu et al.: Joint retrieval of aerosol and water-leaving radiance When aerosol is a major target of retrieval, waterleaving radiance is often empirically estimated or even neglected in operational algorithms employed by currentgeneration satellite imagers due to their small contribution to TOA signals.For examples, the MODIS Collection 5 algorithm uses zero water-leaving radiance for all but the 550 nm band, where a value of reflectance 0.005 is assumed (ρ 550 nm = 0.005, cf.Remer et al., 2005Remer et al., , 2006)).MISR and POLDER assume zero water-leaving radiance in the red and NIR bands for aerosol retrieval (Kahn et al., 2010;Deuzé et al., 2000).Then observations are searched within a lookup table (LUT), which contains precalculated radiation fields for a limited number of aerosol models.The aerosol model with lowest fitting residue is selected as the solution.Depending on the sensitivity of measurements, different types and combinations of aerosol models have been designed.For example, MODIS LUT has 20 combinations of fine and coarse aerosol models for retrieving aerosol information from multispectral radiance-only observations (Remer et al., 2005(Remer et al., , 2006)); MISR LUT has 74 aerosol mixtures for retrieving multiangle multispectral radiance-only observations (Kahn et al., 2010); POLDER LUT consists of 12 aerosol models for retrieving multiangle, multispectral, and polarimetric observations (Deuzé et al., 2000).
As the main disadvantage of LUT approach, the solutions have to be selected from a finite number of aerosol models which might not be sufficiently representative in the relevant parameter space.New research efforts have been proposed to expand the LUT to cover more aerosol models (e.g., Limbacher and Kahn, 2014).An alternative to the LUT approach is optimization-based retrieval.It involves a direct inversion of the measurements within the context of a parametric description of the aerosol and surface characteristics that govern the radiation field observed at the TOA.The optimization-based retrieval is featured by a more compact and continuous representation of the relevant parameter space.A review of modern aerosol retrieval algorithms used by airborne and satellite-borne passive remote sensing instruments has been recently given by Kokhanovsky (2015).
When water-leaving radiance becomes the major target of retrieval, traditional retrievals decouple the atmosphere and surface using "atmospheric correction" procedures.The Ocean Biology Processing Group (OBPG) uses the atmospheric correction developed by Gordon and Wang (1994) and Gordon (1997) and refined by Ahmad et al. (2010).In this algorithm an aerosol optical property lookup table is built for ten aerosol models and eight relative humidity (RH) values based on the aerosol property statistics from Aerosol Robotic Network (AERONET) observations (Ahmad et al., 2010).Aerosol optical depth (AOD) and type are determined by fitting the observations in two near-infrared bands (e.g., 748 and 869 nm for MODIS), where water-leaving radiance is assumed negligible.The selected aerosol model is then extrapolated to shorter-wavelength visible bands and applied to the measured TOA radiances to retrieve normalized water-leaving radiance (nLw) (Gordon and Wang, 1994;Gordon, 1997).To reduce errors caused by this atmospheric correction procedure and instrumental radiometric uncertainties, empirical gain factors are derived by forcing agreement between retrieved nLw values and in situ measurements obtained at the Marine Optical Buoy (MOBY) site in Lanai, Hawaii (Franz et al., 2007).
For single-angle, nonpolarimetric instruments such as MODIS and the Sea-Viewing Wide Field-of-View Sensor (SeaWiFS), Franz et al. (2007) pointed out that, "the performance of satellite-based ocean color retrieval process is relatively insensitive to the aerosol model assumption . . .at least for open-ocean conditions where maritime aerosols dominate and aerosol concentrations are relatively low (i.e., aerosol optical thickness generally less than 0.3 at 500 nm)".Therefore, the gain factors derived from conditions at the MOBY site can be applied globally to improve the agreement between satellite and in situ nLw over deep (Case 1) waters.In more challenging observing conditions, e.g., in the presence of absorbing aerosols or complex, spatially diverse (Case 2) waters, inaccurate knowledge of the absorbing aerosol optical properties or height distribution can lead to incorrect assumptions regarding CDOM and phytoplankton absorption coefficients (Moulin et al., 2001;Schollaert et al., 2003;Banzon et al., 2009).In addition, the vertical distribution of absorbing aerosols can affect the reflectance of the oceanatmosphere system, resulting in errors in nLw (Duforêt et al., 2007).In coastal regions, where the traditional assumption of zero water-leaving radiance in the near-infrared (NIR) (Gordon, 1997;Siegel et al., 2000) breaks down, backscattering from suspended hydrosol particles (e.g., algae or sediment) can be misinterpreted as aerosols, leading to overestimation of AOD.The resulting overcorrection can lead to underestimated or even negative water-leaving radiances in the blue and green (e.g., Hu et al., 2000;Bailey et al., 2010;He et al., 2012).
The National Aeronautics and Space Administration's Pre-Aerosol, Clouds, and ocean Ecosystem (PACE) mission, with an anticipated launch date early in the next decade, is aimed at expanding upon current satellite ocean color measurements.The PACE payload is envisioned to include an ocean color spectrometer to measure ocean carbon storage and ecosystem function, and possibly a multiangle, multispectral polarimeter to provide advanced data records on clouds and aerosols and to assist with atmospheric correction of the ocean biology measurements.The capability of multiangle polarimetry in characterizing aerosols for the purposes of assessing their climatic or environmental impacts and improving nLw retrievals over turbid waters or in the presence of absorbing (dust or carbon-containing) aerosols motivates supplementing the vicarious calibration and LUT-based atmospheric correction procedures with one that permits simultaneous extraction of AOD, particle properties, and nLw.Inclusion of spectral bands covering the UV, visible, NIR, and shortwave infrared (SWIR), multiple view angles, and polarimetry in the retrieval enables retrieval of aerosol types that may be beyond the capabilities of the LUT and potentially improves accuracy of both the aerosol and ocean water properties.Given that measurements of atmospheric mineral dust and carbonaceous aerosols show a strong spectral dependence of absorption coefficient in the near-UV (e.g., Koven and Fung, 2006;Bergstrom et al., 2007;Russell et al., 2010) and have a spectral signature similar to those of CDOM, accurate modeling of radiative transfer (RT) in the coupled atmosphere-ocean system (CAOS) becomes necessary.
Without using bio-optical models, some RT models for CAOS consider specular reflection by assuming a flat ocean surface (Jin and Stamnes, 1994;Bulgarelli et al., 1999;Chami et al., 2001;Sommersten et al., 2009;Zhai et al., 2009) for simplicity.Better modeling fidelity and accuracy is then achieved by including sea surface roughness into the RT models (Nakajima and Tanaka, 1983;Fischer and Grassl, 1984;Masuda and Takashima, 1986;Kattawar and Adams, 1989;Mobley, 1994;Deuzé, 1989;Jin et al., 2006;Spurr, 2006) and including the water-leaving radiance and/or ocean foam reflection based on a Lambertian or a more general bidirectional reflectance distribution model (Koepke, 1984;Lyapustin and Muldashev, 2001;Mobley et al., 2003;Sayer et al., 2010;Sun and Lukashin, 2013;Gatebe et al., 2005).Though empirical parameterization of water-leaving radiance simplifies the radiative transfer, the relationship between water-leaving radiance and inherent optical properties (IOP) of dissolved or suspended ocean constituents is indirect.Such a weakness can be overcome by using bio-optical models to relate IOP directly to water-leaving radiance.The bio-optical model-based RT methods make it feasible to perform a one-step retrieval of IOP and aerosol optical properties from TOA measurements of radiance and polarization (e.g., Hasekamp et al., 2011), which is a complementary retrieval strategy to the prevailing two-step retrieval that obtains nLw from TOA via atmospheric correction and then determines IOP from nLw (IOCCG, 2006).Various RT solutions involving the use of bio-optical models have been developed and can be used for this purpose.These include the invariant imbedding method adopted by HydroLight (Mobley, 2008) and its faster version EcoLight (Mobley, 2011a) for scalar (intensity only) RT, in addition to the addingdoubling method (Chowdhary et al., 2006) and successiveorder-of-scattering method (Zhai et al., 2010) for polarized RT in the CAOS.
Joint retrieval of aerosol and nLw properties requires supplementing the forward RT calculations with a sophisticated and computationally efficient inverse model to disentangle their contributions to TOA radiometry and polarimetry.Motivated by the development of a multiangle imaging polarimeter at JPL -the Airborne Multiangle SpectroPolarimetric Imager (AirMSPI) (Diner et al., 2013) -this paper describes the development of a coupled aerosol-ocean retrieval methodology.Our method (1) employs a simplified bio-optical model to obtain a reasonable estimate of nLw in the first retrieval step, followed by an empirical refinement in the subsequent step; (2) applies constraints on the spatial smoothness of aerosol and Chl loadings across neighboring image patches and spectral constraints on aerosol optical properties and on nLw across relevant bands to improve the convergence and stability of the algorithm; and (3) models and stores the RT fields in the aerosol/Rayleigh mixed layer, the pure Rayleigh-scattering layers, and the ocean medium separately, then couples them to obtain the radiative field at the sensor -thereby enhancing the Jacobian evaluations by reusing RT fields in the unperturbed layers.The Markov chain and doubling methods are applied to the mixed and uniform layers, respectively, to gain computational efficiency.
The parameters of our retrieval include spectrally dependent real and imaginary parts of aerosol refractive index, aerosol concentrations of different size components, mean height and width of aerosol distribution, nonspherical particle fraction, wind speed over ocean surface, and normalized water-leaving radiance.As auxiliary product, aerosol phase matrix is obtained from the retrieved refractive index and normalized size distribution.Throughout the paper, we use the definition of "exact" normalized water-leaving radiance (nLw) given by Morel et al. (2002).It is consistent with the definition adopted by Franz et al. (2007) and Zibordi et al. (2009) and is related to the remote sensing reflectance (R rs ) by R rs = nLw/F 0 , where F 0 is the extraterrestrial solar irradiance.
The paper is organized as follows.In Sect.2, we introduce our development of the RT model that integrates the Markov chain, doubling and adding methods for CAOS.The multipatch retrieval algorithm is described in Sect.3. In Sect.4, a truth-in/truth-out test is performed to assess the retrieval uncertainties for a variety of synthetic scenarios combined from three types of aerosols, five aerosol loadings, three Chl a concentrations, three solar incidence angles, four viewing geometries, and two types of measurement noise.To test the algorithm with real data, retrievals applied to AirMSPI observations over the USC_SeaPRISM AERONET site and near the La Jolla AERONET site are compared to the independent AERONET results.A summary is presented in Sect. 5.

A integrated radiative transfer model for a coupled atmosphere-ocean system
A five-layer model is established for a CAOS system, which consists (from the bottom up) of the ocean medium, the airwater interface, a pure Rayleigh layer, an aerosol/Rayleigh mixed layer, and a second pure Rayleigh layer (see Fig. 1).All layers are vertically homogeneous except for the mixed layer, where the aerosol has its own vertical distribution profile, different than that of the Rayleigh-scattering molecular atmosphere.Parameterizations of distribution profile, size, and single scattering properties of aerosols in the mixed layer are demonstrated in Appendix A.

Radiative transfer modeling and Jacobian evaluation strategy
The five-layer CAOS model allows the use of different RT methods to model radiative transfer in different layers based on their computational strength.As an example, the Markov chain method (Esposito and House, 1978;Xu et al., 2010Xu et al., , 2011Xu et al., , 2012)), which exhibits high computational efficiency for modeling RT in vertically inhomogeneous media (Esposito, 1979), is adopted in this work for the aerosol/Rayleigh mixed layer (see Appendix B for details).The doubling method (Stokes, 1862;van de Hulst, 1963;Hansen, 1971;de Haan et al., 1987;Evans and Stephens, 1991; among others), which exhibits high efficiency for modeling RT in optically homogeneous media (Esposito, 1979) is used for the two pure Rayleigh layers and the ocean medium (assumed to be homogeneous throughout the paper).Appendix C gives an example of using the doubling method for modeling RT in the ocean medium.The radiative fields from all layers are then coupled using an adding strategy to obtain the TOA fields.
In addition to the benefit of enabling a combination of the strengths of different RT methods, the strategy of separate RT modeling in five layers also makes for an efficient optimization-based retrieval.During the iterative optimization process, Jacobians are calculated to represent how the radiation fields vary as a function of the model parameters.When they are evaluated by perturbing a model parameter within one of the layers, the diffuse RT fields for all other layers are unchanged from the values obtained from the forward RT simulation and thus can be reused.For example, calculation of the Jacobians with respect to surface or ocean bio-optical parameters does not require recomputation of RT in the atmospheric layer because it has already been derived from the previous forward model calculation.Similarly, when evaluating Jacobians with respect to the aerosol parameters, it is unnecessary to repeat the RT computation of the Rayleigh layers and in the ocean or at the air-water interface.Because optimization-based retrievals involve Jacobian evaluations for a large number of parameters at all iterative steps, this strategy significantly improves the retrieval efficiency.

Atmosphere-ocean coupling
For the atmosphere, the diffuse radiative fields for the aerosol/Rayleigh mixed layer and the pure Rayleigh layers are computed by Markov chain and doubling methods, respectively, then coupled to get the diffuse reflection and transmission matrices for the whole atmosphere (see Appendix B).
For the ocean, the radiative field in the bulk medium is computed by the doubling method with the optics of ocean constituents evaluated by use of a simplified biooptical model (see Appendix D), then coupled with the reflection and transmission across the air-water interface, and finally corrected to account for Raman scattering (see Appendix C).In addition to the contribution by water-leaving radiance from the simplified bio-optical model, total light leaving ocean surface also includes polarized specular reflection (R W , see Appendix E), a Lambertian term for depolarizing ocean foam reflection and an empirical Lambertian correction term, a WL , to account for the errors of the single-parameter-based bio-optical model for water-leaving radiance (i.e., departures from the predetermined functional relationships to Chl a concentration).Thus, the overall bidirectional ocean surface reflection matrix R surf is described as , 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 W by f foam = 2.95×10 −6 ×W 3.52 (Koepke, 1984), and R Bio WL is the reflection matrix of the ocean-interface system with Raman-scattering correction (see Appendix C).Note that R Bio WL is a physically based term in which Chl a concentration ([Chl_a]) is an adjustable free parameter.The last two terms of Eq. (1) constitute our water-leaving radiance model.With and without assuming a WL to be 0 the simplified and the empirically adjusted bio-optical models are formulated, respectively.Though the water-leaving radiance model in Eq. ( 1) has angular dependence, to be consistent with the conventional ocean color products we derive from R Bio WL and a WL in Eq. ( 1) the normalized water-leaving radiance by setting the Sun at zenith (θ 0 = 0 • ) and viewing angle to be nadir (θ v = 0 • ): nLw = (2) where d 0 is the Earth-Sun distance at which the solar irradiance 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. ( 1)-( 2) are all spectrally dependent.Once the diffuse reflection and transmission matrices of the atmosphere and reflection from ocean system are individually known, their coupling to get RT field for the full CAOS is implemented by using the adding method.Two operators Q and S are defined to account for the interaction between the ocean and atmosphere via single and higher orders of reflection, respectively: where R surf is the diffuse reflection matrix from ocean surface and R * atmos is the diffuse reflection matrix from atmosphere with light illumination from the bottom of the atmosphere.The matrices for downwelling and upwelling diffuse light at the atmosphere-ocean interface are given by where µ 0 = cos θ 0 .The reflection matrix of the full CAOS is where µ = | cos θ |.For simplicity in describing the conceptual scheme, the superscript "m" that denotes Fourier series order was not shown in the above expression.In actuality, the TOA radiation fields are reconstructed from all orders of Fourier terms: CAOS,31 cos mφ (4c) where the bidirectional reflectance factor BRF tot and tot are used to fit the observation.Since the sunlight is unpolarized, other matrix entries (namely R CAOS,ij , with j ≥ 2) are not involved in Stokes vector calculation for the diffuse light from the reflection matrix.
Note that the above formalism for modeling RT in a CAOS assumes a horizontally homogeneous atmosphere above a uniform surface, which is known as the independent pixel/patch approximation (IPA) in RT theory (Cahalan et al., 1994).In reality, however, aerosol properties and surface reflection vary across the pixels/patches.To reduce the IPA errors, the single-scattering contribution to the total field evaluated by Eq. ( 4) is replaced by an exact evaluation of radiance along the line of sight.Moreover, for simplicity of model demonstration, our five-layer model assumes the sensor to be located at the TOA.For real airborne measurements, however, the sensor is located inside the atmosphere.Therefore to improve the modeling accuracy, the radiative field is actually computed at the sensor location.This is realized by adding an extra Rayleigh layer above the sensor altitude (e.g., h > h AirMSPI = 20 km in our case), then using the U term in the adding method to compute the diffuse upwelling light reaching the sensor.Moreover, ozone correction is made by BRF tot, corr (λ , where τ ozone is the total ozone optical depth and f ozone is the fraction of ozone above the sensor (in our current study f ozone is assumed to be 20 % for h AirMSPI = 20 km).
The integrated RT model established in the current section will be used as the forward model in retrieval, which is to be introduced in the next section.
3 Optimization approach for joint aerosol and water-leaving radiance retrieval Within the framework of optimization-based retrievals for nonlinear problems, various approaches have been proposed to invert passive remote sensing data for aerosol, ocean, and surface properties.Ideally, the solution vector x that con-F.Xu et al.: Joint retrieval of aerosol and water-leaving radiance tains all relevant parameters characterizing aerosol properties, water-leaving radiance, and surface reflection is approached in an iterative way by x k+1 = x k − x k with x k being the solution after k iterations and x k being the increment being obtained by x k = (J T k ) −1 y k , where J k is the Jacobian matrix evaluated with x k , and y k is the difference between model and measurement ( y k = y(x k ) − y meas ).Unfortunately, the determinant of J k is often close to 0 and as a result J k is ill conditioned.Therefore, a stable retrieval that ensures convergence to a physically sensible solution must impose constraints such that det , where C f is the covariance matrix of the measured signals, W k,i denotes the imposed various constraints, γ k is a Lagrange multiplier that assigns a weight to the constraint, and y k incorporates y k and the relevant a priori constraints and Lagrange multipliers.Introduction of various types of constraints and/or an a priori estimate of W, and establishment of a means for determinant γ k are key elements of optimization-based algorithms.Different approaches include the Levenberg-Marquardt algorithm (Levenberg, 1944;Marquardt, 1963), the Phillips-Tikhonov-Twomey algorithm (Phillips, 1962;Tikhonov, 1963;Twomey, 1963Twomey, , 1975)), and the Twomey-Chahine algorithm (Chahine, 1968), as discussed by Dubovik et al. (2004).
To maximize the use of information provided by different remote sensing instruments on aerosol and surface properties, various algorithms have been applied to inverse radiance and polarimetric signals (Kokhanovsky, 2015;Kokhanovsky et al., 2015).For the particular application of AirMSPI aerosol and water-leaving radiance retrievals, an adaptation of the inversion approach of Dubovik (2004) and Dubovik et al. (2008Dubovik et al. ( , 2011) ) is used.This approach considers inversion as a multiterm least square fitting with multiple a priori constraints.Particularly, as suggested by Dubovik et al. (2008Dubovik et al. ( , 2011)), additional constraints on temporal or spatial variability of the retrieved characteristics can be used if the retrieval is performed for a group of observed pixels/patches.In the present application, a smoothness constraint is imposed to constrain spatial variation of aerosol properties and Chl a concentration over a target area of finite size.While the term "multipixel algorithm" is introduced by Dubovik et al. (2011) for POLDER/PARASOL retrievals with pixel data of ∼ 6 km × 7 km resolution at nadir, the term "multipatch algorithm" is used here since the AirMSPI pixel resolution is much finer (10 m × 10 m) and 50 × 50 pixels are merged into a "patch" to reduce IPA errors.Moreover, as an extension of what is meant by multispectral and multiangle, even polarimetric, a multipixel algorithm can be understood as one based on a forward signal model that can predict how radiances escaping from different pixels are physically coupled, which is tantamount to using 3-D RT (cf.Langmore et al., 2013, for a background-aerosol and gas-plume retrieval demonstration).To avoid confusion, we use the terminology "multipatch" here.Note that though accurate forwarding RT modeling with multiple aerosol species is possible, the increased number of free parameters challenges the ability to retrieve a globally optimized solution in an efficient way.Therefore, as described in Appendix A, a single aerosol species is assumed to represent an effective set of aerosol optical properties, size distribution (which may be multimodal), and vertical profile.Five log-normal size distribution components (N sc = 5) are used to represent the aerosol size distribution, with median radii and standard deviations optimally chosen and given in Table 1, and size-independent refractive index are assumed.Retrieval with more than five size components has also been performed and comparison shows that they both retrieve well aerosol optical properties after being optimally set as log-normally shaped (Dubovik et al., 2006).Since fivecomponent-based retrieval is faster, it is adopted in the current study.Nevertheless, our retrieval leaves the option open for adopting more than five components as well as for retrieving size-dependent refractive index when extra constraints or sensitivity from observation in some observation cases are available.
In the next three subsections, we will give some details on the design of a multipatch retrieval algorithm for joint aerosol and water-leaving radiance retrieval.Readers not interested in it could skip over them.

Multipatch retrieval algorithm with smoothness constraints
Imposing smoothness constraints on both the spatial variations of aerosol loading and Chl a concentration and on spectral variations of aerosol optical properties and nLw leads to the minimization of the following cost function in fitting an N-patch image (Dubovik et al., 2011).
where x i is an iterative solution for the set of parameters being retrieved and x * i is an a priori estimate of the solution corresponding to the ith patch, x = [x 1 , x 2 , x 3 , . . .x N ]; f (x i ), s (x i ) and a (x i ) correspond to the residues of fitting observations, the spectral smoothness constraints, and the a priori estimate, respectively; s,i is a smoothness matrix for constraining the spectral variation of aerosol optical properties and water-leaving radiances across the relevant bands; W f and W a are the weighting matrices for measurements and the a priori estimate, respectively; γ denotes the relevant Lagrange multipliers; y i is the difference between the model and measurements for the ith patch [ y i = y(x i ) − y meas ]; and interpatch is the interpatch smoothness matrix constructed for the patches along two orthogonal directions (u and v) of the image, namely where the derivative matrix S (m) is constructed from the mth order difference, and γ u and γ v are the Lagrange multipliers and their values are shown in Table 2 for all retrieval parameters.
The optimal solution is approached in an iterative way so that after k iterations, the solution vector x i,k+1 containing parameters of aerosol and surface properties for the ith patch is updated as where the multiplier t p (0 ≤ t p ≤ 1) is introduced to improve the convergence of the nonlinear numerical algorithm (Orega and Reinboldt, 1970).Solving the following normal system constructed for the N -patches image at the kth iteration gives the increment of solution for each patch ( x i,k ), where the Fisher matrix for the ith patch is a function of Jacobian matrix J i,k and weighting matrix W f,i , and ∇ (x i,k ) is the gradient of the minimized quadratic form: , where y meas contains the measurement data, y k contains the modeled radiance and polarization with x k , W f is the weighting matrix defined as the covariance matrix C f normalized by its first diagonal element namely W f = (1/σ 2 SD,1 )C (with σ SD being the standard deviation), W a is the weighting matrix of the a priori estimate x * , and s is the single-patchbased smoothness matrix containing sub-smoothness matrices for all parameters.The Lagrange multipliers γ s reflects the strength of the smoothness constraints.
As listed in Table 2, the parameters of the retrieval include spectrally dependent real (m r ) and imaginary (m i ) parts of aerosol refractive index, aerosol concentrations of all size components (C v (r m )), mean height (h a ) and half width (σ a ) of aerosol layer, nonspherical particle fraction (f ns ), wind speed over ocean (W ), Chl a concentration ([Chl_a]), and a WL , which adjust the nLw values in the second step of the retrieval.These parameters form the solution vec- where the natural logarithm is used to ensure nonnegativity of the real solution after dynamic positive or negative changes during the iterative optimization process.The term a WL, Const is an offset determined from nLw using [Chl_a] from the first retrieval step to ensure that the adjustment of nLw in logarithmic space is real.Then γ s s is constructed as a block matrix from diagonal concatenation of the spectral smoothness matrices for real and imaginary parts of the refractive index and a λ , namely for all patches: γ s s =diag 0, 0, 0, 0, 0, γ s,m r s,m r , γ s,m i s,m i , 0, 0, (11) 0, 0, 0, 0, γ s,(a WL,Const + a WL ) s,(a WL,Const + a WL ) , where 0 represents a zero submatrix for a parameter not being subject to any smoothness constraints; and the Lagrange multipliers γ s are predetermined and given in Table 2.
a Nonspherical (spheroidal) particle fraction is excluded in truth-in/truth-out test but included in real data retrieval; b The subscript "1" of nLw means the normalized water-leaving radiance determined from Chl a concentration ([Chl_a] 1 ) retrieved at step 1; c Determined with the consideration of constant offset π/F 0 × (d/d 0 ) 2 × nLw 1 (λ); d The normalized water-leaving radiance is not directly retrieved.After the second step retrieval, the updated Chl a concentration [Chl_a] 2 and the adjustment term a WL are used to derive it via Eq.(2).
In our retrieval test, an a priori estimate is assumed unavailable so we set a i,k = a * i,a .Therefore Eq. ( 10) simplifies to When the spectral and spatial smoothness constraints are turned off (namely setting γ s = γ u = γ v = 0), the multipatch algorithm reduces to the traditional Levenberg-Marquardt algorithm (Levenberg, 1944;Marquardt, 1963), which has been used for retrieval tests with MISR synthetic radiances (Diner et al., 2011;Xu et al., 2012).
Ideally, the retrieval is deemed successful when the minimization of the cost function is achieved, such that 2 where N f,i , N s,i , N a,i and N a * ,i are the number of observations, spectral smoothness, number of unknowns, and a priori estimates of parameters corresponding to ith patch, respectively; N interpatch is the number of spatial smoothness constraints; and ε 2 f is the expected variance due to measurement errors.In practice, forward RT modeling error and other unmodeled effects can impede realization of the condition shown in Eq. ( 13).Therefore, the retrieval is also terminated when the relative difference of fitting residues with solutions from two successive iterations drops below a user-specified threshold value, ε 2 c .Namely, is the second criterion to terminate the optimization.

Determination of Lagrange multipliers
Following Dubovik and King (2000), the Lagrange multipliers reflecting the strength of smoothness constraints are defined as where ε 2 f , ε 2 a and ε 2 g are the first diagonal elements of the covariance matrices corresponding to the measurements (C f ), to the a priori estimates (C a ) and to the smoothness constraints (C g , with the subscript "g" indicating the spectral smoothness constraint "s" or spatial smoothness constraint "u" or "v"), respectively.To estimate ε 2 g for a given parameter to be retrieved (x j ), which is a function of t, the most unsmooth-known solution x ns j (t) over the target area is used, namely, where t min and t max specify the lower and upper bound of t.
In practical implementation of our algorithm, however, the Lagrange multipliers are modified in the following way: There are two differences between γ Final ... and γ ... : 1.The multipliers N f /N g and N f /N a are introduced to account for possible redundancy of the measured and a priori data.Considering that ε 2 ... is a variance of the error in a single measured or estimated a priori value, if we have N values of similar kind the total variance increases proportionally to N. Introducing this coefficient ensures that when there are several kinds of data, the data with fewer values are given comparable weight as the data type for which there is a greater number of available values.
2. The multiplier ε 2 f /ε 2 f is introduced with ε 2 f estimated as the dynamic fitting residual during iterations: .
With the multiplier ε 2 f /ε 2 f , the fitting residual is used as an estimate of measurement error variance.As a result, during the first few iterations the contribution of the a priori term is strongest, and its influence decreases as the retrieval progresses.This is done to ensure mostly monotonic convergence, as in the Levenberg-Marquardt procedure (Levenberg, 1944;Marquardt, 1963).However, the Levenberg-Marquardt approach does not specify a particular scheme for introducing these terms, rather it relies on the implementer's intuition.Our algorithm requires the fitting errors in the initial iterations to be dominated by model linearization errors as opposed to random measurement errors.Because at each iterative step the full forward model is replaced by its linear approximation, the errors of linearization decrease as convergence toward the final solution progresses, and they practically disappear so that ε 2 f becomes equal to ε 2 f .As a result of this adjustment of the Lagrange multiplier, the nonlinear iteration becomes significantly more monotonic.

Implementation of two-step retrieval
As water-leaving radiance is a small contribution to TOA signals, opening a large number of parameters for its retrieval increases the risk of obtaining solutions at local minima of the fitting metric and a significant slowdown of the retrieval.To improve retrieval efficiency and reliability, we use a two-step retrieval strategy: obtaining a reasonable estimate of water-leaving radiance (i.e., close to the truth) by using a bio-optical model constrained by a single parameter (namely Chlorophyll a concentration, [Chl_a] which governs the abundance of CDOM and phytoplankton in a prescribed way) during the first step of the retrieval.This is accomplished by setting a WL to zero so that only Chl a concentration (the ocean parameter to which the measurements have the largest information content) is retrieved.Other ocean parameters (e.g., CDOM concentration) are models as dependent on [Chl_a].In light of the possibility that the biooptical model parameterized by Chl a concentration only can have inaccuracies (particularly in Case 2 waters), this constraint is relaxed in a subsequent step so that the nLw retrieval is improved by letting the Chl a concentration and the a WL term be optimized simultaneously ( a λ,WL is allowed to be negative).To mitigate the propagation of instrumental and atmospheric modeling errors to the water-leaving radiance, the second retrieval step (1) allows the adjustment of the bio-optical model-based nLw values only within a confined range (e.g., −15 % ≤ nLw adjust /nLw 1 ≤ +15 %, with nLw 1 being the nLw from the first retrieval step); and (2) imposes a spectral smoothness constraint on nLw(λ).

Validation of optimization algorithm
Technologies to extend the observational capabilities of JPL's Multi-angle Imaging SpectroRadiometer (MISR, Diner et al., 1998) have been developed over the past decade for the purpose of providing additional observational constraints on aerosol and surface properties.These have been incorporated into AirMSPI, as described in Diner et al. (2013).AirMSPI is an ultraviolet-visible-near-infrared imager that has been flying aboard the NASA ER-2 high-altitude aircraft since October 2010.At the heart of the instrument is an 8-band (355, 380, 445, 470, 555, 660, 865, and 935 nm) pushbroom camera mounted on a gimbal to acquire multiangle observations over a ±65 • along-track range.Three of AirMSPI's spectral bands (470, 660, and 865 nm) include measurements of the Q and U Stokes polarization parameters.To validate the retrieval approach, the algorithm was applied to simulated and real AirMSPI data.

Retrievals with simulated AirMSPI observations
Prior to performing retrievals with actual AirMSPI data, truth-in/truth-out tests with simulated data were conducted to assess the accuracy and stability of our optimization approach.The simulation generates modeled TOA radiance and polarization fields based on AirMSPI observations over the USC SeaPRISM AERONET-OC site (118.12• W, 33.56 • N) off the coast of southern California on 6 February 2013.Images of the targeted area were obtained at 9 viewing angles (0, ±29, ±47, ±59, and ±65 • ).At nadir, the imaged area covers 10 km × 11 km swath.The data are mapped to a 10 m spatial grid.Patches comprised of averages of data within 50 pixel × 50 pixel areas were generated, and a total of 102 patches seen at all angles, corresponding to a 5 km × 5 km area, were used simultaneously in the retrievals to take advantage of the multipatch retrieval algorithm.Totally 126 signals per patch are measured, which include radiances at nine angles and eight spectral bands and Q and U at nine angles and three polarimetric bands.Since we use DoLP in retrieval and did not model or make use of AirMSPI's watervapor band at 935 nm, in fact we have 90 signals per patch.Moreover, patch-averaged radiance and degree of linear polarization (DoLP) are used in retrieval.The algorithm tests include three steps: 1. Using the AirMSPI observational characteristics described above, simulated measurements were generated for five different aerosol loadings, three aerosol types, three Chl a concentrations, and nine combinations of Sun illumination and viewing geometries.The five aerosol loadings correspond to AOD of 0.02, 0.1, 0.3, 0.5, and 1.0 in the AirMSPI green band (555 nm).
The three aerosol types include (a) weakly absorbing aerosols from the MODIS/SeaWiFS LUT (Ahmad et al., 2010) with RH = 85 % and fine mode volume fraction = 50 %; (b) moderately absorbing particles from the same LUT with RH = 30 % and fine mode volume fraction = 80 %; and (c) dust aerosols (Sokolik and Toon, 1999).Hygroscopic growth is assumed for the water-soluble and smoke aerosols but is excluded for dust aerosols.The refractive index, size parameters, and vertical profile parameters for these three types of aerosols, and the assumed wind speed, are listed in Table 3.The size distributions of the first two aerosol types were fitted by our five-component aerosol size model.The three Chl a concentrations used were 0.05, 0.2, and 1 mg m −3 .A perturbation of ±10 % was imposed on the water-leaving radiance predicted by the Chl abased bio-optical model to simulate modeling errors and to test the validity of the two-step retrieval strategy.
The wind speed was assumed to be 4 m s −1 .The mean height and half width of the aerosol distribution profile were set to 1 and 0.75 km, respectively.
To cover a wide range of observing geometries, a total of nine scenarios based on the AirMSPI USC_SeaPRISM viewing geometry is used, as illustrated in Fig. 2: the Sun is placed at the original incidence angle θ 0 = 49.1 • as well as at 25 • and overhead Sun (θ 0 = 0 • ).Relative azimuth angles of φ ≈ 50, 95, 140, and 176 • are also modeled.The latter case includes glint.For the case with overhead Sun, only one azimuth angle is necessary.
2. Random noise was added to the simulated radiance and DoLP values.This is a commonly adopted measure to test the impact of measurement errors on retrieval algorithm performance (Dubovik et al., 2011;Hasekamp andLandgraf, 2005, 2007).We added a relative measurement uncertainty of σ L = ±1 % to the radiances and an absolute uncertainty of δ DoLP = ±0.005 to the DoLP.After a random-error test, an extra ±4 % systematic error was added to study the influence of calibration bias.
3. Retrieved aerosol properties and Chl a concentrations were compared to their known (input truth) values.

Influence of aerosol loading and absorption on nLw retrieval
As an example, we use one of the simulated scenarios of AirMSPI observation over USC_SeaPRISM AERONET-OC site (θ 0 = 25 • , φ ≈ 95 • ) as input.Figures 3-6 compare retrieved AOD, SSA, particle size distribution (PSD), and nLw,  igure 2. Simulation geometries based on AirMSPI observations over the AERONET OC-site USC_SeaPRISM on 6 February 2013.The three red dots indicate the Sun's location θ 0 = 49.1 • , the actual value at the time of the AirMSPI overflight, as well as 25 and 0 • .For each incidence angle, four viewing geometries corresponding to the azimuthal angles ≈ 50, 95, 140, and 176 • are simulated, which are marked in different colors: black, blue, dark red, and dark yellow, respectively.Due to symmetry, only one azimuthal plane is necessary to simulate for zenith Sun location.Therefore totally nine geometries are created for truth-in/truth-out test.The viewing angles corresponding to the 9 AirMSPI images form line segments.Each line segment is composed of densely sampled cross-track positions contributed by all patches in the image.For each azimuthal case, a total of nine segments are plotted.respectively, to the "true" values used in the simulation.In all figures, the top, middle, and bottom rows of the panels correspond to Chl a concentrations of 0.05, 0.2, and 1 mg m −3 , respectively (with ±10 % perturbation on water-leaving radiances in different bands).The left, middle, and right panels correspond to weakly absorbing, moderately absorbing, and dust aerosols, respectively.
For all aerosol types, the shapes of AOD, SSA, and nLw, as a function of wavelength and PSD as a function of particle radius, are similar to their true values.Due to the limited contribution of nLw to TOA radiance, the aerosol retrieval accuracy is not significantly affected by the Chl a concentration within the range modeled here.The retrievals over dust are less accurate than for the weakly and moderately absorbing aerosols, due to the fact that dust aerosols are dominated by coarse-mode particles and the extinction is more spectrally neutral, so the information provided by the multispectral measurements between 355 and 865 nm is less effective to constrain the aerosol retrieval.As expected, Fig. 6 shows higher retrieved nLw accuracy at low AOD loading (τ 555 ≤ 0.1) due to greater atmospheric transparency and increased fraction of nLw in the TOA signals.
When the aerosol species changes from weakly absorbing aerosols (corresponding to the three figures in the left column of Fig. 6) to moderately absorbing aerosols (middle column), then to dust (right column), the bias in nLw increases.This is because the water-leaving radiance signal becomes weaker with increased atmospheric absorption and retrieval of absorbing aerosol properties is more uncertain than for nonabsorbing aerosols, and the errors propagate to the water- A more comprehensive view of aerosol retrieval errors is displayed in Fig. 7a-d.Though the absolute error of retrieved AOD increases as the aerosol loading increases (see Fig. 7a), the relative error of AOD (100 × |AOD retrieved − AOD true |/AOD true ) generally decreases as the TOA radiance carries more aerosol information at higher loading (see Fig. 7b).For the same reason, an inverse relationship between aerosol loading and absolute error in single-scattering albedo (|SSA retrieved − SSA true |) is observed, as shown in Fig. 7c.To evaluate the retrieval error for size distribution, the effective radius is used and calculated for fine and coarse modes: where the lower size limit r min = 0.04 µm and the upper size limit r max = 15 µm.Setting r cri to be 0.75 µm for weakly and moderately absorbing aerosols and r cri to be 0.25 µm for dust aerosols to distinguish fine and coarse modes, a generally inverse relationship between aerosol loading and the relative error in effective radius of fine and coarse-mode aerosols (100 × |r eff, retrieved − r eff, true |/r eff, true ) is also observed for all types of aerosols, as shown in Fig. 7d.For τ 555 ≥ 0.3, the maximum retrieval error in AOD is ∼ 2.5, 2.5, and 7 % for weakly, moderately absorbing aerosols, and dust particles, respectively.The maximum retrieval error for SSA ω 0,355 nm is ∼ 0.005, 0.015, and 0.025 for weakly absorbing, moderately absorbing, and dust aerosols, respectively.We find that the maximum error in SSA for the weakly absorbing aerosol appears at red and near-infrared bands (660 and 865 nm) for all aerosol loading cases, suggesting that there is less sensitivity to SSA as the ocean reflectance decreases.For the moderately absorbing aerosols, the maximum error is observed at the two UV bands (355 and 385 nm), indicating higher errors as absorption increases, particularly at low aerosol loading.Moreover, increasing AOD is found to be helpful for constraining the SSA retrieval for both weakly and moderately aerosols.However, for dust aerosols, where SSA spans a larger range as a function of wavelength compared to the weakly and moderately absorbing aerosols, limited improvement on SSA retrieval accuracy is gained by increasing AOD. Figure 7d shows that for weakly and moderately absorbing aerosols the effective radius for coarse-mode aerosols has larger retrieval errors than the fine mode aerosol.We attribute this to the fact that the longest spectral band of AirMSPI used in the retrievals (865 nm) is insufficient to fully constrain the coarse-mode aerosol PSD.
In Fig. 7e-f, which correspond to Chl a concentration to be 0.05, 0.2, and 1.0 mg m −3 (with ±10 % perturbation imposed on the water-leaving radiance), the retrieval error of normalized water-leaving radiance ( nLw = nLw retrieved − nLw true ) is plotted against uncertainty metrics specified by the PACE Science Definition Team (SDT) (Del Castillo et al., 2012), i.e., a relative error of 5 % or an absolute error of 0.001 × F 0 /π (whichever is larger) in the visible, and twice these values in the UV.For weakly and moderately absorbing aerosols, the accuracy of nLw at all visible bands mostly falls within the PACE SDT requirement for all aerosol loadings and Chl a concentrations.The uncertainty in retrieved nLw in the pair of UV bands, however, falls outside the specified bounds when τ 555 > ∼ 0.1.As the TOA signals in the UV are dominated by Rayleigh scattering, accurate retrieval of water-leaving radiance remains challenging even after the interpatch smoothness constraints on aerosol variation and spectral smoothness constraints on aerosol optical properties are imposed.For all Chl a concentrations, errors in nLw are largest for dust aerosols, and fall outside the PACE SDT requirement for τ 555 > ∼ 0.1, even in the visible.These errors can potentially be reduced if an improved bio-optical model can be devised that relates the more accurately determined visible nLw values to the values in the UV.
Figure 7 shows that for all aerosol types, even though the retrieval errors of SSA and AOD at low aerosol loading (τ 555 < 0.1) are relatively larger than at high AOD, these errors do not propagate to the retrieval of nLw.This is because in the single-scattering regime, the path radiance is dominated by scattering optical depth, which is the product of AOD and SSA.This means AOD and SSA errors counteract each other to some extent (i.e., an overestimated AOD is compensated by an underestimated SSA and vice versa) so that scattering optical depth is less biased, leading to a reduced impact on the retrieval of nLw.However, when AOD increases, the fraction of water-leaving radiance in the TOA signal reduces significantly, and accurate separation of its weak contribution in the multiple-scattering regime becomes more difficult.The presence of dust aerosols further complicates the retrievals as the aerosols and CDOM share a similar shape of absorption spectra, namely, increasing absorption at shorter wavelengths (Aurin and Dierssen, 2012;Bergstrom et al., 2007).

Effect of multipatch vs. single-patch retrieval
Taking the case of median loading (τ 555 = 0.3) of weakly absorbing aerosols and median Chl a concentration [Chl_a] = 0.2 mg m −3 as an example, Fig. 8a  servation are highly nonunique subjected to local optimum solutions.Through the imposition of interpatch smoothness constraints on aerosol loading and Chl a concentration, the multipatch retrieval yields results that are closer to the truth.As indicated in Fig. 8b, the multipatch algorithm also shows greater noise resistance in all three quantities (nLw, AOD and SSA) simultaneously.The AOD error in the single-patch retrieval decreases as the level of random noise in intensity increases from 0.5 to 2.0 %, due to that fact that the errors mainly propagate into nLw and SSA.

Comparison to direct water-leaving radiance retrieval
For the same scene, parameters used to compare the singleand multipatch-based retrievals, Fig. 9 compares a retrieval using the bio-optical model and one in which nLw is modeled using unconstrained Lambertian reflectance factors at each wavelength.Using the bio-optical model reduces the parameter space for the water-leaving radiance from seven independent spectral values to a single parameter (Chl a concentration) that establishes the spectral variation of the surface signal.While there is little difference between AOD retrieved with and without the bio-optical model, SSA re- trieval accuracy improves by 0.01 and 0.02 at 350 and 865 nm, respectively.Moreover, a remarkable gain in nLw accuracy by about 6, 11, and 12 %, or 0.07, 0.12, and 0.03 mW cm −2 sr −1 µm −1 in absolute magnitude at 445, 470, 555 nm, respectively, is achieved when the bio-optical model is used.Given that the PACE SDT specification tolerates an uncertainty of ∼ 0.06 mW cm −2 sr −1 µm −1 in these bands, the accuracy gain from using the bio-optical model is significant.

Influence of systematic error
The above truth-in/truth-out tests were performed assuming instrumental errors are completely random.Such an assumption, however, is not applicable to radiometric errors and their band-to-band variations, which represent systematic deviations from the true values due to calibration errors.For a satellite instrument such as MISR, the radiometric uncertainty is 4 % and the band-to-band variations are about 1.5 % (Bruegge et al., 2002).Because the absolute error is larger in magnitude than band-to-band error and represents a systematic bias that applies to all measurements, it can potentially have greater impact on retrieval accuracy than band-to-band errors and random noise.To model its effect, we keep the random noise levels used in the previous analysis and add a ±4 % systematic error to the simulated radiance signals.
The resulting retrieval errors of AOD, SSA, effective radii of fine-and coarse-mode aerosol, nLw, and band-to-band ratio are displayed in Fig. 10a-f, respectively.Comparison of Figs. 7 and 10 shows that systematic errors have a larger impact on retrieval accuracy than random errors, as the latter are suppressed by using a lot of patches for retrieval while the former are not.For AOD and SSA, a negative radiance bias causes larger retrieval errors than a positive bias.Comparison of Figs.10e and 7e shows that errors in nLw due to an intensity bias increase at all AODs: at low aerosol loading the errors propagate to nLw while at high loading the contribution of nLw to the TOA signal is weak, exacerbating errors.On the other hand, comparison of Figs.10f with 7h shows a much smaller effect of systematic errors on nLw(λ)/ nLw(555); in other words, the systematic errors mainly propagate to the overall magnitude of nLw(λ) curve while the relative spectral shape is affected to a much lesser degree.

Retrievals with real AirMSPI observations
Following algorithm validation using the truth-in/truth-out tests, we applied the algorithm to actual AirMSPI ob- servations acquired over the USC_SeaPRISM AERONET-OC site and near the AERONET site in La Jolla.The USC_SeaPRISM and La_Jolla scenes were chosen from a larger set of AirMSPI field campaign images to ensure cloud free conditions.The data were processed with the recently upgraded data processing pipeline, which includes vicarious radiometric calibrations and improved polarimetric calibration making use of onboard polarization sources.Nadir intensity and DoLP images from combinations of different spectral bands for these two target areas are shown in Figs.11a  and 12a.Maps of retrieved AOD and SSA at 555 nm, nLw at 445 and 555 nm spectral bands are displayed in Figs.11b  and 12b.
Selecting the image patch that is closest to the AERONET site, our retrieved AOD, SSA, size distribution, and nLw are compared to the independent AERONET results, as shown in Figs.11c and 12c.We first discuss results from the USC_SeaPRISM retrievals.The AERONET site reported a relatively high 550 nm AOD of 0.30 and 0.26 at 19:08 and 20:08 UTC, respectively, and our retrieval returns an intermediate value of 0.27 from the AirMSPI data acquired at 19:40 UTC.The differences between the AirM-SPI and AERONET AOD and SSA retrievals are within the AERONET SSA retrieval uncertainties (e.g., 0.015 for τ 440 and 0.03 for ω 0,440 at τ 440 > 0.2, see Table 4 of Dubovik et al., 2000).The Generalized Retrieval of Aerosol and Surface Properties (GRASP) algorithm by Dubovik et al. (2011Dubovik et al. ( , 2014) ) was also run, and the difference between the GRASP and JPL algorithms is on the order of ∼ 0.025 for AOD and ∼ 0.008 for SSA in all bands.
As illustrated in the bottom right panel of Fig. 11c, the retrieved nLw also compares favorably to AERONET reported values.After interpolating AERONET nLw in logarithmic space to obtain nLw in the AirMSPI bands, the differences are found to be 0.0396, 0.0118, 0.0198, and 0.0077 mW cm −2 sr −1 µm −1 in the 445, 470, 555, and 660 nm bands, respectively.These differences are within the AERONET-OC uncertainties of 0.0462, 0.0516, 0.0279, and 0.0167 mW cm −2 sr −1 µm −1 in the four bands, obtained by interpolating combined standard uncertainties in validated nLw at various AERONET-OC sites (Gergely and Zibordi, 2014).Note that the nonspherical particle fraction retrieved using both GRASP and JPL algorithm is negligible and the results are not displayed here.
For the second study site, the AirMSPI target area was about 13 km away from the La Jolla AERONET station.In spite of the distance, the differences between the AirM-SPI and AERONET AOD and SSA values are both within AERONET's uncertainty, as observed from the upper two plots of Fig. 12c.Though the difference in PSD in some size bins falls outside the AERONET uncertainty range, the bimodality of the size distribution is identified even at the low aerosol loading for this case (τ 555 ∼ 0.04).Independent surface measurements to validate the nLw retrieval were not available at this site.

Summary and outlook
Accurate retrieval of both aerosol properties and waterleaving radiance is challenging as the latter only accounts for a small fraction of TOA signals and can be easily contaminated by Rayleigh and/or aerosol scattering.To ensure highquality retrievals of the aerosol properties, traditional atmospheric correction schemes, which are focused primarily on retrieval of surface characteristics, may not be sufficient.In light of the additional information provided by multiangular, multispectral, and polarimetric measurements, we tested the concept of simultaneous aerosol and water-leaving radiance retrieval which include spectrally dependent real and   ocean medium.These features are implemented to enhance computational efficiency.
Next, an optimization approach has been developed for joint aerosol and water-leaving radiance retrieval.The algorithm involves a two-step retrieval strategy, first relying on a bio-optical model to retrieve a single parameter (Chl a concentration) that governs nLw, then allowing adjustment of nLw to account for modeling errors.Our optimization algorithm imposes smoothness constraints on the spatial variation of aerosol loading and Chl a concentration and the spectral variation of aerosol optical properties and nLw.We demonstrated that the use of multipatch constraints in conjunction with the bio-optical model improves the retrieval accuracy of aerosol properties and water-leaving radiance and stabilizes the algorithm.Truth-in/truth-out tests assuming random errors 1.0 % (relative) and 0.005 (absolute) for intensity and DoLP, respectively, show that the retrieval accuracy of nLw in the visible bands meet the requirements of the PACE SDT in the presence of weakly and moderately absorbing aerosols of optical depth at 555 nm less than 1 and Chl a concentrations 0.05, 0.2 and 1 mg m −3 , whereas meeting the PACE SDT goals in the UV and for dust is more challenging.Increased aerosol absorption reduces the nLw retrieval accuracy except when AOD is low.The addition of systematic errors leads to biases in the absolute magnitude of nLw at both low and high AOD.Band ratios between visible bands (e.g., nLw(λ)/nLw(555)), which are widely used in ocean color analyses, are less impacted by systematic errors for weakly absorbing aerosols.Case studies of AOD, SSA, size distribution and nLw using real AirMSPI observations over the AERONET USC_SeaPRISM OC site and near the AERONET La Jolla site compare favorably to AERONET's reported values.
In future work, the influence of modeling errors on nLw retrievals will be investigated.Since water-leaving radiance accounts for a small fraction of the TOA signals, small forward modeling errors can translate into large nLw retrieval errors.The modeling error can arise from various sources, e.g., neglect of cirrus cloud contamination, approximate treatment of trace-gas absorption and the atmosphere profile, salinity of sea-water, assumption of plane-parallel atmosphere, retrieval of effective aerosol optical properties from assuming single aerosol species and size-independent refrac-  tive index, δ-truncation of phase matrix, finite stream number and truncated Fourier terms adopted in the RT model, or errors in the solar spectrum.Further considering the potential errors in empirically adjusted bio-optical model for optically complex waters (e.g., coastal shallow water and inland water), the combined effects on nLw accuracy remain to be studied.Development of a fast, yet accurate, CAOS RT model and algorithm validation using a wider set of AirMSPI scenes is also part of our ongoing effort.The aerosol/Rayleigh mixed layer is defined to have the minimum altitude h min and maximum altitude h max .A single aerosol species is assumed to be distributed throughout it with a Gaussian distribution profile characterized by mean height a and standard deviation σ a , characterizing the width of the aerosol layer.Then, the aerosol concentration profile c a is where the normalization factor F norm is used to ensure that h max h min c k (h)dh = 1 and evaluates to where erf(x) is the error function.
Breaking the aerosol volumetric size distribution dV (r)/d ln(r) into a finite number of size components (Dubovik et al., 2011), the total AOD (τ a ) is the sum of all size components: where N sc is the total number of size components; K ext,a,i and C v,i are the extinction coefficient (in units of km −1 ) and column volume concentration (in units of km) of the ith aerosol size component, respectively; C v, tot is the total volume concentration (C v, tot = C v, 1 +C v, 2 +C v, 3 +. ..); and f i is the volume fraction of the ith component ( Moreover, the total aerosol size distribution is constituted as The associated normalized size distribution is Using a log-normal volume weighted size distribution for all size components, dv i (r)/d ln r is dimensionless and is parameterized by a median radius for volume size distribution r m,i and a geometric standard deviation σ i , namely, The mixed layer is subdivided into N sublayers, each bounded by the altitudes h n and h n+1 (h n < h n+1 ).Assuming no trace gases and optical homogeneity of each sublayer, the optical depth ( τ (n) ), single-scattering albedo (SSA, ω (n) 0 ) and phase matrix (P (n) ) of the nth sublayer are contributed by aerosol and Rayleigh molecules only, therefore where P R and P a are the Rayleigh and aerosol phase matrix, respectively; the SSA of aerosol ω 0,a is a function of scattering coefficient (K sca,a ) and extinction coefficient (K ext,a ): ω 0,a = K sca,a /K ext,a ; τ (n) a is the AOD in the nth sublayer and can be evaluated analytically after considering the aerosol distribution profile (Eq.A1) according to the following: R in Eqs.(A7)-(A9) is the Rayleigh optical depth of the nth sublayer and is evaluated assuming the US standard atmosphere profile (Tomasi et al., 2005;Bodhaine et al., 2007).
As functions of aerosol refractive index, shape, and size distribution, the elements of P a and the quantities K ext,a and K sca,a are computed using Mie theory for spherical particles (van de Hulst, 1981) and using T-matrix and geometrical optics methods for nonspherical (spheroidal) particles (Dubovik and King, 2000;Dubovik et al., 2006).During the optimization process, the spectrally dependent refractive index (m r + m i i) and concentrations (C v (r m )) of the aerosol size components are updated dynamically.To avoid inefficient on-the-fly Mie computations, these particle properties are precalculated for all size components and saved on a grid of discrete real and imaginary refractive indices.For an arbitrary combination of real and image refractive indices, interpolation is used to obtain the optical properties.Then, the aerosol phase matrix and scattering and extinction coefficients are updated via linear combination of the contribution of all size components, namely where X represents any Mie property of {P a , K ext,a and K sca,a }.Via Eqs.(A7)-(A9), X is then mixed with Rayleigh scattering to obtain the overall scattering properties of each layer, which are used as inputs to the RT model for the mixed layer.distribution, as airborne and satellite-borne passive remote sensing has low sensitivity to the vertical profile.As a consequence of this assumption, the ocean reflection matrix R ocean , which depends on τ ocean , ω ocean , and P ocean , is computed using the doubling method.
As described in Appendix E, reflection of light from ocean surface and its transmission through an air-ocean interface are evaluated using the model of Cox and Munk (1954a, b) for a wind-roughened ocean surface.The set of reflection and transmission matrices (R W , T W ) corresponding to downwelling incident light (in air) and another set of matrices (R * W , T * W ) corresponding to upwelling incident light (in water) are then determined.In accordance with the adding method, two operators Q and S are defined to account for the interaction between the ocean bulk and its interface with air via single and higher orders of reflection, respectively.
However, unlike a real atmospheric layer that attenuates light during its transmission, the air-water interface is a pseudolayer without any thickness, so all attenuation-related terms should be removed.This leads to a modification of the classical adding-doubling scheme (named the "extended adding-doubling method" in the remainder of the paper) for coupling the transfer of radiation between the ocean bulk medium and the air-water interface: the matrices describing the downwelling and upwelling of diffuse light at the top of the ocean now become The reflection matrix describing the upwelling diffusely reflected light leaving the ocean-air interface is where the superscript NR over R indicates that Raman scattering is not considered at this step (but will be included via a correction introduced in the next section).As a numerical validation, Fig. C1 compares top-ofocean radiance and DoLP computed with the extended adding-doubling method via Eq.(C1) and an independent successive-orders-of-scattering code (Zhai et al., 2010).Chl a concentration was set to 0.30 mg m −3 , solar zenith angle to 60 • , surface wind speed to 7 m s −1 , and ocean optical thickness to 10.Using 40 streams in the half plane of 0 ≤ µ ≤ 1 and 30 Fourier terms, this case study shows that the maximum relative difference in computed intensity is < 0.3 % in magnitude, and the maximum absolute difference in degree of linear polarization (DoLP) is 0.005 in the worst case, but more typically about 0.001.The difference can be even smaller by using more streams and Fourier terms.

C2 Correction for Raman scattering
The RT modeling formulated in Section C1 does not account for inelastic scattering processes including Raman scattering by water and fluorescence by chlorophyll and CDOM.Accurate modeling of these processes is necessary (Mobley, 2008;Zhai et al., 2015) but requires additional inputs and computations that can significantly slow down the retrievals (Mobley, 2011b).To optimize the trade-off between computational efficiency and numerical accuracy, the correction scheme proposed by Lee et al. (2013) is used to quantify the contribution With the FF scattering function, the backscattering efficiency can be obtained analytically (Mobley et al., 2002): where δ 90 is δ evaluated at = 90 • .Mobley (2002) found that the detailed shape of particlescattering function is not critical if a correct backscatter fraction B bp is provided.Characterized by a spectrally neutral backscatter efficiency B bp , Huot et al. (2008) obtained an empirical relationship between Chl a concentration and B bp .The spectrally neutral assumption for the backscattering efficiency also indicates that the refractive index and slope parameter are not independent to each other.Knowing B bp from a given Chl a concentration via Eq.(D10) and further assuming a linear relationship between n p and γ p (Mobley et al., 2002), namely, n p = 1.01 + 0.1542(γ p − 3). (D11) Thus, given Chl a concentration B bp is computed from Eq. (D10).Then Eqs.(D9) and (D11) can be solved to determine n p and γ p -the two model parameters of the FF scattering function.Figure D1 illustrates the resulting relationships between n p and B bp , between γ p and B bp , and between n p and γ p .
The particle-scattering coefficients are evaluated based on the model by Morel and Maritorena (2001): The scattering coefficient for CDOM is treated as zero in the present study.

Figure 1 .
Figure1.Depiction of the 5-layer CAOS model.A Gaussian vertical distribution profile for aerosols in the mixed layer is assumed and the Markov chain model is used for RT in this optically inhomogeneous layer.The ocean medium and the two Rayleigh layers (below and above the mixed layer, respectively) are treated as optically homogeneous and the doubling method is used for the RT computations.Coupling of these layers and inclusion of the airwater interface are completed by use of the adding strategy.The Sun illuminates the top-of-atmosphere with solar zenith angle θ 0 and azimuthal plane φ 0 .We define φ = φ v − φ 0 , where the sensor views the atmosphere at viewing angle θ v and azimuthal angle φ v .

Figure 4 .
Figure 4. Panel layout as in Fig. 3 but for retrieved single-scattering albedo.The black line with dots placed at the AirMSPI wavelengths represents the true SSA.The colored symbols represent retrieved SSA for various values of AOD.

Figure 5 .
Figure 5. Panel layout as in Fig. 3 but for retrieved normalized aerosol size distribution.The black lines correspond to the true size distribution, with dots at discrete values of particle radius.The colored lines represent retrieved size distributions for various values of AOD.

F.Figure 6 .
Figure6.Panel layout as in Fig.3but for retrieved values of nLw (mW cm −2 sr −1 µm −1 ).The black lines correspond to the true nLw, with dots placed at the AirMSPI wavelengths.The colored symbols represent retrieved nLw for various values of AOD.

Figure 7 .
Figure 7. Retrieval errors of (a) AOD; (b) AOD (relative difference); (c) SSA; (d) effective radii for fine and coarse-mode aerosols; (e)-(g) nLw (signed difference) corresponding to Chl a concentrations 0.2, 0.05 and 1.0 mg m −3 , respectively (with ±10 % perturbation imposed on the water-leaving radiance); and (h) band ratios (Rλ, nLw = nLw(λ)/ nLw(555)).The retrieval errors of aerosol properties show similar features for all Chl a concentrations.Therefore the results corresponding to [Chl_a] = 0.2 mg m −3 are displayed in (a)-(d).Via truth-in/truth-out tests, the uncertainties are estimated for AirMSPI multispectral, multiangular, and multipolarimetric observations over a 5 km × 5 km ocean area.The simulation is based on nine combinations of Sun incidence and viewing geometries.Relative random noise of 1.0 % is used for radiance and absolute random noise of 0.005 is used for DoLP.The colors correspond to seven different AirMSPI spectral bands.The maximum water-leaving radiance error target specified by the PACE Science Definition Team (SDT) is plotted as black curves.The uncertainty of nLw at 865 nm is not displayed since the PACE SDT did not specify a requirement on this band.The spread of the error, depicted by the vertical bars, reflects the dependence on illumination and viewing geometries.

-Figure 8 .-Figure 9 .
Figure8.Comparison of single-patch-and multipatch-based retrievals of (a) AOD, SSA, aerosol size distribution, and Chl a concentration for the median AOD (τ 555 = 0.3) of weakly absorbing aerosols and Chl a concentration of 0.2 mg m −3 .The simulation uses the Sun and viewing geometry corresponding to the AirMSPI overflight of the USC_ PRISM AERONET site.Image-averaged Chl a concentrations are 0.29 and 0.22 mg m −3 for the single-and multipatch-based retrievals, respectively.A random error of 1.0 % and 0.005 is added to the simulated intensity and DoLP data, respectively; (b) AOD, SSA, and nLw (mW cm −2 sr −1 µm −1 ) retrieved with different levels of random noise (0.5, 1.0, 2.0, and 3.0 %) added to the simulated BRF while the noise in DoLP is kept at 0.005.The aerosol loading, Chl a concentration, and Sun and viewing geometry are the same as in Fig.8a.

Figure 11 .
Figure 11.(a) Nadir AirMSPI intensity image from spectral combination of 445, 555, and 660 nm bands (left image) and DoLP image from spectral combination of 470, 660, and 865 nm bands (right image).The bright spot inside the white circle marked on the intensity image (dark spot inside the white circle marked on the DoLP image) is the AERONET USC_SeaPRISM ocean color station, located on the Eureka oil platform.AirMSPI observations were acquired at 19:44 UTC on 6 February 2013.The yellow frame bounds the area viewed in common from all nine angles (observations over this area are used for retrieval).(b) Maps of retrieved AOD at 555 nm (top left), SSA at 555 nm (top right), nLw at 445 nm (bottom left), and nLw at 555 nm (bottom right).(c) Comparisons of retrieved spectral AOD (top left), SSA (top right), aerosol size distribution (bottom left), and nLw (bottom right) using the GRASP and JPL algorithms to AERONET reported values.

Figure 12 .
Figure 12.Similar to Fig. 11 but corresponding to the AirMSPI observations near AERONET La Jolla site on 14 January 2013 at 21:09 UTC.The bluish part at the bottom right part of the DoLP image indicates the shallow water area which was not captured by all images and hence excluded in retrieval.

F
. Xu et al.: Joint retrieval of aerosol and water-leaving radiance Appendix A: Parameterizations of distribution profile, size, and single-scattering properties of aerosols

F.Figure C1 .
Figure C1.Comparison of top-of-ocean radiance (L) and degree of linear polarization (DoLP) computed by the extended addingdoubling model (solid lines in the left two panels) and successiveorders-of-scattering (labeled as "SOS", dots in the left two panels) with the bio-optical model described in Appendix D for an ocean system (ocean and air-water interface only, no atmosphere).Shown in the top right panel is the percentage difference of reflectance calculated by 100 × (L EAD − L SOS )/L EAD , where the subscript "EAD" denotes our extended adding-doubling method.Shown in the bottom right panel is the difference of DoLP computed by 100 × (DoLP EAD − DoLP SOS ).The chlorophyll concentration is [Chl_a] = 0.30 mg m −3 and the solar zenith angle is 60 • .The ocean surface is roughened by a wind of speed 7 m s −1 and ocean optical thickness is set to 10.An arbitrary combination of refractive index and slope parameter (n p = 1.05, γ p = 3.71) is chosen to compute the Fourier-Forland phase function.The results are plotted for viewing angles (θ v ) increasing from 0 to 87 • with an angular step of 3 • ; the five azimuthal planes (φ v ) are 0, 45, 90, 135 and 180 • (shown in black, red, blue, green, and cyan, respectively) with respect to the principal plane (namely O-XZ in Fig.1).

FFigure D1 .
Figure D1.Left panel: refractive index (n p ) of phytoplankton and their covariant particles as a function of backscattering efficiency (B bp ); Middle panel: slope parameter (γ p ) as a function of backscattering efficiency (B bp ); and right panel: refractive index (n p ) as a function of the slope parameter (γ p ). B bp is computed from Eq. (D10) as a function of Chl a concentration [Chl_a].The refractive index (n p ) and slope parameter (γ p ) characterizing a Junge size distribution are then determined by solving Eqs.(D9) and (D11) numerically.

Table 2 .
Parameters in ocean retrieval and Lagrange multipliers for smoothness constraints.

Table 3 .
Cases for truth-in/truth-out retrieval tests.
Xu et al.:Joint retrieval of aerosol and water-leaving radiance leaving radiance.As AOD and SSA errors are the largest for dust aerosols, the normalized water-leaving radiance retrieval error also becomes largest in the presence of dust.