Polarized Radiative Transfer Simulations: A Tutorial Review and Upgrades of the Vector Discrete Ordinate Radiative Transfer Computational Tool

We present an overview and several important upgrades to the Vector Discrete Ordinate Radiative Transfer (VDISORT) code. VDISORT is a polarized (vector) radiative transfer code that can be applied to a wide range of research problems including the Earth’s atmosphere and ocean system. First, a solution is developed to the complex algebraic eigenvalue problem resulting when the b₂ component of the Stokes scattering matrix is non-zero. This solution is needed to compute the V component of the Stokes vector $\mathbf{I}={[I_{\parallel },I_{\perp },U,V]}^T$. Second, a significant improvement in computational efficiency is obtained by reducing the dimension of the algebraic eigenvalue by a factor of 2 resulting in a speed increase of about 2³ = 8. Third, an important upgrade of the VDISORT code is obtained by developing and implementing a method to enable output at arbitrary polar angles by the integration of the source function (ISF) method for partially reflecting Lambertian as well as general non-Lambertian surfaces. Fourth, a pseudo-spherical treatment has been implemented to provide important corrections for Earth curvature effects at near horizontal solar zenith and observation (viewing) polar angles. Fifth, a post-processing single-scattering correction procedure has been developed to enhance the accuracy and speed for strongly forward-peaked scattering. With these significant improvements the results from the upgraded version of the VDISORT code match published benchmark results for Rayleigh scattering, Mie scattering, and scattering by non-spherical cirrus particles. The performance of VDISORT for a polarized incident beam source is equally satisfactory. The VDISORT vector radiative transfer code is made public and freely available for use by the growing polarimetric research community including the space-borne polarimeters on the future NASA PACE and AOS missions.

1 Introduction

This paper documents a state-of-the-art numerical code called VDISORT for monochromatic polarized radiative transfer in non-isothermal, vertically inhomogeneous, but horizontally homogeneous media. The physical processes included are Planckian thermal emission, scattering with a general phase matrix, absorption, and surface reflection. The system may be driven by parallel or isotropic unpolarized or polarized radiation incident at the boundaries, as well as by internal thermal sources and thermal emission from the boundaries. The Stokes vector is returned at user-specified angles and levels. The azimuthally-averaged Stokes components are also available. Irradiances and mean radiances can be generated from the azimuthally-averaged radiance, i.e., the first component of the Stokes vector [I,Q,U,V]^T.

VDISORT is based on articles published in the open literature, but its theoretical background and algorithmic developments have not been systematically described in one single document. This paper provides an up-to-date complete description of VDISORT including a self-contained account of its theoretical basis and a discussion of the numerical implementation of the theory.

Like DISORT, the scalar version, VDISORT has been designed to be a good scientific software package and a numerical code of general utility. The VDISORT package takes advantage of robust existing software tools to make it numerically well-conditioned, and user-friendly. A set of test cases have been adopted to verify the numerical code against published results.

1.1 Outline of Paper

This paper is organized as follows. Section 1 provides the motivation for and a brief history of VDISORT, while Section 2 provides theoretical background, introduces the radiative transfer equation for unpolarized and polarized radiation, and discusses the inherent optical properties including the phase matrix. The solution of the vector radiative transfer equation (VRTE) is discussed in Section 3 including the discrete ordinate method, important upgrades of the vector discrete ordinate code (VDISORT), the ISF method, and treatment of polarized reflectance from the lower boundary. Section 4 discusses the merits of the 4 × 4 solution versus the 3 × 3 approximation, while Section 5 is devoted to the single-scattering approximation, which is used in a post-processing step to enhance the accuracy of the computed results for a given number of discrete ordinate streams. Test results are presented in Section 6 including comparisons with published benchmarks. Finally, in Section 7 a brief summary and concluding remarks are provided.

1.2 Brief History

The discrete ordinate radiative transfer algorithm (DISORT) has proven to be an accurate, versatile and reliable method for the solution of the scalar radiative transfer problem in plane-parallel, vertically inhomogeneous media (Stamnes et al., 1988; Lin et al., 2015; Laszlo et al., 2016; Stamnes K. et al., 2017). An extension of the scalar discrete ordinate theory to solve for the complete Stokes vector I = [I,Q,U,V]^T was reported by Weng (1992), who adopted an approach to the solution of the vector problem completely analogous to the scalar case. Hence, the computer code for the vector problem, VDISORT, could rely on the same well-tested routine to obtain the eigenvalues and eigenvectors as the one used in the scalar version (DISORT). Also, the same scaling transformation (Stamnes and Conklin, 1984) could be applied to circumvent the notorious ill-conditioning that occurs when applying boundary and layer interface continuity conditions.

The first version of the FORTRAN code developed by Weng (1992) had a few shortcomings related to the fact that it had been applied exclusively in the microwave region, and thus had not been tested for beam source applications. In addition, the procedure employed to compute the Fourier components of the phase matrix turned out to be both inaccurate and inefficient. To correct these shortcomings an improved version of the code was developed by Schulz et al. (1999). In this new version of the code 1) errors in the numerical implementation were corrected, 2) the procedure used to compute the Fourier components of the phase matrix was replaced by a more accurate and efficient method, 3) the basic performance of the code was tested against benchmark results. However, although the code seemed to have the potential to become an accurate and reliable tool for a variety of applications, no attempt was made to test it in a systematic and comprehensive manner. Also, no attempt was made to document the code thoroughly and extensively, and it was assumed that the homogeneous solution involved only real eigenvalues/eigenvectors, which are sufficient to solve for the I, Q, and U Stokes parameters, but the V component requires complex arithmetic.

The original code provided solutions for the I, Q, and U Stokes parameters at the discrete ordinates (i.e. at the quadrature polar angles). Since the computing time required for the discrete ordinate method increases cubically with the number of quadrature angles, it becomes cost-effective to obtain the solution at a limited number of quadrature angles and then generate the solution at additional angles by using an efficient interpolation scheme. To this end analytic expressions were developed that were shown to be accurate for the I, Q, and U components, but not for the V component, at arbitrary angles and optical depths (Schulz and Stamnes, 2000) in much the same way they were developed for the scalar version (Stamnes, 1982). These analytic expressions obtained by the ISF method (Stamnes, 1982) satisfy not only the radiative transfer equation, but also the boundary and layer-interface continuity conditions at arbitrary polar angles, and they have proven to be superior to standard interpolation schemes such as cubic splines (Stamnes, 1982; Schulz and Stamnes, 2000). To complete this development, an extension of this ISF method to include accurate computation of the V Stokes parameter, which requires complex arithmetic, is described Section 3.5.

2 Theoretical Basis

2.1 Unpolarized Radiation

We consider an inhomogeneous horizontal slab of scattering/absorbing material with inherent optical properties that vary only in the vertical direction z, where z increases upward. The corresponding vertical optical depth is defined by

$\tau \left(z\right)={\int }_z^{\infty }\left[\alpha \left(z^{\prime }\right)+\beta \left(z^{\prime }\right)\right]d{z}^{\prime }(1)$

where α and β are the absorption and scattering coefficients in units of reciprocal length, respectively, and the vertical optical depth is defined to increase downward from τ(z = ∞) = 0 at the top of the slab. The slab is assumed to be in local thermodynamic equilibrium so that it emits radiation according to the local temperature T (τ(z)). The diffuse radiance distribution I (τ, u, ϕ) can be described by the scalar radiative transfer equation (RTE)

$u\frac{dI\left(\tau ,u,\varphi \right)}{d\tau }=I\left(\tau ,u,\varphi \right)-S\left(\tau ,u,\varphi \right)(2)$

where the source function is given by

$S\left(\tau ,u,\varphi \right)=S^{\ast }\left(\tau ,u,\varphi \right)+[1-\stackrel{\text{—}}{\varpi }\left(+\right)]B\left(\tau \right)+\frac{\stackrel{\text{—}}{\varpi }\left(\tau \right)}{4\pi }{\int }_0^{2\pi }d{\varphi }^{\prime }{\int }_{-1}^{1}p\left(\tau ,u^{\prime },{\varphi }^{\prime };u,\varphi \right)I\left(\tau ,u^{\prime },{\varphi }^{\prime }\right)d{u}^{\prime }.(3)$

Here u is the cosine of the polar angle θ, ϕ is the azimuth angle, ϖ(τ) = β(τ)/[α(τ) + β(τ)] is the single-scattering albedo, p (τ, u′, ϕ′; u, ϕ) is the scattering phase function, and B(τ) is the thermal radiation field given by the Planck function. The single-scattering source term is given by

$S^{\ast }\left(\tau ,u,\varphi \right)=\frac{\varpi \left(\tau \right)}{4\pi }p\left(\tau ,-{\mu }_0,{\varphi }_0;u,\varphi \right)S_b{e}^{-\tau /{\mu }_0}(4)$

where S_b is the incident (solar) irradiance and μ₀ is the cosine of the solar zenith angle. The differential vertical optical depth is [see Eq. 1]

$d\tau \left(z\right)=-\left[\alpha \left(\tau \right)+\beta \left(\tau \right)\right]dz(5)$

where the minus sign indicates that τ increases in the downward direction, whereas z increases in the upward direction, as noted above. The scattering angle Θ and the polar and azimuth angles are related by (see Figure 1)

${\widehat{\mathbf{\Omega }}}^{\prime }\cdot \widehat{\mathbf{\Omega }}=\cos \mathrm{\Theta }=\cos \theta \cos {\theta }^{\prime }+\sin {\theta }^{\prime }\sin \theta \cos \left({\varphi }^{\prime }-\varphi \right).$

FIGURE 1

FIGURE 1. Coordinate system for scattering by a volume element at O. The points C, A and B are located on the unit sphere. The incident light beam with Stokes vector ${\mathbf{I}}_{\text{S}}^{\text{inc}}$ is in direction OA (θ′, ϕ′) with unit vector ${\widehat{\mathbf{\Omega }}}^{\prime }$, the scattered beam with Stokes vector ${\mathbf{I}}_{\text{S}}^{\text{sca}}$ is in direction OB(θ, ϕ) with unit vector $\widehat{\mathbf{\Omega }}$ (Hovenier et al., 2004).

Here ${\widehat{\mathrm{\Omega }}}^{\prime }$ is the unit vector of the incident beam direction and $\widehat{\mathrm{\Omega }}$ is the unit vector of the scattered direction. By definition, θ = 180° is directed toward nadir (straight down) and θ = 0° toward zenith (straight up). Thus, u = cos θ varies in the range [ − 1, 1] (from nadir to zenith). For oblique illumination of the slab, ϕ₀ = 180° is defined to be the azimuth angle of the incident light.

2.2 Polarized Radiation

To generalize Eq. 2 to apply to polarized radiation, the scalar source function must be replaced by the appropriate vector version. Hence, the multiple-scattering term $S^{\text{ms}}(\tau ,u,\varphi )=\frac{\varpi (\tau )}{4\pi }{\int }_0^{2\pi }d{\varphi }^{\prime }{\int }_{-1}^{1}d{u}^{\prime }p(\tau ,u^{\prime },{\varphi }^{\prime };u,\varphi )I(\tau ,u^{\prime },{\varphi }^{\prime })$ in Eq. 3 must be replaced by

${\mathbf{S}}^{\text{ms}}\left(\tau ,u,\varphi \right)=\frac{\varpi \left(\tau \right)}{4\pi }{\underset{0}{\int }}^{2\pi }d{\varphi }^{\prime }{\underset{-1}{\int }}^{1}d{u}^{\prime }\mathbf{P}\left(\tau ,u^{\prime },{\varphi }^{\prime };u,\varphi \right)\mathbf{I}\left(\tau ,u^{\prime },{\varphi }^{\prime }\right)(6)$

where I (τ, u′, ϕ′) is the Stokes vector, and P (τ, u′, ϕ′; u, ϕ) is the scattering phase matrix. The first element of the vector S^ms(τ, u, ϕ) represents the energy per unit solid angle, per unit frequency interval, and per unit time that is scattered by a unit volume in the direction (u = cos θ, ϕ). Hence, in a plane-parallel (slab) geometry, the integro-differential vector radiative transfer equation (VRTE) for polarized radiation is expressed in terms of a Stokes vector I (τ, u, ϕ) as

$u\frac{d\mathbf{I}\left(\tau ,u,\varphi \right)}{d\tau }=\mathbf{I}\left(\tau ,u,\varphi \right)-\mathbf{S}\left(\tau ,u,\varphi \right)(7)$

where the vector source function is

$\mathbf{S}\left(\tau ,u,\varphi \right)={\mathbf{S}}^{\text{ms}}\left(\tau ,u,\varphi \right)+\mathbf{Q}\left(\tau ,u,\varphi \right).(8)$

Here S^ms(τ, u, ϕ) is given by Eq. 6 and the source term Q (τ, u, ϕ), due to beam and thermal sources, is given by:

$\mathbf{Q}\left(\tau ,u,\varphi \right)=\frac{\varpi \left(\tau \right)}{4\pi }\mathbf{P}\left(\tau ,-{\mu }_0,{\varphi }_0;u,\varphi \right){\mathbf{S}}_b{e}^{-\tau /{\mu }_0}+\left[1-\varpi \left(\tau \right)\right]{\mathbf{S}}_t\left(\tau \right).(9)$

The first term on the right hand side of Eq. 9 describes the incident beam S_b in direction (−μ₀, ϕ₀), which is attenuated at depth τ by a factor $e^{-\tau /{\mu }_0}$ and undergoes single scattering into the direction (u, ϕ). For an unpolarized incident beam S_b has the form

${\mathbf{S}}_b={\left[I_0/2,I_0/\mathrm{2,0,0}\right]}^T\mathrm{o}\mathrm{r}{\left[I_0,\mathrm{0,0,0}\right]}^T(10)$

where the superscript ^T denotes the transpose, and where the first or second expression corresponds to the choice of Stokes vector representation, ${[I_{\parallel },I_{\perp },U,V]}^T$ or [I,Q,U,V]^T. The second term on the right hand side of Eq. 9 is due to thermal emission, which is unpolarized, and S_t(τ) is given by

${\mathbf{S}}_t\left(\tau \right)={\left[B\left(T\left(\tau \right)\right)/2,B\left(T\left(\tau \right)\right)/\mathrm{2,0,0}\right]}^T\mathrm{o}\mathrm{r}{\left[B\left(T\left(\tau \right)\right),\mathrm{0,0,0}\right]}^T(11)$

where B is the Planck function, and where the first or second expression corresponds to the choice of Stokes vector representation. We have set μ₀ ≡|u₀|≡| cos θ₀|, where θ₀ is the polar angle of the incident light beam.

2.3 Effects of Earth Curvature–The Pseudo-Spherical Approximation

For many applications plane-parallel geometry is adequate. For large solar zenith angles (θ₀ ≥ 70° and for near horizontal polar viewing angles θ), however, the plane-parallel approximation (PPA) provides inaccurate results. Then the Earth curvature must be considered. Large solar zenith angles occur around the times of sunrise and sunset at any location on a planet. Such large solar zenith angles are present, for example, in observations made by instruments onboard geostationary satellites that observe a large part of Earth’s disk throughout the day. Sensors deployed on polar-orbiting satellite platforms also observe at large solar zenith angles at high latitudes.

As discussed by several investigators (see He et al. (2018) and references therein), the so-called pseudo-spherical approximation (PSA) (Dahlback and Stamnes, 1991) represents a very useful correction to the plane-parallel approximation. In the PSA the direct beam single-scattering term, also called the solar pseudo-source term, is treated in spherical geometry while the multiple-scattering term is treated using the PPA. Hence, in the PSA the exponential attenuation in Eqs. 4 and 9 is replaced by the Chapman function, that is, exp (−τ/μ₀) → exp (−τCh(μ₀)), where the Chapman function Ch(μ₀) takes Earth curvature into account, but ignores refraction. As shown by He et al. (2018) the influence of Earth curvature increases rapidly with solar zenith angle, being up to 1, 3, and 12% for solar zenith angles of 75°, 80°, and 85°, respectively.

2.4 Scattering Phase Matrix

The development of vector radiative transfer theory may start with the Stokes vector representation $\mathbf{I}={[I_{\parallel },I_{\perp },U,V]}^T$. In terms of the complex transverse electric field components of the radiation field $E_{\parallel }=|E_{\parallel }|{\mathrm{e}}^{-i{\delta }_{\parallel }}$ and $E_{\perp }=|E_{\perp }|{\mathrm{e}}^{-i{\delta }_{\perp }}$, these Stokes vector components are given by:

$\begin{array}{ccc}\hfill I_{\parallel }& =\hfill & E_{\parallel }{E}_{\parallel }^{\ast }\hfill \\ \hfill I_{\perp }& =\hfill & E_{\perp }{E}_{\perp }^{\ast }\hfill \\ \hfill U& =\hfill & 2|E_{\parallel }||E_{\perp }|\cos \delta \hfill \\ \hfill V& =\hfill & 2|E_{\parallel }||E_{\perp }|\sin \delta \hfill \end{array}(12)$

where the phase difference δ is δ_∥ − δ_⊥. The connection between the $\mathbf{I}={[I_{\parallel },I_{\perp },U,V]}^T$ Stokes vector representation and the more commonly used I_S = [I,Q,U,V]^T representation, where I = I_∥ + I_⊥ and Q = I_∥ − I_⊥, is given by:

${\mathbf{I}}_S=\mathbf{\text{DI}}(13)$

where

$\mathbf{D}=\left(\begin{array}{cccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right),{\mathbf{D}}^{-1}=\frac{1}{2}\left(\begin{array}{cccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 2\hfill \end{array}\right).(14)$

The degree of polarization is defined as

$p={\left[Q^2+U^2+V^2\right]}^{1/2}/I(15)$

so that 0 ≤ p ≤ 1, where p = 1 corresponds to completely polarized light and p = 0 to natural (unpolarized) light. The degree of circular polarization is defined as

$p_c=V/I,(16)$

the degree of linear polarization as

$p_l={\left[Q^2+U^2\right]}^{1/2}/I,(17)$

and alternatively, when U = 0 as

$p_l\left(U=0\right)=-\frac{Q}{I}=\frac{I_{\perp }-I_{\parallel }}{I_{\perp }+I_{\parallel }}.(18)$

The transverse electric field vector ${[E_{\parallel },E_{\perp }]}^T$ of the scattered field can be obtained in terms of the transverse field vector ${[E_{\parallel 0},E_{\text{⊥}0}]}^T$ of the incident field by a linear transformation:

$\left(\begin{array}{c}\hfill E_{\parallel }\hfill \\ \hfill E_{\perp }\hfill \end{array}\right)=\mathbf{A}\left(\begin{array}{c}\hfill E_{\parallel 0}\hfill \\ \hfill E_{\perp 0}\hfill \end{array}\right)(19)$

where A is a 2 × 2 matrix, referred to as the amplitude scattering matrix. The corresponding linear transformation connecting the Stokes vectors of the incident and scattered fields in the scattering plane is called the Mueller matrix (in the case of a single scattering event). For scattering by a small volume containing an ensemble of particles, the ensemble-averaged Mueller matrix is referred to as the Stokes scattering matrix F. Finally, when transforming from the scattering plane to a fixed laboratory frame, the corresponding matrix is referred to as the scattering phase matrix P.

2.4.1 Stokes Vector Representation I_S = [I,Q,U,V]^T

The scattering geometry is illustrated in Figure 1. The plane AOB, defined as the scattering plane, is spanned by the directions of propagation of the incident parallel beam with Stokes vector ${\mathbf{I}}_S^{\text{inc}}$ and the scattered parallel beam with Stokes vector ${\mathbf{I}}_S^{\text{sca}}$. Here the subscript S pertains to the Stokes vector representation I_S = [I,Q,U,V]^T. The scattered radiation, represented by the Stokes vector ${\mathbf{I}}_S^{\text{sca}}$, is related to the incident radiation, represented by the Stokes vector ${\mathbf{I}}_S^{\text{inc}}$, by a 4 × 4 scattering matrix [see Eqs. 20 and 21 below].

If in a small volume of particles any of the following conditions are met (Hovenier and van der Mee, 1983) 1) each particle in the volume element has a plane of symmetry, and the particles are randomly oriented, 2) each volume element contains an equal number of particles and their mirror particles in random orientation, 3) the particles are much smaller than the wavelength of the incident light, then the Stokes scattering matrix in the I_S = [I,Q,U,V]^T representation has the following form

${\mathbf{F}}_S\left(\mathrm{\Theta }\right)=\left[\begin{array}{cccc}\hfill a_1\left(\mathrm{\Theta }\right)\hfill & \hfill b_1\left(\mathrm{\Theta }\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill b_1\left(\mathrm{\Theta }\right)\hfill & \hfill a_2\left(\mathrm{\Theta }\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill a_3\left(\mathrm{\Theta }\right)\hfill & \hfill b_2\left(\mathrm{\Theta }\right)\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -b_2\left(\mathrm{\Theta }\right)\hfill & \hfill a_4\left(\mathrm{\Theta }\right)\hfill \end{array}\right].(20)$

Each of the six independent matrix elements in Eq. 20 depends on the scattering angle Θ, and will in general also depend on the position in the medium. For spherical particles, the matrix in Eq. 20 simplifies, since a₁ = a₂ and a₃ = a₄, so that only four independent elements remain.

Two rotations are required to connect the Stokes vector of the scattered radiation to that of the incident radiation. As illustrated in Figure 1, the first rotation is from the meridian plane OAC, associated with the Stokes vector ${\mathbf{I}}_S^{\text{inc}}$, into the scattering plane OAB, whereas the second rotation is from the scattering plane OAB into the meridian plane OBC, associated with the Stokes vector ${\mathbf{I}}_S^{\text{sca}}$. Hence, the Stokes vector for the scattered radiation is given by (Chandrasekhar, 1960)

${\mathbf{I}}_S^{\text{sca}}={\mathbf{R}}_S\left(\pi -i_2\right){\mathbf{F}}_S\left(\mathrm{\Theta }\right){\mathbf{R}}_S\left(-i_1\right){\mathbf{I}}_S^{\text{inc}}\equiv {\mathbf{P}}_S\left(\mathrm{\Theta }\right){\mathbf{I}}_S^{\text{inc}}.(21)$

Here i₁ and i₂ are the angles between the meridian planes of the incident and the scattered radiation, respectively, and the scattering plane (Figure 1).The Stokes rotation matrix R_S represents a rotation in the clockwise direction with respect to an observer looking into the direction of propagation (Chandrasekhar, 1960). For rotation by an arbitrary angle of ω (0 ≤ ω ≤ 2π) it can be written as

${\mathbf{R}}_S\left(\omega \right)=\left[\begin{array}{cccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \cos \left(2\omega \right)\hfill & \hfill -\sin \left(2\omega \right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \sin \left(2\omega \right)\hfill & \hfill \cos \left(2\omega \right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right].(22)$

Hence, according to Eq. 21, the scattering phase matrix, which connects the Stokes vector of the scattered radiation to that of the incident radiation, is obtained from the Stokes scattering matrix F_S(Θ) in Eq. 20 by

${\mathbf{P}}_S\left({\theta }^{\prime },{\varphi }^{\prime };\theta ,\varphi \right)={\mathbf{R}}_S\left(\pi -i_2\right){\mathbf{F}}_S\left(\mathrm{\Theta }\right){\mathbf{R}}_S\left(-i_1\right)={\mathbf{R}}_S\left(-i_2\right){\mathbf{F}}_S\left(\mathrm{\Theta }\right){\mathbf{R}}_S\left(-i_1\right)(23)$

where R_S (π − i₂) = R_S (−i₂) since the rotation matrix is periodic with a period π.

According to Eq. 21 (see also Figure 1), the Stokes vector ${\mathbf{I}}_S^{\text{inc}}$ of the incident parallel beam must be multiplied by the rotation matrix R_S (−i₁) before it is multiplied by the Stokes scattering matrix F_S(Θ), whereafter it must be multiplied by the rotation matrix R_S (π − i₂). In some radiative transfer (RT) models including Monte Carlo simulations these matrix multiplications are carried out explicitly. In other types of RT models such as the adding-doubling method (De Haan et al., 1987) and the discrete ordinate method (Siewert, 2000; Sommersten et al., 2010; Cohen et al., 2013) they are taken care of implicitly through the expansion of the scattering phase matrix in generalized spherical functions (Siewert, 1981, 1982) as discussed in Section 2.4.3.

Carrying out the matrix multiplications in Eq. 23 one finds:

${\mathbf{P}}_S\left(\mathrm{\Theta }\right)=\left[\begin{array}{cccc}\hfill a_1\hfill & \hfill b_1{C}_1\hfill & \hfill -b_1{S}_1\hfill & \hfill 0\hfill \\ \hfill b_1{C}_2\hfill & \hfill C_2{a}_2{C}_1-S_2{a}_3{S}_1\hfill & \hfill -C_2{a}_2{S}_1-S_2{a}_3{C}_1\hfill & \hfill -b_2{S}_2\hfill \\ \hfill b_1{S}_2\hfill & \hfill S_2{a}_2{C}_1+C_2{a}_3{S}_1\hfill & \hfill -S_2{a}_2{S}_1+C_2{a}_3{C}_1\hfill & \hfill -b_2{C}_2\hfill \\ \hfill 0\hfill & \hfill -b_2{S}_1\hfill & \hfill -b_2{C}_1\hfill & \hfill a_4\hfill \end{array}\right](24)$

where a_j = a_j(Θ), j = 1, … , 4, b_j = b_j(Θ), j = 1, 2, and

$C_1=\cos 2i_1,C_2=\cos 2i_2(25)$

$S_1=\sin 2i_1,S_2=\sin 2i_2.(26)$

A comparison of Eqs. 20 and 24 shows that only the corner elements of F_S(Θ) remain unchanged by the rotations of the reference planes. The (1.1)-element of both the scattering phase matrix P_S(Θ) and the Stokes scattering matrix F_S(Θ) is the scattering phase function. Also, since the (4.4)-element of the scattering phase matrix remains unchanged by the rotations, the state of circular polarization of the incident light does not affect the intensity of the scattered radiation after one scattering event.

To compute P_S(θ′, ϕ′; θ, ϕ) given by Eq. 23 we must relate the angles θ′, ϕ′, θ, and ϕ on the left side with the angles i₁, i₂, and Θ on the right side. Using spherical geometry, we may apply the cosine rule for Θ, θ, and θ′ successively, in Figure 1, to obtain (u = cos θ, u′ = cos θ′) (Hovenier et al., 2004)

$\cos \mathrm{\Theta }=u{u}^{\prime }+{\left(1-u^2\right)}^{1/2}{\left(1-u^{\prime 2}\right)}^{1/2}\cos \left({\varphi }^{\prime }-\varphi \right)(27)$

$\cos i_1=\frac{-u+u^{\prime }\cos \mathrm{\Theta }}{{\left(1-u^{\prime 2}\right)}^{1/2}{\left(1-{\cos }^2\mathrm{\Theta }\right)}^{1/2}}(28)$

$\cos i_2=\frac{-u^{\prime }+u\cos \mathrm{\Theta }}{{\left(1-u^2\right)}^{1/2}{\left(1-{\cos }^2\mathrm{\Theta }\right)}^{1/2}}.(29)$

The trigonometric functions for the double angles can be obtained by using

$\cos 2i=2{\cos }^{2}i-1(30)$

and

$\sin 2i=2\sin i\cos i(31)$

or

$\sin 2i=\left\{\begin{array}{ccc}\hfill 2{\left(1-{\cos }^{2}i\right)}^{1/2}\cos i\hfill & \mathrm{i}\mathrm{f}\hfill & 0<{\varphi }^{\prime }-\varphi <\pi \hfill \\ \hfill -2{\left(1-{\cos }^{2}i\right)}^{1/2}\cos i\hfill & \mathrm{i}\mathrm{f}\hfill & \pi <{\varphi }^{\prime }-\hfill \end{array}\right.(32)$

where i is i₁ or i₂.

Equations 25–32 describe the conventional way to determine the variables C₁, S₁, C₂, and S₂.

A better approach to compute the variables C₁, S₁, C₂, and S₂ appearing in P_S(θ′, ϕ′; θ, ϕ) given by Eq. 24 is described in a recent publication (Berk, 2022). This new approach can be described as follows. Defining

$\mathrm{\Delta }\theta ={\theta }^{\prime }-\theta \text{and}\mathrm{\Delta }\varphi ={\varphi }^{\prime }-\varphi ,(33)$

we have

${\sin }^2\mathrm{\Theta }={\sin }^2\mathrm{\Delta }\theta +4\cos \mathrm{\Delta }\theta \sin \theta \sin {\theta }^{\prime }{\sin }^2\left(\mathrm{\Delta }\varphi /2\right)-4{\sin }^2\theta {\sin }^2{\theta }^{\prime }{\sin }^4\left(\mathrm{\Delta }\varphi /2\right),(34)$

and it can be shown that the variables C₁, S₁, C₂, and S₂ are given by

$C_1=\left\{\begin{array}{cc}\hfill 1\hfill & \mathit{\text{if}}\sin \mathrm{\Theta }=0\hfill \\ \hfill \frac{{\left[2\sin \theta \cos {\theta }^{\prime }{\sin }^2\left(\mathrm{\Delta }\varphi /2\right)+\sin \mathrm{\Delta }\theta \right]}^2-{\sin }^2\theta {\sin }^2\mathrm{\Delta }\varphi }{{\sin }^2\mathrm{\Theta }}\hfill & \mathit{\text{otherwise}}\hfill \end{array}\right.(35)$

$S_1=\left\{\begin{array}{cc}\hfill 0\hfill & \mathit{\text{if}}\sin \mathrm{\Theta }=0\hfill \\ \hfill \frac{\left[-2\sin \theta \sin \mathrm{\Delta }\varphi \right.\left[2\sin \theta \cos {\theta }^{\prime }{\sin }^2\left(\mathrm{\Delta }\varphi /2\right)+\sin \mathrm{\Delta }\theta \right]}{{\sin }^2\mathrm{\Theta }}\hfill & \mathit{\text{otherwise}}\hfill \end{array}\right.(36)$

$C_2=\left\{\begin{array}{cc}\hfill 1\hfill & \mathit{\text{if}}\sin \mathrm{\Theta }=0\hfill \\ \hfill \frac{{\left[2\sin {\theta }^{\prime }\cos \theta {\sin }^2\left(\mathrm{\Delta }\varphi /2\right)-\sin \mathrm{\Delta }\theta \right]}^2-{\sin }^2{\theta }^{\prime }{\sin }^2\mathrm{\Delta }\varphi }{{\sin }^2\mathrm{\Theta }}\hfill & \mathit{\text{otherwise}}\hfill \end{array}\right.(37)$

$S_2=\left\{\begin{array}{cc}\hfill 0\hfill & \mathit{\text{if}}\sin \mathrm{\Theta }=0\hfill \\ \hfill \frac{\left[-2\sin {\theta }^{\prime }\sin \mathrm{\Delta }\varphi \right.\left[2\sin {\theta }^{\prime }\cos \theta {\sin }^2\left(\mathrm{\Delta }\varphi /2\right)-\sin \mathrm{\Delta }\theta \right]}{{\sin }^2\mathrm{\Theta }}\hfill & \mathit{\text{otherwise}}.\hfill \end{array}\right.(38)$

The advantage of using Eqs. 33–38 instead of Eqs. 25–32 is that they eliminate numerical instability issues and the need to treat positive and negative relative azimuth angles as separate cases (Berk, 2022).

We now have all the information needed to compute the scattering phase matrix [see Eq. 24] as a function of the three variables u = cos θ, u′ = cos θ′, and Δϕ = ϕ′ − ϕ. If there is no difference in azimuth (i.e. ϕ′ − ϕ = 0), then the meridian planes of the incident and scattered beams in Fig. 1 coincide with the scattering plane. Hence, there is no need to rotate the reference planes (R( − i₂) and R( − i₁) both reduce to the identity matrix), so that

${\mathbf{P}}_S\left(u^{\prime },u,0\right)={\mathbf{P}}_S\left(u^{\prime },u,\pi \right)={\mathbf{F}}_S\left(\mathrm{\Theta }\right).(39)$

It follows from the cosine law of spherical geometry

$\cos \mathrm{\Theta }={\left.u{u}^{\prime }+{\left(1-u^2\right)}^{1/2}1-u^{\prime 2}\right)}^{1/2}\cos \left({\varphi }^{\prime }-\varphi \right)(40)$

that the phase matrix is invariant to three basic changes in the polar angles u′ and u and azimuthal angles ϕ′ and ϕ which leave the scattering angle unaltered: 1) changing the signs of u and u′ simultaneously: P_S(−u′, − u, ϕ′ − ϕ) = P_S(u′, u, ϕ′ − ϕ), 2) interchange of u and u′: P_S(u′, u, ϕ′ − ϕ) = P_S(u, u′, ϕ′ − ϕ) 3) interchange of ϕ and ϕ′: P_S(u′, u, ϕ′ − ϕ) = P_S(u′, u, ϕ − ϕ′). Also, if the b₂-element in Eq. 24 is zero, then the Stokes parameter V is scattered independently of the others, according to the phase function a₄(Θ), and the remaining part of the scattering phase matrix referring to I, Q, and U becomes a 3 × 3 matrix:

${\mathbf{P}}_S\left(\mathrm{\Theta }\right)=\left[\begin{array}{ccc}\hfill a_1\hfill & \hfill b_1{C}_1\hfill & \hfill -b_1{S}_1\hfill \\ \hfill b_1{C}_2\hfill & \hfill C_2{a}_2{C}_1-S_2{a}_3{S}_1\hfill & \hfill -C_2{a}_2{S}_1-S_2{a}_3{C}_1\hfill \\ \hfill b_1{S}_2\hfill & \hfill S_2{a}_2{C}_1+C_2{a}_3{S}_1\hfill & \hfill -S_2{a}_2{S}_1+C_2{a}_3{C}_1\hfill \end{array}\right].(41)$

Finally, in a plane-parallel or slab geometry, there is no azimuth-dependence for light beams traveling in directions perpendicular to the slab (either up or down). Thus, if either the incident or the scattered beam travels in a perpendicular direction, we may use the meridian plane of the other beam as a reference plane for both beams. Since this plane coincides with the scattering plane, Eq. 39 applies in this situation too.

For Rayleigh scattering with parameter $f=\frac{1-\rho }{1+\rho }$, where ρ is the depolarization factor defined in Eq. 53, the Stokes scattering matrix in the Stokes vector representation I_S = [I,Q,U,V]^T is given by (Chandrasekhar, 1960; Sommersten et al., 2010)

${\mathbf{F}}_S\left(\mathrm{\Theta }\right)=\frac{3}{3+f}\left[\begin{array}{cccc}\hfill 1+f{\cos }^2\mathrm{\Theta }\hfill & \hfill -f{\sin }^2\mathrm{\Theta }\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -f{\sin }^2\mathrm{\Theta }\hfill & \hfill f\left(1+{\cos }^2\mathrm{\Theta }\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 2f\cos \mathrm{\Theta }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \left(3f-1\right)\cos \mathrm{\Theta }\hfill \end{array}\right].(42)$

For the first scattering event of the scalar RTE, only the (1.1)-element of Eq. 42 matters, and leads to the scattering phase function given by

$p_{\text{Ray}}\left(\mathrm{\Theta }\right)=\frac{3}{3+f}\left(1+f{\cos }^2\mathrm{\Theta }\right).(43)$

2.4.2 Stokes Vector Representation $\mathbf{I}={[I_{\Vert },I_{\perp },U,V]}^T$

The Stokes vector $\mathbf{I}={[I_{\Vert },I_{\perp },U,V]}^T$ is related to I_S = [I,Q,U,V]^T by

${\mathbf{I}}_S=\mathbf{D}\mathbf{I}(44)$

where D is given by Eq. 14, so that I = I_‖ + I_⊥, and Q = I_‖ − I_⊥. Denoting the Stokes vector obtained after a rotation by

${\mathbf{I}}_S^{\prime }={\mathbf{R}}_S\left(\omega \right){\mathbf{I}}_S(45)$

we find

${\mathbf{I}}^{\prime }={\mathbf{D}}^{-1}{\mathbf{I}}_S^{\prime }={\mathbf{D}}^{-1}{\mathbf{R}}_S\left(\omega \right){\mathbf{I}}_S={\mathbf{D}}^{-1}{\mathbf{R}}_S\left(\omega \right)\mathbf{D}\mathbf{I}=\mathbf{R}\left(\omega \right)\mathbf{I}.(46)$

Hence, the rotation matrix for the Stokes vector in the representation $\mathbf{I}={[I_{\Vert },I_{\perp },U,V]}^T$ becomes:

$\mathbf{R}\left(\omega \right)={\mathbf{D}}^{-1}{\mathbf{R}}_S\left(\omega \right)\mathbf{D}=\left[\begin{array}{cccc}\hfill {\cos }^2\omega \hfill & \hfill {\sin }^2\omega \hfill & \hfill -\frac{1}{2}\sin \left(2\omega \right)\hfill & \hfill 0\hfill \\ \hfill {\sin }^2\omega \hfill & \hfill {\cos }^2\omega \hfill & \hfill \frac{1}{2}\sin \left(2\omega \right)\hfill & \hfill 0\hfill \\ \hfill \sin \left(2\omega \right)\hfill & \hfill -\sin \left(2\omega \right)\hfill & \hfill \cos \left(2\omega \right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right].(47)$

The scattering phase matrix P(Θ) in the Stokes vector representation $\mathbf{I}={[I_{\Vert },I_{\perp },U,V]}^T$ is related to the scattering phase matrix P_S(Θ) in the Stokes vector representation I_S = [I,Q,U,V]^T by

$\mathbf{P}\left(\mathrm{\Theta }\right)={\mathbf{D}}^{-1}{\mathbf{P}}_S\left(\mathrm{\Theta }\right)\mathbf{D}.(48)$

Similarly, the Stokes scattering matrix F(Θ) associated with the Stokes vector representation $\mathbf{I}={[I_{\Vert },I_{\perp },U,V]}^T$ is related to the Stokes scattering matrix F_S(Θ) in Eq. 20 by

$\mathbf{F}\left(\mathrm{\Theta }\right)={\mathbf{D}}^{-1}{\mathbf{F}}_S\left(\mathrm{\Theta }\right)\mathbf{D}=\left(\begin{array}{cccc}\hfill \frac{1}{2}\left(a_1+a_2+2b_1\right)\hfill & \hfill \frac{1}{2}\left(a_1-a_2\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \frac{1}{2}\left(a_1-a_2\right)\hfill & \hfill \frac{1}{2}\left(a_1+a_2-2b_1\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill a_3\hfill & \hfill b_2\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -b_2\hfill & \hfill a_4\hfill \end{array}\right).(49)$

In the Stokes vector representation $\mathbf{I}={[I_{\Vert },I_{\perp },U,V]}^T$, the Stokes scattering matrix for Rayleigh scattering becomes (using Eqs. 42 and 49 (Chandrasekhar, 1960)):

$\mathbf{F}\left(\mathrm{\Theta }\right)=\frac{3}{2\left(1+2\zeta \right)}\left(\begin{array}{cccc}\hfill {\cos }^2\mathrm{\Theta }+\zeta {\sin }^2\mathrm{\Theta }\hfill & \hfill \zeta \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \zeta \hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \left(1-\zeta \right)\cos \mathrm{\Theta }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \left(1-3\zeta \right)\cos \mathrm{\Theta }\hfill \end{array}\right)(50)$

where $\zeta =\rho /(2-\rho )=\frac{1-f}{1+3f}$.

From Eq. 50 we see that for an incident beam of natural unpolarized light given by ${\mathbf{I}}^{\text{inc}}={[I_{\parallel }^{\text{inc}},I_{\perp }^{\text{inc}},U^{\text{inc}},V^{\text{inc}}]}^T={[\frac{1}{2}{I}^{\text{inc}},\frac{1}{2}{I}^{\text{inc}},\mathrm{0,0}]}^T$, the scattered intensities in the plane parallel and perpendicular to the scattering plane are obtained by carrying out the multiplication I^sca = F(Θ)I^inc:

$I_{\parallel }^{\text{sca}}\propto \frac{3}{4\left(1+2\zeta \right)}\left[2\zeta +\left(1-\zeta \right){\cos }^2\mathrm{\Theta }\right]I^{\text{inc}}(51)$

$I_{\perp }^{\text{sca}}\propto \frac{3}{4\left(1+2\zeta \right)}\left[\left(1+\zeta \right)\right]I^{\text{inc}}.(52)$

Thus, for unpolarized incident light, the scattered light at right angles (Θ = 90°) to the direction of incidence defines the depolarization ratio:

$\rho \equiv {\left(\frac{I_{\parallel }^{\text{sca}}}{I_{\perp }^{\text{sca}}}\right)}_{\mathrm{\Theta }=90°}=\frac{2\zeta }{1+\zeta }(53)$

whereas the degree of linear polarization becomes [Eq. 18]:

$p_l=\frac{I_{\perp }-I_{\parallel }}{I_{\perp }+I_{\parallel }}=\frac{\left(1-\zeta \right)\left(1-{\cos }^2\mathrm{\Theta }\right)}{1+3\zeta +\left(1-\zeta \right){\cos }^2\mathrm{\Theta }}\to \frac{1-\zeta }{1+3\zeta }=\frac{1-\rho }{1+\rho }=f\mathrm{a}\mathrm{s}\mathrm{\Theta }\to 90°.(54)$

2.4.3 Generalized Spherical Functions–The Greek Constants

For the scalar RTE, only the a₁(Θ) element of the Stokes scattering matrix Eq. 20 is relevant, and this element is the scattering phase function given by Eq. 69 in general, and by Eq. 43 for Rayleigh scattering. The scattering phase function can be expanded in Legendre polynomials [see Eq. 69], which enables expression as a Fourier cosine series.

In a similar manner, the scattering phase matrix can be expanded in generalized spherical functions. In the Stokes vector representation I_S = [I,Q,U,V]^T, the scattering phase matrix is P_S(Θ) = P_S(u′, u; ϕ′ − ϕ) with u = cos θ, θ being the polar angle after scattering, and u′ = cos θ′, θ′ being the polar angle prior to scattering. Similarly, ϕ and ϕ′ are the azimuth angles after and prior to scattering, respectively. To expand in generalized spherical functions, the scattering phase matrix is first expanded in a (M + 1)-term Fourier series in the azimuth angle difference (Δϕ′ = ϕ′ − ϕ):

${\mathbf{P}}_S\left(u^{\prime },u;\mathrm{\Delta }\varphi \right)=\sum _{m=0}^M\left\{{\mathbf{P}}_c^m\left(u^{\prime },u\right)\cos m\left(\mathrm{\Delta }{\varphi }^{\prime }\right)+{\mathbf{P}}_s^m\left(u^{\prime },u\right)\sin m\left(\mathrm{\Delta }{\varphi }^{\prime }\right)\right\}(55)$

where ${\mathbf{P}}_c^m(u^{\prime },u)$ and ${\mathbf{P}}_s^m(u^{\prime },u)$ are the coefficient matrices of the cosine and sine terms, respectively, of the Fourier series.

We use an addition theorem for the generalized spherical functions to express the Fourier expansion coefficient matrices directly in terms of the expansion coefficients of the Stokes scattering matrix F_S(Θ) [see Eq. 20] as follows (Siewert, 1981; Siewert, 1982; Mishchenko, 1991):

${\mathbf{P}}_c^m\left(u^{\prime },u\right)={\mathbf{A}}^m\left(u^{\prime },u\right)+{\mathbf{\Delta }}_{\mathrm{3,4}}{\mathbf{A}}^m\left(u^{\prime },u\right){\mathbf{\Delta }}_{\mathrm{3,4}}(56)$

${\mathbf{P}}_s^m\left(u^{\prime },u\right)={\mathbf{A}}^m\left(u^{\prime },u\right){\mathbf{\Delta }}_{\mathrm{3,4}}-{\mathbf{\Delta }}_{\mathrm{3,4}}{\mathbf{A}}^m\left(u^{\prime },u\right)(57)$

where Δ_3,4 = diag (1, 1, −1, 1). The matrix A^m (u′, u) is given by:

${\mathbf{A}}^m\left(u^{\prime },u\right)=\sum _{\ell =m}^M{\mathbf{P}}_{\ell }^m\left(u\right){\mathbf{\Lambda }}_{\ell }{\mathbf{P}}_{\ell }^m\left(u^{\prime }\right).(58)$

The matrix ${\mathbf{P}}_{\ell }^m(u)$ is given by:

${\mathbf{P}}_{\ell }^m\left(u\right)=\left(\begin{array}{cccc}\hfill P_{\ell }^{m,0}\left(u\right)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill P_{\ell }^{m,+}\left(u\right)\hfill & \hfill P_{\ell }^{m,-}\left(u\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill P_{\ell }^{m,-}\left(u\right)\hfill & \hfill P_{\ell }^{m,+}\left(u\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill P_{\ell }^{m,0}\left(u\right)\hfill \end{array}\right)(59)$

where

$P_{\ell }^{m,\pm }\left(u\right)=\frac{1}{2}\left[P_{\ell }^{m,-2}\left(u\right)\pm P_{\ell }^{m,2}\left(u\right)\right](60)$

and the functions $P_{\ell }^{m,0}(u)$ and $P_{\ell }^{m,\pm 2}(u)$ are the generalized spherical functions (Hovenier et al., 2004). The matrix Λ_ℓ in Eq. 58 is

${\mathbf{\Lambda }}_{\ell }=\left(\begin{array}{cccc}\hfill {\alpha }_{1,\ell }\hfill & \hfill {\beta }_{1,\ell }\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\beta }_{1,\ell }\hfill & \hfill {\alpha }_{2,\ell }\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\alpha }_{3,\ell }\hfill & \hfill {\beta }_{2,\ell }\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -{\beta }_{2,\ell }\hfill & \hfill {\alpha }_{4,\ell }\hfill \end{array}\right)(61)$

and

$a_1\left(\mathrm{\Theta }\right)=\sum _{\ell =0}^M{\alpha }_{1,\ell }{P}_{\ell }^{\mathrm{0,0}}\left(\cos \mathrm{\Theta }\right)(62)$

$a_2\left(\mathrm{\Theta }\right)+a_3\left(\mathrm{\Theta }\right)=\sum _{\ell =2}^M\left({\alpha }_{2,\ell }+{\alpha }_{3,\ell }\right)P_{\ell }^{\mathrm{2,2}}\left(\cos \mathrm{\Theta }\right)(63)$

$a_2\left(\mathrm{\Theta }\right)-a_3\left(\mathrm{\Theta }\right)=\sum _{\ell =2}^M\left({\alpha }_{2,\ell }-{\alpha }_{3,\ell }\right)P_{\ell }^{2,-2}\left(\cos \mathrm{\Theta }\right)(64)$

$a_4\left(\mathrm{\Theta }\right)=\sum _{\ell =0}^M{\alpha }_{4,\ell }{P}_{\ell }^{\mathrm{0,0}}\left(\cos \mathrm{\Theta }\right)(65)$

$b_1\left(\mathrm{\Theta }\right)=\sum _{\ell =2}^M{\beta }_{1,\ell }{P}_{\ell }^{\mathrm{0,2}}\left(\cos \mathrm{\Theta }\right)(66)$

$b_2\left(\mathrm{\Theta }\right)=\sum _{\ell =2}^M{\beta }_{2,\ell }{P}_{\ell }^{\mathrm{0,2}}\left(\cos \mathrm{\Theta }\right).(67)$

Here the Greek constants α_j,ℓ and β_j,ℓ are expansion coefficients, and a_j(Θ) and b_j(Θ) are the elements of the Stokes scattering matrix F_S(Θ) in Eq. 20. An example of Greek constants for Rayleigh scattering is provided in Table 1 (see Mishchenko and Travis (1997)) where

$c=\frac{2\left(1-\rho \right)}{2+\rho }d=\frac{2\left(1-2\rho \right)}{2+\rho }(68)$

TABLE 1

TABLE 1. Expansion coefficients for Rayleigh scattering.

and ρ is the depolarization ratio given by Eq. 53.

We note that in the scalar (unpolarized) case all components of the Stokes scattering matrix F_S(Θ) [see Eq. 20] are zero except for a₁(Θ), and:

$a_1\left(\mathrm{\Theta }\right)=\sum _{\ell =0}^M{\alpha }_{1,\ell }\left(\tau \right)P_{\ell }^{\mathrm{0,0}}\left(\cos \mathrm{\Theta }\right)\equiv p\left(\tau ,\cos \mathrm{\Theta }\right)\approx \sum _{\ell =0}^M\left(2\ell +1\right){\chi }_{\ell }\left(\tau \right)P_{\ell }\left(\cos \mathrm{\Theta }\right)(69)$

since $P_{\ell }^{\mathrm{0,0}}(\cos \mathrm{\Theta })\equiv P_{\ell }(\cos \mathrm{\Theta })$, where P_ℓ(cos Θ) is the Legendre polynomial of order ℓ, and α_1,ℓ(τ) ≡ (2ℓ + 1)χ_ℓ(τ). Here the coefficients χ_ℓ(τ) are the moments of the phase function expanded in Legendre polynomials. Note also that the expansion coefficients given above [Eq. 61] are for the scattering phase matrix P_S(Θ), which relates the incident and scattered Stokes vectors in the representation I_S = [I,Q,U,V]^T.

3 Solution of the Vector Radiative Transfer Equation

3.1 Isolation of Azimuth Dependence

We start from the scattering phase matrix expanded in a Fourier series (see Eq. 55) (Δϕ′ = ϕ′ − ϕ):

$\mathbf{P}\left(u^{\prime },u;\mathrm{\Delta }{\varphi }^{\prime }\right)=\sum _{m=0}^M\left\{{\mathbf{P}}_{\text{c}}^m\left(u^{\prime },u\right)\cos m\mathrm{\Delta }{\varphi }^{\prime }+{\mathbf{P}}_{\text{s}}^m\left(u^{\prime },u\right)\sin m\mathrm{\Delta }{\varphi }^{\prime }\right\}.(70)$

To isolate the azimuth dependence of the radiation field we expand the Stokes vector I (τ, u, ϕ) in the VRTE [Eq. 7] and the source term Q (τ, u, ϕ) in Eq. 9 in a Fourier series in a manner similar to the expansion of the scattering phase matrix in Eq. 70 (Δϕ₀ = ϕ₀ − ϕ):

$\mathbf{I}\left(\tau ,u,\varphi \right)=\sum _{m=0}^M\left\{{\mathbf{I}}_{\text{c}}^m\left(\tau ,u\right)\cos m\mathrm{\Delta }{\varphi }_0+{\mathbf{I}}_{\text{s}}^m\left(\tau ,u\right)\sin m\mathrm{\Delta }{\varphi }_0\right\}(71)$

$\begin{array}{ccc}\hfill & \hfill & \hfill \\ \hfill \mathbf{Q}\left(\tau ,u,\varphi \right)& =\hfill & \sum _{m=0}^M\left\{{\mathbf{Q}}_{\text{c\hspace{0.17em}}}^m\left(\tau ,u\right)\cos m\mathrm{\Delta }{\varphi }_0+{\mathbf{Q}}_{\text{s\hspace{0.17em}}}^m\left(\tau ,u\right)\sin m\mathrm{\Delta }{\varphi }_0\right\}\hfill \end{array}(72)$

where the subscript s or c denotes sine or cosine mode. Using these expansions, we obtain the following equations for the Fourier components of the VRTE (see Stamnes and Stamnes (2015) for details)

$\begin{array}{ccc}\hfill u\frac{d{\mathbf{I}}_{\text{c}}^m\left(\tau ,u\right)}{d\tau }& =\hfill & {\mathbf{I}}_{\text{c}}^m\left(\tau ,u\right)-\frac{\varpi \left(\tau \right)}{4}{\int }_{-1}^{1}d{u}^{\prime }\left\{{\mathbf{P}}_{\text{c}}^m\left(\tau ,u^{\prime },u\right){\mathbf{I}}_{\text{c}}^m\left(\tau ,u^{\prime }\right)\left(1+{\delta }_{0m}\right)\right.\hfill \\ \hfill & \hfill & \left.-{\mathbf{P}}_{\text{s}}^m\left(\tau ,u^{\prime },u\right){\mathbf{I}}_{\text{s}}^m\left(\tau ,u^{\prime }\right)\right\}-{\mathbf{Q}}_{\text{c}}^m\left(\tau ,u\right)\hfill \end{array}(73)$

$\begin{array}{ccc}\hfill u\frac{d{\mathbf{I}}_{\text{s}}^m\left(\tau ,u\right)}{d\tau }& =\hfill & {\mathbf{I}}_{\text{s}}^m\left(\tau ,u\right)-\frac{\varpi \left(\tau \right)}{4}{\int }_{-1}^{1}d{u}^{\prime }\left\{{\mathbf{P}}_{\text{c}}^m\left(\tau ,u^{\prime },u\right)\right.{\mathbf{I}}_{\text{s}}^m\left(\tau ,u^{\prime }\right)\hfill \\ \hfill & \hfill & \left.+{\mathbf{P}}_{\text{s}}^m\left(\tau ,u^{\prime },u\right){\mathbf{I}}_{\text{c}}^m\left(\tau ,u^{\prime }\right)\right\}-{\mathbf{Q}}_{\text{s}}^m\left(\tau ,u\right).\hfill \end{array}(74)$

For scattering by randomly oriented particles, the Fourier coefficient matrix ${\mathbf{P}}_{\text{c}}^m(u,u^{\prime })$ has two (2 × 2) zero submatrices, one in the upper right corner and one in the lower left corner, and the matrix ${\mathbf{P}}_{\text{s}}^m(u,u^{\prime })$ has a (2 × 2) zero submatrix in the upper left corner and one in the lower right corner (Hovenier and van der Mee, 1983). Hence

${\mathbf{P}}_{\text{c}}^m\left(u,u^{\prime }\right)=\left(\begin{array}{cccc}\hfill C_{11}^m\hfill & \hfill C_{12}^m\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill C_{21}^m\hfill & \hfill C_{22}^m\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill C_{33}^m\hfill & \hfill C_{34}^m\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill C_{43}^m\hfill & \hfill C_{44}^m\hfill \end{array}\right){\mathbf{P}}_{\text{s}}^m\left(u,u^{\prime }\right)=\left(\begin{array}{cccc}\hfill 0\hfill & \hfill 0\hfill & \hfill S_{13}^m\hfill & \hfill S_{14}^m\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill S_{23}^m\hfill & \hfill S_{24}^m\hfill \\ \hfill S_{31}^m\hfill & \hfill S_{32}^m\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill S_{41}^m\hfill & \hfill S_{42}^m\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right)(75)$

where we have defined ${\mathbf{P}}_{\text{c}}^m{(u,u^{\prime })}_{i,j}\equiv {\mathbf{C}}^m{(u,u^{\prime })}_{i,j}$ and ${\mathbf{P}}_{\text{s}}^m{(u,u^{\prime })}_{i,j}\equiv {\mathbf{S}}^m{(u,u^{\prime })}_{i,j}$ to simplify the notation. Therefore, the homogeneous VRTE [Q = 0 in Eqs. 73 and 74] may be rewritten as (m ∈ [0, 1, 2, … , 2N])

$u\frac{d{\mathbf{I}}_{\text{c}}^m\left(\tau ,u\right)}{d\tau }={\mathbf{I}}_{\text{c}}^m\left(\tau ,u\right)-\frac{\varpi \left(\tau \right)}{4}{\int }_{-1}^1\left[{\mathbf{C}}^m\left(u,u^{\prime }\right){\mathbf{I}}_{\text{c}}^m\left(\tau ,u^{\prime }\right)-\left(1-{\delta }_{0m}\right){\mathbf{S}}^m\left(u,u^{\prime }\right){\mathbf{I}}_{\text{s}}^m\left(\tau ,u^{\prime }\right)\right]d{u}^{\prime }(76)$

$u\frac{d{\mathbf{I}}_{\text{s}}^m\left(\tau ,u\right)}{d\tau }={\mathbf{I}}_{\text{s}}^m\left(\tau ,u\right)-\frac{\varpi \left(\tau \right)}{4}\left(1-{\delta }_{0m}\right){\int }_{-1}^1\left[{\mathbf{C}}^m\left(u,u^{\prime }\right){\mathbf{I}}_{\text{s}}^m\left(\tau ,u^{\prime }\right)+{\mathbf{S}}^m\left(u,u^{\prime }\right){\mathbf{I}}_{\text{c}}^m\left(\tau ,u^{\prime }\right)\right]d{u}^{\prime }.(77)$

3.1.1 Vector Radiative Transfer Equation for the Combined Mode

We note that the ${\mathbf{I}}_{\text{c}}^m(\tau ,u)$ and ${\mathbf{I}}_{\text{s}}^m(\tau ,u)$ components in Eqs. 76 and 77 are still coupled. To produce a pair of independent differential equations, we define combined cosine and sine modes as

${\tilde{\mathbf{I}}}_{\text{c}}^m\left(\tau ,u\right)\equiv \left(\begin{array}{c}\hfill I_{\Vert c}^m\left(\tau ,u\right)\hfill \\ \hfill I_{\text{⊥}c}^m\left(\tau ,u\right)\hfill \\ \hfill U_s^m\left(\tau ,u\right)\hfill \\ \hfill V_s^m\left(\tau ,u\right)\hfill \end{array}\right){\tilde{\mathbf{I}}}_s^m\left(\tau ,u\right)\equiv \left(\begin{array}{c}\hfill I_{\Vert s}^m\left(\tau ,u\right)\hfill \\ \hfill I_{\text{⊥}s}^m\left(\tau ,u\right)\hfill \\ \hfill U_c^m\left(\tau ,u\right)\hfill \\ \hfill V_c^m\left(\tau ,u\right)\hfill \end{array}\right).(78)$

After substantial manipulations, one can show that the mth term of the homogeneous VRTE of the combined modes can now be written as (m ∈ [0, 1, 2, … , 2N]):

$u\frac{d{\tilde{\mathbf{I}}}_c^m\left(\tau ,u\right)}{d\tau }={\tilde{\mathbf{I}}}_c^m\left(\tau ,u\right)-\frac{\varpi \left(\tau \right)}{2}{\int }_{-1}^1{\tilde{\mathbf{P}}}_c^m\left(u,u^{\prime }\right){\tilde{\mathbf{I}}}_c^m\left(\tau ,u^{\prime }\right)d{u}^{\prime }(79)$

$u\frac{d{\tilde{\mathbf{I}}}_s^m\left(\tau ,u\right)}{d\tau }={\tilde{\mathbf{I}}}_s^m\left(\tau ,u\right)-\frac{\varpi \left(\tau \right)}{2}{\int }_{-1}^1{\tilde{\mathbf{P}}}_s^m\left(u,u^{\prime }\right){\mathbf{I}}_s^m\left(\tau ,u^{\prime }\right)d{u}^{\prime }(80)$

where the combined scattering phase matrices are defined as:

${\tilde{\mathbf{P}}}_c^m\left(u,u^{\prime }\right)=\left(\begin{array}{cccc}\hfill C_{11}^m\hfill & \hfill C_{12}^m\hfill & \hfill -S_{13}^m\hfill & \hfill -S_{14}^m\hfill \\ \hfill C_{21}^m\hfill & \hfill C_{22}^m\hfill & \hfill -S_{23}^m\hfill & \hfill -S_{24}^m\hfill \\ \hfill S_{31}^m\hfill & \hfill S_{32}^m\hfill & \hfill C_{33}^m\hfill & \hfill C_{34}^m\hfill \\ \hfill S_{41}^m\hfill & \hfill S_{42}^m\hfill & \hfill C_{43}^m\hfill & \hfill C_{44}^m\hfill \end{array}\right){\tilde{\mathbf{P}}}_s^m\left(u,u^{\prime }\right)=\left(\begin{array}{cccc}\hfill C_{11}^m\hfill & \hfill C_{12}^m\hfill & \hfill S_{13}^m\hfill & \hfill S_{14}^m\hfill \\ \hfill C_{21}^m\hfill & \hfill C_{22}^m\hfill & \hfill S_{23}^m\hfill & \hfill S_{24}^m\hfill \\ \hfill -S_{31}^m\hfill & \hfill -S_{32}^m\hfill & \hfill C_{33}^m\hfill & \hfill C_{34}^m\hfill \\ \hfill -S_{41}^m\hfill & \hfill -S_{42}^m\hfill & \hfill C_{43}^m\hfill & \hfill C_{44}^m\hfill \end{array}\right).(81)$

For the special m = 0 case, we have

${\tilde{\mathbf{P}}}_c^0\left(u,u^{\prime }\right)=\left(\begin{array}{cccc}\hfill C_{11}^0\hfill & \hfill C_{12}^0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill C_{21}^0\hfill & \hfill C_{22}^0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right){\tilde{\mathbf{P}}}_s^0\left(u,u^{\prime }\right)=\left(\begin{array}{cccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill C_{33}^0\hfill & \hfill C_{34}^0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill C_{43}^0\hfill & \hfill C_{44}^0\hfill \end{array}\right).(82)$

Equations 79 and 80 are two independent differential equations to be solved.

3.2 The Discrete Ordinate Method

The discrete ordinate method consists of replacing the integration over u′ in Eqs. 79 and 80 by a discrete sum by introducing the Gaussian quadrature points u_j (the discrete ordinates) and corresponding weights w_j. For each Fourier component one obtains (i = ±1, ±2, … , ±2N):

$u_i\frac{d{\tilde{\mathbf{I}}}_{\text{c}}^m\left(\tau ,u_i\right)}{d\tau }={\tilde{\mathbf{I}}}_{\text{c}}^m\left(\tau ,u_i\right)-\frac{\varpi \left(\tau \right)}{2}\sum _{\begin{array}{c}j=-N\\ j\ne 0\end{array}}^N{\omega }_j{\tilde{\mathbf{P}}}_{\text{c}}^m\left(u_i,u_j\right){\tilde{\mathbf{I}}}_{\text{c}}^m\left(\tau ,u_j\right)(83)$

$u_i\frac{d{\tilde{\mathbf{I}}}_{\text{s}}^m\left(\tau ,u_i\right)}{d\tau }={\tilde{\mathbf{I}}}_{\text{s}}^m\left(\tau ,u_i\right)-\frac{\varpi \left(\tau \right)}{2}\sum _{\begin{array}{c}j=-N\\ j\ne 0\end{array}}^N{\omega }_j{\tilde{\mathbf{P}}}_{\text{s}}^m\left(u_i,u_j\right){\mathbf{I}}_{\text{s}}^m\left(\tau ,u_j\right)(84)$

The convention for the indices of the quadrature points is such that u_j < 0 for j < 0, and u_j > 0 for j > 0. These points are distributed symmetrically about zero, i.e., u_−j = −u_j, and the corresponding weights are equal, i.e. w_−j = w_j.

The solution of the discrete ordinate approximation to the VRTE Eqs. 83 and 84 is analogous to that of the scalar RTE. Detailed derivations, including the removal of the notorious ill-conditioning problem, can be found elsewhere (Schulz et al., 1999; Siewert, 2000; Stamnes and Stamnes, 2015), and will not be repeated here.

3.3 Discrete Ordinate Radiative Transfer Upgrades

Below we describe some new important upgrades of VDISORT:

• First, a new algorithm to handle the complex eigenvalue/eigenvector problem is developed and implemented in VDISORT to give an accurate computation of the V component of the Stokes vector.

• Second, a reduction of the dimension of the complex eigenvalue problem is developed to reduce the computational burden, and the boundary condition for the complex eigenvalue/vector case is discussed.

• Third, to obtain solutions at arbitrary polar angles, we have developed and implemented an accurate ISF approach that works for both real and complex eigensolutions and an arbitrary bidirectional reflectance distribution matrix is added to compute the polarized reflectance at the lower boundary.

• Fourth, a pseudo-spherical treatment has been implemented to provide important corrections for Earth curvature effects at large (near horizontal) solar zenith and observation polar angles.

• Finally, a single-scattering solution is developed and used to enhance the accuracy and speed for strongly forward-peaked scattering, such as by large water droplets or ice particles.

3.3.1 Complex Eigenvalues/Eigenvectors

The appearance of complex eigenvalues/eigenvectors in the vector radiative transfer problem stems from the asymmetric structure of the Stokes scattering matrix [ ± b₂(Θ) in Eq. 20], or equivalently, from the Greek constant matrix [ ± β₂ in Eq. 61]. In the scalar DISORT model, this 4 × 4 scattering matrix degenerates into the scalar scattering phase function Eq. 69 and the complex eigenvalue problem does not occur.

3.3.2 The Complex Homogeneous Solution

By seeking exponential solutions to the homogeneous VRTE’s, Eqs. 83 and 84, one obtains a standard algebraic eigenvalue problem. Writing the VRTE separately for the upward (u > 0) and downward (u < 0) hemispheres, and defining μ = |u| and α = c, s, the homogeneous version of the Fourier components of the VRTE becomes (i = 1, 2, … , 2N):

$\begin{array}{ll}\hfill +{\mu }_i\frac{d{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_i\right)}{d\tau }={\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_i\right)& -\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,-{\mu }_j\right){\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)\hfill \\ \hfill & -\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,+{\mu }_j\right){\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_j\right)\hfill \end{array}(85)$

$\begin{array}{ll}\hfill -{\mu }_i\frac{d{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_i\right)}{dd\tau }={\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_i\right)& -\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(-{\mu }_i,-{\mu }_j\right){\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)\hfill \\ \hfill & -\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(-{\mu }_i,+{\mu }_j\right){\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_j\right).\hfill \end{array}(86)$

Seeking solutions to Eqs. 85 and 86 of the form ${\tilde{\mathbf{I}}}_{\alpha }^m(\tau ,+{\mu }_i)=\mathbf{g}({\mu }_i)\exp (-k\tau )$, we obtain an algebraic eigenvalue problem: Ag = kg, where the 8N × 8N eigenmatrix A can be written as (here the italic I is the identity matrix):

$\left(\begin{array}{cccccc}\hfill \frac{-\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}^m\left({\mu }_1,{\mu }_1\right)}{2{\mu }_1}+\frac{\mathit{I}}{{\mu }_1}\hfill & \hfill \cdots \hfill & \hfill \frac{-\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}^m\left({\mu }_1,{\mu }_N\right)}{2{\mu }_1}\hfill & \hfill \frac{-\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}^m\left({\mu }_1,-{\mu }_1\right)}{2{\mu }_1}\hfill & \hfill \cdots \hfill & \hfill \frac{-\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}^m\left({\mu }_1,-{\mu }_N\right)}{2{\mu }_1}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill \frac{-\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}^m\left({\mu }_N,{\mu }_1\right)}{2{\mu }_N}\hfill & \hfill \cdots \hfill & \hfill \frac{-\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}^m\left({\mu }_N,{\mu }_N\right)}{2{\mu }_N}+\frac{\mathit{I}}{{\mu }_N}\hfill & \hfill \frac{-\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}^m\left({\mu }_N,-{\mu }_1\right)}{2{\mu }_N}\hfill & \hfill \cdots \hfill & \hfill \frac{-\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}^m\left({\mu }_N,-{\mu }_N\right)}{2{\mu }_N}\hfill \\ \hfill \frac{-\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}^m\left(-{\mu }_1,{\mu }_1\right)}{-2{\mu }_1}\hfill & \hfill \cdots \hfill & \hfill \frac{-\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}^m\left(-{\mu }_1,{\mu }_N\right)}{-2{\mu }_1}\hfill & \hfill \frac{-\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}^m\left(-{\mu }_1,-{\mu }_1\right)}{-2{\mu }_1}+\frac{\mathit{I}}{-{\mu }_1}\hfill & \hfill \cdots \hfill & \hfill \frac{-\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}^m\left(-{\mu }_1,-{\mu }_N\right)}{-2{\mu }_1}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill \frac{-\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}^m\left(-{\mu }_N,{\mu }_1\right)}{-2{\mu }_N}\hfill & \hfill \cdots \hfill & \hfill \frac{-\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}^m\left(-{\mu }_N,{\mu }_N\right)}{-2{\mu }_N}\hfill & \hfill \frac{-\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}^m\left(-{\mu }_N,-{\mu }_1\right)}{-2{\mu }_N}\hfill & \hfill \cdots \hfill & \hfill \frac{-\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}^m\left(-{\mu }_N,-{\mu }_N\right)}{-2{\mu }_N}+\frac{\mathit{I}}{-{\mu }_N}\hfill \end{array}\right)(87)$

In the scalar DISORT case real eigenvalues ± k_j and eigenvectors g_±j (±μ_j) are obtained, so that the homogeneous solution can be written as a linear combination of the eigensolutions:

$I_{\text{homogenous}}\left(\tau ,\pm {\mu }_i\right)=\sum _{j=1}^N{C}_{-j}{g}_{-j}\left(\pm {\mu }_i\right)e^{k_j\tau }+\sum _{j=1}^N{C}_j{g}_j\left(\pm {\mu }_i\right)e^{-k_j\tau }.(88)$

For the 4 × 4 VDISORT problem the g_±j (±μ_i) eigenvector in the scalar case must be replaced by a 4 × 1 eigenvector g_±j (±μ_i) for ${[I_{\parallel },I_{\perp },U,V]}^T$. Since we have 2N quadrature angles ± μ₁, ±μ₂, … , ±μ_N, the dimension of the full eigenvector g_±j is 8N × 1. A matrix with real elements has eigenvalue/eigenvector solutions that either are real or occur in complex conjugate pairs. Therefore, defining $k_c,k_c^{\ast }$ as a complex conjugate pair of eigenvalues, ${\mathbf{g}}_c,{\mathbf{g}}_c^{\ast }$ as a complex conjugate pair of eigenvectors, and C _{± 1}, C _{± 2} as arbitrary coefficients, we may, in analogy with the scalar case with only real solutions, write the complex solution to the homogeneous VRTE as:

$\begin{array}{ll}\hfill {\mathbf{I}}_{\text{homo}}\left(\tau ,\pm {\mu }_i\right)=& \underset{{\mathbf{I}}_{\text{REAL}}}{\underbrace{\sum _{j=1}^{N_r}{C}_j{\mathbf{g}}_j\left(\pm {\mu }_i\right)e^{-k_j\tau }+\sum _{j=1}^{N_r}{C}_{-j}{\mathbf{g}}_{-j}\left(\pm {\mu }_i\right)e^{k_j\tau }}}\hfill \\ \hfill & +\sum _{j=1}^{N_i}\left[C_{1j}{\mathbf{g}}_{cj}\left(\pm {\mu }_i\right)e^{-k_{cj}\tau }+C_{2j}{\mathbf{g}}_{cj}^{\ast }\left(\pm {\mu }_i\right)e^{-k_{cj}^{\ast }\tau }\right]\hfill \\ \hfill & +\underset{{\mathbf{I}}_{\text{COMPLEX}}}{\underbrace{\sum _{j=1}^{N_i}\left[C_{-1j}{\mathbf{g}}_{-cj}\left(\pm {\mu }_i\right)e^{k_{cj}^{\ast }\tau }+C_{-2j}{\mathbf{g}}_{-cj}^{\ast }\left(\pm {\mu }_i\right)e^{k_{cj}\tau }\right]}}\hfill \end{array}(89)$

where N_r is the number of real solutions, N_i is the number of complex conjugate pair solutions, and N = N_r + 2N_i is the number of streams in each hemisphere.

The homogeneous solution in Eq. 89 contains complex numbers that must be converted into real values before solving for the coefficients. Since linear combinations of ${\mathbf{g}}_{\pm c}{e}^{\pm k_c\tau }$ and ${\mathbf{g}}_{\text{⊥}c}^{*}{e}^{\pm k_c^{*}\tau }$ are also a solution of Eqs. 85 and 86, we can separate the real and imaginary parts. As shown in Section 8 (Supplementary Appendix A1) the complex homogeneous solutions may be converted into the following real solutions:

$\begin{array}{ll}\hfill {\mathbf{I}}_{\text{homo}}\left(\tau ,\pm {\mu }_i\right)=& \underset{{\mathbf{I}}_{\text{REAL}}}{\underbrace{\sum _{j=1}^{N_r}{C}_j{\mathbf{g}}_j\left(\pm {\mu }_i\right)e^{-k_j\tau }+\sum _{j=1}^{N_r}{C}_{-j}{\mathbf{g}}_{-j}\left(\pm {\mu }_i\right)e^{k_j\tau }}}\hfill \\ \hfill +& \sum _{j=1}^{N_i}\left[C_{1j}{\widehat{\mathbf{g}}}_{1j}\left(\tau ,\pm {\mu }_i\right)+C_{2j}{\widehat{\mathbf{g}}}_{2j}\left(\tau ,\pm {\mu }_i\right)\right]e^{-k_{rj}\tau }\hfill \\ \hfill +& \underset{{\mathbf{I}}_{\text{COMPLEX}}}{\underbrace{\sum _{j=1}^{N_i}\left[C_{-1j}{\widehat{\mathbf{g}}}_{-1j}\left(\tau ,\pm {\mu }_i\right)+C_{-2j}{\widehat{\mathbf{g}}}_{-2j}\left(\tau ,\pm {\mu }_i\right)\right]e^{k_{rj}\tau }.}}\hfill \end{array}(90)$

The new eigenvalue/vector pairs are real numbers, and the coefficients C _{± 1j}, C _{± 2j} will be determined by the top/lower boundary as well as the layer continuity conditions. The beauty of Eq. 90 is that (after the conversion) it has the same form as that of the scalar homogeneous solution, with all of the differences incorporated in the new (4 × 1) eigenvectors ${\widehat{\mathbf{g}}}_{\pm 1}(\tau ,\pm {\mu }_i)$ and ${\widehat{\mathbf{g}}}_{\pm 2}(\tau ,\pm {\mu }_i)$.

3.3.3 Reduction of the Dimension of the Algebraic Eigenvalue Problem

A reduction of the dimension of the eigenmatrix by a factor of 2 saves a significant amount of computing time. In the scalar DISORT code, this reduction is based on the symmetry p^m (μ_i, μ_j) = p^m (−μ_i, − μ_j) of the phase function as explained in detail elsewhere (Stamnes K. et al., 2017).

However, in the vector case, the phase matrix ${\mathbf{p}}_{\alpha }^m({\mu }_i,{\mu }_j)\ne {\mathbf{p}}_{\alpha }^m(-{\mu }_i,-{\mu }_j)$, implying that the reduction of dimension can not be applied directly, and that a transformation must be done first as follows:

${\mathbf{P}}_{\alpha }^m\left(-\mu ,-{\mu }^{\prime }\right)=\mathbf{D}{\mathbf{P}}_{\alpha }^m\left(+\mu ,+{\mu }^{\prime }\right)\mathbf{D}(91)$

where the matrix D = diag (1, 1, −1, −1).

This result just tells us that we need to introduce a new phase matrix $\mathbf{D}{\mathbf{P}}_{\alpha }^m({\mu }_i,{\mu }_j)\mathbf{D}$ in the VRTE to restore the special symmetry structure of the eigenmatrix. For the downward VRTE in Eq. 86, we multiply it by D on both sides, and then add D ⋅D = I on the RHS, so that Eq. 86 becomes:

$\begin{array}{ccc}\hfill -{\mu }_i\frac{d{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)}{d\tau }\mathbf{D}& =\hfill & \mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)-\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j\mathbf{D}\cdot {\tilde{\mathbf{P}}}_{\alpha }^m\left(-{\mu }_i,-{\mu }_j\right)\mathbf{D}\cdot \mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)\hfill \\ \hfill & \hfill & -\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j\mathbf{D}\cdot {\tilde{\mathbf{P}}}_{\alpha }^m\left(-{\mu }_i,+{\mu }_j\right)\mathbf{D}\cdot \mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_j\right)\hfill \\ \hfill & =\hfill & \mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)-\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,+{\mu }_j\right)\mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)\hfill \\ \hfill & \hfill & -\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,-{\mu }_j\right)\mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_j\right).\hfill \end{array}(92)$

In Eq. 85, we simply just multiply by D ⋅D = I in front of the term ${\tilde{\mathbf{I}}}_{\alpha }^m(\tau ,-{\mu }_j)$ on the RHS, so that it becomes

$\begin{array}{ccc}\hfill +{\mu }_i\frac{d{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_i\right)}{d\tau }& =\hfill & {\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_i\right)-\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,-{\mu }_j\right)\mathbf{D}\cdot \mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)\hfill \\ \hfill & \hfill & -\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,+{\mu }_j\right){\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_j\right).\hfill \end{array}(93)$

Now we may rewrite the VRTE with the proper symmetry structure, which is suitable for the reduction of dimension.

$\begin{array}{ll}\hfill +{\mu }_i\frac{d{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_i\right)}{d\tau }& ={\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_i\right)-\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,-{\mu }_j\right)\mathbf{D}\cdot \mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)\hfill \\ \hfill & -\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,+{\mu }_j\right){\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_j\right)\hfill \end{array}(94)$

$\begin{array}{ll}\hfill -{\mu }_i\frac{d\mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)}{d\tau }& =\mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)-\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}^m\left(+{\mu }_i,+{\mu }_j\right)\mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)\hfill \\ \hfill & -\frac{\varpi \left(\tau \right)}{2}\sum _{j=1}^N{w}_j{\tilde{\mathbf{P}}}^m\left(+{\mu }_i,-{\mu }_j\right)\mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_j\right).\hfill \end{array}(95)$

Equations 94 and 95 are identical with Eqs. 85 and 86 if we make the following connections:

$\begin{array}{ll}\hfill {\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_i\right)& \to {\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+{\mu }_i\right)\left(\text{\hspace{0.17em}unchanged}\right)\hfill \\ \hfill {\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)& \to \left[\mathbf{D}{\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-{\mu }_j\right)\right]\hfill \\ \hfill {\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,+{\mu }_j\right)& \to {\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,+{\mu }_j\right)\left(\text{\hspace{0.17em}unchanged}\right)\hfill \\ \hfill {\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,-{\mu }_j\right)& \to \left[{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,-{\mu }_j\right)\mathbf{D}\right].\hfill \end{array}$

To accomplish the reduction of dimension, we define eigenvectors

${\tilde{\mathbf{g}}}_{\alpha ,+}={\mathbf{g}}_{\alpha ,+}=\left(\begin{array}{c}\hfill {\mathbf{g}}_{\alpha ,+1}\hfill \\ \hfill ⋮\hfill \\ \hfill {\mathbf{g}}_{\alpha ,+N}\hfill \end{array}\right)\left(\text{\hspace{0.17em}unchanged}\right){\tilde{\mathbf{g}}}_{\alpha ,-}=\left(\begin{array}{c}\hfill {\mathbf{D}\mathbf{g}}_{\alpha ,-1}\hfill \\ \hfill ⋮\hfill \\ \hfill {\mathbf{D}\mathbf{g}}_{\alpha ,-N}\hfill \end{array}\right)$

as well as two 4N × 4N matrices as we did in the scalar case (Lin et al., 2015):

${\mathbf{E}}_{\alpha }=\left(\begin{array}{ccc}\hfill \frac{\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}_{\alpha }^m\left({\mu }_1,{\mu }_1\right)}{2{\mu }_1}-\frac{\mathit{I}}{{\mu }_1}\hfill & \hfill \cdots \hfill & \hfill \frac{\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}_{\alpha }^m\left({\mu }_1,{\mu }_N\right)}{2{\mu }_1}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill \frac{\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}_{\alpha }^m\left({\mu }_N,{\mu }_1\right)}{2{\mu }_N}\hfill & \hfill \cdots \hfill & \hfill \frac{\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}_{\alpha }^m\left({\mu }_N,{\mu }_N\right)}{2{\mu }_N}-\frac{\mathit{I}}{{\mu }_N}\hfill \end{array}\right)$

${\mathbf{F}}_{\alpha }=\left(\begin{array}{ccc}\hfill \frac{\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}_{\alpha }^m\left({\mu }_1,-{\mu }_1\right)\mathbf{D}}{2{\mu }_1}\hfill & \hfill \cdots \hfill & \hfill \frac{\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}_{\alpha }^m\left({\mu }_1,-{\mu }_N\right)\mathit{D}}{2{\mu }_1}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill \frac{\varpi \left(\tau \right)w_1{\tilde{\mathbf{P}}}_{\alpha }^m\left({\mu }_N,-{\mu }_1\right)\mathit{D}}{2{\mu }_N}\hfill & \hfill \cdots \hfill & \hfill \frac{\varpi \left(\tau \right)w_N{\tilde{\mathbf{P}}}_{\alpha }^m\left({\mu }_N,-{\mu }_N\right)\mathbf{D}}{2{\mu }_N}\hfill \end{array}\right)$

to obtain

$\left({\mathbf{E}}_{\alpha }+{\mathbf{F}}_{\alpha }\right)\left({\mathbf{E}}_{\alpha }-{\mathbf{F}}_{\alpha }\right)=k_{\alpha }^2\left({\tilde{\mathbf{g}}}_{\alpha ,+}+{\tilde{\mathbf{g}}}_{\alpha ,-}\right).(96)$

Equation 96 is similar to the one obtained in the scalar case, so we may proceed exactly as in that case (Lin et al., 2015).

3.4 Discrete-Ordinate Approximation of the Vector Source Function

In VDISORT, the source function is a 4 × 1 vector, and the discrete-ordinate approximation of the mth Fourier component of the vector source function may be written as

$\begin{array}{ccc}\hfill {\mathbf{S}}_{\alpha }^m\left(\tau ,\pm \mu \right)& =\hfill & \frac{\varpi }{2}\sum _{i=1}^N{\omega }_i{\tilde{\mathbf{P}}}_{\alpha }^m\left(-{\mu }_i,\pm \mu \right){\mathbf{I}}_{\alpha }^m\left(\tau ,-{\mu }_i\right)\hfill \\ \hfill & \hfill & +\frac{\varpi }{2}\sum _{i=1}^N{\omega }_i{\tilde{\mathbf{P}}}_{\alpha }^m\left(+{\mu }_i,\pm \mu \right){\mathbf{I}}_{\alpha }^m\left(\tau ,+{\mu }_i\right)\hfill \\ \hfill & \hfill & +{\mathbf{X}}_0^m\left(\pm \mu \right)e^{-\tau /{\mu }_0}\hfill \end{array}(97)$

$\begin{array}{ccc}\hfill {\mathbf{X}}_0^m\left(\pm \mu \right)& =\hfill & \left(2-{\delta }_{0m}\right)\frac{\varpi }{4\pi }{\tilde{\mathbf{P}}}_{\alpha }^m\left(-{\mu }_0,\pm \mu \right){\mathbf{S}}_{\mathbf{b}}.\hfill \end{array}(98)$

In Eq. 97, the Stokes vector ${\mathbf{I}}_{\alpha }^m(\tau ,-{\mu }_i)$ is given at the quadrature angles where i = 1, 2, … , 2N. For the complex eigenvalue/vector case, according to Eq. 168, the general solution can be written as (ignoring the thermal source, which can be treated in a similar manner):

$\begin{array}{ll}\hfill {\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,\pm {\mu }_i\right)=& \underset{{\mathbf{I}}_{\text{REAL}}}{\underbrace{\sum _{j=1}^{N_r}{C}_j{\mathbf{g}}_j\left(\pm {\mu }_i\right)e^{-k_j\tau }+\sum _{j=1}^{N_r}{C}_{-j}{\mathbf{g}}_{-j}\left(\pm {\mu }_i\right)e^{k_j\tau }}}+\underset{{\mathbf{I}}_{\text{PARTICULAR}}}{\underbrace{{\mathbf{Z}}_{\mathbf{0}}\left(\pm {\mu }_i\right)e^{-\tau /{\mu }_0}}}\hfill \\ \hfill +& \sum _{j=1}^{N_i}\left[C_{1j}{\widehat{\mathbf{g}}}_{1j}\left(\tau ,\pm {\mu }_i\right)+C_{2j}{\widehat{\mathbf{g}}}_{2j}\left(\tau ,\pm {\mu }_i\right)\right]e^{-k_{rj}\tau }\hfill \\ \hfill +& \underset{{\mathbf{I}}_{\text{COMPLEX}}}{\underbrace{\sum _{j=1}^{N_i}\left[C_{-1j}{\widehat{\mathbf{g}}}_{-1j}\left(\tau ,\pm {\mu }_i\right)+C_{-2j}{\widehat{\mathbf{g}}}_{-2j}\left(\tau ,\pm {\mu }_i\right)\right]e^{k_{rj}\tau }}}\hfill \end{array}(99)$

where the first two summation terms on the right hand side are the homogeneous solution for real eigenvalues/vectors, the following term is the particular solution, and the last two summation terms are the homogeneous solution for complex conjugate eigenvalues/vectors. The real vectors ${\widehat{\mathbf{g}}}_{1j}(\tau ,\pm {\mu }_i)$ and ${\widehat{\mathbf{g}}}_{2j}(\tau ,\pm {\mu }_i)$ are defined in Eqs. 163–167.

Substituting Eq. 99 into the vector source function given by Eq. 97, we may write the vector source function [Eq. 97] in a compact form similar to the general solutions [Eq. 99]:

$\begin{array}{ll}\hfill {\mathbf{S}}_{\alpha }^m\left(\tau ,\pm \mu \right)=& \sum _{j=1}^{N_r}{C}_j{\tilde{\mathbf{g}}}_j\left(\pm \mu \right)e^{-k_j\tau }+\sum _{j=1}^{N_r}{C}_{-j}{\tilde{\mathbf{g}}}_{-j}\left(\pm \mu \right)e^{k_j\tau }+{\tilde{\mathbf{Z}}}_0\left(\pm \mu \right)e^{-\tau /{\mu }_0}\hfill \\ \hfill +& \sum _{j=1}^{N_i}\left[C_{1j}{\widetilde{\widehat{\mathbf{g}}}}_{1j}\left(\tau ,\pm \mu \right)+C_{2j}{\widetilde{\widehat{\mathbf{g}}}}_{2j}\left(\tau ,\pm \mu \right)\right]e^{-k_{rj}\tau }\hfill \\ \hfill +& \sum _{j=1}^{N_i}\left[C_{-1j}{\widetilde{\widehat{\mathbf{g}}}}_{-1j}\left(\tau ,\pm \mu \right)+C_{-2j}{\widetilde{\widehat{\mathbf{g}}}}_{-2j}\left(\tau ,\pm \mu \right)\right]e^{k_{rj}\tau }.\hfill \end{array}(100)$

In Eq. 100 the following expressions.

${\tilde{\mathbf{g}}}_{\pm j}\left(\pm \mu \right)=\frac{\varpi }{2}\sum _{i=1}^{N_r}\left\{{\omega }_i{\mathbf{P}}_{\alpha }^m\left(-{\mu }_i,\pm \mu \right){\mathbf{g}}_{\pm j}\left(-{\mu }_i\right)+{\omega }_i{\mathbf{P}}_{\alpha }^m\left({\mu }_i,\pm \mu \right){\mathbf{g}}_{\pm j}\left(+{\mu }_i\right)\right\}(101)$

${\widetilde{\widehat{\mathbf{g}}}}_{\pm 1j}\left(\tau ,\pm \mu \right)=\frac{\varpi }{2}\sum _{i=1}^{N_i}\left\{{\omega }_i{\mathbf{P}}_{\alpha }^m\left(-{\mu }_i,\pm \mu \right){\widehat{\mathbf{g}}}_{\pm 1j}\left(\tau ,-{\mu }_i\right)+{\omega }_i{\mathbf{P}}_{\alpha }^m\left({\mu }_i,\pm \mu \right){\widehat{\mathbf{g}}}_{\pm 1j}\left(\tau ,{\mu }_i\right)\right\}(102)$

${\widetilde{\widehat{\mathbf{g}}}}_{\pm 2j}\left(\tau ,\pm \mu \right)=\frac{\varpi }{2}\sum _{i=1}^{N_i}\left\{{\omega }_i{\mathbf{P}}_{\alpha }^m\left(-{\mu }_i,\pm \mu \right){\widehat{\mathbf{g}}}_{\pm 2j}\left(\tau ,-{\mu }_i\right)+{\omega }_i{\mathbf{P}}_{\alpha }^m\left({\mu }_i,\pm \mu \right){\widehat{\mathbf{g}}}_{\pm 2j}\left(\tau ,{\mu }_i\right)\right\}(103)$

${\tilde{\mathbf{Z}}}_0\left(\pm \mu \right)=\frac{\varpi }{2}\sum _{i=1}^N\left\{{\omega }_i{\mathbf{P}}_{\alpha }^m\left(-{\mu }_i,\pm \mu \right){\mathbf{Z}}_0\left(-{\mu }_i\right)+{\omega }_i{\mathbf{P}}_{\alpha }^m\left({\mu }_i,\pm \mu \right){\mathbf{Z}}_0\left({\mu }_i\right)+{\mathbf{X}}_0\left(\pm \mu \right)\right\}(104)$

are simply convenient analytic interpolation formulas of g, ${\widehat{\mathbf{g}}}_1$, ${\widehat{\mathbf{g}}}_2$, and Z₀. They clearly reveal the interpolatory nature of Eq. 100 for the vector source function. The fact that they are derived from the basic VRTE to which we seek solutions indicates that these expressions, like the analogous expressions in the scalar case, will be superior to any other interpolation scheme (Stamnes et al., 1988; Schulz and Stamnes, 2000; Lin et al., 2015; Stamnes K. et al., 2017).

3.5 Integration of the Source Function Method–Solutions at User-Desired Polar Angles

The discrete ordinate solutions for the Stokes vector are computed at the quadrature angles as discussed in the previous section. To obtain values at arbitrary angles as desired by the user, an interpolation algorithm has to be implemented. In previous versions of VDISORT, a standard spline interpolation scheme, shown to work well for Rayleigh scattering (Schulz et al., 1999), was used to obtain output at arbitrary polar angles. However, the spline interpolation generally requires a large number of quadrature angles (number of streams), and it may fail when the particles have sharp forward scattering peaks, such as for large cloud droplets or ice crystals, when the Stokes components may change rapidly with polar angle.

A better approach to interpolation is to use the discrete-ordinate solution to derive explicit expressions for the source function that can be integrated analytically. This ISF method is implemented in DISORT (Stamnes et al., 1988; Lin et al., 2015; Stamnes K. et al., 2017) and also in a previous version of VDISORT (Schulz and Stamnes, 2000), but that solution is not valid for the V component of the Stokes vector because the eigenvalues/vectors were assumed to be real. Therefore, the solution must be extended to apply to the general case for which some of the eigenvalues/eigenvectors may be complex. By doing so, we obtain results at arbitrary polar angles for any given number of streams, and we may save computing time by getting accurate results at arbitrary polar angles for a relatively small number of streams. Below we will use the ISF method to derive a new interpolation algorithm that works well also for the complex eigensolutions.

3.5.1 Single-Layer (Homogeneous) Medium

For a slab of thickness τ*, we may solve Eqs. 83 and 84 to obtain.

${\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+\mu \right)={\tilde{\mathbf{I}}}_{\alpha }^m\left({\tau }^{\ast },+\mu \right)e^{-\left({\tau }^{\ast }-\tau \right)/\mu }+{\int }_{\tau }^{{\tau }^{\ast }}\frac{dt}{\mu }{\mathbf{S}}_{\alpha }^m\left(t,+\mu \right)e^{-\left(t-\tau \right)/\mu }(105)$

${\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-\mu \right)={\tilde{\mathbf{I}}}_{\alpha }^m\left(0,-\mu \right)e^{-\tau /\mu }+{\int }_0^{\tau }\frac{dt}{\mu }{\mathbf{S}}_{\alpha }^m\left(t,-\mu \right)e^{-\left(\tau -t\right)/\mu }.(106)$

Using Eq. 100 in Eqs. 105 and 106, we find that for a slab of thickness τ*, the Stokes vectors become

$\begin{array}{ll}\hfill \tilde{\mathbf{I}}\left(\tau ,+\mu \right)=& \tilde{\mathbf{I}}\left({\tau }^{*},+\mu \right)e^{-\frac{{\tau }^{*}-\tau }{\mu }}+\sum _{\begin{array}{c}j=-N_r\\ j\ne 0\end{array}}^{N_r}{C}_j\frac{{\tilde{\mathbf{g}}}_j\left(+\mu \right)}{1+k_j\mu }\left\{e^{-k_j\tau }-e^{-\left[k_j{\tau }^{*}+\left({\tau }^{*}-\tau \right)/\mu \right]}\right\}\hfill \\ \hfill & +\sum _{\begin{array}{c}j=-N_i\\ j\ne 0\end{array}}^{N_i}\left\{{\mathbf{G}}_j\left(\tau ,+\mu \right)e^{-k_{rj}\tau }-{\mathbf{G}}_j\left({\tau }^{*},+\mu \right)e^{-\left[k_{rj}{\tau }^{*}+\left({\tau }^{*}-\tau \right)/\mu \right]}\right\}\hfill \\ \hfill & +\frac{{\tilde{\mathbf{Z}}}_0\left(+\mu \right)}{1+\mu /{\mu }_0}\left\{e^{-\tau /{\mu }_0}-e^{-\left[{\tau }^{*}/{\mu }_0+\left({\tau }^{*}-\tau \right)/\mu \right]}\right\}\hfill \end{array}(107)$

$\begin{array}{ll}\hfill \tilde{\mathbf{I}}\left(\tau ,-\mu \right)=& \tilde{\mathbf{I}}\left(0,-\mu \right)e^{-\frac{\tau }{\mu }}+\sum _{\begin{array}{c}j=-N_r\\ j\ne 0\end{array}}^{N_r}{C}_j\frac{{\tilde{\mathbf{g}}}_j\left(-\mu \right)}{1-k_j\mu }\left\{e^{-k_j\tau }-e^{-\tau /\mu }\right\}\hfill \\ \hfill & +\sum _{\begin{array}{c}j=-N_i\\ j\ne 0\end{array}}^{N_i}\left\{{\mathbf{G}}_j\left(\tau ,-\mu \right)e^{-k_{rj}\tau }-{\mathbf{G}}_j\left(0,-\mu \right)e^{-\tau /\mu }\right\}\hfill \\ \hfill & +\frac{{\tilde{\mathbf{Z}}}_0\left(-\mu \right)}{1-\mu /{\mu }_0}\left\{e^{-\tau /{\mu }_0}-e^{-\tau /\mu }\right\}\hfill \end{array}(108)$

where

$\begin{array}{ccc}\hfill {\mathbf{G}}_j\left(\tau ,\pm \mu \right)& =\hfill & \frac{1}{{\left(\mu k_{ij}\right)}^2+{\left(1\pm \mu k_{rj}\right)}^2}\hfill \\ \hfill & \hfill & \times \left\{\left(1\pm \mu k_{rj}\right)\left[C_{1j}{\widetilde{\widehat{\mathbf{g}}}}_{1j}\left(\tau ,\pm \mu \right)+C_{2j}{\widetilde{\widehat{\mathbf{g}}}}_{2j}\left(\tau ,\pm \mu \right)\right]\right.\hfill \\ \hfill & \hfill & \left.\pm \mu k_{ij}\left[C_{1j}{\widetilde{\widehat{\mathbf{g}}}}_{2j}\left(\tau ,\pm \mu \right)-C_{2j}{\widetilde{\widehat{\mathbf{g}}}}_{1j}\left(\tau ,\pm \mu \right)\right]\right\}.\hfill \end{array}(109)$

Here we assumed k_−j = −k_j, k_−rj = −k_rj and k_−ij = −k_ij and thus k_−cj = −k_cj = −k_rj − ik_ij.

In Eqs. 107 and 108, the first term on the RHS is due to beam attenuation of the Stokes vector. The first of the following summation terms stems from the integration of the homogeneous solution with real eigensolutions in the vector source function [Eq. 100], whereas the second summation is a new term needed for the complex eigensolutions in the vector source function [Eq. 100]. Because ${\widetilde{\widehat{\mathbf{g}}}}_{1j}(\tau ,+\mu )$ and ${\widetilde{\widehat{\mathbf{g}}}}_{2j}(\tau ,+\mu )$ are functions of ${\widehat{g}}_{1j}(\tau ,+\mu )$ and ${\widehat{g}}_{2j}(\tau ,+\mu )$ that depend on cos (k_iτ) and sin (k_iτ), we used the following equations to handle the integration of the complex eigenvector/vectors:

$\int \cos ax{e}^{bx}dx=\frac{1}{a^2+b^2}{e}^{bx}\left(a\sin ax+b\cos ax\right)(110)$

$\int \sin ax{e}^{bx}dx=\frac{1}{a^2+b^2}{e}^{bx}\left(b\sin ax-a\cos ax\right).(111)$

3.5.2 Multi-Layer (Inhomogenous) Medium

The single-layer case can be extended into a multi-layer medium, for which we need to evaluate the integral by integrating layer-by-layer as follows.

$\begin{array}{ll}\hfill {\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+\mu \right)=& {\tilde{\mathbf{I}}}_{\alpha }^m\left({\tau }_L,+\mu \right)e^{-\left({\tau }_L-\tau \right)/\mu }+{\int }_{\tau }^{{\tau }_p}\frac{dt}{\mu }{\mathbf{S}}_{\alpha ,p}^m\left(t,+\mu \right)e^{-\left(t-\tau \right)/\mu }\hfill \\ \hfill & +\sum _{n=p+1}^L{\int }_{{\tau }_{n-1}}^{{\tau }_n}\frac{dt}{\mu }{\mathbf{S}}_{\alpha ,n}^m\left(t,+\mu \right)e^{-\left(t-\tau \right)/\mu }\hfill \end{array}(112)$

$\begin{array}{ll}\hfill {\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-\mu \right)=& {\tilde{\mathbf{I}}}_{\alpha }^m\left(0,-\mu \right)e^{-\tau /\mu }+{\int }_{{\tau }_{p-1}}^{\tau }\frac{dt}{\mu }{\mathbf{S}}_{\alpha ,p}^m\left(t,-\mu \right)e^{-\left(\tau -t\right)/\mu }\hfill \\ \hfill & +\sum _{n=1}^{p-1}{\int }_{{\tau }_{n-1}}^{{\tau }_n}\frac{dt}{\mu }{\mathbf{S}}_{\alpha ,n}^m\left(t,-\mu \right)e^{-\left(\tau -t\right)/\mu }\hfill \end{array}(113)$

Using Eq. 100 for ${\mathbf{S}}_{\alpha ,n}^m(t,-\mu )$ in the nth layer, Eqs. 112 and 113 become:

$\begin{array}{ll}\hfill & {\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,+\mu \right)={\tilde{\mathbf{I}}}_{\alpha }^m\left({\tau }_L,+\mu \right)e^{-\frac{{\tau }_L-\tau }{\mu }}\hfill \\ \hfill & +\sum _{n=p}^L\sum _{\begin{array}{c}j=-N\\ j\ne 0\end{array}}^N{C}_{jn}\frac{{\tilde{\mathbf{g}}}_{jn}\left(+\mu \right)}{1+k_{jn}\mu }\left\{e^{-k_{jn}{\tau }_{n-1}+\left({\tau }_{n-1}-\tau \right)/\mu }-e^{-\left[k_{jn}{\tau }_n+\left({\tau }_n-\tau \right)/\mu \right]}\right\}\hfill \\ \hfill & +\sum _{n=p}^L\sum _{\begin{array}{c}j=-N\\ j\ne 0\end{array}}^N\left\{{\mathbf{G}}_{jn}\left({\tau }_{n-1},+\mu \right)e^{-k_{\mathit{rjn}}{\tau }_{n-1}+\left({\tau }_{n-1}-\tau \right)/\mu }-{\mathbf{G}}_{jn}\left({\tau }_n,+\mu \right)e^{-\left[k_{jn}{\tau }_n+\left({\tau }_n-\tau \right)/\mu \right]}\right\}\hfill \\ \hfill & +\sum _{n=p}^L\frac{{\tilde{\mathbf{Z}}}_{0n}\left(+\mu \right)}{1+\mu /{\mu }_0}\left\{e^{-\left[{\tau }_{n-1}/{\mu }_0+\left({\tau }_{n-1}-\tau \right)/\mu \right]}-e^{-\left[{\tau }_n/{\mu }_0+\left({\tau }_n-\tau \right)/\mu \right]}\right\}\hfill \end{array}(114)$

with τ_n−1 replaced by τ for n = p, and

$\begin{array}{ll}\hfill & {\tilde{\mathbf{I}}}_{\alpha }^m\left(\tau ,-\mu \right)={\tilde{\mathbf{I}}}_{\alpha }^m\left(0,-\mu \right)e^{-\frac{\tau }{\mu }}\hfill \\ \hfill & +\sum _{n=1}^p\sum _{\begin{array}{c}j=-N\\ j\ne 0\end{array}}^N{C}_{jn}\frac{{\tilde{\mathbf{g}}}_{jn}\left(-\mu \right)}{1-k_{jn}\mu }\left\{e^{-\left[k_{jn}{\tau }_n+\left(\tau -{\tau }_n\right)/\mu \right]}-e^{-\left[k_{jn}{\tau }_{n-1}+\left(\tau -{\tau }_{n-1}\right)/\mu \right]}\right\}\hfill \\ \hfill +& \sum _{n=1}^p\sum _{\begin{array}{c}j=-N\\ j\ne 0\end{array}}^N\left\{{\mathbf{G}}_{jn}\left({\tau }_n,-\mu \right)e^{-\left[k_{\mathit{rjn}}{\tau }_n+\left(\tau -{\tau }_n\right)/\mu \right]}-{\mathbf{G}}_{jn}\left({\tau }_{n-1},-\mu \right)e^{-\left[k_{\mathit{rjn}}{\tau }_{n-1}+\left(\tau -{\tau }_{n-1}\right)/\mu \right]}\right\}\hfill \\ \hfill +& \sum _{n=1}^p\frac{{\tilde{\mathbf{Z}}}_{0n}\left(-\mu \right)}{1-\mu /{\mu }_0}\left\{e^{-\left[{\tau }_n/{\mu }_0+\left(\tau -{\tau }_n\right)/\mu \right]}-e^{-\left[{\tau }_{n-1}/m{u}_0+\left(\tau -{\tau }_{n-1}\right)/\mu \right]}\right\}\hfill \end{array}(115)$

with τ_n replaced by τ for n = p, and where

$\begin{array}{ccc}\hfill {\mathbf{G}}_{jn}\left(\tau ,\pm \mu \right)& =\hfill & \frac{1}{{\left(\mu k_{\mathit{ijn}}\right)}^2+{\left(1\pm \mu k_{\mathit{rjn}}\right)}^2}\hfill \\ \hfill & \hfill & \times \left\{\left(1\pm \mu k_{\mathit{rjn}}\right)\left[C_{\mathit{1jn}}{\widetilde{\widehat{\mathbf{g}}}}_{\mathit{1jn}}\left(\tau ,\pm \mu \right)+C_{\mathit{2jn}}{\widetilde{\widehat{\mathbf{g}}}}_{\mathit{2jn}}\left(\tau ,\pm \mu \right)\right]\right.\hfill \\ \hfill & \hfill & \left.\pm \mu k_{\mathit{ijn}}\left[C_{\mathit{1jn}}{\widetilde{\widehat{\mathbf{g}}}}_{\mathit{2jn}}\left(\tau ,\pm \mu \right)-C_{\mathit{2jn}}{\widetilde{\widehat{\mathbf{g}}}}_{\mathit{1jn}}\left(\tau ,\pm \mu \right)\right]\right\}.\hfill \end{array}(116)$

It can be verified that for a single layer τ_n−1 = τ, τ_n = τ_L = τ* in Eq. 114; τ_n = τ, τ_n−1 = 0 in Eq. 115, they are reduced to Eqs. 107 and 108, as they should.

3.5.3 Numerical Example

Many practical problems including remote sensing applications require Stokes vector components at sensor observing angles. Even though one could envision generating results at sensor observing angles by interpolating or extrapolating the quadrature values to such angles, accurate interpolation of these values is difficult, particularly for optical depths close to zero (upper boundary) and for extrapolation to angles close to μ = 0 and μ = 1.0. For example, Figure 2 shows an example of inaccurate spline interpolation for a benchmark case (Garcia and Siewert, 1989) that has been reproduced by VDISORT.We note that the spline interpolation produces large oscillations, whereas the ISF method, see Eqs. 105 and 106 yields accurate analytic results that agree with the benchmark values at the quadrature angles.

FIGURE 2

FIGURE 2. Comparison of output for (i) quadrature angles (μ = cos θ, where θ is the polar angle); (ii) spline interpolation; and (iii) analytic results at arbitrary output angles.

We emphasize that the ISF method allows one to compute analytically the radiation field at arbitrary angles and optical depths from the discrete ordinate solutions, both for the scalar radiative transfer problem (Lin et al., 2015; Stamnes K. et al., 2017) and, as shown here, for the vector problem including the complete Stokes vector I_S = [I, Q, U, V]^T. This capability is one of the unique advantages of the discrete ordinate method. Also, the ISF solutions satisfy the boundary and continuity conditions not only at the discrete set of quadrature angles, but at arbitrary angles μ. In fact, we have found that the ISF method improves the accuracy of the discrete ordinate solution (Schulz and Stamnes, 2000).

3.6 Polarized Reflectance at the Lower Boundary

In Eq. 99 the homogeneous solution contains unknown constants C_±j, C _{± 1j}, and C _{± 2j}, which are to be determined by the boundary conditions at the top and bottom of the medium and continuity conditions at layer interfaces.

3.6.1 Top Boundary

In VDISORT, we assume that the light at the top of the atmosphere is a direct beam such that the diffuse light contribution is zero:

${\tilde{\mathbf{I}}}_{\alpha }^m\left(0,-\mu \right)={\left[\mathrm{0,0,0,0}\right]}^T(117)$

for m = 0, … , 2N − 1.

3.6.2 Layer Interfaces

At layer interfaces the Stokes vector must satisfy continuity conditions, because by assumption there is no change in the refractive index between layers. Assuming that τ_n,bottom is the optical depth at the bottom of the nth layer, and that τ_n+1,top is the optical depth at the top of the (n + 1)^th layer, we have.

${\tau }_{n,\text{bottom}}={\tau }_{n+1,\text{top}}(118)$

${\tilde{\mathbf{I}}}_{\alpha }^m\left({\tau }_{n,\text{bottom}},\pm \mu \right)={\tilde{\mathbf{I}}}_{\alpha }^m\left({\tau }_{n+1,\text{top}},\pm \mu \right)(119)$

where the Stokes vector ${\tilde{\mathbf{I}}}_{\alpha }(\tau ,\pm \mu )$ is given by the sum of homogeneous and particular solutions in Eq. 99.The continuity conditions are applied layer by layer for n = 1, … , L − 1. One important difference from the scalar DISORT model is that the new real eigenvectors associated with the complex solutions depend on the optical depth τ in each layer [Eq. 90]. Therefore, for jth eigenvector in the nth layer, we have.

${\widehat{\mathbf{g}}}_{1j}\left({\tau }_{n,\text{top}},\pm {\mu }_i\right)\ne {\widehat{\mathbf{g}}}_{1j}\left({\tau }_{n,\text{bottom}},\pm {\mu }_i\right)(120)$

${\widehat{\mathbf{g}}}_{2j}\left({\tau }_{n,\text{top}},\pm {\mu }_i\right)\ne {\widehat{\mathbf{g}}}_{2j}\left({\tau }_{n,\text{bottom}},\pm {\mu }_i\right).(121)$

3.6.3 Lower Boundary

For the RTE, the lower boundary is determined by the bidirectional reflectance distribution function/matrix (BRDF). We introduce the matrix or polarized BRDF as the 4 × 4 reflection matrix R by writing:

$\begin{array}{ccc}\hfill \mathbf{I}\left(\mu ,\varphi ,{\tau }_L\right)& =\hfill & {\int }_0^{2\pi }d{\varphi }^{\prime }{\int }_0^{1}d{\mu }^{\prime }\mathbf{R}\left(-{\mu }^{\prime },\mu ,\varphi -{\varphi }^{\prime }\right)\mathbf{I}\left(-{\mu }^{\prime },{\varphi }^{\prime },{\tau }_L\right)\hfill \\ \hfill & \hfill & +\mathbf{R}\left(-{\mu }_0,\mu ,\varphi -{\varphi }_0\right){\mathbf{S}}_{\mathbf{b}}{e}^{-{\tau }_L/{\mu }_0}.\hfill \end{array}(122)$

The reflection matrix R (−μ′, μ, ϕ − ϕ′) is then expanded in a Fourier series:

$\begin{array}{ccc}\hfill \mathbf{R}\left(-{\mu }^{\prime },\mu ,\varphi -{\varphi }^{\prime }\right)& =\hfill & \sum _{m=1}^L{\mathbf{R}}_c^m\left(-{\mu }^{\prime },\mu \right)\cos m\left(\varphi -{\varphi }^{\prime }\right)\hfill \\ \hfill & \hfill & +\sum _{m=1}^L{\mathbf{R}}_s^m\left(-{\mu }^{\prime },\mu \right)\sin m\left(\varphi -{\varphi }^{\prime }\right).\hfill \end{array}(123)$

The Fourier mode ${\tilde{\mathbf{I}}}_{\alpha }^m({\tau }_L,\mu )$ of the Stokes vector at the polar angle μ just above the lower boundary can then be expressed in terms of $\tilde{\mathbf{R}}(\mu ,{\mu }^{\prime })$ by (detailed derivation omitted for brevity):

${\tilde{\mathbf{I}}}_{\alpha }^m\left({\tau }_L,\mu \right)={\tilde{\mathbf{R}}}_{\alpha ,\text{beam}}^m\left(-{\mu }_0,\mu \right){\mathbf{S}}_{\mathbf{b}}{e}^{-{\tau }_L/{\mu }_0}+\pi \sum _{j=1}^N{\omega }_j{\tilde{\mathbf{R}}}_{\alpha }^m\left(-{\mu }_j,\mu \right){\tilde{\mathbf{I}}}_{\alpha }^m\left({\tau }_L,-{\mu }_j\right)(124)$

where.

${\tilde{\mathbf{R}}}_{c,\text{beam}}^m\left(-{\mu }_0,\mu \right)=\left(\begin{array}{cccc}\hfill R_{\mathit{c11}}^m\hfill & \hfill R_{\mathit{c12}}^m\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill R_{\mathit{c21}}^m\hfill & \hfill R_{\mathit{c22}}^m\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill R_{\mathit{s31}}^m\hfill & \hfill R_{\mathit{s32}}^m\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill R_{\mathit{s41}}^m\hfill & \hfill R_{\mathit{s42}}^m\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right)(125)$

${\tilde{\mathbf{R}}}_{s,\text{beam}}^m\left(-{\mu }_0,\mu \right)=\left(\begin{array}{cccc}\hfill 0\hfill & \hfill 0\hfill & \hfill R_{\mathit{s11}}^m\hfill & \hfill R_{\mathit{s12}}^m\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill R_{\mathit{s21}}^m\hfill & \hfill R_{\mathit{s22}}^m\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill R_{\mathit{c31}}^m\hfill & \hfill R_{\mathit{c32}}^m\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill R_{\mathit{c41}}^m\hfill & \hfill R_{\mathit{c42}}^m\hfill \end{array}\right)(126)$

${\tilde{\mathbf{R}}}_c^m\left(-{\mu }_j,\mu \right)=\left(\begin{array}{cccc}\hfill {\mathrm{\Delta }}_{0m}^{+}{R}_{\mathit{c11}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{+}{R}_{\mathit{c12}}\hfill & \hfill -{\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s13}}\hfill & \hfill -{\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s14}}\hfill \\ \hfill {\mathrm{\Delta }}_{0m}^{+}{R}_{\mathit{c21}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{+}{R}_{\mathit{c22}}\hfill & \hfill -{\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s23}}\hfill & \hfill -{\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s24}}\hfill \\ \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s31}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s32}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{c33}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{c34}}\hfill \\ \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s41}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s42}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{c43}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{c44}}\hfill \end{array}\right)(127)$

${\tilde{\mathbf{R}}}_s^m\left(-{\mu }_j,\mu \right)=\left(\begin{array}{cccc}\hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{c11}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{c12}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s13}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s14}}\hfill \\ \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{c21}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{c22}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s23}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s24}}\hfill \\ \hfill -{\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s31}}\hfill & \hfill -{\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s32}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{+}{R}_{\mathit{c33}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{+}{R}_{\mathit{c34}}\hfill \\ \hfill -{\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s41}}\hfill & \hfill -{\mathrm{\Delta }}_{0m}^{-}{R}_{\mathit{s42}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{+}{R}_{\mathit{c43}}\hfill & \hfill {\mathrm{\Delta }}_{0m}^{+}{R}_{\mathit{c44}}\hfill \end{array}\right)(128)$

where ${\mathrm{\Delta }}_{0m}^{\pm }\equiv 1\pm {\delta }_{0m}$.

In VDISORT, the implementation of the polarized BRDF has been tested by assuming a vacuum layer above a rough ocean surface (Tsang et al., 1985) and by matching the values of the Stokes vector just above the surface (in the upward direction) with analytic values of the BRDF computed without the Fourier expansion. Figure 3 shows an example of a comparison between VDISORT upward reflected Stokes parameters and analytic Cox-Munk BRDF results for a 1-D Gaussian surface with a wind speed of 5 m/s. The V component vanishes because the reflection from an unpolarized direct beam does not produce circular polarization. The results show that VDISORT reproduced the analytic results except for backscattering (ϕ = 180°) of the U component where the value of zero is given as numerical noise $(\sim 10^{-10})$ by VDISORT.

FIGURE 3

FIGURE 3. Upper panel: The upward reflection of Stokes parameter I (μ, τ_L, ϕ) just above the lower boundary for a 1-D Gaussian surface with a ‘Cox-Munk’ slope variance distribution for a wind speed of 5 m/s. Middle panel: Same as upper panel except for the Stokes parameter Q (μ, τ_L, ϕ). Lower panel: Same as upper panel except for the Stokes parameter U (μ, τ_L, ϕ).

4 The 4 × 4 Solution Versus the 3 × 3 Approximation

In the discrete ordinate method of radiative transfer, we need to determine homogenous and particular solutions to arrive at the general solution. The particular solution is formulated as a set of linear equations Ax = b that can quickly be solved using standard techniques of linear algebra, for example Gaussian elimination. The most time-consuming step is solving the homogenous problem, which is formulated as a standard algebraic eigenvalue problem (A− λ)x = 0. In the case of the 4 × 4 solution, this eigenvalue problem involves matrices of size 4N × 4N, where N is the number of quadrature points or “streams” in the upper and lower hemispheres, since a reduction of dimension step (Section 3.3.3) is used to reduce the matrix dimension from 8N × 8N–4N × 4N. This step is completely analogous to that used in DISORT (Stamnes et al., 1988) to reduce a matrix of dimension 2N × 2N to N × N. Since the presence of the sign on b₂ in Eq. 20 leads to a matrix A in the algebraic eigenvalue problem that cannot be made symmetric, the eigenvalues and corresponding eigenvectors in the 4 × 4 representation occur in complex conjugate pairs. In the 3 × 3 approximation, the 4th row and column in Eq. 20 are simply omitted, and the resulting impact on the scattering phase matrix is obtained by setting b₂ = 0 in Eq. 24 which results in Eq. 41. Then the matrix A is symmetric implying that the resulting eigenvalues and eigenvectors are real. Note that setting b₂ = 0 in the 4 × 4 representation leads to I, Q, and U parameters identical to those obtained in the 3 × 3 approximation.

Since setting b₂ = 0 decouples the V component from I, Q, and U, the eigenvalue problem required to solve the homogeneous system in N discrete ordinates is reduced from solving a 4N × 4N system to solving a 3N × 3N system. The computational burden of solving an eigenvalue problem scales like n³ where n is the dimension of the matrix. Therefore, the 3 × 3 approximation reduces the computational burden by a theoretical factor of 4³/3³ ≈ 2.37. Since the resulting eigenvectors and eigenvalues for the 3N × 3N system are real, significant further computational savings are obtained by using an eigensolver, such as ASYMTX available in the DISORT package (Stamnes et al., 1988) that avoids unnecessary complex arithmetic. In the 4 × 4 case some of the eigenvectors and eigenvalues occur in complex conjugate pairs implying that complex arithmetic must be considered to obtain accurate solutions when b₂ ≠ 0.

A comparison of results produced by VDISORT and by a doubling-adding method (VLBLE) is provided in Figure 4 (Stamnes S. et al., 2017), where the reflected components are plotted against the polar angle θ_R, where θ_R = 0° is the zenith direction, and θ_R = 90° is the horizon. The transmitted components are plotted against the polar angle θ_T, where θ_T = 0° is the nadir direction, and θ_T = 90° is the horizon. We note that the results for I and Q produced by the 3 × 3 approximation are essentially identical to the more computationally demanding 4 × 4 results.

FIGURE 4

FIGURE 4. Top panels: reflected I and Q components at polar angles θ_R; bottom panels: transmitted I and Q components at polar angles θ_T. The incident Stokes vector is I_S =[1,0,0,0]. The actual reflection and transmission for I and Q are scaled by a factor $\frac{\pi }{{\mu }_0{F}_0}$ where F₀=1 and μ₀ is the cosine of the solar zenith angle θ₀. NSTR = 120. The input parameters (phase matrix) for this ice (cirrus) cloud consisting on non-spherical ice crystals are described elsewhere (Stamnes S. et al., 2017).

The results shown in Figure 4 pertain to a cirrus cloud consisting of non-spherical ice crystals. The adequacy of the 3 × 3 approximation was investigated for spherical particles by Hansen (1971) who concluded that it is usually adequate to work with 3 × 3 matrices to compute multiple-scattering polarization properties. Hansen (1971) investigated errors only of the reflected radiance and the degree of polarization; errors for the individual Stokes components Q and U, and transmittances were not considered. The results shown in Figure 4 suggest that the 3 × 3 approximation holds not only for water clouds, but also for non-spherical ice crystals, and that it applies not just to the reflected radiation, but also to the transmitted radiation and to the Stokes parameters Q and U (not shown) (Stamnes S. et al., 2017).

5 Single-Scattering Solution

Introducing the half-range Stokes vectors (the ^± signs denote the upper (⁺) and the lower (^-) hemispheres, respectively).

${\mathbf{I}}^{+}\left(\tau ,\theta ,\varphi \right)\equiv {\mathbf{I}}^{+}\left(\tau ,\theta \le \pi /2,\varphi \right)={\mathbf{I}}^{+}\left(\tau ,\mu ,\varphi \right)(129)$

${\mathbf{I}}^{-}\left(\tau ,\theta ,\varphi \right)\equiv {\mathbf{I}}^{-}\left(\tau ,\theta >\pi /2,\varphi \right)={\mathbf{I}}^{-}\left(\tau ,\mu ,\varphi \right)(130)$

and similar definitions for S^±(τ, μ, ϕ), Eq. 7 may be re-written as (μ ≡|u|).

$\mu \frac{d{\mathbf{I}}^{+}\left(\tau ,\mu ,\varphi \right)}{d\tau }={\mathbf{I}}^{+}\left(\tau ,\mu ,\varphi \right)-{\mathbf{S}}^{+}\left(\tau ,\mu ,\varphi \right)(131)$

$-\mu \frac{d{\mathbf{I}}^{-}\left(\tau ,\mu ,\varphi \right)}{d\tau }={\mathbf{I}}^{-}\left(\tau ,\mu ,\varphi \right)-{\mathbf{S}}^{-}\left(\tau ,\mu ,\varphi \right).(132)$

5.1 Single-Layer (Homogeneous) Medium

For a homogeneous slab, we adopt the following notation:

1) τ₀ is the optical depth at the upper boundary (top of the slab);

2) τ_b is the optical depth at the lower boundary (bottom of the slab).

Integrating Eq. 131 from τ_b to τ and Eq. 132 from τ₀ to τ (τ₀ ≤ τ ≤ τ_b), and solving for I^±(τ, μ, ϕ) we obtain.

${\mathbf{I}}^{+}\left(\tau ,\mu ,\varphi \right)={\mathbf{I}}^{+}\left({\tau }_b,\mu ,\varphi \right)e^{-\left({\tau }_b-\tau \right)/\mu }+{\int }_{\tau }^{{\tau }_b}\frac{dt}{\mu }{\mathbf{S}}^{+}\left(t,\mu ,\varphi \right)e^{-\frac{\left(t-\tau \right)}{\mu }}(133)$

${\mathbf{I}}^{-}\left(\tau ,\mu ,\varphi \right)={\mathbf{I}}^{-}\left({\tau }_0,\mu ,\varphi \right)e^{-\left(\tau -{\tau }_0\right)/\mu }+{\int }_{{\tau }_0}^{\tau }\frac{dt}{\mu }{\mathbf{S}}^{-}\left(t,\mu ,\varphi \right)e^{-\frac{\left(\tau -t\right)}{\mu }}.(134)$

Equations 133 and 134 show that if the source functions S^±(t, μ, ϕ) are known, we can obtain a solution to the radiative transfer problem by integration (numerically or analytically).

In the single-scattering approximation, we omit the multiple-scattering term in Eq. 8, so that the source function S^±(τ, μ, ϕ) in Eq. 9 simply becomes:

${\mathbf{S}}^{\pm }\left(\tau ,\mu ,\varphi \right)\equiv {\mathbf{Q}}^{\pm }\left(\tau ,\mu ,\varphi \right)=\frac{\varpi \left(\tau \right)}{4\pi }\mathbf{P}\left(\tau ,-{\mu }_0,{\varphi }_0;u,\varphi \right){\mathbf{S}}_b{e}^{-\tau /{\mu }_0}+\left[1-\varpi \left(\tau \right)\right]{\mathbf{S}}_t\left(\tau \right).(135)$

The first term on the RHS of Eq. 135 is proportional to the incident beam S_b which for an unpolarized incident beam is given by Eq. 10, while the second term is due to thermal emission, which is unpolarized, and S_t(τ) and is given by Eq. 11.

5.2 Multi-Layer (Inhomogeneous) Medium

The vertical variation of the inherent optical properties (IOPs) in a slab may be dealt with by dividing it into a number of adjacent, horizontal layers in which the IOPs are taken to be constant within each layer, but allowed to vary from layer to layer. The number of layers should be large enough to resolve the vertical variation in the IOPs. In such a multi-layered medium, consisting of a total of L layers, we may evaluate the integrals in Eqs. 133 and 134 by integrating layer by layer as follows (τ_p−1 ≤ τ ≤ τ_p and μ > 0, τ_b = τ_L, τ₀ = 0) ignoring the boundary terms (setting I⁺(τ_b, μ, ϕ) = 0 and I⁻(τ₀, μ, ϕ) = 0):

${\mathbf{I}}^{+}\left(\tau ,\mu ,\varphi \right)={\int }_{\tau }^{{\tau }_p}\frac{dt}{\mu }{\mathbf{S}}_p^{+}\left(t,\mu ,\varphi \right)e^{-\left(t-\tau \right)/\mu }+\sum _{n=p+1}^L{\int }_{{\tau }_{n-1}}^{{\tau }_n}\frac{dt}{\mu }{\mathbf{S}}_n^{+}\left(t,\mu ,\varphi \right)e^{-\left(t-\tau \right)/\mu }(136)$

${\mathbf{I}}^{-}\left(\tau ,\mu ,\varphi \right)=\sum _{n=1}^{p-1}{\int }_{{\tau }_{n-1}}^{{\tau }_n}\frac{dt}{\mu }{\mathbf{S}}_n^{-}\left(t,\mu ,\varphi \right)e^{-\left(\tau -t\right)/\mu }+{\int }_{{\tau }_{p-1}}^{\tau }\frac{dt}{\mu }{\mathbf{S}}_p^{-}\left(t,\mu ,\varphi \right)e^{-\left(\tau -t\right)/\mu }.(137)$

We can evaluate the integrals in Eqs. 136 and 137 either numerically or analytically if the source function ${\mathbf{S}}_i^{\pm }(t,\mu ,\varphi )$ in a layer denoted by subscript i = n, or p is known. Explicit solutions obtained in the single-scattering approximation are provided in Section 9 (Supplementary Appendix A2).

5.3 Alternative Single Scattering Solution

Another way to understand the one-layer single-scattering solution is to consider the output I^±(τ, μ, ϕ) as coming from two sources: 1) the attenuated incident radiation from the layer boundary and in the same direction (μ^±, ϕ), and 2) the source term contribution from the direct beam scattering. As seen in Section 5.1 [Eqs. 133 and 134], the attenuated boundary contribution is the first term on the right, while the source term is the second term on the right.

Next, we consider the single-scattering solution at layer boundaries τ_n with n ∈ [0, 1, 2, ⋯L] for a multi-layer medium. Since there is no diffuse radiation at the following two boundaries: 1) TOA (Top-Of-Atmosphere) downward τ₀; 2) BOA (Bottom-Of-Atmosphere) upward τ_L, we have1:

${\mathbf{I}}^{-}\left({\tau }_0,\mu ,\varphi \right)=0(138)$

${\mathbf{I}}^{+}\left({\tau }_L,\mu ,\varphi \right)=0.(139)$

Once we have specified the boundary conditions, the one-layer solution can be called consecutively to create the layer boundary radiation of the next layer. Hence, for n = 1, 2, … , L, we may recursively compute the downward radiation as:

${\mathbf{I}}^{-}\left({\tau }_n,\mu ,\varphi \right)={\mathbf{I}}^{-}\left({\tau }_{n-1},\mu ,\varphi \right)e^{-\left({\tau }_n-{\tau }_{n-1}\right)/\mu }+{\int }_{{\tau }_{n-1}}^{{\tau }_n}\frac{dt}{\mu }{\mathbf{S}}_n^{-}\left(t,\mu ,\varphi \right)e^{-\frac{{\tau }_n-t}{\mu }}.(140)$

Similarly, for n = L, … , 2, 1, we may recursively compute the upward radiation as:

${\mathbf{I}}^{+}\left({\tau }_{n-1},\mu ,\varphi \right)={\mathbf{I}}^{+}\left({\tau }_n,\mu ,\varphi \right)e^{-\left({\tau }_n-{\tau }_{n-1}\right)/\mu }+{\int }_{{\tau }_n}^{{\tau }_{n-1}}\frac{dt}{\mu }{\mathbf{S}}_n^{+}\left(t,\mu ,\varphi \right)e^{-\frac{t-{\tau }_{n-1}}{\mu }}.(141)$

Having obtained all layer boundary contributions in this manner, we may simply apply the single-layer solution again to get the I^±(τ, μ, ϕ) at arbitrary optical depth τ.

5.4 Single Scattering Correction

The single-scattering correction is a post-processing step that further improves the accuracy of radiance output by correcting the single-scattering term. It was developed by Nakajima and Tanaka Nakajima and Tanaka (1988) and can be used together with several phase matrix truncation methods (Wiscombe, 1977; Hu et al., 2000; Lin et al., 2018).

In many radiative transfer models, the phase function/matrix is expanded in generalized spherical functions (see Section 2.4.3). However, for strongly forward-peaked scattering a large number of phase element expansion coefficients are needed for accurate representation of the phase function/matrix. Due to the computational burden incurred by use of such a large number expansion coefficients, a truncation method is commonly applied so that only the first 2N elements are used, where 2N is set to be equal to the number of streams. This method greatly improves the computational efficiency, but also introduces radiance errors due to the approximate phase matrix representation. Since the truncation reduces the scattering cross-section and thereby enhances the direct beam contribution, the differential optical depth dτ and the single-scattering albedo ϖ are both being scaled as follows.

$\widehat{d\tau }=\left(1-f\varpi \right)d\tau (142)$

$\widehat{\varpi }=\frac{\left(1-f\right)\varpi }{1-f\varpi }(143)$

where the factor f depends on the particular truncation method used (Wiscombe, 1977; Hu et al., 2000; Lin et al., 2018).

The single-scattering correction method (Nakajima and Tanaka, 1988) was designed to decrease the error incurred by the truncation. To this end we replace the approximate single-scattering solution obtained by use of the truncated phase matrix by the correct single-scattering solution obtained from the accurate phase matrix as described in Section 5.3. Denoting P and P* as the original and truncated phase matrix, and I* as the singly-scattered radiance, we may write the single-scattering correction algorithm as follows (Nakajima and Tanaka, 1988):

${\mathbf{I}}_{\text{corrected}}\left(\widehat{\tau },\widehat{\varpi },P^{*}\right)=\mathbf{I}\left(\widehat{\tau },\widehat{\varpi },P^{*}\right)-{\mathbf{I}}^{*}\left(\widehat{\tau },\widehat{\varpi },P^{*}\right)+{\mathbf{I}}^{*}\left(\widehat{\tau },\varpi /\left(1-f\varpi \right),P\right).(144)$

On the right hand side of Eq. 144, the first term $\mathbf{I}(\widehat{\tau },\widehat{\varpi },P^{\ast })$ is the uncorrected radiance computed with the truncated phase matrix P*, the scaled optical depth $\widehat{\tau }$ and the scaled single-scattering albedo $\widehat{\varpi }$. The second (subtracted) term ${\mathbf{I}}^{\ast }(\widehat{\tau },\widehat{\varpi },P^{\ast })$ is the uncorrected singly-scattered radiance obtained by use of the truncated phase matrix P^∗, the scaled single-scattering albedo $\widehat{\varpi }$ and the scaled optical depth $\widehat{\tau }$ (as in the first term). The (added) third term ${\mathbf{I}}^{\ast }(\widehat{\tau },\varpi /(1-f\varpi ),P)$ is the accurate singly-scattered radiance obtained by use of the accurate phase matrix P, the scaled $\widehat{\tau }$, and another scaled single-scattering albedo ϖ/(1 − fϖ).

6 Discrete Ordinate Radiative Transfer Test Results

Schulz et al. (1999) tested a previous of VDISORT against benchmark results provided by Garcia and Siewert (1989). The result of the first three Stokes components were reproduced, but the V component was not considered. In our new version, this issue is fixed by implementing the complex eigensolutions as already discussed. The current version of VDISORT has been tested against benchmark results provided by (Garcia and Siewert, 1989) and Siewert (2000), and excellent agreement was found (Lin, 2016). Here we provide comparisons with benchmark results provided Kokhanovsky et al. (2010) for more challenging phase matrices that require more than 100 terms in the phase matrix expansion.

6.1 Comparison With Benchmark Results

The Kokhanovsky et al. (2010) benchmark results were provided for the Stokes parameters of both reflected and transmitted light in the case of molecular, aerosol, and cloudy multiple-scattering media at the wavelength λ = 412 nm. A black underlying surface for three values of the relative azimuth angles ϕ − ϕ₀ = 0, 90, 180° were considered and the solar zenith angle was set to 60°. The optical thickness was set to 0.3262 for all three single layer cases. Since the Rayleigh test is simple and has been well tested in previous versions of VDISORT, we will focus on the more challenging aerosol and cloud cases.

The phase matrix elements for all three cases are shown in Figure 5. Because aerosol and cloud particles are much larger than molecules, their phase matrices were calculated using Mie theory (Kokhanovsky et al., 2010). In contrast to Rayleigh scattering, aerosol and cloud particles both have a strong forward-scattering peak. There are also two peaks around 137° in scattering angle that correspond to the primary and secondary rainbows.

FIGURE 5

FIGURE 5. (A) Normalized phase function a₁; (B) − b₁/a₁; (C) a₄/a₁; (D) b₂/a₁ (Kokhanovsky et al., 2010).

The sharp forward-peaked scattering of aerosol and cloud particles shown in Figure 5 implies that a large number of terms in the Fourier expansion is required for accurate representation of the phase matrix. In fact, Kokhanovsky et al. (2010) provide about 1,000 terms of Greek constants, and below we use the first 148 terms for both the aerosol and cloud cases in VDISORT in an attempt to match the benchmark results.

6.1.1 Aerosol Case

The normalized log-normal density distribution of aerosol particles considered in the benchmark computations of Kokhanovsky et al. (2010) had mode radius r_g = 0.3 μm corresponding to a mode size parameter of $\sim 2.3$ at 412 nm, and a standard deviation σ_g = 0.92. The size distribution was integrated from r₁ = 0.005 μm to 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. Figure 6 shows reflected (top) and transmitted (bottom) results for the Stokes parameter I (τ, μ, ϕ) for solar beam incidence at 60° solar zenith angle on a homogeneous slab of optical thickness 0.3262 overlying a black surface. Similar results for the Stokes parameter Q (τ, μ, ϕ) are shown in Figure 7 and for the Stokes parameters U (τ, μ, ϕ) and V (τ, μ, ϕ) in Figure 8.

FIGURE 6

FIGURE 6. The Stokes parameter I (τ, μ, ϕ) for the aerosol scattering case for reflected (top) and transmitted (bottom) light. Number of discrete ordinate streams, NSTR = 148.

FIGURE 7

FIGURE 7. The Stokes parameter Q (τ, μ, ϕ) for the aerosol scattering case for reflected (top) and transmitted (bottom) light. NSTR = 148.

FIGURE 8

FIGURE 8. The Stokes parameter U (τ, μ, ϕ) (left) and V (τ, μ, ϕ) (right) for the aerosol scattering case for reflected (top) and transmitted (bottom) light. The U and V components vanish for ϕ =0° and ϕ =180°. NSTR = 148.

6.1.2 Cloud Case

Benchmark results for a homogeneous slab of optical thickness 5 consisting of a log-normal distribution of cloud particles with r_g = 5 μm (mode size parameter of $\sim 38$ at 412 nm), and σ_g = 0.4 were also provided by Kokhanovsky et al. (2010). The smallest and largest particle radii were selected to be r₁ = 0.005 μm and r₂ = 100 μm, and the refractive index was set to m = 1.339. These choices yield a single-scattering albedo of 1.0, and an asymmetry factor g = 0.86114 for this ensemble of cloud particles. Figures 9, 10 show that VDISORT simulations with NSTR = 148 yield good agreement with the benchmark results for all four Stokes parameters.

FIGURE 9

FIGURE 9. Upper panel: The Stokes parameter I (τ, μ, ϕ) for the cloud scattering case for reflected (top) and transmitted (bottom) light. NSTR = 148. Lower panel: Same as the upper panel except for the Stokes parameter Q (τ, μ, ϕ).

FIGURE 10

FIGURE 10. The Stokes parameter U (τ, μ, ϕ) (left) and V (τ, μ, ϕ) (right) for the cloud scattering case for reflected (top) and transmitted (bottom) light. The U and V components vanish for ϕ =0° and ϕ =180°. NSTR = 148.

6.2 The Bidirectional Polarized Reflectance Distribution Function (BP_rDF)

The BP_rDF is the vector equivalent that corresponds to the Bidirectional Reflectance Distribution Function (BRDF) in scalar radiative problems for which only the first component of the Stokes vector, the radiance I, is considered. In this version of VDISORT, BP_rDFs for two surface types were implemented: a Lambertian surface and a rough surface with a Gaussian distribution of surface slopes, which is frequently used to model scattering from wind-roughened water surfaces.

In general, both the diffuse light and the direct beam are reflected by BP_rDF shown as below:

${\mathbf{I}}_{\text{refl}}\left({\tau }^{\ast },\mu ,\varphi \right)={\int }_0^{1}d{\mu }^{\prime }{\int }_0^{2\pi }d{\varphi }^{\prime }{\mu }^{\prime }\widetilde{\mathbf{\text{R}}}\left({\mu }^{\prime },\mu ,\varphi -{\varphi }^{\prime }\right){\mathbf{I}}_{\text{inc}}\left({\tau }^{*},{\mu }^{\prime },{\varphi }^{\prime }\right)+\frac{{\mu }_0}{4\pi }\tilde{\mathbf{R}}\left({\mu }_0,\mu ,\varphi -{\varphi }_0\right){\mathbf{S}}_b{e}^{-{\tau }^{*}/{\mu }_0}(145)$

where μ and μ′ are cosines of the polar angles θ and θ′ and ϕ′ and ϕ are the corresponding azimuth angles. The downward Stokes vector at the surface is denoted I_inc (τ*, μ′, ϕ′), and ${\mathbf{S}}_b={[I_0/2,I_0/\mathrm{2,0,0}]}^T$ is the TOA direct beam illumination assumed to consist of unpolarized light of irradiance I₀. The 4 × 4 reflection matrix $\tilde{\mathbf{R}}$ depends on the surface properties. It is expanded into a Fourier series to isolate the azimuth dependence and that expansion is consistent with the expansion of the phase matrix.

A Lambertian surface is a special surface that, regardless of the state of polarization of the incident radiation, gives rise to reflected radiation that is uniform, i.e. isotropic over the upward hemisphere, and unpolarized. Therefore, only the first (m = 0) term in the Fourier expansion contributes. The reflection matrix is given by

${\tilde{\mathbf{R}}}^{m=0}={\rho }_{\text{L}}\left(\begin{array}{cccc}\hfill 0.5\hfill & \hfill 0.5\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0.5\hfill & \hfill 0.5\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right)(146)$

where ρ_L is the surface albedo.

For a wind-roughened water surface, an explicit expression for the reflectance matrix $\tilde{\mathbf{R}}$ is given by

$\tilde{\mathbf{R}}\left({\mu }^{\prime },\mu ,\varphi -{\varphi }^{\prime }\right)=\frac{1}{4{\mu }^{\prime }\mu {\mu }_n^4}p\left({\mu }_n\right)\cdot {\mathbf{C}}_{\text{rs}}^r\left({\mu }^{\prime },\mu ;{\varphi }^{\prime },\varphi \right)\cdot \mathbf{S}\left({\mu }^{\prime },\mu ,\sigma \right).(147)$

In Eq. 147, the matrix ${\mathbf{C}}_{\text{rs}}^r$ is determined by the relative refractive index m and is derived from the Fresnel reflectance, with details described in Supplementary Appendix A3 of Stamnes and Stamnes (2015). p (μ_n) is the rough surface slope probability approximated by a one-dimensional Gaussian distribution (Cox and Munk, 1954).

$p\left({\mu }_n\right)=\frac{1}{\pi {\sigma }^2}\exp \left(-\frac{1-{\mu }_n^2}{{\sigma }^2{\mu }_n^2}\right)(148)$

${\mu }_n=\frac{\mu +{\mu }^{\prime }}{\sqrt{2\left(1-\cos \mathrm{\Theta }\right)}}(149)$

$\cos \mathrm{\Theta }=-{\mu }^{\prime }\mu +\sqrt{1-{\mu }^{\prime 2}}\sqrt{1-{\mu }^2}\cos \left({\varphi }^{\prime }-\varphi \right)(150)$

${\sigma }^2=0.003+0.00512\cdot w.(151)$

σ² is the mean square surface slope determined by the water surface wind speed w in m s⁻¹. S (μ, μ′, σ) is the shadow term that is only important for the conditions of large wind speeds and large viewing zenith angles (Tsang et al., 1985). To demonstrate how the surface reflectance is applied, a simple example of an atmospheric simulation is considered with a layer of non-absorbing molecules of scattering optical thickness 0.32 and a depolarization factor of 0.04 at 412 nm (Bodhaine et al., 1999) overlying wind-roughened water surfaces with different wind speeds (2, 5, 10 m s⁻¹). The solar zenith angle is set to 30°. Figure 11 shows the upward I, Q, and U components at the top of atmosphere. A glint pattern is clearly evident in the principal plane (Δϕ = (ϕ′ − ϕ) = 0°) for all wind speeds.

FIGURE 11

FIGURE 11. TOA polarization components of sunlight reflected from a molecular atmosphere overlying a wind-roughened water surface for three different wind speeds.

6.3 The Single-Scattering Correction

To demonstrate the efficiency and accuracy gained by use of the single-scattering correction (SSC), we provide some examples in this Section. For strongly forward-peaked scattering occurring for particles that are large compared to the light wavelength, use of the delta-M scaling transformation is very useful. Also, use of the ISF method discussed in Section 3.5 helps producing accurate results at polar angles other than the quadrature angles. Figure 12 shows that accurate results for the total polarized radiance (the I Stokes parameter) can be obtained with as little as 16 streams when the SSC is applied in addition to the delta-M and ISF methods. Similar results for the Q and U Stokes parameters are shown in Figures 13, 14.

FIGURE 12

FIGURE 12. The Stokes parameter I (τ, μ, ϕ) for the aerosol scattering case. Top row: Reflected light. Same as top row of Figure 6, but using NSTR = 16 (NSTR = 148 in Figure 6) and with and without SSC. Bottom row: Transmitted light; same as bottom row of Figure 6, but using NSTR = 16 (NSTR = 148 in Figure 6) and with and without SSC.

FIGURE 13

FIGURE 13. The Stokes parameter Q (τ, μ, ϕ) for the aerosol scattering case. Top row: Reflected light. Same as top row of Figure 7, but using NSTR = 16 (NSTR = 148 in Figure 7) and with and without SSC. Bottom row: Transmitted light; same as bottom row of Figure 7, but using NSTR = 16 (NSTR = 148 in Figure 7) and with and without SSC.

FIGURE 14

FIGURE 14. The Stokes parameter U (τ, μ, ϕ) for the aerosol scattering case. Top row: Reflected light. Same as top row of Figure 8, but using NSTR = 16 (NSTR = 148 in Figure 8) and with and without SSC. Bottom row: Transmitted light; same as bottom row of Figure 8, but using NSTR = 16 (NSTR = 148 in Figure 8) and with and without SSC.

6.4 Polarized Beam Incidence

So far we have assumed that the incident beam consisted of natural (unpolarized) light (like sunlight). For some applications, like laser (lidar) or lunar beam illumination, the source would be a polarized beam. VDISORT is capable of handling also the general case of an arbitrarily polarized incident beam source. To test the performance of VDISORT for polarized beam incidence, a polarized beam source term I_inc = [I,Q,U,V]^T = π[1.0,−0.4.0.2.0.05]^T (corresponding to ${\mathbf{I}}_{\text{inc}}={[I_{\parallel },I_{\perp },U,V]}^T=\pi {[0.3,0.7,0.2,0.05]}^T$) was chosen as input for the so-called L = 13 case reported by Garcia and Siewert (1989). In this test case the optical thickness of the slab was assumed to be 1.0, the single-scattering albedo was taken to be 0.99, and the surface was assumed to be a Lambertian reflector with albedo 0.1. We calculated the Stokes parameters I, Q, and U and the degree of polarization (DOP) for this test case and reproduced Figure 5 of Schulz et al. (1999) as shown in Figure 15. Schulz et al. (1999) verified their results by comparing with output produced by the accurate General Adding Program (GAP) described by De Haan et al. (1987).

FIGURE 15

FIGURE 15. The Stokes parameters for polarized beam incidence with ${\mathbf{I}}_{\text{inc}}={[I_{\parallel },I_{\perp },U,V]}^T=\pi {[0.3,0.7,0.2,0.05]}^T$. Same as Figure 5 of Schulz et al. (1999).

7 Concluding Remarks

A review is provided of the vector discrete ordinate (VDISORT) method of solution to the radiative transfer equation pertinent for polarized transfer of radiation in a vertically stratified medium. Several new features are described and discussed including how to 1) deal with the complex solutions required to compute the V component of the Stokes vector I = [I Q U V]^T, 2) obtain accurate radiances at any desired polar observation angles by use of the ISF method, 3) deal with polarized beam incidence at the top of the atmosphere as well as polarized reflectance by the lower boundary, 4) use a pseudo-spherical treatment to correct for Earth curvature effects, and 5) use the single-scattering correction to enhance the efficiency of the method without sacrificing accuracy. Comparisons with benchmark results are provided to demonstrate the versatility of the VDISORT computer code to provide reliable solutions for aerosol and cloud cases including non-spherical ice cloud phase matrices. In particular, it has been shown that as few as 2N = 16 discrete-ordinate streams are sufficient to compute accurate polarized radiances for phase matrices appropriate for ensembles of aerosol particles. We encourage future users to help us improve this freely available tool by 1) reporting on bugs found and how they were fixed, 2) making suggestions for how this tool can be improved, and 3) help make it known to friends and co-workers in need of such a resource that this tool is available.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

KS wrote the first draft of the paper based on notes and material provided by ZL. AB helped with the formulation of the equations in Section 2.4.1. All authors provided feedback on the draft including corrections and constructive suggestions. The numerical coding was done primarily by ZL with significant input from SS. Inspection/testing of the numerical VDISORT code was done by WL and IL, who provided substantial suggestions for improvement.

Funding

Partial support for the work described in this paper was provided by grants to Stevens Institute of Technology from the Air Force Research Laboratory as well as NASA’s Remote Sensing Theory Program.

Conflict of Interest

AB is employed by Spectral Sciences.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsen.2022.880768/full#supplementary-material

Footnotes

¹We are ignoring the contribution ${\mathbf{S}}_b{e}^{-{\tau }_L/{\mu }_0}\cdot \rho ({\mu }_0,\mu ,\varphi -{\varphi }_0)$ from the lower boundary, which is typically small.

References

Berk, A. (2022). A Note on Computing Stokes Phase Matrices for Macroscopically Isotropic Mirror Symmetric Scattering Medium Stokes Scattering Matrices. J. Quantitative Spectrosc. Radiat. Transf. 286, 108220. doi:10.1016/j.jqsrt.2022.108220

CrossRef Full Text | Google Scholar

Bodhaine, B. A., Wood, N. B., Dutton, E. G., and Slusser, J. R. (1999). On Rayleigh Optical Depth Calculations. J. Atmos. Ocean. Technol. 16, 1854–1861. doi:10.1175/1520-0426(1999)016<1854:orodc>2.0.co;2

CrossRef Full Text | Google Scholar

Chandrasekhar, S. (1960). Radiative Transfer. New York: Dover Publications.

Google Scholar

Cohen, D., Stamnes, S., Tanikawa, T., Sommersten, E. R., Stamnes, J. J., Lotsberg, J. K., et al. (2013). Comparison of Discrete Ordinate and Monte Carlo Simulations of Polarized Radiative Transfer in Two Coupled Slabs with Different Refractive Indices. Opt. Expr. 21, 9592–9614. doi:10.1364/oe.21.009592

PubMed Abstract | CrossRef Full Text | Google Scholar

Cox, C., and Munk, W. (1954). Measurement of the Roughness of the Sea Surface from Photographs of the Sun’s Glitter. Josa 44, 838–850. doi:10.1364/josa.44.000838

CrossRef Full Text | Google Scholar

Dahlback, A., and Stamnes, K. (1991). A New Spherical Model for Computing the Radiation Field Available for Photolysis and Heating at Twilight. Planet. Space Sci. 39, 671–683. doi:10.1016/0032-0633(91)90061-e

CrossRef Full Text | Google Scholar

De Haan, J., Bosma, P., and Hovenier, J. (1987). The Adding Method for Multiple Scattering Calculations of Polarized Light. Astronomy Astrophysics 183, 371–391.

Google Scholar

Garcia, R., and Siewert, C. (1989). The Fn Method for Radiative Transfer Models that Include Polarization Effects. J. Quantitative Spectrosc. Radiat. Transf. 41, 117–145. doi:10.1016/0022-4073(89)90133-7

CrossRef Full Text | Google Scholar

Hansen, J. E. (1971). Multiple Scattering of Polarized Light in Planetary Atmospheres Part Ii. Sunlight Reflected by Terrestrial Water Clouds. J. Atmos. Sci. 28, 1400–1426. doi:10.1175/1520-0469(1971)028<1400:msopli>2.0.co;2

CrossRef Full Text | Google Scholar

He, X., Stamnes, K., Bai, Y., Li, W., and Wang, D. (2018). Effects of Earth Curvature on Atmospheric Correction for Ocean Color Remote Sensing. Remote Sens. Environ. 209, 118–133. doi:10.1016/j.rse.2018.02.042

CrossRef Full Text | Google Scholar

Hovenier, J. W., der Mee, C. D. V., and Domke, H. (2004). Transfer of Polarized Light in Planetary Atmospheres. Dordrecht, Netherlands: Kluwer Academic Publsihers.

Google Scholar

Hovenier, J. W., and van der Mee, C. V. M. (1983). Fundamental Relationships Relevant to the Transfer of Polarized Light in a Scattering Atmosphere. Astron. Astrophys. 128, 1–16.

Google Scholar

Hu, Y.-X., Wielicki, B., Lin, B., Gibson, G., Tsay, S.-C., Stamnes, K., et al. (2000). δ-Fit: A Fast and Accurate Treatment of Particle Scattering Phase Functions with Weighted Singular-Value Decomposition Least-Squares Fitting. J. Quantitative Spectrosc. Radiat. Transf. 65, 681–690. doi:10.1016/s0022-4073(99)00147-8

CrossRef Full Text | Google Scholar

Kokhanovsky, A. A., Budak, V. P., Cornet, C., Duan, M., Emde, C., Katsev, I. L., et al. (2010). Benchmark Results in Vector Atmospheric Radiative Transfer. J. Quantitative Spectrosc. Radiat. Transf. 111, 1931–1946. doi:10.1016/j.jqsrt.2010.03.005

CrossRef Full Text | Google Scholar

Laszlo, I., Stamnes, K., Wiscombe, W. J., and Tsay, S.-C. (2016). “The Discrete Ordinate Algorithm, Disort for Radiative Transfer,” in Light Scattering Reviews (Berlin Heidelberg: Springer), Vol. 11, 3–65. doi:10.1007/978-3-662-49538-4_1

CrossRef Full Text | Google Scholar

Lin, Z., Chen, N., Fan, Y., Li, W., Stamnes, K., and Stamnes, S. (2018). New Treatment of Strongly Anisotropic Scattering Phase Functions: the Delta-M+ Method. J. Atmos. Sci. 75, 327–336. doi:10.1175/jas-d-17-0233.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, Z. (2016). Development and Applications of Radiative Transfer Models for Unpolarized and Polarized Light (Hoboken, NJ: Stevens Institute of Technology). PhD Thesis.

Lin, Z., Stamnes, S., Jin, Z., Laszlo, I., Tsay, S.-C., Wiscombe, W., et al. (2015). Improved Discrete Ordinate Solutions in the Presence of an Anisotropically Reflecting Lower Boundary: Upgrades of the Disort Computational Tool. J. Quantitative Spectrosc. Radiat. Transf. 157, 119–134. doi:10.1016/j.jqsrt.2015.02.014

CrossRef Full Text | Google Scholar

Mishchenko, M. I. (1991). Light Scattering by Randomly Oriented Rotationally Symmetric Particles. J. Opt. Soc. Am. A 8, 871–882. doi:10.1364/josaa.8.000871

CrossRef Full Text | Google Scholar

Mishchenko, M. I., and Travis, L. D. (1997). Satellite Retrieval of Aerosol Properties over the Ocean Using Polarization as Well as Intensity of Reflected Sunlight. J. Geophys. Res. 102, 16989–17013. doi:10.1029/96jd02425

CrossRef Full Text | Google Scholar

Nakajima, T., and Tanaka, M. (1988). Algorithms for Radiative Intensity Calculations in Moderately Thick Atmospheres Using a Truncation Approximation. J. Quantitative Spectrosc. Radiat. Transf. 40, 51–69. doi:10.1016/0022-4073(88)90031-3

CrossRef Full Text | Google Scholar

Schulz, F., and Stamnes, K. (2000). Angular Distribution of the Stokes Vector in a Plane-Parallel, Vertically Inhomogeneous Medium in the Vector Discrete Ordinate Radiative Transfer (Vdisort) Model. J. Quantitative Spectrosc. Radiat. Transf. 65, 609–620. doi:10.1016/s0022-4073(99)00115-6

CrossRef Full Text | Google Scholar

Schulz, F., Stamnes, K., and Weng, F. (1999). Vdisort: an Improved and Generalized Discrete Ordinate Method for Polarized (Vector) Radiative Transfer. J. quantitative Spectrosc. Radiat. Transf. 61, 105–122. doi:10.1016/s0022-4073(97)00215-x

CrossRef Full Text | Google Scholar

Siewert, C. (2000). A Discrete-Ordinates Solution for Radiative-Transfer Models that Include Polarization Effects. J. Quantitative Spectrosc. Radiat. Transf. 64, 227–254. doi:10.1016/s0022-4073(99)00006-0

CrossRef Full Text | Google Scholar

Siewert, C. E. (1981). On the Equation of Transfer Relevant to the Scattering of Polarized Light. Astrophys. J. 245, 1080–1086. doi:10.1086/158884

CrossRef Full Text | Google Scholar

Siewert, C. E. (1982). On the Phase Matrix Basic to the Scattering of Polarized Light. Astron. Astrophys. 109, 195–200.

Google Scholar

Sommersten, E. R., Lotsberg, J. K., Stamnes, K., and Stamnes, J. J. (2010). 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. Transf. 111, 616–633. doi:10.1016/j.jqsrt.2009.10.021

CrossRef Full Text | Google Scholar

Stamnes, K., and Conklin, P. (1984). A New Multi-Layer Discrete Ordinate Approach to Radiative Transfer in Vertically Inhomogeneous Atmospheres. J. Quantitative Spectrosc. Radiat. Transf. 31, 273–282. doi:10.1016/0022-4073(84)90031-1

CrossRef Full Text | Google Scholar

Stamnes, K. (1982). On the Computation of Angular Distributions of Radiation in Planetary Atmospheres. J. Quantitative Spectrosc. Radiat. Transf. 28, 47–51. doi:10.1016/0022-4073(82)90096-6

CrossRef Full Text | Google Scholar

Stamnes, K., and Stamnes, J. J. (2015). Radiative Transfer in Coupled Environmental Systems. Hoboken, NJ: Wiley VCH.

Google Scholar

Stamnes, K., Thomas, G. E., and Stamnes, J. J. (2017a). Radiative Transfer in the Atmosphere and Ocean. 2 edn. New York, NY: Cambridge University Press.

Google Scholar

Stamnes, K., Tsay, S.-C., Wiscombe, W., and Jayaweera, K. (1988). Numerically Stable Algorithm for Discrete-Ordinate-Method Radiative Transfer in Multiple Scattering and Emitting Layered Media. Appl. Opt. 27, 2502–2509. doi:10.1364/ao.27.002502

PubMed Abstract | CrossRef Full Text | Google Scholar

Stamnes, S., Ou, S., Lin, Z., Takano, Y., Tsay, S., Liou, K., et al. (2017b). Polarized Radiative Transfer of a Cirrus Cloud Consisting of Randomly Oriented Hexagonal Ice Crystals: The 3× 3 Approximation for Non-spherical Particles. J. Quantitative Spectrosc. Radiat. Transf. 193, 57–68. doi:10.1016/j.jqsrt.2016.07.001

CrossRef Full Text | Google Scholar

Tsang, L., Kong, J. A., and Shin, R. T. (1985). Theory of Microwave Remote Sensing. Hoboken, NJ: A Wiley-Interscience Publication.

Google Scholar

Weng, F. (1992). A Multi-Layer Discrete-Ordinate Method for Vector Radiative Transfer in a Vertically-Inhomogeneous, Emitting and Scattering Atmosphere – I. Theory. J. Quantitative Spectrosc. Radiat. Transf. 47, 19–33. doi:10.1016/0022-4073(92)90076-g

CrossRef Full Text | Google Scholar

Wiscombe, W. (1977). The Delta–M Method: Rapid yet Accurate Radiative Flux Calculations for Strongly Asymmetric Phase Functions. J. Atmos. Sci. 34, 1408–1422. doi:10.1175/1520-0469(1977)034<1408:tdmrya>2.0.co;2

CrossRef Full Text | Google Scholar