## Abstract

Various aspects of infrared radiative transfer through clouds are investigated. First, three solutions to the IR radiative transfer equation are presented and assessed, each corresponding to a different approximation for the Planck function. It is shown that the differences in results between solutions with linear and exponential dependence of the Planck source function are small for typical vertical resolutions in climate models. Second, a new perturbation-based approach to solving the IR radiative transfer equation with the inclusion of cloud scattering is presented. This scheme follows the standard perturbation method, and allows one to identify the zeroth-order equation with the absorption approximation and the first-order equation as including IR scattering effects. This enables one solution to accurately treat cloudy layers in which cloud scattering is included, and allows for an improved and consistent treatment of absorbing aerosol layers in the absence of cloud by using the zeroth-order equation. This new scheme is more simple and efficient compared to previous perturbation method work for treating infrared absorption and scattering. Last, a general method is devised for calculating the random, maximum, and slantwise overlap of cloud layers, which conveniently integrates into the two-stream radiative transfer solution in this work. For several random and maximum (or slantwise) overlap cloud cases with a wide variation of cloud fractions, the error in the cooling rate is generally less than 1 K day^{−1} and the error in the radiative flux is generally less than 3 W m^{−2}.

## 1. Introduction

In the last few decades, the treatment of atmospheric radiation in climate models has achieved remarkable progress: the application of the line-by-line technology and the correlated *k*-distribution method for use in transmission calculations (Clough et al. 1989; Lacis and Oinas 1991; Ramaswamy and Freidenreich 1991; Fu and Liou 1992; Mlawer et al. 1997; Zhong and Haigh 2001, Chou et al. 2001 and others), the application of the Eddington approximation and other newly developed methods to address the radiative transfer process (Shettle and Weinman 1970; Coakley and Chýlek 1975; Joseph et al. 1976; Liou et al. 1988; Stamnes et al. 1988; Li and Ramaswamy 1996, and others), and the application of Mie and other more computationally intensive scattering methods in solving for nonspherical cloud and aerosol particle optical properties and their corresponding parameterizations (e.g., Slingo 1989; Fu et al. 1998; Liou et al. 1999; Dobbie et al. 1999; Lindner and Li 2000). However, treatment of the interaction of radiation with clouds in climate models [mostly the general circulation models (GCMs)] has received comparably less attention and is still far from satisfactory.

Although it is realized that clouds play an important role in the radiation process, there have been few studies that systematically discuss the cloud treatment in general circulation models. The last review of the treatment of radiation in GCMs was presented by Stephens (1984), which is almost two decades old now. This work focuses on addressing some current issues in radiative transfer through clouds for a one-dimensional column model with direct applicability to GCMs. This paper seeks to provide a more recent update on the treatment of radiation in GCMs.

There are two parts to this paper. Part I discusses various aspects of the interaction of IR radiation with clouds. New solutions for the treatment of radiative transfer in GCMs is presented, including such areas as cloud scattering and cloud overlap. In Part II (Li and Barker 2002), the compounding problem of cloud horizontal inhomogeneity in GCMs is addressed.

This work begins with an analysis of the relevant physics pertaining to the absorption approximation in solving the infrared radiative transfer equation. The solutions corresponding to different assumptions of thermal source function (Planck function) are then given. There is focus on three kinds of source functions, one with linear dependence on the optical depth, another with exponential dependence on optical depth, and another that is an isothermal constant source function. The accuracy for each solution is shown by comparing with the results from very high vertical resolution calculations.

In this work, a new perturbation method for cloud infrared scattering is proposed. This method is simpler than the method proposed in Li and Fu (2000), and can handle the radiative transfer through aerosols outside cloud layers more correctly with the zeroth-order equation implicit in absorption approximation method. The full solution is invoked in regions that contain cloud so as to correctly account for the absorption and scattering processes. The method allows for a consistent treatment of cloud and clear layers and is computationally efficient. Moreover, the perturbation method is easy to match with the cloud overlap scheme based on a two-stream (or higher stream) approximation for radiative transfer.

Since individual clouds are often subgrid in horizontal extent in GCMs, it is vital to properly account for cloud overlap. For the radiative treatment of cloud overlap in the infrared, traditionally the cloud matrix method is widely used. This method, however, is computationally expensive both in time and in the use of computer memory. Furthermore, the matrix method cannot be integrated into the two-stream radiative transfer method used in this paper. In view of these limitations, a conceptually new scheme to handle the cloud overlap problem is presented here. This scheme is tested using one-dimensional Monte Carlo calculations.

## 2. Solution to the absorption approximation

The plane-parallel, homogeneous radiative transfer equation for an azimuthally averaged diffuse infrared radiance *I*(*τ, **μ*) is

where *μ* = cos*θ, **θ* is the zenith angle, *τ* is the optical depth, *ω̃* is the single-scattering albedo, *B*[*T*(*τ*)] is the Planck function for a temperature *T, **T* is assumed to satisfy a certain distribution in *τ* (see appendix A), and *P*(*μ, **μ*′) is the azimuthal-independent scattering phase function.

The phase function with the *δ*-Eddington approximation (Joseph et al. 1976) is

where *f* is the fractional scattering into the forward peak, *g* is the asymmetry factor. Since in the infrared, absorption is strong and scattering is weak, the scattering process can be simplified by setting *f* = 1. Thus only the forward-scattering peak (scattering angle is 0°) is kept (and enhanced), and scattering in all other directions is completely neglected.

Under this assumption, the multiple-scattering source term in (1) is rewritten as

Equation (4) is referred to as the absorption approximation (AA). Note that the scattering has not been simply neglected. If it was, then *P*(*μ, **μ*′) = 0 [see (18)] and the corresponding results for cloudy atmospheres would be inaccurate.

As mentioned in Li (2000b), there are two types of solutions for (4). The first one, which is widely used at present in GCMs (e.g., Fels and Schwarzkopf 1981; Morcrette 1990), is based on the direct numerical interpretation of the solution of (4). The solutions of (4) are

where *F*^{+}(*F*^{−}) is the upwelling (downwelling) flux at pressure *p, **p*_{s} is the surface pressure, *T*_{s} is the surface temperature, *T*′ is the temperature at pressure *p*′, and Ψ[*κ*(*p, **p*′)] is the flux transmittance defined in terms of the absorptance depth *κ*(*p, **p*′) for a slab of atmosphere between *p* and *p*′,

with

where *τ*(*p, **p*′) is the optical depth between *p* and *p*′. For clear sky conditions, *ω*(*p*) = 0 and *κ*(*p, **p*′) = *τ*(*p, **p*′).

In (5), the upwelling (downwelling) flux at a level is directly determined by thermal emission and absorption from all layers below (above) that level. Figure 1a shows the transfer process for (5). Since the contributions exchanged between all layers is considered explicitly, the computation time is quadratically proportional to the number of model layers. In the second approach to solving (4), all the exchange contributions from outside of a layer are represented by the incoming flux at the boundary of the layer. The transfer process is shown in Fig. 1b. The calculation of flux is done individually layer by layer. The computation time is therefore linearly proportional to the number of model layers. Li (2000b) shows that these two types of solution are mathematically equivalent. Because of its computational advantage, the second type of solution is used hereafter. Also, it is important to point out that it is much easier to handle the cloud scattering and cloud overlap by using the second type of solution to (4).

Assuming that level 1 is the top of the atmosphere (TOA), level *N* is the surface, and layer *i* is between levels *i*th and (*i* + 1)th there are two types of source functions, isothermal and nonisothermal. For a nonisothermal layer, an optical depth dependence on the Planck function has to be assumed. The dependence of the Planck function on the optical depth for the nonisothermal case is discussed in appendix A.

A commonly used nonisothermal source is the Planck function with a linear optical depth dependence (the L source scheme), which was first proposed by Schuster (Wiscombe 1976). It is assumed that the Planck function within the *i*th layer is

where *α*_{i} = *B*_{i+1} − *B*_{i}, with *B*_{i} and *B*_{i+1} being the Planck functions for the temperature at levels *i* and *i* + 1, respectively, and *τ*_{i} is the optical depth of the layer *i.* From the general solution (B5) we obtain the corresponding solutions of (4) for the upwelling and downwelling radiance at levels *i* and *i* + 1,

where *I*_{i}(±*μ*) and *I*_{i+1}(±*μ*) are the upwelling (downwelling) radiance at level *i* and level *i* + 1, respectively, *κ*_{i} = *ε*_{i}*τ*_{i} is the absorptance optical depth for layer *i,* and *ε*_{i} = 1 − *ω̃*_{i} is the single-scattering coalbedo. The corresponding upwelling and downwelling fluxes, *F*^{+}_{i} and *F*^{−}_{i+1}, can be obtained through angular integration via Gaussian quadrature (Li 2000b). Considering only the two-stream case (one-node Gaussian quadrature), we have

where *ζ*_{i} = *πα*_{i}*μ*_{1}/*κ*_{i}, *B̃*_{i} = *πB*_{i}, and 1/*μ*_{1} = *e*^{1/2} = 1.648 721 3 is the diffusivity factor. This diffusivity factor is slightly different from the widely used value of 1.66. Our value, however, has a rigorous mathematical base, as derived from Gaussian quadrature theory (Li 2000b). The general case for higher-order streams is also shown in Li (2000b).

It is unnecessary to list the downwelling flux, since the physical process is the same for the upwelling and downwelling flux. The expression of *F*^{−}_{i+1} should be the same as that of *F*^{+}_{i} , except that the physical quantities related to “up” and “down” should be exchanged. The general transformations for the “up” and “down” are

The physical quantities unrelated to the transfer direction like *ω*_{i}, *g*_{i}, etc., are unchanged. A good example is that the upwelling and downwelling radiance in (8) follow the up and down exchange rule shown above exactly.

Though there is no singularity in the solution of (9), calculations show that the results are generally unstable in the stratosphere and above since the value of *κ*_{i} could be very small, especially for a correlated *k*-distribution scheme of gaseous transmission. We therefore have to put the following limit on (9)

Another commonly used source function is the Planck function with exponential optical depth dependence (the E source scheme). For the layer *i,* the Planck function is assumed as (e.g., Fu et al. 1998; Li and Fu 2000)

where *β*_{i} = ln(*B*_{i+1}/*B*_{i}).

where *η*_{i}(*μ*) = 1/(1 + *β*_{i}*μ*/*ε*_{i}).

The upwelling flux is

In the solution of (13), there could exist singularity for *η*_{i}(−*μ*_{1}) → ∞, as *κ*_{i}/*μ*_{1} → *β*_{i}. However, through the algebraic derivation we obtain

Therefore as *κ*_{i}/*μ*_{1} → *β*_{i}, accurate results can be achieved by keeping the first few terms in the expansion series, thus the singularity is removed. Similarly, for the downwelling flux the singularity could occur as *η*_{i}(*μ*_{1}) → ∞, at *κ*_{i}/*μ*_{1} → −*β*_{i}. Since *κ*_{i} is always positive, *κ*_{i}/*μ*_{1} → −*β*_{i} occurs in the region of temperature inversions (i.e., temperature increase with height), which could occur within the atmospheric boundary. So it is necessary to eliminate the singularity for the downwelling flux in the same way as for (14).

In the above we have discussed the AA solution for nonisothermal sources. For an isothermal source, that is, the constant Planck function (the C source scheme), the corresponding solution can be obtained directly from (9) or (13) by taking *B*_{i} = *B*_{i+1} = *B*_{i+1/2} and *α*_{i} = 0 (or *β*_{i} = 0), where *B*_{i+1/2} is the Planck function determined by *T* at the center of layer *i.* We have

The differences in cooling rate and flux among these three kinds of solution should be small for the high vertical resolution, since the temperature difference inside each layer is small. The differences between these solutions is supposed to increase as the model vertical resolution decreases. Therefore, the vertical resolution dependence for different solutions for an atmospheric profile has to be shown.

We use the radiation model used in CCC GCM written by Li, which is an algorithm that employs the correlated *k*-distribution for gaseous transmission by resolving the whole spectrum from 0 to 2500 cm^{−1} into a total of 46 intervals in the cumulative probability function space. The cloud and aerosol optical properties are resolved into 9 bands for the same spectrum range of 0–2500 cm^{−1}, and each band contains 2 to 9 intervals in the cumulative probability function space. The continuum scheme is based on numerical calculations from the algorithm of LBLRTM (Clough et al. 1989; Mlawer et al. 1997).

First for the clear sky case, the comparisons for the upwelling flux at TOA and the downwelling flux at the surface are shown in Fig. 2. The middle latitude summer atmospheric profile with 60-km height (McClatchey et al. 1972) is used. The model is set with equal layer thickness. Benchmark results are obtained by using the model with extremely high resolution (layer thickness of 10 m). In Fig. 2, as expected, the errors for both the upwelling flux at the TOA and the downwelling flux at the surface increase as the layer thickness increases. For the L and E source schemes, the differences in both the upwelling flux at TOA and the downwelling flux at the surface are small; while the C source scheme performs completely different from the other two schemes. For the upwelling flux at the TOA the error of the C source scheme is smaller than that of the L source scheme or that of the E source scheme and with opposite sign. But for the downwelling flux at the surface, the error for C source scheme could be as large as 25 W m^{−2} for a model layer thickness close to 2 km. Note, although the Planck source function assumptions for the L and E source schemes are very different, the differences in the upwelling flux at the TOA and the downwelling flux at the surface are less than 0.05 W m^{−2} for layer thickness less than 1 km. This is because in both schemes the Planck functions at the upper and lower boundaries for a layer are the same. Because of this, the difference in emission from a layer inside seems not to be crucial to the modeling situation provided that the layer thickness is not too large. Also, even for layer thicknesses of 1 km, the temperature difference across these layers is usually less than 10 K, since the typical measured atmospheric lapse rate is about 6–7 K km^{−1} in the troposphere. As shown in appendix A, for temperature differences less than 10 K, both the L and E source schemes produce very accurate results compared to the results obtained from the linear interpretation of temperature inside the layer.

For the L and E source schemes, if the model vertical solution is less than 0.25 km the results would be very accurate with error in flux at TOA and surface less than 0.1 W m^{−2}. However this resolution is not usually present in current GCM simulations above the boundary layer. In Fig. 3 the vertical profiles of error in fluxes is presented for the case with layer thickness of 2 km. We can see that the errors in fluxes are mostly produced in the lower atmosphere region, below 12 km. Equation (9) [or (15)] shows that the upwelling flux at level *i* is composed of two parts. The first part is the term *F*^{+}_{i+1 }*e*^{−κi/μ1}, which represents the flux transfer from the lower level; the second part is the remained of (9) [or (15)] and represents the emission from layer *i.* In the upper atmosphere, the first term dominates. This can be seen clearly in the limit result in (10). As the optical depth approaches zero, the contribution of the second part is negligible. The first term is not sensitive to the layer thickness at all. However, in the lower atmosphere, the layer optical depth becomes larger and the emission contribution becomes more important. The emission part depends on how the Planck function is represented in the layer and is sensitive to the model layer thickness.

The above analysis shows that vertical resolution should be high for the lower atmosphere, but can be relatively low for the higher atmosphere. This is true for most of the atmospheric profiles used in climate models. In Fig. 4 an atmospheric profile with two different layer thicknesses is shown. Layer thicknesses of 2 km above 12 km and 0.25 km below 12 km. In comparison with Fig. 3, the errors in flux dramatically decreases for all the three source schemes. This verifies the point that a relatively coarse resolution can be used for atmospheric profiles in high-altitude regions.

For the cloudy case, we use the same middle latitude summer atmospheric profile with vertical layers of 0.25 km. We consider three cloud cases of low cloud, middle cloud, and high cloud. A low cloud is positioned from 1.0 to 2.0 km with effective radius (*r*_{e}) of 5.98 *μ*m. A middle cloud is positioned from 4.0 to 5.0 km with effective radius (*r*_{e}) of 6.2*μ*m. And a high cloud is positioned from 10 to 12 km with general mean effective size (*D*_{ge}) of 30 *μ*m. For the low and middle clouds, the liquid water content (LWC) varies from 0.001 to 1 g m^{−3}; and the parameterization for the water cloud optical properties given by Lindner and Li (2000) is used. For the high cloud, the ice water content (IWC) varies from 0.00001 to 0.01 g m^{−3}, and the parameterization for ice cloud optical properties of Fu et al. (1998) is used.

The upper panels in Fig. 5 show the results of the difference in the upwelling flux at TOA for three different schemes against the high-resolution calculations. The middle panels are the corresponding results for the downwelling fluxes at the surface. It is shown that for the C source scheme, the differences in the upwelling flux at TOA and the downwelling flux at the surface increase as LWC increases. The error for the downwelling flux at the surface is about 2 W m^{−2} for very small LWC, which is consistent with the clear sky results (see Fig. 2). For the high cloud, the change in flux difference is less than 0.1 W m^{−2} with increases in IWC because the optical thickness of the high cloud is small and its influence on flux is limited. For the L and E source schemes, the differences in the upwelling flux at TOA and the downwelling flux at the surface show no obvious increase with increases of LWC or IWC. Figure 5 indicates that the L and E source schemes are reliable for the cloudy sky case.

It is difficult to compare the cooling rate profiles between the high-resolution calculations and the lower-resolution calculations since the cloud cooling profile strongly relies on the model vertical resolution. Therefore, instead of comparing the cloud-cooling rates, a comparison is made of the net energy deposited in a cloud, that is, the net flux into the cloud.

We use the same cloud cases as in the upper and middle panels of Fig. 5. The corresponding errors for net flux into the cloud layers are shown in the lower panels of Fig. 5. Also the benchmark results are from the high vertical resolution calculations with a layer thickness of 10 m. It is shown that for L and E source schemes the errors in deposited energy are very small with a maximum of about 0.1 W m^{−2} for the considered range of LWC (IWC). The same as flux, the C source scheme produces much larger errors in comparison with L and E source schemes.

One advantage for the solution shown in this paper is that the limit case of cooling rate as the optical depth approaches zero can be expressed explicitly. This is difficult to be achieved in the method of (5). In the limit case of an *i*th layer optical depth *τ*_{i} → 0, the cooling rate for layer *i* by the L source scheme becomes

where *c*_{p} is the specific heat at constant pressure, *ρ* is the ambient air density, and *k*_{i} is the cloud extinction coefficient for the layer *i.* The physics of (16a) is very clear in that the cooling rate is determined by the net incoming flux, *F*^{−}_{i} − *F*^{+}_{i+1} , minus the net outgoing (emitted) flux, *B̃*_{i+1} + *B̃*_{i}. For the E source scheme, under the same limitation,

Both (16a) and (16b) reduce to the result of the C source scheme as *α*_{i} = 0 (*β*_{i} → 0) and *B*_{i} = *B*_{i+1} = *B*_{i+1/2}.

Equation (16) could be useful to evaluate the cooling rates for very high-altitude (mesosphere or higher) regions of the atmosphere. Since the flux difference for a layer becomes very small in the high-altitude region, instability in the numerical calculation could occur by using the flux difference method to evaluate the cooling rate. Equation (16) could provide a more stable solution to this problem.

## 3. Scattering

In (9) or (15), the upwelling and downwelling fluxes are decoupled (i.e., independent of each other). This makes flux calculations extremely simple. However, the simple treatment of scattering in (3) is usually insufficient and can lead to cooling rate errors of several K day^{−1} (Toon et al. 1989; Fu et al. 1997). To counter this but at the same time maintain numerical efficiency, Li and Fu (2000) proposed a simple perturbation method to deal with the infrared scattering by cloud droplets. In this method, the zero order approximations of upwelling radiance and downwelling radiance at each model level are computed without scattering effects. Then, the zero-order radiance is used to recalculate the first-order correction in radiance with the proper scattering effects included as a perturbation term. It is shown that this approach could produce accurate results and is computationally efficient.

Based on the method of Li and Fu (2000) the upwelling flux for the two-stream case including scattering effect for the E source scheme is

with

where *ψ*(*μ, **μ*_{1}) = 1 + 3*g*_{i}*μμ*_{1}, *g*_{i} is the asymmetry factor for layer *i, **F*^{(0)±}_{i} is the upwelling (downwelling) flux obtained from the zero order equation.

The advantage for the perturbation method is that the scattering effect calculation of *S*^{+}_{i} only has to be done within the layers that have cloud. This is why the perturbation method is very efficient, usually about 1 to 2 orders faster than the matrix inverse method [see the comparison in Li and Fu (2000)].

Based on (17), four calculation paths are required to obtain *F*^{(0)−}_{i}, *F*^{(0)+}_{i}, *F*^{−}_{i}, and *F*^{+}_{i} (see Fig. 6a). This process can be simplified. First, AA is used to calculate the downwelling flux of *F*^{(0)−}_{i} (*i* = 1, 2, … , *N*). Second, the upwelling flux *F*^{+}_{i} (*i* = 1, 2, … , *N*) is calculated with the proper scattering effect considered by perturbation method shown as (17), but we use *F*^{+}_{i+1} instead of *F*^{(0)+}_{i+1} in calculating *S*^{+}_{i}, that is,

Finally, we calculate the upwelling flux again with the scattering effect included, in this process we use *F*^{+}_{i} instead of *F*^{(0)+}_{i} in calculating *S*^{+}_{i}. Thus, one calculation path for *F*^{(0)+}_{i} (*i* = 1, 2, … , *N*) is saved. The calculation processes for the three calculation path scheme is sketched in Fig. 6b. Since the scattering effect is small in the infrared, the differences in the cooling rates and the downwelling fluxes are very small by using the *F*^{(0)+}_{i+1} or *F*^{+}_{i+1} in calculating the scattering effect. We found that the difference in the cooling rate by using this approximation is not noticeable. The differences in the upwelling fluxes at TOA and in the downwelling fluxes at the surface are less than 0.1 W m^{−2}.

As shown in Li and Fu (2000), the difference in downwelling flux between AA and the perturbation method is always much smaller than that for the upwelling flux. As we mentioned, in AA the forward-scattering peak is enhanced and backward-scattering effect is neglected. For the upwelling flux at the level *i,* the contribution from the forward scattering is proportional to *F*^{(0)+}_{i+1} and the contribution from the backward scattering is proportional to *F*^{(0)−}_{i}. Since, generally, *F*^{(0)+}_{i+1} is larger than *F*^{(0)−}_{i}, the enhancement of the forward-scattering peak and the omission of the backward scattering would produce a larger upwelling flux. For the downwelling flux at level *i* + 1, the situation is the opposite. The forward scattering is proportional to *F*^{(0)−}_{i} and the backward scattering is proportional to *F*^{(0)+}_{i+1}. The reduction of flux from omission of the backward scattering could largely compensate the enhancement of the forward scattering. Thus the error in the downwelling flux for AA is always much smaller compared to that for upwelling flux. The above analysis helps us to understand the error caused by AA.

However, in Li and Fu (2000) the perturbation term is

so the zero order equation is

which is different from the AA in (4), but we use the AA solution as the zero order equation for the perturbation method. This does not follow the standard perturbation method. Also this inconsistency leads to a relatively complicated solution of (17). In (17) there are two types of transmission function *e*^{−τi/μ1} and *e*^{−κi/μ1}. In addition, as mentioned before, the zero-order equation of (18) is not suitable for handling absorbing particles outside cloud layers in comparison with AA (also see discussion in appendix B).

We therefore proposed an improved perturbation method in appendix B, in which the zero-order equation is AA. The solution for the E source scheme for the two-stream case is

with

This solution is simpler than that of (17) and contains only one type of transmission function, *e*^{−κi/μ1}. The physical interpretation for the three terms in *S*^{+}_{i} is straightforward. They represent the effects of the forward scattering, backward scattering, and the internal anisotropic scattering, respectively. The solution of (19) can be extended to any higher order stream following Li and Fu (2000).

Similarly, we can obtain the perturbation solution for the L source scheme case. The upwelling flux is

with

The results for the C source scheme can be obtained by letting *α*_{i} = 0 in (19) [or *β*_{i} = 0 in (20)] and *B*_{i} = *B*_{i+1} = *B*_{i+1/2}.

In Fig. 7, the comparisons in the cloud-cooling rate among (19) and (20) are presented, The results are all presented versus the results of (17), since the accuracy of (17) has been systematically investigated in Li and Fu (2000). A low cloud is positioned from 1.0 to 2.0 km with *r*_{e} = 5.98 *μ*m and LWC = 0.22 g m^{−3}. A middle cloud is positioned from 4.0 to 5.0 km with *r*_{e} = 6.2 *μ*m and LWC = 0.28 g m^{−3}. And a high cloud is positioned from 10 to 12 km with *D*_{ge} = 30 *μ*m and IWC = 0.048 g m^{−3}. It is shown that (19) generates the slightly less cooling rate for the low cloud and the middle cloud in comparison with (17). The difference between results from (17) and (19) are small (relative error <1%), smaller than the difference between (17) and D2S [matrix inverse method, see Li and Fu (2000)].

Table 1 shows the comparisons of the upwelling flux at TOA and the downwelling flux at the surface. It is noted that the differences among the three methods of (17), (19), and (20) are very small.

Li and Fu (2000) shows that the cloud infrared scattering effect has been well addressed by the perturbation method of (17). For example, the error in cooling rate is up to about 4 K day^{−1} for the middle cloud case in Fig. 7 by AA method; while the error is reduced to 0.35 K day^{−1} by using (17). The new perturbation approach of (19) and (20) can produce similar accuracy results, but with more rigorous mathematical and physical basis, as discussed above.

## 4. Cloud overlap

In most of the previous works on infrared radiative transfer, the upwelling (downwelling) flux at one level is determined by thermal emission and absorption from all layers below (above) that level. In other words, the exchange of photons between any two layers is directly considered in the radiative transfer process. Corresponding to the radiative transfer method, a cloud matrix method was proposed to deal with the random and maximum cloud overlap (Manabe and Strickler 1964; Washington and Williamson 1979; Harshvardhan et al. 1987; Liang and Wang 1997; Raisanen 1998; Morcrette and Jakob 2000; Chou et al. 2001). In the cloud matrix method, the matrix elements are used to describe the mutual cloud relation between any two levels. For example, if two cloud layers with cloud fractions *c*_{i} and *c*_{i+n} are located in the layers *i* and *i* + *n,* the two nonadjacent clouds are assumed as random and total cloud coverage between level *i* and level *i* + *n* + 1 is described by the matrix elements

This cloud matrix method has been summarized by Li (2000a). In particular, Li (2000a) found a simple and correct method to calculate the effective cloud emissivity, which has not been correctly handled before. However, based on the two-stream (or higher) radiative transfer method discussed in the previous section, all exchange contributions from the outside of one layer are represented by the boundary conditions of this layer. Therefore, the cloud matrix method for the cloud overlap is no longer appropriate for the two-stream (or higher) infrared radiative transfer scheme. Also the cloud matrix method is expensive in its use of computer memory (Li 2000a). A new scheme has to be developed to handle the cloud overlap in consistency with the two-stream (or higher) radiative transfer method discussed above.

For a single-layer cloud located in the layer *i* with the cloud fraction *c*_{i} (see Fig. 8a), the downwelling flux at the level *i, **F*^{−}_{i}, is obtained by (9) or (13) for cloud-free layers. At the level *i* + 1 there are two downwelling flux paths, one from the clear sky portion and another from the cloudy sky portion. The flux for the clear sky portion is

where *τ*_{i} is the *i*th layer optical depth for the clear sky portion, and the flux for the cloudy sky portion is

where *τ*_{i} is the *i*th layer optical depth for the cloudy sky portion. Here we only use the E source scheme to illustrate the radiative transfer process. The same is for the L source scheme.

Therefore, the downwelling flux at level *i* + 1 is the cloud fractional weighted mean,

Thus, below cloud *c*_{i} the fluxes from the clear and cloudy portions are well mixed.

Considering a random overlap cloud system, we assume there is another cloud located at layer *c*_{i+k}, (*k* = 2, 3, …) (see Fig. 8a). Usually cloud *c*_{i+k} is assumed to be in random overlap to cloud *c*_{i}, if the two clouds are separated. Between level *i* + 1 and level *i* + *k,* there is no cloud, we can obtain the downwelling flux at level *i* + *k* based on the clear sky radiative transfer calculation. Then again using the approach of (22)–(24), we can obtain the downwelling flux at level *i* + *k* + 1. This is the principle to handle the radiative transfer through random clouds by using two-stream (or higher) scheme. The crucial point for the method is the well mixing assumption for flux in the cloud free layers. This results in no correlation between the upper and lower clouds.

Now we consider an adjacent cloud block occupying several model layers with different cloud fractions in each layer (Fig. 8b or 8c). Similar to (22)–(24), we can obtain the fluxes of ℳ^{−}_{i+1}, 𝒩;th^{−}_{i+1}, and *F*^{−}_{i+1}. The downwelling fluxes at level *i* + 1 for the clear portion and the cloud portion are dependent on the relationship of *c*_{i+1} and *c*_{i}. Figures 8b and 8c illustrate that the downwelling flux at level *i* + 1 for the clear sky portion is ℳ^{−}_{i+1} (for *c*_{i+1} > *c*_{i}) or [(*c*_{i} − *c*_{i+1})𝒩;th^{−}_{i+1} + (1 − *c*_{i})ℳ^{−}_{i+1}]/(1 − *c*_{i+1}) (for *c*_{i+1} ≤ *c*_{i}). Therefore,

where *τ*_{i+1} is the (*i* + 1)th layer optical depth for the clear sky portion. Similarly, the flux for the cloudy portion is

where *τ*_{i+1} is the (*i* + 1)th layer optical depth for the cloudy sky portion. Also we have

By the similar process we can obtain *F*^{−}_{i+k}(*k* = 3, 4, 5, …) for the lower adjacent layers.

In the above, the discussion is based on the random and maximum overlap assumption. The observations and cloud-resolving models suggest that the cloud could be in slantwise overlap (Barker et al. 1999). Figure 8d shows two cloud layers in slantwise overlap. The overlap region for the two adjacent cloud layers is smaller than that for the maximum overlap. For example, let the overlap of *c*_{i} and *c*_{i+1} be *o*_{i,i+1}, *o*_{i,i+1} < Max(*c*_{i}, *c*_{i+1}). Following the scheme discussed above, we can easily deal with the radiative transfer for the slantwise overlap case. The downwelling flux at level *i* + 1 is the same as (22)–(24). For level *i* + 2, we have,

where *τ*_{i+1} is the (*i* + 1)th layer optical depth for the clear sky portion, and

where *τ*_{i+1} is the (*i* + 1)th layer optical depth for the cloudy sky portion. Note, if *o*_{i,i+1} = *c*_{i} (i.e., *c*_{i} < *c*_{i+1}) or *o*_{i,i+1} = *c*_{i+1} (i.e., *c*_{i} > *c*_{i+1}), (28)–(29) become to identical to (25)–(26). Therefore, (28)–(29) are the more general results.

First we test the accuracy of the proposed scheme for a random cloud system. For the same cloud cases as Fig. 6, we consider the low, middle, and high clouds in a column together. Each cloud block is assigned a cloud fraction, *c.* Within a contiguous block, clouds are treated assuming the maximum overlap. For the individual cloud blocks, the random overlap is assumed.

For the benchmark result, we consider *N* calculations by the 1D column model. In each calculation the cloud fraction is randomly set as either 0 or 1 for each of the three cloud blocks (the clouds are identical within a contiguous block). For a cloud block, the chance of the cloud fraction being 1 in *N* calculations is equal to the mean cloud fraction, *c,* for that cloud block. The radiative transfer code is then applied to each column (either clear or overcast) for *N* times and the results are averaged to produce the fluxes and profile of heating rates for the random cloud system. This is the statistical simulation of the random cloud based on independent column approximation (ICA) model method. This method is equivalent to a 1D Monte Carlo simulation, since the results are a statistical average of a random set of clouds. The standard deviation for the 1D Monte Carlo simulation should be proportional to *N*^{−1/2}. Therefore a relatively large number of *N* is required in order to obtain a sufficiently small standard deviation.

Figure 9 shows results for three cloud cases (ran1, ran2, and ran3) where the ICA results based on *N* = 10 000 calculations. The cloud fractions are shown in the figure as the shaded areas. The top panels in Fig. 9 show the cooling rate profiles for three random cloud cases from the benchmark ICA and model calculations. The bottom panels shows the difference in cooling rate between the model results and the benchmark ICA results. It is found that the errors are extremely small, less than 0.1 K day^{−1} for most of cases.

The errors in the upwelling flux at TOA and downwelling flux at the surface are listed in Table 2. The maximum error is 0.5 W m^{−2} (relative error about 0.2%). Figure 9 and Table 2 demonstrates that the physical consideration for dealing with the radiative transfer through clouds with random overlap is correct.

The results in Fig. 9 and Table 2 are very accurate. The domain mean radiation results can be achieved by only a single calculation regardless how complicated the cloud field is. This is the main point of this paper. We will further explore this point in the following discussion for adjacent clouds and especially the cloud subgrid variability in Part II.

For the maximum cloud overlap, we consider five cloud overlap cases shown in Fig. 10 (the shaded areas indicate the cloud configurations). Cases A–D are cloud configurations with the maximum overlap; case E is a cloud configuration with the slantwise overlap.

Figure 10 shows the comparison results of the cooling rate. All results are against the mean results of ICA. In a maximum overlap case, only a small number ICA calculations are needed, corresponding to different regions of the cloud configuration. It is found that this method can handle the radiative transfer through maximum overlapping fractional clouds accurately. In Fig. 10, the cooling rate profiles could be totally different for different cloud overlap configurations, even the total cloud fractions are the same. This indicates the way of dealing with cloud overlap is a very important issue.

In Fig. 10 the worst cases always occur in cloud case D. For the low and middle clouds, the errors in the cooling rate are larger in comparison with other cloud cases. In cloud case D, the cloud configuration is concave. Inside the concave region the downwelling (upwelling) flux from the cloud portion would not well mix with the downwelling (upwelling) flux from clear sky portion. This causes the relatively large errors in cooling rate, since well mixing is the basic assumption for this method. However, in the real atmosphere, such a concave cloud geometry seldom happens. Just for testing the validation of our method, we choose cloud configurations with extremes in variation.

Table 3 listed the results for the upwelling flux at the TOA and the downwelling flux at the surface. Generally the errors in flux are very small, less than 3 W m^{−2} for most cases. The relatively large errors also occur for cloud case D. The error could be 12.4 W m^{−2}, about 3% in relative error.

## 5. Discussion and summary

GCM intercomparisons (e.g., Cess et al. 1989) show that the results of radiative transfer calculations for clear sky cases for different GCMs are always very consistent with one another. However, large differences among the models occur for cloudy sky cases. This indicates that the results are very sensitive to the radiative transfer treatment and that there may be significant differences between models.

By comparison with solar radiative transfer, the infrared radiative transfer is more simple. This is because cloud scattering is a second-order effect in the IR. Based on AA, the upwelling and downwelling fluxes are decoupled, and the upwelling and downwelling calculation paths can be treated separately. This dramatically simplifies the calculation process and makes it possible for us to handle the radiative transfer through clouds accurately in a GCM.

In this work, the infrared radiative transfer solution for a multilayer system is solved in a way that all of the radiation contributions above and below a given layer are represented as boundary conditions for that layer. It is found that this simplifies the calculation process. Based on this approach, we have systematically investigated three different solutions for AA corresponding to three different assumptions for the Planck source function. It is shown in appendix A that the E source scheme matches better with the linear temperature distribution compared with the L source scheme. However, we found that if the model layer thickness is not too large (e.g., <0.25 km), then the differences in cooling rates and fluxes are very small for the L and E source schemes. Both schemes could provide accurate results in the clear sky and the cloudy sky conditions. It is also noted that the calculations and analysis show that a relatively coarse resolution can be applied to the higher regions of an atmospheric profile. We find that the error associated with the change in flux across a cloud layer is less than 0.1 W m^{−2} for the L and E schemes regardless of the values of LWC or IWC.

A new approach to the perturbation method for solving cloud IR radiative transfer including scattering effects has been proposed. It is more simple in its form and allows for a more consistent treatment of radiation through cloudy and clear layers. In addition, aerosols in cloud-free layers could be handled more accurately with the zero-order equation implicit in AA. The perturbation method for cloud scattering is very efficient and only one more calculation path is needed in comparison with the AA (which does not account for scattering effects).

Cloud overlap is a significant problem because it can have a large impact on the radiative transfer solution for a column. In this work, a new approach for the radiative treatment for random and maximum cloud overlap is presented. This cloud overlap treatment is designed to be consistent with and easy to implement for the two path calculation process (three paths if perturbation method is used for scattering) typical in GCM column calculations. Compared with the matrix method, this new method is physically intuitive and computationally more efficient. To test the scheme, Monte Carlo calculations have been used. The Monte Carlo method has been used to create random cloud fields by applying large a number of individual 1D column model calculations with the cloud layer fraction 0 or 1. This enables us to generate the benchmark results for comparison of radiative transfer for random clouds. We also compared our model with ICA for the results of maximum overlap and slantwise overlap for different cloud configurations. The error in the cooling rates and radiative fluxes are generally less than 1 K day^{−1} and 3 W m^{−2}, respectively. All discussions in this paper are based on the two-stream method; however, it is straightforward to extend them to the higher stream case by Li (2000b).

Although accurate results have been achieved concerning the cloud scattering and cloud overlap problems, the clouds treated are ideal in that they are assumed to be horizontally homogeneous within the cloud. This is a “classical” way to describe cloud in climate models. This simplified cloud treatment can be traced back to the early work of Manabe and Strickler (1964). In Part II, we will present a way to treat the within-cloud horizontal variability of cloud water path for radiation calculations.

## Acknowledgments

The author would like to thank Drs. H. W. Barker, S. Dobbie, and K. v. Salzen, and two anonymous reviewers for their helpful comments. The author would especially like to thank Dr. H. W. Barker for his suggestion of using a 1D Monte Carlo method to create the random cloud fields.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

_{2}15μ band cooling rates.

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Evaluation of Thermal Source Function

Wiscombe (1976) presents a relationship between the L source function and the ambient temperature distribution. In the following, an intuitive picture will be presented showing the temperature deviation from the assumed linear interpretation of temperature in a model layer for both L and E source schemes.

The Planck function is

where *λ* is the wavelenth, *T* is the temperature. Here *c*_{1} = 2*πc*^{2}*h, **c*_{2} = *hc*/*b, **c* is the speed of light, *b* is the Boltzmann constant, *h* is the Planck constant.

Wiscombe (1976) assumes that the extinction coefficient exponential in height. Similarly Chou and Arking (1980) parameterize the extinction coefficient as

where *k* and *k*_{0} are the extinction coefficient at pressure *p* and surface *p*_{s}, respectively. Chou and Arking (1980) shows that the above relation holds true at least for troposphere and lower stratosphere. As shown in section 2 of this paper, the radiative transfer is not sensitive to the source function for the higher atmosphere. Based on the above assumption for extinction coefficient we derived

where *H* is the scale hight, *z*_{i} (*τ*_{i}) is the thickness (optical depth) for layer *i,* and *z* (*τ*) is the distance (optical depth) starting from level *i.* For the L source scheme, the Planck function inside layer *i* is

where *B*_{i} and *B*_{i+1} are the Planck functions for the temperatures at the top and the bottom of layer *i.* From (A1) and (A3),

where *T*_{i} and *T*_{i+1} are the temperatures at levels *i* and *i* + 1, respectively. The upper panels in Fig. A1 show the temperature variation within layer *i* by (A4) for *λ* = 4, 10, and 100 *μ*m. Here 4 and 100 *μ*m are close to the endpoints of the infrared spectrum, while 10 *μ*m is within the atmospheric window region. The dotted lines in Fig. A1 show the results for a linear distribution of temperature inside the layer. It is found that for a large value of *λ,* the variation of temperature by (A4) matches with the linear distribution results. This can be verified, since for a large value of *λ, *(A4) becomes

Since *z*/*H* and *z*_{i}/*H* are small, we approximately have

This is a linear distribution in temperature.

For a small value of *λ, *(A4) becomes

The exponential function yields results that deviate from the linear distribution for large differences in the two-level temperatures. Therefore, if the temperature is supposed to be linear distributed inside the layer, this assumption holds true for a large value of *λ.* For small values of *λ,* the results are accurate only for temperature differences less than 10 K (see Fig. A1). Also, by comparing the top and the third up from the bottom panels in Fig. A1, we find that the deviation of temperature from the linear distribution becomes larger as the mean temperature of the layer decreases.

For the E source scheme,

by (A1),

The second level and the bottom panels in Fig. A1 show the temperature variation within layer *i* for *λ* = 4, 8, and 100 *μ*m. In comparison with the upper panels, generally the E source scheme produces layer temperature profiles very close to the linear distribution results. This can be verified, since for a large value of *λ, *(A9) becomes

which reduces to the linear distribution of (A5) when the difference between *T*_{i} and *T*_{i+1} is small. For a small value of *λ,*

In (A11), this function is linear in inverse temperature. Usually, (*T*_{i+1} − *T*_{i})/*T*_{i+1} is small, and the first order of the expansion of (A11) is

This is why the E source scheme can always match with the linear temperature results for all of the spectrum range.

In summary, a certain source function should correspond to a certain distribution of temperature in a layer. If a linear distribution in temperature is assumed, generally the E source scheme matches better with the assumption as compared with the L source scheme, especially for small values of *λ.* Both the L and E source schemes match well with the linear distribution for a layer temperature difference of less than 10 K.

Actually *B*[*T*(*τ*)] should be determined by the distribution of temperature directly. However, the form of *B*[*T*(*τ*)] would be very complicated even for the linear distribution of (A6), and the corresponding analytical solution does not exist. We therefore have to deal with the problem in a reverse way by assuming *B*[*T*(*τ*)] to be a very simple form, with which the analytical solution of AA exists. This appendix shows the reasonableness for the assumptions of L and E source schemes judged by the linear temperature distribution within a layer.

### APPENDIX B

#### Perturbation Method for Scattering

We can write the radiative transfer equation as

since the cloud absorption is usually strong in the infrared, and the single-scattering albedo *ω̃* is small, we can expend the radiance in powers of *ω̃*, *I*(*τ, **μ*) = *I*^{(0)}(*τ, **μ*) + *ω̃**I*^{(1)}(*τ, **μ*) + ···. Thus a sequence of equations for the outer expansion is

We only considered equations up to the first order. This is the standard perturbation method. Actually we can expand in a power of a smaller coefficient, like *ω̃*/2, to make the solution more accurate.

The first-order correction is supposed to only apply to layers with cloud present (for saving computing time). For aerosols outside of cloud layers, the infrared scattering effect is very small; however, the absorption can be strong. The zeroth-order equation of (B2a) is not suitable for dealing with the radiative transfer for an absorbing medium since it contains no single-scattering coalbedo factor. However, AA is able to handle the radiative transfer through aerosol properly (Li and Min 2001). Therefore, it is better to expand the equation with the zero-order equation as in AA. If we let the transfer equation be

with *I*(*τ, **μ*) = *I*^{(0)}(*τ, **μ*) + *ω̃**I*^{(1)}(*τ, **μ*), we have

where *ε* = 1 − *ω̃*. Equation (B4a) is now is the AA equation.

For the L source scheme the upwelling radiance and downwelling radiance between levels *i* and *i* + 1 for the zero-order equation (B4a) are

where *I*^{(0)}_{i+1}(*μ*) and *I*^{(0)}_{i}(−*μ*) are the upwelling radiance and downwelling radiance at levels *i* + 1 and *i,* respectively; *α*_{i} = *B*_{i+1} − *B*_{i}, with *B*_{i} and *B*_{i+1} being the Planck functions for the temperature at levels *i* and *i* + 1, respectively; *τ* is the optical depth starting from level *i*; *τ*_{i} is the optical depth between levels *i* and *i* + 1; and *κ*_{i} = *ε*_{i}*τ*_{i}. Similarly we can obtain the solution for the E source scheme.

We do not use the phase function of (2) directly, instead, we use

with scalings for *τ, **ω,* and *g* following Joseph et al. (1976).

In the two-stream case, we have (Li and Fu 2000)

where 1/*μ*_{1} = 1.647 213.

## Footnotes

*Corresponding author address:* Dr. Jiangnan Li, Canadian Center for Climate Modeling and Analysis, Meteorological Service of Canada, P.O. Box 1700, University of Victoria, Victoria, BC, Canada V8P 2Y2. Email: Jiangnan.Li@ec.gc.ca