Vector Radiative Transfer in a Vertically Inhomogeneous Scattering and Emitting Atmosphere. Part I: A New Discrete Ordinate Method

1.   Introduction

Radiative transfer models (RTMs) have been widely utilized in space sensor simulations, satellite data assimilation, and remote sensing algorithms (Weng, 2007; Boukabara et al., 2011; Kan et al., 2020). The doubling and adding method is developed to solve the radiative transfer equation (Evans and Stephens, 1991) and applied in the fast radiative transfer models such as the Community Radiative Transfer Model (CRTM; Liu and Weng, 2006) and Advanced Radiative Transfer Modeling System (ARMS; Weng et al., 2020). The successive order of scattering (SOS; Min and Duan, 2004; Duan et al., 2010) approach and successive-order-of-interaction (SOI; Heidinger et al., 2006; O’Dell et al., 2006) method are also used in RTM. As an example of SOS, Second Simulation of a Satellite Signal in the Solar Spectrum (6S) is developed and widely used for radiative transfer in visible wavelengths (Vermote et al., 1997). There are also some other methods to solve the radiative transfer equation, such as the Monte Carlo method (Lu and Hsu, 2005), discrete ordinate adding approximation method (Li et al., 2022), matrix-operator method (Liu and Weng, 2013), and others (Emde et al., 2018; Wang et al., 2020).

As a powerful technique, the discrete ordinate method is used for solving the vector radiative transfer equation. Originally, the discrete ordinate radiative transfer (DISORT) model was developed to solve the radiance at the top of atmosphere (Stamnes et al., 1988a). Then, a vector discrete ordinate radiative transfer (VDISORT) model was first developed by Weng (1992) and takes into account the Stokes vector of radiation. Based on VDISORT, a few studies are performed on the radiative transfer simulations and satellite applications (Barichello and Siewert, 1999), including the integration of the source function method to enable arbitrary output (Stamnes and Swanson, 1981; Schulz and Stamnes, 2000), a vector linearized discrete ordinate radiative transfer (VLIDORT) model for retrieval applications (Spurr, 2004, 2006), and so on. Several approximate methods can also be performed together with the discrete ordinate method to enhance the efficiency of VDISORT (Liou et al., 2005; Wang and Wu, 2011). In general, the Fourier expansion and the Gaussian quadrature schemes are used to approximate the double integrations of multiple scattering terms in VDISORT (Marshak and Davis, 2005). The Fourier expansion separates the Stokes elements [i.e., $\left({I}_l,{I}_r,{I}_u,{I}_v\right)^{\mathrm{T}}$] into the cosinusoidal and sinusoidal modes with respect to the relative azimuthal angle. For polarized RTMs, the cosinusoidal mode of the first two components in the Stokes vector is combined with the sinusoidal modes of the 3rd and 4th components as a new vector [i.e., $\left({I}_{_l}^{\mathrm{c}},{I_r^{\mathrm{c}}},{I}_u^{\mathrm{s}},{I}_v^{\mathrm{s}}\right)^{\mathrm{T}}$]. The cosinusoidal mode of an ordinary differential RTM is then formed for a solution. This combination works well with the properties of the scattering phase matrix of spherical particles. The other elements form the Stokes vector [i.e., $\left({I}_{_l}^{\mathrm{s}},{I}_r^{\mathrm{s}},{I}_u^{\mathrm{c}},{I}_v^{\mathrm{c}}\right)\mathrm{^T}$] in the sinusoidal mode of RTM. In the past, many studies on RTM focused on the solution of the cosinusoidal mode of the combined Stokes vector (Liu and Weng, 2006, 2013; Natraj et al., 2023). It performs well for the radiative transfer issues with nonpolarized sources in Mie or Rayleigh scattering atmosphere (Spurr, 2004, 2006; Kokhanovsky et al., 2010; Shi et al., 2021). The contribution from the combined sine-mode Stokes vector to the total radiance also needs further research, especially for the cases of radiative transfer with the polarized sources. Considering the sinusoidal mode together with the cosinusoidal mode of the Stokes vector, a typical case with Mie scattering was tested under a polarized source to derive the total Stokes radiance (Schulz et al., 1999; Lin et al., 2022).

Based on similar principles, several RTMs have been used to investigate the radiative transfer problems, including RTMs based on the discrete ordinate, matrix-operator, or adding-doubling method (He et al., 2010; Lin et al., 2022). Besides, by source function integration (Schulz and Stamnes, 2000; Spurr, 2006) or simple adjustments on the integral scheme (Evans and Stephens, 1990; Heidinger et al., 2006), RTM can maintain the arbitrary zenith output and therefore can calculate the Stokes signal at arbitrary viewing geometry. Additional efforts have been made to deal with the complex eigenvalue problem in the homogeneous solution (Stamnes et al., 1988b; Lin et al., 2022). However, these RTMs all rely on an assumption that the atmospheric scattering phase matrix and the earth-surface reflective matrix have the same azimuthal symmetric characteristics as Mie scattering phase matrix. However, in some circumstances, the scattering phase matrix or the boundary reflection matrix does not obey this symmetric relationship. The actual sea-surface reflection should be more related to the surface winds rather than the relative azimuthal differences to the incident radiance (Cox and Munk, 1954; Yueh, 1997). Atmospheric scattering by non-spherical particles with a fixed orientation also has different characteristics from Mie scattering (Mishchenko, 2000). Since the non-spherical particle scattering theory develops and the polarized bidirectional reflection distribution function (pBRDF) has already been developed (Yueh, 1997; Fan et al., 2010; Liu et al., 2011; He and Weng, 2023), it is necessary to expand the discrete ordinate radiative transfer theory to the normal cases of the atmospheric scattering and the boundary reflection without relying on the azimuthal symmetric assumption.

In this study, an updated VDISORT model is derived for solving the radiative transfer equation (RTE) in vertically inhomogeneous scattering and emitting atmospheres. The model extends the discrete ordinate radiative transfer theory applicable for the general atmospheric and surface scattering and emitting medium. The VDISORT combines both the cosinusoidal and sinusoidal modes of all four Stokes components and outputs the complete Stokes signal. The complex eigenvalue is used in the VDISORT codes. Using the integral scheme mentioned by Heidinger et al. (2006), the model can also output the Stokes radiance at arbitrary viewing geometry. In this paper, Section 2 presents the theoretical basis of VDISORT. In Section 3, the results are validated against benchmarks for Rayleigh, Mie, and random-oriented non-spherical particle scattering atmospheres (Coulson, 1960; Garcia and Siewert, 1986, 1989; Siewert, 2000). With a polarized source, the significance of the sinusoidal mode information of RTM is also discussed in Section 3. Section 4 presents the summary and conclusions.

2.   Theory and methodology

2.1   Radiative transfer equation

In a plane-parallel atmosphere, radiative transfer equation (RTE) is expressed as (Weng, 1992)

where ${\boldsymbol{I}}(\tau ,\mu ,\varphi )$ denotes the Stokes vector $\left({I}_l,{I}_r,{I}_u,{I}_v\right)^{\mathrm{T}}$, which is the function of optical depth $\tau$, cosine of zenith angle $\mu{\text{ }}$and azimuthal angle $\phi$; $\mu{\text{ }}$is positive for the upwelling radiation and negative for downwelling radiation; $\omega (\tau )$ is the single scattering albedo at the corresponding atmosphere layer; ${\boldsymbol{M}}(\tau ,\mu ,\phi ;\mu ',\phi ')$ is the scattering phase matrix, where $(\mu ',\phi ')$ and $(\mu ,\phi )$ represent the incident and scattering direction, respectively; $\boldsymbol{S}_{\mathrm{b}}$ is the internal beam source radiation from the direction $( - {\mu _0},{\phi _0})$ at the top of atmosphere; $\boldsymbol{S}_{\mathrm{t}}(\tau)$ is the thermal source of radiation from atmosphere. Taking $B[T(\tau )]$ as the Planck function at temperature $T(\tau )$, the thermal source is assumed to be unpolarized and is expressed as $\boldsymbol{S}\mathrm{_t}\left(\tau\right)=\left(\dfrac{B[T(\tau)]}{2}, \dfrac{B[T(\tau)]}{2},0,0\right)^{\mathrm{T}}$.

Expanding ${\boldsymbol{I}}(\tau ,\mu ,\varphi )$ and ${\boldsymbol{M}}(\tau ,\mu ,\phi ;\mu ',\phi ')$ into Fourier series as follows:

Considering $\boldsymbol{S}_{\mathrm{b}}$ as a constant, ${\boldsymbol{Q}}(\tau ,\mu ,\phi )$ can be written by

where ${\delta _{0m}} = \left\{ \begin{gathered} 1,{\text{ }}m = 0 \\ 0,{\text{ }}m \ne 0 \\ \end{gathered} \right.$. Set ${\boldsymbol{M}}\left( {\tau ,\mu ,\phi ;\mu ',\phi '} \right) =( {m_{pq}}; p,$ q = 1,2,3,4). For spherical particle scattering, ${\boldsymbol M}^{\mathrm{Mie}}$$(\tau,\mu,\phi;$$\mu',\phi')$, the four elements in the upper left and lower right corner are even functions of $\phi ' - \phi$, with other elements odd functions of $\phi ' - \phi$. It means the Fourier expanded item of $\boldsymbol{M}\mathrm{^{Mie}}\left(\tau,\mu,\phi;\mu',\phi'\right)$ can be written as

where ${\boldsymbol{M}}_{m,pq}^{\rm{{c,s}}} = \left( {\begin{array}{*{20}{c}} {m_{m,\left( {2p - 1} \right)\left( {2q - 1} \right)}^{\rm{{c,s}}}}&{m_{m,\left( {2p - 1} \right)\left( {2q} \right)}^{\rm{{c,s}}}} \\ {m_{m,\left( {2p} \right)\left( {2q - 1} \right)}^{\rm{{c,s}}}}&{m_{m,\left( {2p} \right)\left( {2q} \right)}^{\rm{{c,s}}}} \end{array}} \right), p,q = 1,2$.

Using the orthogonality of trigonometric functions and applying the Gaussian or double-Gaussian integration method to replace the zenith integral, the cosinusoidal and the sinusoidal modes of RTE at each Fourier harmonic component was traditionally derived as (Schulz et al., 1999; Lin et al., 2022)

where $\boldsymbol{I}_m^{\mathrm{c,s}}(\tau,\mu)=\left[{I}_{m,l}^{\mathrm{c,s}},{I}_{m,r}^{\mathrm{c,s}},{I}_{m,u}^{\mathrm{c,s}}, {I}_{m,v}^{\mathrm{c,s}}\right]\mathrm{^T}$; ${\mu _j},{w_j}$ are the Gaussian or double-Gaussian integral points and corresponding integral weights, respectively; N is the number of the Gaussian integral points for double Gaussian method in the range $\left[ {0,1} \right]$ for $\mu$; ${\mu _{ - j}} = - {\mu _j},{\text{ }}{w_{ - j}} = {{{w}}_j}$ are used to derive the Gaussian integral points and weights for $\mu {\text{ }}$in $\left[ { - 1,0} \right]$. The cosinusoidal mode equation can be derived from Eqs. (5a) and (5b) as:

where $2N$ is the number of streams, $\boldsymbol{I}_{m,lr}^{\mathrm{c,s}}(\tau,\mu)=\left({I}_{m,l}^{\mathrm{c,s}},{I}_{m,r}^{\mathrm{c,s}}\right)\mathrm{^T},\text{ }\boldsymbol{I}_{m,uv}^{\mathrm{c,s}}(\tau,\mu)=\left({I}_{m,u}^{\mathrm{c,s}},{I}_{m,v}^{\mathrm{c,s}}\right)\mathrm{^T}$, and so is $\boldsymbol{Q}_m^{\mathrm{c,s}}$.

However, Eqs. (5) and (6) clearly rely on the properties of the Mie scattering phase matrix in Eq. (4). It cannot be used under all circumstances of atmospheric scattering. For the general scattering case, the Fourier harmonics of the phase matrix should be

Thus, the new combinations should take the form as follows:

Similarly, the sinusoidal mode equation of RTE should be

2.2   Combination of the cosinusoidal and sinusoidal modes

Combining the cosinusoidal and sinusoidal modes, the equation of RTM at each Fourier harmonic is expressed as

To make the solution complete and universal, the sinusoidal and cosinusoidal modes of the radiation can be solved at the same time through Eq. (10). In Eq. (10), it is obvious that there are eight ordinary differential equations at a specified μ. Setting $\mu = {\mu _{ \pm j}}$, $16N$ equations are then formed. With a known phase matrix ${\boldsymbol{M}}\left( {\tau ,\mu ,\phi ;\mu ',\phi '} \right)$, single scattering albedo $\omega (\tau )$, input source $\boldsymbol{S}_{\mathrm{b}}$, and thermal source $\boldsymbol{S}\mathrm{_t}(\tau)$, there are $16N$ unknown Stokes elements in the different $16N$ equations for each Fourier harmonic, making RTE a closure for solving all the unknowns.

Compared to the previous Eqs. (5) and (6), which was used in Schulz et al. (1999) and Lin et al. (2022), Eqs. (8)–(10) are more generalized for all the atmospheric scattering cases without dependence on the properties of Mie scattering phase matrix. With the similar forms taken into account the boundary reflection, the updated VDISORT in this paper is well-coupled for the real surface reflection without azimuthal symmetric assumptions. It is also a great step for considering the fixed-oriented non-spherical atmospheric scattering, although there are still the vector characteristics of the extinction coefficients under this circumstance that needs consideration.

2.3   Calculation of radiances at arbitrary zenith angle

If the incident and scattering zenith angles are set in the same Gaussian quadrature scheme, RTE can be solved only at the Gaussian points. An interpolation is often used for calculating the radiances at the angles other than Gaussian points (Weng, 1992; Schulz et al., 1999). The interpolation can be problematic when $N$ is small. In this paper, two additional points are added to the Gaussian quadrature schemes with their weights being set to zero (Heidinger et al., 2006). The revised Gaussian quadrature scheme is:

Since the weights of two added points are zero, it will have no impact on Eq. (10) if adding these two points into the equation. It should be noted that all the points must satisfy ${\mu _j} \ne - {\mu _0}$. A small increment is added to $\mu$ when $\mu = - {\mu _0}$. The same processing is also operated when $\mu = 0$. After taking these two points into Eq. (10), the new set of equations is expressed as

where the radiance at the specified zenith angle is now available through the two extra Gaussian points. When it comes to the downwelling radiation, the solution at ${\mu _{ - \left( {N + 1} \right)}}$ is used. The solution at ${\mu _{N + 1}}$ is for the upwelling radiation. By adding two Gaussian points with zero weights, RTE can be solved at arbitrary zenith angles without interpolation.

2.4   VDISORT solution

If the total atmosphere is treated as a combination of a number of the homogeneous atmospheric layers, RTE should be solved for each homogeneous layer to derive the solution at arbitrary optical depth. Figure 1 shows a schematic diagram of the overall inhomogeneous atmosphere containing $L$ homogeneous atmospheric layers. In each layer, the optical parameters such as the single scattering albedo ${\omega _i}$ and the phase matrix ${{\boldsymbol{M}}_i}(\tau ,\mu ,\phi ;\mu ',\phi ')$ are unique. The theoretical approach to solve Eq. (12) in each homogenous layer is similar to Weng (1992). The differences are the size of the matrix and the number of the equations. To better illustrate the boundary condition for updated VDISORT, the solving process should be introduced briefly. Equation (12) can be simplified into

Download

Larger image

PowerPoint slide

Fig  1.  Schematic diagram of the homogeneous layers in the inhomogeneous atmosphere, with the phase matrix ${{\boldsymbol{M}}_i}(\tau ,\mu ,\phi ;\mu ',\phi ')$ and the single scattering albedo ${\omega _i}$ specified for each layer.

where

and for each pair of ${\mu _s}$ and ${\mu _j}$:

and

Given that one side of Eq. (13) has been fixed in Eqs. (14c) and (14d), the other parameters (e.g., ${{\boldsymbol{A}}_m}$) can be derived by each equation and the definitions of these parameters, such as Eq. (14b). Equation (13) is a set of ordinary differential equations and can be solved through seeking for the general and particular solutions. Apparently, the general solution should be in the form (Weng, 1992)

where

The coefficients ${C_{jmp}}$ are determined by the boundary conditions in Section 2.5; ${\lambda _{jmp}}$ and ${{\boldsymbol{g}}_{jmp}}$ are the j-th eigenvalue and eigenvector by solving the homogeneous differential equation

Substituting Eq. (15) into Eq. (16), the general solution can be derived by solving the eigenvalue issue:

The eigenvalues and the eigenvectors are the forms of complex variables with the LAPACK Fortran library. The potential complex eigenvalue problems can be taken into consideration by this approach.

The particular solution is related to the single scattering item from the beam source and thermal radiation $\boldsymbol{P}_{mp}(\tau,\mu)=\boldsymbol{P}_{\mathrm{b}mp}(\tau,\mu)+\boldsymbol{P}_{\mathrm{t}mp}(\tau,\mu)$. Since the thermal source is nonpolarized, $\boldsymbol{P}_{\mathrm{t}mp}(\tau,\mu)$ is nonzero only for the zero-order cosinusoidal component of the first two Stokes elements. Adopting the first-order polynomial approximation for $\boldsymbol{P}_{\mathrm{t}mp}(\tau,\mu)$ and the linear assumption for the thermal emission to simplify the calculation in the homogeneous layer, the differential equation and the approximate method are derived as (Weng, 1992)

and

The constants ${c_{0p}}$ and ${c_{1p}}$ should be derived by the linear assumption. Assuming that the thermal emission varies linearly between the upper and the lower boundary of each atmospheric layer, ${c_{0p}}$ and ${c_{1p}}$ can be calculated by the thermal radiation at the upper and the lower interface of the layer. When $m = 0$, Eq. (18) can be simplified as

Then, the thermal particular solution $\boldsymbol{P}_{\mathrm{t}mp}(\tau,\mu)=\delta_{0m} \left[\boldsymbol{b}_{0p}+\boldsymbol{b}_{1p}\tau\right]$ can be obtained.

The equation for the particular solution corresponding to the beam source is

where $\boldsymbol{F}_{\mathrm{b}m}$ in Eq. (21) should be:

The particular solution takes the form (Weng, 1992)

Substituting Eq. (23) into Eq. (21), the solution can be derived through

Hence, the solution for RTE should be

where

The only unknown parameters for RTE’s solution in Eq. (25) are the coefficients ${{\boldsymbol{C}}_{jmp}}$ for the general solution, which can be resolved through the continuity conditions in the lower and upper boundary, and the interfaces between the internal layers.

2.5   Boundary conditions

2.5.1   Upper boundary conditions

The upper boundary is set at the top of the atmosphere, and the downward solution for RTM is equal to the incident radiation, which mainly results from the cosmic background radiation (Weng, 1992):

Then, the upper boundary condition is

2.5.2   Lower boundary conditions

At the bottom boundary, the upwelling solution should be equal to the sum of the boundary thermal emission, the reflection of the downwelling radiation and the reflected beam source, which is expressed as

The radiative equation at the lower boundary is (Weng, 1992)

where ${\boldsymbol{B}}( + \mu ,\phi ; - \mu ',\phi ')$ is the pBRDF matrix, and ${\boldsymbol{E}}\left( { + \mu ,\phi } \right)$ is the emissivity matrix at the lower boundary. Similar to RTE, it can be rewritten by using the Fourier expansion and the Gaussian integral method as

where

Unlike the previous VDISORT versions (Weng, 1992; Schulz et al., 1999; Lin et al., 2022), Eq. (30) is derived more generally without assumption on the properties of the pBRDF matrix. Moreover, owing to the different Fourier expansion bases, Eq. (30) can be in a better coincidence with the physical implications of pBRDF. Note that pBRDF is more related to ${\phi _w}$ rather than ${\phi _0}$ (Cox and Munk, 1954; Yueh, 1997; He and Weng, 2023). If ${\boldsymbol{E}}\left( { + \mu ,\phi } \right)$ is decomposed based on $\left( {\phi - {\phi _w}} \right)$ rather than $\left( {{\phi _0} - \phi } \right)$, its Fourier harmonics can be transformed by:

The reflected downward radiation also results from the solution in Eq. (25). After re-arranging the coefficients from the general solution and the lower boundary conditions, the continuity condition at the lower boundary satisfies the following equation

where for each $+ \mu$

and for each pair of $+ \mu$ and $- {\mu _j}$

2.5.3   Continuity conditions at the internal interfaces

For the interface between atmospheric layers, the solution of the upper layer should be equal to the solution of the lower layer when it comes to the interface, implying

Substituting Eq. (25) into Eq. (36) results in

2.5.4   Coefficient processing

Based on Eqs. (27), (33), and (37), the coefficients can be solved simultaneously. Before the solving process, the scale transformation should be performed on the coefficients to ensure the stability of RTM as the previous version does (Stamnes and Conklin, 1984; Weng, 1992). The scale transformation takes the form

After the scale transformation, the coefficients ${C_{jmp}}$ will be replaced by the new ${\tilde C_{jmp}}$.

It should be noticed although the equations here is similar to those equations in Schulz et al. (1999) to some extent, the complete form of the equations has been derived first time in this study. Traditionally, RTMs based on the equations like Eqs. (5) and (6) or those equations in Schulz et al. (1999) are fine to the cases when their scattering phase matrix and boundary reflection matrix have the assuming azimuthal symmetric relationship similar to the Mie phase matrix. The equations derived here are complete form for the common case without the symmetric properties or assumptions. For example, the additional third items on the right side of Eqs. (8) and (9) are obvious evidences for its improvements, compared to Eq. (6) or Eqs. (A11) and (A12) in the earlier study by Schulz et al. (1999).

3.   Results

3.1   Rayleigh case

To validate the accuracy of this new VDISORT model, a Rayleigh case is first tested without thermal sources in Eq. (1). The surface type is Lambertian having an albedo ${\lambda _0}$ of 0.25 and other parameters are set as follows: $\boldsymbol{S}_{\rm{b}}=\left(0.5\pi,0.5\pi,0,0\right)^{\rm{T}},\text{ }\omega=1,\text{ }\mu_0=0.8,\text{ }\phi_0=0^{\circ},\rm{\ and}\text{ }\tau_L=1$. The Rayleigh phase matrix at each Fourier harmonic is calculated by the formulas in Chandrasekhar (1960). The scattering matrix and phase matrix of Rayleigh takes the form of:

where $\boldsymbol{M}_{LRUV}=\boldsymbol{L}\left(\pi-i_2\right)\boldsymbol{S}_{LRUV}\boldsymbol{L}\left(-i_1\right)$, $\boldsymbol{L}$ is the rotational matrix, and ${i_1}{\text{ and }}{i_2}$ are two rotational angles. The details of $\boldsymbol{M}_{LRUV}^{\rm{Rayleigh}}$ and ${\boldsymbol{L}}$ can be found on pages 37–42 of Chandrasekhar (1960). After the azimuthal Fourier expansion, the scattering phase matrix elements only have the zeroth to the second Fourier harmonics.

With only one atmospheric layer, Fig. 2 shows the zenith distributions of I, Q, U, and P at the lower boundary when $\phi = 90^\circ$ with the Gaussian stream $2N = 28$. Instead of the last Stokes element ${{V}}{\text{ }}$, the degree of polarization ${P}=\sqrt{{Q}^2+{U}^2+{V}^2}/{I}$ is used here, where ${{I}} = {{{I}}_l} + {{{I}}_r}, {\text{ }}{{Q}} = {{{I}}_l} - {{{I}}_r},{\text{ }}{{U}} = {{{I}}_u}$. Note that the negative sign was added to ${U}$ in the benchmark (Coulson et al., 1960) for a consistent definition with other results by Evans and Stephens (1991). Given that the definition of ${Q}\text{ }$ is different, ${Q}\text{ }$ in the benchmark is also multiplied by $- 1$. From the Rayleigh phase matrix in Eq. (39), it is obvious that no ${V}\text{ }$ component exists for unpolarized source and Lambertian surface in Rayleigh scattering.

Figure 2 illustrates the angular distribution of VDISORT radiance for Rayleigh case. Overall, the simulations by VDISORT agree well with the benchmark. The parameters ${I}$ and ${P}\text{ }$ increase first and then decrease with $\mu {\text{ }}$ ranging from −1 to 0, while ${U}$ decreases first and then increases. They have similar maximum or minimum around $\mu = - 0.3$. The parameter Q always decreases within the same range of $\mu {\text{ }}$. The maximum of the absolute bias is less than 0.0042, 0.002, 0.00003, and 0.0065 for ${I},\ {Q},\ {U},{\rm{\ and}}\ {P}$, respectively. Overall, the updated VDISORT reproduces the benchmark results for the Rayleigh case as the earlier VDISORT version (Schulz et al., 1999).

Download

Larger image

PowerPoint slide

Fig  2.  Comparisons of the zenith distributions for I, Q, U, and P at the bottom boundary for Rayleigh case when $\phi = 90^\circ$ with ${{\boldsymbol{S}}_{\rm{b}}} = {\left( {0.5\pi, 0.5\pi, 0, 0} \right)^{\rm{T}}}, {\text{ }}\omega = 1, {\text{ }}{\lambda _0} = 0.25, {\text{ }}{\mu _0} = 0.8, {\text{ }}{\phi _0} = 0^\circ, {\text{ }}{\tau _L} = 1$, between the benchmarks (circles) and the simulations of the updated VDISORT (solid line).

3.2   L13 case with spherical particle scattering

A case with a group of spherical particle scattering based on the Legendre expansion coefficients of 0th–13th orders of the phase matrix (hereinafter L13) was studied by Garcia and Siewert (1986, 1989) for a Mie scattering atmosphere with a Gamma-size distribution of spherical particles and the results are formed as another benchmark. The gamma distribution of spherical particles is with an effective radius of 0.2 μm, the effective variance of 0.07 and the index of refraction n = 1.44. The form of the atmospheric scattering matrix for this case is

With $\boldsymbol{S}_{LRUV}=\boldsymbol{D}^{-1}\boldsymbol{S}_{IQUV}\boldsymbol{D}$, $\boldsymbol{S}_{LRUV}$ can be:

Based on ${{\boldsymbol{M}}_{LRUV}} = {\boldsymbol{L}}\left( {\pi - {i_2}} \right){{\boldsymbol{S}}_{LRUV}}{\boldsymbol{L}}\left( { - {i_1}} \right)$, it is clear that the phase matrix for L13 case has all 16 elements with the symmetric relationship in Eq. (5). The expansion coefficients published in Vestrucci and Siewert (1984) for L13 case are used for calculating the Mie phase matrix from 0th to 13th Fourier harmonics at an electromagnetic wavelength of 0.951 μm, with the method described in De Haan et al., (1987). Set the Lambert surface’s albedo ${\lambda _0} = 0.1$ and $\boldsymbol{S}\mathrm{_b}=\left(0.5\pi,0.5\pi,0,0\right)\mathrm{^T}, \text{ }\omega=0.99,\text{ }\mu_0=0.2,\text{ }\phi_0=0^{\circ},\text{ }\tau_L=1$. No thermal source is included in the L13 case. Figure 3 shows the zenith distributions of ${I},\ {Q},\ {U},\ \mathrm{and}\ {V}$ at the lower boundary for $\phi = 90^\circ$ with the double-Gaussian stream $2N = 16$. Clearly, the updated VDISORT also matches well with the benchmark in Garcia and Siewert (1986, 1989).

Download

Larger image

PowerPoint slide

Fig  3.  Comparisons of the zenith distribution for I, Q, U, and V at the bottom boundary for L13 case when $\phi = 90^\circ$ with $\boldsymbol{S}_{\mathrm{b}}=\left(0.5\pi,0.5\pi,0,0\right)\mathrm{^T},\text{ }\omega=0.99,\text{ }\lambda_0=0.1,\text{ }\mu_0=0.2,\text{ }\phi_0=0^{\circ},\text{ }\tau_L=1$, between the benchmarks (circles) and the simulations of the updated VDISORT (solid line).

At the lower boundary, ${I}\ \mathrm{and}\ {Q}$ increase first and then decrease with $\mu {\text{ }}$ in $\left[ { - 1,0} \right]$ when $\phi = 90^\circ$. The maximum turns up near $\mu = - 0.5$; ${U}$ shows an opposite trend and decreases first and then increases, with its minimum near $\mu = - 0.4$; ${V}$ decreases before $\mu = - 0.9$, and then increases but finally turns back to decrease around $\mu = - 0.4$. The maximum absolute bias is less than 0.0003, 0.0006, 0.00003, and 0.0000003 for I, Q, U, and V, respectively. Considering the symmetric relationship in Eq. (4) for Mie scattering, Eq. (14b) can be rewritten as:

With an unpolarized source with the constant first two Stokes components $\left({I}_l,{I}_r,0,0\right)\mathrm{^T}$, only $\left({I}_{_l}^{\mathrm{c}},{I}_r^{\mathrm{c}},{I}_u^{\mathrm{s}},{I}\mathrm{_{\mathit{v}}^s}\right)\mathrm{^T}$ is generated in the multi-scattering process. The single scattering process is similar to the multi-scattering process. The reflection and emission of the Lambertian surface are also assumed unpolarized. Hence, the traditional VDISORT with only cosinusoidal mode is enough for accurate simulations in this case. It is clear that the updated VDISORT can also simulate the cases as the previous VDISORT version does.

3.3   Polarized source and sinusoidal mode

To investigate the effects of the sinusoidal mode of RTE, L13 case with a polarized source is tested. Referring to Schulz et al. (1999), a polarized beam source $\boldsymbol{S}_{\mathrm{b}}=\left(0.7\pi,0.3\pi,0.2\pi,0.05\pi\right)^{\mathrm{T}}$ is considered. The other parameters and methods are all the same with L13 case in Section 3.2. Although the overall ${I},\ {Q},\ {U,}$ and ${P}$ was solved in Schulz et al. (1999), all the cosinusoidal and the sinusoidal modes of the four Stokes elements ${I},\ {Q},\ {U},\ \mathrm{and}\ {V}$, instead of ${P}$, along with the final total Stokes elements derived by updated VDISORT are shown in Fig. 4. The similar trends of overall ${I},\ {Q},\ \mathrm{and}\ {U}$ are found with those in Schulz et al. (1999), which can also validate the accuracy of updated VDISORT in this paper. It also provides us with detailed information of the cosinusoidal solution $\left({I}_{_l}^{\mathrm{c}},{I}_r^{\mathrm{c}},{I}_u^{\mathrm{s}},{I}_v^{\mathrm{s}}\right)^{\mathrm{T}}$ and the sinusoidal solution $\left({I}_{_l}^{\mathrm{s}},{I}_r^{\mathrm{s}},{I}_u^{\mathrm{c}},{I}_v^{\mathrm{c}}\right)^{\mathrm{T}}$ of RTE at the same time. The stream of the double-Gaussian scheme is $2N = 16$.

Download

Larger image

PowerPoint slide

Fig  4.  Zenith distributions of I, Q, U, and V (solid line), along with their cosinusoidal (dashed line) and sinusoidal (dotted line) components, at the bottom boundary for L13 case when ϕ = 90° under a polarized source $\boldsymbol{S}\mathrm{_b}=\left(0.7\pi,0.3\pi,0.2\pi,0.05\pi\right)^{\mathrm{T}}$, with other parameters same as in Fig. 3.

Considering Eq. (42) with a polarized source, the output radiance should have both the cosinusoidal and the sinusoidal components of all four Stokes elements. It is due to the non-zero ${U}$ and ${V}$ of the source radiation, different from the L13 case in Section 3.2. The polarized $\boldsymbol{S}\mathrm{_b}$ in this case produces the complete cosinusoidal and sinusoidal components during the scattering process, which enables us to estimate the significance of the sinusoidal mode of RTE on the total radiance. As is shown in Fig. 4, the zenith distributions of the total Stokes elements all differ from those of their cosinusoidal or sinusoidal components. The variables ${I}$ and ${Q}\text{ }$, along with their cosinusoidal and sinusoidal components, increase first and then decrease with $\mu$ ranging from –1 to 0. However, it is obvious that the positions of their maximum differ from those in their cosinusoidal and sinusoidal components, especially for ${Q}$. This combined feature results from the different maximum positions of the cosinusoidal and sinusoidal components. Moreover, although ${U}$ and ${V}$ share the similar zenith trends with their sinusoidal components, their cosinusoidal components have the different zenith trends with their sinusoidal components and themselves. ${U}$ and its sinusoidal component decrease first and then increase with $\mu$, while ${V}$ and its sinusoidal component have the opposite trends. The cosinusoidal component of ${U}$ increases first and then decreases slowly. The cosinusoidal component of $V$ decreases with $\mu$ in $\left[ { - 1,0} \right]$. However, it must lead to the changes in the positions of their extremums due to their different characteristics of zenith distributions, which is obvious for all the four Stokes elements in Fig. 4. In addition, the total $V$ has significant differences in value with its sinusoidal component in the cosinusoidal mode of RTE. All these analyses imply that with the polarized source input, the sinusoidal mode equation should have significant impacts on the total radiation, especially for $V$ in this case.

Figure 4 proves that with the polarized source in radiative transfer, the sinusoidal solution of RTE is generated and cannot be ignored. The values of $I, Q,\; {\rm{and}}\; U$ agree with those in Schulz et al. (1999) and prove the effectiveness of updated VDISORT. Since the reflection of the earth’s surface in reality is polarized and some active radars with polarized emission are deployed (Yueh, 1997; Raney et al., 2021), RTMs based on a complete radiative transfer theory for the issues with the general boundary reflection, polarized sources or the general atmospheric scattering are required, especially for the case with non-0 elements in Eq. (42). The updated VDISORT in this paper can be an effective solver to complement the existing RTMs.

3.4   L11 case with random-oriented non-spherical particle scattering

To better validate the accuracy and the operation of the updated VDISORT, a case with random-oriented oblate particle scattering based on the Legendre expansion coefficients of 0th–11th orders of the phase matrix (hereinafter L11) is used for further investigations. The scattering matrix and phase matrix for random-oriented non-spherical particles is the same as Eqs. (40) and (41). The expansion coefficients were publically available in the model 2 of Kuik et al. (1992). The corresponding simulated Stokes radiance has been public in Siewert (2000), as the benchmark here for validation. Figure 5 shows the Stokes radiance at the bottom boundary for $\phi = 90^\circ$, with the double-Gaussian stream $2N = 16$ and $\boldsymbol{S}_{\mathrm{b}}=\left(0.5\pi, 0.5\pi,0,0\right)\mathrm{^T},\text{ }\omega=0.973527,\text{ }\lambda_0=0,\text{ }\mu_0=0.6,\text{ }\phi_0=0^{\circ},\text{ }\tau_L=1$. Also, no thermal radiance is considered in L11 case. Overall, the simulated radiance also matches well with the benchmark. Generally, $I$ and $Q{\text{ }}$ decreases with $\mu$ in $\left[ { - 1,0} \right]$. $U$ and $V$ share the same trend at the beginning but turn to increase around $\mu = - 0.6$ and $\mu = - 0.9$, respectively. Unlike $U$, $V$ turns back to decrease around $\mu = - 0.3$. The maximum of the absolute bias are less than 0.0008, 0.00005, 0.00006, 0.000003 for ${I},\ {Q},\ {U},\ {\rm{and}}\ {V}$, respectively. Similar to L13 case in Section 3.2, the Stokes radiance in L11 case with unpolarized $\boldsymbol{S}_{\mathrm{b}}$ and Lambertian reflection results from $\left({I}_{_l}^{\mathrm{c}},{I}_r^{\mathrm{c}},{I}_u^{\mathrm{s}},{I}_v^{\mathrm{s}}\right)\mathrm{^T}$, due to the similar symmetric relationship with Eq. (4) of Mie scattering and consequent Eq. (42). Although Figs. 2–5 can illustrate some zenith characteristics of the radiation, the most important thing is that they all prove the accuracy of the updated VDISORT model. The model can achieve similar simulations with the former versions and benchmarks.

Download

Larger image

PowerPoint slide

Fig  5.  Comparisons of the zenith distributions of I, Q, U, and V for L11 case when $\phi = 90^\circ$ with $\boldsymbol{S}_{\mathrm{b}}=\left(0.5\pi,0.5\pi,0,0\right)^{\mathrm{T}},\text{ }\omega=0.973527,$ λ₀ = 0, µ₀= 0.6, ϕ₀ = 0°, τ_L = 1, between the benchmarks (circles) and the simulations of the updated VDISORT (solid line).

4.   Summary and conclusions

This paper updates the discrete ordinate radiative transfer theory to better fit the general cases of atmospheric scattering and boundary reflection, and develops a new version of the discrete ordinate radiative transfer (VDISORT) model. Theoretically, the model does not need the assumption that the atmospheric scattering phase matrix and boundary reflection matrix should have the same azimuthal symmetric characteristics on the relative azimuthal distribution with Mie scattering phase matrix in Eq. (4). Meanwhile, it can simulate the complete Stokes information with both its sinusoidal and cosinusoidal modes. With two added zero-weighted Gaussian points, VDISORT can output the radiance at arbitrary zenith angle without interpolation. Through the complex variable and subroutines in Fortran codes, the complex eigenvalue problem is also taken into consideration.

Several tests illustrate that the updated VDISORT simulates the accurate results with the benchmarks for Rayleigh, L13 case with Mie scattering atmosphere for a group of spherical particles, L13 with a polarized source, and L11 case for a non-spherical particle scattering atmosphere, as the previous RTMs did. The maximum of the absolute bias for Rayleigh scattering case, unpolarized L13 and L11 are less than 0.0065, 0.0006, and 0.0008, respectively. With a polarized source input, it is clear that the sine-mode radiative transfer equation along with its Stokes information has obvious impacts on the total Stokes radiance, which proves the significance of VDISORT in preserving the full polarized radiation.

The main purpose of this paper is to update the discrete ordinate radiative transfer theory and to validate the accuracy of the updated VDISORT. The updated VDISORT does not require the properties of Mie scattering phase matrix on the atmospheric scattering and the boundary reflection. It can do more research with a pBRDF surface reflection without the azimuthal symmetric assumption, which will be tested in part 2 of this series of papers. It is believed that we have made a great step in radiative transfer simulation toward the oriented non-spherical particle scattering. With all advantages mentioned above, the updated VDISORT has more application potentials in radiative transfer sciences compared to the previous versions. The updated VDISORT will also be integrated into ARMS to support satellite data assimilation.