## Abstract

Proper quantification of the solar radiation budget and its transfer within the atmosphere is of utmost importance in climate modeling. The delta-four-stream (DFS) approximation has been demonstrated to offer a more accurate computational method of quantifying the budget than the simple two-stream approximations widely used in general circulation models (GCMs) for radiative-transfer computations. Based on this method, the relative improvement in the accuracy of solar flux computations is investigated in the simulations of the third-generation Canadian Climate Center atmosphere GCM. Relative to the computations of the DFS-modified radiation scheme, the GCM original-scheme whole-sky fluxes at the top of the atmosphere (TOA) show the largest underestimations at high latitudes of a winter hemisphere on the order of 4%–6% (monthly means), while the largest overestimations of the same order are found over equatorial regions. At the surface, even higher overestimations are found, exceeding 20% at subpolar regions of a winter hemisphere. Flux differences between original and DFS schemes are largest in the tropics and at high latitudes, where the monthly zonal means and their dispersions are within 5 W m^{−2} at the TOA and 10 W m^{−2} at the surface in whole sky, but differences may be as large as 20 and −40 W m^{−2}. In clear sky, monthly zonal means and their dispersions remain within 2 W m^{−2}, but may be as large as 25 and −12 W m^{−2}. Such differences are found to be mostly determined by variations in cloud optical depth and solar zenith angle, and by aerosol loading in a clear sky.

## 1. Introduction

The earth climate system and all processes within it are principally governed by the balance between incoming solar radiation and outgoing terrestrial radiation. Insofar as the climate system may be regarded as one massive thermodynamic engine fuelled by solar radiation, proper quantification of the solar radiation budget, particularly the transfer of this radiation within the atmosphere, is of utmost importance in climate modeling. The importance of precisely computing the radiative fluxes in climate models can hardly be overemphasized, given that a radiative inaccuracy of even a few W m^{−2} can seriously compromise the evaluation of heating rates and climate forcings (Fouquart et al. 1991).

Several attempts have been made to evaluate the accuracy of the various solar radiation transfer schemes used in climate models. The Intercomparison of Radiation Codes in Climate Models (ICRCCM) international program (Fouquart et al. 1991) showed that, with the same atmospheric profiles, the range of irradiances estimated by the different codes often exceeded 20% (with 4%–10% typical root-mean-square differences relative to medians), with the largest disparities occurring in aerosol-laden and cloudy atmospheres. Barker et al. (2003) assessed the performance of solar radiative transfer codes used both for research and in weather and climate models, with emphasis on handling unresolved clouds. Barker et al.’s assessment shows that even though most codes since the ICRCCM have been modified to better resolve the solar spectrum, spectroscopic databases have been updated, and new gaseous transmittance and cloud optical property parameterizations have been implemented, most codes still underestimate atmospheric absorption of solar radiation even for clear-sky and homogeneous overcast cloud cases. While mainly assessing the performance of research-type solar radiative transfer codes, Halthore et al. (2005) point out that recent studies have shown considerable discrepancy between modeled and measured solar radiation components in both clear and cloudy skies, indicating that large uncertainties may be unaccounted for in measurements or that possible absorption mechanisms may be either underrepresented or perhaps completely absent in radiative transfer models.

Quantifying the shortwave (SW) radiation budget in climate models is tantamount to solving a problem of radiative transfer in planetary atmosphere. Given that the exact solution of the radiative transfer equation in a scattering and absorbing medium is difficult to obtain in a computationally efficient manner even under a plane-parallel assumption, approximate methods are necessary (Li and Ramaswamy 1996). In radiative transfer approximations, the higher the order of stream expansions, the more accurate the approximation becomes. However, such a gain in accuracy comes at the expense of a more computationally intensive scheme. In general circulation models (GCMs), where radiative transfer calculations are performed over all grid cells, vertical levels, radiation bands, and clear/cloudy conditions, computational economy is a must. Consequently, shortwave radiative transfer schemes in GCMs often resort to the least numerically expensive methods, with the two-stream and Eddington approximations being the most common choices, as they provide simple yet accurate analytical solutions that are feasible to incorporate into the GCM radiative schemes.

More accurate methods for calculation of radiative fluxes in planetary atmosphere have been developed based on four-stream approximations (Cuzzi et al. 1982; Li and Ramaswamy 1996; Li and Dobbie 1998; Liou 1974; Liou et al. 1988). Liou et al. (1988) formulated a solution of the delta-four-stream (DFS) approximation in analytic form, which makes it suitable for computation of solar radiative fluxes with high accuracy and computational economy. They demonstrated the superiority of DFS against the delta-two-stream (DTS) method through single-layer computations involving an exact method with a wide range of input radiative parameters. The accuracy of this method was investigated by Shibata and Uchiyama (1992) and by Chou (1992), who compared the DFS flux computations to “exact” values in a multilayered atmospheric-column model that included absorbing gases, aerosols, and clouds. Their results show that fluxes and heating rates are accurately computed to within a few percent by the DFS method. Kay et al. (2001) compared the accuracy of the two-stream, the DFS, and the matrix inversion methods relative to the exact discrete ordinates (32-stream) method by using a multilayered atmospheric-column model with different surface/aerosol systems and cloud layer. Their results show that the DFS approximation is the technique that stands out in terms of accuracy and efficiency. Furthermore, the DFS method has been used by Fu and Liou (1993) to provide accurate computations of the radiative properties and effects of cirrus clouds, while Christopher et al. (2000) have included the DFS method in a radiative transfer model to estimate the radiative forcing of biomass-burning aerosols.

A solar radiation transfer code has been developed for use in the third-generation Canadian Climate Center Atmosphere GCM (CCC GCM III), which is based on the DFS method. While the “absolute” accuracy of this method has been assessed for typical model atmospheres, the main objective of this study is to assess the relative improvement in the accuracy of solar radiation computations in GCM simulations using the DFS method. The model used and the simulation setup are described in section 2. Results of the simulations performed for this study are presented in section 3, ending with some discussion and conclusions in section 4. Specific aspects pertaining to the implementation of the DFS method in CCC GCM III are discussed in the appendix.

## 2. Model and simulation description

As the DFS method seems to offer a compromise between accuracy and computational efficiency, a shortwave transfer code based on this method has been developed and integrated into the solar radiation scheme of the third-generation Canadian Climate Center Atmosphere GCM (CCC GCM III). In many aspects, the solar radiation scheme of the third-generation GCM is still similar to that of the second-generation GCM (McFarlane et al. 1992). In this scheme, solar radiation computations are done every hour over four bands in the visible and near-infrared region (0.25, 0.69, 1.19, 2.38, and 4.00 *μ*m), accounting for the absorption and scattering by gases, aerosols, and clouds. Given the resultant optical properties of each model layer, the reflectance and transmittance are computed using two different approaches: (i) the two-stream approximation (Coakley and Chylek 1975) for optical mixtures of aerosols and gases, and (ii) the delta-Eddington approximation (Joseph et al. 1976) for optical mixtures of gases, aerosols, and clouds. The first calculation (two stream) is always done for purposes of computing the cloud radiative forcing. In case clouds are present in a layer, the overall reflectance and transmittance are computed as a cloud-coverage-weighted mean between the clear and the cloudy fractions of each layer. These reflectances and transmittances are then combined upward (from the layer in question to the top of the atmosphere) and downward (from the layer in question to the surface), resulting in the final reflectances and transmittances of each model layer. These are used to evaluate the solar fluxes across model levels by multiplication with the top-of-atmosphere (TOA) solar flux and the ozone transmission function. Absorption and scattering by principal aerosol species (sea salt, sulfate, black and organic carbon, and soil dust) are accounted for by coupling to the CCC GCM the Canadian Aerosol Module (CAM; Gong et al. 2003) and using the aerosol optical parameters module described by Ayash et al. (2008).

The original shortwave radiation scheme in CCC GCM III is modified by using the DFS code for both clear- and whole-sky computations. The results presented here were obtained from a two-year model run by making parallel calls (at every model hour) to the shortwave radiation routine in the same model run: first with the original codes (two stream and delta Eddington) with fluxes fed back into the model, and then with the DFS modified code with fluxes only for diagnosis. This approach is adopted to compare original and modified radiative calculations against the exact same input.

## 3. Results

The changes in the simulated flux brought about by using the DFS code are examined here. Monthly zonal means of the differences between original and modified daily mean fluxes, along with their limits and their standard deviations, are presented in Figs. 1 and 2 for the months of January and July. Given that the DFS method represents an approximation that is closer to the exact method than the delta-Eddington and two-stream methods, it would be useful to examine not only the flux differences but also the differences in flux computations relative to the “more exact” one by defining the percentage difference (PD) as

Flux differences and PD are examined for the SW fluxes outgoing (upward) at the TOA and downwelling (downward) at the surface, since they signify the net effect of all computations performed within the model’s atmospheric vertical levels and are the principal parameters considered in any radiative climate analysis. The sign convention followed here is that a flux is always positive, regardless of its direction.

According to several studies on radiative transfer (Cuzzi et al. 1982; Liou et al. 1988; Stephens et al. 2001), the accuracies of parameterizations are most sensitive to the parameters: optical depth, single-scattering albedo (SSA), and cosine of solar zenith angle (CSZA). Radiative properties of the atmosphere are largely influenced by the presence of clouds, while, in a clear sky, atmospheric radiative properties are influenced by variations in the aerosol loading. Therefore, the monthly zonal means of the column-integrated optical depths of clouds and aerosols simulated by our models are included in Figs. 1 and 2, respectively. Another good measure of the spatial distribution of clouds—the cloud radiative forcing (CRF), which is computed as the clear-sky minus whole-sky outgoing SW fluxes at the TOA (Harrison et al. 1990)—is shown in Fig. 5.

### a. Whole-sky fluxes

We find that our results may be characteristically explained by Liou et al.’s (1988) calculations of the accuracies of the DTS and the DFS approximations relative to the exact adding method. For an SSA of 0.8, their single-layer computation results show that

for the layer reflection, the DTS exhibits higher underestimations with decreasing CSZA (up to 30%) compared to the DFS (up to 10%) below a CSZA between 0.3 and 0.5, but otherwise shows higher overestimations with increasing optical depth (up to 15%) compared to the DFS (up to 5%);

for the layer total transmission, the DTS shows higher overestimations with decreasing CSZA (up to 30%) compared to the DFS (up to 10%) below an optical depth of about 5.0, but otherwise shows much higher underestimations with increasing optical depth (up to 90%) compared to the DFS (within 5%).

At the TOA, upward flux is determined by the resultant reflection across model layers. In line with Liou et al.’s (1988) findings for reflection, lower CSZA values at high latitudes in both hemispheres that remain within 0.75 result in average underestimations of upward flux, with mean differences and standard deviations within 5 W m^{−2} in the respective summer hemispheres and with PD values on the order of 4%–6%. At equatorial and tropical latitudes, a combination of higher CSZA values (within unity) and maximum cloud optical depths leads to average overestimation of upward flux with about the same magnitude (but with opposite sign) of flux underestimation at high latitudes. This is further confirmed by the spatial correlation between areas of higher mean PD values and those of high CRF, such as the regions of South America and South Africa in January, and the tropical west Pacific, East Asia, and India in July (monsoon season).

At all latitudes, however, TOA flux differences may attain higher positive and negative values, with magnitudes as high as 15 W m^{−2}. Flux differences at polar latitudes differ between the hemispheres in that these differences come to a null average at the highest latitudes in the SH, while consistently maintaining an average underestimation in the NH in July. This may be explained by the significantly lower cloud optical depths at the SH highest latitudes compared to the NH ones, which reflect the typically dry atmosphere of the Antarctic. Under such conditions, GCM fluxes would be obtained mostly by clear-sky computations, to which the two-stream approximation is well suited.

At the surface, average differences in downward flux are largely positive, with the largest overestimations occurring at high latitudes in both hemispheres. With the downward fluxes determined by the resultant transmission across model layers, such findings are also in line with Liou et al.’s (1988) findings for layer transmission, given that high latitudes are associated with lower CSZA values. While the largest average differences and standard deviations are found at the high latitudes of the respective summer hemispheres, the opposite holds for the highest PD values. This is because CSZA values are lowest over the winter hemisphere and hence exhibit the largest relative differences and PD values, while solar flux is larger over the summer hemisphere, leading to larger differences even with lower PD values. Means and standard deviations of the flux differences at the surface are within 10 W m^{−2}, but these differences may be as high as 20 W m^{−2}, and may also attain large negative values with high cloud optical depth and solar flux, reaching −40 W m^{−2} at northern midlatitudes for the month of April (not shown here). Such large negative differences are also in line with Liou et al.’s (1988) findings for high optical depths. On the other hand, surface PD attains particularly large values, exceeding 20% at subpolar regions.

### b. Clear-sky fluxes

Differences in simulated fluxes have also been investigated for clear-sky conditions, with and without the optical properties of aerosols. With aerosols excluded, negligible PD values are found, remaining within 1% at all times and locations. Such a finding demonstrates that the DFS computations do in fact converge to those of the two-stream method under thin-atmosphere conditions, where the latter method’s superiority to other two-stream formulations has been demonstrated by Coakley and Chylek (1975). The inclusion of aerosols is found to lead to significant changes in PD values, with mean overestimations on the order of 4%–6% at the TOA and underestimations exceeding 10% at the surface over the subpolar regions in January and July. The effect of aerosols is illustrated in Fig. 2, which shows that at northern tropical and subtropical latitudes in April, the higher aerosol optical depths that are largely caused by the high soil-dust loadings of the desert regions at these latitudes lead to high flux differences, both at the TOA and at the surface. While average differences and standard deviations remain within 2 W m^{−2} for all months, these differences may be as high as 25 W m^{−2} at the surface and −12 W m^{−2} at the TOA.

Such large differences in clear-sky fluxes may be attributed in part to the delta adjustment in the DFS method, which is absent in the two-stream method. Because of the fact that scattering by atmospheric particulates is highly peaked in the forward direction, such that a highly peaked diffraction pattern is typical for atmospheric particles, it has been shown that a delta-function adjustment offers significant enhancement in flux computations (Liou 2002). Another possible reason may lie in the fact that the solution of the two-stream method employed in the GCM radiation code is that of a thin-atmosphere approximation, which is mostly valid for typical background aerosol profiles. However, for certain conditions of exceptionally high aerosol loadings, such as the dust-storm episodes that are simulated by CAM, the thin-atmosphere assumption may become a poor approximation for flux calculations.

Finally, we note that variations in the monthly zonal mean of cloud SSA (ranging from 0.4 to 0.9) and aerosols (found between 0.9 and 1.0) could not help explain the variations in flux differences. Therefore, our results suggest that, for monthly averages in a global climate model, the accuracy of solar flux computation is largely determined by values of optical depth and CSZA.

### c. Computational time

The computational time of using the DFS code within the CCC GCM III running platform was measured. It was found that the ratio of the time required by the modified SW radiation routine relative to the original routine is anywhere between two to three, depending primarily on the extent of cloudiness within a model vertical column. It should be noted that in the original routine the clear-sky (two stream) computations are performed only once, and the resulting reflectance and transmittance are transferred to the whole-sky computations by using “equivalenced” arrays. In implementing the DFS code, however, maintaining such a scheme led to unexpected numerical errors, which could be avoided only by repeating the clear-sky DFS computations for whole-sky computations. This probably increased the computational time demand by the modified routine, which otherwise could have been less than the figure reported here.

## 4. Discussion and conclusions

The accuracies of four-stream approximations have been investigated in several studies, each looking differently into atmospheric and radiative properties and employing different four-stream formulations. The formulation we use here is that of the DFS by Liou et al. (1988), who suggest that the DFS significantly improves the accuracy of single-layer computations of reflection and transmission compared to the delta-two-stream approximation. Our results here for global-model flux computations are in line with these findings in showing that computations by two-stream and delta-Eddington approximations may significantly differ relative to DFS computations, on the order of 6% at the TOA and over 20% at the surface. Using the same DFS formulation in a multilayered atmospheric-column model, Kay et al. (2001) show that net-flux computations by the two-stream method (with a delta-function adjustment) may differ from those of the DFS by 6–10 W m^{−2} for clear-sky conditions including oceanic and terrestrial aerosols, and by 5–9 W m^{−2} when a single cloud layer is added. While our results for monthly zonal averages of flux differences are within these numbers, our model simulations show, however, that flux differences may also be significantly larger, which may be attributed to the distribution of clouds in multiple layers in actual GCM simulations and the occurrence of aerosol loadings that are higher than the average profiles considered by Kay et al. (2001).

On the other hand, contrary to findings by this and other studies, Stephens et al. (2001) find that the use of a four-stream method offers little gain over two-stream methods. One possible reason for such a finding is that a four-stream doubling model was used, but with no delta adjustment (G. L. Stephens 2007, personal communication). The inclusion of such an adjustment in a four-stream model has been shown to offer some improvement in the overall consistency and accuracy of flux computations, and the four-stream model without a delta adjustment can be less accurate than the two-stream model with a delta adjustment, such as for reflection at low optical depths and small CSZA (Cuzzi et al. 1982). In any case, Stephens et al. (2001) base their findings on single-layer computations of reflection and transmission, which do offer the advantage of classifying atmospheric radiative conditions based on the optical depth and SSA encountered in actual GCM simulations, but do not examine flux computations in a multilayered GCM atmosphere.

Our analysis is based only on relative differences in flux computations between two-stream/delta-Eddington and DFS methods, since it would be practically impossible to include a numerically expensive exact computational method in long-term GCM runs. However, Chou (1992), who examined the accuracy of Liou et al.’s (1988) DFS formulation against detailed exact calculations in a multilayered atmospheric-column model, found that the accuracy of DFS computations is within 7.5 W m^{−2} over ultraviolet and visible bands and is within 4 W m^{−2} over near-infrared bands.

Overall, our results show that the use of the more accurate DFS method in the CCC GCM III solar radiation scheme yields significant changes in the radiative flux computations. Such changes are most significant in the tropics and at high latitudes where the monthly zonal means of flux differences and their dispersions for whole-sky computations are within 5 W m^{−2} at the TOA and 10 W m^{−2} at the surface (although differences may be as large as 20 and −40 W m^{−2}). For clear-sky computations, monthly zonal means of flux differences and their dispersions remain within 2 W m^{−2}, but may be as large as 25 and −12 W m^{−2}. Such differences are found to be mostly determined by variations in cloud optical depth and solar zenith angle, and by aerosol loading in a clear sky. The implications of such significant changes in solar flux computations for the climate dynamics simulated by the GCM would be worth investigating. Furthermore, the improvement of the overall accuracy of GCM flux simulations by using the DFS approximation is also worth investigating against observational data, which would involve closure experiments between modeled fluxes and satellite flux measurements and accounting for radiatively significant variables, such as cloud cover and surface albedo.

## Acknowledgments

We thank J. Li for his valuable input.

## REFERENCES

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### DFS Implementation in CCC GCM III

The DFS code output should be equivalent to that of the delta-Eddington and two-stream codes in the shortwave radiation routines of CCC GCM III.

The delta-Eddington computations provide each model layer’s reflectance and transmittance, with and without reflection from the lower boundary (which could be the actual ground surface or the upper boundary of the underlying layer). Therefore, by setting the reflectivity of the surface underlying the layer first to the underlying albedo and then to zero, the boundary value and intensity computations are performed twice, resulting in two sets of the layers’ reflectance and transmittance calculated as

where

Note that with the diffuse intensity terms being proportional to *F*_{Θ} this term is set to 1 in Eq. (A1).

The two-stream computations result in different sets of output, which are the layer’s reflectance and transmittance with multiple reflections and with single reflection from the lower boundary. While the former set of output corresponds to the DFS’s output with reflection, only the latter’s transmittance corresponds to the DFS’s transmittance without reflection. Therefore, an additional quantity needs to be computed, which is the layer’s reflectance with only single reflection from below. This may be realized by considering the two-stream equations for a layer’s reflectance and transmittance with multiple reflections from below (Liou 2002), as shown here:

where

The required quantity, that is, reflectance with single reflection, is given by

The spherical transmission (*γ*) is not part of the DFS’s computations and, hence, needs to be expressed in terms of other quantities that are. Rearrangement of Eqs. (A2) and (A3) yields

which gives the desired quantity in terms of computed ones in the DFS code.

## Footnotes

*Corresponding author address:* Tarek Ayash, Dept. of Chemical Engineering and Applied Chemistry, University of Toronto, 200 College Street, Toronto, ON M5S 3E5, Canada. Email: tarek.ayash@utoronto.ca