# Joint retrieval of aerosol and water-leaving radiance from multispectral, multiangular and polarimetric measurements over ocean

Feng Xu, Oleg Dubovik, Peng-Wang Zhai, David J. Diner, Olga V. Kalashnikova, Felix C. Seidel, Pavel Litvinov, Andrii Bovchaliuk, Michael J. Garay, Gerard van Harten, et al.

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

## ► To cite this version:

Feng Xu, Oleg Dubovik, Peng-Wang Zhai, David J. Diner, Olga V. Kalashnikova, et al. Joint retrieval of aerosol and water-leaving radiance from multispectral, multiangular and polarimetric measurements over ocean. Atmospheric Measurement Techniques, 2016, 9, pp.2877-2907. 10.5194/amt-9-2877-2016. insu-03686204

HAL Id: insu-03686204

https://insu.hal.science/insu-03686204

Submitted on 3 Jun 2022

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

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

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Joint retrieval of aerosol and water-leaving radiance from multispectral, multiangular and polarimetric measurements over ocean

Feng Xu<fnref n="1" />, Oleg Dubovik<fnref n="2" />, Peng-Wang Zhai<fnref n="3" />, David J. Diner<fnref n="1" />, Olga V. Kalashnikova<fnref n="1" />, Felix C. Seidel<fnref n="1" />, Pavel Litvinov<fnref n="2" />, Andrii Bovchaliuk<fnref n="2" />, Michael J. Garay<fnref n="1" />, Gerard van Harten<fnref n="1" />, and Anthony B. Davis<fnref n="1" />

*Correspondence to:* Feng Xu (feng.xu@jpl.nasa.gov)

Received: 15 December 2015 – Published in Atmos. Meas. Tech. Discuss.: 18 January 2016

Revised: 27 May 2016 – Accepted: 2 June 2016 – Published: 8 July 2016

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

## 1 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.

<footnote n="1">¹Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, USA</footnote>

<footnote n="2">²Laboratoire d'Optique Atmosphérique, UMR8518, CNRS/Université Lille-1, Villeneuve d'Ascq, France</footnote>

<footnote n="3">³Department of Physics, University of Maryland, Baltimore County, Baltimore, Maryland, USA</footnote>
When aerosol is a major target of retrieval, water-leaving radiance is often empirically estimated or even neglected in operational algorithms employed by current-generation 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\text{nm}} = 0.005$, cf. Remer et al., 2005, 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, 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 ocean-atmosphere 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 adding-doubling method (Chowdhary et al., 2006) and successive-order-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.

## 2 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 air-water 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,
<figure><img src="image_4.png" /><figcaption>Figure 1. 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 air-water interface are completed by use of the adding strategy. The Sun illuminates the top-of-atmosphere with solar zenith angle $\theta_0$ and azimuthal plane $\phi_0$. We define $\phi = \phi_v - \phi_0$, where the sensor views the atmosphere at viewing angle $\theta_v$ and azimuthal angle $\phi_v$.</figcaption></figure>

and single scattering properties of aerosols in the mixed layer are demonstrated in Appendix A.

## 2.1 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., 2010, 2011, 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.

## 2.2 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 bio-optical 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, $\Delta 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

$$ \pi \mathbf{R}_{\text{surf}} = f_{\text{foam}} a_{\text{foam}} \mathbf{D}_0 + (1 - f_{\text{foam}}) \mathbf{R}_W + (1 - f_{\text{foam}}) \mathbf{R}_{\text{WL}}^{\text{Bio}} + (1 - f_{\text{foam}}) \Delta a_{\text{WL}} \mathbf{D}_0, \quad (1) $$

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 \times 10^{-6} \times W^{3.52}$ (Koepeke, 1984), and $R_{WL}^{\text{Bio}}$ is the reflection matrix of the ocean-interface system with Raman-scattering correction (see Appendix C). Note that $R_{WL}^{\text{Bio}}$ 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 $\Delta 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_{WL}^{Bio}$ and $\Delta a_{WL}$ in Eq. (1) the normalized water-leaving radiance by setting the Sun at zenith ($\theta_0 = 0^\circ$) and viewing angle to be nadir ($\theta_v = 0^\circ$):

$$nLw = \frac{F_0}{\pi} \left(\frac{d_0}{d}\right)^2 \left[ R_{WL,11}^{Bio}(\theta_v = 0^\circ; \theta_0 = 0^\circ; [\text{Chl}_a]) + \Delta a_{WL} \right], \quad (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_{WL}^{Bio}$, $a_{foam}$, $\Delta 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:

$$Q_1 = R_{\text{atmos}}^* R_{\text{surf}} \quad (3a)$$

$$Q_n = Q_1 Q_{n-1} \quad (3b)$$

$$S = \sum_{n=1}^{\infty} Q_n, \quad (3c)$$

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

$$D = T_{\text{atmos}} + S \exp\left(-\frac{\tau_{\text{atmos}}}{\mu_0}\right) + S T_{\text{atmos}} \quad (3d)$$

$$U = R_{\text{surf}} \exp\left(-\frac{\tau_{\text{atmos}}}{\mu_0}\right) + R_{\text{surf}} D, \quad (3e)$$

where $\mu_0 = \cos\theta_0$. The reflection matrix of the full CAOS is

$$R_{\text{CAOS}} = R_{\text{atmos}} \exp\left(-\frac{\tau_{\text{atmos}}}{\mu}\right) U + T_{\text{atmos}}^* U, \quad (3f)$$

where $\mu = |\cos\theta|$. 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:

$$\text{BRF}_{\text{tot}} = \pi \sum_{m=0}^{\infty} (2 - \delta_{0m}) R_{\text{CAOS},11}^{(m)} \cos m\phi \quad (4a)$$

$$\text{qBRF}_{\text{tot}} = \pi \sum_{m=0}^{\infty} (2 - \delta_{0m}) R_{\text{CAOS},21}^{(m)} \cos m\phi \quad (4b)$$

$$\text{uBRF}_{\text{tot}} = \pi \sum_{m=0}^{\infty} (2 - \delta_{0m}) R_{\text{CAOS},31}^{(m)} \cos m\phi \quad (4c)$$

$$\text{vBRF}_{\text{tot}} = \pi \sum_{m=0}^{\infty} (2 - \delta_{0m}) R_{\text{CAOS},41}^{(m)} \cos m\phi, \quad (4d)$$

where the bidirectional reflectance factor $BRF_{tot}$ and $DoLP = \frac{\sqrt{qBRF_{tot}^2 + uBRF_{tot}^2 + vBRF_{tot}^2}}{BRF_{tot}^2}$ are used to fit the observation. Since the sunlight is unpolarized, other matrix entries (namely $R_{CAOS,ij}$, with $j \ge 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_{\text{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}(\lambda) = BRF_{tot}(\lambda) \exp[-\tau_{ozone}(\lambda)(1/\cos\theta_0 + f_{ozone}/\cos\theta_v)]$, where $\tau_{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_{\text{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-
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 - \Delta x_k$ with $x_k$ being the solution after $k$ iterations and $\Delta x_k$ being the increment being obtained by $\Delta x_k = (J_k^T)^{-1}\Delta y_k$, where $J_k$ is the Jacobian matrix evaluated with $x_k$, and $\Delta y_k$ is the difference between model and measurement ($\Delta y_k = y(x_k) - y_{\text{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[J_k^T(C_f)^{-1}J_k + \gamma_{k,1}W_{k,1} + \gamma_{k,2}W_{k,2} + \dots] > 0$ and $\Delta x_k = [J_k^T(C_f)^{-1}J_k + \gamma_{k,1}W_{k,1} + \gamma_{k,2}W_{k,2} + \dots]^{-1}\Delta y'_k$, where $C_f$ is the covariance matrix of the measured signals, $W_{k,i}$ denotes the imposed various constraints, $\gamma_k$ is a Lagrange multiplier that assigns a weight to the constraint, and $\Delta y'_k$ incorporates $\Delta 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 $\gamma_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, 1963, 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. (2008, 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. (2008, 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

<table><caption>**Table 1.** Median radius ($r_m$) and standard deviation ($\sigma$) of $N_{sc} = 5$ volume weighted log-normal size components, namely $dv_i(r)/d\ln r$ in Eqs. (A4)–(A5).</caption><thead><tr><th>Component number</th><th>Median radius ($r_m$, µm)</th><th>Standard deviation ($\sigma$)</th></tr></thead><tbody><tr><td>1</td><td>0.1</td><td>0.35</td></tr><tr><td>2</td><td>0.1732</td><td>0.35</td></tr><tr><td>3</td><td>0.3</td><td>0.35</td></tr><tr><td>4</td><td>1</td><td>0.5</td></tr><tr><td>5</td><td>2.9</td><td>1</td></tr></tbody></table>

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 five-component-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.

## 3.1 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).

$$
\begin{align}
\mathbf{C}(\mathbf{x}) &= \sum_{i=1}^{N} \Psi(\mathbf{x}_i) + \frac{1}{2} \mathbf{x}^\mathrm{T} \mathbf{\Omega}_{\mathrm{interpatch}} \mathbf{x} \tag{5} \\
&= \sum_{i=1}^{N} [\Psi_{\mathrm{f}}(\mathbf{x}_i) + \Psi_{\mathrm{s}}(\mathbf{x}_i) + \Psi_{\mathrm{a}}(\mathbf{x}_i)] + \frac{1}{2} \mathbf{x}^\mathrm{T} \mathbf{\Omega}_{\mathrm{interpatch}} \mathbf{x} \nonumber \\
&= \frac{1}{2} \sum_{i=1}^{N} \left[ \Delta \mathbf{y}_i^\mathrm{T} \mathbf{W}_{\mathrm{f},i}^{-1} \Delta \mathbf{y}_i + \gamma_s \mathbf{x}_i^\mathrm{T} \mathbf{\Omega}_{\mathrm{s},i} \mathbf{x}_i \right. \nonumber \\
&\qquad \left. + \gamma_a (\mathbf{x}_i - \mathbf{x}_i^*)^\mathrm{T} \mathbf{W}_{\mathrm{a},i}^{-1} (\mathbf{x}_i - \mathbf{x}_i^*) \right] + \frac{1}{2} \mathbf{x}^\mathrm{T} \mathbf{\Omega}_{\mathrm{interpatch}} \mathbf{x}, \nonumber
\end{align}
$$

where $\mathbf{x}_i$ is an iterative solution for the set of parameters being retrieved and $\mathbf{x}_i^*$ is an a priori estimate of the solution corresponding to the $i$th patch, $\mathbf{x} = [\mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3, \dots, \mathbf{x}_N]$; $\Psi_{\mathrm{f}}(\mathbf{x}_i)$, $\Psi_{\mathrm{s}}(\mathbf{x}_i)$ and $\Psi_{\mathrm{a}}(\mathbf{x}_i)$ correspond to the residues of fitting observations, the spectral smoothness constraints, and the a priori estimate, respectively; $\mathbf{\Omega}_{\mathrm{s},i}$ is a smoothness matrix for constraining the spectral variation of aerosol optical properties and water-leaving radiances across the relevant bands; $\mathbf{W}_{\mathrm{f}}$ and $\mathbf{W}_{\mathrm{a}}$ are the weighting matrices for measurements and the a priori estimate, respectively; $\gamma$ denotes the relevant Lagrange multipliers; $\Delta\mathbf{y}_i$ is the difference between the model and measurements for the $i$th patch [$\Delta\mathbf{y}_i = \mathbf{y}(\mathbf{x}_i) - \mathbf{y}_{\mathrm{meas}}$]; and $\mathbf{\Omega}_{\mathrm{interpatch}}$ is the interpatch smoothness matrix constructed for the patches along two orthogonal directions ($u$ and $v$) of the image, namely

$$
\mathbf{\Omega}_{\text{interpatch}} = \gamma_u \mathbf{S}^{(m_u)} \cdot \mathbf{T} \mathbf{S}^{(m_u)} + \gamma_v \mathbf{S}^{(m_v)} \cdot \mathbf{T} \mathbf{S}^{(m_v)}, \quad (6)
$$

where the derivative matrix $\mathbf{S}^{(m)}$ is constructed from the $m$th order difference, and $\gamma_u$ and $\gamma_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 $\mathbf{x}_{i,k+1}$ containing parameters of aerosol and surface properties for the $i$th patch is updated as

$$
\mathbf{x}_{i,k+1} = \mathbf{x}_{i,k} - t_p \Delta \mathbf{x}_{i,k}, \quad (7)
$$

where the multiplier $t_p (0 \le t_p \le 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 $k$th iteration gives

the increment of solution for each patch ($\Delta\mathbf{x}_{i,k}$,

$$
\begin{align}
& \left[ \begin{pmatrix} \mathbf{A}_{1,k} & 0 & \dots & 0 \\ 0 & \mathbf{A}_{2,k} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \mathbf{A}_{N,k} \end{pmatrix} + \mathbf{\Omega}_{\text{interpatch},k} \right] \begin{pmatrix} \Delta \mathbf{x}_{1,k} \\ \Delta \mathbf{x}_{2,k} \\ \vdots \\ \Delta \mathbf{x}_{N,k} \end{pmatrix} \nonumber \\
& = \left[ \begin{pmatrix} \nabla \Psi(\mathbf{x}_{1,k}) \\ \nabla \Psi(\mathbf{x}_{2,k}) \\ \vdots \\ \nabla \Psi(\mathbf{x}_{N,k}) \end{pmatrix} + \mathbf{\Omega}_{\text{interpatch},k} \begin{pmatrix} \mathbf{x}_{1,k} \\ \mathbf{x}_{2,k} \\ \vdots \\ \mathbf{x}_{N,k} \end{pmatrix} \right], \tag{8}
\end{align}
$$

where the Fisher matrix for the $i$th patch is a function of Jacobian matrix $\mathbf{J}_{i,k}$ and weighting matrix $\mathbf{W}_{f,i}$,

$$
\mathbf{A}_{i,k} = \mathbf{J}_{i,k}^{\mathrm{T}} \mathbf{W}_{f,i}^{-1} \mathbf{J}_{i,k} + \gamma_{s,i,k} \mathbf{\Omega}_{\Delta,i} + \gamma_{a,i,k} \mathbf{W}_{a,i}^{-1}, \quad (9)
$$

and $\nabla\Psi(\mathbf{x}_{i,k})$ is the gradient of the minimized quadratic form:

$$
\nabla \Psi(\mathbf{x}_{i,k}) = \mathbf{J}_{i,k}^{\mathrm{T}} \mathbf{W}_{f,i}^{-1} (\mathbf{y}_{i,k} - \mathbf{y}_{i,\mathrm{meas}}) + \gamma_{s,i,k} \mathbf{\Omega}_{\mathrm{s},i} \mathbf{x}_{i,k} \quad (10) \\
\qquad + \gamma_{a,i,k} \mathbf{W}_{a,i}^{-1} (\mathbf{x}_{i,k} - \mathbf{x}_{i}^{*}),
$$

where $\mathbf{y}_{\mathrm{meas}}$ contains the measurement data, $\mathbf{y}_k$ contains the modeled radiance and polarization with $\mathbf{x}_k$, $\mathbf{W}_f$ is the weighting matrix defined as the covariance matrix $\mathbf{C}_f$ normalized by its first diagonal element namely $\mathbf{W}_f = (1/\sigma_{\mathrm{SD},1}^2)\mathbf{C}$ (with $\sigma_{\mathrm{SD}}$ being the standard deviation), $\mathbf{W}_a$ is the weighting matrix of the a priori estimate $\mathbf{x}^*$, and $\mathbf{\Omega}_s$ is the single-patch-based smoothness matrix containing sub-smoothness matrices for all parameters. The Lagrange multipliers $\gamma_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 ($\sigma_a$) of aerosol layer, nonspherical particle fraction ($f_{ns}$), wind speed over ocean ($W$), Chl $a$ concentration ([Chl_a]), and $\Delta a_{WL}$, which adjust the nLw values in the second step of the retrieval. These parameters form the solution vector $\mathbf{x} = \log[m_r(\lambda), m_i(\lambda), C_v(r_m), h_a, \sigma_a, f_{ns}, W, \text{Chl}_a, a_{WL, \text{Const}}(\lambda) + \Delta a_{WL}(\lambda)]^\mathrm{T}$, 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, \text{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 $\gamma_s \mathbf{\Omega}_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 $\Delta a_\lambda$, namely for all patches:

$$
\gamma_s \mathbf{\Omega}_s = \text{diag} \left\{ \begin{matrix} \mathbf{0}, \mathbf{0}, \mathbf{0}, \mathbf{0}, \mathbf{0}, \gamma_{s,m_r} \mathbf{\Omega}_{s,m_r}, \gamma_{s,m_i} \mathbf{\Omega}_{s,m_i}, \mathbf{0}, \mathbf{0}, \end{matrix} \begin{matrix} \mathbf{0}, \mathbf{0}, \mathbf{0}, \mathbf{0}, \gamma_{s,(a_{WL,Const}+\Delta a_{WL})} \mathbf{\Omega}_{s,(a_{WL,Const}+\Delta a_{WL})} \end{matrix} \right\}, \quad (11)
$$

where $\mathbf{0}$ represents a zero submatrix for a parameter not being subject to any smoothness constraints; and the Lagrange multipliers $\gamma_s$ are predetermined and given in Table 2.
<table><caption>**Table 2.** Parameters in ocean retrieval and Lagrange multipliers for smoothness constraints.</caption><thead><tr><th></th><th>Range</th><th>Order of<br>finite<br>difference for<br>spectral<br>smoothness<br>constraints<br>(m<sub>s</sub>)</th><th>Lagrange<br>regularization<br>factor (γ<sub>s</sub>)</th><th>Order of finite<br>difference for<br>interpatch<br>smoothness<br>constraints<br>(m<sub>(u,v)</sub>)</th><th>Lagrange<br>regularization<br>factor γ<sub>(u,v)</sub></th></tr></thead><tbody><tr><td rowspan="6">Aerosol parameters</td><td>Volume concentration of size<br>components<br>(C<sub>v</sub>(r<sub>m</sub>), µm<sup>3</sup> µm<sup>-2</sup>)</td><td>[1.0 × 10<sup>-6</sup>, 5]</td><td>—</td><td>1</td><td>1</td></tr><tr><td>Mean height of aerosol<br>distribution profile (h<sub>a</sub>, km)</td><td>[0.05, 10]</td><td>—</td><td>1</td><td>0.01</td></tr><tr><td>Width of aerosol distribution<br>profile (σ<sub>a</sub>)</td><td>[0.5, 2.5]</td><td>—</td><td>1</td><td>0.01</td></tr><tr><td>Refr. index (real part: n<sub>r</sub>(λ))</td><td>[1.33, 1.60]</td><td>1</td><td>0.1</td><td>10</td></tr><tr><td>Refr. index (imag. part: n<sub>i</sub>(λ))</td><td>[5 × 10<sup>-7</sup>, 5 × 10<sup>-1</sup>]</td><td>2</td><td>0.01</td><td>1</td></tr><tr><td>Nonspherical particle fraction<br>(f<sub>ns</sub>)<sup>a</sup></td><td>[1 × 10<sup>-3</sup>, 1]</td><td>—</td><td>—</td><td>0.1</td></tr><tr><td rowspan="6">Surface parameters</td><td>Adjustment term<br>(Δa<sub>WL</sub>(λ),<br>mW cm<sup>-2</sup> sr<sup>-1</sup> µm<sup>-1</sup>)</td><td>π/F<sub>0</sub> × (d/d<sub>0</sub>)<sup>2</sup> × nLw<sub>1</sub> ×<br>[-15%, +15%]<sup>b</sup></td><td>3<sup>c</sup></td><td>0.1<sup>c</sup></td><td>0.1<sup>c</sup></td></tr><tr><td>Chlorophyll a concentration<br>([Chl_a], mg m<sup>-3</sup>)</td><td>Step 1 for [Chl_a]<sub>1</sub>:<br>[0.02, 15]<br>Step 2 for [Chl_a]<sub>2</sub>:<br>[0.85, 1.15] × [Chl_a]<sub>1</sub></td><td>—</td><td>—</td><td>1</td></tr><tr><td>Normalized water-leaving<br>radiance<sup>d</sup> (nLw (λ),<br>mW cm<sup>-2</sup> sr<sup>-1</sup> µm<sup>-1</sup>)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Surface wind speed (W, m s<sup>-1</sup>)</td><td>[1, 30]</td><td>—</td><td>—</td><td>1</td></tr></tbody></table>

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

$$ \nabla \Psi(\mathbf{x}_{i,k}) = \mathbf{J}_{i,k}^{\mathrm{T}} \mathbf{W}_{f,i}^{-1} (\mathbf{y}_{i,k} - \mathbf{y}_{i,\mathrm{meas}}) + \gamma_{s,i,k} \mathbf{\Omega}_{s,i} \mathbf{x}_{i,k}. \quad (12) $$

When the spectral and spatial smoothness constraints are turned off (namely setting $\gamma_s = \gamma_u = \gamma_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

$$ \begin{aligned} 2 \sum_{i=1}^{N_{\mathrm{patch}}} \Psi(\mathbf{x}_{k,i}) + \mathbf{x}_k \mathbf{\Omega}_{\mathrm{interpatch},k} \mathbf{x}_k^{\mathrm{T}} &\le N_{\mathrm{interpatch}} \varepsilon_f^2 \\ + \sum_{i=1}^{N_{\mathrm{patch}}} (N_{\mathrm{f},i} + N_{\mathrm{s},i} + N_{\mathrm{a},i}^* - N_{\mathrm{a},i}) \varepsilon_f^2, \end{aligned} \quad (13) $$

where $N_{\mathrm{f},i}$, $N_{\mathrm{s},i}$, $N_{\mathrm{a},i}$ and $N_{\mathrm{a},i}^*$ are the number of observations, spectral smoothness, number of unknowns, and a priori estimates of parameters corresponding to $i$th patch, respectively; $N_{\mathrm{interpatch}}$ is the number of spatial smoothness constraints; and $\varepsilon_f^2$ 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

<footnote><sup>a</sup> Nonspherical (spheroidal) particle fraction is excluded in truth-in/truth-out test but included in real data retrieval; <sup>b</sup> The subscript “1” of nLw means the normalized water-leaving radiance determined from Chl a concentration ([Chl a]<sub>1</sub>) retrieved at step 1; <sup>c</sup> Determined with the consideration of constant offset π/F<sub>0</sub> × (d/d<sub>0</sub>)<sup>2</sup> × nLw<sub>1</sub>(λ); <sup>d</sup> The normalized water-leaving radiance is not directly retrieved. After the second step retrieval, the updated Chl a concentration [Chl a]<sub>2</sub> and the adjustment term Δa<sub>WL</sub> are used to derive it via Eq. (2).</footnote>
threshold value, $\epsilon_c^2$. Namely,

$$ \epsilon_c^2 \geq \frac{\left[ 2 \sum_{i=1}^{N_{\text{patch}}} \Psi(\mathbf{x}_{k+1,i}) + \mathbf{x}_{k+1} \mathbf{\Omega}_{\text{interpatch},k+1} \mathbf{x}_k^T \right] - \left[ 2 \sum_{i=1}^{N_{\text{patch}}} \Psi(\mathbf{x}_{k,i}) + \mathbf{x}_k \mathbf{\Omega}_{\text{interpatch},k} \mathbf{x}_k^T \right]}{2 \sum_{i=1}^{N_{\text{nxel}}} \Psi(\mathbf{x}_{k,i}) + \mathbf{x}_k \mathbf{\Omega}_{\text{interpatch},k} \mathbf{x}_k^T}, \quad (14) $$

is the second criterion to terminate the optimization.

## 3.2 Determination of Lagrange multipliers

Following Dubovik and King (2000), the Lagrange multipliers reflecting the strength of smoothness constraints are defined as

$$ \gamma_g = \epsilon_f^2 / \epsilon_g^2 \quad \text{and} \quad \gamma_a = \epsilon_f^2 / \epsilon_a^2, \quad (15) $$

where $\epsilon_f^2$, $\epsilon_a^2$ and $\epsilon_g^2$ 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 $\epsilon_g^2$ for a given parameter to be retrieved ($x_j$), which is a function of $t$, the most unsmooth-known solution $x_j^{\text{ns}}(t)$ over the target area is used, namely,

$$ \epsilon_g^2 = \int_{t_{\text{min}}}^{t_{\text{max}}} \left( \frac{d^m[x_j^{\text{us}}(t)]}{d^m t} \right)^2 dt, \quad (16) $$

where $t_{\text{min}}$ and $t_{\text{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:

$$ \gamma_g^{\text{Final}} = \frac{N_f \tilde{\epsilon}_f^2}{N_g \epsilon_f^2} \gamma_g \quad \text{and} \quad \gamma_a^{\text{Final}} = \frac{N_f \tilde{\epsilon}_f^2}{N_a \epsilon_f^2} \gamma_a. \quad (17) $$

There are two differences between $\gamma_{...}^{\text{Final}}$ and $\gamma_{...}$:

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 $\epsilon_{...}^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 $\tilde{\epsilon}_f^2/\epsilon_f^2$ is introduced with $\tilde{\epsilon}_f^2$ estimated as the dynamic fitting residual during iterations:

$$ \tilde{\epsilon}_f^2(\mathbf{x}_k) \approx \frac{2 \sum_{i=1}^{N_{\text{patch}}} \Psi(\mathbf{x}_{k,i}) + \mathbf{x}_k \mathbf{\Omega}_{\text{interpatch},k} \mathbf{x}_k^T}{N_{\text{interpatch}} + \sum_{i=1}^{N_{\text{patch}}} (N_{f,i} + N_{s,i} + N_{a^*,i} - N_{a,i})}. \quad (18) $$

With the multiplier $\tilde{\epsilon}_f^2/\epsilon_f^2$, 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 $\tilde{\epsilon}_f^2$ becomes equal to $\epsilon_f^2$. As a result of this adjustment of the Lagrange multiplier, the nonlinear iteration becomes significantly more monotonic.

## 3.3 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 $\Delta 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 bio-optical 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 $\Delta a_{WL}$ term be optimized simultaneously ($\Delta a_{\lambda, 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<sub>adjust</sub>/nLw<sub>1</sub> ≤ +15%, with nLw<sub>1</sub> being the nLw from the first retrieval step); and (2) imposes a spectral smoothness constraint on nLw(λ).

# 4 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.

## 4.1 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 water-vapor 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⁻³. A perturbation of ±10% was imposed on the water-leaving radiance predicted by the Chl *a*-based 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⁻¹. 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, \text{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 and Landgraf, 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.

### 4.1.1 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,
<table><caption>**Table 3.** Cases for truth-in/truth-out retrieval tests.</caption><thead><tr><th></th><th>Weakly absorbing aerosol</th><th>Moderately absorbing aerosol</th><th>Dust aerosol</th></tr></thead><tbody><tr><td>Targeted AOT at 555 nm</td><td></td><td>0.02, 0.1, 0.3, 0.5, 1.0</td><td></td></tr><tr><td>Volume fractions (f<sub>v,1-5</sub>)</td><td>4, 32, 20, 4, 40%</td><td>16, 56, 6, 6, 16%</td><td>2, 8, 1, 24, 65%</td></tr><tr><td>Aerosol</td><td></td><td></td><td></td></tr><tr><td>Mean height of aerosol distribution profile (h<sub>a</sub>, km)</td><td></td><td>1</td><td></td></tr><tr><td>Half width of aerosol distribution profile (σ<sub>a</sub>, km)</td><td></td><td>0.75</td><td></td></tr><tr><td>Refractive index (mean of real part n<sub>r</sub>(λ))</td><td>1.388<sup>a</sup></td><td>1.522<sup>a</sup></td><td>1.497</td></tr><tr><td>Refractive index (mean of imag. part: n<sub>i</sub>(λ))</td><td>1.98 × 10<sup>-3</sup><sup>a</sup></td><td>1.32 × 10<sup>-2</sup><sup>a</sup></td><td><sup>b</sup></td></tr><tr><td>Surface</td><td></td><td></td><td></td></tr><tr><td>Chlorophyll a ([Chl_a], mg m<sup>-3</sup>)</td><td></td><td>0.05, 0.2, 1.0</td><td></td></tr><tr><td>Adjustment term (Δa<sub>WL</sub>(λ), mW cm<sup>-2</sup> sr<sup>-1</sup> µm<sup>-1</sup>)</td><td></td><td>corresponding to ±10% perturbation on bio-optical model simulated nLw at AirMSPI 355, 385, 445, 475, and 550, 660 and 865 nm spectral bands</td><td></td></tr><tr><td>Surface wind speed (W, m s<sup>-1</sup>)</td><td></td><td>4</td><td></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 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 θ<sub>0</sub> = 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.</figcaption></figure>

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<sup>-3</sup>, 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 (τ<sub>555</sub> ≤ 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-

<footnote><sup>a</sup> slightly dependent on wavelength but average values are listed here; <sup>b</sup> 8.04 × 10<sup>-3</sup> (355 nm), 7.74 × 10<sup>-3</sup> (380 nm), 4.98 × 10<sup>-3</sup> (445 nm), 4.10 × 10<sup>-3</sup> (470 nm), 2.01 × 10<sup>-3</sup> (555 nm), 4.28 × 10<sup>-4</sup> (660 nm), and 3.27 × 10<sup>-4</sup> (865 nm).</footnote>
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.

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<sub>retrieved</sub> – AOD<sub>true</sub>|/AOD<sub>true</sub>) 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<sub>retrieved</sub> – SSA<sub>true</sub>|) 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:

$$r_{\text{eff,fine}} = \left[ \int_{r_{\text{min}}}^{r_{\text{cri}}} \frac{dv(r)}{d\ln r} d\ln r \right] \left[ \int_{r_{\text{min}}}^{r_{\text{cri}}} \frac{1}{r} \frac{dv(r)}{d\ln r} d\ln r \right]^{-1} \quad (19)$$

$$r_{\text{eff,coarse}} = \left[ \int_{r_{\text{cri}}}^{r_{\text{max}}} \frac{dv(r)}{d\ln r} d\ln r \right] \left[ \int_{r_{\text{cri}}}^{r_{\text{max}}} \frac{1}{r} \frac{dv(r)}{d\ln r} d\ln r \right]^{-1}, \quad (20)$$

where the lower size limit $r_{\text{min}} = 0.04$ µm and the upper size limit $r_{\text{max}} = 15$ µm. Setting $r_{\text{cri}}$ to be 0.75 µm for weakly and moderately absorbing aerosols and $r_{\text{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_{\text{eff, retrieved}}$ – $r_{\text{eff, true}}$|/$r_{\text{eff, true}}$) is also observed for all types of aerosols, as shown in Fig. 7d. For $\tau_{555} \ge 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 $\omega_{0,355\text{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<sup>-3</sup> (with ±10% perturbation imposed on the water-leaving radiance), the retrieval error of normalized water-leaving radiance ($\Delta nLw = nLw_{\text{retrieved}} – nLw_{\text{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/\pi$ (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 $\tau_{555} > \sim 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 $\tau_{555} > \sim 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 ($\tau_{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).

### 4.1.2 Effect of multipatch vs. single-patch retrieval

Taking the case of median loading ($\tau_{555} = 0.3$) of weakly absorbing aerosols and median Chl a concentration [Chl_a] = 0.2 mg m<sup>-3</sup> as an example, Fig. 8a compares simulated single-patch- and multipatch-based retrievals of AOD, SSA, PSD, and Chl a concentration. The Sun illumination and AirMSPI viewing geometry at the USC_SeaPRISM AERONET-OC site on 6 February 2013
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Simulated true and retrieved spectral AOD for different scene conditions. Left column of three panels: weakly absorbing aerosol. Middle column of three panels: moderately absorbing aerosol. Right column of three panels: dust aerosol. AOD is retrieved for three values of Chl *a* concentration: 0.05 (top row of panels), 0.2 (middle row of panels), and 1.0 mg m⁻³ (bottom row of panels), with ±10% perturbation of water-leaving radiance. Five aerosol loadings, corresponding to τ₅₅₅ = 0.02, 0.10, 0.30, 0.50, and 1.0, are plotted in dark blue, dark red, dark yellow, purple, and green, respectively. The lines with crosses at the AirMSPI wavelengths represent the true AODs, while the open circles correspond to the retrieved values. The synthetic data are from one of the simulated scenarios of AirMSPI observation over USC_SeaPRISM AERONET-OC site (θ₀ = 25°, φ ≈ 95°). Though not plotted, the spatial variation of the retrieved AOD across the whole image is less than 1% for all spectral bands.</figcaption></figure>

is used. While the single patch-based retrieval leads to spa-
tially highly variable Chl *a* concentrations with a mean value
of 0.26 mg m⁻³ (namely 30% retrieval error), the multi-
patch algorithm yields a more stable and accurate value of
0.21 mg m⁻³, which is within 5% of the true value. Cor-
respondingly, the accuracy of the nLw retrieval improves
by 0.04, 0.03, and 0.01 mW cm⁻² sr⁻¹ µm⁻¹ at 445, 470,

and 555 nm, respectively, which is a nonnegligible amount
compared to the PACE-tolerated uncertainty 0.07, 0.06, and
0.05 mW cm⁻² sr⁻¹ µm⁻¹ at these bands; the AOD accuracy
at 355, 555, 865 nm improves by 3.4, 6.0, and 6.4%, respec-
tively; and the SSA accuracy improves by 0.008, 0.013, and
0.019. For the single patch-based approach, combinations
of aerosol type, amount, and nLw that fit the simulated ob-
<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

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.

### 4.1.3 Comparison to direct water-leaving radiance retrieval

For the same scene, parameters used to compare the single- and 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-
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

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⁻² sr⁻¹ µm⁻¹ 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⁻² sr⁻¹ µm⁻¹ in these bands, the accuracy gain from using the bio-optical model is significant.

### 4.1.4 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
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Panel layout as in Fig. 3 but 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.</figcaption></figure>

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 $\Delta$nLw($\lambda$)/nLw(555); in other words, the systematic errors mainly propagate to the overall magnitude of nLw($\lambda$) curve while the relative spectral shape is affected to a much lesser degree.

## 4.2 Retrievals with real AirMSPI observations

Following algorithm validation using the truth-in/truth-out tests, we applied the algorithm to actual AirMSPI ob-
<figure><img src="image_10.png" /><figcaption>Figure 7.</figcaption></figure>

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 AirMSPI and AERONET AOD and SSA retrievals are within the

<figure><img src="image_11.png" /><figcaption>Figure 7.</figcaption></figure>

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. (2011, 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,
<figure><img src="image_12.png" /><figcaption>Figure 7.</figcaption></figure>

the differences are found to be 0.0396, 0.0118, 0.0198, and 0.0077 mW cm⁻² sr⁻¹ µm⁻¹ 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⁻² sr⁻¹ µm⁻¹ 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 AirMSPI 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} \sim 0.04$). Independent

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

**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⁻³, 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⁻³ are displayed in (a)–(d). Via truth-in/truth-out tests, the uncertainties are estimated for AirMSPI multispectral, multiangular, and multipolar-metric 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.

surface measurements to validate the nLw retrieval were not available at this site.

## 5 Summary and outlook

Accurate retrieval of both aerosol properties and water-leaving 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 high-quality 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
<figure><img src="image_14.png" /></figure>

**Figure 8.** Comparison of single-patch- and multipatch-based retrievals of (a) AOD, SSA, aerosol size distribution, and Chl *a* concentration for the median AOD ($\tau_{555} = 0.3$) of weakly absorbing aerosols and Chl *a* concentration of $0.2 \text{ mg m}^{-3}$. The simulation uses the Sun and viewing geometry corresponding to the AirMSPI overflight of the USCPRISM AERONET site. Image-averaged Chl *a* concentrations are $0.29$ and $0.22 \text{ 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.

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. An efficient RT modeling strategy has been developed that couples separate runs for modeling RT in two Rayleigh layers, an aerosol/Rayleigh mixed layer, and an ocean medium. Repeated, time-consuming RT computations for layers whose properties are not perturbed during Jaco-

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

**Figure 9.** AOD, SSA, aerosol size distribution, and nLw retrieved using the bio-optical model compared to retrievals in which water-leaving radiance is modeled simply as Lambertian with arbitrary albedo.

<figure><img src="image_16.png" /><figcaption>**Figure 10.**</figcaption></figure>

bian evaluations are avoided. The Markov chain method is used for modeling RT in the mixed layer and the doubling method is used to model RT in the pure Rayleigh layer and
<figure><img src="image_17.png" /><figcaption>Figure 10.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 10.</figcaption></figure>

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

<figure><img src="image_19.png" /><figcaption>**Figure 10.** Similar to Fig. 7a–e and h for [Chl *a*] = 0.2 mg m⁻³ but with additional systematic error of +4% (open squares) and -4% (closed squares) included in the truth-in/truth-out retrieval tests.</figcaption></figure>

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⁻³, 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-
<figure><img src="image_20.png" /><figcaption>**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.</figcaption></figure>

<caption>**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.</caption>

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 our 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.
## Appendix A: Parameterizations of distribution profile, size, and single-scattering properties of aerosols

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 $h_a$ and standard deviation $\sigma_a$, characterizing the width of the aerosol layer. Then, the aerosol concentration profile $c_a$ is

$$c_a(h) = F_{\text{norm}}^{-1} \exp \left[ - \frac{(h - h_a)^2}{\sigma_a^2} \right], \quad \text{(A1)}$$

where the normalization factor $F_{norm}$ is used to ensure that $\int_{h_{min}}^{h_{max}} c_k(h)dh = 1$ and evaluates to

$$F_{\text{norm}} = \frac{\sqrt{\pi \sigma_a}}{2} \left[ \text{erf} \left( \frac{h_{\text{max}} - h_a}{\sigma_a} \right) - \text{erf} \left( \frac{h_{\text{min}} - h_a}{\sigma_a} \right) \right], \quad \text{(A2)}$$

where $\text{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 ($\tau_a$) is the sum of all size components:

$$\tau_a = \sum_{i=1}^{N_{\text{sc}}} C_{v,i} K_{\text{ext,a},i} = C_{v, \text{tot}} \sum_{i=1}^{N_{\text{sc}}} f_i K_{\text{ext,a},i}, \quad \text{(A3)}$$

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 ($f_i = C_{v,i}/C_{v,tot}$).

Moreover, the total aerosol size distribution is constituted as

$$\frac{dV(r)}{d\ln r} = \sum_{i=1}^{N_{\text{sc}}} \frac{dV_i(r)}{d\ln r} = \sum_{i=1}^{N_{\text{sc}}} C_{v,i} \frac{dv_i(r)}{d\ln r}. \quad \text{(A4)}$$

The associated normalized size distribution is

$$\frac{dv(r)}{d\ln r} = \sum_{i=1}^{N_{\text{sc}}} f_i \frac{dv_i(r)}{d\ln r}. \quad \text{(A5)}$$

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 $\sigma_i$, namely,

$$\frac{dv_i(r)}{d\ln r} = \frac{1}{\sqrt{2\pi}\sigma_i} \exp \left[ - \frac{(\ln r - \ln r_{m,i})^2}{2\sigma_i} \right]. \quad \text{(A6)}$$

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 ($\Delta\tau^{(n)}$), single-scattering albedo (SSA, $\omega_0^{(n)}$) and phase matrix ($\mathbf{P}^{(n)}$) of the nth sublayer are contributed by aerosol and Rayleigh molecules only, therefore

$$\Delta\tau^{(n)} = \Delta\tau_a^{(n)} + \Delta\tau_R^{(n)} \quad \text{(A7)}$$

$$\omega_0^{(n)} = \frac{\Delta\tau_R^{(n)} + \omega_{0,a}^{(n)} \Delta\tau_a^{(n)}}{\Delta\tau_R^{(n)} + \Delta\tau_a^{(n)}} \quad \text{(A8)}$$

$$\mathbf{P}^{(n)}(\Theta) = \frac{\Delta\tau_R^{(n)} \mathbf{P}_R^{(n)}(\Theta) + \omega_{0,a}^{(n)} \Delta\tau_a^{(n)} \mathbf{P}_a^{(n)}(\Theta)}{\Delta\tau_R^{(n)} + \omega_{0,a}^{(n)} \Delta\tau_a^{(n)}}, \quad \text{(A9)}$$

where $\mathbf{P}_R$ and $\mathbf{P}_a$ are the Rayleigh and aerosol phase matrix, respectively; the SSA of aerosol $\omega_{0,a}$ is a function of scattering coefficient ($K_{sca,a}$) and extinction coefficient ($K_{ext,a}$): $\omega_{0,a} = K_{sca,a}/K_{ext,a}$; $\Delta\tau_a^{(n)}$ is the AOD in the nth sublayer and can be evaluated analytically after considering the aerosol distribution profile (Eq. A1) according to the following:

$$\Delta\tau_a^{(n)} = \tau_a \left[ \text{erf}\left(\frac{h^{(n+1)} - h_a}{\sigma_a}\right) - \text{erf}\left(\frac{h^{(n)} - h_a}{\sigma_a}\right) \right] \quad \text{(A10)} \\ \left[ \text{erf}\left(\frac{h_{\text{max}} - h_a}{\sigma_a}\right) - \text{erf}\left(\frac{h_{\text{min}} - h_a}{\sigma_a}\right) \right]^{-1}.$$

$\Delta\tau_R^{(n)}$ 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 $\mathbf{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$) 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

$$\mathbf{X}_{\text{a,ext/sca}} = \sum_{i=1}^{N_{\text{sc}}} f_i \mathbf{X}_{\text{a,ext/sca},i}, \quad \text{(A11)}$$

where $\mathbf{X}$ represents any Mie property of $\{\mathbf{P}_a, K_{ext,a}$ and $K_{sca,a}\}$. Via Eqs. (A7)–(A9), $\mathbf{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.

## Appendix B: Modeling radiative transfer in atmosphere system

### B1 Markov chain method for RT in aerosol/Rayleigh mixed layer

The light propagation direction in the mixed layer is discretized into a finite number of angles over the range $0 \le \mu \le 1$, where $\mu = |u| = |\cos\theta|$, and $\theta$ is the angle of propagation relative to the downward normal. Within the framework of the Markov chain method, the probability of a photon to transition from one state $(n, u_i)$ to another $(n', u_j)$ is given by the transition matrices $T_{\text{Refl}}$ and $T_{\text{Trans}}$ for difusely reflected and transmitted light, respectively. The transition probability from state $(n', u_j)$ to emergence from the top and bottom of the mixed layer in direction $u_e$ is given by the extinction matrices $E_{\text{Refl}}$ and $E_{\text{Trans}}$, respectively. Given the initial distribution of photons in all states ($\Pi_0$) from the single-scattering computations, the multiple-scattering (indicated by subscript “M”) contributions to the reflection and transmission matrices of the whole aerosol/Rayleigh mixed layer (“AR”) are expressed as a sequence of matrix operations for each azimuthal component $m$ (Xu et al., 2010):

$$ \left\{ \begin{aligned} (2 - \delta_{0m}) \mathbf{R}_{\text{M,AR}}^{(m)} &= \mathbf{E}_{\text{Refl}}^{(m)} [\mathbf{I}_d - \mathbf{T}_{\text{Refl}}^{(m)}]^{-1} \mathbf{\Pi}_0^{(m)} \\ (2 - \delta_{0m}) \mathbf{T}_{\text{M,AR}}^{(m)} &= \mathbf{E}_{\text{Trans}}^{(m)} [\mathbf{I}_d - \mathbf{T}_{\text{Trans}}^{(m)}]^{-1} \mathbf{\Pi}_0^{(m)} \end{aligned} \right. , \quad (B1) $$

where $\delta_{0m}$ is the Kronecker delta, $\mathbf{I}_d$ is the identity matrix, and $\mathbf{R}_{\text{M,AR}}^{(m)}$ and $\mathbf{T}_{\text{M,AR}}^{(m)}$ are the $m$th Fourier sine and cosine components of the mixed layer reflection and transmission matrices, respectively, namely $\mathbf{R}_{\text{M,AR}}^{(m)} = [\mathbf{R}_{\text{M,AR,c}}^{(m)}, \mathbf{R}_{\text{M,AR,s}}^{(m)}]^{\text{T}}$ and $\mathbf{T}_{\text{M,AR}}^{(m)} = [\mathbf{T}_{\text{M,AR,c}}^{(m)}, \mathbf{T}_{\text{M,AR,s}}^{(m)}]^{\text{T}}$. Analytical expressions for $\mathbf{\Pi}_0^{(m)}$, $\mathbf{E}_{\text{Refl}}^{(m)}$, $\mathbf{E}_{\text{Trans}}^{(m)}$, $\mathbf{R}_{\text{Refl}}^{(m)}$ and $\mathbf{T}_{\text{Trans}}^{(m)}$ have been given by Xu et al. (2010) as a function of optical depth, phase matrix, and SSA for mixed Rayleigh and aerosol scattering (Eqs. A7–A9). Including the contributions of single scattering, $\mathbf{R}_{\text{S,AR}}^{(m)}$ and $\mathbf{T}_{\text{S,AR}}^{(m)}$ give the total reflection and transmission matrices of the mixed layer, namely

$$ \mathbf{R}_{\text{AR}}^{(m)} = \mathbf{R}_{\text{M,AR}}^{(m)} + \mathbf{R}_{\text{S,AR}}^{(m)} \quad (B2) $$

$$ \mathbf{T}_{\text{AR}}^{(m)} = \mathbf{T}_{\text{M,AR}}^{(m)} + \mathbf{T}_{\text{S,AR}}^{(m)}. \quad (B3) $$

Equation (B1) is the basic form of the Markov chain method. The majority of computational time is spent in computing the matrix inverse $[\mathbf{I}_d - \mathbf{X}^{(m)}]^{-1}$, with $\mathbf{X}$ being $T_{\text{Refl}}$ or $T_{\text{Trans}}$. To gain computational efficiency, the chain-to-chain adding strategy is applied to reduce the matrix dimension via subgrouping the layers (Esposito, 1979), and a truncated

Neumann series expansion is applied to approximate the matrix inverse.

$$ [\mathbf{I}_d - \mathbf{X}^{(m)}]^{-1} \approx \mathbf{I}_d + \sum_{n=1}^{N_{\text{max}}} \prod_{i=1}^{n} \mathbf{X}_i^{(m)} \quad (B4) $$

Setting 3–4 sublayers for each subgroup, fast convergence and accuracy of matrix inverse computation is usually achieved by using the first 3–4 series terms of Eq. (B4) (namely $N_{\text{max}} = 3$ or 4).

### B2 Coupling with doubling method for RT in atmosphere system

The reflection and transmission matrices of the two Rayleigh-scattering layers above and below the mixed layer, $(\mathbf{R}_R, \mathbf{T}_R)$, are computed using the doubling method (Hansen, 1971). Together with the reflection matrix of the mixed layer $(\mathbf{R}_{\text{AR}})$ computed from the Markov chain, a set of reflection and transmission matrices $(\mathbf{R}_{\text{atmos}}, \mathbf{T}_{\text{atmos}})$ for TOA illumination is obtained by applying the adding method twice (e.g., using Eq. 3 of Lacis and Hansen, 1974): two adjacent layers each time. In a similar way, another set of reflection and transmission matrices $(\mathbf{R}_{\text{atmos}}^*, \mathbf{T}_{\text{atmos}}^*)$ corresponding to illuminations from bottom of the atmosphere is evaluated by switching the location of the illumination sources from the top to the bottom of the mixed layer to evaluate $(\mathbf{R}_{\text{AR}}^*, \mathbf{T}_{\text{AR}}^*)$ and get $(\mathbf{R}_R^*, \mathbf{T}_R^*)$ from $(\mathbf{R}_R, \mathbf{T}_R)$ using the symmetric relationship (Hansen, 1970), then using the adding method to couple them (e.g., Eq. 4 of Lacis and Hansen 1974).

## Appendix C: Modeling radiative transfer in ocean system

In the five-layer CAOS system illustrated in Fig. 1, the ocean system is composed of the ocean medium and the air-water interface. The diffuse reflection matrix of the ocean medium and the reflection and transmission matrices of the air-water interface need to be known before they are coupled to evaluate the diffuse field at the top of ocean.

### C1 Extended adding-doubling method

Evaluation of the reflection matrix of the ocean system follows a similar methodology as the atmosphere system. However, instead of considering the contributions by molecules and aerosols, RT in the ocean involves scattering and absorption by sea water, CDOM, and phytoplankton and their covariant particles. Evaluation of the IOPs of these components relies on a simplified bio-optical model described in Appendix D, which determines absorption and scattering of CDOM and phytoplankton particles, then bulk optical depth $\tau_{\text{ocean}}$, phase matrix $\mathbf{P}_{\text{ocean}}$, and single-scattering albedo $\omega_{\text{ocean}}$ as a function of Chl a concentration. We also assume that the ocean components have a uniform vertical
<figure><img src="image_21.png" /><figcaption>**Figure C1.** Comparison of top-of-ocean radiance (L) and degree of linear polarization (DoLP) computed by the extended adding-doubling model (solid lines in the left two panels) and successive-orders-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 \times (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 \times (DoLP_{EAD} - DoLP_{SOS})$. The chlorophyll concentration is $[\text{Chl}_a] = 0.30 \text{ mg m}^{-3}$ and the solar zenith angle is $60^{\circ}$. The ocean surface is roughened by a wind of speed $7 \text{ m s}^{-1}$ and ocean optical thickness is set to 10. An arbitrary combination of refractive index and slope parameter ($n_p = 1.05, \gamma_p = 3.71$) is chosen to compute the Fourier-Forland phase function. The results are plotted for viewing angles ($\theta_v$) increasing from 0 to $87^{\circ}$ with an angular step of $3^{\circ}$; the five azimuthal planes ($\phi_v$) are 0, 45, 90, 135 and $180^{\circ}$ (shown in black, red, blue, green, and cyan, respectively) with respect to the principal plane (namely O-XZ in Fig. 1).</figcaption></figure>

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 $\mathbf{R}_{\text{ocean}}$, which depends on $\tau_{\text{ocean}}$, $\omega_{\text{ocean}}$, and $\mathbf{P}_{\text{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 ($\mathbf{R}_W, \mathbf{T}_W$) corresponding to downwelling incident light (in air) and another set of matrices ($\mathbf{R}_W^*, \mathbf{T}_W^*$) corresponding to upwelling incident light (in water) are then determined. In accordance with the adding method, two operators $\mathbf{Q}$ and $\mathbf{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.

$$Q_1 = \mathbf{R}_W^* \mathbf{R}_{\text{ocean}} \quad \text{(C1a)}$$

$$Q_n = Q_1 Q_{n-1} \quad \text{(C1b)}$$

$$S = \sum_{n=1}^{\infty} Q_n \quad \text{(C1c)}$$

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

$$\mathbf{D} = \mathbf{T}_W + \mathbf{S} \mathbf{T}_W \quad \text{(C1d)}$$

$$\mathbf{U} = \mathbf{R}_{\text{ocean}} \mathbf{D}. \quad \text{(C1e)}$$

The reflection matrix describing the upwelling diffusely reflected light leaving the ocean-air interface is

$$\mathbf{R}_{\text{OS}}^{\text{Bio, NR}} = \mathbf{T}_W^* \mathbf{U}, \quad \text{(C1f)}$$

where the superscript NR over $\mathbf{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-of-ocean 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 \text{ mg m}^{-3}$, solar zenith angle to $60^{\circ}$, surface wind speed to $7 \text{ m s}^{-1}$, and ocean optical thickness to 10. Using 40 streams in the half plane of $0 \le \mu \le 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
by Raman scattering, namely,

$$ \frac{R_{\text{rs}}^{\text{Raman}}}{R_{\text{rs}}^{\text{NR}}} = \zeta(\lambda) \frac{R_{\text{rs}}^{\text{Total}}(440)}{R_{\text{rs}}^{\text{Total}}(550)} + \xi_1(\lambda) [R_{\text{rs}}^{\text{Total}}(550)]^{\xi_2(\lambda)}, \quad (\text{C2}) $$

where $R_{\text{rs}}^{\text{Total}}$ is the total remote sensing reflectance as a sum of Raman scattering ($R_{\text{rs}}^{\text{Raman}}$) and non-Raman scattering ($R_{\text{rs}}^{\text{NR}}$), and $\zeta, \xi_1$, and $\xi_2$ are model parameters for empirical correction. Assuming an isotropic distribution of the radiance contributed by Raman scattering, the corrected reflection matrix of ocean and air-water interface system for use by Eq. (1) is

$$ \mathbf{R}_{\text{WL}}^{\text{Bio}} = \pi \left[ \mathbf{R}_{\text{OS}}^{\text{Bio, NR}}(\theta_v, \phi_v; \theta_0, \phi_0) + \mathbf{R}_{\text{WL,11}}^{\text{Bio,Raman}} \mathbf{D}_0 \right]. \quad (\text{C3}) $$

Since the two reference spectral bands at 440 and 550 nm in Eq. (C2) are close to the AirMSPI bands at 445 and 555 nm, $R_{\text{rs}}^{\text{Total}}(440)$ and $R_{\text{rs}}^{\text{Total}}(550)$ are directly replaced by $R_{\text{rs}}^{\text{Total}}(445)$ and $R_{\text{rs}}^{\text{Total}}(555)$ in our calculation. The parameters $\zeta, \xi_1$, and $\xi_2$ for the other AirMSPI bands are obtained by interpolating the values listed for the SeaWiFS bands in Lee et al. (2013). Fluorescence is neglected in our RT model due to its tiny contribution to TOA signals over open ocean, though it is known to have some impact on nLw at 685 nm (Gordon, 1979).

## Appendix D: A simplified bio-optical model

As indicated in the last two terms of Eq. (1), our water-leaving radiance model consists of two parts. The first part ($R_{\text{WL}}^{\text{Bio}}$) is a physically based term, which is dependent on a single parameter (namely Chl a concentration, or [Chl a]). The absorption and scattering properties of colored dissolved organic matter (CDOM or yellow substance) and phytoplankton and their covariant particles are dependent on this single parameter in a prescribed way. To deal with effects not captured by this model, a second empirical term ($\Delta a_{\text{WL}}$), represented as Lambertian water-leaving radiance adjustment with arbitrary spectral albedo, is added. This appendix describes the computation of ocean bulk optical properties as a function of Chl a concentration, which are then used as input to obtain $R_{\text{WL}}^{\text{Bio, NR}}$ via radiative transfer modeling, followed by a Raman-scattering correction to compute $R_{\text{WL}}^{\text{Bio}}$ (see Appendix C).

Pure sea water, CDOM, and phytoplankton and their covariant particles are considered to be the primary contributors to the oceanic absorption and scattering.

### D1 Pure sea water

The absorption coefficients of water ($\alpha_w$) are taken from the tabulated experimental data by Pope and Fry (1997). The scattering phase function of pure seawater (Morel 1974) is

$$ F_{w,11}(\Theta) = 4\pi \times 0.06225 \times (1 + 0.835\cos^2\Theta), \quad (\text{D1}) $$

where $\Theta$ is scattering angle. To obtain the other entries of the 4 × 4 phase matrix, we use ratios defined by Rayleigh scattering.

$$ F_{w,ij}(\Theta) = F_{w,ij}(\Theta) \times \frac{F_{R,ij}(\Theta)}{F_{R,11}(\Theta)} \quad (\text{D2}) $$
for $i \neq 1$ and $j \neq 1$

The depolarization factor of sea water is currently set to zero.

Invoking the Einstein–Smoluchowski theory of fluctuation scattering provides $\beta_w$ (Mobley, 1994), and the scattering coefficient for pure seawater is given by

$$ \beta_w = 0.00193(550/\lambda)^{4.32}. \quad (\text{D3}) $$

Due to the symmetry of scattering function of seawater around 90°, the backscattering coefficient $\beta_{bw}$ for the sea water is

$$ \beta_{bw} = 1/2\beta_w. \quad (\text{D4}) $$

### D2 Phytoplankton and their covariant particles

Phytoplankton and their covariant particles are assumed to conform to the hyperbolic (Junge) size-distribution, namely,

$$ n(r) = \frac{C}{r^{\gamma_p}}, \quad (\text{D5}) $$

with $n(r)dr$ being the number of particles per unit volume with radius between $r$ and $r+dr$ and $C$ is included to ensure proper normalization after integrating $n(r)$ over all sizes, namely,

$$ \int_0^\infty n(r)dr = 1. \quad (\text{D6}) $$

Knowing the real refractive index of particles ($n_p$) and the slope parameter ($\gamma_p$) of the hyperbolic size distribution, the Fournier–Forland (FF) scattering function which is a Mie-theory-based analytical approximation to the real scattering function of an ensemble of particles, can be determined (Fournier and Forland, 1994; Fournier and Jonasz, 1999), namely,

$$ F_{\text{FF}}(\Theta) = \frac{1}{(1 - \delta)^2 \delta^v} \left\{ v(1 - \delta) - (1 - \delta^v) + [\delta(1 - \delta^v) - v(1 - \delta)] \sin^{-2}(\Theta/2) \right\} + \frac{1 - \delta_{180}^v}{4(\delta_{180} - 1)\delta_{180}^v} (3\cos^2\Theta - 1) \quad (\text{D7}) $$

where $\Theta$ is the scattering angle, $\delta_{180}$ is the value of $\delta$ at $\Theta = 180^\circ$, and $\delta$ and $v$ are expressed as

$$ v = (3 - \gamma_p)/2 \text{ and } \delta = \frac{4}{3(n_p - 1)^2} \sin^2(\Theta/2). \quad (\text{D8}) $$
<figure><img src="image_22.png" /><figcaption>**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 ($\gamma_p$) as a function of backscattering efficiency ($B_{bp}$); and right panel: refractive index ($n_p$) as a function of the slope parameter ($\gamma_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 ($\gamma_p$) characterizing a Junge size distribution are then determined by solving Eqs. (D9) and (D11) numerically.</figcaption></figure>

With the FF scattering function, the backscattering efficiency can be obtained analytically (Mobley et al., 2002):

$$B_{bp} = 1 - \frac{1 - \delta_{90}^{v+1} - 0.5(1 - \delta_{90}^{v})}{(1 - \delta_{90}^{v})\delta_{90}^{v}}, \quad (D9)$$

where $\delta_{90}$ is $\delta$ evaluated at $\Theta = 90^{\circ}$.

Mobley (2002) found that the detailed shape of particle-scattering 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}$.

$$B_{bp} = \frac{1}{4\pi} \int_{\pi/2}^{\pi} F_{p}(\Theta) \sin \Theta d\Theta \quad (D10) \\ = 0.002 + 0.01(0.5 - 0.25\log_{10}[\text{Chl}_a])$$

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 $\gamma_p$ (Mobley et al., 2002), namely,

$$n_p = 1.01 + 0.1542(\gamma_p - 3). \quad (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 $\gamma_p$ – the two model parameters of the FF scattering function. Figure D1 illustrates the resulting relationships between $n_p$ and $B_{bp}$, between $\gamma_p$ and $B_{bp}$, and between $n_p$ and $\gamma_p$.

The absorption coefficients of phytoplankton and their covariant particles for $400 \le \lambda \le 700$ nm are parameterized by Bricaud et al. (1998) as

$$\alpha_p = A_p(\lambda)[\text{Chl}_a]^{E_p(\lambda)}. \quad (D12)$$

Integrated with Vasilkov et al.'s (2005) $A_p(\lambda)$ and $E_p(\lambda)$ spectra for $300 \le \lambda \le 400$ nm from coastal California water measurements, Morrison and Nelson's $A_p(\lambda)$ spectra for $300 \le \lambda \le 750$ nm from Bermuda Atlantic Time Series (BATS) site measurements (Morrison and Nelson, 2004), and setting $A_p(\lambda)$ and $E_p(\lambda)$ to 0 beyond 720 nm, the $A_p(\lambda)$ and $E_p(\lambda)$ spectra for $300 \le \lambda \le 1000$ nm are available from http://www.oceanopticsbook.info and adopted here.

The particle-scattering coefficients are evaluated based on the model by Morel and Maritorena (2001):

$$\beta_p = \beta_p(\lambda_0)(\lambda/\lambda_0)^\kappa, \quad (D13)$$

where $\beta_p(\lambda_0)$ is the scattering coefficient at the reference wavelength $\lambda_0$. Following Huot et al. (2008), we use $\lambda_0 = 660$ nm and

$$\beta_p(660) = 0.347[\text{Chl}_a]^{0.766}, \quad \text{with} \quad (D14)$$

$$\kappa = 0.5(\log_{10}[\text{Chl}_a] - 0.3), \quad (D15)$$

$$0.02 < [\text{Chl}_a] < 2 \text{ mg m}^{-3}$$

$$\kappa = 0, [\text{Chl}_a] < 0.02 \text{ mg m}^{-3}. \quad (D16)$$

## D3 CDOM

Absorption of CDOM ($\alpha_{CDOM}$) is estimated using the model of Bricaud et al. (1981):

$$\alpha_{CDOM}(\lambda) = \alpha_{CDOM}(\lambda_0) \exp[-S(\lambda - \lambda_0)], \quad (D17)$$

where for the reference wavelength $\lambda_0 = 440$ nm, $S = 0.014$ and, according to Bricaud et al. (1998),

$$\alpha_{CDOM}(440) = 0.2[\alpha_w(440) + \alpha_p(440)]. \quad (D18)$$

The scattering coefficient for CDOM is treated as zero in the present study.
## D4 Total inherent optical properties of sea water

Summarizing the contribution of all components gives the total absorption coefficient of ocean bulk ($\alpha_{ocean}$, cf. Zhai et al., 2010; Chowdhary et al., 2012) and the total scattering coefficient:

$$ \alpha_{ocean} = \alpha_w + \alpha_{CDOM} + \alpha_p \quad (D19) $$

$$ \beta_{ocean} = \beta_w + \beta_p. \quad (D20) $$

The total scattering function for sea water is

$$ P_{ocean,11}(\Theta) = [\beta_w F_w(\Theta) + \beta_p F_{FF}(\Theta)]/\beta_{ocean}. \quad (D21) $$

Polarized radiative transfer computations require the full phase matrix of bulk ocean scattering. To this purpose, we construct other phase matrix entries ($i \neq 1$ and $j \neq 1$) by using the ratio of measured sea water phase matrix entries, namely,

$$ P_{ocean,ij}(\Theta) = P_{ocean,11}(\Theta) \times [F_{VF,ij}(\Theta)/F_{VF,11}(\Theta)], \quad (D22) $$

where the ratio "$F_{VF,ij}(\Theta)/F_{VF,11}(\Theta)$" is taken from the averaged experimental measurements of Voss and Fry (1984).

Taking the geometric thickness of ocean as $\Delta H$, the total ocean optical thickness is then obtained from

$$ \tau_{ocean} = [\alpha_{CDOM} + (\alpha_w + \beta_w) + (\alpha_p + \beta_p)]\Delta H \quad (D23) \\ = (\alpha_{ocean} + \beta_{ocean})\Delta H. $$

The ocean single-scattering albedo is

$$ \omega_{ocean} = \beta_{ocean}/(\alpha_{ocean} + \beta_{ocean}). \quad (D24) $$

With $\tau_{ocean}$, $\omega_{ocean}$, and $P_{ocean}$, the reflection matrix of ocean and air-water interface system $R_{\lambda, WL}^{Bio, NR}$ is determined from radiative transfer modeling (see Appendix C). Further inclusion of a Raman-scattering correction via Eq. (C3) yields $R_{\lambda, WL}^{Bio}$ for the bio-optical model-based water-leaving radiances. As [Chl_a] is the only independent parameter in the simplified model, modeling errors are unavoidable. To account for them, the water-leaving radiances are adjusted in the second retrieval step by allowing $\Delta a_{WL} \neq 0$ in Eq. (1).

## Appendix E: Reflection and transmission matrices of the air-water interface

With the micro-facet assumption of oceanic wave structure, the polarized ocean surface reflectance (Tsang, 1985; Mishchenko, 1997) is modeled as

$$ \mathbf{R}_W = \frac{\pi P_S(z_x, z_y) S_h(\cos \theta_v, \cos \theta_i)}{4 \cos^4 \beta \cos \theta_i \cos \theta_v} \quad (E1) \\ \mathbf{r}(\pi - i_2)\mathbf{F}_r(n_w, \theta_i)\mathbf{r}(-i_1), $$

where $F_r$ is the Fresnel matrix for reflection as a function of the refractive index of water ($n_w$) and incidence angle $\theta_i$; the rotation matrices $\mathbf{r}(\pi - i_2)$ and $\mathbf{r}(i_1)$ are dependent on the angles $i_1$ and $i_2$, which account for the Stokes vector rotations between the meridian and reflection planes (Hovenier, 1969); $\Theta$ is the scattering angle; $\beta$ is the tilt angle of the facet surface normal; $S_h(\mu, \mu_0)$ is a shadowing function (Smith, 1967; Sancer, 1969; Brown, 1980); and $z_x$ and $z_y$ are the two components of surface slope:

$$ z_x = \frac{-\sin \theta_v \sin \phi}{\cos \theta_0 + \cos \theta_v} \quad (E2) $$

$$ z_y = \frac{\sin \theta_0 + \sin \theta_v \sin \phi}{\cos \theta_0 + \cos \theta_v}, \quad (E3) $$

where $\theta_0$ and $\theta_v$ are solar incidence and viewing angles, respectively, and $\phi$ is the relative azimuth angle. Without consideration of the wind direction, the wave slope probability distribution conforms to Cox and Munk's model (1954a, b):

$$ P_S(z_x, z_y) = \frac{1}{2\pi\sigma^2} \exp\left(-\frac{\tan^2\beta}{2\sigma^2}\right) \quad (E4) \\ \text{with } \tan^2\beta = z_x^2 + z_y^2, $$

where the slope variance is related to the wind speed ($W$) by $\sigma^2 = (0.003 + 0.00512W)/2$.

For the downwelling light, the transmission matrix (Zhai et al., 2010) is

$$ \mathbf{T}_W = \left[ \frac{n_w^2 \cos \theta_t \cos \theta_i}{(n_w \cos \theta_t - n_i \cos \theta_i)^2} \right] \frac{\pi P_S(z_x, z_y) S_h(\cos \theta_v, \cos \theta_i)}{4 \cos^4 \beta \cos \theta_i \cos \theta_v} \quad (E5) \\ \mathbf{r}(\pi - i_2)\mathbf{F}_t(n_w, \theta_i)\mathbf{r}(-i_1), $$

in which, compared to the reflection matrix, the Fresnel matrix for transmission $F_t$ is used and the extra term in the bracket accounts for the beam convergence when the light transmits from air through the air-water interface.

The Eqs. (E1) and (E5) also apply to the evaluation of reflection and transmission matrices $R_W^*$ and $T_W^*$ for the incidence of upwelling light from water into air after substituting $n_w$ with the new relative refractive index $1/n_w$.
*Acknowledgements.* The authors are grateful to Zia Ahmad at NASA Goddard Space Flight Center for providing the information on aerosol models used in MODIS ocean color retrieval and Jianwei Wei at Optical Oceanography Laboratory of University of Massachusetts Boston for discussing the AERONET Ocean Color product of normalized water-leaving radiance. This work was performed at the Jet Propulsion Laboratory, California Institute of Technology under contract with the National Aeronautics and Space Administration.

Edited by: A. Kokhanovsky

Reviewed by: two anonymous referees

## References

Ahmad, Z., Franz, B. A., McClain, C. R., Kwiatkowska, E. J., Werdell, J., Shettle, E. P., and Holben, B. N.: New aerosol models for the retrieval of aerosol optical thickness and normalized water-leaving radiances from the SeaWiFS and MODIS sensors over coastal regions and open oceans, Appl. Opt., 50, 626–626, 2011.

Aurin, D. A. and Dierssen, H. M.: Advantages and limitations of ocean color remote sensing in CDOM-dominated, mineral-rich coastal and estuarine waters, Remote Sens. Environ. 125, 181–197, 2012.

Bailey, S. W., Franz B. A., and Werdell P. J.: Estimation of near-infrared water-leaving reflectance for satellite ocean color data processing, Opt. Express, 18, 7521–7527, 2010.

Banzon, V. F., Gordon, H. R., Kuchinke, C. P., Antoine, D., Voss, K. J., and Evans, R. H.: Validation of a SeaWiFS dust-correction methodology in the Mediterranean Sea: Identification of an algorithm-switching criterion, Remote Sens. Environ., 113, 2689–2700, 2009.

Bergstrom, R. W., Pilewskie, P., Russell, P. B., Redemann, J., Bond, T. C., Quinn, P. K., and Sierau, B.: Spectral absorption properties of atmospheric aerosols, Atmos. Chem. Phys., 7, 5937–5943, doi:10.5194/acp-7-5937-2007, 2007.

Bodhaine, B. A., Wood, N. B., Dutton, E. G., and Slusser, J. R.: On Rayleigh optical depth calculations, J. Atmos. Oceanic Technol., 16, 1854–1861, 1999.

Bricaud, A., Morel, A., and Prieur, L.: Absorption by dissolved organic matter of the sea (yellow substance) in the UV and visible domains, Limnol. Oceanogr., 26, 43–53, 1981.

Bricaud, A., Morel, A., Babin, M., Allali, K., and Claustre, H.: Variations of light absorption by suspended particles with chlorophyll a concentration in oceanic (case 1) waters: analysis and implications for bio-optical models, J. Geophys. Res., 103, 31033–31044, 1998.

Brown, G. S.: Shadowing by non-Gaussian random surfaces, IEEE Trans. Ant. Prop., 28, 788–790, 1980.

Bruegge, C. J., Chrien, N. L., Ando, R. R., Diner, D. J., Abdou, W. A., Helmlinger, M. C., Pilorz, S. H., and Thome, K. J.: Validation of Multi-angle Imaging SpectroRadiometer (MISR) radiometric data products, IEEE Trans. Geosci. Remote Sens., 40, 1477–1492, 2002.

Bulgarelli, B., Kisselev, V., and Roberti, L.: Radiative transfer in the atmosphere-ocean system: the finite-element method, Appl. Opt., 38, 1530–1542, 1999.

Cahalan, R. F., Ridgway, W., Wiscombe, W. J., Bell, T. L., and Snider, J. B.: The albedo of fractal stratocumulus clouds, J. Atmos. Sci. 51, 2434–2455, 1994.

Chahine, M. T.: Determination of the temperature profile in an atmosphere from its outgoing radiance, J. Opt. Soc. Am., 58, 1634–1637, 1968.

Chami, M., Santer, R., and Dilligeard, E.: Radiative transfer model for the computation of radiance and polarization in an ocean-atmosphere system: polarization properties of suspended matter for remote sensing, Appl. Opt., 40, 2398–2416, 2001.

Chowdhary, J., Cairns, B., and Travis, L. D.: Contribution of water-leaving radiances to multiangle multispectral polarimetric observations over the open ocean: bio-optical model results for case 1 waters, Appl. Opt., 45, 5542–5567, 2006.

Chowdhary, J., Cairns, B., Waquet, F., Knobelspiesse, K., Ottaviani, M., Redemann, J., Travis, L., and Mishchenko M.: Sensitivity of multiangle, multispectral polarimetric remote sensing over open oceans to water-leaving radiance: Analyses of RSP data acquired during the MILAGRO campaign, Remote Sens. Environ., 118, 284–308, 2012.

Cox, C. and Munk, W.: Measurement of the roughness of the sea surface from photographs of the sun’s glitter, J. Opt. Soc. Amer., 44, 838–850, 1954.

Cox, C. and Munk, W.: Statistics of sea surface derived from sun glitter, J. Mar. Res., 13, 198–227, 1954.

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

Del Castillo, C. E., Platnick, S., Antoine, D., Balch, B., Behrenfeld, M., Boss, E. S., Cairns, B., Chowdhary, J., DaSilva, A., Diner, D., Dubovik, O., Franz, B., Frouin, R., Gregg, W., Huemmrich, K. F., Kahn, R., Marshak, A., Massie, S., McClain, C., McNaughton, C., Meister, G., Mitchell, G., Muller-Karger, F., Puschell, J., Riedi, J., Siegel, D., Wang, M., and Werdell, J.: Pre-Aerosol, Clouds, and Ocean Ecosystem (PACE) Mission, Science Definition Team Report, NASA, 2012.

Deuzé, J. L., Herman, M., and Santer, R.: Fourier series expansion of the transfer equation in the atmosphere-ocean system, J. Quant. Spectrosc. Radiat. Transfer, 41, 483–494, 1989.

Deuzé, J. L., Goloub, P., Herman, M., Marchand, A. G., Perry, Susana, S., and Tanré, D.: Estimate of the Aerosol Properties over the Ocean with POLDER, J. Geophys. Res., 105, 15329–15346, 2000.

Diner, D. J., Beckert, J. C., Reilly, T. H., Bruegge, C. J., Conel, J. E., Kahn, R. A., Martonchik, J. V., Ackerman, T. P., Davies, R., Gerstl, S. A. W., Gordon, H. R., Muller, J.-P., Myneni, R. B., Sellers, P. J., Pinty, B., and Verstraete, M.: Multi-angle Imaging SpectroRadiometer (MISR) instrument description and experiment overview, IEEE Trans. Geosci. Remote Sens., 36, 1072–1087, 1998.

Diner, D. J., Hodos, R. A., Davis, A. B., Garay, M. J., Martonchik, J. V., Sanghavi, S. V., von Allmen, P., Kokhanovsky, A. A., and Zhai, P. W.: An optimized approach for aerosol retrievals using simulated MISR radiances, Atmos. Res., 116, 1–14, 2011.

Diner, D. J., Xu, F., Garay, M. J., Martonchik, J. V., Rheingans, B. E., Geier, S., Davis, A., Hancock, B. R., Jovanovic, V. M., Bull, M. A., Capraro, K., Chipman, R. A., and McClain, S. C.: The Airborne Multiangle SpectroPolarimetric Imager (AirMSPI): a
new tool for aerosol and cloud remote sensing, Atmos. Meas. Tech., 6, 2007–2025, doi:10.5194/amt-6-2007-2013, 2013.

Dubovik, O. and King, M. D.: A flexible inversion algorithm for retrieval of aerosol optical properties from Sun and sky radiance measurements, J. Geophys. Res., 105, 20673–20696, 2000.

Dubovik, O., Smirnov, A., Holben, B. N., King, M. D., Kaufman, Y. J., Eck, T. F., and Slutsker, I.: Accuracy assessments of aerosol optical properties retrieved from aerosol robotic network (AERONET) Sun and sky radiance measurements, J. Geophys. Res., 105, 9791–806, 2000.

Dubovik, O.: Optimization of Numerical Inversion in Photopolari- metric Remote Sensing, in: Photopolarimetry in Remote Sens- ing, edited by: Videen, G., Yatskiv, Y., and Mishchenko, M., Kluwer Academic Publishers, Dordrecht, the Netherlands, 65– 106, 2004.

Dubovik, O., Sinyuk, A., Lapyonok, T., Holben, B. N., Mishchenko, M., Yang, P., Eck, T. F., Volten, H., Muñoz, O., Veihelmann, B., van der Zande, W. J., Leon, J.-F., Sorokin, M., and Slutsker, I.: Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust, J. Geophys. Res., 111, D11208, doi:10.1029/2005JD006619, 2006.

Dubovik, O., Lapyonok, T., Kaufman, Y. J., Chin, M., Ginoux, P., Kahn, R. A., and Sinyuk, A.: Retrieving global aerosol sources from satellites using inverse modeling, Atmos. Chem. Phys., 8, 209–250, doi:10.5194/acp-8-209-2008, 2008.

Dubovik, O., Herman, M., Holdak, A., Lapyonok, T., Tanré, D., Deuzé, J. L., Ducos, F., Sinyuk, A., and Lopatin, A.: Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite obser- vations, Atmos. Meas. Tech., 4, 975–1018, doi:10.5194/amt-4- 975-2011, 2011.

Dubovik, O., Lapyonok, T., Litvinov, P., Herman, M., Fuertes, D., Ducos, F., Lopatin, A., Chaikovsky, A., Torres, B., Derimian, Y., Huang, X., Aspetsberger, M., and Federspiel, C.: GRASP: a ver- satile algorithm for characterizing the atmosphere, SPIE: News- room, Sep., 2014.

Duforêt, L., Frouin R., and Dubuisson, P.: Importance and estima- tion of aerosol vertical structure in satellite ocean-colour remote sensing, Appl. Opt., 46, 1107–1119, 2007.

Esposito, L. W. and House, L. L.: Radiative transfer calculated by a Markov chain formalism, Astrophys. J. 219, 1058–1067, 1978.

Esposito, L. W.: An “adding” algorithm for the Markov chain for- malism for radiation transfer, Astrophys. J., 233, 661–663, 1979.

Evans, K. F. and Stephens, G. L.: A new polarized atmospheric ra- diative transfer model, J. Quant. Spectrosc. Radiat. Transfer, 46, 413–423, 1991.

Fischer, J. and Grassl, H.: Radiative transfer in an atmosphere- ocean system: an azimuthally dependent matrix-operator ap- proach, Appl. Opt., 23, 1032–1039, 1984.

Fourier, G. and Forand J. L.: Analytic phase function for ocean wa- ter, in: Ocean optics XII, edited by: Jaffe, J. S., Proceedings of the SPIE, 2258, 194–201, 1994.

Fournier, G. and Jonasz, M.: Computer-based underwater imaging analysis, in: Airborne and in-water underwater imaging, edited by: Gilbert, G., Proceedings of the SPIE, 3761, 62–77, 1999.

Franz, B. A., Bailey, S. W., Werdell, P. J., and McClain, C. R.: Sensor-independent approach to the vicarious calibration of satellite ocean color radiometry, Appl. Opt., 46, 5068–5082, 2007.

Gatebe, C., King, M. D., Lyapustin, A., Arnold, G. T., and Rede- mann, J.: Airborne Spectral Measurements of Ocean Directional Reflectance, J. Atmos. Sci., 62, 1072–1092, 2005.

Gergely M. and Zibordi, G.: Assessment of AERONET LWN un- certainties, Metrologia, 51, 40–47, 2014.

Gordon, H. R.: Diffuse reflectance of the ocean: the theory of its augmentation by chlorophyll a fluorescence at 685 nm, Appl. Opt., 18, 1161–1166, 1979.

Gordon, H. R. and Wang, M.: Retrieval of water-leaving radiance and aerosol optical thickness over the oceans with SeaWiFS: A preliminary algorithm, Appl. Opt., 33, 443–452, 1994.

Gordon, H. R.: Atmospheric correction of ocean color imagery in the Earth observing era, J. Geophys. Res., 102, 17081–17106, 1997.

Hansen, J. E.: Multiple scattering of polarized light in planetary at- mospheres – Part I: The doubling method, J. Atmos. Sci., 28, 120–125, 1971.

Hasekamp, O. P. and Landgraf, J.: Linearization of vector ra- diative transfer with respect to aerosol properties and its use in satellite remote sensing, J. Geophys. Res., 110, D04203, doi:10.1029/2004JD005260, 2005.

Hasekamp, O. P. and Landgraf, J.: Retrieval of aerosol properties over land surfaces: capabilities of multiple-viewing-angle inten- sity and polarization measurements, Appl. Opt., 46, 3332–3344, 2007.

Hasekamp, O., Litvinov, P., and Butz, A.: Aerosol proper- ties over the ocean from PARASOL multi-angle photopo- larimetric measurements, J. Geophys. Res., 116, D14204, doi:10.1029/2010JD015469, 2011.

He, X., Bai, Y., Pan, D., Tang, J., and Wang, D.: Atmospheric correction of satellite ocean color imagery using the ultraviolet wavelength for highly turbid waters, Opt. Express, 20, 20754– 20770, 2012.

Hovenier, J. W.: Symmetry relationships for scattering of polarized light in a slab of randomly oriented particles, J. Atmos. Sci., 26, 488–499, 1969.

Hu, C., Carder, K. L., and Muller-Karger, F. E.: Atmospheric correc- tion of SeaWiFS imagery over turbid coastal waters: A practical method, Rem. Sens. Environ., 74, 195–206, 2000.

Huot, Y., Morel, A., Twardowski, M. S., Stramski, D., and Reynolds, R. A.: Particle optical backscattering along a chloro- phyll gradient in the upper layer of the eastern South Pa- cific Ocean, Biogeosciences, 5, 495–507, doi:10.5194/bg-5-495- 2008, 2008.

IOCCG, Remote Sensing of Inherent Optical Properties: Funda- mentals, Tests of Algorithms, and Applications, in: Reports of the International Ocean-Colour Coordinating Group, No. 5, edited by: Lee, Z. P., p. 126, IOCCG, Dartmouth, Canada, 2006.

Jin, Z. and Stamnes, K.: Radiative transfer in nonuniformly refract- ing layered media: atmosphere-ocean system, Appl. Opt., 33, 431–442, 1994.

Jin, Z., Charlock, T. P., Rutledge, K., Stamnes, K., and Wang, Y.: Analytical solution of radiative transfer in the coupled atmosphere-ocean system with a rough surface, Appl. Opt., 45, 7443–7455, 2006.

Kahn, R. A., Gaitley, B. J., Garay, M. J., Diner, D. J., Eck, T. F., Smirnov, A., and Holben, B. N.: Multiangle Imaging SpectroRa- diometer global aerosol product assessment by comparison with
the Aerosol Robotic Network, J. Geophys. Res.-Atmos., 115, D23209, doi:10.1029/2010JD014601, 2010.

Kattawar, G. W. and Adams, C. N.: Stokes vector calculations of the submarine light field in an atmosphere-ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on the radiance and polarization, Limnol. Oceanogr., 34, 1463–1472, 1989.

Koepke, P.: Effective reflectance of oceanic whitecaps, Appl. Opt., 23, 1816–1824, 1984.

Kokhanovsky, A. A.: The modern aerosol retrieval algorithms based on the simultaneous measurements of the intensity and polarization of reflected solar light: A review, Front. Environ. Sci. 3, doi:10.3389/fenvs.2015.00004, 2015.

Kokhanovsky, A. A., Davis A. B., Cairns B., Dubovik O., Hasekamp O., Sano I., Mukai S., Rozanov, V. V., Litvinov, P., Lapyonok, T., Kolomiets, I. S., Oberemok, Y. A., Savenkov, S., Martin, W., Wasilewski, A., Di Noia, A., Stap, F. A., Rietjens, J., Xu, F., Natraj, V., Duan, M., Cheng, T., and Munro, R.: Space-based remote sensing of aerosols: The multi-angle spectro-polarimetric frontier, Earth Sci. Rev., 145, 85–116, 2015.

Koven, C. D. and Fung, I.: Inferring dust composition from wavelength-dependent absorption in Aerosol Robotic Network (AERONET) data, J. Geophys. Res., 111, D14205, doi:10.1029/2005JD006678, 2006.

Lacis, A. A. and Hansen, J. E.: A parameterization for the absorption of solar radiation in the Earth's atmosphere, J. Atmos. Sci. 31, 118–133, 1974.

Langmore, I., Davis, A. B., and Bal, G.: Multi-pixel retrieval of structural and optical parameters in a 2D scene with a path-recycling Monte Carlo forward model and a new Bayesian inference engine, IEEE Trans. Geosci. Remote Sens., 51, 2903–2919, 2013.

Lee, Z. P., Hu, C., Shang, S. L., Du, K. P., Lewis, M., Arnone, R., and Brewin, R.: Penetration of UV-Visible solar light in the global oceans: Insights from ocean color remote sensing, J. Geophys. Res., 118, 4241–4255, doi:10.1002/jgrc.20308, 2013.

Levenberg, K.: A method for the solution of certain non-linear problems in least squares, Quart. Appl. Math., 2, 164–168, 1944.

Limbacher, J. A. and Kahn, R. A.: MISR research-aerosol-algorithm refinements for dark water retrievals, Atmos. Meas. Tech., 7, 3989–4007, doi:10.5194/amt-7-3989-2014, 2014.

Lyapustin, A. I. and Muldashev, T. Z.: Generalization of Marshak boundary condition for non-Lambert reflection, J. Quant. Spectrosc. Radiat. Transfer, 67, 457–446, 2001.

Marquardt, D.: An algorithm for least-squares estimation of nonlinear parameters, SIAM J. Appl. Math., 11, 431–441, 1963.

Masuda, K. and Takashima, T.: Computational accuracy of radiation emerging from the ocean surface in the model atmosphere-ocean System, Pap. Meteorol. Geophys., 37, 1–13, 1986.

Mishchenko, M. I. and Travis, L. D.: Satellite retrieval of aerosol properties over the ocean using polarization as well as intensity of reflected sunlight, J. Geophys. Res. 102, 16989–17013, 1997.

Mobley, C. D.: Light and Water: Radiative Transfer in Natural Waters, Academic, 1994.

Mobley, C. D., Sundman, L. K., and Boss, E.: Phase function effects on oceanic light fields, Appl. Opt., 41, 1035–1050, 2002.

Mobley, C. D., Zhang, H., and Voss, K. J.: Effects of optically shallow bottoms on upwelling radiances: bidirectional reflectance

distribution function effects, Limnol. Oceanog., 48, 337–345, 2003.

Mobley, C. D. and Sundman, L. K.: HydroLight Technical Documentation, Sequoia Scientific, Inc., http://www.hydrolight.info, 2008.

Mobley, C. D.: EcoLight-S 1.0 User's Guide and Technical Documentation, Sequoia Scientific, Inc., http://www.hydrolight.info, 2011a.

Mobley, C. D.: Fast light calculations for ocean ecosystem and inverse models, Opt. Express, 19, 18927–18944, 2011b.

Morel, A.: Optical properties of pure water and pure sea water, in: Optical aspects of oceanography, Academic Press, London, 1–24, 1974.

Morel, A. and Maritorena, S.: Bio-optical properties of oceanic waters: a reappraisal, J. Geophys. Res., 106, 7163–7180, 2001.

Morel, A., Antoine, D., and Gentili, B.: Bidirectional reflectance of oceanic waters: accounting for Raman emission and varying particle scattering phase function, Appl. Opt., 41, 6289–6306, 2002.

Morrison, J. R. and Nelson, N. B.: Seasonal cycle of phytoplankton UV absorption and the Bermuda Atlantic Time-series Study (BATS) site, Limnol. Oceanogr., 49, 215–224, 2004.

Moulin, C., Gordon, H. R., Chomoko, R. M., Banzon, V. F., and Evans, R. H.: Atmospheric correction of ocean color imagery through thick layers of Saharan dust, Geophys. Res. Lett., 28, 5–8, 2001.

Nakajima, T. and Tanaka, M.: Effect of wind-generated waves on the transfer of solar radiation in the atmosphere-ocean system, J. Quant. Spectrosc. Radiat. Transfer, 29, 521–537, 1983.

Ortega, J. M. and Rheinboldt, W. C.: Iterative solution of nonlinear equations in several variables, Academic Press, New York-London, 1970.

Phillips, B. L.: A technique for the numerical solution of certain integral equations of the first kind, J. Assoc. Comp. Mach., 9, 84–97, 1962.

Pope, R. M. and Fry, E. S.: Absorption spectrum (380–700 nm) of pure water: II. Integrating cavity measurements, Appl. Opt., 36, 8710–8723, 1997.

Remer, L. A., Kaufman, Y. J., Tanré, D., Mattoo, S., Chu, D. A., Martins, J. V., Li, R.-R., Ichoku, C., Levy, R. C., Kleidman, R. G., Eck, T. F., Vermote, E., and Holben, B. N.: The MODIS aerosol algorithm, products and validation, J. Atmos. Sci., 62, 947–973, 2005.

Remer, L. A., Tanré, D., and Kaufman, Y. J.: Algorithm for remote sensing of tropospheric aerosol from MODIS: Collection5, MODIS Algorithm Theoretical Basis Document, http://modis-atmos.gsfc.nasa.gov/_docs/MOD04-MYD04_ATBD_C005.pdf, 2006.

Russell, P. B., Bergstrom, R. W., Shinozuka, Y., Clarke, A. D., DeCarlo, P. F., Jimenez, J. L., Livingston, J. M., Redemann, J., Dubovik, O., and Strawa, A.: Absorption Angstrom Exponent in AERONET and related data as an indicator of aerosol composition, Atmos. Chem. Phys., 10, 1155–1169, doi:10.5194/acp-10-1155-2010, 2010.

Sancer, M. I.: Shadow-corrected electromagnetic scattering from a randomly rough surface, IEEE Trans. Ant. Prop., 17, 577–585, 1969.

Sayer, A. M., Thomas, G. E., and Grainger, R. G.: A sea surface reflectance model for (A)ATSR, and application to aerosol re-
trievals, Atmos. Meas. Tech., 3, 813–838, doi:10.5194/amt-3-813-2010, 2010.

Sokolik, I. N. and Toon, O. B.: Incorporation of mineralogical com-
position into models of the radiative properties of mineral aerosol
from UV to IR wavelengths, J. Geophys. Res., 104, 9423–9444,
1999.

Schollaert, S. E., Yoder, J. A., O'Reilly, J. E., and Westphal,
D. L.: Influence of dust and sulfate aerosols on ocean color
spectra and chlorophyll a concentrations derived from Sea-
WiFS off the U.S. east coast, J. Geophys. Res., 108, 3191,
doi:10.1029/2000JC000555, 2003.

Siegel, D. A., Wang, M., Maritorena, S., and Robinson, W.: At-
mospheric correction of satellite ocean color imagery: the black
pixel assumption, Appl. Opt., 39, 3582–3591, 2000.

Smith, B. G.: Geometrical shadowing of a random rough surface,
IEEE Trans. Ant. Prop., 17, 15668–15671, 1967.

Sommersten, E. R., Lotsberg, J. K., Stamnes, K., and Stamnes, J. J.: Discrete ordinate and Monte Carlo simulations for polarized radiative transfer in a coupled system consisting of two medium with different refractive indices, J. Quant. Spectrosc. Radiat. Transfer, 111, 616–633, 2010.

Spurr, R. J. D.: VLIDORT: A linearized pseudo-spherical vector
discrete ordinate radiative transfer code for forward model and
retrieval studies in multilayer multiple scattering media, J. Quant.
Spectrosc. Radiat. Transfer, 102, 316–342, 2006.

Stokes, G. G.: On the intensity of the light reflected from or trans-
mitted through a pile of plates, Proc. Roy. Soc. London, 11, 545–
556, 1862.

Sun, W. and Lukashin, C.: Modeling polarized solar radia-
tion from the ocean-atmosphere system for CLARREO inter-
calibration applications, Atmos. Chem. Phys., 13, 10303–10324,
doi:10.5194/acp-13-10303-2013, 2013.

Tikhonov, A. N.: On the solution of incorrectly formulated prob-
lems and the regularization method, Dokl. Akad. Nauk SSSR,
151, 501–504, 1963.

Tomasi, C., Vitale, V., Petkov, B., Lupi, A., and Cacciari, A.: Im-
proved algorithm for calculations of Rayleigh-scattering optical
depth in standard atmospheres, Appl. Opt., 44, 3320–3341, 2005.

Tsang, L., Kong, J. A., and Shin, R. T.: Theory of Microwave Re-
mote Sensing, John Wiley, New York, 1985.

Twomey, S.: On the numerical solution of Fredholm integral equa-
tions of the first kind by the inversion of the linear system pro-
duced by quadrature, J. Assoc. Comp. Mach., 10, 97–101, 1963.

Twomey, S.: Comparison of constrained linear inversion and an it-
erative nonlinear algorithm applied to the indirect estimation of
particle size distributions, J. Comput. Phys., 18, 188–200, 1975.

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

van de Hulst, H. C.: Light Scattering by Small Particles, Dover, New York, 1981.

Voss, K. J. and Fry, E. S.: Measurement of the Mueller matrix for
ocean water, Appl. Opt., 23, 4427–4439 1984.

Xu, F., Davis, A. B., West, R. A., and Esposito, L. W.: Markov
Chain formalism for polarized light transfer in plane-parallel
atmospheres, with numerical comparison to the Monte Carlo
method, Opt. Express, 19, 946–967, 2010.

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

Xu, F., Davis, A. B., Suniti, S. V., Martonchik, J. V., and Diner,
D. J.: Linearization of Markov chain formalism for vector radia-
tive transfer in a plane-parallel atmosphere/surface system, Appl.
Opt., 51, 3491–3507, 2012.

Zhai, P. W., Hu, Y., Trepte, C. R., and Lucker, P. L.: A vector radia-
tive transfer model for coupled atmosphere and ocean systems
based on successive order of scattering method, Opt. Express,
17, 2057–2079, 2009.

Zhai, P. W., Hu, Y., Chowdhary, J., Trepte, C. R., Lucker, P. L., and Josset, D. B.: A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface, J. Quant. Spectrosc. Radiat. Transfer, 111, 1025–1040, 2010.

Zhai, P. W., Hu, Y., Winker, D. M., Franz, B. A., and Boss, E.: Con-
tribution of Raman scattering to polarized radiation field in ocean
waters, Opt. Express, 23, 23582–23596, 2015.

Zibordi, G., Mélin, F., Berthon, J-F, Holben, B., Slutsker, I., Giles, D., D'Alimonte, D., Vandemark, D., Feng, H., Schuster, G., Fabbri, B. E., Kaitala, S., and Seppälä, J.: AERONET-OC: A Network for the Validation of Ocean Color Primary Products, J. Atmos. Oceanic Technol., 26, 1634–1651, 2009.