Atmospheric Radiative Transfer Parameterizations

Abstract

Various radiative transfer (RT) schemes are presented in the chapter including four-stream discrete ordinates adding method (4DDA), four-stream harmonic expansion approximation (4SDA) for the solar spectra and absorption approximation (AA), variational iteration method (VIM) for the infrared spectra. 4DDA uses Gaussian quadrature method to deal with the integration in the RT equation. 4SDA considers four-order spherical harmonic expansion in radiative intensity. VIM allows the zeroth-order solution to be identified as AA, and the scattering effect is included in the first-order iteration. By applying 4DDA/4SDA to a realistic atmospheric profile with gaseous transmission considered, it is found that the accuracy of 4DDA/4SDA is superior to two stream spherical harmonic (Eddington approximation) adding method (2SDA) and two-stream discrete ordinates adding method (2DDA), especially for the cloudy conditions. It is shown that the relative errors of 4DDA/4SDA are generally less than 1% in heating rate, while the relative errors of both 2SDA and 2DDA are over 6%. By applying VIM to a full radiation algorithm a gaseous gaseous transmission, it is found that VIM is generally more accurate than the discrete ordinates method (DOM). Computationally, VIM is slightly faster than DOM in the two-stream case but more than twice as fast in the four-stream case. In view of its high overall accuracy and computational efficiency, 4DDA, 4SDA, as well as VIM are well suited in solving radiative transfer in climate models.

1. Introduction

Solving the radiative transfer (RT) equation is a key issue when dealing with solar/infrared radiative processes related to climate simulations (e.g., radiative flux and heating rate calculation at each level of a climate model). In recent decades, numerous approximation techniques have been proposed to solve the RT equation [1, 2, 3, 4, 5].

For solar radiation, analytical solutions of various two-stream approximations [6] have been widely used in climate models [7]. More accurate methods of four-stream approximation and six-stream approximation have also been proposed (e.g., [8, 9, 10, 11, 12, 13]). Based on the invariance principle [14], an adding method of four-stream discrete ordinates method (DOM) (4DDA) has been proposed. The advantage of the adding method is that it can handle partial cloud in climate models (e.g., [15, 16, 17]). 4DDA is much more accurate than the adding method of two-stream DOM (2DDA) in flux and heating rate calculation, especially under cloudy-sky conditions, where cloud plays an important role in radiative transfer with different dynamics and microphysics [18, 19]. Zhang and Li [20] proposed a new solar RT parameterization (4SDA) which is based on four-stream harmonic expansion to simulate RT in climate modeling.

Since scattering is much weaker for infrared radiation than for solar radiation, an absorption approximation (AA) is used in most current climate models [7]. In AA, the scattering phase function is replaced by a δ-function, meaning that scattering is neglected in all but the forward direction. The advantage of AA is that the upward and downward transfer processes are completely separated. If scattering is considered, the upward and downward intensities are coupled. Under this condition, Fu [21] derived an inverse matrix method to deal with multiple scattering in the atmosphere by using the two-stream DOM (2DOM) and the four-stream DOM (4DOM). It is found that 4DOM is more accurate than 2DOM, especially for thin optical depths. Even if scattering is included, TOA errors for 2DOM can still be >2 W m⁻² for cirrus clouds [22]. The 4DOM results tend to be very accurate, but the calculation process is complicated. Zhang et al. [23] apply the four-stream adding method to resolve multilayer infrared RT.

In the past two decades, the variational iteration method (VIM) has been used to deal with many nonlinear problems [24, 25, 26, 27, 28, 29, 30, 31]. In addition, VIM has been shown to converge faster to the exact solution than other methods do. An accurate solution can be found for most nonlinear problems in only one or two iterations for small solution domains in VIM [27, 31]. Because of its flexibility, convenience and efficiency, VIM has been applied to various nonlinear equations ([27, 28, 29] and many others). In VIM, a nonlinear problem is separated into linear and nonlinear terms, where the latter are usually difficult to deal with and are initially approximated as a first guess. Subsequently, a correction function is constructed by using a general Lagrange multiplier that can be identified optimally via variational theory. Although the infrared RT equation is not a nonlinear equation, it contains an integral term for the scattering, and this is very complicated to deal with. Therefore, VIM is applied to solve the infrared RT equation.

The objective of the chapter is to introduce 4DDA/4SDA scheme [15] for solar RT, AA, and VIM [32] schemes and apply them to resolve solar/infrared radiative transfer in RT models. The simulation results and their comparison to results from the 128-stream discrete ordinates method radiative transfer scheme [33] are shown in Section 4. Finally, a summary is given in Section 5.

2. Parameterizations for solar radiative transfer

The azimuthally averaged solar RT equation is

$$\mu \frac{\mathit{dI}\left(\tau ,\mu \right)}{\mathit{d\tau }}=I\left(\tau ,\mu \right)-\frac{\omega }{2}{\int }_{\mathrm{-1}}^{1}I\left(\tau ,{\mu }^{\text{'}}\right)P\left(\mu ,{\mu }^{\text{'}}\right)d{\mu }^{\text{'}}-\frac{\omega }{4\pi }{F}_0{e}^{-\tau /{\mu }_0}P\left(\mu ,-{\mu }_0\right)$$ E1

where $\mu$ is the cosine of the zenith angle (we use the positive and negative $\mu$ to denote the upward and downward light beams, respectively), $\tau$ is optical depth, $\omega$ is single-scattering albedo, and $P\left(\mu ,{\mu }^{\text{'}}\right)$ is the azimuthally averaged scattering phase function, defining the light incidence at ${\mu }^{\text{'}}$ , which is scattered in the direction $\mu$ . $F_0$ is the incoming solar flux and ${\mu }_0$ is the solar zenith angle. The $P\left(\mu ,{\mu }^{\text{'}}\right)$ can be written as

$$P\left(\mu ,{\mu }^{\text{'}}\right)=\sum _{l=0}^N{\omega }_l{P}_l\left(\mu \right)P_l\left({\mu }^{\text{'}}\right)$$ E2

where $P_l$ is the Legendre function. The ${\omega }_l$ can be determined from the orthogonal property of Legendre polynomials: ${\omega }_l=\left[\left(2l+1\right)/2\right]\int P\left(\text{cos}\mathrm{\Theta }\right)P_l\left(\text{cos}\mathrm{\Theta }\right)d\text{cos}\mathrm{\Theta }$ . In our notations, ${\omega }_0=1$ and ${\omega }_1=3\mathrm{g}$ . $g$ is the asymmetry factor.

2.1. Four-stream discrete ordinates scheme

The four-stream discrete ordinates scheme (4DDA) use the Gaussian quadrature method to deal with the integration in Eq. (1). Eq. (1) yields

$${\mu }_i\frac{\mathit{dI}\left(\tau ,{\mu }_i\right)}{\mathit{d\tau }}=I\left(\tau ,{\mu }_i\right)-\frac{\omega }{2}\sum _{l=0}^3{\omega }_l{P}_l\left({\mu }_i\right)\sum _{j=-2}^{2}I\left(\tau ,{\mu }_j\right)P_l\left({\mu }_j\right)a_j-\frac{\omega }{4\pi }{F}_0{e}^{-\tau /{\mu }_0}\sum _{l=0}^3{\omega }_l{P}_l\left({\mu }_i\right)P_l\left(-{\mu }_0\right)$$ E3

where the quadrature point ${\mu }_{-i}=-{\mu }_i$ , the weight $a_{-i}=a_i$ with $i=\pm 1,\pm 2$ , ${\mu }_1=0.2113248$ , ${\mu }_2=0.7886752$ and $a_1=a_2=0.5$ .

After a lengthy and laborious derivation, the solution of (3) is given by

$$\left[\begin{array}{l}I\left(\tau ,{\mu }_2\right)\\ I\left(\tau ,{\mu }_1\right)\\ I\left(\tau ,{\mu }_{-1}\right)\\ I\left(\tau ,{\mu }_{-2}\right)\end{array}\right]=\left[\begin{array}{llll}{\phi }_2^{+}{e}_2& {\phi }_1^{+}{e}_1& {\phi }_1^{-}{e}_3& {\phi }_2^{-}{e}_4\\ {\phi }_2^{+}{e}_2& {\phi }_1^{+}{e}_1& {\phi }_1^{-}{e}_3& {\phi }_2^{-}{e}_4\\ {\phi }_2^{-}{e}_2& {\phi }_1^{-}{e}_1& {\phi }_1^{+}{e}_3& {\phi }_2^{+}{e}_4\\ {\phi }_2^{-}{e}_2& {\phi }_1^{-}{e}_1& {\phi }_1^{+}{e}_3& {\phi }_2^{+}{e}_4\end{array}\right]G+\left[\begin{array}{l}{Z}_2^{+}\\ Z_1^{+}\\ Z_1^{-}\\ Z_2^{-}\end{array}\right]\text{exp}\left(-\tau /{\mu }_0\right)$$ E4

where G is determined by the boundary conditions, and values of other parameters are shown in [13].

The adding method which is used to deal with multilayer RT is shown in [13].

2.2. Four-stream spherical harmonic expansion scheme

The purpose of the four-stream spherical harmonic expansion (4SDA) is to separate out the angle-dependent intensity by assuming

$$I\left(\tau ,\mu \right)=\sum _{l=0}^3{I}_l\left(\tau \right)P_l\left(\mu \right)$$ E5

By substituting Eqs. (2) and (5) into Eq. (1) and using the orthogonality relation of the Legendre function in $-1\le \mu \le 1$ , we obtain

$$\frac{{\mathit{dI}}_1}{\mathit{d\tau }}=a_0{I}_0-b_0{e}^{-\tau /{\mu }_0}$$ E6

$$2\frac{{\mathit{dI}}_2}{\mathit{d\tau }}+\frac{{\mathit{dI}}_0}{\mathit{d\tau }}=a_1{I}_1-b_1{e}^{-\tau /{\mu }_0}$$ E7

$$3\frac{{\mathit{dI}}_3}{\mathit{d\tau }}+2\frac{{\mathit{dI}}_1}{\mathit{d\tau }}=a_2{I}_2-b_2{e}^{-\tau /{\mu }_0}$$ E8

$$3\frac{{\mathit{dI}}_2}{\mathit{d\tau }}=a_3{I}_3-b_3{e}^{-\tau /{\mu }_0}$$ E9

The solution of Eqs. (6)–(9) is

$$\left[\begin{array}{l}{I}_0\\ I_1\\ I_2\\ I_3\end{array}\right]=\left[\begin{array}{llll}{e}_1& e_3& e_2& e_4\\ R_1{e}_1& -R_1{e}_3& R_2{e}_2& -R_2{e}_4\\ {\widehat{P}}_1{e}_1& {\widehat{P}}_1{e}_3& {\widehat{P}}_2{e}_2& {\widehat{P}}_2{e}_4\\ Q_1{e}_1& -Q_1{e}_3& Q_2{e}_2& -Q_2{e}_4\end{array}\right]G+\left[\begin{array}{l}{\eta }_0\\ {\eta }_1\\ {\eta }_2\\ {\eta }_3\end{array}\right]\text{exp}\left(-\tau /{\mu }_0\right)$$ E10

where $G$ is determined by the boundary conditions, and values of other parameters can be found at [20].

The 4SDA can also be applied to solve multilayer solar RT. The detailed process can be found at [20].

3. Parameterizations for infrared radiative transfer

The infrared RT equation for intensity $I\left(\tau ,\mu \right)$ is

$$\mu \frac{\mathit{dI}\left(\tau ,\mu \right)}{\mathit{d\tau }}=I\left(\tau ,\mu \right)-\left(1-\omega \right)B\left(\tau \right)-\frac{\omega }{2}{\int }_{\mathrm{-1}}^{1}I\left(\tau ,{\mu }^{\text{'}}\right)P\left(\mu ,{\mu }^{\text{'}}\right)d{\mu }^{\text{'}}$$ E11

where $\mu$ , $\tau$ , $\omega$ , and $P$ are the same as one in Eq. (1). The Planck function is approximated exponentially in optical depth [17, 18] as $B\left(\tau \right)=B_0{e}^{-\mathit{\beta \tau }}$ where $\beta =-\left[\text{ln}\left(B_1/B_0\right)\right]/{\tau }_1$ . ${\tau }_1$ is the optical depth of the considered layer. Planck functions $B_0$ and $B_1$ are evaluated by the temperature of the top ( $\tau =0$ ) and bottom ( $\tau ={\tau }_1$ ) of the layer, respectively.

3.1. Absorption approximation

In absorption approximation (AA), the scattering phase function $P\left(\mu ,{\mu }^{\text{'}}\right)$ is simplified as a $\delta \left(\mu ,{\mu }^{\text{'}}\right)$ function [22], and the infrared RT equation becomes

$$\mu \frac{{\mathit{dI}}^{\mathit{0}}\left(\tau ,\mu \right)}{\mathit{d\tau }}=\left(1-\omega \right)I^0\left(\tau ,\mu \right)-\left(1-\omega \right)B_0{e}^{-\mathit{\beta \tau }}$$ E12

We use $I^0$ as the radiative intensity from the absorption approximation. The solution of Eq. (12) is

$$I^0\left(\tau ,\pm \mu \right)=\left[I^0\left({\tau }_s,\pm \mu \right)\mp \frac{1-\omega }{\mathit{\mu \beta }\pm \left(1-\omega \right)}{B}_0{e}^{-{\mathit{\beta \tau }}_s}\right]e^{\pm \left(1-\omega \right)\left(\tau -{\tau }_s\right)/\mu }\pm \frac{1-\omega }{\mathit{\mu \beta }\pm \left(1-\omega \right)}{B}_0{e}^{-\mathit{\beta \tau }}$$ E13

where $\pm \mu$ is corresponding to the upward (downward) path and ${\tau }_s$ is the initial point at optical depth ${\tau }_s=0$ for the downward direction and at ${\tau }_s={\tau }_1$ for the upward direction. Generally, AA is applied to the two-stream case (δ-2AA) and the four-stream case (δ-4AA). In the following subsection, the AA solution is used as the initial approximate solution of the VIM.

3.2. Variational iteration method

A general nonlinear system is used to illustrate the basic idea of variational iteration method (VIM) [24, 25, 26, 27, 28, 29, 30, 31]:

$$\mathit{L\upsilon }\left(x\right)+\mathit{N\upsilon }\left(x\right)=h\left(x\right)$$ E14

where $\upsilon \left(x\right)$ is the function to be solved and $h\left(x\right)$ , $L$ , and $N$ are an inhomogeneous term, a linear operator, and a nonlinear operator, respectively. If $\upsilon \left(x\right)$ is found in the (n)th iteration, the (n + 1)th-order functional solution is

$${\upsilon }^{n+1}\left(x\right)={\upsilon }^n\left(x\right)+{\int }_0^x\lambda \left(\zeta \right)\left[{\mathit{L\upsilon }}^n\left(\zeta \right)+N{\tilde{\upsilon }}^n\left(\zeta \right)-h\left(\zeta \right)\right]\mathit{d\zeta }$$ E15

Here, ${\tilde{\upsilon }}^n\left(\zeta \right)$ is considered as a restricted variation [i.e., ${\widehat{\delta }\tilde{\upsilon }}^n\left(\zeta \right)=0$ ], where $\widehat{\delta }$ represents the variation. The term $\lambda \left(\zeta \right)$ is a general Lagrange multiplier and should be identified optimally by using variational theory.

Based on VIM, the functional reiteration of the infrared RT equation can be deduced as

$$\begin{array}{l}{I}^{n+1}\left(\tau ,\pm \mu \right)=I^n\left(\tau ,\pm \mu \right)+{\int }_{{\tau }_s}^{\tau }{\lambda }^{\pm }\left(s\right)[\frac{{\mathit{dI}}^n\left(s,\pm \mu \right)}{\mathit{ds}}\mp \frac{I^n\left(s,\pm \mu \right)}{\mu }\pm \frac{\omega }{2\mu }{\int }_{\mathrm{-1}}^1{\tilde{I}}^n\left(s,{\mu }^{\text{'}}\right)P\left(\pm \mu ,{\mu }^{\text{'}}\right)d{\mu }^{\text{'}}\\ \pm \left(1-\omega \right)B_0{e}^{-\mathit{\beta s}/\mu }]\mathit{ds}\end{array}$$ E16

According to VIM theory, under the conditions that ${\left.\widehat{\delta }{I}^n\left(s,\pm \mu \right)\right|}_{s={\tau }_s}=0$ , $\widehat{\delta }{\int }_{{\tau }_s}^{\tau }\frac{\left(1-\omega \right)B_0{e}^{-\mathit{\beta s}}}{\mu }\mathit{ds}=0$ and the restricted variation ${\widehat{\delta }\tilde{I}}^n\left(s,{\mu }^{\text{'}}\right)=0$ , Eq. (16) yields

$${\left.\widehat{\delta }{I}^{n+1}\left(\tau ,\pm \mu \right)=\left[1+{\lambda }^{\pm }\left(s\right)\right]\widehat{\delta }{I}^n(s,\pm \mu )\right|}_{s=\tau }+\widehat{\delta }{\int }_{{\tau }_s}^{\tau }\left[-\frac{{\mathit{d\lambda }}^{\pm }\left(s\right)}{\mathit{ds}}\mp \frac{{\lambda }^{\pm }\left(s\right)}{\mu }\right]I^n\left(s,\pm \mu \right)\mathit{ds}$$ E17

For the above correction functional to be stationary [i.e., $\widehat{\delta }{I}^{n+1}\left(\tau ,\pm \mu \right)=0$ ], we require the following stationary conditions:

$${\left.-\frac{{\mathit{d\lambda }}^{\pm }\left(s\right)}{\mathit{ds}}\mp \frac{{\lambda }^{\pm }\left(s\right)}{\mu }=0\text{and}{\lambda }^{\pm }\left(s\right)\right|}_{s=\tau }=-1$$ E18

Therefore the Lagrangian multiplier ${\lambda }^{\pm }\left(s\right)$ can be identified as

$${\lambda }^{\pm }\left(s\right)=-e^{\pm \left(\tau -s\right)/\mu }$$ E19

By substituting Eq. (19) into Eq. (16), we obtain

$$\begin{array}{l}{I}^{n+1}\left(\tau ,\pm \mu \right)=I^n\left(\tau ,\pm \mu \right)-{\int }_{{\tau }_s}^{\tau }{e}^{\pm \left(\tau -s\right)/\mu }[\frac{{\mathit{dI}}^n\left(s,\pm \mu \right)}{\mathit{ds}}\mp \frac{I^n\left(s,\pm \mu \right)}{\mu }\pm \frac{\omega }{2\mu }{\int }_{\mathrm{-1}}^1{I}^n\left(s,{\mu }^{\text{'}}\right)P\left(\pm \mu ,{\mu }^{\text{'}}\right){\mathit{d\mu }}^{\text{'}}\\ \pm \left(1-\omega \right)B_0{e}^{-\mathit{\beta s}/\mu }]\mathit{ds}\end{array}$$E20

In two-stream VIM (2VIM), the term of ${\int }_{\mathrm{-1}}^1{I}^n\left(s,{\mu }^{\text{'}}\right)P\left(\pm \mu ,{\mu }^{\text{'}}\right){\mathit{d\mu }}^{\text{'}}$ in Eq. (22) is simplified as

$${\int }_{\mathrm{-1}}^1{I}^n\left(s,{\mu }^{\text{'}}\right)P\left(\pm \mu ,{\mu }^{\text{'}}\right){\mathit{d\mu }}^{\text{'}}=\left(1+3{\mathit{g\mu \mu }}_{\mathit{1}}\right)I^n\left(s,{\mu }_1\right)+\left(1+3{\mathit{g\mu \mu }}_{-\mathit{1}}\right)I^n\left(s,{\mu }_{-1}\right)$$ E21

where ${\mu }_1=-{\mu }_{-1}=1/1.66$ . The AA expression shown in Eq. (15) is used as the initial zeroth-order solution. By substituting Eq. (15) into Eq. (22), we obtain the first-order solution of upward intensity at $\tau =0$ and downward intensity at $\tau ={\tau }_1$ . Detailed process can be found at [30]. The upward and downward fluxes are written as

$$F^{+}\left(0\right)={\mathit{\pi I}}^1\left(0,{\mu }_1\right)\text{and}{F}^{-}\left({\tau }_1\right)={\mathit{\pi I}}^1\left({\tau }_1,{\mu }_1\right)$$ E22

We consider a surface boundary condition with a surface emissivity ${\epsilon }_s$ as follows in 2VIM

$$I_s\left({\mu }_1\right)=\left(1-{\epsilon }_s\right)I_s\left({\mu }_{-1}\right)+{\epsilon }_{s}B\left(T_s\right)$$ E23

where $T_s$ , $I_s\left({\mu }_1\right)$ , and $I_s\left({\mu }_1\right)$ are the surface temperature, upward intensity, and downward intensity, respectively, at the surface.

In four-stream VIM (4VIM), the term of ${\int }_{\mathrm{-1}}^1{I}^n\left(s,{\mu }^{\text{'}}\right)P\left(\pm \mu ,{\mu }^{\text{'}}\right){\mathit{d\mu }}^{\text{'}}$ in Eq. (22) is simplified as

$${\int }_{\mathrm{-1}}^1{I}^n\left(s,{\mu }^{\text{'}}\right)P\left({\mu }_i,{\mu }^{\text{'}}\right){\mathit{d\mu }}^{\text{'}}=\frac{1}{2}\sum _{j=-2}^2\left(1+3{\mathit{g\mu }}_i{\mu }_j\right)I^n\left(s,{\mu }_j\right)$$ E24

where ${\mu }_1=-{\mu }_{-1}=0.2113248$ and ${\mu }_2=-{\mu }_{-2}=0.7886752$ . The AA expression shown in Eq. (15) is used as the initial zeroth-order solution. Detailed process can be found at [30]. The upward and downward fluxes are written as

$$F^{+}\left(0\right)=2\pi \left[a_1{\mu }_1{I}^1\left(0,{\mu }_1\right)+a_2{\mu }_2{I}^1\left(0,{\mu }_2\right)\right]$$ E25

$$F^{-}\left({\tau }_1\right)=2\pi \left[a_1{\mu }_1{I}^1\left({\tau }_1,{\mu }_{-1}\right)+a_2{\mu }_2{I}^1\left({\tau }_1,{\mu }_{-2}\right)\right]$$ E26

The surface boundary condition with an emissivity ${\epsilon }_s$ is given by

$$I_s\left({\mu }_i\right)=2\left(1-{\epsilon }_s\right)\sum _{j=1}^2{a}_j{\mu }_j{I}_s\left(-{\mu }_j\right)+{\epsilon }_{s}B\left(T_s\right)$$ E27

To enhance the accuracy of the radiative schemes, a δ-function adjustment is used to adjust the optical parameters following Wiscombe [34]. We refer to the VIM solutions with the δ-function applied as δ-2VIM and δ-4VIM.

4. Results and discussion

4.1. Solar spectra

RT in the atmosphere is a complicated process. It depends not only on the single-layer direct reflection and transmission but also on the diffuse results and the gaseous transmission, cloud-aerosol scattering and absorption, etc. It is important to evaluate errors in radiation under a variety of atmospheric conditions by using a radiation algorithm. The Fu-Liou radiation model [35] is used in this study. This model adopts the correlated-k distribution method for gaseous transmission, including five major greenhouse gases, H₂O, CO₂, O₃, NO₂, and CH₄.

In benchmark calculations, the discrete ordinates model [33] with a 128-stream scheme (128S) is incorporated with the gaseous transmission scheme of the Fu-Loiu radiation model [35] by replacing the existing radiative transfer algorithm in the model. Also the two-stream discrete ordinates adding method (2DDA), Eddington adding method (2SDA), four-stream discrete ordinates adding method (4DDA), and four stream spherical harmonic adding method (4SDA) schemes are incorporated with the same gaseous transmission scheme. The atmosphere was vertically divided into 280 layers, each of thickness 0.25 km. The mid-latitude winter of atmospheric profile [36] is used with a surface albedo of 0.2 for each band. For a water cloud the optical properties are parameterized in terms of liquid water content (LWC) and effective radius (re) [36]. Two calculations are performed: (1) a low cloud ( $\mathrm{LWC}=0.22{\mathrm{gm}}^{‐3}$ , $r_e=5.89\mathrm{\mu m}$ ) positioned from 1.0 to 2.0 km and (2) a middle cloud ( $\mathrm{LWC}=0.28{\mathrm{gm}}^{‐3}$ , $r_e=6.19\mathrm{\mu m}$ ) positioned from 4.0 to 5.0 km. Three solar zenith angles of ${\mu }_0=1$ , ${\mu }_0=0.5$ , and ${\mu }_0=0.25$ are generally considered.

In Figure 1, the benchmark results of heating rate are shown for three solar zenith angles under low/middle cloud condition (Figure 1a–c for low cloud condition, Figure 1g–i for middle cloud condition) as well as the absolute errors of 2DDA, 2SDA, 4DDA, and 4SDA against the benchmark results (Figure 1d–f for cloud condition, Figure 1j–l for middle cloud condition). When ${\mu }_0=1$ , the absolute errors of 2DDA are up to about 1.5 and 2.3 K day⁻¹ (relative errors about 6%) for the low and middle clouds. The error of 2SDA is smaller than 2DDA. When ${\mu }_0=0.25$ , the result becomes opposite as the error of 2SDA becomes larger than 2DDA. By using 4DDA and 4SDA, relative errors are much suppressed. In general, the relative error becomes less than 1% for the three solar zenith angles. This indicates both 4DDA and 4SDA are accurate enough to obtain cloud-top solar heating.

Figure 1.

Heating rate profiles computed from 128S and the error profiles from 2DDA, 2SDA, 4DDA, and 4SDA. Three solar zenith angles ${\mu }_0=1$ , ${\mu }_0=0.5$ , and ${\mu }_0=0.25$ are considered.

4.2. Infrared spectra

For the infrared spectra, the accuracy and efficiency of δ-2AA, δ-4AA, δ-2VIM, δ-4VIM, δ-2DOM, and δ-4DOM will be systematically compared. In addition, the discrete ordinates model [28] with a 128-stream scheme (δ-128DOM) is used as the benchmark model.

A radiation model [17] is used to study the accuracy of the VIM scheme for multiple layers within a model atmosphere. A correlation-k distribution scheme is used to simulate the gaseous transmission with profiles for H₂O, CH₄, CO₂, NO₂, O₃, and CFCs. This model is reasonably efficient because it neglects to scatter for certain intervals with very large absorption coefficients and water vapor continuum at high altitudes. Nine infrared bands are adopted in this model in wavenumber ranges 0–340, 340–540, 540–800, 800–980, 980–1100, 1100–1400, 1400–1900, 1900–2200, and 2200–2500 cm⁻¹. The optical properties of ice and water clouds are calculated based on the radiative property parameterization of [37, 38, 39], respectively. The mid-latitude winter atmospheric profiles [36] are used. The atmospheric profile is divided into homogeneous layers with a geometrical thickness of 0.25 km.

In this model, a low cloud with an effective radius $r_e=5.89\mathrm{\mu m}$ and liquid water content $\mathrm{LWC}=0.22{\mathrm{gm}}^{‐3}$ is located at 1.0–2.0 km, a middle cloud with $r_e=5.89\mathrm{\mu m}$ and $\mathrm{LWC}=0.22{\mathrm{gm}}^{-3}$ is located at 4.0–5.0 km, and a high cloud with a mean effective size $D_e=41.1\mathrm{\mu m}$ and an ice water content $\mathrm{IWC}=0.0048{\mathrm{gm}}^{-3}$ is located at 9.0–11.0 km. The surface emissivity ${\epsilon }_s$ is set to 1.

In the left column of Figure 2, the benchmark heating rates calculated by δ-128DOM are given under conditions of low clouds (Figure 2a), middle clouds (Figure 2d), high clouds (Figure 2g) and the all-sky case containing a combination of low, middle, and high clouds (Figure 2j). The middle column of Figure 2 shows the errors in the calculated heating rates of δ-2AA, δ-2DOM, and δ-2VIM compared to those of δ-128DOM. Furthermore, the right column of Figure 2 shows the errors in the calculated heating rates of δ-4AA, δ-4DOM, and δ-4VIM compared to those of δ-128DOM.

Figure 2.

Heating rate profiles calculated by δ-128DOM and the error profiles produced by δ-2AA, δ-4AA, δ-2DOM, δ-4DOM, δ-2VIM, and δ-4VIM in the mid-latitude winter atmospheric profile with surface emissivity ${\epsilon }_s=1$ for (a)–(c) low cloud, (d)–(f) middle cloud, (g)–(i) high cloud, and (j)–(l) all three cloud types.

For the low-cloud case, the maximum errors of δ-2AA, δ-2DOM, and δ-2VIM are 20.8, 10.5, and 10.3 K day⁻¹, respectively, at a height corresponding to the top of the cloud (Figure 2b). The maximum errors of δ-4AA, δ-4DOM, and δ-4VIM are 20.9, 10.1, and 10.03 K day⁻¹, respectively, at the same height (Figure 2c). This shows that, for low-level water clouds, δ-2VIM (δ-4VIM) is more accurate than δ-2AA (δ-4AA) and δ-2DOM (δ-4DOM).

For the middle cloud, the maximum errors of δ-2AA, δ-2DOM, and δ-2VIM are approximately 20.9, 11.1, and 10.8 K day⁻¹, respectively, at a height corresponding to the top of the cloud. At the bottom of the cloud, the errors are 10.6 K day⁻¹ for δ-2AA and just 20.05 and 10.05 K day⁻¹ for δ-2DOM and δ-2VIM, respectively (Figure 2e). At the top of the cloud, δ-4AA, δ-4DOM, and δ-4VIM produce biases of approximately 11.1, 10.7, and 10.5 K day⁻¹, respectively (Figure 2f).

For the high-cloud case, the optical depth of the ice cloud is much smaller than that of the water cloud. Based on the results shown in Figure 2, the accuracies of δ-2DOM and δ-2VIM are similar and are better than that of δ-2AA. Furthermore, δ-4AA, δ-4DOM, and δ-4VIM have very similar accuracies. As seen in the third row in Figure 2, δ-2AA produces an error of approximately 10.6 K day⁻¹ at a height corresponding to the top of a high cloud. The accuracies of δ-2DOM and δ-2VIM are similar; the maximum errors of both are approximately 10.5 K day⁻¹ (Figure 2h). The maximum errors of δ-4AA, δ-4DOM, and δ-4VIM are all about 10.1 K day⁻¹(Figure 2i). However, the difference is that the maximum errors occur at top of the high cloud for δ-4AA but on the bottom for δ-4DOM and δ-4VIM.

In the case of all three clouds together, δ-2VIM is more accurate than δ-2AA and δ-2DOM for the low and middle clouds (Figure 2k). In general, δ-4DOM and δ-4VIM are comparable in accuracy for cloud heating rate (Figure 2l).

The results of upward (downward) flux at the TOA (surface) for the six schemes are presented in Table 1 for mid-latitude winter profiles with the surface emissivity ${\epsilon }_s=1$ . For the low and middle clouds, the δ-2AA overestimates the upward flux at TOA by 1.8 and 2.1 W m⁻², respectively. The errors of the δ-2DOM are up to 21.0 and 21.2 W m⁻², while that of the δ-2VIM is limited to 20.7 W m⁻². Also, the δ-4VIM is generally more accurate than the δ-4DOM. However, the results for the high clouds are reversed, with the δ-2VIM errors becoming larger than that of δ-2DOM. For the downward flux at the surface, the AA schemes are generally more accurate than the corresponding DOM and VIM schemes except in the middle-cloud case. We also calculate the mean square error of each scheme. For upward flux at TOA, the mean square error of VIM is generally smaller than one of DOM and AA. However, the mean square error of 2VIM and 4VIM is 0.10 W m⁻² and 0.17 W m⁻², respectively, while the mean square error of 2AA is limited to 0.08 W m⁻².

Atmospheric condition	128DOM	2AA	2DOM	2VIM	4AA	4DOM	4VIM
Upward flux (TOA)
Low clouds	209.5326	211.2735; (1.7409)	208.5206; (−1.012)	209.0425; (−0.4901)	210.9585; (1.4259)	209.1645; (−0.3681)	209.5458; (0.0132)
Middle clouds	183.5646	185.7200; (2.1554)	182.4027; (−1.1619)	182.9060; (−0.6586)	182.4002; (−1.1644)	183.1096; (−0.4550)	183.5357; (−0.0289)
High clouds	198.1357	200.3545; (2.2188)	198.6394; (0.5037)	198.8373; (0.7016)	199.2441; (1.1084)	197.2351; (−0.9006)	197.4354; (−0.7003)
Low, middle, and high clouds	176.3821	179.7768; (3.3947)	175.6421; (−0.7400)	176.2356; (−0.1465)	178.8788; (2.4967)	175.5891; (−0.7930)	175.8856; (−0.4965)
Mean square error		6.0309	0.7939	0.2969	2.7128	0.4456	0.1845
Downward flux (surface)
Low clouds	302.7236	302.7267; (0.0031)	302.9533; (0.2297)	302.9384; (0.2148)	302.7429; (0.0193)	302.8662; (0.1426)	302.8449; (0.1213)
Middle clouds	287.4458	286.9155; (−0.5303)	287.7658; (0.3200)	287.6141; (0.1683)	287.7658; (0.3200)	287.5678; (0.1220)	287.5418; (0.096)
High clouds	247.8497	248.0313; (0.1816)	248.3309; (0.4812)	248.3303; (0.4806)	248.2968; (0.4471)	248.5608; (0.7111)	248.6452; (0.7955)
Low, middle, and high clouds	302.1410	302.0497; (−0.0913)	302.4413; (0.3003)	302.4227; (0.2817)	302.0189; (−0.1221)	302.3161; (0.1751)	302.3151; (0.1741)
Mean square error		0.0806	0.1192	0.0962	0.0794	0.1429	0.1718

Table 1.

Comparison of δ-2AA, δ-4AA, δ-2DOM, δ-4DOM, δ-2VIM, and δ-4VIM for flux (W m⁻²) at the top and surface by using the mid-latitude winter atmospheric profile with ${\epsilon }_s=1$ .

The flux differences between the six approximate schemes and δ-128DOM are listed in parentheses. All δ symbols are neglected in the table.

For climate modeling, the efficiency of radiative transfer parameterization is also very important. Table 2 lists the computing times required for δ-2AA, δ-4AA, δ-2DOM, δ-4DOM, δ-2VIM, and δ-4VIM, which were computed by HP EliteDesk 880 G1 TWR with 8 Intel(R) Core(TM) i7–4790 CPUs, 32-bit operating system, and 8GB memory. The results are normalized to that of δ-2AA. The computational efficiency of δ-2VIM is slightly better than that of δ-2DOM, which took more than double the time of δ-2AA. However, δ-4VIM is more than twice as fast as δ-4DOM for the radiation algorithm alone and the radiation model.

	2AA	4AA	2DOM	4DOM	2VIM	4VIM
Algorithm only	1.0	1.4412	2.2223	14.8195	2.0365	5.3615
Radiation model	1.0	1.1095	1.4114	5.8899	1.3759	2.2739

Table 2.

Computing times of δ-2AA, δ-4AA, δ-2DOM, δ-4DOM, δ-2DDA, δ-4DDA, δ-2VIM, and δ-4VIM (normalized by the computing time of δ-2AA).

All δ symbols are neglected in the table.

5. Summary and conclusions

The objective of the paper is to introduce 4DDA/4SDA [13, 20] for solar RT, AA, and VIM [32] for infrared RT and applied them to radiative transfer models.

4DDA use Gaussian quadrature method to deal with the integration in the RT equation. 4SDA is based on four-stream harmonic expansion in radiative intensity. By applying 4DDA/4SDA to a realistic atmospheric profile with gaseous transmission considered, it is found that the accuracy of 4DDA/4SDA is superior to Eddington adding method (2SDA) and two-stream discrete ordinates adding method (2DDA), especially for the cloudy conditions. It is shown that the relative errors of 4SDA are generally less than 1% in heating rate, while the relative errors of both 2SDA and 2DDA are over 6%.

VIM differs from other analytical methods for solving nonlinear differential equations in that it requires neither linearization nor small perturbations. The optimal result is constructed through variation by a Lagrange multiplier. It was found that the scattering term in the infrared RT equation could be dealt with as a nonlinear operator in VIM. By taking the AA solution as the zeroth-order solution, the scattering effect was properly included in the first-order iterative solution.

The six schemes of δ-2AA, δ-4AA, δ-2DOM, δ-4DOM, δ-2VIM, and δ-4VIM were compared systematically against the benchmark results provided by δ-128DOM for a multilayer case. For a multilayer atmosphere, VIM gave more accurate results than those of DOM and AA for the low and middle clouds in both the two- and four-stream cases. However, the errors from VIM for high clouds were similar to those from AA.

Computationally, δ-2VIM and δ-2DDA were slightly faster than δ-2DOM, which took more than double time of δ-2AA. However, δ-4VIM/δ-4DDA was more than twice as fast as δ-4DOM for both the pure radiation algorithm and the radiation model in a layered, cloudy atmosphere. In general, the main benefit of δ-2VIM was an improved accuracy with a computational time similar to that of δ-2DOM. The main benefit of δ-4VIM/δ-4DDA was improved computational efficiency with accuracy similar to that of δ-4DOM.

In view of their overall high accuracy and computational efficiency, the conclusion is that 4DDA, 4SDA, δ-2VIM, and δ-4VIM are well suited for parameterizing the solar/infrared RT in climate models.

Acknowledgments

The work is supported by the National Key R&D Program of China (2018YFC1507002) and National Natural Science Foundation of China (41675003 and 41675056).