Separation of Rapidly-Varying and Slowly-Varying Processes and Development of Diffraction Decomposition Order Method in Radiative Transfer

Abstract

Single scattering in radiative transfer is separated into rapidly-varying and slowly-varying processes, where the rapidly-varying process (RVP) is mainly contributed by the diffraction effect. Accordingly, the diffraction decomposition order (DDO) method is developed to solve the vector radiative transfer equation (VRTE). Instead of directly solving the original VRTE, we decompose it into a series of order equations, where the zeroth-order equation replaces the RVP with a $\delta$-function while the high-order equations are the same as the zeroth-order one, except that the high decomposition orders of the RVP are used as driven sources. In this study, the DDO method is numerically realized using the successive order of the scattering method. The DDO is computationally efficient and accurate. More importantly, all physical processes in the VRTE are fully decomposed due to the order decomposition of the RVP and can be straightforwardly discussed.

1. Introduction

The radiative transfer process of electromagnetic energy is significant in atmospheric science, and is managed by the vector radiative transfer equation (VRTE) [1,2,3]. Many algorithms are developed to effectively solve the VRTE, such as the Discrete-Ordinate, adding-doubling, successive order of scattering, and Monte Carlo methods [4]. The radiative transfer methods are comprehensively discussed, and also compared to generate benchmark results (e.g., [5]).

Multiple scattering is an important process in radiative transfer. Scattering distribution in every scattering has a strong but steep forward peak due to the diffraction effect [6,7,8]. Direct calculation of the VRTE is usually time-consuming due to the strong diffraction peak. As a result, truncation techniques are developed to promote computational efficiencies, such as $\delta$-M, $\delta$-fit, and small-angle approximation methods [9,10,11]. Even though computational efficiency is significantly improved, computational accuracy is reduced due to the truncation with respect to the diffraction peak (e.g., [12,13,14]).

Other than diffraction, the ordinary scattering in every scattering is a slowly-varying process. Although the diffraction peak accounts for a large percentage of scattering energy, it is limited to narrow scattering angles and a rapidly-varying process, where an exemplary scattering phase function is shown in Figure 1. Based on the recognition, the diffraction peak can be decomposed into an order expansion, where the zeroth-order term is a $\delta$-function and the high-order terms are associated with the high-order differences between the diffraction peak and the $\delta$-function. Consequently, the original VRTE is decomposed into a series of order equations, called the diffraction decomposition order method (DDO). The rapidly-varying process is separated from the single scattering process and acts as driven sources for the high-order equations in the DDO method. The DDO method is efficient due to the direct $\delta$-function replacement of the diffraction peak at the zeroth-order equation and accurate due to the successive consideration of the diffraction peak at the high-order equations. Moreover, the processes of successive order of scattering and diffraction decomposition are both reflected and discussed in this study.

2. Method

Only the radiative transfer process in the solar radiation is considered since the diffraction peak is weakened in the longwave radiation. Defining $u=cos\theta$ and $\mu =\left|u\right|$, the VRTE in the plane-parallel assumption can be written as

$$u\frac{\partial \mathbf{I}(\Omega )}{\partial \tau }=\mathbf{I}(\Omega )-\frac{\omega }{4\pi }{\int \int }_{4\pi }\mathbf{P}(\Omega ,{\Omega }^{\prime })\mathbf{I}\left({\Omega }^{\prime }\right)d{\Omega }^{\prime },$$

(1)

and the boundary condition at the top of the atmosphere is

$$\left\{\begin{array}{cc}\hfill \mathbf{I}(\tau =0;\Omega )=& {\mathbf{E}}_0\delta (\Omega -{\Omega }_0),\hfill \\ \hfill & {\mathbf{E}}_0={\left(\begin{array}{cccc}{\mathrm{F}}_0/{\mu }_0& 0& 0& 0\end{array}\right)}^T,\hfill \end{array}\right.$$

(2)

where $\Omega$ and ${\Omega }_0$ denote the viewing $(u,\phi )$ and solar $(-{\mu }_0,{\phi }_0)$ solid angles, respectively; $d\Omega =dud\phi$; $\delta (\Omega -{\Omega }_0)=\delta (u+{\mu }_0)\delta (\phi -{\phi }_0)$; $\mathbf{I}$, $\mathbf{P}$, and $\omega$ are the Stokes vector, scattering phase matrix, and single scattering albedo, respectively; ${\mathrm{F}}_0$ is the solar irradiance. The argument $\tau$ of the physical quantities is suppressed in the equations. ${\Omega }_u$ and ${\Omega }_d$ are used to denote the upward $(\mu ,\phi )$ and downward $(-\mu ,\phi )$ solid angles.

In the DDO method, the scattering phase matrix $\mathbf{P}$ can be separated into the rapidly-varying process (RVP) (denoted by subscript ‘r’) and the slowly-varying process (SVP) (denoted by subscript ‘s’) as

$$\mathbf{P}=f{\mathbf{P}}_r+(1-f){\mathbf{P}}_s,$$

(3)

where f is the proportion factor of the RVP, and is usually a predefined parameter. We define the difference symbol between the RVP and the $\delta$-function as

$${\mathbf{P}}_{r\delta }={\mathbf{P}}_r-4\pi \mathrm{Diag}(1,1,1,1)\delta (\Omega -{\Omega }^{\prime }).$$

(4)

Equation (1) is reduced into the zeroth-order equation only associated with the SVP when ${\mathbf{P}}_{r\delta }=\mathbf{0}$ as

$$\left\{\begin{array}{c}u\frac{\partial {\mathbf{I}}_0(\Omega )}{\partial {\tau }_e}={\mathbf{I}}_0(\Omega )-\frac{{\omega }_e}{4\pi }{\int \int }_{4\pi }{\mathbf{P}}_s(\Omega ,{\Omega }^{\prime }){\mathbf{I}}_0\left({\Omega }^{\prime }\right)d{\Omega }^{\prime },\hfill \\ {\mathbf{I}}_0({\tau }_e=0;\Omega )={\mathbf{E}}_0\delta (\Omega -{\Omega }_0),\hfill \end{array}\right.$$

(5)

where the effective optical depth and single scattering albedo are defined as

$${\tau }_e=(1-\omega f)\tau ,{\omega }_e=\frac{\omega (1-f)}{1-\omega f}.$$

(6)

The residual equation of Equation (1) after the zeroth-order reduction is

$$\begin{array}{cc}\hfill u\frac{\partial {\mathbf{I}}_{res}(\Omega )}{\partial \tau }={\mathbf{I}}_{res}(\Omega )& -\frac{\omega }{4\pi }{\int \int }_{4\pi }\mathbf{P}(\Omega ,{\Omega }^{\prime }){\mathbf{I}}_{res}\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }\hfill \\ \hfill & -\frac{\omega f}{4\pi }{\int \int }_{4\pi }{\mathbf{P}}_{r\delta }(\Omega ,{\Omega }^{\prime }){\mathbf{I}}_0\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }.\hfill \end{array}$$

(7)

Repeating the same operations as Equation (1) reduced to Equations (5)–(7), Equation (7) is successively decomposed into a series of order equations, as follows:

$$\begin{array}{cc}\hfill u\frac{\partial {\mathbf{I}}_i(\Omega )}{\partial {\tau }_e}={\mathbf{I}}_i(\Omega )& -\frac{{\omega }_e}{4\pi }{\int \int }_{4\pi }{\mathbf{P}}_s(\Omega ,{\Omega }^{\prime }){\mathbf{I}}_i\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }\hfill \\ \hfill & -\frac{{\omega }_r}{4\pi }{\int \int }_{4\pi }{\mathbf{P}}_{r\delta }(\Omega ,{\Omega }^{\prime }){\mathbf{I}}_{i-1}\left({\Omega }^{\prime }\right)d{\Omega }^{\prime },\hfill \end{array}$$

(8)

where $i=1,\dots ,\infty$, and the single scattering albedo associated with the RVP is

$${\omega }_r=\frac{\omega f}{1-\omega f}.$$

(9)

The solution of Equation (5) is physically composed of direct and diffuse terms

$${\mathbf{I}}_0={\mathbf{I}}_{0,dir}+{\mathbf{I}}_{0,dif}=exp\left(-\frac{{\tau }_e}{{\mu }_0}\right){\mathbf{E}}_0\delta (\Omega -{\Omega }_0)+{\mathbf{I}}_{0,dif}.$$

(10)

Substituting Equation (10) into Equation (5), the zeroth-order equation is simplified as

$$\begin{array}{cc}\hfill u\frac{\partial {\mathbf{I}}_{0,dif}(\Omega )}{\partial {\tau }_e}={\mathbf{I}}_{0,dif}(\Omega )& -\frac{{\omega }_e}{4\pi }{\int \int }_{4\pi }{\mathbf{P}}_s(\Omega ,{\Omega }^{\prime }){\mathbf{I}}_{0,dif}\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }\hfill \\ \hfill & -\frac{{\omega }_e}{4\pi }exp\left(-\frac{{\tau }_e}{{\mu }_0}\right){\mathbf{P}}_s(\Omega ,{\Omega }_0){\mathbf{E}}_0.\hfill \end{array}$$

(11)

For consistency, the notation ${\mathbf{I}}_0$ is used to replace the diffuse term ${\mathbf{I}}_{0,dif}$ hereafter.

Substituting Equation (10) into Equation (8) at $i=1$, Equations (8), (10) and (11) can be consistently organized as follows:

$$u\frac{\partial {\mathbf{I}}_i(\Omega )}{\partial {\tau }_e}={\mathbf{I}}_i(\Omega )-\frac{{\omega }_e}{4\pi }{\int \int }_{4\pi }{\mathbf{P}}_s(\Omega ,{\Omega }^{\prime }){\mathbf{I}}_i\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }-{\mathbf{S}}_i(\Omega ),$$

(12)

where $i=0,\dots ,\infty$, and the driven source functions ${\mathbf{S}}_i$ are explicitly represented as

$$\left\{\begin{array}{cc}\hfill {\mathbf{S}}_0(\Omega )=& \frac{{\omega }_e}{4\pi }exp\left(-\frac{{\tau }_e}{{\mu }_0}\right){\mathbf{P}}_s(\Omega ,{\Omega }_0){\mathbf{E}}_0,\hfill \\ \hfill {\mathbf{S}}_1(\Omega )=& \frac{{\omega }_r}{4\pi }{\int \int }_{4\pi }{\mathbf{P}}_{r\delta }(\Omega ,{\Omega }^{\prime }){\mathbf{I}}_0\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }+{\mathbf{S}}_{1\delta }(\Omega ),\hfill \\ \hfill {\mathbf{S}}_i(\Omega )=& \frac{{\omega }_r}{4\pi }{\int \int }_{4\pi }{\mathbf{P}}_{r\delta }(\Omega ,{\Omega }^{\prime }){\mathbf{I}}_{i-1}\left({\Omega }^{\prime }\right)d{\Omega }^{\prime },i=2,\dots ,\infty ,\hfill \\ \hfill {\mathbf{S}}_{1\delta }(\Omega )=& \frac{{\omega }_r}{4\pi }exp\left(-\frac{\tau }{{\mu }_0}\right){\mathbf{P}}_r(\Omega ,{\Omega }_0){\mathbf{E}}_0,\hfill \\ \hfill -& \frac{{\omega }_e}{4\pi }\left[exp\left(-\frac{{\tau }_e}{{\mu }_0}\right)-exp\left(-\frac{\tau }{{\mu }_0}\right)\right]{\mathbf{P}}_s(\Omega ,{\Omega }_0){\mathbf{E}}_0,\hfill \end{array}\right.$$

(13)

and the corresponding direct term is

$${\mathbf{I}}_{dir}=exp\left(-\frac{\tau }{{\mu }_0}\right){\mathbf{E}}_0\delta (\Omega -{\Omega }_0).$$

(14)

The solution of the VRTE in Equation (1) can be obtained by

$$\mathbf{I}={\mathbf{I}}_{dir}+\sum _{i=0}^{\infty }{\mathbf{I}}_i.$$

(15)

The source function ${\mathbf{S}}_{1\delta }$ in the first-order equation is explicitly from the direct term of the zeroth-order equation of Equation (10). It is reduced to a $\delta$-function if the RVP is reduced to a $\delta$-function, that is, when ${\mathbf{P}}_{r\delta }=\mathbf{0}$, ${\mathbf{S}}_{1\delta }=\left[exp\left(-{\tau }_e/{\mu }_0\right)-exp\left(-\tau /{\mu }_0\right)\right]{\mathbf{E}}_0\delta (\Omega -{\Omega }_0)$, which is exactly the direct term difference between Equations (10) and (15). The RVP stays in the driven source functions of the high-order equations while the SVP is involved in the multiple scattering processes, and the two processes are fully separated. The contribution of the diffraction effect is successively considered in the high-order equations. Equations (12)–(15) are the specific expressions of the DDO method.

The total optical depth is denoted as b and the effective total optical depth is $b_e=(1-\omega f)b$, as in Equation (6). Formally, Equation (12) can be written as upward and downward components:

$$\begin{array}{cc}\hfill {\mathbf{I}}_i({\tau }_e;{\Omega }_u)=& {\Gamma }_{\mu }\left({\tau }_e,b_e\right){\mathbf{I}}_i(b_e;{\Omega }_u)\hfill \\ \hfill & +{\int }_{{\tau }_e}^{b_e}\frac{d{\tau }_e^{\prime }}{\mu }{\Gamma }_{\mu }\left({\tau }_e,{\tau }_e^{\prime }\right)\frac{{\omega }_e}{4\pi }{\int \int }_{4\pi }d{\Omega }^{\prime }{\mathbf{P}}_s({\tau }_e^{\prime };{\Omega }_u,{\Omega }^{\prime }){\mathbf{I}}_i({\tau }_e^{\prime };{\Omega }^{\prime })\hfill \\ \hfill & +{\int }_{{\tau }_e}^{b_e}\frac{d{\tau }_e^{\prime }}{\mu }{\Gamma }_{\mu }\left({\tau }_e,{\tau }_e^{\prime }\right){\mathbf{S}}_i({\tau }_e^{\prime };{\Omega }_u),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\mathbf{I}}_i({\tau }_e;{\Omega }_d)=& {\Gamma }_{\mu }\left(0,{\tau }_e\right){\mathbf{I}}_i(0;{\Omega }_d)\hfill \\ \hfill & +{\int }_0^{{\tau }_e}\frac{d{\tau }_e^{\prime }}{\mu }{\Gamma }_{\mu }\left({\tau }_e^{\prime },{\tau }_e\right)\frac{{\omega }_e}{4\pi }{\int \int }_{4\pi }d{\Omega }^{\prime }{\mathbf{P}}_s({\tau }_e^{\prime };{\Omega }_d,{\Omega }^{\prime }){\mathbf{I}}_i({\tau }_e^{\prime };{\Omega }^{\prime })\hfill \\ \hfill & +{\int }_0^{{\tau }_e}\frac{d{\tau }_e^{\prime }}{\mu }{\Gamma }_{\mu }\left({\tau }_e^{\prime },{\tau }_e\right){\mathbf{S}}_i({\tau }_e^{\prime };{\Omega }_d),\hfill \end{array}$$

(17)

where the optical depth dependence is explicitly given for clarity and the transmittance function ${\Gamma }_{\mu }$ is defined as

$${\Gamma }_{\mu }\left({\tau }_1,{\tau }_2\right)=exp\left(-\frac{{\tau }_2-{\tau }_1}{\mu }\right).$$

(18)

Equations (16) and (17) can be straightforwardly solved using the successive order of the scattering (SOS) process. One can define the multiple scattering term as

$${\mathbf{M}}_{ij}({\tau }_e;\Omega )=\frac{{\omega }_e}{4\pi }{\int \int }_{4\pi }d{\Omega }^{\prime }{\mathbf{P}}_s({\tau }_e;\Omega ,{\Omega }^{\prime }){\mathbf{I}}_{i(j-1)}({\tau }_e;{\Omega }^{\prime }),$$

(19)

where j enumerates the scattering order and $j=2,\dots ,\infty$. Double integration with respect to solid angles in the equations is reduced to only zenith integration after using Fourier order expansion associated with azimuthal dependence, which is described in details [15]. For simplification, the multiple and single scattering notations are merged into one notation as

$${\mathbf{L}}_{ij}({\tau }_e;\Omega )=\left\{\begin{array}{cc}\hfill {\mathbf{S}}_i({\tau }_e;\Omega ),& j=1,\hfill \\ \hfill {\mathbf{M}}_{ij}({\tau }_e;\Omega ),& j=2,\dots ,\infty .\hfill \end{array}\right.$$

(20)

Correspondingly, the upward Stokes vector in Equation (16) can be expanded as successive orders

$${\mathbf{I}}_{ij}({\tau }_e;{\Omega }_u)={\int }_{{\tau }_e}^{b_e}{\Gamma }_{\mu }\left({\tau }_e,{\tau }_e^{\prime }\right){\mathbf{L}}_{ij}({\tau }_e^{\prime };{\Omega }_u)\frac{d{\tau }_e^{\prime }}{\mu },$$

(21)

$${\mathbf{I}}_i({\tau }_e;{\Omega }_u)={\Gamma }_{\mu }\left({\tau }_e,b_e\right){\mathbf{I}}_i(b_e;{\Omega }_u)+\sum _{j=1}^{\infty }{\mathbf{I}}_{ij}({\tau }_e;{\Omega }_u).$$

(22)

Similarly, the downward Stokes vector in Equation (17) can be expanded as

$${\mathbf{I}}_{ij}({\tau }_e;{\Omega }_d)={\int }_0^{{\tau }_e}{\Gamma }_{\mu }\left({\tau }_e^{\prime },{\tau }_e\right){\mathbf{L}}_{ij}({\tau }_e^{\prime };{\Omega }_d)\frac{d{\tau }_e^{\prime }}{\mu },$$

(23)

$${\mathbf{I}}_i({\tau }_e;{\Omega }_d)={\Gamma }_{\mu }\left(0,{\tau }_e\right){\mathbf{I}}_i(0;{\Omega }_d)+\sum _{j=1}^{\infty }{\mathbf{I}}_{ij}({\tau }_e;{\Omega }_d).$$

(24)

The flow-chart of the diffraction decomposition order algorithm using the successive order of scattering method is shown in Figure 2. Each block represents the process solving the corresponding order equation using the SOS method.

3. Numerical Results

The benchmark results are given by Kokhanovsky et al. [5]. Figure 3 shows the reflected radiance and corresponding relative difference in the case of a predefined aerosol layer. The results show excellent consistency with the benchmark, with errors within 0.5%. Large errors are shown when the viewing zenith angles are close to 90° because the effective optical depth is close to infinity. The zenith angles and azimuthal angles for the solar, upward and downward directions are illustrated in Figure 4.

The scattering phase matrix, and corresponding rapidly-varying and slowly-varying components in Figure 5 are used to exemplify the DDO method. The rapidly-varying scattering phase matrix is truncated at ${4.5}^{\circ }$ and the corresponding proportion factor f is $0.452$ by calculation. The scattering phase function $P_{11}$ from $0^{\circ }$ to ${4.5}^{\circ }$ is further shown in the inlet. For the rapidly-varying component ${\mathbf{P}}_r$, the elements $P_{22}$, $P_{33}$, and $P_{44}$ are extremely close to $P_{11}$ and the elements $P_{12}$ and $P_{34}$ are close to 0, which confirm that the zeroth-order of ${\mathbf{P}}_{r\delta }$ equal to 0 is reasonable. Only one homogeneous scattering layer and black surface are used in the following results to discuss the DDO method, that is, ${\mathbf{I}}_i(b_e;{\Omega }_u)=0$ in Equation (22) and ${\mathbf{I}}_i(0;{\Omega }_d)=0$ in Equation (24). The total optical depth b and single-scattering albedo $\omega$ are 5 and 1, respectively. The upright direction is set to be the positive z-direction for the zenith angle and the projection direction of solar incidence is set to be the positive x-direction for the azimuthal angle since only relative azimuthal angle (RAA) $\Delta \phi =(\phi -{\phi }_0)$ matters. The solar zenith angle (SZA) is fixed to be ${120}^{\circ }$, that is, ${\mu }_0=0.5,{\phi }_0=0^{\circ }$. All scattering results are calculated by the straightforward successive order of scattering method (denoted as SOS in the figure legends) and the DDO method with different orders.

Stokes vectors are highly influenced by the RVP when the viewing azimuthal angle (VAA) aligns with the solar azimuthal angle. The upward and downward Stokes vectors and their differences with $\phi =0^{\circ }$ at the layer top, middle, and bottom are shown in Figure 6, respectively. The legend DDO00 denotes the results from the zeroth-order equation, while the other DDO legends denote results from accumulated high-order contributions. The viewing zenith angles (VZA) are all from $0^{\circ }$ to ${180}^{\circ }$, but only non-zero parts are plotted in the figures as the downward ones in Figure 6a–d and the upward ones in Figure 6i–l. Each scattering of the RVP is limited to small angles. Consequently, the influence on upward directions is much less than on downward directions, and the influence on small upward VZAs is much less than on large upward VZAs. The VZAs close to the solar zenith angle are sensitive to the RVP, as shown in Figure 6e–l. The differences between the original and zeroth-order equations can be gradually reduced by the successive consideration of the diffraction orders. The order number reflects the degree of influence from the RVP.

The calculation process is almost the same if the VZAs are the same since the azimuthal dependence is expanded using Fourier order. The final results are obtained by summarizing the Fourier orders using specific VAAs. The results at layer middle are shown in Figure 7 when the VAA is ${60}^{\circ }$. The viewing directions are far away from the direction of solar incidence so that the influence from the RVP is small, which can be verified from the zeroth-order results. Only one diffraction order is enough to improve the results. However, many Fourier orders are necessary to mutually cancel the effect from the RVP when we compare the radiance values in Figure 6e–h and Figure 7.

The difference features at each decomposition order are shown in Figure 8.The differences between the original and zeroth-order equations are apparent in the direction of solar incidence. As the diffraction order increases, the difference can be gradually reduced. The order number reflects the extent of influence from the RVP.

4. Conclusions

In conclusion, the diffraction decomposition order method is developed in this study by separating the rapidly-varying process from a single scattering process. In the DDO method, the original radiative transfer equation is expanded into a series of order equations, where the zeroth-order equation is fully composed of the SVP and the RVP only stays in the driven source of the high-order equations. The RVP and SVP are completely separated in the DDO method, and the influence from the RVP can be iteratively considered through high-order equations. The influence from the RVP is limited to small scattering angles so the contribution from the RVP is predictable. The radiative transfer process can be efficiently and accurately reproduced by using different diffraction orders for different degrees of influence.