Markov chain formalism for polarized light transfer in plane-parallel atmospheres, with numerical comparison to the Monte Carlo method

Feng Xu; Anthony B. Davis; Robert A. West; Larry W. Esposito

doi:10.1364/OE.19.000946

1. Introduction and overview

The polarization state of atmospheric radiation potentially contains abundant information about aerosol properties. A computationally efficient and accurate method for vector radiative transfer is therefore important for practical retrieval algorithms. For this purpose, the Markov chain method was developed to calculate the diffusely reflected intensity from a vertically homogeneous and inhomogeneous atmosphere [1–2]. Application to Venus’ atmosphere indicates its higher computation efficiency than that of adding/doubling method [3–4] by one magnitude [2]. Moreover, the Markov chain method retains (i) the advantage of the adding/doubling method of easy physical interpretation and high numerical accuracy in calculating the atmosphere of large optical depth, (ii) the advantage of the successive orders of scattering method in separating the contribution of a specific scattering order numerically [5–7] or analytically [8], and (iii) the advantage of the discrete source method [9–11] for the internal radiation field calculation. These merits form the motivation of the current research to extend the Markov chain method to polarized radiative transfer in a vertically inhomogeneous atmosphere.

This study is organized as follows. In order to ease the derivation of vector radiative transfer formalism, Section 2 gives an introduction to the Markov chain formalism for the scalar case in a physical way. The formalism differs to that in Ref [1]. in two respects: (1) inclusion of transmission, and (2) consideration of a vertically inhomogeneous atmosphere. Then as a generalization, Section 3 extends the Markov chain formalism to vector radiative transfer in a vertically inhomogeneous plane-parallel atmosphere. A comparison of numerical results calculated by the Markov chain method to the counterparts from the Monte Carlo method is presented in Section 4, followed by a summary and outlook in Section 5.

2. Scalar radiative transfer

Assume that a vertically inhomogeneous plane-parallel atmosphere of optical depth τ ₀ is illuminated by a monochromatic solar flux πF ₀. The atmosphere is divided into N sublayers so that each one has its optical thickness Δτ_n = τ_n ₊₁−τ_n (1 ≤ n ≤ N with τ_N ₊₁ = τ ₀), single scattering albedo ω₀(n) and phase function P(cosΘ,n) which is normalized in such a way that its integral over all directions is 4π. The scattering angle Θ is defined by

\cos Θ = u u_{0} + {(1 - u^{2})}^{1 / 2} {(1 - u_{0}^{2})}^{1 / 2} \cos (ϕ - ϕ_{0}),

where u ₀ = cosθ ₀ and u = cosθ denote the direction of the incident and scattered light, respectively, with respect to the downward increasing optical depth,

ϕ_{0}

and ϕ describe the azimuthal plane of the incident and diffuse light.

When polarization is neglected, the phase function is even function about ( $ϕ - ϕ_{0}$ ) and thus can be expanded into the Fourier series components (P ⁽ ^m ⁾) for cosine modes only:

P (\cos Θ, n) = P (u, u_{0}, ϕ - ϕ_{0}; n) = P^{(0)} (u, u_{0}; n) + 2 \sum_{m = 1}^{\infty} P^{(m)} (u, u_{0}; n) \cos m (ϕ - ϕ_{0}) .

To further specify the scattering type to be reflection or transmission according to the sign of u and specify the incidence to be upwelling or downwelling according to the sign of u ₀, we define μ = |u| and μ ₀ = |u ₀|, getting the following four regimes [3]:

P_{r} (μ, μ_{0}, ϕ - ϕ_{0}; n) = P (- μ, μ_{0}, ϕ - ϕ_{0}; n),

P_{t} (μ, μ_{0}, ϕ - ϕ_{0}; n) = P (μ, μ_{0}, ϕ - ϕ_{0}; n),

P_{r}^{*} (μ, μ_{0}, ϕ - ϕ_{0}; n) = P (μ, - μ_{0}, ϕ - ϕ_{0}; n),

P_{t}^{*} (μ, μ_{0}, ϕ - ϕ_{0}; n) = P (- μ, - μ_{0}, ϕ - ϕ_{0}; n),

where the subscripts “r” and “t” of the phase function denote the reflection and transmission, respectively for the illumination from above while the asterisk (*) in the superscript denote the illumination from below. Accordingly, the notations for the Fourier series components become

P_{r}^{(m)} (μ, μ_{0}; n) = P^{(m)} (- μ, μ_{0}; n),

P_{t}^{(m)} (μ, μ_{0}; n) = P^{(m)} (μ, μ_{0}; n),

P_{r}^{(m), *} (μ, μ_{0}; n) = P^{(m)} (μ, - μ_{0}; n),

P_{t}^{(m), *} (μ, μ_{0}; n) = P^{(m)} (- μ, - μ_{0}; n) .

From Eq. (1), we have $P_{r / t}^{*} (μ, μ_{0}, ϕ - ϕ_{0}; n)$ = $P_{r / t} (μ, μ_{0}, ϕ - ϕ_{0}; n)$ for scalar radiative transfer so that $P_{r / t}^{(m), *} (μ, μ_{0}; n)$ = $P_{r / t}^{(m)} (μ, μ_{0}; n)$ . In the same way the radiation field can be expanded as a series:

I (τ, μ, ϕ) = I^{(0)} (τ, μ) + 2 \sum_{m = 1}^{\infty} I^{(m)} (τ, μ) \cos m (ϕ - ϕ_{0}) .

Both the diffusely reflected field I _r(τ = 0, μ, ϕ) and transmitted field I _t(τ = τ ₀, μ, ϕ) are resolved once I ⁽⁰⁾(τ, μ) and I ⁽ ^m ⁾(τ, μ) are determined.

2.1 Initial scattering in the atmosphere

We begin the calculation of multiple scattering with the intensity of the first order scattered light. This is the initial state for the Markov chain. Taking into account the exponential attenuation of the incident solar light (direction cosine μ ₀) and the outgoing scattered light (direction cosine μ_i), the contribution of first-order scattering to the diffusely reflected intensity from the upper boundary of the n ^th layer (τ = τ_n) has the Fourier series $U_{(n, i), 0, r}^{(m)}$ /2, where

\begin{array}{l} U_{(n, i), 0, r}^{(m)} = w_{i} F_{0} \exp (- \frac{τ_{n}}{μ_{0}}) \int_{0}^{Δ τ_{n}} \frac{d x}{μ_{i}} \exp (- \frac{x}{μ_{0}}) ω_{0} (n) \frac{P_{r}^{(m)} (μ_{i}, μ_{0}; n)}{2} \exp (- \frac{x}{μ_{i}}) \\ = w_{i} F_{0} \exp (- \frac{τ_{n}}{μ_{0}}) \frac{μ_{0}}{μ_{i} + μ_{0}} ω_{0} (n) \frac{P_{r}^{(m)} (μ_{i}, μ_{0}; n)}{2} {1 - \exp [- Δ τ_{n} (\frac{1}{μ_{0}} + \frac{1}{μ_{i}})]} \end{array}

and the single-scattering contribution to the diffused transmitted intensity through the lower boundary of the n^th layer (τ = τ_n ₊₁) has the Fourier series

U_{(n, i), 0, t}^{(m)}

/2, where

\begin{array}{l} U_{(n, i), 0, t}^{(m)} = w_{i} F_{0} \exp (- \frac{τ_{n}}{μ_{0}}) \int_{0}^{Δ τ_{n}} \frac{d x}{μ_{i}} \exp (- \frac{x}{μ_{0}}) ω_{0} (n) \frac{P_{t}^{(m)} (μ_{i}, μ_{0}; n)}{2} \exp (- \frac{Δ τ_{n} - x}{μ_{i}}) \\ = w_{i} F_{0} \exp (- \frac{τ_{n}}{μ_{0}}) \frac{μ_{0}}{μ_{i} - μ_{0}} ω_{0} (n) \frac{P_{t}^{(m)} (μ_{i}, μ_{0}; n)}{2} {\exp (- \frac{Δ τ_{n}}{μ_{i}}) - \exp (- \frac{Δ τ_{n}}{μ_{0}})} . \end{array}

Note that for integrating the contribution of the incident light in all directions μ_i to the scattered intensity in one direction in the next order of scattering, a quadrature weight w_i is introduced for μ_i on the interval 0 ≤ μ_i ≤ 1.

Assuming the emergent single scattering intensity from the boundary to be produced by uniform source distribution in each sublayer, the source intensity can be calculated through dividing Eqs. (12)-(13) by an average attenuation factor c ₍ _n,i ₎ along the path Δτ_n/μ_i, namely,

Π_{(n, i), 0, r / t}^{(m)} = U_{(n, i), 0, r / t}^{(m)} / c_{(n, i)},

where

c_{(n, i)} = \frac{1}{Δ τ_{n}} \int_{0}^{Δ τ_{n}} \exp (- \frac{x}{μ_{i}}) d x = \frac{μ_{i}}{Δ τ_{n}} [1 - \exp (- \frac{Δ τ_{n}}{μ_{i}})] .

2.2 Intermediate scattering in the atmosphere

To analyze the multiple scattering, we assume the light last scattered from the layer n in the direction μ_i travels upward to the layer n' for the next order of scattering. The reduced flux for scattering in n' ^th layer is:

W_{(n^{'}, i), (n, i)} = [\frac{1}{Δ τ_{n}} \int_{0}^{Δ τ_{n}} \exp (- \frac{x}{μ_{i}}) d x] \exp (- \frac{τ_{n} - τ_{n^{'} + 1}}{μ_{i}}) [\int_{0}^{Δ τ_{n^{'}}} \exp (- \frac{x}{μ_{i}}) \frac{d x}{μ_{i}}],

where on the basis of the uniform source distribution in each layer the first term stands for the intensity of the light leaving the upper boundary of n ^th layer of optical thickness Δτ_n, the second term indicates the remaining light intensity reaching the lower boundary of n' ^th layer after traveling an optical depth (τ_n−τ_n' ₊₁), and the last term indicates the attenuated light flux in the n' ^th layer. Equation (16) has the following analytical form:

W_{(n^{'}, i), (n, i)} = \frac{μ_{i}}{Δ τ_{n}} [1 - \exp (- \frac{Δ τ_{n}}{μ_{i}})] \exp (- \frac{τ_{n} - τ_{n^{'} + 1}}{μ_{i}}) [1 - \exp (- \frac{Δ τ_{n^{'}}}{μ_{i}})] .

In a similar way, we can calculate the attenuated light flux in the n' ^th layer for the downward travelling light (n < n')

W_{(n^{'}, i), (n, i)} = \frac{μ_{i}}{Δ τ_{n}} [1 - \exp (- \frac{Δ τ_{n}}{μ_{i}})] \exp (- \frac{τ_{n^{'}} - τ_{n + 1}}{μ_{i}}) [1 - \exp (- \frac{Δ τ_{n^{'}}}{μ_{i}})] .

In addition to the upward and downward propagation, partial light flux remains in the same layer (n = n') for scattering:

W_{(n, i), (n, i)} = 1 - \frac{1}{Δ τ_{n}} \int_{0}^{Δ τ_{n}} \exp (- \frac{x}{μ_{i}}) d x = 1 - \frac{μ_{i}}{Δ τ_{n}} [1 - \exp (- \frac{Δ τ_{n}}{μ_{i}})] .

Knowing the reduced flux, the scattered intensity from n' ^th layer in all directions can be calculated from its Fourier series which satisfies

(2 - δ_{0 m}) I_{}^{(m)} = \frac{1}{2} A_{}^{(m)} B_{}^{(m)},

where the B-matrix

B_{(n, i)}^{(m)}

is the intensity of the incident light from the n ^th layer in direction μ_i (for the second order of scattering

B_{}^{(m)}

= [

Π_{r, 0}^{(m)}, Π_{t, 0}^{(m)}

]^T) and the A-matrix

A_{(n^{'}, j), (n, i)}^{(m)}

describes the scattered intensity from the n'-th layer and direction μ_j, which is contributed by the incident light from n ^th layer and direction μ_i,

A_{(n^{'}, j), (n, i)}^{(m)} = ω_{0} (n^{'}) \frac{P^{(m)} (μ_{j}, μ_{i}; n^{'})}{2} W_{(n^{'}, i), (n, i)} .

For upwelling light

P_{r}^{(m)}

and

P_{t}^{(m)}

are used as

P^{(m)}

in above equation while for the downwelling light

P_{r}^{(m), *}

and

P_{t}^{(m), *}

are used. The scattered light becomes the incident light in the next order of scattering. For integrating the contribution of the incident light in all directions to the next order of scattering, a Gaussian quadrature weight w_j is introduced to define a transition factor

Q_{(n^{'}, j), (n, i)}^{(m)}

:

Q_{(n^{'}, j), (n, i)}^{(m)} = w_{j} A_{(n^{'}, i), (n, i)}^{(m)} .

Note that Q represents the transition of light between two states due to one order of scattering. It forms the basis for a later description of all intermediate multiple scattering processes.

2.3 Scattering to the boundary

The emergent intensities from the top of the atmosphere (TOA, or τ = 0) and from the bottom of the atmosphere (BOA, or τ = τ ₀) are contributed by both downward- and upward-traveling light. Assuming they come from layer n' in direction j, and are redirected to emerge in direction μ _e through scattering, we then analyze their contribution to the radiation fields at the BOA and TOA.

2.3.1 Downwelling light

Downwelling light from n' ^th layer partially leaves the TOA through diffuse reflection and partially leaves the BOA through diffuse transmission. The diffuse scattering occurs in all layers n with n ≥ n' so that the Fourier series of the emergent intensities are:

R_{e, (n^{'}, j)}^{(m)} = \sum_{n = n^{'}}^{N} R_{(n, e), (n^{'}, j)}^{(m)}

for diffuse reflection at the TOA and

T_{e, (n^{'}, j)}^{(m)} = \sum_{n = n^{'}}^{N} T_{(n, e), (n^{'}, j)}^{(m)}

for diffuse transmission at the BOA, where the subscripts of R and T on the left hand side of the above two equations mean the light leaving the atmosphere in the direction μ _e is contributed by the scattering of the light coming from n' ^th layer in direction μ_j, and the subscript (n, e) of each series term on the right hand side of the two equations means the scattering into the direction μ _e occurs in n ^th layer.

The contribution of scattering in n ^th layer to the emergent intensities is calculated by integrating over τ bounded by τ_n and τ_n ₊₁. In terms of the Fourier series we have

R_{(n, e), (n^{'}, j)}^{(m)} = \int_{τ_{n}}^{τ_{n + 1}} I_{^{r, (n, μ_{e}), (n^{'}, μ_{j})}}^{(m), d} (τ) \frac{d τ}{μ_{e}},

T_{(n, e), (n^{'}, j)}^{(m)} = \int_{τ_{n}}^{τ_{n + 1}} I_{^{t, (n, μ_{e}), (n^{'}, μ_{j})}}^{(m), d} (τ) \frac{d τ}{μ_{e}},

where the kernels of the above integrals

I_{r, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), d}

and

I_{t, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), d}

represent the contribution of the scattering of the downwelling (denoted by the superscript “d”) incident light from n' ^th layer to the diffusely emergent light through reflection (denoted by the subscript “r”) and transmission (denoted by the subscript “t”) from/through n ^th layer, respectively. Denoting

I_{0}^{(m), d}

to be the Fourier series vector for all states of the incident light intensity from n' ^th layer, multiplication with matrices

I_{r}^{(m), d} / 2

and

I_{t}^{(m), d} / 2

gives the Fourier series of the reflected and transmitted fields, respectively, where

I_{r}^{(m), d}

and

I_{t}^{(m), d}

are evaluated by

\begin{array}{l} I_{r, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), d} (τ) = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} \int_{τ_{n^{'}}}^{τ_{n^{'} + 1}} \frac{d x}{μ_{j}} \exp (- \frac{τ}{μ_{e}}) ω_{0} (n) \frac{P_{r}^{(m)} (μ_{e}, μ_{j}; n)}{2} \exp (- \frac{τ - x}{μ_{j}}) \\ = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} ω_{0} (n) \frac{P_{r}^{(m)} (μ_{e}, μ_{j}; n)}{2} {\exp [(\frac{τ_{n^{'} + 1} - τ}{μ_{j}} - \frac{τ}{μ_{e}})] - \exp [(\frac{τ_{n^{'}} - τ}{μ_{j}} - \frac{τ}{μ_{e}})]}, \end{array}

\begin{array}{l} I_{t, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), d} (τ) = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} \int_{τ_{n^{'}}}^{τ_{n^{'} + 1}} \frac{d x}{μ_{j}} \exp (- \frac{τ_{0} - τ}{μ_{e}}) ω_{0} (n) \frac{P_{t}^{(m)} (μ_{e}, μ_{j}; n)}{2} \exp (- \frac{τ - x}{μ_{j}}) \\ = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} ω_{0} (n) \frac{P_{t}^{(m)} (μ_{e}, μ_{j}; n)}{2} {\exp [(\frac{τ_{n^{'} + 1} - τ}{μ_{j}} - \frac{τ_{0} - τ}{μ_{e}})] - \exp [(\frac{τ_{n^{'}} - τ}{μ_{j}} - \frac{τ_{0} - τ}{μ_{e}})]} . \end{array}

Substitution of Eq. (27) and (28) into Eq. (25) and (26) respectively gives the analytical solution for

R_{(n, e), (n^{'}, j)}^{(m)}

and

T_{(n, e), (n^{'}, j)}^{(m)}

:

\begin{array}{l} R_{(n, e), (n^{'}, j)}^{(m)} = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} \frac{μ_{j}}{μ_{e} + μ_{j}} ω_{0} (n) \frac{P_{r}^{(m)} (μ_{e}, μ_{j}; n)}{2} {[\exp (\frac{τ_{n^{'} + 1} - τ_{n}}{μ_{j}} - \frac{τ_{n}}{μ_{e}}) - \exp (\frac{τ_{n^{'}} - τ_{n}}{μ_{j}} - \frac{τ_{n}}{μ_{e}})] \\ - [\exp (\frac{τ_{n^{'} + 1} - τ_{n + 1}}{μ_{j}} - \frac{τ_{n + 1}}{μ_{e}}) - \exp (\frac{τ_{n^{'}} - τ_{n + 1}}{μ_{j}} - \frac{τ_{n + 1}}{μ_{e}})]}, \end{array}

\begin{array}{l} T_{(n, e), (n^{'}, j)}^{(m)} = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} \exp (- \frac{τ_{0}}{μ_{e}}) \frac{μ_{j}}{μ_{e} - μ_{j}} ω_{0} (n) \frac{P_{t}^{(m)} (μ_{e}, μ_{j}; n)}{2} {[\exp (\frac{τ_{n^{'} + 1} - τ_{n}}{μ_{j}} + \frac{τ_{n}}{μ_{e}}) - \exp (\frac{τ_{n^{'}} - τ_{n}}{μ_{j}} + \frac{τ_{n}}{μ_{e}})] \\ - [\exp (\frac{τ_{n^{'} + 1} - τ_{n + 1}}{μ_{j}} + \frac{τ_{n + 1}}{μ_{e}}) - \exp (\frac{τ_{n^{'}} - τ_{n + 1}}{μ_{j}} + \frac{τ_{n + 1}}{μ_{e}})]} . \end{array}

In the particular case that diffuse scattering occurs in the layer n = n', Eq. (25) and (26) become

R_{(n, e), (n, j)}^{(m)} = \int_{0}^{Δ τ_{n}} I_{r, (n, μ_{e}), (n, μ_{j})}^{(m), d} (τ) \exp (- \frac{τ_{n + 1} - τ}{μ_{e}}) \frac{d τ}{μ_{e}}

and

T_{(n, e), (n, j)}^{(m)} = \int_{0}^{Δ τ_{n}} I_{t, (n, μ_{e}), (n, μ_{j})}^{(m), d} (τ) \exp (- \frac{τ_{0} - τ_{n} - τ}{μ_{e}}) \frac{d τ}{μ_{e}},

respectively, where

\begin{array}{l} I_{r, (n, μ_{e}), (n, μ_{j})}^{(m), d} (τ) = \frac{1}{2} \frac{1}{Δ τ_{n}} \int_{0}^{τ} \frac{d x}{μ_{e}} \exp (- \frac{τ}{μ_{e}}) ω_{0} (n) \frac{P_{r}^{(m)} (μ_{e}, μ_{j}; n)}{2} \exp (- \frac{x}{μ_{j}}) \\ = \frac{1}{2} \frac{1}{Δ τ_{n}} \frac{μ_{j}}{μ_{j} + μ_{e}} ω_{0} (n) \frac{P_{r}^{(m)} (μ_{e}, μ_{j}; n)}{2} {1 - \exp [- τ (\frac{1}{μ_{j}} + \frac{1}{μ_{e}})]} \end{array}

and

\begin{array}{l} I_{t, (n, μ_{e}), (n, μ_{j})}^{(m), d} (τ) = \frac{1}{2} \frac{1}{Δ τ_{n}} \int_{0}^{τ} \frac{d x}{μ_{e}} \exp (- \frac{τ - x}{μ_{e}}) ω_{0} (n) \frac{P_{t}^{(m)} (μ_{e}, μ_{j}; n)}{2} \exp (- \frac{x}{μ_{j}}) \\ = \frac{1}{2} \frac{1}{Δ τ_{n}} \frac{μ_{j}}{μ_{j} - μ_{e}} ω_{0} (n) \frac{P_{t}^{(m)} (μ_{e}, μ_{j}; n)}{2} [\exp (- \frac{τ}{μ_{j}}) - \exp (- \frac{τ}{μ_{e}})] . \end{array}

Substitution of Eqs. (33) and (34) into Eqs. (31) and (32) gives analytical expressions for diffuse reflection at the TOA and diffuse transmission at the BOA respectively, namely,

\begin{array}{l} R_{(n, e), (n, j)}^{(m)} = \\ \frac{1}{2} \frac{1}{Δ τ_{n}} \frac{1}{μ_{j} + μ_{e}} ω_{0} (n) \frac{P_{r}^{(m)} (μ_{e}, μ_{j}; n)}{2} {μ_{e} [\exp (- \frac{τ_{n}}{μ_{e}}) - \exp (- \frac{τ_{n + 1}}{μ_{e}})] - μ_{j} \exp (- \frac{τ_{n + 1}}{μ_{e}}) [1 - \exp (- \frac{Δ τ_{n}}{μ_{j}})]}, \end{array}

\begin{array}{l} T_{(n, e), (n, j)}^{(m)} = \\ \frac{1}{2} \frac{1}{Δ τ_{n}} \frac{1}{μ_{j} - μ_{e}} ω_{0} (n) \frac{P_{t}^{(m)} (μ_{e}, μ_{j}; n)}{2} \exp (- \frac{τ_{0} - τ_{n}}{μ_{e}}) {\frac{μ_{j} μ_{e}}{μ_{j} - μ_{e}} [\exp (- \frac{μ_{j} - μ_{e}}{μ_{j} μ_{e}} Δ τ_{n}) - 1] - Δ τ_{n}} . \end{array}

2.3.2 Upwelling light

Upwelling light from n' ^th layer leaves the TOA through diffuse transmission and partially leaves the BOA through diffuse reflection. The diffuse scattering occurs in all layers n with n ≤ n' so that the Fourier series of the emergent intensities are,

R_{e, (n^{'}, j)}^{(m)} = \sum_{n = 1}^{n'} R_{(n, e), (n^{'}, j)}^{(m)}

and

T_{e, (n^{'}, j)}^{(m)} = \sum_{n = 1}^{n'} T_{(n, e), (n^{'}, j)}^{(m)} .

Integrating the scattering over τ of the n ^th layer bounded by τ_n and τ_n ₊₁ gives

R_{(n, e), (n^{'}, j)}^{(m)} = \int_{τ_{n}}^{τ_{n + 1}} I_{t, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), u} (τ) \frac{d τ}{μ_{e}}

and

T_{(n, e), (n^{'}, j)}^{(m)} = \int_{τ_{n}}^{τ_{n + 1}} I_{r, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), u} (τ) \frac{d τ}{μ_{e}},

where the kernels of the above integrals

I_{r, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), u}

and

I_{t, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), u}

represent the contribution of the scattering of the upwelling (denoted by the superscript “u”) incident light from n' ^th layer to the diffusely emergent light through transmission (denoted by the subscript “t”) and reflection (denoted by the subscript “r), from/through the n ^th layer, respectively. Denoting

I_{0}^{(m), u}

to be the Fourier series vector for all states of the incident light intensity from n' ^th layer, multiplication with the matrices

I_{r}^{(m), u} / 2

and

I_{t}^{(m), u} / 2

gives the Fourier series of the reflected and transmitted fields, respectively, where

I_{r}^{(m), u}

and

I_{t}^{(m), u}

are evaluated by

\begin{array}{l} I_{t, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), u} (τ) = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} \int_{τ_{n^{'}}}^{τ_{n^{'} + 1}} \frac{d x}{μ_{j}} \exp (- \frac{τ}{μ_{e}}) ω_{0} (n) \frac{P_{t}^{(m), *} (μ_{e}, μ_{j}; n)}{2} \exp (- \frac{x - τ}{μ_{j}}) \\ = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} ω_{0} (n) \frac{P_{t}^{(m), *} (μ_{e}, μ_{j}; n)}{2} [\exp (\frac{τ - τ_{n^{'}}}{μ_{j}} - \frac{τ}{μ_{e}}) - \exp (\frac{τ - τ_{n^{'} + 1}}{μ_{j}} - \frac{τ}{μ_{e}})] \end{array}

\begin{array}{l} I_{r, (n, μ_{e}), (n^{'}, μ_{j})}^{(m), u} (τ) = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} \int_{τ_{n^{'}}}^{τ_{n^{'} + 1}} \frac{d x}{μ_{j}} \exp (- \frac{τ_{0} - τ}{μ_{e}}) ω_{0} (n) \frac{P_{r}^{(m), *} (μ_{e}, μ_{j}; n)}{2} \exp (- \frac{x - τ}{μ_{j}}) \\ = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} ω_{0} (n) \frac{P_{r}^{(m), *} (μ_{e}, μ_{j}; n)}{2} [\exp (\frac{τ - τ_{n^{'}}}{μ_{j}} - \frac{τ_{0} - τ}{μ_{e}}) - \exp (\frac{τ - τ_{n^{'} + 1}}{μ_{j}} - \frac{τ_{0} - τ}{μ_{e}})] . \end{array}

Substitution of Eq. (41) and (42) into Eq. (39) and (40) respectively gives the analytical solution for

R_{(n, e), (n^{'}, j)}^{(m)}

and

T_{(n, e), (n^{'}, j)}^{(m)}

:

\begin{array}{l} R_{(n, e), (n^{'}, j)}^{(m)} = \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} \frac{μ_{j}}{μ_{e} - μ_{j}} ω_{0} (n) \frac{P_{t}^{(m), *} (μ_{e}, μ_{j}; n)}{2} {[\exp (\frac{τ_{n + 1} - τ_{n^{'}}}{μ_{j}} - \frac{τ_{n + 1}}{μ_{e}}) - \exp (\frac{τ_{n + 1} - τ_{n^{'} + 1}}{μ_{j}} - \frac{τ_{n + 1}}{μ_{e}})] \\ - [\exp (\frac{τ_{n} - τ_{n^{'}}}{μ_{j}} - \frac{τ_{n}}{μ_{e}}) - \exp (\frac{τ_{n} - τ_{n^{'} + 1}}{μ_{j}} - \frac{τ_{n}}{μ_{e}})]} \end{array}

\begin{array}{l} T_{(n, e), (n^{'}, j)}^{(m)} = \\ \frac{1}{2} \frac{1}{Δ τ_{n^{'}}} \frac{μ_{j}}{μ_{e} + μ_{j}} \exp (- \frac{τ_{0}}{μ_{e}}) ω_{0} (n) \frac{P_{r}^{(m), *} (μ_{e}, μ_{j}; n)}{2} {[\exp (\frac{τ_{n + 1} - τ_{n^{'}}}{μ_{j}} + \frac{τ_{n + 1}}{μ_{e}}) - \exp (\frac{τ_{n + 1} - τ_{n^{'} + 1}}{μ_{j}} + \frac{τ_{n + 1}}{μ_{e}})] \\ - [\exp (\frac{τ_{n} - τ_{n^{'}}}{μ_{j}} + \frac{τ_{n}}{μ_{e}}) - \exp (\frac{τ_{n} - τ_{n^{'} + 1}}{μ_{j}} + \frac{τ_{n}}{μ_{e}})]} . \end{array}

In the case where diffuse scattering occurs in the layer (n = n'), Eq. (39) and (40) become

R_{(n, e), (n, j)}^{(m)} = \int_{0}^{Δ τ_{n}} I_{t, (n, μ_{e}), (n, μ_{j})}^{(m), d} (τ) \exp (- \frac{τ_{n}}{μ_{e}}) \frac{d τ}{μ_{e}}

and

T_{(n, e), (n, j)}^{(m)} = \int_{0}^{Δ τ_{n}} I_{r, (n, μ_{e}), (n, μ_{j})}^{(m), d} (τ) \exp (- \frac{τ_{0} - τ_{n + 1}}{μ_{e}}) \frac{d τ}{μ_{e}},

respectively. Substitution of Eq. (33) into Eq. (46) and Eq. (34) into Eq. (45) gives

\begin{array}{l} R_{(n, e), (n, j)}^{(m)} = \\ \frac{1}{2} \frac{1}{Δ τ_{n}} \frac{1}{μ_{j} - μ_{e}} ω_{0} (n) \frac{P_{t}^{(m), *} (μ_{e}, μ_{j}; n)}{2} {μ_{j} \exp (- \frac{τ_{n}}{μ_{e}}) [1 - \exp (- \frac{Δ τ_{n}}{μ_{j}})] - μ_{e} [\exp (- \frac{τ_{n}}{μ_{e}}) - \exp (- \frac{τ_{n + 1}}{μ_{e}})]} \end{array}

and

\begin{array}{l} T_{(n, e), (n, j)}^{(m)} = \\ \frac{1}{2} \frac{1}{Δ τ_{n}} \frac{μ_{j}}{μ_{j} + μ_{e}} ω_{0} (n) \frac{P_{r}^{(m), *} (μ_{e}, μ_{j}; n)}{2} \exp (- \frac{τ_{0} - τ_{n + 1}}{μ_{e}}) (\frac{Δ τ_{n}}{μ_{j}} + \frac{μ_{e}}{μ_{j} + μ_{e}} {\exp [- (\frac{1}{μ_{j}} + \frac{1}{μ_{e}}) Δ τ_{n}] - 1}) . \end{array}

2.4 Radiative transfer through multiple scattering

On the basis of the Π vector and the Q, R, and T matrices, the radiative transfer through different orders of scattering can be achieved in matrix form after invoking the orthogonality relation [12] for the trigonometric functions, namely,

(1 + δ_{0 m}) I_{r}^{(m)} (τ = 0) = R_{eI}^{(m)} Π_{I0}^{(m)} + R_{eJ}^{(m)} Q_{JI}^{(m)} Π_{I0}^{(m)} + R_{eK}^{(m)} Q_{KJ}^{(m)} Q_{JI}^{(m)} Π_{I0}^{(m)} + R_{eL}^{(m)} Q_{LK}^{(m)} Q_{KJ}^{(m)} Q_{JI}^{(m)} Π_{I0}^{(m)} + \dots

(1 + δ_{0 m}) I_{t}^{(m)} (τ = τ_{0}) = T_{eI}^{(m)} Π_{I0}^{(m)} + T_{eJ}^{(m)} Q_{JI}^{(m)} Π_{I0}^{(m)} + T_{eK}^{(m)} Q_{KJ}^{(m)} Q_{JI}^{(m)} Π_{I0}^{(m)} + T_{eL}^{(m)} Q_{LK}^{(m)} Q_{KJ}^{(m)} Q_{JI}^{(m)} Π_{I0}^{(m)} + \dots

where “I”, “J”, “K”, and “L”, in the subscripts denote the state of the light described by the layer number (n) and the propagating direction (μ _i); “0” denotes the initial light state in the atmosphere and “e” denotes the light leaving the atmosphere in direction μ _e. Change of the state of light from “I” to “J” is represented in the sequence “JI”. The first term on the right side of the above equation represents the diffusely emergent light experiencing two scattering events, the second term represents the emergent light experiencing three scattering events, and so on. In the form of an infinite series summation, all orders of multiple scattering are accounted for. If a specific order of scattering is not explicitly expressed and counted, Eqs. (49)-(50) have an equivalent but more compact form:

(2 - δ_{0 m}) I_{r}^{(m)} (τ = 0) = R^{(m)} {(E - Q^{(m)})}^{- 1} Π_{0}^{(m)},

(2 - δ_{0 m}) I_{t}^{(m)} (τ = τ_{0}) = T^{(m)} {(E - Q^{(m)})}^{- 1} Π_{0}^{(m)},

where E is the identity matrix.

(2.5) Total radiation field

The above description has not included the contribution of single-scattered light to the total emergent intensity. For a vertically inhomogeneous atmosphere, the contribution of single scattering to diffuse reflection and transmission can be analytically calculated by summarizing the contribution of scattering of incident solar light in all sublayers, namely,

I_{r / t, s} = \sum_{n = 1}^{N} I_{r / t, s, n},

where the contribution to the diffuse reflection by each layer is

\begin{array}{l} I_{r, s, n} (τ = 0, μ_{e}, ϕ_{e}) = \frac{ω_{0} (n)}{4} \int_{τ_{n}}^{τ_{n + 1}} \frac{d x}{μ_{e}} \exp (- \frac{x}{μ_{e}}) P_{r} (μ_{e}, ϕ_{e}; μ_{0}, ϕ_{0}, n) \exp (- \frac{x}{μ_{0}}) \\ = \frac{ω_{0} (n)}{4} \frac{μ_{0}}{μ_{e} + μ_{0}} P_{r} (μ_{e}, ϕ_{e}; μ_{0}, ϕ_{0}, n) {\exp [- τ_{n} (\frac{1}{μ_{e}} + \frac{1}{μ_{0}})] - \exp [- τ_{n + 1} (\frac{1}{μ_{e}} + \frac{1}{μ_{0}})]}, \end{array}

and the contribution to the diffuse transmission by each layer is

\begin{array}{l} I_{t, s, n} (τ = τ_{0}, μ_{e}, ϕ_{e}) = \frac{ω_{0} (n)}{4} \int_{τ_{n}}^{τ_{n + 1}} \frac{d x}{μ_{e}} \exp (- \frac{τ_{0} - x}{μ_{e}}) P_{t} (μ_{e}, ϕ_{e}; μ_{0}, ϕ_{0}, n) \exp (- \frac{x}{μ_{0}}) \\ = \frac{ω_{0} (n)}{4} \frac{μ_{0}}{μ_{0} - μ_{e}} P_{t} (μ_{e}, ϕ_{e}; μ_{0}, ϕ_{0}, n) [\exp (- \frac{τ_{0} - τ_{n + 1}}{μ_{e}}) \exp (- \frac{τ_{n + 1}}{μ_{0}}) - \exp (- \frac{τ_{0} - τ_{n}}{μ_{e}}) \exp (- \frac{τ_{n}}{μ_{0}})] . \end{array}

Knowing the contribution of single and multiple scattering, the total intensities of the diffusely reflected light from the TOA and the transmitted light through the BOA are obtained:

I_{r / s, d i f} = I_{r / t, s} + I_{r / t}^{(0)} + 2 \sum_{m = 1}^{\infty} I_{r / t}^{(m)} .

Recall that this expression is for the diffuse (single and multiply scattered) light. In the case of transmission through the BOA, the directly transmitted component in the direction (µ ₀, ϕ ₀), which can be estimated from F ₀, τ₀, and µ₀, has to be included for the total field. Our development in this section is parallel to that of Refs [1]-[2], reproduces their results, and extends their formalism to diffuse transmission through a vertically inhomogeneous atmosphere.

3. Vector radiative transfer

In analogy to the phase function used in scalar radiative transfer, the radiative transfer of polarized light is described by a 4 × 4 Mueller matrix $P (u, u_{0}, ϕ - ϕ_{0}, n)$ which is defined as the multiplication of the phase matrix P(cosΘ,n) with two rotational matrices about angles i ₁ and i ₂ [13], namely

P (u, u_{0}, ϕ - ϕ_{0}, n) = L (π - i_{2}) P (\cos Θ, n) L (- i_{1}) .

To eliminate the azimuthal angle dependence in the radiative transfer calculation, we use the analytical form of Fourier series expansion of the Mueller matrix [14],

\begin{array}{l} P (u, u_{0}, ϕ - ϕ_{0}, n) = \\ P_{c}^{(0)} (u, u_{0}) + 2 \sum_{m = 1}^{\infty} [P_{c}^{(m)} (u, u_{0}, n) \cos m (ϕ - ϕ_{0}) + P_{s}^{(m)} (u, u_{0}, n) \sin m (ϕ - ϕ_{0})], \end{array}

where, for particles with a plane of symmetry [15]

P_{c}^{(m)} (μ, μ_{0}) = [\begin{matrix} P_{11, c}^{(m)} & P_{12, c}^{(m)} & 0 & 0 \\ P_{21, c}^{(m)} & P_{22, c}^{(m)} & 0 & 0 \\ 0 & 0 & P_{33, c}^{(m)} & P_{34, c}^{(m)} \\ 0 & 0 & P_{43, c}^{(m)} & P_{44, c}^{(m)} \end{matrix}]

and

P_{s}^{(m)} (μ, μ_{0}) = [\begin{matrix} 0 & 0 & P_{13, s}^{(m)} & P_{14, s}^{(m)} \\ 0 & 0 & P_{23, s}^{(m)} & P_{24, s}^{(m)} \\ P_{31, s}^{(m)} & P_{32, s}^{(m)} & 0 & 0 \\ P_{41, s}^{(m)} & P_{42, s}^{(m)} & 0 & 0 \end{matrix}] .

Again, by use of μ = |u| four scattering regimes in a similar form as Eqs. (3)–(6) and Eqs. (7)–(10) can be distinguished for the Mueller matrix as well as its series, with the only difference of replacing the scalar phase function “P” by the Mueller matrix “P”. In the same manner of expansion, the series form of the polarized light intensity described by the Stokes vector [I, Q, U, V]^T is:

I (τ, μ, ϕ) = I_{c}^{(m)} (τ, μ) + 2 \sum_{m = 1}^{\infty} [I_{c}^{(m)} (τ, μ) \cos m (ϕ - ϕ_{0}) + I_{s}^{(m)} (τ, μ) \sin m (ϕ - ϕ_{0})],

where the m ^th Fourier series component of the stokes vector

I_{c}^{(m)}

=

{[I_{c}^{(m)}, Q_{c}^{(m)}, U_{c}^{(m)}, V_{c}^{(m)}]}^{T}

and

I_{s}^{(m)}

=

{[I_{s}^{(m)}, Q_{s}^{(m)}, U_{s}^{(m)}, V_{s}^{(m)}]}^{T}

.

3.1 Polarized radiative transfer through multiple scattering

Analysis of polarized radiative transfer can proceed in a similar way as for the scalar case so that the polarized radiative transfer through multiple scattering can again be described as matrix operation,

[\begin{array}{l} 2 I_{r,c}^{(m)} (τ = 0) \\ 2 I_{r,s}^{(m)} (τ = 0) \end{array}] = [\begin{matrix} R_{c}^{(m)} & - R_{s}^{(m)} \\ R_{s}^{(m)} & R_{c}^{(m)} \end{matrix}] [\begin{array}{l} Π_{0, c}^{(m)} \\ Π_{0, s}^{(m)} \end{array}] + [\begin{matrix} R_{c}^{(m)} & - R_{s}^{(m)} \\ R_{s}^{(m)} & R_{c}^{(m)} \end{matrix}] [\begin{matrix} Q_{c}^{(m)} & - Q_{s}^{(m)} \\ Q_{s}^{(m)} & Q_{c}^{(m)} \end{matrix}] [\begin{array}{l} Π_{0, c}^{(m)} \\ Π_{0, s}^{(m)} \end{array}] + ...

for diffuse reflection and

[\begin{array}{l} 2 I_{t, c}^{(m)} (τ = τ_{0}) \\ 2 I_{t, s}^{(m)} (τ = τ_{0}) \end{array}] = [\begin{matrix} T_{c}^{(m)} & - T_{s}^{(m)} \\ T_{s}^{(m)} & T_{c}^{(m)} \end{matrix}] [\begin{array}{l} Π_{0, c}^{(m)} \\ Π_{0, s}^{(m)} \end{array}] + [\begin{matrix} T_{c}^{(m)} & - T_{s}^{(m)} \\ T_{s}^{(m)} & T_{c}^{(m)} \end{matrix}] [\begin{matrix} Q_{c}^{(m)} & - Q_{s}^{(m)} \\ Q_{s}^{(m)} & Q_{c}^{(m)} \end{matrix}] [\begin{array}{l} Π_{0, c}^{(m)} \\ Π_{0, s}^{(m)} \end{array}] + ...

for diffuse transmission. For m = 0 only the cosine mode is needed in calculation so that

I_{r,c}^{(m)} (τ = 0) = R_{c}^{(m)} Π_{0, c}^{(m)} + R_{c}^{(m)} Q_{c}^{(m)} Π_{0, c}^{(m)} + ...

I_{t, c}^{(m)} (τ = τ_{0}) = T_{c}^{(m)} Π_{0, c}^{(m)} + T_{c}^{(m)} Q_{c}^{(m)} Π_{0, c}^{(m)} + ...

The n ^th term on the right hand side of Eqs. (62)-(65) means the contribution of diffuse light after experiencing (n + 1) scattering events. Difference of the vector and scalar radiative transfer case in calculating the Π vector and the Q, R, and T matrices is expressed here. For calculating the contribution of the downwelling incident light to the diffuse field at the TOA and BOA, each element of the R-matrix in Eqs. (62) and (64) becomes a 4×4 cell which is calculated by Eq. (35) after replacing the phase function series

P_{r}^{(m)}

by the Mueller matrix series

P_{r}^{(m)}

. And each element of the T-matrix in Eqs. (63) and (65) becomes a 4×4 cell which is calculated by Eq. (36) after replacing

P_{t}^{(m)}

by

P_{t}^{(m)}

. For calculating the contribution of the upwelling incident light to the diffuse field, the elements of R in Eqs. (62) and (64) are calculated by Eq. (43) after replacing phase function series

P_{t}^{(m), *}

by the Mueller matrix series

P_{t}^{(m), *}

and the elements of T of Eqs. (63) and (65) are calculated by Eq. (44) after replacing

P_{r}^{(m), *}

by

P_{r}^{(m), *}

. For the initial one-time scattered light, the elements of the Π vector becomes

[\begin{array}{l} Π_{(n, i), 0, r, c}^{(m)} \\ Π_{(n, i), 0, r, s}^{(m)} \end{array}] = \frac{w_{i}}{c_{(n, i)}} \exp (- \frac{τ_{n}}{μ_{0}}) \frac{μ_{0}}{μ_{i} + μ_{0}} ω_{0} (n) {1 - \exp [- Δ τ_{n} (\frac{1}{μ_{0}} + \frac{1}{μ_{i}})]} [\begin{array}{l} \frac{P_{r, c}^{(m)} (μ_{0}, μ_{e}, n)}{2} \\ \frac{P_{r, s}^{(m)} (μ_{0}, μ_{e}, n)}{2} \end{array}] F_{0}

and

[\begin{array}{l} Π_{(n, i), 0, t, c}^{(m)} \\ Π_{(n, i), 0, t, s}^{(m)} \end{array}] = \frac{w_{i}}{c_{(n, i)}} \exp (- \frac{τ_{n}}{μ_{0}}) \frac{μ_{0}}{μ_{i} - μ_{0}} ω_{0} (n) {\exp (- \frac{Δ τ_{n}}{μ_{i}}) - \exp (- \frac{Δ τ_{n}}{μ_{0}})} [\begin{array}{l} \frac{P_{t, c}^{(m)} (μ_{0}, μ_{e}, n)}{2} \\ \frac{P_{t, s}^{(m)} (μ_{0}, μ_{e}, n)}{2} \end{array}] F_{0} .

Equations (62)-(65) can be written in the compact form of Eqs. (51) and (52) if we denote the single-scattered light distribution by

Π_{0}^{(m)} = [\begin{array}{l} Π_{0, c}^{(m)} \\ (1 - δ_{0 m}) Π_{0, s}^{(m)} \end{array}],

the transition matrix by

Q_{}^{(m)} = [\begin{matrix} Q_{c}^{(m)} & - (1 - δ_{0 m}) Q_{s}^{(m)} \\ (1 - δ_{0 m}) Q_{s}^{(m)} & (1 - δ_{0 m}) Q_{c}^{(m)} \end{matrix}],

the reflection matrix by

R_{}^{(m)} = [\begin{matrix} R_{c}^{(m)} & - (1 - δ_{0 m}) R_{s}^{(m)} \\ (1 - δ_{0 m}) R_{s}^{(m)} & (1 - δ_{0 m}) R_{c}^{(m)} \end{matrix}],

and the transmission matrix by

T_{}^{(m)} = [\begin{matrix} T_{c}^{(m)} & - (1 - δ_{0 m}) T_{s}^{(m)} \\ (1 - δ_{0 m}) T_{s}^{(m)} & (1 - δ_{0 m}) T_{c}^{(m)} \end{matrix}] .

3.2 “Chain-to-chain adding” strategy

In numerical calculations of the diffuse field at the TOA and BOA, most computation time is spent on solving the matrix Eqs. (51) and (52). Extension from scalar radiative transfer to polarized radiative transfer means an increase of matrix size by a factor of 16 since in scalar case the 1×1 phase function is used while in the vector case a 4×4 Mueller matrix is used. If circular polarization can be neglected then the vector case becomes 3×3. Moreover, for the sake of accuracy a sufficiently small sublayer thickness and a fine Gaussian-Legendre ordinate division on the interval 0 ≤ μ ≤ 1 is required, leading to a large set of linear equations. To save computation time, an “adding” algorithm was proposed in Ref [2], with the aim of reducing the matrix dimension. The strategy is to divide the whole set of sublayers into subsets and apply the Markov chain formalism for each group in a sequential way. Starting from the subset of layers at the bottom, the reflected and transmitted fields are calculated through the Markov chain method described in Sections 2 and 3.1. The subset of layers is then compressed into a single sublayer and added as a new chain knot to the next subset of layers. In an upward manner from the bottom layer subset to the top layer subset, the diffusely reflected and transmitted fields at the required optical depth are obtained. This “chain-to-chain adding” strategy differs from the “layer-to-layer adding” concept in “adding/doubling method” [3]. Because of its ability of handling a number of layers at one time, the “chain-to-chain adding” algorithm is particularly effective in dealing with high vertical inhomogeneity. However, there is a cost for gaining computation efficiency: the advantages of the regular Markov chain method in separating the contribution of a specific orders of scattering as well as getting the internal radiation field are lost.

4. Illustrative numerical calculations

The Markov chain vector 1D radiative transfer model described in the previous sections was programmed in MatLab® for ease of matrix and vector manipulations. Results from this prototype were compared with results from 1D vector Monte Carlo models developed simultaneously in FORTRAN, and a selection of such comparisons are presented in this section.

4.1 Implementation of the Monte Carlo method

We use two Monte Carlo models that are described in detail by Davis and Xu [16]. The only difference between them is in the representation of vertical variability. One code takes exactly the same input information as the Markov chain model, i.e., N layers with their respective single scattering albedo and phase function; it then tracks histories in optical units, assuming the layers are internally uniform. This version of the Monte Carlo model was benchmarked extensively against the high-precision values provided by Kokhanovsky et al. [17] for a Rayleigh case, an aerosol case, and a cloud case, all single uniform layers over an absorbing surface.

The other code tracks trajectories, where convenient, in physical units and treats spatial variability continuously but in a specific parameterized representation; for the present study, the atmosphere is composed of an exponential distribution of (Rayleigh scattering) molecules and a uniform layer of particles, e.g., polydisperse (Mie scattering) spheres.

In all of the following, 10⁸ Monte Carlo histories were traced using the code that ingests discrete layers; this leads to relative uncertainties on the order of 10⁻⁴ for well-sampled signals but this precision deteriorates for poorly populated signals (typically, the 3^rd and 4^th components of the Stokes vector). Radiances were computed using the local estimation technique [18, 19] at the same ordinates as used by the Markov chain model in the principal and perpendicular planes.

4.2 Separate orders of scattering in a pure Rayleigh atmosphere

We start our comparisons of the Markov chain and Monte Carlo methods in predicting various orders of scattering by a homogeneous atmosphere of pure Rayleigh scattering with optical depth τ₀ = 0.5 and single scattering albedo ω₀ = 1. The atmosphere is divided into 20 layers. In Markov chain formalism, contributions of single and higher-order scattering are calculated by Eqs. (54)-(55) after replacing the phase function series $P_{r/t}^{(m)}$ by the Mueller matrix series $P_{r/t}^{(m)}$ , and by substituting relevant terms in Eqs. (62)-(65) into Eq. (61), respectively. In the Monte Carlo method, tallying the various orders of scattering is trivial to encode. An illustrative calculation of the first two, three, then all, orders of scattering contributing to the diffusely reflected and transmitted fields in the principal plane (ϕ-ϕ ₀ = 0° and 180°) is given in Figs. 1 -4 , showing agreement to be better than 99%. In these figures, only the I and Q components are plotted since the U and V components vanish in the principal plane. Moreover, positive viewing zenith angles are used for the ϕ-ϕ ₀ = 0° plane while their negative counterparts are used for ϕ-ϕ ₀ = 180°. Note that, although the Markov chain formalism is presented for the vertically inhomogeneous atmosphere in the current paper, by setting the Mueller matrix P(cosΘ,n) and single scattering albedo ω₀(n) to be constant for all sublayers, the formalism reduces to its counterpart for the homogeneous atmosphere [1].

Fig. 1 The first two, three, then all, orders of scattering contributing to the diffuse reflection field (I component) in the principal plane for a homogeneous Rayleigh atmosphere of optical depth τ₀ = 0.5 and unit single scattering albedo. Error bars are based on the standard deviation among 10 realizations of the Monte Carlo simulation, each with 10⁷ histories (cf. Figures 12 and 16 for a situation with degraded precision).

Download Full Size | PDF

Fig. 12 Same as Fig. 9 but the V component of the diffusely reflected field is plotted.

Download Full Size | PDF

Fig. 4 Same as Fig. 1 but the Q component of the diffuse transmission field is plotted.

Download Full Size | PDF

Fig. 2 Same as Fig. 1 but the Q component of the diffuse reflection field is plotted.

Download Full Size | PDF

Fig. 3 Same as Fig. 1 but the I component of the diffuse transmission field is plotted.

Download Full Size | PDF

4.3 Vertically stratified aerosol-and-Rayleigh atmosphere

As another example of numerical calculation, we choose an Earth-like atmosphere of total optical depth τ₀ = 1.22934. The atmosphere is composed of Mie aerosols (optical depth 1.0074) of refractive index n = 1.38+10⁻⁸i and Rayleigh molecules (optical depth 0.2286). For the aerosols, a lognormal distribution with parameters s = 1 and a ₀ = 0.1 μm is assumed:

f (a) = \frac{1}{\sqrt{2 π} (s a)} \exp {- \frac{1}{2} {[\frac{\ln (a / a_{0})}{s}]}^{2}} .

This aerosol model was used in a recent intercomparison of aerosol retrieval algorithms described in Ref [20]. and the aerosol optical depth corresponds to their “Case 12.” The phase matrix for the aerosol distribution is calculated by integrating the contribution of particles over the size interval 0 ≤ a ≤ 30 μm, which is discretized to 40000 quadrature points. To highlight molecular scattering, we selected the MISR blue wavelength of 446.4 nm [21], leading to the stated Rayleigh optical depth.

The atmosphere is divided into 42 sublayers, with the optical thickness of each layer as well as the Rayleigh scattering fraction in each layer are listed in Table 1 . This layering approximates in optical distance space the combination of an exponential decay of the Rayleigh column with a scale height of 8 km with a uniform aerosol layer between 0 and 2 km. The surface is black and solar illumination is at zenith angle θ ₀ = 60°. This atmospheric structure, as well as the lower-boundary and illumination conditions, were prescribed in Ref [20]. for 443 nm albeit at a slightly different wavelength (hence Rayleigh optical depth and aerosol phase matrix) and without molecular depolarization effects. The layers in Table 1 were generated automatically from top to bottom under the sole condition that none of them has an optical thickness exceeding 0.03, which is taken in Markov chain method as a reasonable rule-of-thumb to ensure the single scattering domination in each layer.

Table 1. Optical thickness and Rayleigh fraction of the atmosphere divided into 42 sublayers

View Table

Figures 5 -6 give the I and Q components of the reflection field in the principal plane (ϕ-ϕ ₀ = 0° and 180°) while Figs. 7 -8 give the I and Q components of the transmission field in the same plane. The U and V components vanish identically in the principal plane. In the perpendicular plane ϕ-ϕ ₀ = ±90°, however, all components of the reflected and transmitted fields appear, as demonstrated in the panels of Figs. 9 -16 . Using 45 Gauss-Legendre ordinates on the interval 0 < μ < 1 (equivalently, 90 streams), Figs. 5-8 and Figs. 9-16 show the disagreement between the Markov chain and Monte Carlo methods to be less than 1%. Using finer (coarse) sublayer thickness and more (less) Gauss-Legendre ordinates on μ’s interval can further reduce (increase) the disagreement, but costing (saving) in computation time. For example, the maximum error of Q for diffuse reflection and transmission in the principal and perpendicular planes increases to about 4% when 23 Gauss-Legendre ordinates are used but 75% computation time is saved.

Fig. 5 Diffuse reflection field (I component) in the principal plane of the vertically inhomogeneous atmosphere described in the main text under solar illumination at 60° incidence. As in Fig. 1, the error bars are based on the standard deviation among 10 realizations of the Monte Carlo simulation, each with 10⁷ histories.

Download Full Size | PDF

Fig. 6 Same as Fig. 5 but the Q component of the diffuse reflection field is plotted.

Download Full Size | PDF

Fig. 7 Same as Fig. 5 but the I component of the diffuse transmission field is plotted.

Download Full Size | PDF

Fig. 8 Same as Fig. 5 but the Q component of the diffuse transmission field is plotted.

Download Full Size | PDF

Fig. 9 Same as Fig. 5 but the I component of the diffusely reflected field in the perpendicular plane (ϕ−ϕ ₀ = ±90°) is plotted.

Download Full Size | PDF

Fig. 16 Same as Fig. 9 but the V component of the diffusely transmitted field is plotted.

Download Full Size | PDF

Fig. 10 Same as Fig. 9 but the Q component of the diffusely reflected field is plotted.

Download Full Size | PDF

Fig. 11 Same as Fig. 9 but the U component of the diffusely reflected field is plotted.

Download Full Size | PDF

Fig. 13 Same as Fig. 9 but the I component of the diffusely transmitted field is plotted.

Download Full Size | PDF

Fig. 14 Same as Fig. 9 but the Q component of the diffusely transmitted field is plotted.

Download Full Size | PDF

Fig. 15 Same as Fig. 9 but the U component of the diffusely transmitted field is plotted.

Download Full Size | PDF

5. Conclusion and outlook

Accurate yet efficient modeling of vector (polarized) polarized radiative transfer is becoming increasingly important, largely to keep up with instrument development in several application areas: medical imaging, ocean and atmospheric optics (including remote sensing), planetary missions, etc. Space-based Earth observation sensors with polarization capability include various incarnations of the POLDER (POLarization and Directionality of Earth Reflectances) sensor, most recently on ESA’s PARASOL satellite [22], and the APS (Aerosol Polarimetry Sensor) on NASA’s upcoming Glory mission [23].

Even when polarization is not a focus of the mission, it needs to be accounted for to keep the accuracy of the forward radiative transfer modeling at par with instrument precision. For instance the Orbiting Carbon Observatory (OCO) mission [24] targets global column amounts of CO₂ based on differential optical absorption spectroscopy (DOAS) in the near-IR. This calls for DOAS in both CO₂ and O₂ bands, and the later is significantly affected by aerosol scattering, therefore by polarization [25]. To remedy this situation without sacrificing computational efficiency, Natraj and Spurr [26] developed an approximation model based on two orders of scattering that compensates for polarization effects left out by OCO’s operational scalar code.

Since the Markov chain formalism has inherent computational efficiency in handling all orders of scattering from a vertically inhomogeneous atmosphere, it forms another attractive solution to polarized radiative transfer. Comparison of the diffusely reflected and transmitted radiation fields calculated by Markov chain method to those by Monte-Carlo method shows the deviation to be less than 1%, depending on the optical depth and quadrature mesh. On the basis of the current work, we plan to further extend the Markov chain approach to the cases concerning non-uniform and/or bi-directional polarized surface reflectivity, spherical shell geometry, and more structurally complex 3D atmospheres.

Acknowledgements

The research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. F. Xu is supported by an appointment to the NASA Postdoctoral Program (NPP) at the Jet Propulsion Laboratory (JPL); the NPP is administered by Oak Ridge Associated Universities under contract with NASA. We thank Michael Garay for providing some successive orders of scattering results for benchmarking in the early stage of this project. The whole “(Vector) Radiative Transfer” task force at JPL is also acknowledged for lively discussions on a weekly basis. Copyright 2010. All rights reserved.

References and links

1. L. W. Esposito and L. L. House, “Radiative transfer calculated by a Markov chain formalism,” Astrophys. J. 219, 1058–1067 (1978). [CrossRef]

2. L. W. Esposito, “An ‘adding’ algorithm for the Markov chain formalism for radiation transfer,” Astrophys. J. 233, 661–663 (1979). [CrossRef]

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

4. K. F. Evans and G. L. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant. Spectrosc. Radiat. Transf. 46(5), 413–423 (1991). [CrossRef]

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

6. 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(26), 6762–6774 (2006). [CrossRef] [PubMed]

7. P. W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17(4), 2057–2079 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-4-2057. [CrossRef] [PubMed]

8. J. H. Hannay, “Radiative transfer: exact Rayleigh scattering series and a daylight formula,” J. Opt. Soc. Am. A 26(3), 669–675 (2009). [CrossRef]

9. 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. Transf. 47(1), 19–33 (1992). [CrossRef]

10. R. J. D. Spurr, “VLIDORT: A linearized pseudo-spherical vector discrete ordinate radiative transfer code for forward model and retrieval studies in multilayer multiple scattering media,” J. Quant. Spectrosc. Radiat. Transf. 102(2), 316–342 (2006). [CrossRef]

11. 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. Transf. 61(1), 105–122 (1999). [CrossRef]

12. J. V. Dave, “Intensity and polarization of the radiation emerging from a plane-parallel atmosphere containing monodispersed aerosols,” Appl. Opt. 9(12), 2673–2684 (1970). [CrossRef] [PubMed]

13. J. W. Hovenier, “Symmetry relations for scattering of polarized light in a slab of randomly oriented particles,” J. Atmos. Sci. 26(3), 488–499 (1969). [CrossRef]

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

15. J. W. Hovenier, and C. V. M. van der Mee, “Basic Relationships for Matrices Describing Scattering by Small Particles”, in Light Scattering by Nonspherical Particles, pp. 61–85, M. Mishchenko, J. W. Hovenier and L. Travis, (eds.), Academic Press, San Diego, 2000.

16. A. B. Davis and F. Xu, “Monte Carlo modeling of polarized light transfer in vertically varying plane-parallel atmospheres, with application to lofted aerosol layer detection using O₂ spectroscopy,” J. Quant. Spectrosc. Radiat. Transf.

17. A. A. Kokhanovsky, V. P. Budak, C. Cornet, M. Duan, C. Emde, I. L. Katsev, D. A. Klyukov, S. V. Korkin, L. C-Labonnote, and B. Mayer, “Benchmark results in vector atmospheric radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 111(12-13), 1931–1946 (2010). [CrossRef]

18. G. Marchuk, G. Mikhailov, N. Nazaraliev, R. Darbinjan, B. Kargin, and B. Elepov, The Monte Carlo Methods in Atmospheric Optics, Springer-Verlag, New-York, NY (1980).

19. K. F. Evans, and A. Marshak, “Numerical Methods,” in 3D Radiative Transfer in Cloudy Atmospheres, pp. 243–281, A. Marshak and A. B. Davis (eds.), Springer, Heidelberg, Germany (2005).

20. A. A. Kokhanovsky, J. L. Deuzé, D. J. Diner, O. Dubovik, F. Ducos, C. Emde, M. J. Garay, R. G. Grainger, A. Heckel, M. Herman, I. L. Katsev, J. Keller, R. Levy, P. R. J. North, A. S. Prikhach, V. V. Rozanov, A. M. Sayer, Y. Ota, D. Tanré, G. E. Thomas, E. P. Zege, A. A. Kokhanovsky, J. L. Deuzé, D. J. Diner, O. Dubovik, F. Ducos, C. Emde, M. J. Garay, R. G. Grainger, A. Heckel, M. Herman, I. L. Katsev, J. Keller, R. Levy, P. R. J. North, A. S. Prikhach, V. V. Rozanov, A. M. Sayer, Y. Ota, D. Tanré, G. E. Thomas, and E. P. Zege, “The inter-comparison of major satellite aerosol retrieval algorithms using simulated intensity and polarization characteristics of reflected light,” Atm. Meas. Techn. 3, 909–932 (2010). [CrossRef]

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

22. P. Y. Deschamps, F. M. Breon, M. Leroy, A. Podaire, A. Bricaud, J. C. Buriez, and G. Seze, “The POLDER mission: Instrument characteristics and scientific objectives,” IEEE Trans. Geosci. Rem. Sens. 32(3), 598–615 (1994). [CrossRef]

23. 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. Am. Meteorol. Soc. 88(5), 677–691 (2007). [CrossRef]

24. D. Crisp, R. Atlas, F. Breon, L. Brown, J. Burrows, P. Ciais, B. Connor, S. Doney, I. Fung, and D. Jacob, “The Orbiting Carbon Observatory (OCO) mission,” Adv. Space Res. 34(4), 700–709 (2004). [CrossRef]

25. V. Natraj, R. J. D. Spurr, H. Boesch, Y. Jiang, and Y. Yung, “Evaluation of errors in neglecting polarization in the forward modeling of O2 A band measurements from space, with relevance to CO2 column retrieval from polarization sensitive instruments,” J. Quant. Spectrosc. Radiat. Transf. 103(2), 245–259 (2007). [CrossRef]

26. V. Natraj and R. J. D. Spurr, “A fast linearized pseudo-spherical two orders of scattering model to account for polarization in vertically inhomogeneous scattering-absorbing media,” J. Quant. Spectrosc. Radiat. Transf. 107(2), 263–293 (2007). [CrossRef]

Layer Number	Optical Thickness	Rayleigh Fraction
1	3.0000E-02	1.0000E+00
2	3.0000E-02	1.0000E+00
3	3.0000E-02	1.0000E+00
4	3.0000E-02	1.0000E+00
5	3.0000E-02	1.0000E+00
6	2.8034E-02	1.0000E+00
7	2.9646E-02	4.2727E-02
8	2.9655E-02	4.3017E-02
9	2.9664E-02	4.3310E-02
10	2.9673E-02	4.3605E-02
11	2.9682E-02	4.3902E-02
12	2.9692E-02	4.4200E-02
13	2.9701E-02	4.4501E-02
14	2.9710E-02	4.4803E-02
15	2.9720E-02	4.5107E-02
16	2.9729E-02	4.5414E-02
17	2.9739E-02	4.5722E-02
18	2.9749E-02	4.6032E-02
19	2.9758E-02	4.6345E-02
20	2.9768E-02	4.6659E-02
21	2.9778E-02	4.6976E-02
22	2.9788E-02	4.7294E-02
23	2.9798E-02	4.7614E-02
24	2.9808E-02	4.7937E-02
25	2.9818E-02	4.8262E-02
26	2.9829E-02	4.8588E-02
27	2.9839E-02	4.8917E-02
28	2.9849E-02	4.9248E-02
29	2.9860E-02	4.9581E-02
30	2.9870E-02	4.9916E-02
31	2.9881E-02	5.0253E-02
32	2.9892E-02	5.0593E-02
33	2.9902E-02	5.0934E-02
34	2.9913E-02	5.1278E-02
35	2.9924E-02	5.1624E-02
36	2.9935E-02	5.1972E-02
37	2.9946E-02	5.2323E-02
38	2.9957E-02	5.2675E-02
39	2.9969E-02	5.3030E-02
40	2.9980E-02	5.3388E-02
41	2.9991E-02	5.3747E-02
42	7.8896E-03	5.3975E-02

Markov chain formalism for polarized light transfer in plane-parallel atmospheres, with numerical comparison to the Monte Carlo method

Abstract

1. Introduction and overview

2. Scalar radiative transfer

2.1 Initial scattering in the atmosphere

2.2 Intermediate scattering in the atmosphere

2.3 Scattering to the boundary

2.3.1 Downwelling light

2.3.2 Upwelling light

2.4 Radiative transfer through multiple scattering

(2.5) Total radiation field

3. Vector radiative transfer

3.1 Polarized radiative transfer through multiple scattering

3.2 “Chain-to-chain adding” strategy

4. Illustrative numerical calculations

4.1 Implementation of the Monte Carlo method

4.2 Separate orders of scattering in a pure Rayleigh atmosphere

4.3 Vertically stratified aerosol-and-Rayleigh atmosphere

5. Conclusion and outlook

Acknowledgements

References and links

Cited By

Figures (16)

Tables (1)

Equations (72)

Optics Express