Comparison of discrete ordinate and Monte Carlo simulations of polarized radiative transfer in two coupled slabs with different refractive indices

D. Cohen; S. Stamnes; T. Tanikawa; E. R. Sommersten; J. J. Stamnes; J. K. Lotsberg; K. Stamnes

doi:10.1364/OE.21.009592

1. Introduction

In remote sensing of the Earth from space one goal is to retrieve atmospheric and surface parameters from measurements of the radiation emerging at the top-of-the-atmosphere (TOA) at a number of wavelengths [1, 2]. These retrieval parameters (RPs), such as aerosol type and loading, and concentrations of aquatic constituents in an open ocean or coastal water area, depend on the inherent optical properties (IOPs) of the atmosphere and the water. If there is a model providing a link between the RPs and the IOPs, a forward radiative transfer (RT) model can be used to compute how the measured TOA radiation field should respond to changes in the RPs, and an inverse RT problem can be formulated and solved to derive information about the RPs [3, 4]. A forward RT model, employing IOPs that describe how atmospheric and aquatic constituents absorb and scatter light, can be used to compute the multiply scattered light field in any particular direction (with specified polar and azimuth angles) at any particular depth level (including the TOA) in a vertically stratified medium, such as a coupled atmosphere-water system [5,6]. In order to solve the inverse RT problem it is important to have an accurate and efficient forward RT model. Accuracy is important in order to obtain reliable and robust retrievals, and efficiency is an issue because standard iterative solutions of the nonlinear inverse RT problem require running the forward RT model repeatedly to compute the radiation field and its partial derivatives with respect to the RPs (the Jacobians) [3,4].

Recently RT theory has been employed in the detection of skin diseases. Letting light at different wavelengths illuminate the skin and analyzing the reflected light using an RT model for a coupled air-tissue system, one may detect differences between the IOPs of normal and malignant skin tissues [7–9].

While solutions to the scalar RT equation (RTE), which involves only the first component of the Stokes vector (the radiance or intensity), are well developed, modern RT models that solve the vector RTE are capable of calculating also polarization effects described by the second, third, and fourth component of the Stokes vector. Even if one’s interest would lie primarily in the radiance, it is important to realize that solutions of the scalar RTE, which ignores polarization effects, introduce errors in the computed radiances [10–12].

As mentioned above, in scalar RT theory only the first component of the Stokes vector (the radiance) is considered. Numerical realizations using the discrete ordinate method developed by Chandrasekhar and others in the 1940’s [13] for the scalar problem, were first explored in the early 1970’s [14]. The scalar RT method was further developed in the 1980’s [15–17], and eventually implemented numerically in a robust, well-tested and freely available computer code (DISORT) in 1988 [18]. The DISORT code computes radiances and irradiances in a vertically inhomogenous medium that is divided into homogenous layers to resolve the vertical variation in the IOPs. The method was extended to be applicable to two adjacent media with different refractive indices, such as a coupled atmosphere-water system [5] or a coupled air-tissue system [19], and the performance of the resulting C-DISORT codes was tested against Monte Carlo (MC) computations [20,21].

The DISORT and MC RT models have different advantages and drawbacks, and should be treated as complementary. The DISORT RT model applies only to stratified slabs, but is efficient and requires much less computer time than the MC RT model, which on the other hand can easily be extended to solve RT problems including a 3D geometry. Since these two RT models are based on different numerical techniques, a comparison between them can be used to test their performance. The DISORT RT model relies on the replacement of the integral in the RTE by a discrete sum and employs methods of linear algebra to solve the RTE, whereas the MC RT model relies on the use of probability concepts and statistics to solve the RT problem by simulating the random walk of photons in a host medium.

For a layered medium with constant refractive index, solutions of the vector RTE [see Eq. (1)] were developed and implemented in the VDISORT code [22–25]. For a coupled system consisting of two media with different refractive indices, the Coupled Vector Discrete Ordinate (C-VDISORT) RTM [26, 27] was developed and tested for pure Rayleigh scattering [28]. It may seem trivial to go from a Rayleigh scattering medium to a Mie scattering medium. However, we found that this extension required a lot of attention to detail, and also further development of both the coupled VDISORT (C-VDISORT) and Monte Carlo (C-PMC) codes in order to get them to produce correct results. Also, we are not aware of detailed comparisons between two such radically different RT formulations and numerical results produced by the corresponding computational codes for the coupled case. The comparisons with benchmark results are also a new feature and this effort proved to be very helpful in checking the results. These benchmark results were not available when the results for Rayleigh scattering [28] were produced. For that special case, it was easy to distinguish between the direct and diffuse components of the Stokes vector in the C-PMC code. The direct component was simply separated in post-processing [28]. For a medium consisting of strongly forward-scattering particles we needed to add a counter to keep track of scattering events in order to separate the diffuse component from the direct one. The quadrature of the phase matrix also had to be increased and the number of azimuthal and polar detector bins had to be significantly increased to get accurate results from the C-PMC code. For the C-VDISORT code we had to create a new subroutine to handle a generalized phase matrix, and the δ-M technique [29] was needed to get accurate results without using a large number of streams which adversely affects the computational efficiency.

In this paper, we extend the C-VDISORT code to be applicable also to non-Rayleigh scattering in a coupled, stratified system and test it against Coupled Polarized Monte Carlo (C-PMC) RT model simulations [30–37]. We compare both radiances as well as other components of the Stokes vector for non-Rayleigh scattering in a coupled, stratified system by solving the vector RTE for the Stokes vector using both the C-VDISORT and the C-PCM RT models. We start by comparing the results produced by the two codes with recently published benchmark results for a single, homogeneous slab with a constant refractive index that has strong forward scattering due to aerosol and cloud particles. We then compare output from these two RT models at 3 different levels in a stratified coupled system: (i) just below the interface between the two slabs (with different refractive indices), (ii) just above the interface, and (iii) at the top of the upper slab. For an atmosphere-water system these levels correspond to locations just below the water surface, just above the water surface, and at the top-of-the-atmosphere (TOA). We focus on particle scattering in the Mie regime in both media because the previous study comparing C-VDISORT and C-PMC results [28] considered pure Rayleigh scattering only.

2. Radiative transfer theory

2.1. General radiative transfer

In a plane-parallel (slab) geometry, the integro-differential equation for polarized radiation, given in terms of the Stokes vector I⃗ = [I_||, I_⊥, U, V]^T, where the superscript T denotes the transpose, is [13]

\begin{array}{l} u \frac{d \vec{I} (τ, u, ϕ)}{d τ} & = & \vec{I} (τ, u, ϕ) + \vec{Q} (τ, u, ϕ) \\ - \frac{a (τ)}{4 π} \int_{0}^{2 π} d ϕ^{'} \int_{- 1}^{1} d u^{'} M (τ, u, ϕ; u^{'}, ϕ^{'}) \vec{I} (τ, u^{'}, ϕ^{'}) . \end{array}

Here I_|| and I_⊥ are the intensity components that are parallel and perpendicular, respectively, to the scattering plane, which is spanned by the directions of the incident and the scattered radiation, U is the degree of linear polarization in the 45°/135° plane, and V is the degree of circular polarization. This Stokes vector is related to the more common one, I⃗_S = [I, Q, U, V]^T, by I = I_|| + I_⊥ and Q = I_|| − I_⊥.

In Eq. (1), dτ = (α + σ)dz is the differential optical depth, α and σ are the absorption and scattering coefficients, respectively, dz is the differential vertical path length, u′ = cosθ′, θ′ and ϕ′ are the polar angle and the azimuth angle, respectively, of the direction of an incident parallel beam, u = cosθ, θ and ϕ are the polar angle and the azimuth angle, respectively, of the observation direction, a(τ) = σ/(α + σ) is the single-scattering albedo, M is the phase matrix describing the scattering properties of the medium [38–40], and Q⃗ is the source term, which has different expressions in the two media of different refractive indices comprising a coupled system [5, 28]. For the atmosphere-water system, the source terms Q⃗_atm and Q⃗_w are given by Eqs. (10) and (11) in [28]. Explicit expressions for the phase matrices used as input to the C-VDISORT and C-PMC RT models are given in Appendix A.

At the interface between the two media (such as the atmosphere-water or air-tissue interface), the reflection law, Snell’s law, and Fresnel’s equations are used to account for the reflection, refraction, and transmission of radiation due to the change in the refractive index. In this study, we consider a smooth planar interface between the two media.

2.2. Types of scattering & input/output

The type of scattering event occurring depends on the size of the particle compared to the wavelength of the radiation. For a spherical particle, the size parameter is defined as

x = \frac{2 π r}{λ}

where r is the radius of the particle and λ is the wavelength of the incident light [41,42]. When x << 1, the scattering is said to be in the Rayleigh regime, which corresponds to scattering by an electric dipole. When 0.1 < x < 50, the scattering is called intermediate or non-Rayleigh scattering. For spherical particles, this regime is referred to as Mie scattering, because solutions to the scattering problem applicable to dielectric spheres published by Gustav Mie in 1908 can be applied [42]. For larger particles, the Mie series solution converges very slowly, and raytracing methods based on geometrical optics are better suited to solve the scattering problem [42,43].

For the simulations considered in this paper we used a collection of spheres with sizes described by a normalized log-normal density distribution, n(r), the number of particles with radius r in a volume $\frac{4 π r^{3}}{3}$ so that the unit is m⁻¹, given by [44]:

n (r) = \frac{1}{\sqrt{2 π} σ_{g}^{2} r} exp [- {(ln r - ln r_{g})}^{2} / (2 σ_{g}^{2})], \int_{0}^{\infty} n (r) d r = 1

where

ln r_{g} = \int_{0}^{\infty} ln r n (r) dr

and σ_g, the dimensionless standard deviation of the size distribution, is given by

σ_{g}^{2} = \int_{0}^{\infty} {(ln r - ln r_{g})}^{2} n (r) d r .

The effective radius r_eff and the effective variance ν_eff are defined as (

G = \int_{r_{1}}^{r_{2}} π r^{2} n (r) d r

)

r_{eff} = \frac{1}{G} \int_{r_{1}}^{r_{2}} r π r^{2} n (r) d r = \frac{r_{g}}{{(1 + ν_{eff})}^{5 / 2}}

ν_{eff} = \frac{1}{G r_{eff}^{2}} \int_{r_{1}}^{r_{2}} {(r - r_{eff})}^{2} π r^{2} n (r) d r = exp [σ_{g}^{2}] - 1 .

Here r₁ and r₂ are respectively the radii of the smallest and largest particles in the distribution, and Eqs. (3)–(5) have been used in the last steps of Eqs. (6) and (7) [44].

Aerosol particles

In Section 4 we compare results produced by C-VDISORT and C-PMC with benchmark computations for a homogeneous slab of aerosol particles provided by Kokhanovsky et al. [45]. The aerosol particles considered in the benchmark computations had r_g = 0.3 μm, and σ_g = 0.92, and the smallest and largest particle radii were selected to be r₁ = 0.005 μm and r₂ = 30 μm. The refractive index of the aerosol particles was set to m = 1.385, which yields a single-scattering albedo of 1.0, and an asymmetry factor g = 0.79275 [45] [see Eq. (8) below].

Cloud particles

Kokhanovsky et al. [45] also provided benchmark results for a homogeneous slab of cloud particles, for which they used r_g = 5 μm, and σ_g = 0.4, and the smallest and largest particle radii were selected to be r₁ = 0.005 μm and r₂ = 100 μm. The refractive index was set to m = 1.339, which yields a single-scattering albedo of 1.0, and an asymmetry factor g = 0.86114 for the cloud particles.

The phase matrix determines how the Stokes vector changes as the result of a scattering event. The element in the first row and first column of the phase matrix is the scattering phase function p(cosΘ) ≡ a₁(Θ), which is the only element that matters if polarization effects are ignored. The first moment of the scattering phase function:

g = \frac{1}{4 π} \int_{4 π} d ω p (cos Θ) cos Θ = \frac{1}{2} \int_{- 1}^{1} d (cos Θ) p (cos Θ) cos Θ

is called the asymmetry factor of the scattering phase function. It is positive for scattering in the forward direction, negative for scattering in the backward direction, and zero for scattering phase functions that are symmetric about Θ = π/2.

To determine the IOPs for non-Rayleigh scattering, we may use a Mie-Lorentz solution for a collection of spheres or a T-matrix code for scattering by a collection of non-spherical particles [46–48]. For a given size distribution of the particles and refractive index relative to the medium in which they are embedded, the IOPs are the absorption coefficient, the scattering coefficient, and the scattering phase matrix. For the benchmark comparisons we used IOPs provided by Kokhanovsky et al. [45].

2.3. The interface between the two media

At the interface between the two media (such as the air-water interface), the difference in the refractive index must be taken into account. The reflection and refraction at the interface are described by reflection and transmission matrices, both for illumination from the upper medium and the lower medium. In the latter case, if the upper slab consists of air and the lower slab of water, the critical angle must be considered to determine whether or not total internal reflection takes place at the water-air interface.

The reflection law θ_r = θ_i, where θ_i and θ_r are the angles of incidence and reflection, respectively, applies at the interface. The angle of refraction (or transmission) θ_t is determined by Snell’s law n₁ sinθ_i = n₂ sinθ_t, where n₁ and n₂ are the real parts of the refractive indices in medium 1 (containing the incident beam) and medium 2 (containing the transmitted beam), respectively. The reflection and transmission coefficients for parallel and perpendicular polarization follow from Fresnel’s equations [42, 49], and are given by

R_{| |} = \frac{E_{| | r}}{E_{| | i}} = \frac{m_{1} cos θ_{t} - m_{2} cos θ_{i}}{m_{1} cos θ_{t} + m_{2} cos θ_{i}}

T_{| |} = \frac{E_{| | t}}{E_{| | i}} = \frac{2 m_{1} cos θ_{i}}{m_{1} cos θ_{t} + m_{2} cos θ_{i}}

for parallel polarization, and by

R_{⊥} = \frac{E_{⊥ r}}{E_{⊥ i}} = \frac{m_{1} cos θ_{i} - m_{2} cos θ_{t}}{m_{1} cos θ_{i} + m_{2} cos θ_{t}}

T_{⊥} = \frac{E_{⊥ t}}{E_{⊥ i}} = \frac{2 m_{1} cos θ_{i}}{m_{1} cos θ_{i} + m_{2} cos θ_{t}}

for perpendicular polarization, where m₁ = n₁ + in′₁ and m₂ = n₂ +in′₂ are the complex refractive indices, E_|| is the component of electric field that is parallel to the plane of incidence (the plane spanned by the unit vector in the direction of the incident beam and the interface normal), while E_⊥ is the component perpendicular to the plane of incidence, R_|| and T_|| are the reflection and transmission coefficients for the parallel polarization component, while R_⊥ and T_⊥ are the reflection and transmission coefficients for the perpendicular polarization component.

For reflection of a collimated beam incident on the interface, one multiplies the incident Stokes vector I⃗_inc = [I_||i, I_⊥i, U_i, V_i]^T by the reflection matrix to obtain the Stokes vector of the reflected beam [42]

(\begin{matrix} I_{| | r} \\ I_{⊥ r} \\ U_{r} \\ V_{r} \end{matrix}) = (\begin{matrix} S_{11}^{r} & 0 & 0 & 0 \\ 0 & S_{22}^{r} & 0 & 0 \\ 0 & 0 & S_{33}^{r} & S_{34}^{r} \\ 0 & 0 & S_{43}^{r} & S_{44}^{r} \end{matrix}) (\begin{matrix} I_{| | i} \\ I_{⊥ i} \\ U_{i} \\ V_{i} \end{matrix})

where the elements of the reflection matrix are given by [42]

S_{11}^{r} = R_{| |} R_{| |}^{*} = {| R_{| |} |}^{2}

S_{22}^{r} = R_{⊥} R_{⊥}^{*} = {| R_{⊥} |}^{2}

S_{33}^{r} = S_{44}^{r} = Re (R_{| |} R_{⊥}^{*})

S_{34}^{r} = Im (R_{| |} R_{⊥}^{*})

S_{43}^{r} = - S_{34}^{r} = Im (R_{| |}^{*} R_{⊥}) .

In the case of total internal reflection, the reflection coefficients in Eqs. (9) and (11) are complex numbers even for the case in which m₁ = n₁ and m₂ = n₂ are real numbers (n′₁ = n′₂ = 0). Then the absolute value of each of the reflection coefficients is 1, as expected for total internal reflection. Total internal reflection couples the two Stokes vector elements U and V since the off-diagonal elements

S_{34}^{r}

and

S_{43}^{r}

in Eq. (13) are non-vanishing.

The Stokes vector of the transmitted collimated beam is related to the Stokes vector of the incident collimated beam by [42]

(\begin{matrix} I_{| | t} \\ I_{⊥ t} \\ U_{t} \\ V_{t} \end{matrix}) = (\begin{matrix} S_{11}^{t} & 0 & 0 & 0 \\ 0 & S_{22}^{t} & 0 & 0 \\ 0 & 0 & S_{33}^{t} & S_{34}^{t} \\ 0 & 0 & S_{43}^{t} & S_{44}^{t} \end{matrix}) (\begin{matrix} I_{| | i} \\ I_{⊥ i} \\ U_{i} \\ V_{i} \end{matrix})

where the elements of the transmission matrix are given by [42]

S_{11}^{t} = {| K_{t} |}^{2} T_{| |} T_{| |}^{*} = {| K_{t} |}^{2} {| T_{| |} |}^{2}

S_{22}^{t} = {| K_{t} |}^{2} T_{⊥} T_{⊥}^{*} = {| K_{t} |}^{2} {| T_{⊥} |}^{2}

S_{33}^{t} = S_{44}^{t} = {| K_{t} |}^{2} Re (T_{| |} T_{⊥}^{*})

S_{34}^{t} = {| K_{t} |}^{2} Im (T_{| |} T_{⊥}^{*})

S_{43}^{t} = - S_{34}^{t} = {| K_{t} |}^{2} Im (T_{| |}^{*} T_{⊥}) .

The factor |K_t|² accounts for the different propagation directions of the incident and transmitted beams relative to the interface normal due to the refractive index difference between the two media. Through considerations involving energy conservation and ray bending at the interface, one finds [28]

{| K_{t} |}^{2} = m_{rel}^{3} \frac{μ_{t}}{μ_{i}}

where m_rel = n₂/n₁.

3. Radiative transfer models

Two different methods were used to obtain solutions of the RTE in Eq. (1), both being capable of handling coupling between two slabs of different refractive indices, absorption, and multiple scattering. In order to deal with vertical inhomogeneity each of the two slabs constituting the coupled system was divided into a number of homogenous horizontal layers. The IOPs, including the phase matrix, for each layer as well as the direction of the incident beam, given by the solar angles (θ₀, ϕ₀) for a coupled system, were provided as input to the polarized RT model, which calculated the Stokes vector at user-specified optical depths and observation directions, each given by its polar angle θ and azimuth angle ϕ.

3.1. C-VDISORT

The C-VDISORT RT model is based on solving the RTE in Eq. (1) by use of the discrete ordinate method, implying that the integral in Eq. (1) is replaced by a discrete sum, which leads to a system of coupled ordinary differential equations that is solved using techniques of linear algebra. Boundary conditions are applied at the top of the first medium (TOA for a coupled atmosphere-water system) and at the bottom of the second medium (the water). At the interface between the two media where the refractive index changes abruptly, Snell’s law and Fresnel’s equations are employed to determine the reflection, refraction, and change in polarization of the radiation field. Furthermore, it is required that the radiation field be continuous across each of the interfaces between horizontal layers in each of the two media. The integral term in Eq. (1) contains the phase matrix M, which can be calculated as explained in Appendix A.

3.2. C-PMC

C-PMC RT simulations are based on statistics and the use of Monte Carlo simulations to determine a random walk of a large number of photons that are incident at the top of the first medium (TOA for a coupled atmosphere-water system) and that proceed until each of the photons is either absorbed or scattered back to space. For each photon incident at the TOA, a random number ρ between 0 and 1 is generated to determine its initial optical path length τ. To compute the new path length τ one generates a new random number ρ between 0 and 1, and τ is computed from the formula

τ = \sum_{j = 1}^{N} d_{j} k_{e, j} = ln ρ

where d_j is the geometric path length and k_e,j is the attenuation coefficient. The summation is carried out until the sum becomes larger than ln ρ, and then d_N is adjusted in layer N so that the sum in Eq. (25) becomes equal to ln ρ. This procedure determines the new vertical position of the beam. Once the new vertical position of the beam has been calculated, the attenuation coefficient, given by k_e = α + σ, where α and σ are the absorption and scattering coefficients, respectively, determines whether or not the beam is scattered or absorbed. If

ρ > \frac{σ}{k_{e}}

, the beam is absorbed; otherwise it is scattered. Next, a new random number ρ between 0 and 1 is generated and used to determine whether the scattering is of the Rayleigh or non-Rayleigh type, each defined by a user-provided Mueller matrix.

At the interface between the two media, where the refractive index changes abruptly, a random number ρ is used to determine whether the beam is reflected or refracted by comparing it to the reflectance R for unpolarized light given by

R = \frac{1}{2} [\frac{{sin}^{2} (θ_{i} - θ_{t})}{{sin}^{2} (θ_{i} + θ_{t})} + \frac{{tan}^{2} (θ_{i} - θ_{t})}{{tan}^{2} (θ_{i} + θ_{t})}] .

If ρ < R the beam is reflected; otherwise the beam is refracted according to Snell’s law, for details see [30, 31]. In Eq. (26),

R = (R_{| |}^{2} + R_{⊥}^{2}) / 2

where R_|| and R_⊥ are the Fresnel reflection coefficients given by Eqs. (9) and (11).

In C-PMC simulations, each beam contains information about the full Stokes vector, but two rotations are required to relate the beam’s initial and final Stokes vectors, as discussed in Appendix A.

4. Comparison of C-VDISORT and C-PMC RT simulations

In the simulations presented in this paper, identical inputs to the C-VDISORT and C-PMC RT codes were used. In Figs. 1–5 we focus on comparisons with benchmark results for reflected and transmitted polarized radiation for a single slab with constant refractive index. Next, in Figs. 7–9 we compare C-VDISORT and C-PMC results for Stokes vector components I, Q, and U, and the degree of linear polarization $P = \sqrt{Q^{2} + U^{2}} / I$ at 3 specific levels of observation in a two-slab system where the index of refraction is set to 1.0 in the upper slab, and to 1.338 in the lower slab. The 3 levels of observation are: (i) just below the interface between the two slabs, (ii) just above the interface, and (iii) at the top of the upper slab. These results are shown as a function of the cosine of the polar angle of observation. All benchmark results (Figs. 1–4) are for a collimated (solar) beam incident at a solar zenith angle of 60°.

Fig. 1 Left four panels: Comparison of reflected polarized radiation components computed by C-VDISORT (dashed curves) with benchmark results [45] (solid curves) for a homogeneous layer of non-absorbing aerosol particles, simulated by putting half of the aerosol particles (optical depth = 0.1631) in each slab, and setting the index of refraction to 1.0 in both slabs. From top to bottom: Stokes parameter I, Q, U, and Degree of linear polarization. Right four panels: Same as left four panels except that the comparison is between C-VDISORT (dashed curves) and C-PMC (solid curves). Red: Δϕ = 0°; green: Δϕ = 90°; blue: Δϕ = 180°. 200 discrete ordinate streams, 10⁹ photons.

Download Full Size | PDF

Fig. 2 Similar results as in Fig. 1 but for the transmitted polarized radiation components.

Download Full Size | PDF

Fig. 3 Left three panels: Relative difference between reflected polarized radiation components computed by C-VDISORT and benchmark results for the results shown in Fig. 1. Right three panels: Relative difference between transmitted polarized radiation components computed by C-VDISORT and benchmark results for the results shown in Fig. 2. The color code is the same as in Fig. 1. Note that the U parameter computed by the C-VDISORT code is close to zero for both Δϕ = 0° and Δϕ = 180° (in agreement with the benchmark) which explains why the red and blue curves are missing for the U parameter.

Download Full Size | PDF

Fig. 4 Left four panels: Comparisons of reflected polarized radiation components computed by C-VDISORT (dashed curves) with benchmark results [45] (solid curves) for a homogeneous layer of non-absorbing cloud particles, simulated by setting the refractive index equal to 1.0 in both slabs, and putting half of the cloud particles (optical depth = 2.5) in each slab. From top to bottom: Stokes parameter I, Q, U and Degree of linear polarization. Right four panels: Same as left four panels except that the comparisons are between C-VDISORT (dashed curves) and C-PMC (solid curves). Red: Δϕ = 0°; green: Δϕ = 90°; blue: Δϕ = 180°. 200 discrete ordinate streams, 10⁹ photons.

Download Full Size | PDF

Fig. 5 Results similar to those in Fig. 4, but for the transmitted polarized radiation components. 200 discrete ordinate streams, 10⁹ photons.

Download Full Size | PDF

Fig. 6 Left three panels: Relative difference between reflected polarized radiation components computed by C-VDISORT and benchmark results for the cloud case shown in Fig. 4. Right three panels: Relative difference between transmitted polarized radiation components computed by C-VDISORT and benchmark results for the cloud case shown in Fig. 5. The color code is the same as in Fig. 1.

Download Full Size | PDF

Fig. 7 Comparisons between polarized radiation components computed with C-VDISORT (dashed curves) and C-PMC (solid curves) for the same physical situation as in Fig. 1 except that the refractive index was set equal 1.0 in the upper slab and equal to 1.338 in the lower slab. Upper four panels: just above interface. Lower four panels: just below the interface. Red: Δϕ = 45°; green: Δϕ = 90°; blue: Δϕ = 135°. 64 discrete ordinate streams in upper slab, 96 streams in lower slab, 10⁹ photons.

Download Full Size | PDF

Fig. 8 Similar to Fig. 7 (aerosol particles with an optical depth of 0.1631 in each slab) except that we put non-absorbing aerosol particles of optical depth 0.3262 in the upper slab, and cloud-like particles with single-scattering albedo of 0.9, and optical depth of 5.0 in the lower slab. 96 discrete ordinate streams in upper slab, 128 streams in lower slab, 10⁹ photons.

Download Full Size | PDF

Fig. 9 Upper four panels: Similar to Fig. 7, but at TOA for non-absorbing aerosol particles with an optical depth 0.1631 in each slab. Lower four panels: Similar to Fig. 8, but at TOA for non-absorbing aerosol particles with an optical depth of 0.3262 in the upper slab, and cloud particles with an optical depth of 5.0 and single-scattering albedo of 0.9 in the lower slab.

Download Full Size | PDF

4.1. Mie Scattering in the Rayleigh limit for a coupled system

Results for Rayleigh scattering both in the upper slab (atmosphere) and the lower slab (water) were provided by Sommersten et al. [28], who showed that C-VDISORT and C-PMC results were in good agreement. Since our RT codes encompass scattering by Rayleigh particles as well as non-Rayleigh particles, they should give results that are essentially the same as those for pure Rayleigh scattering provided the size of the particles is much smaller than the wavelength of the incident light [i.e. $\frac{2 π r}{λ} ≪ 1$ , see Eq. (2)]. To test this requirement, we used a Mie-Lorentz code [46–48] to prepare phase matrix input for the C-VDISORT and the C-PMC RT codes for a collection of spherical particles in the Rayleigh regime with radii between r₁ = 14 nm and r₂ = 15 nm [see Eqs. (3)–(7)]. We found that the C-VDISORT and C-PMC RT models produced results in very good agreement with one another, and also with the results obtained for pure Rayleigh scattering [28], thus verifying that both RT models produce correct results in the Rayleigh limit.

4.2. C-VDISORT and C-PMC vs benchmark – aerosol layer – reflection

In the left four panels of Fig. 1, we compare C-VDISORT results (dotted curves) with benchmark results (solid curves) for a layer of aerosol particles [45] with mean radius r_g = 0.3 μm, standard deviation σ_g = 0.92, and refractive index m = 1.385, which yields a single-scattering albedo of 1.0, and an asymmetry factor g = 0.79275 as explained in §2.2. In the simulations using the C-VDISORT and C-PMC RT models, we set the index of refraction to 1.0 in both slabs, and put half of the aerosol particles (optical depth = 0.1631) in each slab. The C-VDISORT results are seen to be indistinguishable from the benchmark results; in fact the dotted and solid curves lie on top of each other in the left panels of Fig. 1. The right four panels of Fig. 1 show comparisons between reflected polarized radiation components computed by the C-VDISORT model (same results as in the four left panels) and the C-PMC RT model for the benchmark case [45]. Except for statistical fluctuations in the C-PMC results, the agreement is seen to be very good. A quantitative assessment of these results is provided in §4.3.

4.3. C-VDISORT and C-PMC vs benchmark – aerosol layer – transmission

The left four panels of Fig. 2 show comparisons between C-VDISORT and benchmark results for the transmitted radiation for the same aerosol benchmark case [45] as in Fig. 1. Except for artifacts in the I, Q, and U parameters around the forward direction (cosθ = −0.5), the agreement between C-VDISORT and the benchmark results [45] is excellent. These artifacts are most likely due to the treatment of the forward scattering peak in the C-VDISORT RT model, which uses the δ-M approximation [29] to truncate the forward scattering peak of the scattering phase function. A comparison between C-VDISORT and C-PMC results for transmitted polarized radiation components are shown in the right four panels of Fig. 2. The C-VDISORT results are the same as those shown in the left four panels. The C-VDISORT and C-PMC RT models give very similar results except for small differences due to statistical fluctuations in the C-PMC results and artifacts around the forward scattering direction due to the use of the δ-M approximation in C-VDISORT.

To quantify the difference between the C-VDISORT and the benchmark results shown in Figs. 1 and 2 we define the relative difference as

DIFF = {[S_{C - VDISORT} - S_{BENCH}] / | max {S_{BENCH}} |} \times 100

where S stands for the I, Q, or U component of the Stokes vector. Note that we are using the absolute value of the maximum magnitude of the Stokes parameter [max{S_BENCH}] in the interval μ ≡ cosθ ∈ [−1, 1] as the point of reference to avoid dividing by values close to zero for the Q and U components. This difference between the C-VDISORT and the benchmark results [Eq. (27)] is provided in Fig. 3 for the aerosol case. The relative difference for the I parameter shown in the upper left panel of Fig. 3 is a fraction of one percent for all viewing angles. The difference in the Q parameter is of similar magnitude except around the forward direction where it lies between 1 and 2%. The lower left panel shows that the difference in the U parameter is a fraction of one percent.

Similar results for the transmission are displayed in the right panels of Fig. 3, which show that the difference is small except around the forward direction where it can be as large as 20–40%.

4.4. C-VDISORT and C-PMC vs benchmark – cloud layer – reflection

For cloud particles, Kokhanovsky et al. [45] used a mean radius of r_g = 5 μm, and a standard deviation of σ_g = 0.4, and the smallest and largest particle radii were selected to be r₁ = 0.005 μm and r₂ = 100 μm. The refractive index was set to m = 1.339, which yields a single-scattering albedo of 1.0, and an asymmetry factor g = 0.86114 for the cloud particles. In the simulations with C-VDISORT and C-PMC, we set the refractive index equal to 1.0 in both slabs, and put half of the cloud particles (optical depth = 2.5) in each slab. The left four panels of Fig. 4 show comparisons of C-VDISORT and benchmark results [45] for the reflected polarized radiation components for this two-layer system of cloud particles with a total optical depth of 5.0. The agreement is seen to be very good. The right four panels of Fig. 4 show comparisons between reflected polarized radiation components computed by the C-VDISORT and C-PMC RT models for the cloud benchmark case. The C-VDISORT and C-PMC RT models give almost indistinguishable results except for small differences due to statistical fluctuations in the C-PMC results and artifacts in the I and Q parameters for Δϕ = 180° around the forward scattering direction due to the use of the δ-M approximation. A quantitative assessment of these results is provided in §4.5.

4.5. C-VDISORT and C-PMC vs benchmark – cloud layer – transmission

The left four panels of Fig. 5 show comparisons between C-VDISORT and benchmark results for the transmitted radiation for the same cloud benchmark case as in Fig. 4. The agreement between C-VDISORT and the benchmark results [45] is excellent. Comparisons between C-VDISORT and C-PMC results for the transmitted polarized radiation components are shown in the right four panels of Fig. 5. Disregarding small differences due to statistical fluctuations in the C-PMC results, we see that the agreement between the C-VDISORT and the C-PMC results is generally very good, except in the forward direction where the C-PMC result is smaller than the C-VDISORT result (which agrees with the benchmark), presumably due to a sampling problem in the Monte Carlo algorithm. For example, if the resolution of detector bins in polar or azimuth-space is insufficiently sampled then a beam scattered in the forward direction will not be detected at exactly the forward direction, but rather at angles close to the forward direction. Or if the quadrature of the phase matrix is too coarse, then a forward scattered beam is actually not scattered precisely in the forward direction, but at an angle close to the forward direction. These problems could presumably be solved by increasing the resolution of the detector bins about the forward direction, and by increasing the resolution of the phase matrix quadrature about the 0° scattering direction. A change in the resolution of the detector bins will also need to be accounted for when averaging over the bins to reduce noise. The C-VDISORT results are the same as those shown in the left four panels.

To quantify the difference between the C-VDISORT and the benchmark results shown in Figs. 4 and 5 for the cloud case we proceed as in the aerosol case discussed above and use Eq. (27) for this purpose. The resulting differences are shown in Fig. 6. For the reflected components (left panels) these differences are similar to those for the aerosol case (Fig. 3) except that the differences around the forward direction are more pronounced. On the other hand, the right panels show that the differences for the transmitted components are less pronounced in the forward direction for the cloud case than for the aerosol case (right panels of Fig. 3).

4.6. C-VDISORT vs C-PMC – aerosol particles – coupled case

Next, we consider results for a case similar to that in Figs. 1 and 2 (right panels) except that the refractive indices are different in the two slabs. To simulate the atmosphere-water system we set the refractive index equal to 1.0 in the upper slab and equal to 1.338 in the lower slab, and as before we put half of the particles (optical depth = 0.1631) in each slab. Comparisons between polarized radiation components computed by the C-VDISORT and C-PMC RT models are shown in Fig. 7 just above the interface (upper four panels), and just below the interface (lower four panels). The two codes give very similar results except for small discrepancies in the Q and U components for directions close to the horizon at Δϕ = 45°. The corresponding results at the TOA are shown in Fig. 9 (upper panel); there is generally very good agreement also at the TOA.

4.7. C-VDISORT vs C-PMC – aerosol/cloud particles – coupled case

Next we consider results for a case similar to that in Figs. 4 and 5 except that we set the refractive index equal to 1.0 in the upper slab and equal to 1.338 in the lower slab, and we put aerosol particles in the upper slab (atmosphere, optical depth = 0.3262) and cloud particles in the lower slab of optical depth = 5.0 to mimic large particles in the ocean. This choice may seem artificial, but our purpose is simply to compare the performance of the two codes for indentical input. Comparisons between polarized radiation components computed by the C-VDISORT and C-PMC RT models are shown in Fig. 8 just above the interface, just below the interface, and the corresponding results at the TOA are shown in Fig. 9 (lower panel). The two codes give very similar results. A plot of the differences between the C-VDISORT and C-PMC results (not shown) indicate that the patterns and magnitudes of the discrepancies are similar to but somewhat larger than those shown in Figs. 3 and 6. The increased difference is likely due to statistical noise and could be dealt with by smoothing using simple running mean averaging over the detector bins, as well as by using more photons. In Monte Carlo simulations there is a trade-off between the number of detector bins one adopts (an increase in the number of bins increases both accuracy and noise) and the number of photons one uses (an increase in the number of photons tends to reduce noise but increase computational time.)

5. Conclusions

We have compared results from two vector RT models for a coupled system consisting of two adjacent slabs with different refractive indices. The two RT models are based on the coupled vector discrete ordinate (C-VDISORT) method and the coupled polarized Monte Carlo (C-PMC) method. Both RT models were extended to apply to particle scattering beyond the Rayleigh limit, and emphasis was placed on testing the RT models for non-Rayleigh scattering by particles with sizes comparable to or larger than the wavelength of light (Mie scattering). If the two adjacent slabs are vertically stratified, they will apply to a coupled atmosphere-water system. The comparisons were carried out for small particles in the Rayleigh scattering limit as well as for large particles in the Mie scattering regime, and for the same or different refractive indices in the two slabs with identical input IOPs to the two RT models. For a number of different cases we found the two RTMs to give very similar results for identical input IOPs. Since comparisons of results from C-PMC and C-VDISORT RT models previously were provided only for pure Rayleigh scattering [28], our main focus was to extend the treatment in [28] to non-Rayleigh scattering media.

The results presented here show that our extension of the two RT models to include particle scattering in the Mie regime gives results that agree with benchmark results in the limit of no coupling or change in the refractive index between the two media. Also, good agreement was obtained between the two RT models for cases of coupling when the refractive indices in the two slabs were different, indicating that the scattering phase matrix as well as the interface conditions have been correctly implemented in the two RT models.

The C-VDISORT RT model solves the integro-differential radiative transfer equation for polarized radiation using a combination of analytic and numerical techniques based on linear algebra, while the C-PMC RT model is based on probabilistic and statistical concepts. Since these two radically different RT models give essentially the same results, we can confidently apply either RT model in various applications. The advantage of C-VDISORT RT model is that it is very fast compared to C-PMC RT model. The disadvantage of the current C-VDISORT RT model is that it is limited to plane parallel geometry, while the C-PMC RT model can readily be extended to apply to any desired geometry including 3-D RT simulations.

6. Appendix A

6.1. Stokes vector representation I⃗_S = [I, Q, U, V]^T used in C-PMC

The scattering plane is spanned by the directions of propagation of the incident parallel beam and the scattered parallel beam. We define two unit vectors r̂_⊥ and r̂_||, which are perpendicular and parallel, respectively, to the scattering plane, such that r̂_|| × r̂_⊥ is a unit vector in the direction of propagation of the incident beam. In this (r̂_||, r̂_⊥) reference plane, the Stokes vector ${\vec{I}}_{S}^{sca}$ of the radiation scattered by a particle is related to that of the Stokes vector ${\vec{I}}_{S}^{inc}$ of the incident radiation by a 4×4 Mueller matrix. If we rotate the reference plane about the direction of propagation through an angle ω (0 ≤ ω ≤ 2π) in the clockwise direction when looking in the direction of propagation of the incident beam, then a Stokes vector I⃗′_S in the rotated reference plane is related to the Stokes vector I⃗_S in the (r̂_||, r̂_⊥) plane by

{\vec{I}}_{S}^{'} = R_{S} (ω) {\vec{I}}_{S}

where

R_{S} (ω) = [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & cos (2 ω) & - sin (2 ω) & 0 \\ 0 & sin (2 ω) & cos (2 ω) & 0 \\ 0 & 0 & 0 & 1 \end{matrix}]

is called the Stokes rotation matrix.

For scattering by a small volume of particles, the ensemble-averaged Mueller matrix is referred to as the Stokes scattering matrix F_S. If any of the following conditions are fulfilled [51]:

each particle in the volume element has a plane of symmetry, and the particles are randomly oriented,
each volume element contains an equal number of particles and their mirror particles in random orientation,
the particles are much smaller than the wavelength of the incident light,

then the Stokes scattering matrix pertaining to the representation I⃗_S = [I, Q, U, V]^T has the following form

F_{S} (Θ) = [\begin{matrix} a_{1} & b_{1} & 0 & 0 \\ b_{1} & a_{2} & 0 & 0 \\ 0 & 0 & a_{3} & b_{2} \\ 0 & 0 & - b_{2} & a_{4} \end{matrix}] .

Each of the six independent matrix elements in Eq. (30) depends on the scattering angle Θ, but are independent of the optical depth within a given horizontal layer of the medium. For spherical particles, the matrix in Eq. (30) simplifies, since a₁ = a₂ and a₃ = a₄, so that only four independent elements remain.

The phase matrix used in C-PMC is obtained from the Stokes scattering matrix F_S(Θ) in Eq. (30) by

M_{S} (Θ) = R_{S} (Ψ) F_{S} (Θ) R_{S} (Φ)

where R_S is the rotation matrix described in Eq. (29) with ω replaced by Ψ or Φ where Ψ and Φ are angles describing the rotation of the incident and scattered parallel beams from the frame of the scattering plane to the laboratory frame [13].

In C-PMC simulations, two rotations are required to obtain the Stokes vector ${\vec{I}}_{S}^{sca}$ of the scattered parallel beam of radiation from the Stokes vector ${\vec{I}}_{S}^{inc}$ of the incident parallel beam of radiation:

{\vec{I}}_{S}^{sca} = M_{S} (Θ) {\vec{I}}_{S}^{inc} = R_{S} (Ψ) F_{S} (Θ) R_{S} (Φ) {\vec{I}}_{S}^{inc} .

According to Eq. (32), the Stokes vector

{\vec{I}}_{S}^{inc}

of the incident parallel beam must be multiplied by the rotation matrix R_S(Φ) before it is multiplied by the Stokes scattering matrix F_S(Θ), whereafter it must be multiplied by the rotation matrix R_S(Ψ). These matrix multiplications are carried out explicitly in C-PMC, while they are automatically taken care of in C-VDISORT.

6.2. Stokes vector representation I⃗ = [I_||, I_⊥, U, V]^T used in C-VDISORT

The Stokes vector representation used in C-VDISORT, I⃗ = [I_||, I_⊥, U, V]^T, is related to I⃗_S = [I, Q, U, V]^T used in C-PMC by

{\vec{I}}_{S} = D \vec{I}

where

D = [\begin{matrix} 1 & 1 & 0 & 0 \\ 1 & - 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}]

so that I = I_|| + I_⊥, and Q = I_|| − I_⊥. Since [see Eq.(28)]

{\vec{I}}^{'} = D^{- 1} {\vec{I}}_{S}^{'} = D^{- 1} R_{S} (ω) {\vec{I}}_{S} = D^{- 1} R_{S} (ω) D \vec{I}

the rotation matrix for the Stokes vector in the representation I⃗ = [I_||, I_⊥, U, V]^T becomes:

R (ω) = D^{- 1} R_{S} (ω) D = [\begin{matrix} {cos}^{2} ω & {sin}^{2} ω & - \frac{1}{2} sin 2 ω & 0 \\ {sin}^{2} ω & {cos}^{2} ω & \frac{1}{2} sin (2 ω) & 0 \\ sin (2 ω) & - sin (2 ω) & cos 2 ω & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] .

The phase matrix M in the Stokes vector representation I⃗ = [I_||, I_⊥, U, V]^T is related to phase matrix M_S in the Stokes vector representation I⃗_S = [I, Q, U, V]^T by

M (Θ) = D^{- 1} M_{S} (Θ) D .

Similarly, the Stokes scattering matrix F(Θ) associated with the Stokes vector representation I⃗ = [I_||, I_⊥, U, V]^T used in Eq. (1) is related to the Stokes scattering matrix F_S(Θ) in Eq. (30) by

F (Θ) = D^{- 1} F_{S} (Θ) D = (\begin{matrix} \frac{1}{2} (a_{1} + a_{2} + 2 b_{1}) & \frac{1}{2} (a_{1} - a_{2}) & 0 & 0 \\ \frac{1}{2} (a_{1} - a_{2}) & \frac{1}{2} (a_{1} + a_{2} - 2 b_{1}) & 0 & 0 \\ 0 & 0 & a_{3} & b_{2} \\ 0 & 0 & - b_{2} & a_{4} \end{matrix}) .

For Rayleigh scattering with depolarization factor ρ the Stokes scattering matrix in the Stokes vector representation I⃗_S = [I, Q, U, V]^T used in C-PMC is given by [13,28]

F_{S} (Θ) = \frac{3}{4} \frac{(1 - ρ)}{(1 + ρ / 2)} [\begin{matrix} \frac{1 + ρ}{1 - ρ} + {cos}^{2} Θ & - {sin}^{2} Θ & 0 & 0 \\ - {sin}^{2} Θ & 1 + {cos}^{2} Θ & 0 & 0 \\ 0 & 0 & 2 cos Θ & 0 \\ 0 & 0 & 0 & \frac{2 - 4 ρ}{1 - ρ} cos Θ \end{matrix}] .

In the Stokes vector representation I⃗ = [I_||, I_⊥, U, V]^T, used in C-VDISORT, the corresponding Stokes scattering matrix for Rayleigh scattering becomes (using Eq. (37) [13]:

F (Θ) = \frac{3}{2 (1 + 2 γ)} (\begin{matrix} {cos}^{2} Θ + γ {sin}^{2} Θ & γ & 0 & 0 \\ γ & 1 & 0 & 0 \\ 0 & 0 & (1 - γ) cos Θ & 0 \\ 0 & 0 & 0 & (1 - 3 γ) cos Θ \end{matrix})

where ρ = 2γ/(1 + γ) or γ = ρ/(2 − ρ).

6.3. Generalized Spherical Functions

In the C-VDISORT RT model, the number of streams is equal to the number of terms in the discrete sum that replaces the integral in the solution of the RT problem. Each element of the phase matrix M in Eq. (1) is expanded in generalized spherical functions with particular coefficients. To that end, the phase matrix M is first expanded in a Fourier series in the azimuth angles using 2N coefficients as shown in the following equation:

M (u, ϕ; u^{'}, ϕ^{'}) = \sum_{m = 0}^{2 N - 1} {M_{m}^{c} (u, u^{'}) cos m (ϕ^{'} - ϕ) + M_{m}^{s} (u, u^{'}) sin m (ϕ^{'} - ϕ)} .

Here

M_{m}^{c}

is the coefficient of the cosine term of the Fourier series, while

M_{m}^{s}

is the coefficient of the sine term of the Fourier series, u = cosθ, θ being the polar angle after scattering, while u′ = cosθ′, θ′ being the polar angle prior to scattering, ϕ is the azimuth angle after scattering, and ϕ′ is the azimuth angle prior to scattering. An addition theorem for the generalized functions can be used to express the Fourier expansion coefficients directly in terms of the expansion coefficients of the Stokes scattering matrix [38,39,46]:

M_{m}^{c} (u, u^{'}) = A^{m} (u, u^{'}) + \tilde{D} A^{m} (u, u^{'}) \tilde{D}

M_{m}^{s} (u, u^{'}) = A^{m} (u, u^{'}) \tilde{D} - \tilde{D} A^{m} (u, u^{'})

where D̃ = diag{1, 1, −1, −1}. The matrix A^m(u,u′) is given by:

A^{m} (u, u^{'}) = \sum_{ℓ = m}^{2 N - 1} P_{m}^{ℓ} (u) Λ^{ℓ} P_{m}^{ℓ} (u^{'})

where

Λ^{ℓ} = (\begin{matrix} α_{1}^{ℓ} & β_{1}^{ℓ} & 0 & 0 \\ β_{1}^{ℓ} & α_{2}^{ℓ} & 0 & 0 \\ 0 & 0 & α_{3}^{ℓ} & β_{2}^{ℓ} \\ 0 & 0 & - β_{2}^{ℓ} & α_{4}^{ℓ} \end{matrix})

and

a_{1} (Θ) = \sum_{ℓ = 0}^{2 N - 1} α_{1}^{ℓ} P_{0, 0}^{ℓ} (cos Θ)

a_{2} (Θ) + a_{3} (Θ) = \sum_{ℓ = 2}^{2 N - 1} (α_{2}^{ℓ} + α_{3}^{ℓ}) P_{2, 2}^{ℓ} (cos Θ)

a_{2} (Θ) - a_{3} (Θ) = \sum_{ℓ = 2}^{2 N - 1} (α_{2}^{ℓ} - α_{3}^{ℓ}) P_{2, - 2}^{ℓ} (cos Θ)

a_{4} (Θ) = \sum_{ℓ = 0}^{2 N - 1} α_{4}^{ℓ} P_{0, 0}^{ℓ} (cos Θ)

b_{1} (Θ) = \sum_{ℓ = 2}^{2 N - 1} β_{1}^{ℓ} P_{0, 2}^{ℓ} (cos Θ)

b_{2} (Θ) = \sum_{ℓ = 2}^{2 N - 1} β_{2}^{ℓ} P_{0, 2}^{ℓ} (cos Θ) .

Here

α_{j}^{ℓ}

and

β_{j}^{ℓ}

are expansion coefficients, and a_j(Θ) and b_j(Θ) are the elements of the Stokes scattering matrix in Eq. (30).

The matrix $P_{m}^{ℓ} (u)$ occurring in Eq. (43) is defined as:

P_{m}^{ℓ} (u) = (\begin{matrix} P_{m, 0}^{ℓ} (u) & 0 & 0 & 0 \\ 0 & P_{m, +}^{ℓ} (u) & P_{m, -}^{ℓ} (u) & 0 \\ 0 & P_{m, -}^{ℓ} (u) & P_{m, +}^{ℓ} (u) & 0 \\ 0 & 0 & 0 & P_{m, 0}^{ℓ} (u) \end{matrix}) .

where

P_{m, \pm}^{ℓ} (u) = \frac{1}{2} [P_{m, - 2}^{ℓ} (u) \pm P_{m, 2}^{ℓ} (u)] .

The functions

P_{m, 0}^{ℓ} (u)

and

P_{m, \pm 2}^{ℓ} (u)

are the generalized spherical functions.

Acknowledgments

Partial support of this research was provided by the United States National Aeronautics and Space Administration (NASA) through a grant to Stevens Institute of Technology. Constructive criticisms and suggestions by three anonymous reviewers, which helped improve the quality of this paper, are greatly appreciated.

References and links

1. H. R. Gordon, “Atmospheric correction of ocean color imagery in the Earth Observation System era,” J. Geophys. Res. 102, 17081–17106 (1997) [CrossRef] .

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

3. C. Rodgers, Inverse Methods for Atmospheric Sounding (World Scientific, 2000).

4. W. Li, K. Stamnes, R. Spurr, and J. J. Stamnes, “Simultaneous retrieval of aerosols and ocean properties: A classic inverse modeling approach. II. SeaWiFS case study for the Santa Barbara channel,” Int. J. Rem. Sens. 29, 5689–5698 (2008). [CrossRef]

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

6. R. Spurr, K. Stamnes, H. Eide, W. Li, K. Zhang, and J. J. Stamnes, “Simultaneous retrieval of aerosol and ocean properties: A classic inverse modeling approach: I. Analytic Jacobians from the linearized CAO-DISORT model,” J. Quant. Spectrosc. Radiat. Transfer 104, 428–449 (2007) [CrossRef] .

7. K. Nielsen, L. Zhao, G. A. Ryzhikov, M. S. Biryulina, E. R. Sommersten, J. J. Stamnes, K. Stamnes, and J. Moan, “Retrieval of the physiological state of human skin from UV-VIS reflectance spectra: A feasibility study,” J. Photochem. Photobiol. B 93, 23–31 (2008) [CrossRef] [PubMed] .

8. D. L. Swanson, S. D. Laman, M. Biryulina, K. P. Nielsen, G. Ryzhikov, J. J. Stamnes, B. Hamre, L. Zhao, F. S. Castellana, and K. Stamnes, “Optical transfer diagnosis of pigmented lesions: a pilot study,” Skin Res. Technol. 15, 330–337 (2009). [CrossRef] [PubMed]

9. D. L. Swanson, S. D. Laman, M. Biryulina, K. P. Nielsen, G. Ryzhikov, J. J. Stamnes, B. Hamre, L. Zhao, E. Sommersten, F. S. Castellana, and K. Stamnes, “Optical transfer diagnosis of pigmented lesions,” Dermatol. Surg. 36, 1–8 (2010). [CrossRef]

10. C. N. Adams and G. W. Kattawar, “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,” Limm. Ocean. 34, 1453–1472 (1989) [CrossRef] .

11. M. I. Mishchenko, A.A. Lacis, and L. D. Travis, “Errors due to the neglect of polarization in radiance calculations for Rayleigh-scattering atmospheres,” J. Quant Spectrosc. Radiat Transfer 51, 491–510 (1994) [CrossRef] .

12. A. A. Lacis, J. Chowdhary, M. I. Mishchenko, and B. Cairns, “Modeling errors in diffuse-sky radiation: vector vs. scalar treatment,” Geophys. Res. Lett. 25, 135–138 (1998) [CrossRef] .

13. S. Chandrasekhar, Radiative Transfer (Dover Publications, 1960).

14. K. N. Liou, “A numerical experiment on Chandrasekhar’s discrete-ordinate method for radiative transfer,” J. Atmos. Sci. 30, 1303–1326 (1973) [CrossRef] .

15. K. Stamnes and R. A. Swanson, “A new look at the discrete ordinate method for radiative transfer calculations in anisotropically scattering atmospheres,” J. Atmos. Sci. 38, 387–399 (1981) [CrossRef] .

16. K. Stamnes and H. Dale, “A new look at the discrete ordinate method for radiative transfer calculations in anisotropically scattering atmospheres. II. Intensity computations,” J. Atmos. Sci. 38, 2696–2706 (1981) [CrossRef] .

17. K. Stamnes and P. Conklin, “A new multi-layer discrete ordinate approach to radiative transfer in vertically inhomogeneous atmospheres,” J. Quant. Spectrosc. Radiat. Transfer 31, 273–282 (1984) [CrossRef] .

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

19. K. P. Nielsen, L. Zhao, P. Juzenas, K. Stamnes, J. J. Stamnes, and J. Moan, “Reflectance spectra of pigmented and non-pigmented skin in the UV spectral region,” Photochem. Photobiol. 80, 450–455 (2004) [PubMed] .

20. K. I. Gjerstad, J. J. Samnes, B. Hamre, J. K. Lotsberg, B. Yan, and K. Stamnes, “Monte Carlo and discrete-ordinate simulations of irradiances in the coupled atmosphere-ocean system,” Appl. Opt. 42, 2609–2622 (2003) [CrossRef] [PubMed] .

21. K. Hestenes, K. P. Nielsen, L. Zhao, J. J. Stamnes, and K. Stamnes, “Monte Carlo and discrete-ordinate simulations of spectral radiances in the coupled air-tissue system,” Appl. Opt. 46, 2333–2350 (2007) [CrossRef] [PubMed] .

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

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

24. F. M. Schulz, K. Stamnes, 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] .

25. F. M. Schulz and K. Stamnes, “Angular distribution of the Stokes vector in a plane parallel, vertically inhomogeneous medium in the vector discrete ordinate radiative transfer (VDISORT) model,” J. Quant. Spectrosc. Radiat. Transfer 65, 609–620 (2000) [CrossRef] .

26. S. Jiang, “Radiative Transfer in the Coupled Atmosphere-Sea Ice-Ocean System with Applications in Remote Sensing,” PhD thesis, Stevens Institute of Technology (2003).

27. E. R. Sommersten, “CAO-VDISORT: A discrete ordinate method for polarized (vector) radiative transfer in a coupled system consisting of two media with different indices of refraction,” MSc thesis, University of Bergen (2005).

28. E. R. Sommersten, J. K. Lotsberg, K. Stamnes, and J. J. Stamnes, “Discrete ordinate and Monte Carlo simulations for polarized radiative transfer in a coupled system consisting of two media with different refractive indices,” J. Quant. Spectrosc. Radiat. Transfer 111, 616–633 (2010) [CrossRef] .

29. 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] .

30. J. K. Lotsberg, “Impact of particulate oceanic composition on the radiance and polarization of the natural light field,” PhD thesis, University of Bergen (2005).

31. J. K. Lotsberg and J. J. Stamnes, “Impact of particulate oceanic composition on the radiance and polarization of underwater and backscattered light,” Opt. Expr. 18, 10432–10445 (2010) [CrossRef] .

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

33. 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, 1980) [CrossRef] .

34. 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] .

35. H. Ishimoto and K. Masuda, “A Monte Carlo approach for the calculation of polarized light,” J. Quant. Spectrosc. Radiat. Transfer 72, 467–483 (2002) [CrossRef] .

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

37. P. W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Lucker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Expr. 17, 2057–2079 (2009) [CrossRef] .

38. C. E. Siewert, “On the equation of transfer relevant to the scattering of polarized light,” Astrophys. J. 245, 1080–1086 (1981) [CrossRef] .

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

40. J. F. de Haan, P. B. Bosma, and J. W. Hovenier, “The adding method for multiple scattering calculations of polarized light,” Astron. Astrophys. 183, 371–391 (1987).

41. J. M. Wallace and P. V. Hobbs, Atmospheric Science: Introductory Survey (Academic, 1977).

42. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley, 1998) [CrossRef] .

43. X. Zhou, S. Li, and K. Stamnes, “Geometrical-optics code for computing optical properties of large dielectric spheres,” Appl. Opt. 42, 4295–4306 (2003) [CrossRef] [PubMed] .

44. J. E. Hansen and L. D. Travis, “Light scattering in planetary atmospheres,” Space Sci. Rev. 16, 527–610 (1974) [CrossRef] .

45. A. A. Kokhanovsky, C. Cornet, M. Duan, C. Emde, I. L. Katsev, L. C-Labonnote, Q. Min, T. Nakajima, Y. Ota, and A.P. Prikhach, and others, “Benchmark results in vector radiative transfer,” J. Quant. Spectrosc. Radiat. Transfer 111, 1931–1946 (2010). [CrossRef] .

46. M. I. Mishchenko, “Light scattering by randomly oriented rotationally symmetric particles,” J. Opt. Soc. Am. A 8, 871–882 (1991) [CrossRef] .

47. M. I. Mischenko and L. D. Travis, “Capabilities and Limitations of a Current FORTRAN Implementation of the T-Matrix Method for Randomly Oriented, Rotationally Symmetric Scatterers,” J. Quant. Spectrosc. Radiat. Transfer 60, 309–324 (1998) [CrossRef] .

48. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

49. M. Born and E. Wolf, Principles of Optics: 7th edition (Cambridge University, 2002).

50. J. F. de Haan, P. B. Bosma, and J. W. Hovenier, “The adding method for multiple scattering calculations of polarized light,” Astron. Astrophys. 183, 371–391 (1987).

51. 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).

Comparison of discrete ordinate and Monte Carlo simulations of polarized radiative transfer in two coupled slabs with different refractive indices

Abstract

1. Introduction

2. Radiative transfer theory

2.1. General radiative transfer

2.2. Types of scattering & input/output

Aerosol particles

Cloud particles

2.3. The interface between the two media

3. Radiative transfer models

3.1. C-VDISORT

3.2. C-PMC

4. Comparison of C-VDISORT and C-PMC RT simulations

4.1. Mie Scattering in the Rayleigh limit for a coupled system

4.2. C-VDISORT and C-PMC vs benchmark – aerosol layer – reflection

4.3. C-VDISORT and C-PMC vs benchmark – aerosol layer – transmission

4.4. C-VDISORT and C-PMC vs benchmark – cloud layer – reflection

4.5. C-VDISORT and C-PMC vs benchmark – cloud layer – transmission

4.6. C-VDISORT vs C-PMC – aerosol particles – coupled case

4.7. C-VDISORT vs C-PMC – aerosol/cloud particles – coupled case

5. Conclusions

6. Appendix A

6.1. Stokes vector representation I⃗_S = [I, Q, U, V]^T used in C-PMC

6.2. Stokes vector representation I⃗ = [I_||, I_⊥, U, V]^T used in C-VDISORT

6.3. Generalized Spherical Functions

Acknowledgments

References and links

Cited By

Figures (9)

Equations (54)

Optics Express

Abstract

1. Introduction

2. Radiative transfer theory

2.1. General radiative transfer

2.2. Types of scattering & input/output

Aerosol particles

Cloud particles

2.3. The interface between the two media

3. Radiative transfer models

3.1. C-VDISORT

3.2. C-PMC

4. Comparison of C-VDISORT and C-PMC RT simulations

4.1. Mie Scattering in the Rayleigh limit for a coupled system

4.2. C-VDISORT and C-PMC vs benchmark – aerosol layer – reflection

4.3. C-VDISORT and C-PMC vs benchmark – aerosol layer – transmission

4.4. C-VDISORT and C-PMC vs benchmark – cloud layer – reflection

4.5. C-VDISORT and C-PMC vs benchmark – cloud layer – transmission

4.6. C-VDISORT vs C-PMC – aerosol particles – coupled case

4.7. C-VDISORT vs C-PMC – aerosol/cloud particles – coupled case

5. Conclusions

6. Appendix A

6.1. Stokes vector representation I⃗S = [I, Q, U, V]T used in C-PMC

6.2. Stokes vector representation I⃗ = [I||, I⊥, U, V]T used in C-VDISORT

6.3. Generalized Spherical Functions

Acknowledgments

References and links

Cited By

Figures (9)

Equations (54)

Optics Express

6.1. Stokes vector representation I⃗_S = [I, Q, U, V]^T used in C-PMC

6.2. Stokes vector representation I⃗ = [I_||, I_⊥, U, V]^T used in C-VDISORT