## Abstract

A vector radiative transfer model has been developed for coupled atmosphere and ocean systems based on the Successive Order of Scattering (SOS) Method. The emphasis of this study is to make the model easy-to-use and computationally efficient. This model provides the full Stokes vector at arbitrary locations which can be conveniently specified by users. The model is capable of tracking and labeling different sources of the photons that are measured, e.g. water leaving radiances and reflected sky lights. This model also has the capability to separate florescence from multi-scattered sunlight. The *δ* - fit technique has been adopted to reduce computational time associated with the strongly forward-peaked scattering phase matrices. The exponential - linear approximation has been used to reduce the number of discretized vertical layers while maintaining the accuracy. This model is developed to serve the remote sensing community in harvesting physical parameters from multi-platform, multi-sensor measurements that target different components of the atmosphere-oceanic system.

©2009 Optical Society of America

## 1. Introduction

The radiative transfer equation (RTE) describes the propagation of radiation in turbid media phenomenologically [1, 2, 3, 4, 5, 6]. Recently, the rederivation of the RTE from the electromagnetic wave theory establishes a new interpretation of the radiative transfer theory [7]. The RTE is a differential-integral equation which generally has no analytical solution for realistic medium configurations and boundary conditions. Therefore, efficient numerical solutions are crucial for the applications of the radiative transfer theory. A non-complete list of these methods includes the discrete-ordinate (DISORT) method[1, 8, 9, 10, 11, 12, 13, 14, 15], the invariant embedding method [1, 16, 17, 18, 19, 20], the adding and doubling method [3, 21, 22], the matrix operator method [23, 24, 25, 26, 27, 28, 29], the spherical harmonics method [30, 31], the multicomponent method [32, 33], the spherical harmonics discrete ordinate method [34, 35], the FN method [36, 37], the Successive Order of Scattering (SOS) method [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48], and the Monte Carlo method [49, 50, 51, 52, 53, 54, 55, 56].

The solution for the RTE in a coupled Atmosphere - Ocean System (AOS) is particularly of interest in the ocean color remote sensing, the general circulation model, and many other studies. There are a few radiative transfer codes currently available for an AOS [26, 27, 32, 40, 57, 58, 59, 60, 61, 62, 63, 64, 65]. Based on the matrix operator method, [26] and [27] solve the scalar RTE in an AOS with wind induced rough ocean surface. The multi-components method developed by Zege et al. [32, 33] solves the RTE in an AOS including polarization. The model by Chami [40] is based on the SOS method which also includes polarization. A flat ocean surface is assumed and the circularly polarized component is ignored in this model. [60, 63] solve the scalar RTE using the DISORT method. Bulgarelli et al. has developed finite element method for solving the scalar RTE in an AOS with a smooth ocean surface [61]. Another important work is the Hydrolight code based on the invariant embedding method by Mobley [4, 59]. Polarization is not included in the Hydrolight. The model by Chowdhary et al. solves the vector radiative transfer equation in an AOS based on the adding method [64]. The rest of the list are mainly based on the Monte Carlo method [54, 57, 58, 62, 65]. The Monte Carlo technique is slow for remote sensing purposes because it requires an accumulation of statistics. In our classification, the models in [20, 45, 46, 47] do not solve the RTE in an AOS because they either treat the ocean surface as bidirectional reflectance surface or use a semi-analytical model for the ocean [66], which is not a real simulation of light propagation in ocean.

Most discrete ordinate methods (excluding the Monte Carlo method) apply the Fourier transformation to the azimuth dependence of the Stokes vector in RTE in order to reduce the dimension of differential equations, as the computation time for eigen solution is proportional to the third power of the dimension. For each Fourier component, Gaussian quadrature scheme is used to discretize the cosine of the zenith angle *μ* = cos *θ* for numerical integration of the scattering term in RTE. At an ocean surface, the refracted light in the ocean from the atmosphere is confined within the Fresnel cone, which only cover part of the downward solid angle. To match the boundary condition at the ocean surface, two sets of quadrature points are usually selected for ocean [32, 60, 63]. One corresponds to the refracted directions of the quadrature points in the atmosphere. The other covers the region outside the Fresnel cone. While the two sets of Gaussian quadratures in water are required in most models that treat the RTE as eigenvalue problem, it is not entirely clear whether that is a neccessity for the successive order of scattering approach. While less accurate, relaxing to one set of Gaussian quadrature may help reducing computational time and memory requirement.

For optically thin media (such as clear sky with moderate aerosol loading) and for highly absorbing media (such as ocean), the SOS method has advantages when compared with other methods. First, each order of scattering has clear physics interpretation. This results in efficient programing. For example, because of the fact that the source term associated with the refracted light into the water can be formulated analytically, two sets of quadrature points for the ocean are not mandatory. The same quadrature points can be used in the ocean just as in the atmosphere. In other words, the quadrature points in the ocean do not have to match Snel’s law [67] with the quadrature points in the atmosphere, as long as the angle change associated with refraction is properly represented in the scattering source term (scattered light at arbitrary quadrature angle). The use of one set of Gaussian quadrature points could decrease the memory requirements for the ocean by a rough factor of 1/4 because a major part of the memory in the computer code is allocated to store the phase matrices, whose number of elements are proportional to *P*
^{2}, where *P* is the number of quadrature points. The CPU time reduction, however, depends on the specific program implementation.

This paper documents the physics principle and detailed mathematical derivation of a Fortran 90 computing package for solving the vector RTE based on the SOS method. The current code is developed for a flat ocean surface. However, the technique can be extended to rough wave surface model [68]. Both the single and double ocean quadrature versions of the code are available to the public science community by contacting the corresponding author.

Section 2 outlines the vector RTE, while section 3 describes the SOS formulas for solving the vector RTE. Section 4 presents the treatment of the air-sea interface boundary condition. Section 5 shows the numerical procedure used in the code. Section 6 provides the validation and simulation results. The final section, section 7, presents the conclusions.

## 2. The vector radiative transfer equation

A beam light of arbitrary polarization may be represented by the Stokes vector **L** = (*I*,*Q*,*U*,*V*)*T*, where the four elements *I*, *Q*, *U*, and *V* are the Stokes parameters; the superscript *T* stands for the transpose of the vector. A formal definition of the Stokes parameters in a medium may be found in [7] (see Eq. 2.6.4). The mathematical definitions of the Stokes parameters depends on propagation direction, electric vector components, as well as the dielectric properties of the medium. The vector RTE which describes multiple light scattering in a plane-parallel medium can be written as :

where *τ* is the optical depth defined as the product of the extinction coefficient *β _{e}* and the path lengh

*I*;

*π*= cos(

*θ*) and

*θ*is the zenith angle;

*ϕ*is the azimuthal angle. All the zenith angles in this paper are defined in terms of the upward direction.

*θ*< 90° means that the direction is pointing upward and

*θ*> 90° when pointing downward. Both

*μ*and

*ϕ*are pertained to the propagation direction of radiance

**L**.

**S**is the source matrix:

where *ω* = *β _{s}*/

*β*is the single scattering albedo and

_{e}*β*is the scattering coefficient. Mishchenko et al. have shown that Eq. (2) is only valid in a macroscopically isotropic and mirror symmetric medium [7]. In this equation and hereafter the dot product is assumed between the matrix

_{s}**P**and vector

**L**. The same assumption applies if matrices are written next to each other.

**P**is the phase matrix defined as:

The rotation matrix **R** rotates the reference plane of the Stokes vector between the meridian plane and the scattering plane; *i*
_{1} and *i*
_{2} are the rotation angles [1, 69]; **F** is the scattering matrix for the turbid medium; the scattering angle Θ is a function of variables *μ*, *ϕ*, μ′, and *ϕ*′. Equation (1) does not contain thermal emissions.

The integral form of Eq.(1) is needed to apply the SOS method. Standard manipulation leads to [1, 5]:

where *τ _{l}* and

*τ*are the lower and upper limits of the medium under consideration, namely

_{u}*τ*<

_{l}*τ*<

*τ*. For the atmosphere,

_{u}*τ*= 0 and

_{l}*τ*=

_{u}*τ*

_{a}^{*}are the top and bottom of the atmosphere, respectively. For the ocean,

*τ*=

_{l}*τ*and

_{o}*τ*=

_{u}*τ*

_{o}^{*}are the top and bottom of the ocean, respectively. Please note that

*τ*=

_{o}*τ*

_{a}^{*}+ ϵ is just below the ocean surface, where e is an infinitesimal number.

In other words, the ocean surface is assumed to be a degenerated layer with no geometric thickness. The only impact of the ocean surface to the system is through the refraction and reflection of the dielectric interface, which is governed by Snel’s law and Fresnel formulas [57, 65]. **L**(0, *μ* < 0,*ϕ*) is 0; and **L**(τ^{*}
_{a},μ > 0,ϕ) includes the light reflected by the dielectric interface from the downward light in the atmosphere and the light transmitted through the interface from the ocean. **L**(*τ*
_{0}, *μ* < 0,*ϕ*) includes the light reflected by the dielectric interface from the upward light in the ocean and the refracted light through the dielectric interface from the atmosphere; and **L** (*τ _{o}*

^{*},

*μ*> 0,

*ϕ*) is the reflection of the ocean bottom.

## 3. The successive orders of scattering method

In the SOS method, the total radiance vector is expressed as a summation of contributions from photons scattered a number of times, from 0 to a maximum number of *N* [39, 47, 70, 71],

Substitution of Eq. (5) into Eqs.(4a) and (4b) leads to:

where

$$\phantom{\rule{7.6em}{0ex}}+{e}^{\left(2{\tau}_{a}^{*}-\tau \right)/{\mu}_{0}}\mathbf{P}\left(\tau ,\mu ,\varphi ,-{\mu}_{0},{\varphi}_{0}\right)\mathbf{r}\left(\pi -{\theta}_{0}\right)\}{\mathbf{E}}_{0},$$

It should be noted that the direct solar light has been taken out of Eqs. (6). The source terms caused by the direct solar light have been integrated and written out explicitly as **S**
_{1} in Eqs. (7). In other words, Eqs. (5) and (6) only contains the diffuse light in the turbid medium. In Eqs.(7), **E**
_{0} = (*E*
_{0},0,0,0)^{T} is the irradiance vector for the solar radiation; The matrices **r** and **t** are the Fresnel reflection and transmission matrices, respectively. The readers are refered to [57] for the exact expressions (see Eqs. (12–18) of [57]). *μ*
_{0} = cos(*θ*
_{0}) and *φ*
_{0} are the cosine of the solar zenith angle and the solar azimuthal angle, respectively. *μ*
_{0} < 0 because *θ*
_{0} > 90°. *μ*′_{0} = cos(*θ*′_{0}). *θ*′_{0} is the refracted angle for the solar incident angle *θ*
_{0}, i.e. sin(*π* - *θ*
_{0}) = *n _{w}* sin(

*π*-

*θ*′

_{0}) and

*n*is the refractive index of the water relative to air. The factor of |

_{w}*μ*

_{0}/

*μ*′

_{0}| takes the different projection areas of the solar beam at the two sides of the ocean surface into account. The radiance vectors at the boundaries are:

where **r**
^{*}(*μ*,*φ*,*μ*′,*φ*′) is the bidirectional reflectance matrix for the ocean bottom. For a Lambertian surface, **r**
^{*}
_{1,1} = *ρ* is a constant, and **r**
^{*}
_{i,j} = 0 for *i* ≠ 1 or *j* ≠ 1; *ρ* is the single scattering albedo of the reflective surface. *μ*′ = cos(*θ*′) and sin(*θ*) = *n _{w}* sin(

*θ*′) in Eqs. (8c) and (8e). The factors of

*n*

_{w}^{2}in Eqs.(8c) and (8e) are due to the n-square law of the radiance transmission through a dielectric interface [4].

The next step is to expand the radiance vectors and the phase matrices into Fourier series of the azimuth. The expansion of the radiance vectors can be written as:

with

$${\mathbf{L}}_{n,\mathrm{sin}}^{m}(\tau ,\mu )={({\mathrm{0,0},U}_{n}^{m},{V}_{n}^{m})}^{T}.$$

The expansion of the phase matrices can be written as:

$$\phantom{\rule{6.0em}{0ex}}=\sum _{m=0}^{M}\left(2-{\delta}_{0m}\right)\mathrm{cos}\left[m\left(\varphi -\varphi \prime \right)\right]{\mathbf{P}}_{\mathrm{cos}}^{m}(\tau ,\mu ,\mu \prime ),$$

$$\phantom{\rule{6.0em}{0ex}}=\sum _{m=0}^{M}\left(2-{\delta}_{0m}\right)\mathrm{sin}\left[m\left(\varphi -\varphi \prime \right)\right]{\mathbf{P}}_{\mathrm{sin}}^{m}(\tau ,\mu ,\mu \prime )$$

In the pioneer works [72, 73, 74], the analytical expressions of **P**
* ^{m}_{sin}* and

**P**

*have been found in terms of the generalized spherical functions (GSF) [75]. Lenoble et al. have recently used the so-called generalized Legendre functions as the expansion basis, which are real by definitions [47]. Not trying to create our own convention, we find it is more convenient to use the Wigner*

^{m}_{cos}*d*functions [76] as the expansion basis. The link between the generalized spherical functions and the Wigner

*d*functions has already been suggested by Hovenier and van de Mee [74]. Previously, the Wigner

*d*functions have been used as the basis to expand the scattering matrices by Mishchenko et al. [77]. The Wigner

*d*functions are also real and better documented than the generalized Legendre functions. More importantly, they can be calculated via stable up recurrence relations. Interested readers may refer to Appendix B of [77] for detailed properties of the Wigner

*d*functions. The derivation of

**P**

*and*

^{m}_{sin}**P**

*in terms of the Wigner*

^{m}_{cos}*d*functions is given in the Appendix A.

The standard procedure includes substituting Eqs.(9) and (10) into Eqs. (7) and (8), carrying out the integration over *ϕ*′, and identifying the same terms containing cos[*m*(*ϕ* - *ϕ*
_{0})] and sin[*m*(*ϕ* - *ϕ*
_{0})] at both sides of the equations. The dependence of *ϕ* is then reduced to a set of equations with variables of *τ* and *μ* and order number *m*. The equations with different order *m* are completely decoupled and can be solved separately. For details, please refer to [47, 73, 74, 78, 79]. After lengthy manipulations, the resultant formulas are:

where

$$\phantom{\rule{4.5em}{0ex}}+{e}^{\left(2{\tau}_{a}^{*}-\tau \right)/{\mu}_{0}}{\mathbf{P}}^{m}(\tau ,\mu ,{-\mu}_{0})\mathbf{r}\left(\pi -{\theta}_{0}\right)\}{\mathbf{E}}_{0},$$

and

In the previous equations, **L**
* ^{m}_{n}* =

**L**

^{m}

_{n,cos}+

**L**

*m*

_{n,sin}and

**P**

^{m}=

**P**

^{m}

_{cos}+

**P**

^{m}

_{sin}

**D**with

**D**=

*diag*{1,1,-1.-1}. The purpose of

**D**is to flip signs for the upper right 2 × 2 sub-matrix of

**P**

^{m}

_{sin}[47, 73, 74].

**r**

^{*m}has been treated similarly as

**P**

^{m}. For the Lambertian surface, its (1,1) element is reduced to a constant, i.e., the single scattering albedo; other elements are zero. The dummy integration variable

*μ*″ is used for later convenience.

Equations (11), (12), and (13) are the equations to be used in the SOS method. From a mathematical point of view, the radiances with the factor of *n*
^{2}
_{w} in Eqs. (13) need to be treated with extreme care. At the ocean side, these terms are the contribution transmitted from the atmosphere, confined within the Fresnel cone. Thus extra integration quadrature points are generally needed to cover the region outside the Fresnel cone, which increases both the memory requirement and the CPU time for the computation. In the next section, a single quadrature scheme is presented based on the separation of directly transmitted sky light and the diffuse light in ocean. The double quadrature scheme is also implemented in the SOS code. Numerical comparison is made in section 6.

## 4. The boundary condition at the ocean surface

The starting point is to evaluate the source function Eq.(12c) with the boundary conditions of Eqs. (13). In the following description, we assume that the source function in the ocean is to be found. Substituting Eq.(13e) in Eq. (11a) leads to:

where the first term in Eq.(14a) is the direct contribution from photons transmitted into the ocean from the atmosphere. This part of the radiance field is confined within the Fresnel cone. The total radiance field as shown in Eq. (14a) is not continuous because the first term jumps from a finite value at the edge of the Fresnel cone to zero just beyond the Fresnel cone. In the case of wind roughed ocean surface, this discontinuity is partly smoothed out but the rate of change of the radiance field near the critical angle will still be large. The term **L**
* ^{m}_{n,diff}* (

*τ*,

*μ*) represents the diffuse radiance field. The factor

**r**(

*π*-

*θ*) also changes abruptly around the edge of the Fresnel cone because of the total internal reflection for incident angles larger than the critical angle. Both these two parts contribute to the discontinuity around the critical angle. The quadrature techniques, however, requires the function to be smooth in order to get an accurate integration. This property has prompted people to use two sets of quadrature points to do the integration in Eq. (12c). One set of

*μ*covers the Fresnel cone and the other covers the rest of the region.

Further analysis finds that the integration of the first term of Eq. (14a) can be done by transforming the integral variables *μ*″, which is below the ocean surface, to the variable *μ*′, which is above the ocean surface. Specifically, sin(*π* - *θ*′) = *n _{w}* sin(

*π*-

*θ*″) and cos(

*θ*′)

*dθ*′ =

*n*cos(

_{w}*θ*″)

*dθ*″ lead to the following equation:

$$\phantom{\rule{2.0em}{0ex}}=-{n}_{w}\mathrm{sin}\left({\theta \prime}^{}\right)\frac{\mathrm{cos}(\theta \prime )d\theta \prime}{{n}_{w}\mathrm{cos}(\theta ")}$$

$$\phantom{\rule{2.0em}{0ex}}=\frac{{\mu}^{\prime}}{{\mu}^{"}}d{\mu}^{\prime}.$$

This relation is substituted into the source term caused by the direct sky light transmission from the atmosphere:

$$\phantom{\rule{4.8em}{0ex}}=\frac{\omega \left(\tau \right)}{2}{\int}_{-1}^{0}{\mathbf{P}}^{m}(\tau ,{\mu}_{p},\mu "){\mathbf{t}(\pi -\theta \prime )\mathbf{L}}_{n-1}^{m}({\tau}_{a}^{*},\mu \prime ){e}^{-\left({\tau}_{l}-\tau \right)/\mu "}\frac{\mu \prime}{\mu "}d\mu \text{\'}\text{}$$

$$\phantom{\rule{4.8em}{0ex}}=\frac{\omega \left(\tau \right)}{2}\sum _{q=1}^{P/2}{w}_{q}{\mathbf{P}}^{m}(\tau ,{\mu}_{p},{\mu}_{q}"){\mathbf{t}(\pi -{\theta}_{q}\prime )\mathbf{L}}_{n-1}^{m}({\tau}_{a}^{*},{\mu}_{q}\prime ){e}^{-\left({\tau}_{l}-\tau \right)/{\mu}_{q}"}\frac{{\mu}_{q}\prime}{{\mu}_{q}"},$$

where *μ _{c}* = cos(

*θ*) and

_{c}*θ*=

_{c}*π*- arcsin(1/

*n*) is the critical angle which defines the edge of Fresnel cone;

_{w}*μ*′ and

_{q}*w*are the quadrature points and weights in the atmosphere, respectively. The summation in the last step is over half of the quadrature points, which covers the range of -1 ≤

_{q}*μ*′ ≤ 0. In the simulation, only

**P**

^{m}(

*τ*,

*μ*,

_{p}*μ*) will be pre-calculated and stored in memory, where

_{q}*μ*

_{p,q}are the quadrature points. In the single quadrature scheme, this means

**P**

^{m}(

*τ*,

*μ*p,

*μ*″

_{q}) will not be found because

*μ*″

_{q}is not one of the quadrature points in the ocean. The solution is simply to use linear interpolation between the two nearby quadrature points. Now it is assumed that the discontinuity caused by

**r**(

*π*-

*θ*) in the diffuse light is insignificant compared to the total radiance field. Thus the integration of the diffuse light in Eq. (14b) can be done using a single quadrature scheme:

$$\phantom{\rule{4.5em}{0ex}}=\frac{\omega \left(\tau \right)}{2}\sum _{q=1}^{p}{w}_{q}{\mathbf{P}}^{m}(\tau ,{\mu}_{p},{\mu}_{q}){\mathbf{L}}_{n-1,\mathrm{diff}}^{m}(\tau ,{\mu}_{q}).$$

where *μ*
_{p,q} and *w _{q}* are the quadrature points and the associated weights;

*P*is the order of the quadrature scheme. The total source term is:

Equations (16) and (17) only use one set of Gaussian quadrature points in the ocean, which leads to less memory usage and CPU time. By assuming insignificant contribution of **r**(*π* - *θ*) in the diffuse light, the question remains as to how large the error will be. We have also implemented the double quadrature scheme in the ocean. Numerical comparison shows that the
difference between the two schemes is negligible. The source terms by the direct transmitted light from the atmosphere is taken care by Eq. (16). If the radiance field in the ocean needs to be exported at the end of the simulation, a detector array is defined to record the first term Eq.(14) separately. The separation of the direct and diffuse contribution to the detector offers another benefit. If only the second term in Eq.(13c) is substituted in Eq.(11b), and the optical depth integration term is ignored, the resultant equation is the water leaving radiance:

where the subscript *wlr* stands for the water leaving radiance. The water leaving radiance is important for the retrieval of the inherent optical properties (IOP) of ocean waters. Thus an efficient and accurate model for water leaving radiance simulation is very useful. The paper by Chami et al. has tried to estimate the water leaving radiance by removing the sky light reflected from the ocean surface [40]. That procedure is inefficient and less accurate in the sense that it ignores the multiple interaction between the diffuse light in the atmosphere and ocean. Equation (19), however, provides an exact solution for the water leaving radiance without causing any extra computational burden.

## 5. Numerical procedure

The atmosphere and ocean are discretized into *K _{a}* and

*K*homogeneous layers, respectively. The optical thickness of each layer needs to be small enough to meet the required numerical accuracy. For the atmosphere, the optical depth at each layer is denoted as

_{o}*τ*with

_{k}*k*= 0,1,2,…,

*K*.

_{a}*τ*

_{0}= 0,

*τ*

_{1}=

*τ*

_{0}+

*δτ*,…,

*τ*

*k*=

*τ*

_{k-1}+

*δτ*,…

*τ*

*K*=

_{a}*τ*

^{*}

_{a}. This means that the

*k*th layer is between

*τ*

_{k-1}and

*τ*

_{k}. Practically,

*δτ*does not have to be equal for each layer and can be adjusted for convenience. For the ocean,

*k*=

*K*+1,

_{a}*K*+2, …,

_{a}*K*+

_{a}*K*+1 with

_{o}*τ*

_{Ka+1}=

*τ*

_{o},

*τ*

_{Ka+2}=

*τ*

_{Ka+1}+

*δτ*, … ,

*τ*

_{Ka+Ko+1}=

*τ*

^{*}

_{o}.

The zenith angle cosine *μ* is discretized into the Gaussian quadrature of even order *P*. *μ _{p}* and

*w*are the quadrature points and weights with

_{p}*p*= 1,2, …,

*P*, respectively.

**P**

^{m}(

*τ*,

*μ*,

_{p}*μ*) for each

_{q}*p*and

*q*are calculated and stored in memory before the SOS iteration. Even

**P**

^{m}is a function of

*τ*while many layers with different

*τ*may share the same

**P**

^{m}. Thus an index is created for each layer. The layers with the same scattering properties are labeled with the same index numbers. These index numbers will be used to refer to the correct

**P**

^{m}when needed.

The source functions **S**
^{m}
_{1} (*τ*
_{k},*μ _{p}*) for the atmosphere and ocean are initialized according to Eqs.(12a) and (12b), respectively. The SOS iteration begins with Eqs. (11) with known

**S**

^{m}

_{1}(

*τ*k,

*μ*). At this point, the first terms in Eqs.(11) are neglected temporarily. The integration for

_{p}*μ*< 0 is done in the following way:

$$\phantom{\rule{4.3em}{0ex}}=-{\int}_{{\tau}_{l}}^{{\tau}_{k-1}}{e}^{-(\tau \prime -{\tau}_{k})/{\mu}_{p}}{\mathbf{S}}_{n}^{m}(\tau \prime ,{\mu}_{p})\mathrm{d\tau}\prime /{\mu}_{p}+{\int}_{{\tau}_{k-1}}^{{\tau}_{k}}{e}^{-(\tau \prime -{\tau}_{k})/{\mu}_{p}}{\mathbf{S}}_{n}^{m}(\tau \prime ,{\mu}_{p})\mathrm{d\tau}\prime /{\mu}_{p},$$

$$\phantom{\rule{4.3em}{0ex}}={e}^{\mathrm{\delta \tau}/{\mu}_{p}}{\mathbf{L}}_{n}^{m}({\tau}_{k-1},{\mu}_{p}<0)-{\int}_{{\tau}_{k-1}}^{{\tau}_{k}}{e}^{-(\tau \prime -{\tau}_{k})/{\mu}_{p}}{\mathbf{S}}_{n}^{m}(\tau \prime ,{\mu}_{p})\mathrm{d\tau}\prime /{\mu}_{p}.$$

The exponential - linear approximation has been implemented to evaluate the integration between *τ*
_{k-1} and *τ*k [80]. Equation (20) shows that **L**
* ^{m}_{n}*(

*τ*,

_{k}*μ*< 0) can be updated successively from

_{p}*k*= 0 to

*K*for the atmosphere and from

_{a}*k*=

*K*+ 1 to

_{a}*K*+

_{a}*K*+ 1 for the ocean. For

_{o}*μ*> 0, the procedure is:

_{p}$$\phantom{\rule{4.3em}{0ex}}={\int}_{{\tau}_{k+1}}^{{\tau}_{u}}{e}^{-(\tau \prime -{\tau}_{k})/{\mu}_{p}}{\mathbf{S}}_{n}^{m}(\tau \prime ,{\mu}_{p})\mathrm{d\tau}\prime /{\mu}_{p}+{\int}_{{\tau}_{k}}^{{\tau}_{k+1}}{e}^{-(\tau \prime -{\tau}_{k})/{\mu}_{p}}{\mathbf{S}}_{n}^{m}(\tau \prime ,{\mu}_{p})\mathrm{d\tau}\prime /{\mu}_{p},$$

$$\phantom{\rule{4.3em}{0ex}}={e}^{\mathrm{-\delta \tau}/{\mu}_{p}}{\mathbf{L}}_{n}^{m}({\tau}_{k+1},{\mu}_{p}>0){+\int}_{{\tau}_{k}}^{{\tau}_{k+1}}{e}^{-(\tau \prime -{\tau}_{k})/{\mu}_{p}}{\mathbf{S}}_{n}^{m}(\tau \prime ,{\mu}_{p})\mathrm{d\tau}\prime /{\mu}_{p}.$$

**L**
* ^{m}_{n}* (

*τ*k,

*μ*> 0) will be updated successively from

_{p}*τ*

_{Ka+Ko+1}to

*τ*

_{Ka+1}for the ocean and from

*τ*

_{Ka}to

*τ*

_{0}for the atmosphere. Again the exponential - linear approximation has been implemented to evaluate the integration between

*τ*and

_{k}*τ*+ 1 [80]. After this step, all of the source terms should be reset to zero to prepare the source terms for the next order of scattering.

_{k}The next step is to consider the boundary conditions. The bottom reflection can be taken into account by adding Eq. (13f) to **L**
_{mn} (*τ _{k}*,

*μ*> 0) according to Eq. (11b) for the ocean. For the interface, only the first terms in Eqs. (13c) and (13e) are added to

_{p}**L**

*(*

^{m}_{n}*τ*

_{k},

*μ*) for the purpose of source integration in Eqs.(17). That is:

_{p}where ⇐ is understood to evaluate the right hand side and give the resultant value to the left hand side; *θ _{p}* = arccos(

*μ*). After Eqs. (22), the single quadrature scheme uses Eqs.(17) to get the source function due to ocean diffuse light. The double quadrature scheme, however, integrates the source function in two separate regions, -

_{p}*μ*<

_{c}*μ*<

*μ*, and |

_{c}*μ*| >

*μ*. The integration of the second terms in Eqs. (13c) and (13e) involving

_{c}*n*

^{2}

_{w}in both schemes is done using Eq.(16) or its equivalent form for the transmission from the ocean to the atmosphere. The only difference is that the double quadrature scheme does not involve the interpolation of

**P**

^{m}, because

*μ*″

_{q}in Eq.(16) is one of the predefined quadrature points in the ocean. A separately allocated array

**L**

*(*

^{m}_{n,det}*τ*,

_{k}*μ*) is used to store the radiance field for the final output:

_{p}$$\phantom{\rule{7.2em}{0ex}}+{e}^{-\left({\tau}_{o}-{\tau}_{k}\right)/{\mu}_{p}}{n}_{w}^{2}\mathbf{t}\left(\pi -{\theta}_{p}^{\prime}\right){\mathbf{L}}_{n-1}^{m}({\tau}_{a}^{*},{\mu}_{p}^{\prime}),$$

$$\phantom{\rule{7.2em}{0ex}}+\frac{{e}^{\left({\tau}_{k}-{\tau}_{a}^{*}\right)/{\mu}_{p}}}{{n}_{w}^{2}}\mathbf{t}\left(\pi -{\theta}_{p}^{\prime}\right){\mathbf{L}}_{n-1}^{m}({\tau}_{o},{\mu}_{p}^{\prime}).$$

Equations (23) are done only at levels where the output is desired. Again, in the single quadrature scheme **L**
^{m}
_{n} (*τ _{k}*,

*μ*′) is not directly available because

_{p}*μ*′ is not one of the quadrature points. The solution is also to interpolate the values at the two near points. This is also physically sound because

_{p}**L**

*(*

^{m}_{n}*τ*,

_{k}*μ*) is generally smooth before interacting with the ocean interface.

_{p}In addition, Eq. (13g) needs to be evaluated and stored to prepare the next order of scattering. After this step, **L**
^{m}
_{n} (*τ _{k}*,

*μ*) are set to zero. The steps from Eq.(20) to here are repeated

_{p}*N*times, where

*N*is a preset value for the maximum order of scattering based on sensitivity studies. The

*N*could also be determined based on convergence criterias during the run time [40]. However, it is more attractive to use a preset value in terms of efficiency. After the iteration of the procedure for

*N*times, the remainder is taken into account by a geometric series:

Duan and Min have shown that this relation is satisfied with *N* ≈ 40 for most, if not all plane parallel systems [42]. This requirement is smaller for turbid media with small asymmetry factors
or single scattering albedos. This relation is assumed for all the four Stokes parameters unless *L*
^{m}
_{n-1}(*τ _{k}*,

*μ*) = 0 for one of the elements

_{p}*Q*,

*U*, and

*V*. Then the total radiance field is given by:

## 6. Validation and results

The validation of the SOS code is done with the Monte Carlo code by Zhai et al. [65]. This Monte Carlo code uses biased techniques and estimation schemes to increase the efficiency and accuracy. It can deal with both one-dimensional and three dimensional vector radiative transfer problems. Even it is well known that the Monte Carlo method converges slow, the error associated with the Monte Carlo method is easy to understand, i.e. ∝ *N*
^{-½}
_{photon}, where *N _{photon}* is the number of the photon bundles. As long as

*N*is large enough, the code can produce reliable results. In the validation case, the atmosphere is chosen to be a conservative Rayleigh medium with optical depth of

_{photon}*τ*= 0.1. The Rayleigh scattering matrix is:

_{a}The ocean phase function is a Henyey-Greenstein (HG) scattering phase function [81]:

where *μ* = cos(*θ*) and *g* is the asymmetry parameter defined by:

The range of the asymmetry parameter is [-1,1]. *g* ≈ 1 sugguests the scattering function is strongly peaked in the forward direction. *g* ≈ 0 suggests the scattering function is nearly isotropic. In our calculation, we take *g* = 0.9185, which means the scattering function has a large forward peak. The HG phase function with *g* = 0.9185 has the same integrated backscattering coefficients *b _{b}* as the Petzold phase function [82]. The other scattering matrix elements for the ocean is determined by:

where the superscript *HG* shows the matrix elements are with the HG scattering phase function while *RAY* shows that they are with the Rayleigh scattering phase function. In other words, the reduced scattering matrix elements of the ocean water are the same as those of the Rayleight scattering matrix [83]. *δ* - fit technique for the phase function truncation has been extended to the full scattering matrix to fulfill the polarization calculation requirement in this code. The details are provided in Appendix B. The optical depth is 10.0 for the ocean. The single scattering albedo in ocean depends on the wavelength with a possible range of highly absorptive (*ω* ≈ 0.1) to nearly conservative (*ω* ≈ 0.9). The single scattering albedo has been chosen to be 0.5 as a representative value for the primary study in this paper. The bottom is a Lambertian surface with single scattering albedo of 0.5. The refractive index of the ocean water is 1.338.

Figures 1 and 2 show the Stokes parameters *I*, *Q*/*I*,*U*/*I*, and *V*/*I* as functions of observation angles for the system at the top and bottom of the atmosphere, respectively. The observation angle is 0° for the upward (from the bottom of the ocean to the top of the atmosphere (TOA) direction and 180° for the downward (from the TOA to the bottom of the ocean) direction. The principle plane of the observation angles shown in the figures is obtained by rotating the solar incident plane by an azimuthal angle of 90° clock-wise if looking in the zenith direction. The solid lines are calculated by the Monte Carlo method. The circle and cross symbols are computed by the SOS method with the single Gaussian (SG) quadrature and double Gaussian (DG) quadrature schemes, respectively. A solar source with an incident angle of 120° is used. The principle plane for the observation angles is at an azimuthal angle of 90° from the solar incident plane. For the Monte Carlo code, the total photon number used is 5000000. The Gaussian quadrature order of *P _{a}* = 36 is used in the atmosphere. The ocean Gaussian quadrature order of

*P*= 72 is used for the DG case. The total number of scattering orders is

_{o}*N*= 10.

*δτ*= 0.01 and 0.2 have been used for the optical depth integration of the ocean.

*δτ*= 0.01 has been used in the atmosphere if

*δτ*= 0.01 in the ocean.

*δτ*= 0.2, however, cannot be used for the optical depth integration of the atmosphere because it is larger than 0.1, the optical depth of the atmosphere. Thus the value

*δτ*= 0.05 is selected for the atmosphere in this case, which is the largest value one can select because there have to be at least two layers (

*K*= 3) in the atmosphere to solve for the three unknown coefficients in the exponential - linear approximation. The order of the Fourier expansion used is

_{a}*M*= 8 in Eq. (9). The maximum order of the HG phase function expansion is

*L*= 36 in Eq.(A-13). The

*δ*- fit technique is generalized to the full scattering matrix for the polynomial expansion. The SOS code was run using an I-Mac with a 2.8 Ghz Intel processor. The fortran compiler is GFORTRAN with the -

*O*3 compile option. Table 1 shows the average CPU times for both the SG and DG schemes for the cases shown in Figs. 1 and 2.

Both the SG and DG quadrature schemes results agree well with the Monte Carlo results for all the four Stokes parameters at both the top and bottom of the atmosphere. Note that only the Stokes parameters within observation angles of 0-90° are shown at the TOA because they are zeros for observation angles larger than 90°. For both *δτ* = 0.01 and *δτ* = 0.2, the differences
between the SG and DG are negligible, for the TOA and bottom of the atmosphere. For the TOA, the errors of the SG and DG scheme radiances with *δτ* =0.2 vary between 0.25% at 0° to 1.0% at 82.6°, if using the DG scheme with *δτ* =0.01 as a standard. The error at 87.5° is -2.5%, which is slightly larger than smaller angles. The linear polarization components *Q*/*I* and *U*/*I* also have larger errors at larger observation angles. The differences of the circular polarization components *V*/*I* between the SOS code and Monte Carlo code are larger compared to those of *I*, *Q*/*I*, and *U*/*I*. The DG scheme has slightly better agreement with the Monte Carlo method. For the bottom of the atmosphere, the differences between the SG, DG, and the Monte Carlo results are smaller than those of the TOA. The errors of the SG and DG cases with *δτ* = 0.2 are from -0.5% to 0.5%. Another interesting feature is that the ratios *Q*/*I* and *U*/*I* have smaller errors than the absolute radiance. Considering the large forward peak of the phase function, the results in Figs.1 and 2 are quite impressive.

Figures 3 show the water leaving radiance for the same system as in Fig.1 at the bottom of the atmosphere (just above the ocean surface). The radiance percentage differences between
the SOS results and the Monte Carlo results are around 2%, which are larger than those of the total radiance plots. The reason is that the ocean has a large asymmetry factor *g* = 0.9185 and an optical depth of *τ* = 10, which is rather difficult for the Monte Carlo method to achieve convergence. The SOS method, however, virtually takes all order of scattering into account (see Eq. (25)). The better agreements between the two methods for the total radiance is due to the contribution of the atmospheric Rayleigh layer, which is thin and the phase function is symmetric around the scattering angle of 90°. It is obvious that the polarization components *Q*/*I*, *U*/*I*, and *V*/*I* agree well between the two methods. Another very interesting feature is that the polarized components converge before the absolute radiance for both methods. The small value of the circular component makes it more sensitive to the optical depth integration step *δτ*. Again, the difference between the SG and DG results are negligible. The water leaving radiance for *δτ* = 0.2 has relative errors from -0.6% to -1.5%, if comparing to the results with *δτ* = 0.01.

The next case to be studied is designed to answer the following question. For a clear sky condition with an aerosol optical depth up to 0.3, how many orders of scattering is needed to achieve a radiance error smaller than 0.5% in an AOS- For this purpose, a layer of aerosol is added at the bottom of the atmosphere in the first case studied. The aerosol consists in a fine mode of spherical particles with log-normal size distribution [84]. The modal radius is *r _{m}* =0.08

*μm*and the variance is σ =0.46. The refractive index is

*m*=1.45. The zero imaginary part of the refractive index leads to a single scattering albedo of 1. The incident wavelength is 0.443

*μm*. The conservative Rayleigh layer of optical depth 0.1 is unchanged, located at the top of the aerosol layer. The aerosol layer at the atmosphere bottom does not mix with the Rayleigh layer. To make sure the standard converged results are obtained before the sensitivity study, an optical depth of

*τ*= 0.3 is given to the aerosol layer. The ocean is also kept the same as the previous case. Figure 4 shows the effects of different Gaussian quadrature orders

_{sol}*P*and Fourier series orders

*M*. The radiance and

*Q*/

*I*shown in the figure is at the TOA. A total scattering order of

*N*= 40 is used for all the SOS calculation. The azimuthal angle for Fig. 4 is 180°, because the Stokes parameter at this azimuthal angle contains minimum amount of sun glint, which is very important in the satellite remote sensing. The Gaussian order

*P*shown in Fig.4 refers to the Gaussian quadrature order used in the atmosphere. For the DG case, the corresponding Gaussian order used in the ocean is 2

*P*. In all the cases, the differences between the SG and DG case are negligible. It is observed that the results with the Gaussian quadrature order of

*P*= 18 and the Fourier series order of

*M*= 8 are already converged to the results with

*P*= 36 and

*M*= 36. The Gaussian order

*P*= 8 produces larger errors because

*P*= 8 is much smaller than the ocean Legender expansion order of

*L*= 36. The following sensitivity studies of scattering

*N*use the DG scheme with

*P*= 36 and

*M*= 8. The results with the scattering orders of

*N*= 40 are regarded as the true solution.

Table 2 shows the order of scattering needed for convergence. A set of aerosol optical depth *τ*sol are used to test the sensitivity of scattering orders *N*. Other system configurations are the same as the previous cases. The radiance at the TOA are studied and the observation azimuthal angle is 180°. Both the total radiance and water leaving radiance are tested. A tolerance level of ±0.5% is used in the calculation of Table 2. It is shown that scattering orders of *N* = 4 to *N* = 8 are needed to achieve errors smaller than ±0.5% for the total radiance. The scattering orders are from *N* = 8 to *N* = 10 are needed for the water leaving radiance calculation. The requirement for scattering orders *N* is higher for the water leaving radiance than the total radiance. The reason is that the total radiance consists of the contributions from both the atmosphere and the ocean. The atmosphere is a thin conservative medium in which a few order of scattering are enough for convergence. Thus the total radiance is easier to achieve the ±0.5% error tolerance level. However, the water leaving radiance comes only from the ocean. The ocean is a optically thick medium, which leads to more orders of scattering. The results with the scattering order
of *N* = 10 have errors smaller than 0.1% for all the total radiance calculation. The maximum errors with *N* = 10 associated with the water leaving radiance are around 0.25%.

It should be emphasize that the SOS method will converge even for optically thick and near conservative media given a large enough order of scattering. Min and Duan have presented that the SOS method can converge for a turbid medium of optical depth up to 50, with a single scattering albedo of 0.99 [42]. The limitations of the SOS code may be its computational speed, which can only be found by crosscheck with other methods. This will be an interesting work in the future. The last case we use 0.5% as a standard to study the how many orders of scattering are needed for converging a realistic atmosphere-ocean system. Applications of radiative transfer models for aerosol retrievals may require the polarized radiance measurement as accurate as 0.1% [85, 86]. As stated in the previous paragraph, the total polarized radiance converges within 0.1% at the TOA with the scattering order of *N* = 10 for the AOS studied. A higher single scattering albedo may lead to larger scattering orders.

## 7. Conclusion

This paper presents a vector radiative transfer model for an AOS. The model is based on the SOS Method. The SOS formulas have apparent physical significances and the errors are easy to be analyzed. The source terms in the RTE involves numerical integration with Gaussian quadratures. Previously published models normally use two sets of Gaussian quadrature points in ocean and one set in atmosphere to treat the interface boundary condition correctly. This paper has presented a Single Gaussian (SG) quadrature scheme in ocean in comparison with the Double Gaussian (DG) quadrature scheme. Although two sets of Gaussian quadratures is essential for discrete ordinate methods, the differences between the SG and DG schemes can be small for SOS in a typical atmosphere and ocean configuration. Both memory requirements and CPU time are smaller for the SG scheme, which can be helpful for remote sensing applications. The current SOS model (both SG and DG schemes) can differentiate water leaving radiance from light scattered by the atmosphere and by sea surface reflection of the sky light. The model results compare well against the Monte Carlo simulations. For a realistic case involving an ocean of optical depth *τ* = 10 and single scattering albedo *ω* = 0.5, the SOS converges with the order of scattering *N* = 10 and the Gaussian quadrature order *P* = 36. The *δ* - fit technique is generalized to the whole 4×4 scattering matrix to fulfill the polarization calculation requirement. This work is primarily of interest in aerosol and ocean color remote sensing, and other relevant scientific disciplines. Both the SG and DG scheme computing packages are available to the public science community by contacting the corresponding author.

## Appendix A: the expression of P^{m}

First a few properties of the Wigner *d* functions are rephrased here for convenience. The generalized spherical functions are related to the Wigner *d* functions in a way of [74, 75, 77]:

where ξ = cos(*θ*). In the case of **n** = ±2,

Substituting Eq. (A-2) in Eqs. (14a) and (14b) of [73] leads to:

with the notations of:

If *m* = *n* = 0, the Wigner *d* functions are reduced to the Legendre polynomials,

For *n* = 0, the Wigner *d* functions are related to the associated Legendre function as follows:

Setting *m* = 2 in Eq.(A-6) leads to

For the real derivation, the starting point is the scattering matrix: (Eq. (3) in [73]):

The elements *a*
_{1}, *a*
_{2}, *a*
_{3}, *a*
_{4}, *b*
_{1}, and *b*
_{2} are expanded in terms of the generalized spherical
functions. Our goal is to express everything in terms of the Wigner *d* functions. Substituting Eq. (A-5) in Eqs. (19a) and (19b) of [73], the coefficients *β _{l}* and

*δ*are:

_{l}The coefficients *γ _{l}* and

*ϵ*are found by substituting Eq. (A-7) into Eqs. (19c) and (19d) of [73]:

_{l}The coefficients *α _{l}* and

*ζ*can be obtained by setting

_{l}*m*= 2 in Eq. (A-3) and plugging in Eqs. (23a) and (23b) of [73]:

Note that these coefficients are consistent with those given in Page 103-104 of [77] even different notations have been used. The relations are:

where the order number *s* in [77] has been replaced with *l*; and the subscripts *M* are added to the coefficients in [77] for clarification.

In the visible range, scattering particles are often much larger than the incident wavelength. In turn, the scattering phase function *a*
_{1} usually has a strong forward peak, which leads to hundreds or thousands of expansion terms *β _{l}*. In this case the truncation techniques are normally
used to reduce the number of fitting coefficients significantly [87, 88]. In this paper the

*δ*- fit technique is generalized to fit the scattering matrix [88], which is presented in Appendix B. If no truncation is needed, the method in [77] is used to calculate the expansion coefficients as in (A-12).

Equations (10), which are also used in [47], are different from Eq. (6) of [73] by a factor
of 2. Careful check finds that *C ^{m}* = 2

*P*

^{m}_{cos}and

*S*= 2

^{m}*P*

^{m}_{sin}, where

*C*and

^{m}*S*are the notations used in [73]. These two factors of 2 cancel out and the whole formulas remain consistent. Now define

^{m}**P**

^{m}=

**P**

^{m}

_{cos}+

**P**

^{m}

_{sin}

**D**with

**D**=

*diag*{1,1,-1.-1}. The expression of

**P**

^{m}can be found by substituting Eqs. (A-6) and (A-3) into Eqs. (7a) and (7b) of [73], with the knowledge of which a factor of 2 has been canceled out:

where *L* is the maximum order of *l* used for the expansion; *d ^{l}*

_{m0}=

*d*

^{l}_{m0}(

*μ*),

*d*

^{l,+}

_{m2}=

*d*

^{l,+}

_{m2}(

*μ*), and

*d*

^{l,-}

_{m2}=

*d*

^{l,-}

_{m2}(

*μ*); similarly,

*d*′

^{l}

_{m0}=

*d*

^{l}_{m0}(

*μ*′),

*d*′

^{l,+}

_{m2}=

*d*

^{l,+}

_{m2}(

*μ*′), and

*d*′

^{l,-}

_{m2}=

*d*

^{l,-}

_{m2}(

*μ*′).

**P**
^{m} can be precalculated and stored in the computer memory before the SOS iteration, which will speed up the code significantly. The computation aspects and the recurrence relations of theWigner *d* functions are documented in Appendix B of [77], which will not be repeated here.

## Appendix B: generalized δ - fit for the expansion of the scattering matrix

Both the Delta-M [87] and the *δ* - fit [88] techniques approximate the strong forward peak of the scattering phase function by a delta function.

The difference between the two techniques is how to determine the truncation factor *f*. The Delta M method sets *f* = *β*
_{L+1}, and rescales the coefficients by

This selection is accurate for irradiance simulation. However, it causes large errors in the radiance simulation [89]. Therefore, *δ* - fit uses a different strategy, which is to find the expansion
coefficients that minimize the relative difference *ϵ* between the approximate phase function *a*
^{*}
_{1} (*θ*
_{k}) and the actual phase function *a*
_{1}(*θ _{k}*):

in which

*θ _{k}* is the scattering angle,

*w*is the weight for each scattering angle. By setting the weights for the forward-scattering angles to zero, the forward peak will be automatically truncated. The singular-value decomposition method is used to solve the equations ∂

_{k}*ϵ*/∂

*β*

^{*}

_{l}= 0. After a set of coefficients

*β*

^{*}

_{l}are found, the renormalization factor is

*f*= 1-

*β*

^{*}

_{0}. All the coefficients are divided by

*β*

^{*}

_{0}to renormalize the fitted phase function. The

*δ*- fit is proven to be accurate for radiance simulations [88].

This work deals with the polarized radiative transfer equation. Thus the δ - fit needs to be generalized to all the scattering matrix elements. To conserve the polarization properties of the radiance field, it is important to keep the reduced scattering matrix elements unchanged. Therefore the following procedure is proposed.

First, the expansion coefficients for the phase function are found by using the least squares fitting. The coefficients are substituted into Eq.(B-4) to produce the fitted phase function. The fitted scattering matrix elements are then obtained by:

where *F _{ij}*(

*θ*) are the scattering matrix elements other than

_{k}*a*

_{1}(

*θ*) =

_{k}*F*

_{11}(

*θ*) in Eq. (A-8). For instance,

_{k}*F*

_{12}(

*θ*k) =

*b*

_{1}(

*θ*),

_{k}*F*

_{22}(

*θ*) =

_{k}*a*

_{2}(

*θ*), …, etc. Finally the expansion coefficients

_{k}*c*

^{l*}

_{ij}for the elements

*F*

^{*}

_{ij}(

*θ*) are found by solving:

For *mn* = 00 and (*i*, *j*) = (4,4), *c*
^{l*}
_{ij} correspond to the coefficients *δ _{l}* in Eq.(A-9b). For

*mn*= 20 and (

*i*,

*j*) = (1,2), (4,3),

*c*

^{l*}

_{ij}correspond to

*γ*, and

_{l}*ϵ*in Eq.(A-10b), respectively. For

_{l}*mn*= 22 and (

*i*,

*j*) = (2,2), (3,3)

*c*

^{l*}

_{ij}correspond to

*ζ*and

_{l}*α*in Eq.(A-11b), respectively. Note that the weight is no longer needed because the truncation is already done by the least squares fitting of the phase function. To use the

_{l}*δ*- fit coefficients, the optical depth and single scattering albedo for the layer have to be adjusted similarly as the Delta-M method:

## Acknowledgment

This research was supported by Peng-Wang Zhai’s appointment to the NASA Postdoctoral Program at the NASA Langley Research Center administered by Oak Ridge Associated Universities through a contract with NASA. He also thanks the Science Systems & Applications, Inc. (SSAI) in Hampton, VA for providing the office space and computer support. This study is also supported by Dr. Hal Maring of NASA Radiation Science program and Dr. Paula Bontempi of NASA Biogeochemistry program.

## References and links

**1. **S. Chandrasekhar, *Radiative Transfer*(Dover, New York, 1960).

**2. **R. W. Preisendorfer, *Radiative Transfer on Discrete Spaces*(Pergamon Press, Oxford, 1965).

**3. **H. C. van de Hulst, *Multiple Light Scattering: Tables, Formulas, and Applications***1** and **2** (Academic Press, New York, 1980).

**4. **C. D. Mobley, *Light and Water: Radiative Transfer in Natural Waters*(Academic, San Diego, Calif., 1994).

**5. **K.N. Liou, *An Introduction to Atmospheric Radiation* 2nd edition, (Academic, San Diego2002).

**6. **A. Marshak and A. B. Davis (Eds.), *3D Radiative Transfer in Cloudy Atmospheres*(Springer, 2005). [CrossRef]

**7. **M. I. Mishchenko, L. D. Travis, and A. A. Lacis, *Multiple Scattering of Light by Particles* (Cambridge University Press, New York, 2006).

**8. **K. N. Liou, “A numerical experiment on Chandrasekhar’s discrete-ordinate method for radiative transfer: Application to cloudy and hazy atmospheres,” J. Atmos. Sci. **30**, 1303–1326 (1973). [CrossRef]

**9. **K. Stamnes, S.-C. Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinate method radiative trasfer in multiple scattering and emitting layered media,” Appl. Opt. **27**, 2502–2509 (1988). [CrossRef] [PubMed]

**10. **F. Weng, “A multi-layer discrete-ordinate method for vector radiative transfer in a vertically-inhomogeneous, emitting and scattering atmosphere-I. theory,” J. Quant. Spectrosc. Radiat. Transfer **47**, 19–33 (1992). [CrossRef]

**11. **F. Weng, “A multi-layer discrete-ordinate method for vector radiative transfer in a vertically-inhomogeneous, emitting and scattering atmosphere-II. Application,” J. Quant. Spectrosc. Radiat. Transfer **47**, 35–42 (1992). [CrossRef]

**12. **A. Sánchez, T. F. Smith, and W. F. Krajewski, “A three-dimensional atmospheric radiative transfer model based on the discrete-ordinates method,” Atmos. Research **33**, 283–308 (1994). [CrossRef]

**13. **J. L. Haferman, T. F. Smith, and W. F. Krajewski, “A multi-dimensional discrete-ordinates method for polarized radiative transfer. Part I: validation for randomly oriented axisymmetric particles,” J. Quant. Spectrosc. Radiat. Transfer **58**, 379–398 (1997). [CrossRef]

**14. **F. M. Schulz, K. Stanes, and F. Weng, “VDISORT: An improved and generalized discrete ordinate method for polarized vector radiative transfer,” J. Quant. Spectrosc. Radiat. Transfer **61**, 105–122 (1999). [CrossRef]

**15. **C. E. Siewert, “A discrete-ordinates solution for radiative-transfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transfer **64**, 227–254 (2000). [CrossRef]

**16. **V. A. Ambartzumian, “A new method for computing light scattering in turbid media,” Izv. Akad. Nauk SSSR, Ser. Geogr. Geofiz. **3**, 97–104 (1942).

**17. **V. A. Ambartzumian, *Theoretical Astrophysics*(Pergamon Press, New York, 1958).

**18. **C. N. Adams and G. W. Kattawar, “Solutions of the equation of radiative transfer by an invariant imbedding approach,” J. Quant. Spectrosc. Radiat. Transfer **10**, 341–366 (1970). [CrossRef]

**19. **M. I. Mishchenko, “Reflection of polarized light by plane-parallel slabs containing randomly-oriented, nonspherical particles,” J. Quant. Spectrosc. Radiat. Transfer **46**, 171–181 (1991). [CrossRef]

**20. **M. I. Mishchenko and L. D. Travis, “Satellite retrieval of aerosol properties over the ocean using polarization as well as intensity of reflected sunlight,” J. Geophys. Res. **102**, 16989–17013 (1997). [CrossRef]

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

**22. **J. E. Hansen, “Multiple scattering of polarized light in planetary atmospheres. part II. Sunlight reflected by terrestrial water clouds,” J. Atmos. Sci. **28**, 1400–1426 (1971). [CrossRef]

**23. **G. N. Plass, G. W. Kattawar, and F. E. Catchings, “Matrix operator theory of radiative transfer. 1: Rayleigh-scattering,” Appl. Opt. **12**, 314–329 (1973). [CrossRef] [PubMed]

**24. **G. W. Kattawar, G. N. Plass, and F. E. Catchings, “Matrix operator theory of radiative transfer. 2: scattering from marine haze, Appl. Opt. **12**, 1071–1084 (1973). [CrossRef] [PubMed]

**25. **P. C. Waterman, “Matrix-exponential description of radiative transfer, J. Opt. Soc. Am. **71**, 410–422 (1981).

**26. **T. Nakajima and M. Tanaka, “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). [CrossRef]

**27. **J. Fischer and H. Grassl, “Radiative transfer in an atmosphere-ocean system: an azimuthally dependent matrix operator approach,” Appl. Opt. **23**, 1032–1039 (1984). [CrossRef] [PubMed]

**28. **Q. Liu and E. Ruprecht, “Radiative transfer model: matrix operator method,” Appl. Opt. **35**, 4229–4237 (1996). [CrossRef] [PubMed]

**29. **P. Zhai, G. W. Kattawar, and P. Yang, “Impulse response solution to the three-dimensional vector radiative transfer equation in atmosphere-ocean systems. II. The hybrid matrix operator-Monte Carlo method,” Appl. Opt. **47**, 1063–1071 (2008), http://www.opticsinfobase.org/abstract.cfm?URI=ao-47-8-1063. [CrossRef] [PubMed]

**30. **V. Kourganoff, *Basic Methods in Transfer Problems*(Clarendon Press, London, 1952).

**31. **R. D. M. Garcia and C. E. Siewert, “A generalized spherical harmonics solution for radiative transfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transfer **36**, 401–423 (1986). [CrossRef]

**32. **E. P. Zege, I. L. Katsev, and I. N. Polonsky, “Multicomponent approach to light propagation in clouds and mists,” Appl. Opt. **32**, 2803–2812 (1993). [CrossRef] [PubMed]

**33. **E. P. Zege and L. I. Chaikovskaya, “New approach to the polarized radiative transfer problem,” J. Quant. Spec-trosc. Radiat. Transfer **55**, 19–31 (1996). [CrossRef]

**34. **K. F. Evans, “The spherical harmonic discrete ordinate method for three-dimensional atmospheric radiative transfer,” J. Atmos. Sci. **55**, 429–446 (1998). [CrossRef]

**35. **K. F. Evans, “SHDOMPPDA: A radiative transfer model for cloudy sky data assimilation”. J. Atmos. Sci. **64**, 3858–3868 (2007). [CrossRef]

**36. **R. D. M. Garcia and C. E. Siewert, “The FN method for radiative transfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transfer **41**, 117–145 (1989). [CrossRef]

**37. **W. M. F. Wauben and J. W. Hovenier, “Polarized radiation of an atmosphere containing randomly-oriented spheroids,” J. Quant. Spectrosc. Radiat. Transfer **47**, 491–504 (1992). [CrossRef]

**38. **R.B. Myneni, G. Asrar, and E. T. Kanemasu, “Light scattering in plant canopies: the method of Successive Orders of Scattering Approximations (SOSA),” Agric. For. Meteorol. **39**, 1–12 (1987). [CrossRef]

**39. **J. Lenoble, *Atmospheric radiative transfer*(A. Deepak Publishing1993).

**40. **M. Chami, R. Santer, and E. Dilligeard, “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), http://www.opticsinfobase.org/abstract.cfm?URI=ao-40-15-2398. [CrossRef]

**41. **Q. Min and M. Duan, “A successive order of scattering model for solving vector radiative transfer in the atmosphere,” J. Quant. Spectrosc. Radiat. Transfer **87**, 243–259 (2004). [CrossRef]

**42. **M. Duan and Q. Min, “A semi-analytic technique to speed up successive order of scattering model for optically thick media,” J. Quant. Spectrosc. Radiat. Transfer **95**, 21–32 (2005). [CrossRef]

**43. **T. Greenwald, R. Bennartz, C. O’Dell, and A. Heidinger, “Fast computation of microwave radiances for data assimilation using the Successive Order of Scattering method,” J. Appl. Meteorol. **44**, 960–966 (2005). [CrossRef]

**44. **A. K. Heidinger, C. O’Dell, R. Bennartz, and T. Greenwald, “The successive-order-of-interaction radiative transfer model. part I: model development,” J. Appl. Meteorol. Climatol. **45**, 1388–1402 (2005). [CrossRef]

**45. **S. Y. Kotchenova, E. F. Vermote, R. Matarrese, and F. J. Klemm Jr., “Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part I: Path radiance,” Appl. Opt. **45**, 6762–6774 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=ao-45-26-6762. [CrossRef] [PubMed]

**46. **S. Y. Kotchenova and E. F. Vermote, “Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part II. Homogeneous Lambertian and anisotropic surfaces,” Appl. Opt. **46**, 4455–4464 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=ao-46-20-4455. [CrossRef] [PubMed]

**47. **J. Lenoble, M. Herman, J.L. Deuzé, B. Lafrance, R. Santer, and D. Tanré, “A successive order of scattering code for solving the vector equation of transfer in the earths atmosphere with aerosols,” J. Quant. Spectrosc. Radiat. Transfer **107**, 479–507 (2007). [CrossRef]

**48. **T. Suzuki, T. Nakajima, and M. Tanaka, “Scaling algorithms for the calculation of solar radiative fluxes,” J. Quant. Spectrosc. Radiat. Transfer **107**, 458–469 (2007). [CrossRef]

**49. **G. N. Plass and G. W. Kattawar, “Monte Carlo calculations of light scattering from clouds,” Appl. Opt. **7**, 415–419 (1968). [CrossRef] [PubMed]

**50. **G. W. Kattawar and G. N. Plass, “Radiance and polarization of multiple scattered light from haze and clouds,” Appl. Opt. **7**, 1519–1527 (1968). [CrossRef] [PubMed]

**51. **G. I. Marchuk, G. A. Mikhailov, M. A. Nazaraliev, R. A. Darbinjan, B. A. Kargin, and B. S. Elepov, *the Monte Carlo Methods in Atmospheric Optics*(Springer-Verlag, Berlin, 1980).

**52. **G. W. Kattawar, editor, *Selected Papers on Multiple Scattering in Plane Parallel Atmospheres and Oceans: Methods* SPIE milestone seres, (MS42, Bellingham, WA, 1991).

**53. **L. Roberti and C. Kummerow, “Monte Carlo calculations of polarized microwave radiation emerging from cloud structures,” J. Geophy. Res. **104**, 2093–2104 (1999). [CrossRef]

**54. **H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,” Appl. Opt. **40**, 400–412 (2001). [CrossRef]

**55. **D. M. O’Brien, “Accelerated quasi Monte Carlo integration of the radiative transfer equation,” J. Quant. Spec-trosc. Radiat. Transfer **48**, 41–59 (1992). [CrossRef]

**56. **
I3RC group, “I3RC Monte Carlo community model of 3D radiative transfer,” http://i3rc.gsfc.nasa.gov/I3RC-intro.html.

**57. **G. W. Kattawar and C. N. Adams, “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 radiance and polarization,” Limnol. Oceanogr. **34**, 1453–1472 (1989). [CrossRef]

**58. **H. R. Gordon, “Ship perturbation of irradiance measurements at sea. 1: Monte Carlo simulations,” Appl. Opt. **24**, 4172–4182 (1985). [CrossRef] [PubMed]

**59. **C. D. Mobley, B. Gentili, H. R. Gordon, Z. Jin, G. W. Kattawar, A. Morel, P. Reinersman, K. Stamnes, and R. H. Stavn, “Comparison of numerical models for computing underwater light fields,” Appl. Opt. **32**, 7484–7504 (1993). [CrossRef] [PubMed]

**60. **Z. Jin and K. Stamnes, “Radiative transfer in nonuniformly refracting layered media: atmosphere-ocean system,” Appl. Opt. **33**, 431–442 (1994). [CrossRef] [PubMed]

**61. **B. Bulgarelli, V. B. Kisselev, and L. Roberti, “Radiative Transfer in the Atmosphere-Ocean System: The Finite-Element Method,” Appl. Opt. **38**, 1530–1542 (1999), http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-38-9-1530. [CrossRef]

**62. **P. N. Reinersman and K. L. Carder, “Monte Carlo simulation of the atmospheric point-spread function with an application to correction for the adjacency effect,” Appl. Opt. **34**, 4453–4471 (1995). [CrossRef] [PubMed]

**63. **Z. Jin, T. P. Charlock, K. Rutledge, K. Stamnes, and Y. Wang, “Analytical solution of radiative transfer in the coupled atmosphereocean system with a rough surface,” Appl. Opt. **45**, 7443–7455 (2006). [CrossRef] [PubMed]

**64. **J. Chowdhary, B. Cairns, and L. D. Travis, “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), http://www.opticsinfobase.org/abstract.cfm?URI=ao-45-22-5542. [CrossRef] [PubMed]

**65. **P. Zhai, G. W. Kattawar, and P. Yang, “Impulse response solution to the three-dimensional vector radiative transfer equation in atmosphere-ocean systems. I. Monte Carlo method,” Appl. Opt. **47**, 1037–1047 (2008), http://www.opticsinfobase.org/abstract.cfm?URI=ao-47-8-1037. [CrossRef] [PubMed]

**66. **A. Morel, “Optical modeling of the upper ocean in relation to its biogenous matter content (Case I Waters),” J. Geophys. Res. **93**, 10479–10768 (1988). [CrossRef]

**67. **
Concise Dictionary of Scientific Biography, (Scribner, New York, 1981), p. 643. Willebrord Snel von Royen used only one l in his last name.

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

**69. **Y. Hu, D. Winker, P. Yang, B. Baum, L. Poole, and L. Vann, “Identification of cloud phase from PICASSO-CENA lidar depolarization: a multiple scattering sensitivity study,” J. Quant. Spectrosc. Radiat. Transfer **70**, 569–579 (2001). [CrossRef]

**70. **A. Hammad and S. Chapman, “The primary and secondary scattering of sun light in a plane-stratified atmosphere of uniform composition,” Philos. Mag. **28**, 99–110 (1939).

**71. **J. Lenoble, editor. *Radiative transfer in scattering and absorbing atmospheres: standard computational procedures* (A. Deepak Publishing, 1985).

**72. **I. Kuščer and M. Ribarič, “Matrix formalism in the theory of diffusion of light,” Opt. Acta **6**, 42–51 (1959). [CrossRef]

**73. **C. E. Siewert, “On the phase matrix basic to the scattering of polarized light,” Astron. Astrophys. **109**, 195–200 (1982).

**74. **J. W. Hovenier and C. V. M. van der Mee, “Fundamental relationships relevant to the transfer of polarized light in a scattering atmosphere,” Astron. Astrophys. **128**, 1–16 (1983).

**75. **I. M. Gel’fand and Z. Y. Sapiro, “Representation of the group of rotations of 3-dimensional space and their applications,” Amer. Math. Soc. Transl. **2**, 207–316 (1956).

**76. **E. P. Wigner, *Group theoy and its application to the quantum mechanics of atomic spectra*(Academic Press, New York, 1959).

**77. **M. I. Mishchenko, L. D. Travis, and A. A. Lacis, *Scattering, Absorption, and Emission of Light by Small Particles* (Cambridge University Press, Cambridge2002).

**78. **J. L. Deuzé, M. Herman, and R. Santer, “Fourier series expansion of the transfer equation in the oceanatmosphere system,” J. Quant. Spectrosc. Radiat. Transfer **41**, 483–494 (1989). [CrossRef]

**79. **C. V. M. Van der Mee and J. W. Hovenier, “Expansion coefcients in polarized light transfer,” Astron. Astrophys. **228**, 559–568 (1990).

**80. **A. Kylling and K. Stamnes, “Efficient yet accurate solution of the linear transport equation in the presence of internal sources: the exponential-linear-in depth approximation,” J. Comput. Phys. **102**, 265–276 (1992). [CrossRef]

**81. **L. C. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. **93**, 70–83 (1941). [CrossRef]

**82. **T. J. Petzold, *Volume Scattering Functions for Selected Ocean Waters* (Scripps Institution of Oceanography, 1977).

**83. **K. J. Voss and E. S. Fry, “Measurement of the Mueller matrix for ocean water,” Appl. Opt. **23**, 4427–4439 (1984). [CrossRef] [PubMed]

**84. **M. Herman, J.L. Deuzé, A. Marchand, B. Roger, and P. Lallart, “Aerosol remote sensing from POLDER/ADEOS over the ocean: improved retrieval using a nonspherical particle model,” J. Geophys. Res. **110**, D10S02 (2005). [CrossRef]

**85. **B. Cairns, E. E. Russell, J. D. LaVeigne, and P. M. W. Tennant. “Research scanning polarimeter and airborne usage for remote sensing of aerosols,” Proc SPIE **5158**, 33–44 (2003). [CrossRef]

**86. **M. I. Mishchenko, B. Cairns, G. Kopp, C. F. Schueler, B. A. Fafaul, J. E. Hansen, R. J. Hooker, T. Itchkawich, H. B. Maring, and L. D. Travis, “Accurate monitoring of terrestrial aerosols and total solar irradiance: introducing the Glory Mission,” Bull. Amer. Meteorol. Soc. **88**, 677–691 (2007). [CrossRef]

**87. **W. J. Wiscombe, “The Delta-M method: rapid yet accurate radiative flux calculations for strongly asymmetric phase functions,” J. Atmos. Sci. **34**, 1408–1422 (1977). [CrossRef]

**88. **Y. -X. Hu, B. Wielicki, B. Lin, G. Gibson, S. -C. Tsay, K. Stamnes, and T. Wong, “δ-Fit: A fast and accurate treatment of particle scattering phase functions with weighted singular-value decomposition least-squares fitting,” J. Quant. Spectrosc. Radiat. Transfer **65**, 681–690 (2000). [CrossRef]

**89. **T. Nakajima and M. Tanaka, “Algorithms for radiative intensity calculations in moderately thick atmospheres using a truncation approximation,” J. Quant. Spectrosc. Radiat. Transfer **40**, 51–69 (1988). [CrossRef]