## Abstract

Solar flux densities and heating rates predicted by a broadband, multilayer δ-Eddington two-stream approximation are compared to estimates from a Monte Carlo model that uses detailed descriptions of cloud particle phase functions and facilitates locally nonzero net horizontal flux densities. Results are presented as domain averages for 256-km sections of cloudy atmospheres inferred from A-Train satellite data: 32 632 samples for January 2007 between 70°S and 70°N with total cloud fraction *C* > 0.05. The domains are meant to represent grid cells of a conventional global climate model and consist of columns of infinite width across track and Δ*x* ≈ 1 km along track. The δ-Eddington was applied in independent column approximation (ICA) mode, while the Monte Carlo was applied using both Δ*x* → ∞ (i.e., ICA) and Δ*x* ≈ 1 km. Mean-bias errors due to the δ-Eddington’s neglect of phase function details and horizontal transfer, as functions of cosine of solar zenith angle μ_{0}, are comparable in magnitude and have the same signs.

With minor dependence on cloud particle sizes, the δ-Eddington over- and underestimates top-of-atmosphere reflected flux density for the cloudy portion of domains by ~10 W m^{−2} for μ_{0} > 0.9 and −3 W m^{−2} for μ_{0} < 0.2; full domain averages are ~8 and −2 W m^{−2}, respectively, given mean *C* > 0.75 for all μ_{0}. These errors are reversed in sign, but slightly larger, for net surface flux densities. The δ-Eddington underestimates total atmospheric absorption by ~2.5 W m^{−2} on average. Hence, δ-Eddington mean-bias errors for domain-averaged layer heating rates are usually negative but can be positive. Rarely do they exceed ±10% of the mean heating rate; the largest errors are when the sides of liquid clouds are irradiated by direct beams.

## 1. Introduction

With near unanimity, global climate models (GCMs) employ multilayer two-stream approximations to compute solar radiative fluxes (e.g., Coakley and Chýlek 1975; Meador and Weaver 1980; Zdunkowski et al. 1980).^{1} The varied abilities of two streams to handle scattering by cloud droplets and ice crystals (see King and Harshvardhan 1986) has led to those that employ delta scaling of optical properties (Joseph et al. 1976; Wiscombe 1977) being the models of choice for GCMs. Paralleling over three decades of widespread usage of two streams in GCMs, numerous studies have demonstrated the shortcomings of 1D solutions of the solar radiative transfer equation in general, with particular emphases on common two-stream approximations. The focus has been on demonstrating their limitations at representing either details of particulate optical properties (e.g., King and Harshvardhan 1986; Li and Ramaswamy 1996; Räisänen 2002) or their neglect of horizontal radiative transfer for media characterized by fluctuations at scales less than a few kilometers (Marshak and Davis 2005b; Hogan and Shonk 2013). As yet, however, these two classes of errors have not been assessed on a head-to-head basis.

The purpose of this study is to provide estimates of flux and heating rate biases that can be expected from two-stream-based multilayer broadband solar transfer codes, as used in GCMs, due to their necessary neglect of both details in scattering phase functions for cloud droplets and ice crystals and horizontal transport of photons at spatial scales below the resolution of GCMs. This was achieved using cloud properties inferred from *CloudSat* (Stephens et al. 2002) and *CALIPSO* (Winker et al. 2003) data, along with ECMWF state variables. The two-stream solar code employed here was that used in the Canadian Centre for Climate Modelling and Analysis (CCCma) GCM (von Salzen et al. 2013). A 3D Monte Carlo photon transport algorithm (Barker et al. 2012), which is effectively an infinite-stream approximation, provided benchmark estimates of fluxes and heating rates. Both transfer models used the CCCma’s clear-sky optical properties, as well as common descriptions of liquid and ice cloud optical properties, derived from the Lorenz–Mie scattering theory (Wiscombe 1980) and a combination of techniques (Yang et al. 2013), respectively.

*CloudSat*–*CALIPSO* cross sections, which consist of columns taken to be infinitely wide across track and Δ*x* = 1 km along track, were partitioned into 256-km domains to approximately represent GCM cells. The CCCma two-stream model was applied to each 1-km column and averaged over 256 km, thus affecting the independent column approximation (ICA) (Stephens et al. 1991). Each 256-km domain was then acted on by the Monte Carlo model assuming various phase functions, as well as both Δ*x* → ∞ (i.e., ICA mode) and Δ*x* = 1 km. While the two-stream approximation could utilize only the asymmetry parameter of the particle scattering phase functions, the Monte Carlo used the exact phase functions. It is noted here that, because horizontal variability of cloud in the across-track direction was neglected, errors due to the ICA are underestimated relative to use of fully 3D clouds, potentially by up to ±30% for some clouds under some illumination conditions (cf. Barker 1996; Pincus et al. 2005). Domain-average fluxes and heating rate profiles predicted by the models were compared for data from January 2007.

In its native configuration, the CCCma two-stream approximation executes in the Monte Carlo ICA (McICA) framework (Pincus et al. 2003), as it lacks access to explicit distributions of unresolved clouds. Errors from the McICA method arise through limitations in stochastic generators used to conjure up subgrid-scale clouds (e.g., Räisänen et al. 2004). These errors were not addressed here, so this study was restricted to assessment of limitations in two-stream approximations.

The following section provides brief overviews of the radiative transfer models and cloud optical properties used in this study. The third section shows results for single-layer clouds irradiated by monochromatic, collimated radiation. This helps set the stage for the more complicated broadband results for A-Train atmospheres, which are shown in section 4. A summary and conclusions are drawn in the final section.

## 2. Radiative transfer models

The primary purpose of this study was to furnish estimates of bias errors in solar fluxes computed by two-stream approximations as a result of their neglect of details in cloud scattering phase functions and cloud geometry. The two-stream approximation was exemplified here by the popular δ-Eddington approximation (Joseph et al. 1976), which is used in the CCCma GCM (von Salzen et al. 2013). Zdunkowski et al.’s (1980) practical improved flux method was also used, but differences between it and the δ-Eddington were negligible. The two-stream ICA and the Monte Carlo model both used the same clear-sky optical properties as produced by the CCCma solar radiation code’s correlated *k*-distribution method, which has 31 quadrature points in cumulative probability space (Li and Barker 2005). They also used the same cloud optical properties. This section’s three subsections recapitulate the δ-Eddington and Monte Carlo algorithms as well as cloud optical properties.

### a. δ-Eddington approximation

Let be the azimuthally averaged, normalized phase function that describes the probability of radiation incident on a particle from direction being scattered into direction . By the addition theorem for spherical harmonics,

where are expansion coefficients, and are Legendre polynomials of order *l*. In addition to a layer’s optical depth and single-scattering albedo , many two-stream approximations also utilize the phase function’s asymmetry parameter, which is defined as

The δ-Eddington approximation (Joseph et al. 1976) approximates with

where is the asymmetry parameter of , *f* is the fraction of radiation taken to be scattered directly forward, and is the integrand of the Dirac function. By demanding equality between the second moments of and the that it is approximating, the δ-Eddington approximation ends up transforming layer optical properties as

and

Often the Henyey and Greenstein (1941) function, defined as

Using (4), (5), and (6) in the generalized two-stream solutions (Meador and Weaver 1980) with appropriate coefficients [see King and Harshvardhan’s (1986) Table 2] leads to expressions for layer reflectance and transmittance for both isotropic irradiance and direct-beam incident at solar zenith angle . Linking these expressions vertically for multilayered atmospheres follows well-established methods (e.g., Wiscombe 1977; Coakley et al. 1983; Oreopoulos and Barker 1999).

### b. 3D Monte Carlo model

The 3D Monte Carlo algorithm used here has been described elsewhere (Barker and Davies 1992; Barker et al. 1998, 2003, 2012). It employs cyclic horizontal boundary conditions and can use either simple phase functions, such as or , or tabulated ones obtained from detailed calculation (as discussed in section 2c).

For domains consisting of arrays of columns, photons are injected randomly onto the top of each column, making sure that all columns get sampled. Photon trajectories are determined by generating random distances between scattering events based on local extinction, with updated directions at scattering events based on the phase function of the particle doing the scattering (see Barker et al. 2003). Photons are tracked until they are either absorbed by a particle or the surface or exit the top of the domain (Evans and Marshak 2005). As alluded to in section 1, Monte Carlo solutions of the radiative transfer equation (RTE) can be thought of as infinite-stream models. This is because simulated photons are free to move in any direction.

Let *x* be domain-averaged flux, normalized by top-of-atmosphere (TOA) irradiance, as simulated by a Monte Carlo model using photons. By treating each photon as a Bernoulli trial and using the central limit theorem, it can be shown (Evans and Marshak 2005) that as the number of experiments approaches infinity, the standard deviation of *x* about its expectation value is

This is commonly, and aptly, referred to as Monte Carlo uncertainty.

### c. Cloud optical properties

Clouds were represented fairly simply in this study. Droplet sizes for liquid clouds were assumed to follow gamma distributions with effective variance (Hansen and Travis 1974) and effective radii of either 13 or 15 *μ*m. Their size-integrated optical properties were computed by Wiscombe’s (1980) scattering routines. Spectral integrations were over the CCCma solar wave bands, which are defined by wavelengths *λ*: 0.25–0.69 *μ*m; 0.69*–*1.19 *μ*m; 1.19*–*2.38 *μ*m; and 2.38*–*5.0 *μ*m. The spectral density function was downwelling flux, predicted by Mlawer et al.’s (2000) Code for High-Resolution Accelerated Radiative Transfer (CHARTS), at the midlatitude summer’s tropopause (McClatchey et al. 1972) for . Cumulative phase functions, defined as

where is cosine of scattering angle, were tabulated at 1800 equal increments of . The resulting tables were used in the Monte Carlo model; at each scattering event a uniform random number between 0 and 1 gets generated to represent , which yields from the appropriate table.

The optical properties of ice clouds are known to be important for various applications (Bi et al. 2014; Yang et al. 2015). In this study, the general habit mixture ice cloud model was adopted (Baum et al. 2011, 2014). Spectral optical properties for ice crystals were derived from a database developed by Yang et al. (2013), which covered wide ranges of crystal habits, discrete sizes, wavelengths, and surface roughness conditions. The Amsterdam discrete dipole approximation (ADDA) (Yurkin and Hoekstra 2011) was used for small size parameters, and the improved geometric-optics method (IGOM; Yang and Liou 1996), with inclusion of the edge effect (Bi et al. 2010, 2014), was used for moderate and large size parameters (see Yang et al. 2013). Particle surface roughness was not considered for small size parameters (ADDA calculations); however, moderately and severely rough conditions were considered fully for large size parameters (IGOM simulations).

Yang et al.’s (2013) full dataset was adapted for the CCCma GCM by averaging over *λ*, using the TOA spectrum as a density function, for distributions of habit, size, and roughness (Baum et al. 2011, 2014; Yi et al. 2013). Resulting phase functions and optical properties are functions of effective diameter defined as

where *V*, *A*, and *D* are geometric volume, orientation-averaged projected area, and maximum dimension of ice particle, respectively; *n*(*D*) denotes crystal size distribution; and *f*_{i} indicates the percentage of each ice particle habit and roughness (Yang et al. 2007). For this study, only and 80 *μ*m were used. While Yang et al.’s (2007) tables are for 498 variable angular ranges, their cumulative versions were tabulated, using (9), into 18 000 equal-area bins.

Figure 1 shows size- and spectrally integrated phase functions for the spectral interval 0.25–0.69 *μ*m. All functions exhibit very prominent peaks for with precipitous reductions by . The Lorenz–Mie phase functions show their characteristic minima near , while the ice functions decline almost monotonically back to . Despite the Lorenz–Mie functions having less-pronounced forward peaks than the crystal phase functions, *g* for these distributions of droplets are ~0.867, while for the crystals they are 0.786 and 0.802 for and 80 *μ*m, respectively. Phase functions for this visible wave band were shown since they are responsible for the majority of cloud radiative effects (cf. Räisänen and Barker 2004).

## 3. Results I: Narrowband fluxes for single-layer plane-parallel clouds

Errors in the δ-Eddington approximation that arise from use of minimal information pertaining to cloud particle phase functions can be difficult to explain when looking at broadband fluxes for multilayered inhomogeneous cloud systems. Hence, narrowband results are shown in this section for single homogeneous layers that contain cloud only (see King and Harshvardhan 1986; Li and Ramaswamy 1996). Section 4 addresses the δ-Eddington’s neglect of cloud geometry at the 1-km horizontal scale.

All Monte Carlo results reported in this section used = 400 000 photons per combination of and ; 10 values of from 0.1 to 1.0 in increments of 0.1; and 25 values of in equal logarithmic intervals from 0.1 to 100. Hence, using (8), Monte Carlo uncertainties maximize at for *x* = 0.5, which is ~0.16% uncertainty. Clouds were represented as isolated homogeneous slabs irradiated from above by direct beam only with black underlying surface.

The most basic test of the conventional, analytic δ-Eddington approximation is to compare it to a Monte Carlo model that is equipped with both and optical properties scaled as in (4), (5), and (6). Figures 2a,b and 3a,b show relative differences in single-layer reflectance *r* and transmittance *t* between the δ-Eddington two-stream approximation and a δ-Eddington Monte Carlo model as functions of and for liquid clouds with μm and and ice clouds with μm. Results for μm and μm are very similar to those for μm and μm and so are not shown. Generally speaking, the two-stream aspect of the δ-Eddington is a minor issue, as relative errors for these phase functions are typically less than ~5%. Corresponding results for more absorbing bands are very similar; the notable exception being huge, but unimportant, relative errors for very small values of *t* at ≳ 20 (cf. Fig. 10 from King and Harshvardhan 1986).

Figures 2c,d and 3c,d show the next obvious experiments: comparisons of relative differences in single-layer *r* and *t* between the δ-Eddington two-stream approximation and the Monte Carlo when the latter uses exact phase functions (see Fig. 1) and the former uses just their corresponding *g*. The plot for droplets is reminiscent of those shown by King and Harsvardhan (1986) and Li and Ramaswamy (1996). The δ-Eddington reflects too little and too much at small and large , respectively. Differences, as a function of , are proportionally greatest as ; the single-scattering regime. This can be appreciated by looking at direct-beam backscattered fractions (Wiscombe and Grams 1976) which are defined, using (1), as

Figure 4 shows for droplets with μm and and for crystals with . The value of for the δ-Eddington exceeds the Lorenz–Mie function’s values for , which, as shown in Fig. 2, is where the δ-Eddington reflects too much for small ; the opposite is true for . Results for the ice crystal phase functions are similar: the alternating albedo errors seen in Fig. 3 track closely errors in . As increases out of the single-scattering regime, multiple scattering reduces the importance of differences in , but, as Fig. 4 shows, for a given the sign of the bias remains.

The third band of the CCCma radiation model, which spans *λ* from 1.19 to 2.38 *μ*m, is responsible for most of the atmospheric cloud radiative effects (see Räisänen and Barker 2004). Figure 5 shows relative differences in *r*, *t*, and absorptance *a,* for this band, between δ-Eddington and Monte Carlo estimates. Albedo and transmittance errors retain much of the forms seen in Figs. 2 and 3. Absorptance errors for liquid and ice resemble one another with tendencies to underestimate for most conditions, especially at middling values of and small-to-intermediate . Clearly, the δ-Eddington has a pronounced penchant to underestimate numbers of effective scattering events.

Results shown in Figs. 2, 3, and 5 reveal some of the most basic issues involving two-stream models, but they provide only part of the story pertaining to dynamical modeling, which necessarily centers on errors associated with broadband integrations over myriad cloud configurations and illumination conditions. Hence, the following section is devoted to examining broadband flux and heating rate differences between the multilayer δ-Eddington and its corresponding Monte Carlo model using numerous cloud configurations retrieved from A-Train satellite data.

## 4. Results II: Broadband fluxes and heating rates for A-Train clouds

Results presented in this section expand upon those shown in section 3. Cloud fields used here were represented as 2D fields of particles imbedded in scattering and absorbing gases over reflecting surfaces. Their spatial distributions and some of their properties were inferred from A-Train satellite data, which are discussed in the following subsection. The second subsection describes the setup of surface–atmosphere conditions used in the experiments. In the third subsection, values of domain-average reflected fluxes at TOA, net surface flux, and net atmospheric absorption, as well as heating rate profiles predicted by the multilayer δ-Eddington (in ICA mode) are compared to their Monte Carlo counterparts, which used both detailed phase functions and finite horizontal grid spacing.

### a. A-Train satellite data

Cloud properties used here were derived from data collected between 1 January and 28 January 2007 by *CloudSat*’s 94-GHz cloud-profiling radar (CPR) and *CALIPSO*’s dual-wavelength lidar (Stephens et al. 2002; Winker et al. 2003). Estimates of cloud water content (CWC) came from *CloudSat*’s 2B Radar–Lidar Geometrical Profiling Product (2B-GEOPROF-lidar) data (Mace 2007; Mace et al. 2007). *CloudSat* columns are effectively 1.1 km along track. The radar’s circular footprints are ~1.4 km long, so integrations for ~1 km were oversampled slightly. Each layer is 0.24 km thick. Ground clutter contamination renders the lowest two or three layers unworkable. *CloudSat*’s ECMWF-AUX files report profiles of pressure, temperature *T*, and mixing ratios of water vapor and ozone for each column.

If retrieved values of CWC were reported in the *CloudSat* files, they were used. When retrievals of CWC failed, however, the following values, based on L’Ecuyer et al. (2007), were used:

for cells with *T* > 273 K and

for cells with *T* < 253 K. Mixed-phase conditions were assumed for .

To affect A-Train samples that are meant to represent GCM grid cells, *CloudSat*–*CALIPSO* data were partitioned into 256-km domains. Near-surface precipitation was removed approximately following Barker (2008). Only domains between latitudes 70°S and 70°N with total (i.e., vertically projected) cloud fraction and diurnal mean for the sun-up portion of days were used. This resulted in 32 632, 256-km domains.

Figure 6 shows a summary of the clouds used here. Values of *C* are zonal averages and resemble closely those reported by others (e.g., King et al. 2013). For these domains, mean *C* was ~0.75 which is larger than passive imager-only estimates but close to those from A-Train data (e.g., Mace et al. 2009; Hagihara et al. 2010; Kato et al. 2010). The partition into liquid and ice differs from what gets reported for passive sensors, for if a column had both liquid and ice, both phases received a count regardless of whether they were exposed directly to space or not. Hence, liquid and ice fractions shown in Fig. 6 sum to more than *C*. Note, however, that when liquid and ice fractions are assumed to overlap at random, resulting total fractions are often close to *C*, especially between 40°S and 20°N.

This way of counting clouds has an impact on estimates of cloud mean . For instance, if a passive imager classifies cloud phase at cloud top as ice, the inferred gets associated with ice cloud only, regardless of what might have been below. With active sensor data, one can attempt to draw the distinction and partition between ice and liquid (e.g., Ceccaldi et al. 2013). Hence, values of shown in Fig. 6, at visible wavelengths, have mean values of 24.8 and 7.2 for liquid and ice, respectively. Corresponding values inferred from MODIS imagery are ~15 and ~21 (King et al. 2013). At least the discrepancies between MODIS values and those shown in Fig. 6 resemble differences reported by Stein et al. (2011); their estimates of for ice clouds, using oblate aggregates to represent crystals, are approximately half those reported by MODIS.

The final plot in Fig. 6 is of zonal-average values of

where and are mean and standard deviation of cloud water paths for 256-km domains. This quantity is documented less well than *C* and (e.g., Barker et al. 1996; Rossow et al. 2002; Cole et al. 2011), but what is shown in Fig. 6 seems reasonable: values of ~4 are often associated with overcast stratiform clouds at high latitudes, while values near 1 are expected for multilayer cumuliform clouds that occur often in the tropics (cf. Table 1 from Shonk et al. 2010).

### b. Experimental setup

Two straightforward, and somewhat extreme, cloud-surface scenarios were considered in order to supply likely bounds on estimates of errors associated with two-stream approximations. Scenario A represents oceanic clouds with relatively large particles: all ice clouds had (Okamoto et al. 2010; Stein et al. 2011; King et al. 2013; Deng et al. 2013); liquid cloud droplets had (King et al. 2013); and surface albedos were set at 0.06 for all *λ*. Scenario B approximates land surfaces with smaller cloud particles: all ice clouds had , which is at the small end of current estimates (see previous references); liquid cloud droplets had (King et al. 2013); and surface albedos were 0.06 for and 0.25 for (e.g., Barker et al. 2002). All surfaces were Lambertian. Table 1 summarizes scenarios A and B.

Aerosols were not considered, but Rayleigh scattering was. For the δ-Eddington, Rayleigh scattering was accounted for by simply using *g* = 0. The Monte Carlo model used the standard Rayleigh phase function of

Both models used the CCCma GCM’s description of Rayleigh extinction as a function of *λ* and layer thickness (Li and Barker 2005).

For the δ-Eddington, layer optical properties for clouds and molecular scattering and absorption were merged via extinction weighting to form effective values. For each grid cell of the Monte Carlo, a cumulative extinction function was defined as a single increment for each type of attenuator. When, in the simulation, a photon ensemble was deemed to undergo an interaction with matter, a uniform random number between 0 and 1 was generated and an attenuating species selected using the cumulative extinction function (see Barker et al. 2003).

The A-Train samples clouds over a very narrow range of time of day. To get a picture of the diurnal range and mean values of potential errors due to the use of the δ-Eddington approximation, for each 256-km A-Train domain was sampled by randomly selecting a time between sunrise and noon, depending on day number and latitude (Barker et al. 1998). In so doing, it is recognized that cloud conditions particular to the A-Train’s crossing times were extrapolated to entire sun-up periods and that this will inflict biases on some regions. It has to be stressed here, however, that when it comes to setting global distributions of cloud vertical structure, the A-Train’s observing system is all that is available. Moreover, the point of this study was to establish approximate errors, not necessarily describe them in full detail.

To summarize, for each 256-km A-Train domain, for both scenarios A and B, domain-average radiative transfer simulations were computed for the following:

Multilayer, analytic δ-Eddington in ICA mode,

1D Monte Carlo using with (reported in the appendix only),

1D Monte Carlo using exact phase functions with , and

3D Monte Carlo using exact phase functions with km.

Each Monte Carlo simulation used *N*_{p} = 100 000 photons. When , photons were injected into the domains parallel to the along-track direction of the A-Train satellites.

### c. Broadband fluxes and heating rates

This section focuses on errors in broadband fluxes and heating rates incurred by the analytic δ-Eddington’s neglect of phase function details and cloud geometry. Figures 7a–c and 8a–c show -dependent distributions of TOA reflected flux, net surface flux, and net atmospheric absorption differences between the analytic δ-Eddington and the Monte Carlo using exact phase functions with for the cloudy portions of 256-km A-Train domains. The characters of the mean and median curves could have been anticipated from the contour plots in Figs. 2, 3, and 5. Basically, the δ-Eddington’s underestimation of reflectance at small stems from neglect of details in both ice and liquid phase functions. At large , however, its overestimation of reflectance is fueled exclusively by liquid clouds, with ice clouds offering some temperance and indeed being responsible for the ~20% of cases at in which the δ-Eddington underestimates reflectance. Similarly, the δ-Eddington’s near-systematic underestimation of atmospheric absorption is due to both liquid and ice clouds (cf. large swaths of negative values across Figs. 5c and 5f).

The center rows of Figs. 7 and 8 show corresponding differences between the Monte Carlo using exact phase functions with (i.e., ICA mode) and km. They show the familiar effects of neglecting net horizontal transport of radiation [see chapters in Marshak and Davis (2005a)]. The largest errors lead from the 1D solution reflecting too much (transmitting too little) radiation at large . Mean-bias errors for are ~7 W m^{−2}, which, when scaled by cloud fractions, become ~5 W m^{−2} for full domain averages. This is close to values reported by Cole et al. (2005) for multiscale modeling clouds, Pincus et al. (2005) for simulated cumulus clouds, and Barker et al. (1999) for simulated tropical convective clouds.

While there are several cases in which the solution with overreflects by >15 W m^{−2} for , there are some cases where it actually reflects less. However, median difference for the latter cases is only roughly −1 W m^{−2}, and typical Monte Carlo uncertainty for flux differences is approximately W m^{−2}, where is defined in (8) and W m^{−2} (Kopp and Lean 2011; Ball et al. 2012). Hence, for many of these cases, it is likely that the 1D solution actually reflected slightly more than the 3D solution.

To get a better sense of the magnitude of Monte Carlo noise present in the results shown in Figs. 7 and 8, let total variance for flux differences in a bin be

where and are variances that arise from atmospheric variations and Monte Carlo noise, respectively. From (8), Monte Carlo variance for the cloudy portion of a domain is approximately

where 〈⋅〉 denotes bin mean. When differences are between a Monte Carlo simulation and the noiseless δ-Eddington, (17) represents all the Monte Carlo noise. For differences between two Monte Carlo simulations, the last term in (16) becomes ~. For estimates of atmospheric absorption, are enhanced further since they were computed as differences between net fluxes at TOA and surface.

Figure 9 shows mean and median differences between Monte Carlo results for scenario A, with and km (a replot of curves in Fig. 7) as well as and . This demonstrates that, for a given range of , the lion’s share of the dispersions indicated in Figs. 7 and 8 stem from variations in cloud conditions, as opposed to Monte Carlo noise. Thus, only small reductions to the quantiles shown in Figs. 7 and 8 for reflected flux at TOA and net flux at surface are needed to approximately remove Monte Carlo noise. Adjustments for atmospheric absorption are largest but are still fairly minor for most .

At small there is a slight tendency for the km solution to reflect most. This is the outcome of the well-documented cloud side-illumination effect. While its magnitude is less than that for idealized clouds (cf. Welch and Wielicki 1985; Barker and Davies 1992; Barker et al. 2003), it resembles closely values reported by Cole et al. (2005). This is because, when integrated over wide ranges of realistic clouds, side illumination of dense low-level clouds is much attenuated by frequently occurring layers of thin upper-level clouds, which increase effective at low altitudes.

As for atmospheric absorption, when km, the atmosphere absorbs ~1.5–2 W m^{−2} more than the ICA for a wide range of . As discussed by Barker et al. (1998) and Barker (2005), this is often because of both illumination of cloud sides at and trapping photons within clouds at larger . As shown later, there is minimal enhancement of clear-sky absorption in the lower layers as a result of radiation emerging from the sides of finite clouds in the downward direction (e.g., Welch and Wielicki 1985).

From Figs. 7a–f and 8a–f, it is obvious that neglecting phase function details by the δ-Eddington and cloud geometry by 1D models generally affect -dependent biases of almost the same sign. Hence, these two approximations work together to enhance overall bias errors, as shown in Figs. 7g–i and 8g–i. Errors in TOA reflectance for are ~12 W m^{−2} for the cloudy portion of domains—or ~9 W m^{−2} for full domain averages—with 90% of the errors between about −5 and 25 W m^{−2}. At low sun, mean-bias errors for TOA reflectance are about −4 W m^{−2}, but ~5% of them are larger than −10 W m^{−2}. Corresponding values for net surface flux are slightly larger but of opposite sign. In about 5% of cases it may be that the δ-Eddington underestimates net surface flux by >25 W m^{−2}. Moreover, for all but the largest the δ-Eddington underestimates atmospheric absorption for cloudy conditions by 2–3 W m^{−2}.

The top grouping of plots in Figs. 10 and 11 show the same data as were shown in Figs. 7 and 8 but sorted and averaged as a function of latitude and . As expected, differences between the δ-Eddington and Monte Carlo using show that the strong -dependent biases seen in Figs. 7 and 8 are quite independent of latitude despite integrations over distributions of occurrences of cloud types. The same can be said, for the most part, for differences between the Monte Carlo models using and km and, hence, for the combined errors. The obvious anomaly is for TOA reflectance and net surface flux for , where there are distinct local maxima near the descending branches of the Hadley cells, centered at ~30°S and ~20°N. This is because of the relative lack of high ice clouds (see Fig. 6) and the strong effects of using improper liquid phase functions and side illumination of broken liquid clouds.

The lower set of plots in Fig. 10 shows standard deviation of flux differences that correspond to the mean-bias errors shown above them. These values were determined by local variability of clouds convolved with the strengths of errors associated with neglect of phase function details and horizontal transport. For instance, areas poleward of ~50° exhibit marked reductions in standard deviations for all relative to the rest of Earth. This is because of the persistence of very cloudy conditions. At ~30°S and ~20°N, standard deviations are maximal as a result of relatively changeable low-level clouds being exposed to direct beam during times of no, or only very thin, ice clouds aloft.

Figures 10 and 11 also show that, in the descending branches of the Hadley cells, near 30°S and 20°N, underestimation of atmospheric absorption by the δ-Eddington is greatest at intermediate values of . Figures 12a and 12b show differences in mean all-sky heating rates for all cloud conditions for scenario A. It is clear from this plot, though not surprising by now, that the weak absorption due to the δ-Eddington’s neglect of phase function details extends throughout the entire lower atmosphere. It also shows that, for low-level liquid clouds, neglect of horizontal transfer, and hence side illumination of low liquid clouds, underestimates heating to the same extent as neglect of phase function details. Upon screening out ice clouds, Figs. 12c and 12d show that effects due to neglect of both phase function details and horizontal transfer enhance underestimation of heating by low clouds; almost by −0.1 K day^{−1}, which is, however, only about a 5% error.

Another area with notable differences in atmospheric absorption is from 40°S to 10°N at : the only place and time that the δ-Eddington overestimates total atmospheric absorption. Figures 13a and 13b show mean all-sky heating rate profiles for all cloud conditions for scenario B. Through much of the troposphere, errors due to neglect of horizontal transfer are very weak. For thin ice clouds above ~12 km, the δ-Eddington predicts too little heating (cf. Fig. 5), but below that it predicts too much; the δ-Eddington allows too much radiation through thin ice clouds where it is then overly absorbed by liquid clouds and water vapor. Figures 13c and 13d show that for cases with no ice cloud, the δ-Eddington now exhibits too little heating near the tops of low liquid clouds, as expected from Fig. 5, but again allows too much transmittance, so heating by near-surface water vapor is overestimated.

While the heating rate bias errors just discussed are clear and explainable, they are nevertheless usually smaller than 3% of the local heating rate. The left plots in Fig. 14 show mean maximum heating rates for layers below 20 km for the δ-Eddington and the Monte Carlo models with and km as functions of and latitude for scenario A. The right plots show mean maximum heating rate differences between the models as listed on the plots. Note that, for any individual domain, the layer with maximum heating rate is not necessarily the layer with maximum difference in heating rate between models.

Figure 14 simply tries to put the largest bias errors into context with the estimated heating rates. The largest mean maximum heating rate differences occur between the δ-Eddington and the Monte Carlo model with in areas of maximum cloudiness (cf. Fig. 6) at intermediate values of . This can be explained by results shown in Fig. 5. The largest mean maximum differences are about −0.5 K day^{−1}, or about 10%–15% of corresponding mean maximum heating rates. Maximum differences due to cloud geometry are only about −0.1 K day^{−1} and occur, as expected, at small . While they are slightly larger for scenario B (not shown), they are again up to 10%–15% of corresponding mean maximum heating rates.

## 5. Summary and conclusions

This study addressed two basic omissions in multilayer two-stream approximations of the solar radiative transfer equation when applied to cloudy atmospheres: namely, their neglect of details in scattering phase functions for cloud particles and horizontal transfer of radiation. A-Train satellite data were partitioned into a large number of 256-km segments so as to approximately represent typical GCM cells (Astin and Di Girolamo 2006). These domains were used to show that biases affected by the two mentioned omissions generally have the same sign as a function of solar zenith angle and that they are approximately additive, thus leading to overall bias errors that could be important when attempting to estimate climate sensitivities due to changes in Earth–atmosphere conditions.

Using single-layer isolated clouds, it was demonstrated that phase functions for typical cloud particle size distributions have sufficient detail that, when only their asymmetry parameters are used by two-stream models (the δ-Eddington in particular), systematic errors in reflectance, transmittance, and absorptance will follow. This general result is not new (cf. King and Harshvardhan 1986), but the provision of broadband estimates integrated over most cloud and illumination conditions is. For high-sun conditions, one can expect the δ-Eddington to underestimate surface flux for the cloudy portion of a GCM’s grid cell by ~10 W m^{−2} because of neglect of phase function details. Similarly, for most , the δ-Eddington underestimates total atmospheric absorption by about 2 W m^{−2} for the cloudy portions of GCM grid cells. Corresponding bias errors due to neglect of horizontal transport, as provided here, are approximately of comparable magnitude yet are underestimates, because i) 2D rather than 3D clouds were used (cf. Barker 1996; Pincus et al. 2005), and ii) horizontal grid spacings were only 1 km. Regarding the last point, however, Cole et al.’s (2005) simulations suggest that, for most cloud types, flux differences for km are already very close to their asymptotic values. Therefore, taken together, the δ-Eddington can be expected to underestimate net surface flux for the cloudy portion of a GCM’s grid cell by ~15 W m^{−2} at large , tapering to zero after reversing sign, as (approximately the opposite of that for reflected fluxes).

Layer heating rate errors, once averaged over all cloud conditions, are largest because of the omission of phase function details. Corresponding errors due to neglect of horizontal transfer are smaller. This is not entirely unexpected (see Cole et al. 2005). The result is that maximum heating rate errors in the lower atmosphere due to neglect of phase function details and horizontal transfer of radiation can be expected to be up to about 10%–15% of the corresponding maximum heating rates. Most errors, however, are much smaller at about ±2%.

It is recognized fully that, while the dataset used and many assumptions made in this study are not perfect, they are probably not catastrophically far from both the best available and reality. Hence, estimates reported here of flux errors associated with the multilayer δ-Eddington model, and two streams in general, as a result of neglect of details pertaining to cloud particle phase functions and horizontal transport of solar radiation are likely to be *quantitatively approximate* but *qualitatively accurate*.

The 256-km domains used in this study were supposed to represent grid cells of conventional GCMs, yet unlike GCMs, they were resolved to about the 1-km scale. Hence, this study sidestepped an additional source of flux error endemic to GCM radiation calculations: generation of unresolved clouds (Pincus et al. 2003). Whether these errors are correlated with those analyzed here remains to be seen.

To address the impacts that biases discussed here would have on GCM simulations, one would have to employ a solar transfer model, like the Monte Carlo used in this study, which is able to utilize important details of cloud particle phase functions and horizontal transport. This could be realized using a so-called superparameterized GCM (e.g., Cole et al. 2005; Randall 2013), but it would be only partially fruitful in conventional GCMs that are, as yet, unable to provide suitable descriptions of clouds at unresolved scales.

## Acknowledgments

This study was supported by a grant from the U.S. Department of Energy Atmospheric Radiation Measurement Program (subcontracted to Environment Canada from The Pennsylvania State University), and a contract from the European Space Agency (subcontracted to Environment Canada from the Free University of Berlin).

### APPENDIX

#### Biases due to the Use of the Henyey–Greenstein Phase Functions

The Henyey–Greenstein function , as defined in (7) and shown in Fig. 1, is used often to approximate exact phase functions, especially when fluxes are the only concern. Use of to simplify a Monte Carlo model does, however, entail biases relative to use of exact phase functions. This short appendix shows, by way of example, that use of may not always be a good choice.

Using the same input data that were used to produce Figs. 7 and 8, Fig. A1 shows differences between Monte Carlo models using and the exact phase functions for scenario A with . For midrange values of , use of leads to TOA reflectances that are ~5 W m^{−2} too large on average. As increases, mean differences are reduced largely because of biases associated with ice crystals with large diffraction peaks (cf. Fig. 1). Differences in net fluxes at the surface steadily increase in magnitude with increasing , maximizing near −5 W m^{−2}. As for atmospheric absorption, use of yields minor underestimates for small and typically ~3 W m^{−2} overestimation at = 1. As such, since errors due to use of , rather than their exact counterparts, are approximately comparable in magnitude to those due to use of the δ-Eddington, if the δ-Eddington ends up being deemed unacceptable for reasons presented in the main text, it stands to reason that use of will suffer the same criticism, too.

## REFERENCES

*3D Cloud Structure and Radiative Transfer*, A. Marshak and A. Davis, Eds., Springer-Verlag, 449–486, doi:.

**33**, 657–676, doi:.

*Three-Dimensional Cloud Structure and Radiative Transfer*. A. Marshak and A. Davis, Eds., Springer-Verlag, 243–281.

*3D Radiative Transfer in Cloudy Atmospheres*. Springer, 686 pp., doi:.

*3D Radiative Transfer in Cloudy Atmospheres*, A. Marshak and A. B. Davis, Eds., Springer, 653–686.

*Lidar Remote Sensing for Industry and Environment Monitoring III*, U. N. Singh, T. Itabe, and Z Liu, Eds., International Society for Optical Engineering (SPIE Proceedings, Vol. 4893),

## Footnotes

^{1}

For readability, the term “flux density” is replaced, hereinafter, by simply “flux,” noting that flux implies units of watts per meter squared.