Abstract

The single-layer solutions using a four-stream discrete ordinates method (DOM) in infrared radiative transfer (IRT) have been obtained. Two types of thermal source assumptions—Planck function exponential and linear dependence on optical depth—are considered. To calculate the IRT in multiple layers with a vertically inhomogeneous atmosphere, an analytical adding algorithm has been developed by applying the infrared invariance principle. The derived adding algorithm of the delta-four-stream DOM (δ-4DDA) can be simplified to work for the delta-two-stream DOM (δ-2DDA).

The accuracy for monochromatic emissivity is investigated for both δ-2DDA and δ-4DDA. The relative error for the downward emissivity can be as high as 15% for δ-2DDA, while the error is bounded by 2% for δ-4DDA. By incorporating δ-4DDA into a radiation model with gaseous transmission, δ-4DDA is much more accurate than δ-2DDA. Also, δ-4DDA is much more efficient, since it is an analytical method. The computing time of δ-4DDA is about one-third of the corresponding inverse matrix method.

1. Introduction

The radiative transfer equation (RTE), which describes the radiative transfer process, is an integro-differential equation. Many approximative methods have been developed to solve RTE (e.g., Chandrasekhar 1950; Liou 1974; Coakley and Chýlek 1975; Wiscombe 1977; Meador and Weaver 1980; Liou et al. 1988; Toon et al. 1989; Shibata and Uchiyama 1992; Li and Ramaswamy 1996; Fu et al. 1997; Zhang et al. 2010; van Oss and Spurr. 2002; Lin et al. 2013; Zhang et al. 2013; Zhang and Li 2013), since there are no exact solutions of RTE.

For solar radiation, analytical solutions of various two-stream approximations (Meador and Weaver 1980) have been widely used in climate models (Oreopoulos et al. 2012). More accurate methods of four-stream approximation and six-stream approximation have also been proposed (e.g., Liou et al. 1988; Li and Ramaswamy 1996; Zhang et al. 2010; van Oss and Spurr 2002; Lin et al. 2013; Zhang et al. 2013; Zhang and Li 2013).

Based on the invariance principle (Chandrasekhar 1950), 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., Chou and Suarez 1999; Nakajima et al. 2000; Li et al. 2005).

For solar radiation, 4DDA is much more accurate than the adding method of two-stream DOM (2DDA) in fluxes and heating rate, especially under cloudy-sky conditions, where cloud plays an important role in radiative transfer with different dynamics and microphysics (Lu et al. 2014, 2016). Zhang et al. 2015 used 4DDA in climate models to estimate aerosol direct radiative effect and forcing, showing that 2DDA underestimates the aerosol shortwave direct radiative effect and heating rate significantly.

In Zhang et al. (2015), the four-stream effect in the infrared has not been included. Compared to solar radiation, the scattering effect is much weaker for infrared radiation. Therefore, in most current climate models the method of absorption approximation is widely used in solving infrared radiative transfer (IRT) (Oreopoulos et al. 2012), in which the scattering effect is neglected. However, studies have shown that neglecting scattering can cause an overestimation of outgoing longwave radiation in climate simulations; the difference can be larger than 4 W m−2 (Fu et al. 1997; Li 2000, 2002).

The infrared scattering effect was considered in Wiscombe (1976) by a pure doubling method. First, the solutions are found for each very thin layer by applying the “thin-layer approximation,” then the doubling process is continually carried out to solve the inhomogeneous source and find the solutions for 2n layers (where n is an integer). This method in principle is a numerical approach. It lacks a corresponding adding algorithm for infrared radiation since the invariance principle in Chandrasekhar (1950) has only been applied to solar radiation. In the commonly used absorption approximation method, the downward intensity and upward intensity are decoupled as a result of no scattering effect. This makes the layer-to-layer connection simple and straightforward. When the multiple scattering is included in an infrared radiative transfer process, the downward and upward intensities are coupled. Therefore, the layer-to-layer connection should be dealt with by an adding algorithm, which follows the physics of invariance principle.

Creating an adding algorithm for infrared radiation is the main purpose of this study. First, an analytical solution of four-stream DOM for infrared radiation will be found that is consistent with the four-stream DOM for solar radiation (Liou 2002). There are two types of thermal source assumptions: Planck function exponential or linear dependence on optical depth. Both types will be considered in solutions. Second, the adding method of four-stream DOM for radiation will be built to solve infrared radiative transfer process through a multilayer atmosphere. In addition, the obtained four-stream adding algorithm can be degenerated by using in two-stream DOM.

In section 2, the single-layer solutions by four-stream DOM for infrared radiation are discussed. In section 3, the four-stream adding algorithm, 4DDA, is formulated based on invariance principle. The corresponding result of an adding algorithm to the two-stream DOM (2DDA) is shown in  appendix C. In addition, the accuracy checks of δ-2DDA and δ-4DDA are performed by comparing with δ-64-stream discrete ordinates calculations in section 4. In section 5, a summary is given.

2. Single-layer solution by four-stream discrete ordinates method

a. Infrared radiation

The infrared radiative transfer equation is (e.g., Chandrasekhar 1950)

 
formula

where , τ, and μ are intensity, the optical depth, and the cosine of the zenith angle, respectively. The single scattering albedo ω is defined as the ratio of the scattering cross section to the extinction cross section. The azimuth-independent phase function can be defined as

 
formula

where is equal to , ϕ is equal to the azimuthal angle, and denotes the Legendre function. Based on the orthogonal property of Legendre polynomials, .

The four-stream DOM has been derived in detail by Liou et al. (1988) and Liou (2002). We obtain the infrared radiative transfer equation by using Gaussian quadrature (GQ) and substituting (2) into (1),

 
formula

where the quadrature point , , and . The two-node Gaussian quadrature is used, in which , , and (Sykes 1951). As usual, for upward and downward radiation we have and , respectively.

For a nonisothermal layer, usually the thermal source (Planck function) is assumed to be exponentially or linearly dependent on optical depth. It is shown that the thermal source schemes match with the temperature distribution in the layer if the layer temperature difference is less than 10 K (Li 2002). For the Planck function with exponential dependence on τ (Fu et al. 1997; Li 2002),

 
formula

where and is the total optical depth of the layer. Planck functions and are evaluated by the temperature of the top and the bottom of the layer, respectively.

By substituting (4) into (3), with a lengthy and laborious derivation, the solution of (3) is given by

 
formula

where , , , and ej (j = 1, 2, 3, 4) are shown in  appendix A. In (5) , with the superscript indicating a matrix transpose, can be calculated from the boundary condition of RTE, which is set to be no radiation from the top and the bottom of the layer; that is, and . Therefore, the intensities of RTE are given by

 
formula

where , , , , ,and .

For the Planck function with linear dependence on τ (Toon et al. 1989; Li 2002),

 
formula

where . By substituting (7) into (3), the homogeneous solution of (3) is the same as (5). However, the particular solution is different. The solution of (3) is

 
formula

where the details in solution and the coefficients of and are shown in  appendix B. The coefficient vector can also be calculated from the above boundary condition of RTE:

 
formula

where , , , and and are the same as the above. Equations (6) and (9) are the single-layer solutions for the four-stream approximation. For using in the next section, the intensity matrices of the single-layer can be defined as

 
formula

b. Diffuse radiation

The solution of RTE for diffuse radiation is also used in the adding method. The RTE is

 
formula
 
formula
 
formula

where denotes the incoming diffused intensity; and are used for specifying the incoming intensity, since two-node Gaussian quadrature for the μ interval is adopted. Therefore, there are two forms of boundary conditions corresponding to : type a: [corresponding to downward flux ] and type b: [corresponding to ]. The above boundary conditions are normalized. Similar to (5) and (8), the solution of (11) can also be derived:

 
formula

Under the type-a boundary condition, we obtain

 
formula

where . The intensities are obtained by

 
formula

To distinguish the solution in (6) and (9), the subscript a is used to denote the solution of RTE under the type-a boundary condition.

Similarly, the solution of RTE under the type-b boundary condition is

 
formula

where .

For diffuse intensity, the reflection and transmission are defined as

 
formula
 
formula
 
formula
 
formula

where . In the above equations, the diffuse transmissions include the direct beam transmission . The reflection and transmission matrices for diffuse radiation are defined as

 
formula
 
formula

Matrices in (17) are used in the next section.

3. Adding method of infrared radiation

In solar radiative transfer, Chandrasekhar (1950) has proposed the principles of invariance. The corresponding principles of invariance in infrared radiative transfer have never been discussed before. In the following, four principles of infrared invariance are proposed. To state it clearly, a combination of two homogeneous layers with optical depth and for the first layer and second layer, respectively, is considered. The reflection and transmission functions for the two homogeneous layers are expressed by and for the first layer and by and for the second layer. The superscript asterisk denotes light beams coming from below. In addition, and for a homogeneous atmosphere. Based on the physical meaning of beam radiative transfer, we state the principles as follows:

  1. The upward intensity at the interface of the two layers comes from two parts: the up-outgoing intensity from the second layer and the reflection of the downward intensity through the second layer.

  2. The downward intensity at the interface comes from two parts: the down-outgoing intensity from the first layer and the reflection of the upward intensity through the first layer.

  3. The outgoing intensity at the top of the double layer comes from two parts: the outgoing intensity from the first layer and the total transmission [including the direct beam transmissions of and ] of the upward intensity through the first layer.

  4. The outgoing intensity at the bottom of the double layer comes from two parts: the outgoing intensity from the second layer and the total transmission [ and ] of downward intensity through the second layer.

The schematic diagram of principles 1–4 is given in Fig. 1. These four principles of invariance are different from the four principles of invariance for solar radiation by Chandrasekhar (Zhang et al. 2013). An important characteristic of the infrared is that the thermal source emission acts everywhere. The thermal source emission is contained in the intensities of and as shown in (5) and (8).

Fig. 1.

Schematic diagram of , , , and .

Fig. 1.

Schematic diagram of , , , and .

Based on the above statement, the four principles of invariance in infrared radiative transfer can be written as

 
formula
 
formula
 
formula
 
formula

In the above equations, μ is positive. We define the upward- and downward-intensity matrices in (10). Similarly, the upward and downward internal intensities at the lower boundary of the first layer (interface) are obtained as

 
formula

Substituting (10), (17), and (19) into (18), and using the two-node GQ to decompose the integrations in (18), we have

 
formula
 
formula
 
formula
 
formula

From (20),

 
formula
 
formula
 
formula
 
formula

where is a 2 × 2 unit matrix, as

 
formula

Under the two forms of boundary condition, the upward and downward diffuse intensities at the interface of the two layers are, respectively,

 
formula
 
formula
 
formula

and

 
formula

In the above equations, includes the direct beam transmission . The result can also be written into matrices as

 
formula

The four principles of invariance can also be applied to the process of diffuse radiation. Therefore, we obtain

 
formula
 
formula
 
formula

and

 
formula

In the above equations, μ is positive. As in (18), we use the two-node Gaussian quadrature to decompose the above integrations. Then by using (17) and (23), we obtain

 
formula
 
formula
 
formula
 
formula

From (25), the diffuse beam incident from above can be obtained:

 
formula
 
formula

Similarly, the diffuse beam, which comes from below, can be obtained:

 
formula
 
formula

By applying (21d) and (27a) into multiple layers from layer 1 to layer n in a downward path calculation, we obtain

 
formula
 
formula

By applying (21c) and (26a) into multiple layers extending from the layer N (surface) to layer n + 1 in an upward path calculation, we obtain

 
formula
 
formula

In (28) and (29), and , where is an angular-dependent surface emissivity and is bidirectional reflectance distribution function. For the Lambert surface, and , with the Lambert surface emissivity. Similar to (21a) and (21b), at level n + 1 (bottom boundary of the layer n) are calculated as

 
formula

and

 
formula

Finally, the upward and downward fluxes at level n + 1 are

 
formula
 
formula

and the upward and downward fluxes at the top of atmosphere (TOA) are

 
formula
 
formula

where is a 1 × 2 matrix, as and is a 2 × 2 matrix, as . In the above discussion, the adding method of four-stream DOM is proposed.

For the two-stream DOM of IRT, the adding method can be easily obtained in a similar way (see  appendix C). For 2DDA, Planck function is only approximated exponentially in τ.

4. Comparison results and discussion

The δ-function adjustment is used to deal with the forward peak contribution in order to enhance the accuracy of radiative schemes in the following calculations (Wiscombe 1977; Liou et al. 1988; Iwabuchi 2006; Zhang et al. 2013; and others). Therefore, 2DDA and 4DDA become δ-2DDA and δ-4DDA, respectively. In this section, the accuracy of δ-2DDA and δ-4DDA will be examined by the benchmark result calculated from the discrete ordinates numerical model (Stamnes et al. 1988). A δ-64-stream scheme (δ-64S) is used in the discrete ordinates numerical model.

a. Double layer

First, we investigate the accuracy of δ-2DDA and δ-4DDA in a case of double layers, in which labels “1” and “2” are for the first layer and second layer, respectively. The upward and downward effective emissivities are defined as

 
formula

A case with the same optical depth and the same Plank function for both layers is considered. The single scattering albedo and asymmetry factors are also set the same (ω1 = ω2 = 0.3637, g1 = g2 = 0.8487), these values are very close to the optical properties of water clouds with an effective radius of 6 μm at the wavelength of 11 μm (Fu et al. 1997). The relative errors for the upward and downward emissivities are shown in Fig. 2. For the upward emissivity, the relative error of δ-2DDA generally exceeds 5% for small optical depths with . The relative error for the downward emissivity can be as high as 10% for small values of optical depth. In general, the relative errors become much smaller by using δ-4DDA. The errors are bounded by 2% for both of the upward and downward emissivities.

Fig. 2.

Relative errors of the (left) δ-2DDA and (right) δ-4DDA for (top) the upward and (bottom) the downward emissivity. Τhe benchmark is δ-64S. The optical depth of the two layers ranges from 0.1 to 100. The Henyey–Greenstein (HG) phase function is used and the asymmetry factors are . The surface emissivity ranges from 0 to 1. The single-scattering albedos are ω1 = ω2 = 0.3637.

Fig. 2.

Relative errors of the (left) δ-2DDA and (right) δ-4DDA for (top) the upward and (bottom) the downward emissivity. Τhe benchmark is δ-64S. The optical depth of the two layers ranges from 0.1 to 100. The Henyey–Greenstein (HG) phase function is used and the asymmetry factors are . The surface emissivity ranges from 0 to 1. The single-scattering albedos are ω1 = ω2 = 0.3637.

The values and are used as an example of cirrus cloud (Fu et al. 1997). The error comparisons are shown in Fig. 3. It is found that the stronger scattering causes larger errors for δ-2DDA in both upward and downward emissivities by comparing with the weaker scattering case shown in Fig. 2. In particular, the downward emissivity shows a very poor result for small optical depths, with relative errors up to 15%. The relative errors of upward emissivity are in the range of 2%–5% in most regions. Again, the relative errors are dramatically reduced by using δ-4DDA. For the upward emissivity and downward emissivity, errors are mostly less than 1%, even in thin-optical-depth regions.

Fig. 3.

As in Fig. 2, but for the single-scattering albedos and asymmetry factors .

Fig. 3.

As in Fig. 2, but for the single-scattering albedos and asymmetry factors .

To illustrate that δ-4DDA performs well in IRT, single- and multilayer results are shown in Table 1. Applying δ-4DOM to a medium with and as a single layer, the emissivity is obtained. Then, we also apply δ-4DDA to the same medium, which is divided into 3, 7, and 10 layers. We consider two total different optical depths (τtotal = 1.0 and 5.0). There are slight differences among the values of emissivity in Table 1 because of numerical round-off errors. Therefore, it presents the result of emissivity as the same for the different layers, which illustrates that the multilayer connection is properly solved by δ-4DDA in the four-stream IRT process.

Table 1.

Emissivities for δ-4DDA vs number of sublayers with a single homogeneous layer, in which the temperature is constant; ω = 0.7771 and g = 0.7720 are considered with two optical depths (τtotal = 1.0 and 5.0).

Emissivities for δ-4DDA vs number of sublayers with a single homogeneous layer, in which the temperature is constant; ω = 0.7771 and g = 0.7720 are considered with two optical depths (τtotal = 1.0 and 5.0).
Emissivities for δ-4DDA vs number of sublayers with a single homogeneous layer, in which the temperature is constant; ω = 0.7771 and g = 0.7720 are considered with two optical depths (τtotal = 1.0 and 5.0).

b. Multilayer atmosphere

For a multilayer atmosphere, the radiation model uses a correlated-k distribution scheme (Li and Barker 2005) for gaseous transmission with H2O, O3, CH4, CO2, and N2O in the spectrum range. This model adopts nine infrared bands 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 parameterization of Yang et al. (2015) and Dobbie et al. (1999), respectively.

The U.S. standard atmospheric profile is used (McClatchey et al. 1972) and surface emissivity is set to 0 or 1.

In the inhomogeneous case, we denote δ-4DDA_E for a Planck function with exponential dependence on medium optical depth and δ-4DDA_L for linear dependence, as solutions of (5) and (8). The benchmark results (the results of δ-64S) of clear-sky heating rate, downward flux, and upward flux are shown in Fig. 4 for (first row) and (third row). The absolute errors of δ-2DDA, δ-4DDA_E, and δ-4DDA_L comparing with the results of δ-64S are shown in the second and fourth rows. For heating rate, δ-2DDA produces relatively large errors for , with an error of up to 0.18 K day−1 near the surface; δ-4DDA_E and δ-4DDA_L yield substantially more accurate results in heating rate, with errors less than 0.06 K day−1. Also δ-4DDA_E and δ-4DDA_L are obviously more accurate than δ-2DDA for downward flux. And the errors of δ-2DDA are up to 0.8 W m−2, whereas the errors of δ-4DDA_E and δ-4DDA_L are less than 0.6 W m−2. For upward flux, the errors of δ-4DDA_E and δ-4DDA_L are generally smaller than δ-2DDA in the case of , especially above 4 km.

Fig. 4.

(top) The profiles of (left) heating rate, (center) downward flux, and (right) upward flux for clear sky with (computed from δ-64S). (third row) As in (top), but with . (second row),(bottom) The absolute errors produced by δ-2DDA, δ-4DDA_E, and δ-4DDA_L.

Fig. 4.

(top) The profiles of (left) heating rate, (center) downward flux, and (right) upward flux for clear sky with (computed from δ-64S). (third row) As in (top), but with . (second row),(bottom) The absolute errors produced by δ-2DDA, δ-4DDA_E, and δ-4DDA_L.

It is interesting to find that the four-stream δ-4DDA can be less accurate than the two-stream δ-2DDA, for upward flux in the case of . Under the clear-sky condition, the distribution of radiative intensity in the atmosphere is relatively isotropic. Because of less angular dependence for δ-2DDA, the accuracy of δ-2DDA could be higher than δ-4DDA. Especially, when the upward radiative intensity is dominated by the surface thermal emission. Generally, a single thermal source is in isotropic distribution. When , the upward radiative intensity depends more on local thermal emissions inside the atmosphere. This makes the upward radiative intensity less isotropic because of many thermal sources from different locations. Therefore, δ-4DDA becomes more accurate as a result of its ability to handle the anisotropic intensity. The same is true for downward radiation, which is entirely the result of local emissions. Thus, the error of δ-2DDA is always larger than that of δ-4DDA regardless of whether or . Because of larger errors in downward flux, δ-2DDA is less accurate than δ-4DDA in net flux (downward minus upward), as shown in the results of heating rate.

A low cloud with effective radius of re = 5.89 μm and liquid water content of LWC = 0.22 g m−3 is positioned from 1.0 to 2.0 km. In Fig. 5, the results for low clouds are presented. The maximum errors in heating rates are about 0.18, 0.1, and less than 0.1 K day−1 for δ-2DDA, δ-4DDA_E, and δ-4DDA_L, respectively. The results are not sensitive to the surface emissivity, which is different from the clear-sky case. For downward flux, the errors of δ-2DDA are generally larger in comparison with δ-4DDA_E and δ-4DDA_L. For upward flux, a large error about 1.2 W m−2 is produced by δ-2DDA, whereas errors are bounded by 0.7 W m−2 for δ-4DDA_E and δ-4DDA_L.

Fig. 5.

As in Fig. 4, but for under low-cloud conditions.

Fig. 5.

As in Fig. 4, but for under low-cloud conditions.

A middle cloud with re = 6.2 μm and LWC = 0.28 g m−3 is positioned from 4.0 to 5.0 km. In Fig. 6, the results under middle-cloud conditions are illustrated. Similar to Fig. 5, the accuracy of δ-4DDA is better than that of δ-2DDA. For , the errors in heating rate are about 1.6, 0.1, and 0.1 K day−1 for δ-2DDA, δ-4DDA_E, and δ-4DDA_L, respectively. For , the errors in heating rate are about 0.55, 0.05, and 0.05 K day−1 for δ-2DDA, δ-4DDA_E, and δ-4DDA_L, respectively.

Fig. 6.

As in Fig. 4, but for under middle-cloud conditions.

Fig. 6.

As in Fig. 4, but for under middle-cloud conditions.

A high cloud with mean effective size De = 41.1 μm and ice water content IWC = 0.0048 g m−3 is positioned from 9.0 to 11.0 km. Figure 7 gives the results for an ice-cloud case. Based on Figs. 1 and 2, the errors of δ-2DDA should increase, since the optical depth of the ice cloud is much smaller than that of the water cloud. In addition, the size of ice crystals is generally much larger than that of water droplets; thus, the scattering effect is much enhanced. Compared to the two-stream method, four-stream method can handle scattering more appropriately because more scattering phase function moments are addressed in a four-stream scheme (Yi et al. 2014). It is shown in Fig. 7 that the maximum errors in heating rates for δ-2DDA is close to 0.35 K day−1. The errors in downward flux by δ-2DDA reach 1.4 W m−2, larger than the result for water clouds. The accuracy in heating rate and flux increases a lot by using the two δ-4DDA schemes.

Fig. 7.

As in Fig. 4, but for under high-cloud conditions.

Fig. 7.

As in Fig. 4, but for under high-cloud conditions.

The results of heating rate and downward and upward fluxes for low, middle, and high clouds are shown in Fig. 8. Again, δ-4DDA yields better results than δ-2DDA.

Fig. 8.

As in Fig. 4, but for the sky containing all three cloud types.

Fig. 8.

As in Fig. 4, but for the sky containing all three cloud types.

From the result of water and ice clouds shown in Figs. 48, generally δ-4DDA_L is slightly more accurate than δ-4DDA_E.

Rapid computations are required for weather and climate modeling. Table 2 illustrates the results of efficiency for δ-2DDA and δ-4DDA. The inverse matrix method used in the Fu–Liou model (Liou et al. 1988; Fu 1991; Fu et al. 1997) is also used to solve layer connection in δ-two-stream DOM scheme (denoted as δ-2DOM) and in δ-four-stream DOM scheme (denoted as δ-4DOM). From Table 2, there is no significant difference between δ-2DDA and δ-2DOM. For a pure radiative transfer calculation without considering gaseous transmission, the computational time of δ-4DDA is only about one-third of δ-4DOM. The efficiency of δ-4DDA is caused by avoiding general eigenvalue decomposition and also caused by the unrolling loop in δ-4DDA to avoid general matrix operations. For a radiative transfer calculation with gaseous transmission and cloud absorption, the computational time of δ-4DDA is about 43% of δ-4DOM. Besides, the advantage of δ-4DDA can be particularly efficient in dealing with the partial cloud in climate modeling and cloud retrieval (Nakajima et al. 2000; Gabriel et al. 2006).

Table 2.

Timing of IRT calculations using various approximation methods (normalized to the δ-2DOM method). Gaseous transmission and cloud absorption are included in “radiation model,” but not in “algorithm only.”

Timing of IRT calculations using various approximation methods (normalized to the δ-2DOM method). Gaseous transmission and cloud absorption are included in “radiation model,” but not in “algorithm only.”
Timing of IRT calculations using various approximation methods (normalized to the δ-2DOM method). Gaseous transmission and cloud absorption are included in “radiation model,” but not in “algorithm only.”

5. Summary

The main purpose of this study is to develop an adding algorithm for infrared radiation based on the invariance principle. For the Planck function with exponential and linear terms of the optical depth, the solutions for four-stream radiative transfer and their corresponding adding schemes of δ-4DDA_E and δ-4DDA_L have been created. The four-stream result can be degenerated to a two-stream result as δ-2DDA.

A wide range of accuracy checks δ-2DDA and δ-4DDA are performed in a double-layer case. It is shown that δ-4DDA is much more accurate than δ-2DDA, especially for small optical depths.

Application of δ-2DDA and δ-4DDA into a realistic atmospheric profile, the errors in heating rate by δ-2DDA are as high as 1.0 K day−1, but less than 0.1 K day−1 by δ-4DDA. For downward flux, the errors by δ-2DDA are up to 1.5 W m−2, while those by δ-4DDA are reduced to 0.8 W m−2. A dramatic improvement occurs in the upward flux, where the errors by δ-2DDA can reach 1.6 W m−2 in the middle-cloud case but are reduced to less than 0.7 W m−2 by δ-4DDA.

The computing time involved by δ-2DDA and δ-4DDA is minimal since both of them are analytical methods. For the pure algorithm, the computational time of δ-4DDA is one-third of the corresponding inverse matrix method used in Fu–Liou model.

In summary, in order to consider the cloud infrared scattering effect, the accurate solution of infrared radiative transfer and its related adding process are needed. This has been achieved by δ-2DDA and δ-4DDA. This methodology of infrared δ-4DDA matches with the methodology of solar δ-4DDA (Zhang et al. 2013); the two parts together form a complete four-stream radiative transfer scheme for climate models. The coding of δ-4DDA is available from the authors upon request.

Acknowledgments

The authors thank three anonymous reviewers for their constructive comments and Professor P. Yang for his editorial effort. The work is supported by the State Key Development Program for Basic Research of China (2015CB452805), National Natural Science Foundation of China (41305004, 91537213, and 41675003), Foundation of State Key Laboratory of Severe Weather (2016LASW-B07), and Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).

APPENDIX A

The Coefficients of Eq. (5)

The coefficients , , and of Eq. (5) are , , , , , , and , where ( and ), , , , , , , and with , ( and ), , , , , , , (), and ().

APPENDIX B

The Solution of Four-Stream Infrared Radiative Transfer with Linear Dependence on Optical Depth

By inserting (7) into (3), the differential equation can be written in a matrix form,

 
formula

where bi,j and bi (i, j = −2, −1, 1, 2) are defined in  appendix A.

To simplify the above equation, by defining the sums and the differences of the upward and downward intensities as new quantities, (B1) can be rewritten as

 
formula

where , (i, j = −2, −1, 1, 2), and are defined in  appendix A.

We assume the particular solution of (B2) in the form of

 
formula

substituting (B3) into (B2), we obtain , , , where , , , with () and . And, c, and are defined in  appendix A.

APPENDIX C

The Adding Process for Two-Stream Infrared Radiative Transfer

The single-layer solution for the two-stream approximation based on DOM has been discussed in detail (Toon et al. 1989; Fu et al. 1997), in which the one-node Gaussian quadrature is used. The adding method for two-stream approximation is similar to 4DDA. However, the matrices , , , , , and of 4DDA are degenerated into scalars as , , , , , and .

The two-stream approximation adding method (2DDA) in IRT is given as follows: From layer 1 to layer n, and can be calculated through a downward path as

 
formula

and

 
formula

and can be calculated through an upward path as

 
formula
 
formula

where . Therefore, upward and downward internal intensities ( and , respectively) can be obtained as

 
formula
 
formula

Therefore, the upward and downward fluxes at level n + 1 are

 
formula
 
formula

and the upward and downward fluxes at TOA are

 
formula
 
formula

REFERENCES

REFERENCES
Chandrasekhar
,
S.
,
1950
:
Radiative Transfer.
Oxford University Press, 393 pp.
Chou
,
M.-D.
, and
M. J.
Suarez
,
1999
: A solar radiation parameterization for atmospheric studies. NASA Tech. Memo. NASA/TM-1999-104606, Vol. 15, 40 pp. [Available online at https://gmao.gsfc.nasa.gov/pubs/docs/Chou136.pdf.]
Coakley
,
J. A. J.
, and
P.
Chýlek
,
1975
:
The two-stream approximation in radiative transfer: Including the angle of the incident radiation
.
J. Atmos. Sci.
,
32
,
409
418
, doi:.
Dobbie
,
J. S.
,
J.
Li
, and
P.
Chýlek
,
1999
:
Two- and four-stream optical properties for water clouds and solar wavelengths
.
J. Geophys. Res.
,
104
,
2067
2079
, doi:.
Fu
,
Q.
,
1991
: Parameterization of radiative processes in vertically nonhomogeneous multiple scattering atmospheres. Ph.D. thesis, University of Utah, 259 pp. [Available from University Microfilm, 305 N. Zeeb Rd., Ann Arbor, MI 48106.]
Fu
,
Q.
,
K. N.
Liou
,
M. C.
Cribb
,
T. P.
Charlock
, and
A.
Grossman
,
1997
:
Multiple scattering parameterization in thermal infrared radiative transfer
.
J. Atmos. Sci.
,
54
,
2799
2812
, doi:.
Gabriel
,
P.
,
M. J.
Christi
, and
G. L.
Stephens
,
2006
:
Calculation of Jacobians for inverse radiative transfer: An efficient hybrid method
.
J. Quant. Spectrosc. Radiat. Transfer
,
97
,
209
227
, doi:.
Iwabuchi
,
H.
,
2006
:
Efficient Monte Carlo methods for radiative transfer modeling
.
J. Atmos. Sci.
,
63
,
2324
2339
, doi:.
Li
,
J.
,
2000
:
Gaussian quadrature and its application to infrared radiative transfer
.
J. Atmos. Sci.
,
57
,
753
765
, doi:.
Li
,
J.
,
2002
:
Accounting for unresolved clouds in a 1D infrared radiative transfer model. Part I: Solution for radiative transfer, including cloud scattering and overlap
.
J. Atmos. Sci.
,
59
,
3302
3320
, doi:.
Li
,
J.
, and
V.
Ramaswamy
,
1996
:
Four-stream spherical harmonic expansion approximation for solar radiative transfer
.
J. Atmos. Sci.
,
53
,
1174
1186
, doi:.
Li
,
J.
, and
H. W.
Barker
,
2005
:
A radiation algorithm with correlated-k distribution. Part I: Local thermal equilibrium
.
J. Atmos. Sci.
,
62
,
286
309
, doi:.
Li
,
J.
,
S.
Dobbie
,
P.
Raisanen
, and
Q.
Min
,
2005
:
Accounting for unresolved cloud in solar radiation
.
Quart. J. Roy. Meteor. Soc.
,
131
,
1607
1629
, doi:.
Lin
,
L.
,
Q.
Fu
,
H.
Zhang
,
J.
Su
,
Q.
Yang
, and
Z.
Sun
,
2013
:
Upward mass fluxes in tropical upper troposphere and lower stratosphere derived from radiative transfer calculations
.
J. Quant. Spectrosc. Radiat. Transfer
,
117
,
114
122
, doi:.
Liou
,
K.-N.
,
1974
:
Analytic two-stream and four-stream solutions for radiative transfer
.
J. Atmos. Sci.
,
53
,
1174
1186
, doi:.
Liou
,
K.-N.
,
2002
:
An Introduction to Atmospheric Radiation.
3rd ed. Academic Press, 583 pp.
Liou
,
K.-N.
,
Q.
Fu
, and
T. P.
Ackerman
,
1988
:
A simple formulation of the delta-four-stream approximation for radiative transfer parameterizations
.
J. Atmos. Sci.
,
45
,
1940
1947
, doi:.
Lu
,
C.
,
Y.
Liu
,
S.
Niu
, and
S.
Endo
,
2014
:
Scale dependence of entrainment-mixing mechanisms in cumulus clouds
.
J. Geophys. Res. Atmos.
,
119
,
13 877
13 890
, doi:.
Lu
,
C.
,
Y.
Liu
,
G. J.
Zhang
,
X.
Wu
,
S.
Endo
,
L.
Cao
,
Y.
Li
, and
X.
Gao
,
2016
:
Improving parameterization of entrainment rate for shallow convection with aircraft measurements and large-eddy simulation
.
J. Atmos. Sci.
,
73
,
761
773
, doi:.
McClatchey
,
R. A.
,
R. W.
Fenn
,
J. E. A.
Selby
,
F. E.
Volz
, and
J. S.
Garing
,
1972
:
Optical Properties of the Atmosphere
. Air Force Rep. AFCRL-71-0279, 85 pp.
Meador
,
W. E.
, and
W. R.
Weaver
,
1980
:
Two-stream approximations to radiative transfer in planetary atmospheres: A unified description of existing methods and a new improvement
.
J. Atmos. Sci.
,
37
,
630
643
, doi:.
Nakajima
,
T.
,
M.
Tsukamoto
,
Y.
Tsushima
,
A.
Numaguti
, and
T.
Kimura
,
2000
:
Modeling of the radiative process in an atmospheric general circulation model
.
Appl. Opt.
,
39
,
4869
4878
, doi:.
Oreopoulos
,
L.
, and Coauthors
,
2012
:
The Continual Intercomparison of Radiation Codes: Results from Phase I
.
J. Geophys. Res.
,
117
,
D06118
, doi:.
Shibata
,
K.
, and
A.
Uchiyama
,
1992
:
Accuracy of the delta-four-stream approximation in inhomogeneous scattering atmospheres
.
J. Meteor. Soc. Japan
,
70
,
1097
1110
.
Stamnes
,
K.
,
S. C.
Tsay
,
W. J.
Wiscombe
, and
K.
Jayaweera
,
1988
:
Numerically stable algorithm for discrete ordinate method radiative transfer in multiple scattering and emitting layered media
.
Appl. Opt.
,
27
,
2502
2509
, doi:.
Sykes
,
J. B.
,
1951
:
Approximate integration of the equation of transfer
.
Mon. Not. Roy. Astron. Soc.
,
111
,
377
386
, doi:.
Toon
,
O. B.
,
C. P.
McKay
, and
T. P.
Ackerman
,
1989
:
Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres
.
J. Geophys. Res.
,
94
,
16 287
16 301
, doi:.
van Oss
,
R. F.
, and
R. J.
Spurr
,
2002
:
Fast and accurate 4 and 6 stream linearized discrete ordinate radiative transfer models for ozone profile retrieval
.
J. Quant. Spectrosc. Radiat. Transfer
,
75
,
177
220
, doi:.
Wiscombe
,
W. J.
,
1976
:
Extension of the doubling method to inhomogeneous sources
.
J. Quant. Spectrosc. Radiat. Transfer
,
16
,
477
489
, doi:.
Wiscombe
,
W. J.
,
1977
:
The delta-M method: Rapid yet accurate radiative flux calculations
.
J. Atmos. Sci.
,
34
,
1408
1422
, doi:.
Yang
,
P.
,
K.
Liou
,
L.
Bi
,
C.
Liu
,
B.
Yi
, and
B.
Baum
,
2015
:
On the radiative properties of ice clouds: Light scattering, remote sensing, and radiation parameterization
.
Adv. Atmos. Sci.
,
32
,
32
63
, doi:.
Yi
,
B.
,
X.
Huang
,
P.
Yang
,
B. A.
Baum
, and
G. W.
Kattawar
,
2014
:
Considering polarization in MODIS-based cloud property retrievals by using a vector radiative transfer code
.
J. Quant. Spectrosc. Radiat. Transfer
,
146
,
540
548
, doi:.
Zhang
,
F.
, and
J.
Li
,
2013
:
Doubling-adding method for delta-four-stream spherical harmonic expansion approximation in radiative transfer parameterization
.
J. Atmos. Sci.
,
70
,
3084
3101
, doi:.
Zhang
,
F.
,
Z.
Shen
,
J.
Li
,
X.
Zhou
, and
L.
Ma
,
2013
:
Analytical delta-four-stream doubling-adding method for radiative transfer parameterizations
.
J. Atmos. Sci.
,
70
,
794
808
, doi:.
Zhang
,
H.
,
F.
Zhang
,
Q.
Fu
,
Z.
Shen
, and
P.
Lu
,
2010
:
Two- and four-stream combination approximations for computation of diffuse actinic fluxes
.
J. Atmos. Sci.
,
67
,
3238
3252
, doi:.
Zhang
,
H.
,
Q.
Chen
, and
B.
Xie
,
2015
:
A new parameterization for ice cloud optical properties use in BCC-RAD and its radiative impact
.
J. Quant. Spectrosc. Radiat. Transfer
,
150
,
76
86
, doi:.

Footnotes

Earth System Modeling Center Contribution Number 124.