A New Parameterization of Canopy Radiative Transfer for Land Surface Radiation Models

Solar radiative transfer (RT) is a process that takes place at the land surface that is important in determining energy, water, moisture vapor and carbon balances, especially for areas covered by a vegetation canopy (Yu et al., 2016). In addition, how to deal with canopy RT is a key issue in remote sensing for the retrieval of land surface information. The RT through the canopy is more complicated than that through the atmosphere. The reflectance of the canopy depends on its structure, the optical properties of the foliage, the viewing geometry, and the soil background reflectance (Ross, 1981). The solution of the canopy RT equation requires knowledge of the phase function as an input. In reality, the leaf phase function is complicated and anisotropic (Biswas, 2007). Though there have been studies on leaf phase function (Biswas, 2007; Otto and Trautmann, 2008), the anisotropic leaf phase function, which depends on the leaf normal distribution and leaf reflection/transmission, has not yet been derived and parameterized. To solve the parameterization problem of leaf phase function and apply it to the canopy RT process is the primary aim of this study.

The RT equation is an integro-differential equation, and thus it is difficult to obtain the exact solution. Therefore, several approximation methods have been proposed for dealing with the atmospheric RT (Liou, 1974; Joseph et al., 1976; Kylling et al., 1995; Nakajima et al., 2000; Lu et al., 2009; Zhang and Li, 2013; Zhang et al., 2013b; Wang et al., 2017). For climate modeling, the widely used methods for RT in the atmosphere are the Eddington approximation, the two-stream discrete ordinate method, and the hemispheric constant method (Coakley and Chylek, 1975). In previous studies on the canopy radiation process (Dickinson, 1983; Sellers, 1985), the hemispheric constant method has been used to obtain the canopy reflectance, transmittance and absorption in canopy models. However, through comparisons with the atmospheric RT method (Joseph et al., 1976; Lu et al., 2009; Zhang and Li, 2013), the Eddington approximation has been found to be more suitable than the hemispheric constant method across a much wider range of optical depths. By using the newly proposed parameterization on leaf phase function, the derived analytical solution for the Eddington approximation is applicable for most canopy radiation problems. To investigate the Eddington approximation for canopy RT is the secondary aim of this study.

In previous studies, the canopy diffuse and direct albedo have been calculated separately (Dickinson, 1983; Sellers, 1985; Dai and Sun, 2006, 2007a, 2007b). The relationship between the diffuse and direct canopy radiation is also discussed in this study.

Following this introduction, in section 2, the parameterization of leaf phase function is proposed and the Eddington approximation for the canopy radiation process is discussed. In section 3, the accuracy of the new method in terms of reflection/transmission is examined through comparison with the hemispheric constant method. A short summary is provided in section 4.

The azimuthally averaged RT equation of the canopy in general form is (e.g., Myneni et al., 1989) \begin{equation} \label{eq1} \mu\frac{dI(L,\mu)}{dL}=I(L,\mu)G(\mu)-\frac{\omega}{2}\int_{-1}^1I(L,\mu')P(\mu,\mu')G(\mu')d\mu' , (1)\end{equation} where I, μ, ω, L and P(μ,μ') are the diffuse intensity, the cosine of the local zenith angle, the single scattering albedo, the cumulative leaf area index, and the azimuthally averaged scattering phase function defining the light incidence at μ' and scattered away at μ, respectively. A canopy is confined between depth zero (L=0) at the top and L=L₀ at the bottom, where L₀ is the total leaf area index. G(μ) is the geometric factor, which can be defined in terms of the leaf normal distribution function \(\hat{g}_L(\mu_L)\), with μ_L=cos(θ_L) and θ_L being the zenith angle of the orientation of the leaf. \(\hat{g}_L(\mu_L)\) a leaf-distribution function, which represents the probability that a single leaf has a normal in the direction about μ_L. It satisfies the following normalization condition: \begin{equation} \label{eq2} \int_0^1\hat{g}_L(\mu_L)d\mu_L=1 . (2)\end{equation} Besides, the azimuthally averaged G(μ) is defined by (Ross, 1981). In many analytic leaf angle distribution models, \(\hat{g}_L(\mu_L)\) and G(μ) have been discussed in detail (Shultis and Myneni, 1988; Biswas, 2007; Otto and Trautmann, 2008). It is well known that G(μ) has a special symmetric property, as G(μ)=G(-μ) (Otto and Trautmann, 2008; Picca and Furfaro, 2013). In this work, we mainly focus on the parameterization of phase function and the solution of the canopy RT equation.

The most difficult part in the calculation of canopy radiation is the phase function for a canopy medium. Using Legendre function expansion, the phase function P(cosΘ) can be written as \begin{equation} \label{eq3} P(\cos\Theta)=\sum_{l=0}^N\omega_lP_l(\cos\Theta) , (3)\end{equation} where P_l is the Legendre function. According to the orthogonal property of Legendre polynomials, we obtain \begin{equation} \omega_l=\frac{2l+1}{2}\int_{-1}^1P(\cos\Theta)P_l(\cos\Theta)d\cos\Theta , (4)\end{equation} where l=0,1,…,N and the cosine of the scattering angle is defined by cosΘ=μμ'+(1-μ²)^1/2(1-μ'²)^1/2cos(φ-φ'), with φ being the azimuthal angle. It defines the light at (μ,φ), and that which scatters away at (μ',φ'). In addition, ω₀ is equal to 1, and ω₁ is equal to 3g, with g being the asymmetry factor. For spherically distributed leaves, the phase function becomes (Ross, 1981) \begin{equation} \label{eq4} P_{\rm sp}(\cos\Theta)=\frac{8}{3\pi}(\sin\Theta-\Theta\cos\Theta)+\frac{8\sigma}{\pi}\cos\Theta , (5)\end{equation} where $$ \sigma=\frac{t_L}{\omega}, $$ and t_L is the leaf transmittance. Based on the conservation of energy, the single albedo ω is equal to r_L+t_L, where r_L is the leaf reflectance. Therefore, the asymmetry factor of spherically distributed leaves is \begin{equation} \label{eq5} g_{\rm sp}=\frac{1}{2}\int_{-1}^1P_{\rm sp}(\cos\Theta)\cos\Theta d\cos\Theta=\frac{8\sigma}{3\pi}-\frac{4}{9} . (6)\end{equation} For arbitrarily distributed leaves, the phase function can not be expressed analytically. In most previous studies on atmospheric RT, the Henyey-Greenstein (HG) scheme has been used (Henyey and Greenstein, 1941). The HG scheme is used in this study to represent the phase function of the canopy: \begin{equation} \label{eq6} P(\cos\Theta)=\frac{1-g^2}{(1+g^2-2g\cos\Theta)^{\frac{3}{2}}} . (7)\end{equation} The asymmetry factor g depends on the leaf scattering characteristics. Therefore, g is difficult to obtain. In the following, we attempt to build a new scheme for the asymmetry factor g. As inferred from the analysis in the appendix of (Norman and Jarvis, 1975), the ratio of the backscattering energy to the total scattering energy is approximated as \begin{equation} \label{eq7} \beta=\frac{1}{2}[1+(1-2\sigma)\langle{\cos^2\theta_L}\rangle] , (8)\end{equation} where \(\langle\cos^2\theta_L\rangle\) indicates a mean of cos²θ_L. In two-stream approximation (Meador and Weaver, 1980), β can be expressed as \begin{equation} \label{eq8} \beta=\frac{1-g}{2} .(9) \end{equation} From Eqs. (8) and (9), \begin{equation} \label{eq9} g=(2\sigma-1)\langle{\cos^2\theta_L}\rangle . (10)\end{equation} For spherically distributed leaves, cos²θ_L is equal to 1/3 (Ross, 1981; De Ridder, 2001). For arbitrarily distributed leaves, the leaf orientation \(\langle\cos^2\theta_L\rangle\) can be parameterized by the leaf inclination index χ_L (De Ridder, 2001), \begin{equation} \label{eq10} \langle{\cos^2\theta_L}\rangle=\left\{ \begin{array}{l@{\quad}l} \dfrac{1}{3}(1+\chi_L) & (\chi_L<0)\\[3mm] \dfrac{1}{3}(1+2\chi_L) & (\chi_L>0) \end{array} \right. . (11)\end{equation} Substituting Eq. (11) into Eq. (10), we finally obtain \begin{equation} \label{eq11} g=\left\{ \begin{array}{l@{\quad}l} \dfrac{1}{3}(2\sigma-1)(1+\chi_L) & (\chi_L<0)\\[3mm] \dfrac{1}{3}(2\sigma-1)(1+2\chi_L) & (\chi_L>0) \end{array} \right. . (12)\end{equation} χ_L is defined by (Ross, 1981): \begin{equation} \label{eq12} \chi_L=\pm\frac{1}{2}\int_0^1|1-\hat{g}(\mu_L)|d\mu_L . (13)\end{equation} A purely horizontal (vertical) leaf angle distribution function (LADF) is defined in χ_L=+1(χ_L=-1), and a uniform LADF is defined in χ_L=0. In the new scheme, g is parameterized by σ and the leaf inclination index χ_L. Figure 1 shows that the value of g varies with σ and χ_L. When \(\sigma\to 1\) and \(\chi_L\to 1\), it yields \(g\to 1\). When \(\sigma\to 1\) and \(\chi_L\to 0\), it yields \(g\to -1\). For isotropic scatters of leaves, we obtain σ =0.5 and g=0.

A simple parameterization of g was also obtained by (Zhou et al., 2009) (hereafter referred to as the "Zhou scheme"): \begin{equation} \label{eq13} g=\frac{1}{2}[1+(1-2\sigma)\langle{\cos^2\theta_L}\rangle] . (14)\end{equation} Figure 2 compares the accuracy of our new scheme with that of the Zhou scheme for spherically distributed leaves. The benchmark result is obtained by Eq. (6). Figure 2 shows that the asymmetry factor obtained by the new scheme is close to the benchmark result, but the Zhou scheme and benchmark calculation differ quite considerably. The result in Fig. 2 indicates that the new scheme can represent the leave phase function properly.

The derived parameterization for phase function can be applied to the two-stream Eddington approximation and higher-order stream RT methods. In the following, only the detailed solution for the Eddington approximation is shown. The extension to a higher-order RT method will be discussed in a subsequent study.

Figure 1.  Asymmetry factor g as a function of single scattering albedo σ and leaf inclination index χ_L.

Figure 2.  Comparison of asymmetry factor g in different schemes, where σ is the single scattering albedo.

Eddington approximation is widely used in atmospheric radiation (Joseph et al., 1976; Li and Ramaswamy, 1996). The RT through the canopy, which considers canopy structure and leaf orientation, is more complicated than the RT through the atmosphere. In the following, the Eddington approximation is applied to the canopy RT for both the direct and diffuse radiation. Furthermore, the relationship between the direct and diffuse radiation is illustrated. It is shown that the direct reflection/transmission can be a function of the variables from the diffuse radiation, which is a new result for the canopy RT.

The incoming diffuse radiance is assumed to be isotropic and the incoming flux is set as1. In addition, no upward diffused flux at the bottom of the canopy is assumed. For diffuse radiation, the canopy RT equation is: \begin{align} \mu\dfrac{dI(L,\mu)}{dL}&=I(L,\mu)G(\mu)-\frac{\omega}{2}\int_{-1}^1I(L,\mu')P(\mu,\mu')G(\mu')d\mu' ;(15a)\\ I(0,\mu)&=\frac{1}{\pi}\quad (-1\le\mu\le 0) ;(15b)\\ I(L_0,\mu)&=0\quad (0\le\mu \le 1) .(15c) \end{align} The upward/downward fluxes are obtained as: \begin{align} F^+(L)&=2\pi\int_0^1I(L,\mu)\mu d\mu ;(16a)\\ F^-(L)&=2\pi\int_0^{-1}I(L,\mu)\mu d\mu . (16b)\end{align} Similarly, the boundary conditions of Eqs. (15b)-(15c) can be written in the form of flux. By multiplying 2π dμ on both sides of Eq. (15) and performing two integrals with intervals of [0,1] and [-1,0], we obtain:

\begin{align} \frac{dF^+}{dL}&=2\pi\int_0^1G(\mu)I(L,\mu)d\mu\nonumber\\ &\quad-\pi\omega\int_0^1\int_{-1}^1G(\mu')I(L,\mu')P(\mu,\mu')d\mu'd\mu ;(17a)\\ -\frac{dF^-}{dL}&=2\pi\int_{-1}^0G(\mu)I(L,\mu)d\mu\nonumber\\ &\quad-\pi\omega\int_{-1}^0\int_{-1}^1G(\mu')I(L,\mu')P(\mu,\mu')d\mu'd\mu ;\ \ (17b) \\ F^-(L&=0)=1 ;(17c)\\ F^+(L&=L_0)=0 . (17d)\end{align} By Eddington approximation, we obtain \begin{equation} \label{eq14} I(L,\mu)=I_0(L)+I_1(L)\mu .(18) \end{equation} Substituting Eq. (18) into Eq. (17), we get: \begin{align} \dfrac{d^2F^+}{dL^2}+(\gamma_4-\gamma_1)\frac{dF^+}{dL}+(\gamma_2\gamma_3-\gamma_4\gamma_1)F^+=0 ; (19a)\\ \dfrac{d^2F^-}{dL^2}+(\gamma_4-\gamma_1)\frac{dF^-}{dL}+(\gamma_2\gamma_3-\gamma_4\gamma_1)F^-=0 ; (19b)\end{align} where γ₁,γ₂,γ₃,γ₄ are listed in the appendix. The corresponding solution of Eq. (16) becomes: \beginsubequations \begin{align} F^+(L)&=A_1{\rm e}^{k_1L}+B_1{\rm e}^{k_2L} ;(20a)\\ F^-(L)&=A_1a{\rm e}^{k_1L}+B_1b{\rm e}^{k_2L} ;(20b) \end{align} where \begin{eqnarray*} k_1&=&\frac{1}{2}\left[\gamma_1-\gamma_4+\sqrt{(\gamma_1+\gamma_4)^2-4\gamma_2\gamma_3}\right] ,\\ k_2&=&\frac{1}{2}\left[\gamma_1-\gamma_4-\sqrt{(\gamma_1+\gamma_4)^2-4\gamma_2\gamma_3}\right] , \end{eqnarray*} a=γ₃/(k₁+γ₄), and b=γ₃/(k₂+γ₄). According to the boundary conditions of Eqs. (17c)-(17d), we get A₁=1/[a-be^{-(k₂-k₁)L₀}] and B₁=1/[b-ae^{(k₂-k₁)L₀}]. Therefore, the reflection/transmission of the diffuse radiation is: \begin{align} \bar{r}&=F^+(0)=\frac{e^{k_2 L_0}-e^{k_1 L_0}}{ae^{k_2L_0}-be^{k_1L_0}} ;(21a)\\ \bar{t}&=F^-(L_0)=\frac{(a-b)e^{(k_1+k_2)L_0}}{ae^{k_2L_0}-be^{k_1L_0}} . (21b)\end{align} \endsubequations For direct radiation, there is a direct solar beam at the top of the considered canopy layer and no upward diffuse radiance from the bottom of the canopy layer. The equation of direct radiation is: \begin{align} \mu\dfrac{dI(L,\mu)}{dL}&=G(\mu)I(L,\mu)-\dfrac{\omega}{2}\int_{-1}^1I(L,\mu')P(\mu,\mu')G(\mu')d\mu' ;(22a)\\ I(0,\mu)&=\delta(\mu,-\mu _0)\frac{F_0}{2\pi}\quad (-1\le\mu \le 0) ;(22b)\\ I(L_0,\mu)&=0\quad (0\le\mu\le 1) ; (22c)\end{align} where μ₀ is the cosine of the solar zenith angle and δ is the delta function. The intensity I(L,μ) can be decomposed into the direct solar beam and diffuse beam (Zhang et al., 2013a): \begin{equation} \label{eq15} I(L,\mu)=I_{\rm dir}(L,\mu)+I_{\rm dif}(L,\mu) . (23)\end{equation} The direct solar beam I _dir(L,μ) and the diffuse beam I _dif(L,μ) are photons without experiencing scattering and photons experiencing scattering, respectively. Therefore, we have: \begin{align} \mu\dfrac{dI_{\rm dir}(L,\mu)}{dL}&=G(\mu)I_{\rm dir}(L,\mu) ; (24a)\\ I_{\rm dir}(0,\mu)&=\delta(\mu,-\mu_0)\frac{F_0}{2\pi}\quad (-1\le\mu\le 0) ;(24b) \\ I_{\rm dir}(L_0,\mu)&=0\quad (0\le\mu \le 1) . (24c)\end{align} The solution of Eq. (24) is \begin{equation} \label{eq16} I_{\rm dir}(L,\mu)=\delta(\mu,-\mu_0)\frac{F_0}{2\pi}e^{LG(\mu)/\mu}\quad (-1\le \mu \le 0) .(25) \end{equation} Substituting Eqs. (23) and (25) into Eq. (22), we get: \beginsubequations \begin{align} \mu\frac{dI_{\rm dif}(\tau,\mu)}{d\tau}&=\!G(\mu)I_{\rm dif}(\tau,\mu)-\!\frac{\omega}{2}\int_{-1}^1I_{\rm dif}(\tau,\mu')P(\mu,\mu')G(\mu')d\mu'\nonumber\\ &\quad-\frac{\omega F_0}{4\pi}e^{-L/\mu_0}G(\mu_0)P(\mu,-\mu_0) ; (26a)\\ I_{\rm dif}(0,\mu)&=\!0\quad (-1\le\mu\le 0) ;(26b)\\ I_{\rm dif}(L_0,\mu)&=\!0\quad (0\le\mu\le 1) . (26c)\end{align} \endsubequations According to the Eddington approximation, Eq. (26) can be written in the format of fluxes F⁺ and F^-: \begin{align} &\frac{dF_{\rm dif}^+}{dL}=\gamma_1F_{\rm dif}^+-\gamma_2F_{\rm dif}^-\omega G(\mu_0)F_0e^{-LG(\mu_0)/\mu_0}\beta_0 ; (27a)\\ &\frac{dF_{\rm dif}^-}{dL}=\gamma_3F_{\rm dif}^+-\gamma_4F_{\rm dif}^-+\omega G(\mu_0)F_0e^{-LG(\mu_0)/\mu_0}(1-\beta_0) ;(27b) \\ &F_{\rm dif}^+(0)=0 ;(27c)\\ &F_{\rm dif}^-(L_0)=0 ; (27d)\end{align} where $$ \beta_0=\int_0^1\frac{P(\mu,-\mu_0)}{2}d\mu . $$ The solution of Eq. (27) is similar to Eq. (20). Consequently, we obtain: \begin{align} F_{\rm dif}^+(L)&=A_2e^{k_1L}+B_2e^{k_2L}+F_0\varsigma e^{-LG(\mu_0)/\mu_0} ;(28a)\\ F_{\rm dif}^-(L)&=A_2ae^{k_1L}+B_2be^{k_2L}+F_0\eta e^{-LG(\mu_0)/\mu_0} ; (28b)\end{align} where \(\varsigma=\omega G(\mu_0)[\gamma_2(\beta_0-1)+d\beta_0]/(cd+\gamma_2\gamma_3)\), η=ω G(μ₀)[c(β₀-1)-γ₃β₀]/(cd+γ₂γ₃), c=G(μ₀)/μ₀+γ₁, and d=G(μ₀)/μ₀-γ₄. According to the boundary conditions of Eqs. (27c) and (27d), we obtain \(A_2=F_0[b\varsigma e^-L_0G(\mu_0)/\mu_0-\) η e^k₂L₀]/(ae^k₂L₀-be^k₁L₀) and \(B_2\!=\!F_0[\eta e^k_1 L_0-a\varsigma e^-L_0G(\mu_0)/\mu_0]/\) (ae^k₂L₀-be^k₁L₀). Therefore, the direct reflection and transmission are: \begin{align} r(\mu_0)&=\frac{F_{\rm dif}^+(0)}{\mu_0F_0}\!=\!\frac{\varsigma(b-a)e^{-L_0G(\mu_0)/\mu_0}+\varsigma(e^{k_1 L_0}-e^{k_2 L_0})}{(ae^{k_2 L_0}-be^{k_1L_0})\mu_0}+\frac{\eta}{\mu_0} ;(29a)\\ t(\mu_0)&=\frac{F_{\rm dif}^-(L_0)}{\mu_0F_0}=\frac{\eta}{\mu_0}e^{-L_0G(\mu_0)/\mu_0}\nonumber\\ &\quad+\frac{\eta(b-a)e^{(k_1+k_2)L_0}+ab\varsigma(e^{k_1L_0}-e^{k_2L_0})e^{-L_0G(\mu_0)/\mu_0}}{(ae^{k_2 L_0}-be^{k_1L_0})\mu_0} .(29b) \end{align} Based on Eq. (21), the direct radiation, Eq. (29) can be represented by the function of the diffuse radiation: \begin{align} r(\mu_0)&=\frac{1}{\mu_0}\{\varsigma-\eta\bar{r}-\varsigma\bar{t}e^{-L_0[G(\mu_0)/\mu_0+k_1+k_2]}\} ;(30a)\\ \tilde{t}(\mu_0)&=t(\mu_0)+e^{-L_0G(\mu_0)/\mu_0}\nonumber\\ &=\frac{\eta-ab\varsigma\bar{r}}{\mu_0}e^{-L_0G(\mu_0)/\mu_0}-\frac{\bar{t}\eta}{\mu_0}+e^{-L_0G(\mu_0)/\mu_0} ; (30b)\end{align} where \(\tilde{t}(\mu_0)=t(\mu_0)+e^-L_0G(\mu_0)/\mu_0\) is the total transmission, including the direct solar beam e^{-L₀G(μ₀)/μ₀}. In previous studies, the direct and diffuse components of canopy radiation have been treated separately (Dickinson, 1983; Sellers, 1985), as have the diffuse and direct albedo calculations. Based on Eq. (30), the relationship between them is revealed.

In the above discussion, the direct and diffuse albedo, transmittance and absorptance are calculated under the boundary condition that the soil is black [both the diffuse reflection \(\bar{r}_\rm s\) and direct reflection r _s(μ₀) of the soil background is 0]. The general solution with a non-zero soil background reflectance can be found from the adding method of (Coakley et al., 1983), which is based on the result of the black soil solution. By assuming the diffuse reflection \(\bar{r}_\rm s\) and direct reflection r _s(μ₀) of the soil background, we obtain the general solution with a non-zero soil background reflectance via the doubling and adding method (Chou, 1992; Tian et al., 2007; Zhang et al., 2013b): \begin{align} r_{\rm co}(\mu_0)&=r(\mu_0)+\bar{t}\frac{r_{\rm s}(\mu_0)e^{-L_0G(\mu_0)/\mu_0}\!+\!\bar{r}_{\rm s}[\tilde{t}(\mu_0)-e^{-L_0G(\mu_0)/\mu_0}]}{1-\bar{r}\bar{r}_{\rm s}} ;(31a)\\ \tilde{t}_{\rm co}(\mu_0)&=\frac{[\tilde{t}(\mu_0)-e^{-L_0G(\mu_0)/\mu_0}] +e^{-L_0G(\mu_0)/\mu_0}r_{\rm s}(\mu_0)\bar{r}}{1-\bar{r}\bar{r}_{\rm s}}\nonumber\\ &\quad +e^{-L_0G(\mu_0)/\mu_0} ;(31b)\\ a_{\rm co}(\mu_0)&=1-r(\mu_0)-e^{-L_0G(\mu_0)/\mu_0}[1-r_{\rm s}(\mu_0)]\nonumber\\ &\quad-[\tilde{t}(\mu_0)-e^{-L_0G(\mu_0)/\mu_0}](1-\bar{r}_{\rm s}) ;(31c)\\ \bar{t}_{\rm co}&=\frac{\bar{t}}{1-\bar{r}\bar{r}_{\rm s}} ;(31d)\\ \bar{r}_{\rm co}&=\bar{r}+\frac{\bar{t}\bar{r}_{\rm s}\bar{t}}{1-\bar{r}\bar{r}_{\rm s}} ;(31e)\\ \bar{a}_{\rm co}&=1-\bar{r}_{\rm co}-\bar{t}(1-\bar{r}_{\rm s}) .(31f) \end{align} In the above equation, r _co(μ₀)/r _co and a _co(μ₀) are the combined albedo and absorption for direct/diffuse radiation; \(\tilde{t}_\rm co(\mu_0)/t_\rm co\) is the combined transmission at the surface for direct/diffuse radiation.

In this section, the accuracy of the Eddington approximation is systematically investigated. Since the hemispheric constant method (Coakley and Chylek, 1975) is widely used in canopy RT, the comparison is also extended to the hemispheric constant method of Dickinson's (1983) model.

The benchmark results are calculated from the successive orders of the scattering approximation (SOSA) with the 128-node Gaussian quadrature (Myneni et al., 1987). We assume that the orientations of leaves are distributed with an equal probability: G(μ)=0.5. The single scattering albedo ω is set as 0.1 and 0.9, which are the typical values for the visible and near-IR spectral regions. For Eddington approximation and the hemispheric constant method, the phase function is P(μ,μ')=1+3gμ'μ. The HG phase function is applied to the calculation of the benchmark results.

In the case of zero surface reflection i.e. \(\bar{r}_\rm s=r_\rm s(\mu_0)=0\), the benchmark results and relative errors of the Eddington approximation and hemispheric constant methods are presented in Fig. 3. When L₀ and μ₀ become larger and smaller, respectively, the reflection and absorption become larger. For ω=0.1, the results for reflection are generally more accurate in the region of small leaf area index for both methods. The reflection of the benchmark results for ω=0.9 is larger than that for ω=0.1 under the same condition. The relative errors are generally less than 2% for the Eddington approximation, which is more accurate than the hemispheric constant method. For transmission, the Eddington approximation is also much more accurate than the hemispheric constant method. For example, when ω=0.1, the error of transmission for the hemispheric constant method is less than 20% for L₀<2.5, but similar error occurs for the Eddington approximation in a much smaller region of L₀ <4.5. For ω=0.9, the relative error of the Eddington approximation is less than 15%, while the error is up to 20% for the hemispheric constant method in the region of L₀ close to 10.

Figure 3.  Benchmark relative errors of the Eddington approximation and hemispheric constant methods in the isotropic case. L₀ is the total leaf area index and μ₀ is the cosine of the solar zenith angle. The soil direct/diffuse reflection r _s(μ₀)/r _s is equal to 0.

Figure 4.  As in Fig. 3 but for r _s(μ₀)=r _s=0.2.

Figure 5.  As in Fig. 3 but for the anisotropic case. The asymmetry factor g is in the range -0.2 to 0; μ₀=0.5; and r _s(μ₀)=r _s=0.

In the case of non-zero soil background reflection i.e. \(\bar{r}_\rm s=r_\rm s(\mu_0)=0.2\), the benchmark results and error comparisons for the two methods are presented in Fig. 4. In this case, the benchmark result for reflection is larger than that for the case of zero surface reflection, since the surface reflection enhances the canopy reflection. For ω=0.1, the region of large relative error (>10%) for reflection is smaller for the Eddington approximation method compared to that of the hemispheric constant method. For transmission, the hemispheric constant method is more accurate at large solar zenith angles (μ₀<0.2), while the Eddington method is more accurate at large leaf area index values (L₀>4). For ω=0.9, the difference in reflection between the two methods becomes very small. However, the hemispheric constant method produces larger errors for transmission compared to the Eddington approximation method.

One of the important new results in this study is the derivation of the asymmetry factor for anisotropic scattering. For the anisotropic case, the HG function is also used to parameterize P(μ,μ'). The asymmetry factor g, as a function of σ and χ_L, is applied to the RT equation directly. In order to simplify the discussion, χ_L is set to 0, then the asymmetry factor, $$ g=\frac{1}{3}(2\sigma-1) . $$ The asymmetry factor g is about -0.2 to 0, while the value range of σ varies from 0.2 to 0.5. A negative value indicates strong backscattering. The single scattering albedo ω is set as 0.1 and 0.9, and both \(\bar{r}_\rm s\) and r _s(μ₀) are set to 0.

In Fig. 5, the accuracy of the two RT methods is investigated under the anisotropic scattering condition. For ω=0.1, the accuracy of the Eddington approximation method is similar to the hemispheric constant method. However, for ω=0.9, the Eddington approximation method becomes more accurate than the hemispheric constant method. In particular, the transmission based on the hemispheric constant method shows very poor results for large cumulative leaf area index, with relative errors up to 20% or higher. This can be greatly improved by using the Eddington approximation method, with the relative error limited to 5% for a large cumulative leaf area index (L₀ up to 10).

The canopy phase function plays a crucial role in the RT through vegetation canopies. In this study, the anisotropic canopy phase function, which is dependent on the leaf normal distribution and leaf reflection/transmission, is derived and a new parameterization for the canopy asymmetry factor is proposed. This new parameterization is much more accurate than the existing scheme for the canopy asymmetry factor. In addition, the analytical solution of the Eddington approximation is obtained. It is interesting to see that there is a relationship between the direct radiation and diffuse radiation. In previous studies, the direct radiation and diffuse radiation have been dealt with separately.

The accuracy in reflection and transmission is examined through comparison with the benchmark result of SOSA. It is shown that the new method, when based on Eddington approximation, can substantially improve the accuracy compared to the previously preferred hemispheric constant method, under both isotropic and anisotropic conditions. Therefore, the canopy albedo can be evaluated more accurately by the analytical solution of non-zero soil background reflection. Given the analytical nature of the solution and its high accuracy, applying this method in land surface radiation models can be carried out very conveniently.