Comparison of δ-2AA, δ-4AA, δ-2DOM, δ-4DOM, δ-2VIM, and δ-4VIM for flux (W m−2) at the top and surface by using the mid-latitude winter atmospheric profile with .
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.
- atmospheric radiative transfer
- solar four-stream discrete ordinates adding method
- solar four-stream harmonic expansion approximation
- infrared variational iteration method
- infrared absorption approximation
- integro-differential equation
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  have been widely used in climate models . 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 , 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  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 . 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  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−2 for cirrus clouds . The 4DOM results tend to be very accurate, but the calculation process is complicated. Zhang et al.  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  for solar RT, AA, and VIM  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  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
where is the cosine of the zenith angle (we use the positive and negative to denote the upward and downward light beams, respectively), is optical depth, is single-scattering albedo, and is the azimuthally averaged scattering phase function, defining the light incidence at , which is scattered in the direction . is the incoming solar flux and is the solar zenith angle. The can be written as
where is the Legendre function. The can be determined from the orthogonal property of Legendre polynomials: . In our notations, and . is the asymmetry factor.
2.1. Four-stream discrete ordinates scheme
where the quadrature point , the weight with , , and .
After a lengthy and laborious derivation, the solution of (3) is given by
where G is determined by the boundary conditions, and values of other parameters are shown in .
The adding method which is used to deal with multilayer RT is shown in .
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
where is determined by the boundary conditions, and values of other parameters can be found at .
The 4SDA can also be applied to solve multilayer solar RT. The detailed process can be found at .
3. Parameterizations for infrared radiative transfer
The infrared RT equation for intensity is
where , , , and are the same as one in Eq. (1). The Planck function is approximated exponentially in optical depth [17, 18] as where . is the optical depth of the considered layer. Planck functions and are evaluated by the temperature of the top () and bottom () of the layer, respectively.
3.1. Absorption approximation
In absorption approximation (AA), the scattering phase function is simplified as a function , and the infrared RT equation becomes
We use as the radiative intensity from the absorption approximation. The solution of Eq. (12) is
where is corresponding to the upward (downward) path and is the initial point at optical depth for the downward direction and at 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
where is the function to be solved and , , and are an inhomogeneous term, a linear operator, and a nonlinear operator, respectively. If is found in the (n)th iteration, the (n + 1)th-order functional solution is
Here, is considered as a restricted variation [i.e., ], where represents the variation. The term 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
According to VIM theory, under the conditions that , and the restricted variation , Eq. (16) yields
For the above correction functional to be stationary [i.e., ], we require the following stationary conditions:
Therefore the Lagrangian multiplier can be identified as
In two-stream VIM (2VIM), the term of in Eq. (22) is simplified as
where . 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 and downward intensity at . Detailed process can be found at . The upward and downward fluxes are written as
We consider a surface boundary condition with a surface emissivity as follows in 2VIM
where , , and are the surface temperature, upward intensity, and downward intensity, respectively, at the surface.
In four-stream VIM (4VIM), the term ofin Eq. (22) is simplified as
The surface boundary condition with an emissivity is given by
To enhance the accuracy of the radiative schemes, a δ-function adjustment is used to adjust the optical parameters following Wiscombe . 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  is used in this study. This model adopts the correlated-k distribution method for gaseous transmission, including five major greenhouse gases, H2O, CO2, O3, NO2, and CH4.
In benchmark calculations, the discrete ordinates model  with a 128-stream scheme (128S) is incorporated with the gaseous transmission scheme of the Fu-Loiu radiation model  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  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) . Two calculations are performed: (1) a low cloud (, ) positioned from 1.0 to 2.0 km and (2) a middle cloud (, ) positioned from 4.0 to 5.0 km. Three solar zenith angles of , , and 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 , the absolute errors of 2DDA are up to about 1.5 and 2.3 K day−1 (relative errors about 6%) for the low and middle clouds. The error of 2SDA is smaller than 2DDA. When , 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.
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  with a 128-stream scheme (δ-128DOM) is used as the benchmark model.
A radiation model  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 H2O, CH4, CO2, NO2, O3, 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−1. 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  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 and liquid water content is located at 1.0–2.0 km, a middle cloud with and is located at 4.0–5.0 km, and a high cloud with a mean effective size and an ice water content is located at 9.0–11.0 km. The surface emissivity 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.
For the low-cloud case, the maximum errors of δ-2AA, δ-2DOM, and δ-2VIM are 20.8, 10.5, and 10.3 K day−1, 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−1, 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−1, respectively, at a height corresponding to the top of the cloud. At the bottom of the cloud, the errors are 10.6 K day−1 for δ-2AA and just 20.05 and 10.05 K day−1 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−1, 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−1 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−1 (Figure 2h). The maximum errors of δ-4AA, δ-4DOM, and δ-4VIM are all about 10.1 K day−1(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 . For the low and middle clouds, the δ-2AA overestimates the upward flux at TOA by 1.8 and 2.1 W m−2, respectively. The errors of the δ-2DOM are up to 21.0 and 21.2 W m−2, while that of the δ-2VIM is limited to 20.7 W m−2. 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−2 and 0.17 W m−2, respectively, while the mean square error of 2AA is limited to 0.08 W m−2.
|Upward flux (TOA)|
|Low, middle, and high clouds||176.3821||179.7768|
|Mean square error||6.0309||0.7939||0.2969||2.7128||0.4456||0.1845|
|Downward flux (surface)|
|Low, middle, and high clouds||302.1410||302.0497|
|Mean square error||0.0806||0.1192||0.0962||0.0794||0.1429||0.1718|
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.
5. Summary and conclusions
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.
The work is supported by the National Key R&D Program of China (2018YFC1507002) and National Natural Science Foundation of China (41675003 and 41675056).