## 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 2^{n} 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)

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

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),

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),

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

where , , , and *e _{j}* (

*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

where , , , , ,and .

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

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

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:

### b. Diffuse radiation

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

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:

Under the type-a boundary condition, we obtain

where . The intensities are obtained by

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

where .

For diffuse intensity, the reflection and transmission are defined as

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

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:

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.

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.

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.

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).

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

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

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

From (20),

where is a 2 × 2 unit matrix, as

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

and

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

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

and

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

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

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

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

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

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

and

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

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

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

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, *g*_{1} = *g*_{2} = 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.

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.

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.

### b. Multilayer atmosphere

For a multilayer atmosphere, the radiation model uses a correlated-*k* distribution scheme (Li and Barker 2005) for gaseous transmission with H_{2}O, O_{3}, CH_{4}, CO_{2}, and N_{2}O 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.

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 *r*_{e} = 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.

A middle cloud with *r*_{e} = 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.

A high cloud with mean effective size *D*_{e} = 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.

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.

From the result of water and ice clouds shown in Figs. 4–8, 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).

## 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

where *b*_{i,j} and *b*_{i} (*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

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

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

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

and

and can be calculated through an upward path as

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

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

and the upward and downward fluxes at TOA are

## REFERENCES

## Footnotes

Earth System Modeling Center Contribution Number 124.