## Abstract

A theoretical study has been conducted on the effects of cloud horizontal inhomogeneity on cloud albedo bias. A two-dimensional (2D) version of the Spherical Harmonic Discrete Ordinate Method (SHDOM) is used to estimate the albedo bias of the plane-parallel (PP–IPA) and independent pixel (IPA–2D) approximations for a wide range of 2D cloud fields obtained from Landsat. They include single-layer trade cumulus, open and closed cell broken stratocumulus, and solid stratocumulus boundary layer cloud fields over ocean. Findings are presented on a variety of averaging scales and are summarized as a function of cloud fraction, mean cloud optical depth, cloud aspect ratio, standard deviation of optical depth, and the gamma function parameter *ν* (a measure of the width of the optical depth distribution). Biases are found to be small for small cloud fraction or mean optical depth, where the cloud fields under study behave linearly. They are large (up to 0.20 for PP–IPA bias, −0.12 for IPA–2D bias) for large *ν*. On a scene-average basis, PP–IPA bias can reach 0.30, while IPA–2D bias reaches its largest magnitude at −0.07. Biases due to horizontal transport (IPA–2D) are much smaller than PP–IPA biases but account for 20% rms of the bias overall.

Limitations of this work include the particular cloud field set used, assumptions of conservative scattering, constant cloud droplet size, no gas absorption or surface reflectance, and restriction to 2D radiative transport. The Landsat data used may also be affected by radiative smoothing.

## 1. Introduction

Recent results by Cahalan et al. (1994a; Cahalan 1994b) have shown that the albedo bias of the independent pixel approximation (IPA-3D) for marine boundary layer clouds is about 1% of the typical albedo, while the plane-parallel (PP-IPA) bias is on the order of 15%. That study included only overcast conditions at one location off the coast of California. This note extends Cahalan’s work to a wide range of marine boundary layer cloud conditions, including broken stratus and trade cumulus cases, found in 45 Landsat scenes of clouds over ocean, ranging from about 20°S to 40°N latitude and from 15° to 150°W longitude. Computations of the radiance field of these scenes using both IPA and 2D assumptions are performed at 0.83 *μ*m. These calculations allow estimation of the error in the PP and IP approximations for a wide range of cloud parameters, complementing Cahalan’s results for overcast marine stratocumulus.

A limitation of the current work is that the IPA is used initially to derive the Landsat optical depth distributions that form the reference field for later testing of the IPA and 2D calculations. Furthermore, since only 2D cloud fields are considered here (effectively, cloud streets), the full effect of realistic cloud fields has not been captured. Horizontal transport would be more important in a 3D field, leading to somewhat higher biases (Chambers 1997). Finally, only one wavelength, with no absorption, has been considered. Nevertheless, the results provide some insight into the effect of broken clouds on albedo bias.

Section 2 briefly describes the two-dimensional cloud fields inferred from Landsat. Section 3 describes the numerical radiation model as well as the radiative quantities studied in this work. Section 4 compares the cloud albedo biases based on the IPA and 2D results. Section 5 presents the major conclusions from this work.

## 2. Landsat-inferred horizontal inhomogeneity

Forty-five Landsat scenes of cloud fields over ocean were available for this study. The criteria that governed the selection of these particular scenes, and their characteristics, are described in detail in Harshvardhan et al. (1994) and Chambers et al. (1996). Landsat provides both better resolution and far more data volume than ground-based measurements of cloud properties. While the Landsat data may suffer from radiative smoothing effects at small scales (e.g., Marshak et al. 1995), test cases have shown that this does not affect results at the averaging scales used in this paper (Chambers et al. 1996, hereafter CWE). Note that seven of the scenes from Harshvardhan et al. (1994) were rejected because they had more than 10% of the pixels saturated in radiance. The final set of 45 scenes contains a majority that are subtropical oceanic boundary layer cloud and therefore should not be interpreted in any way as a globally representative set.

Based on the Landsat reflectance data, interpreted using IPA, the cloud optical depth *τ* at 0.83 *μ*m has been derived in each pixel as described in Harshvardhan et al. (1994). Using an empirical relation from Minnis et al. (1992), an estimate of the physical thickness Δ*z* of the cloud in each pixel can be made. The relation for Δ*z* was derived from simultaneous satellite and surface observations at a 30-km scale, but there is some evidence (CWE) that it is reasonable at the Landsat pixel scale as well. The extinction coefficient *β* can then be calculated for each pixel and distributed through the cloud depth. For this study, the cloud is assumed to have uniform extinction with height in each pixel, that is, *β* = *τ*/Δ*z.* This extinction field is then discretized onto the computational grid, while maintaining total column optical depth at the Landsat value (see CWE). The cloud-top height is held constant, while the cloud base varies according to the calculated Δ*z.* For inversion-capped boundary layer clouds over ocean, this is considered more physically realistic than a constant base height assumption.

Conservative scattering (*ω*_{0} = 1) has been assumed throughout this study. A Mie phase function for a gamma distribution of water droplets with an effective radius, *r*_{e}, of 10 *μ*m and an effective variance of 0.1 is used and is assumed constant throughout the cloud. The mean cosine of the scattering phase function is 0.856, which is close to the value used by Cahalan et al. (1994b) for the Henyey–Greenstein phase function. The effect of atmospheric gases on the radiation within the cloud is neglected. While this may affect the initial retrieval of optical depth (Oreopoulos and Davies 1997, manuscript submitted to *J. Climate*), it has little effect on horizontal radiative transfer characteristics and thus the conclusions of this paper (Chambers 1997).

From each Landsat scene, a series of 20 randomly placed horizontal scan-line samples is taken. The orientation of clouds in the scenes is random enough that horizontal sampling is sufficient. Each scan-line sample is 200 pixels (5.7 km) long. Sensitivity studies confirm this is a sufficient sample length for thin boundary layer clouds. The sampling is required because computing 2D transport in the entire scene is currently beyond the bounds of computer resources in terms of both time and memory. However, PP–IPA bias can be calculated easily for the full scene and compared to the bias found from the sampled scans to test convergence. The correlation between the two is Δ*α*_{sample} = 0.98 × Δ*α*_{full} + 0.004, with a regression coefficient of 0.95 and a standard error of estimate of 0.01, indicating excellent convergence. What differences exist are traceable to a difference in mean properties between the sample and the full scene and are accounted for since all results herein are examined in terms of cloud properties. A cyclical horizontal boundary condition is imposed in the 2D solution, with three logarithmic interpolating points added on the right-hand side of most scans to blend the cloud field from edge to edge. In a few cases, additional interpolating points are added to keep the edge gradient no greater than the largest interior gradient of the scene.

For the scan-line samples taken from the Landsat scenes, mean cloud optical depth ranges from 0.3 to 73.3. Cloud fraction *A*_{c} ranges from 0.0 to 1.0. Aspect ratio (defined in section 3c) ranges from 0.0 to 2.73. Standard deviation of optical depth is between 0.08 and 45. The gamma function parameter (defined in section 3d), measuring the width of the optical depth distribution, ranges from 0.04 to 310.

## 3. Radiative model

### a. SHDOM

The spherical harmonic discrete ordinate method (SHDOM) of Evans is used here in both an independent pixel and two-dimensional solution mode. More details on this method are given in CWE. The algorithm is an earlier version of a 3D method described in Evans (1997, manuscript submitted to *J. Atmos. Sci.*). In brief, it uses both spherical harmonics and discrete ordinates to represent the radiance field during different parts of the solution algorithm. The spherical harmonics are employed for efficiently computing the source function, including the scattering integral. The discrete ordinates are used to integrate the radiative transfer equation through the spatial grid. The solution method is to simply iterate between the source function and radiance field, akin to a successive order of scattering approach.

The two-dimensional (*x* and *z*) cloud fields generated from the Landsat scenes are provided as input to the SHDOM code. A solution is obtained with *L* = 11 zenith terms and *M* = 11 azimuthal terms for the spherical harmonic truncation, and with *N*_{μ} = 12 zenith angles and *N*_{ϕ} = 24 azimuth angles in the discrete ordinate discretization. Using the *δ*-M method (Wiscombe 1977) for the phase function, this truncation has been found to provide good accuracy (mean error about 0.5%) for fluxes, upon which the results of this paper depend. A Monte Carlo method has been used to verify the results for two typical broken cloud scan-line samples, showing agreement within the uncertainty of the Monte Carlo solution using 20 million photons. Along with prior work by one author (Evans), this is considered sufficient verification of the SHDOM numerical model.

For each scan-line sample of each scene, solutions are obtained for solar zenith angles, *θ*_{0}, of 0°, 49°, and 63° (*μ*_{0} = cos *θ*_{0} = 1, 0.65, 0.45) and with both the IPA and 2D solution options.

### b. Albedo and albedo bias

The spectral albedo *α* (*x, **θ*_{0}) is given by

where *F*_{⊚} is the solar flux constant at the appropriate wavelength and *F*^{↑} is the upward spectral flux. All results in this paper are restricted to monochromatic calculations at 0.83 *μ*m. Here, *F*^{↑} is computed by SHDOM as a function of *x* for the 2D scene using both a 2D method and the IPA. To obtain the scan-line sample average IPA and 2D albedos, then, *α*(*x*) is averaged over cloudy pixels in the 200 pixel sample. The PP albedo is retrieved from the average cloud optical depth of the scan (i.e., set *τ*_{c}(*x*) ≡ *τ̄*_{c} for all *x* corresponding to cloudy pixels), using the results for albedo versus optical depth as shown in Fig. 1. Averages at larger scales are obtained analogously.

The IPA–2D albedo bias for cloudy pixels only is defined as

while the PP;t2– IPA bias for cloud areas is

where Φ_{p} = 0 for *τ* = 0, and Φ_{p} = 1 for *τ* > 0. The PP–2D bias is therefore the sum of (2) and (3). Albedo bias for the whole field, applicable to some general circulation models (GCM), is found by multiplying either definition by the cloud fraction. To obtain the relative bias, (2) is divided by *ᾱ*_{IPA}, (3) by *α*_{PP}. These definitions are consistent with the work of Cahalan et al. (1994a, 1994b). As in that work, relative bias is given here in percent, while absolute bias is not.

### c. Aspect ratio

In order to classify the cloud fields in this study, a scan average cloud aspect ratio is computed according to CWE as

where *ρ*(Δ*z*) is the autocorrelation coefficient of the cloud thickness (i.e., the degree to which Δ*z*(*x*) and Δ*z*(*x* + *dx*) are correlated). This definition measures the ratio of the largest vertical cloud extent (Δ*z*_{max}) to a horizontal width that approximates the largest feature width at half-maximum for the sample (see CWE for more details and an example). Scan-line samples that are completely clear over the 5.7-km sample are assigned an aspect ratio of zero. The aspect ratio for completely overcast scan-line samples is computed from the largest vertical cloud extent divided by the full 5.7-km sample width. Other definitions are certainly possible, but this one has been found to be useful.

### d. Gamma function parameter ν

The gamma function has been proposed by Barker et al. (1996) as an analytic approximation for the distribution of optical depths in an actual cloud field. Such an analytic form can enable an approximate IPA solution to be computed for little more than the cost of a PP result. The shape of the distribution is controlled by the parameter *ν*, which is given by

where *τ̄*_{c} is the mean and *σ*_{τ} is the standard deviation of the distribution of optical depths. Small values of *ν* produce optical depths distributed in a monotonically decreasing fashion, characteristic of broken cloud fields. Values of *ν* larger than one produce distributions with a more Gaussian shape and a well-defined mean, characteristic of overcast clouds. This parameter can also be viewed as a measure of the relative variability of a scene (*ν* < 1 ⇒ *σ*_{τ} > *τ̄*_{c} ⇒ highly variable).

## 4. Bias of PP and IPA cloud albedos

At the pixel level, the variability between IPA and 2D albedo is large, due to horizontal transport displacing the peaks and valleys in the radiative flux and thus masking albedo biases. The results are therefore first averaged over the cloudy region of a scan-line sample. The albedo bias over such small horizontal scales is still large and uncorrelated to cloud field properties (as found by Cahalan et al. 1994b), due again to horizontal transport. The behavior of the bias becomes reasonable only on the scale of a Landsat scene (58 km by 58 km) or even larger. The mean cloud albedos for the 20 scan-line samples in each Landsat scene are therefore averaged together in several ways to provide a mesoscale or GCM-scale view. The results are first presented on a scene average basis, in which results for the 20 scan-line samples in each Landsat scene are averaged. The second averaging method classifies each scan-line sample in terms of a cloud parameter, then averages all the samples from any Landsat scene in that cloud parameter bin. Results using the first averaging method still show considerable scatter when plotted against any cloud parameter, due to within-scene variability of the cloud field and sensitivity to variations in the other parameters within the cloud scene. The second approach allows the variability to be reduced using conditional sampling as a function of cloud parameters (cloud fraction *A*_{c}, optical depth *τ*_{c}, aspect ratio, standard deviation *σ*_{τ}, or gamma parameter *ν*) and provides a clearer picture of the trend of albedo bias as a function of these cloud physical properties. The second approach has the same statistical distribution as the scene-average results but provides better sampling for uncertainty estimation. Finally, averaging is done in three broad cloud classes to generate mesoscale and climatological-scale statistics.

### a. Scene-average results

Albedo biases averaged for the 20 scan-line samples in each scene represent the bias at a (58 km)^{2} scale, applicable to GCM modeling. Results are presented as a function of the standard deviation of optical depth, *σ*_{τ}, of the scene. This is the single parameter that best correlates with the PP–IPA bias (aspect ratio best correlates with the IPA–2D bias). In terms of the other cloud physical parameters, the variation is similar to that shown in section 4b, generally falling within the uncertainty (1-sigma) of those results. As shown in Fig. 2a, the PP–IPA bias generally increases with *σ*_{τ}, that is, with the variability of the optical depth field, while the IPA–2D bias is an order of magnitude or more smaller. The PP–IPA bias is strictly positive, due to the shape of the curve in Fig. 1, while for overhead sun the IPA–2D bias is almost evenly distributed among positive and negative values: positive (IPA > 2D) for the broken cloud scenes in which leakage of photons out cloud sides and to the surface dominates, negative (IPA < 2D) for more overcast ones in which “filling in” (i.e., reflected radiance from cloud sides that emerges at the top of a cloud-free area) between closely spaced clouds or cloud cells dominates.

At the higher sun angles shown in Figs. 2b and 2c the results are similar, with somewhat more scatter, particularly in the IPA–2D bias, as *σ*_{τ} increases. The PP–IPA bias remains strictly positive, while the IPA–2D bias now tends to more negative values, particularly for the more variable scenes. This is due to the fact that at slant sun the cloud sides intercept and reflect some of the radiation that impinges on cloud-free regions at the top of the field, making the 2D albedo larger than the IPA albedo (filling-in effect). At *θ*_{0} = 63°, this effect is overwhelmed for some of the more uniform scenes by shadowing of cloud sides in 2D, resulting in a larger IPA albedo and a positive bias.

The importance of horizontal radiative transport to albedo bias can be estimated from the size of the IPA–2D bias relative to the PP–2D bias. For the scene averages, the IPA–2D bias accounts for 20% root-mean-square (rms) of the albedo bias, indicating horizontal transport is generally a small, but not necessarily negligible, factor in the albedo bias.

### b. Parameter space results

Albedo biases averaged over the 5.7-km scan-line samples are binned and averaged to examine the sensitivity of the bias to cloud properties. Results are shown in Figs. 3a–e in terms of five cloud parameters: cloud fraction, mean cloud optical depth, cloud aspect ratio, standard deviation of optical depth, and gamma function parameter *ν*. The cloud fraction bins are those planned for use in the Clouds and the Earth’s Radiant Energy System (CERES) experiment (Baum et al. 1995), while bins for the other cloud properties were chosen after examining the data. The trends shown are not very sensitive to bin boundaries. Note that the PP–IPA bias can theoretically depend only on the frequency distribution of optical depth (i.e., *τ̄*, *σ*_{τ}, and *ν*). The apparent dependence on cloud fraction and aspect ratio comes about because *A*_{c} and aspect ratio are correlated with *τ̄*, *σ*_{τ}, or *ν* in this set of cloud scenes. (For example, small *A*_{c} tends to mean thin clouds in this dataset.)

Two sets of error bars are given in these figures for the 63° solar zenith angle results where the error is generally largest. The error bar through the data points is an estimate of the error in the mean value. It is generated by splitting the cloud scenes into two completely uncorrelated sets (scenes from different years) and taking the difference in the mean value found for each set. This error is quite small, except in bins where there are very few (about 20 or fewer) samples, suggesting good confidence in the mean bias value. The second error bar, slightly offset from the data points, represents the cloud fraction weighted 1-sigma variability of the scan-line samples in that bin.

The variation with cloud fraction is shown in Fig. 3a. The PP–IPA bias is near zero for small cloud fractions, which in these Landsat scenes correspond to thin clouds at the linear limit of the retrieval curve in Fig. 1. The bias increases for intermediate cloud fractions, then falls again in the overcast limit. This behavior is due to the correlation with cloud field variability in this dataset. The parameter *ν* is generally larger, implying less variability, when the cloud fraction is 1.0, compared to the broken cloud fields in this set. The IPA–2D bias has a similar trend, but with a much smaller magnitude and the opposite sign. The trends are similar at all three solar zenith angles, though the bias tends to be smaller at overhead sun (*θ*_{0} = 0°). At slant sun, the IPA–2D bias again is consistent, with added energy incident into the sides of broken clouds leading to a larger 2D albedo.

In terms of optical depth, Fig. 3b, the PP–IPA bias starts out very small as it should in the linear limit. If the optical depth ranged high enough, this bias should drop again at large *τ*_{c} as the retrieval curve (Fig. 1) becomes nearly linear again around *τ* = 50 (unless the variability also increased with *τ*_{c}). Since in this sample the maximum optical depth is below 30, the PP–IPA bias remains quite large. The IPA–2D bias, on the other hand, peaks below *τ* = 5 before dropping again, most obviously at *θ*_{0} = 63°.

The results as a function of aspect ratio are shown in Fig. 3c. The PP–IPA bias increases almost linearly with aspect ratio for all solar zenith angles. The IPA–2D bias has an even more linear trend. Aspect ratio is, in fact, the parameter that best correlates the IPA–2D scan-line sample bias (*r* = 0.56 to 0.79 as *θ*_{0} increases), suggesting that cloud geometry is somewhat more important than the features of the optical depth distribution in controlling horizontal radiative transport.

In terms of the standard deviation of optical depth *σ*_{τ}, results are shown in Fig. 3d. The PP–IPA bias in this case also increases quite linearly with *σ*_{τ}, at all sun angles. This is the parameter that best correlates the PP–IPA scan-line sample bias, with a regression coefficient between 0.71 and 0.88 depending on sun angle. The IPA–2D bias shows a similar behavior but with an increasingly negative bias as the sun drops from overhead for the more variable scenes due to energy intercepted by cloud sides.

Finally, the trends in terms of gamma function parameter *ν* can be observed from Fig. 3e. The biases are nearly constant for large *ν*. As *ν* decreases below one, implying a highly variable cloud field, the biases increase rapidly. In fact, based on the magnitude of the largest biases in this figure, *ν* does the best job of grouping scans with large Δ*α*. This suggests that *ν* is a powerful parameter in determining the influence of horizontal cloud inhomogeneity and should be considered in cloud field classification studies.

Overall, as found by Cahalan et al. (1994a), the albedo bias is less sensitive to mean properties (*A*_{c}, *τ*) than it is to properties of the distribution of *τ* (*σ*_{τ}, *ν*, and aspect ratio).

All the parameter space results show the IPA–2D bias to be smaller, sometimes by orders of magnitude, than the PP–IPA bias, except at small optical depth or cloud fraction where the cloud field approaches the linear limit and at large aspect ratio or small *ν* where the IPA–2D bias can be half the PP–IPA bias at slant sun.

### c. Cloud amount class results

To examine larger-scale results, the Landsat scenes are now composited together in several ways. Results for the whole dataset are summarized in Tables 1a,and 1b, in terms of three cloud amount classes roughly representing scattered, broken, and overcast cloud fields.

Table 1a summarizes the mean and 1-sigma variability of the scene average [(58 km)^{2}] properties for all scenes in a cloud class. [Similar results for spherical albedo bias in Barker et al. (1996) use the same initial dataset but slightly different categorization and scene rejection criteria.] Such results are most applicable to a GCM or mesoscale model. The albedo bias and the variability in that bias are found to be largest for the smallest cloud fractions. This is consistent with the amount of variability in each of the cloud classes, as measured by *σ*_{τ} and *ν*. Cahalan et al. (1994a), in measurements during the First ISCCP (International Satellite Cloud Climatology Project) Regional Experiment off the coast of California, found the largest variability in overcast clouds. In measurements taken during the Atlantic Stratocumulus Experiment, on the other hand, Cahalan et al. (1995) found more variability in broken clouds. More data, obtained from a variety of instruments, may be required to better understand cloud variability on a global basis.

Much of the work of Cahalan et al. (1994a, 1994b) was done using data collected over a period of several days or weeks. To better compare with long-term average data of this sort, Table 1b groups the Landsat scenes in each cloud amount class into composite “scattered,” “broken,” and “overcast” cloud scenes and presents the albedo bias for each [akin to (225 km)^{2} area averages]. It also gives results for a single cloud scene compositing all 45 Landsat scenes [“all scenes,” (400 km)^{2} average]. The results obtained by Cahalan et al. (1994a) correspond to mostly cloudy data between the overcast and all scenes averages, and at 63° solar zenith angle these results agree well with their 15% bias. The 8% bias found in Marshak et al. (1995) is for overcast clouds at an averaging scale between Tables 1a,and 1b and also agrees well with the present results.

Examination of the relative magnitude of the IPA–2D versus PP–2D bias shows that the effect of horizontal radiative transport on albedo bias is largest for scattered clouds and smallest for overcast clouds. IPA–2D bias accounts for 32% (rms over the scenes) of the bias for the scattered cloud class, 11% rms for broken cloud fields, and 7% rms in the overcast clouds. The overcast result is identical to the relative bias contributions found in Cahalan et al. (1994b).

Recall that this dataset is not meant to represent the global distribution of clouds; these bias numbers should therefore not be viewed in any way as global averages. The data should be representative of marine boundary layer clouds, however.

## 5. Conclusions

The albedo bias of the independent pixel approximation (IPA–2D) and plane-parallel method (PP–IPA) for oceanic boundary layer clouds have been assessed using realistic cloud fields from 45 Landsat scenes including broken cloudiness and a two-dimensional SHDOM radiative transport method. The modeling includes conservative scattering, ignores gas absorption, and uses a black surface. Tests indicate that these assumptions have only a minor effect on the results shown.

PP–IPA cloud albedo biases up to near 0.30 were found in this study. They are strictly positive due to the shape of the albedo versus optical depth relationship and are best correlated by the standard deviation of the optical depth distribution *σ*_{τ}. IPA–2D bias, on the other hand, ranges from very small positive values to near −0.12 for extremely variable cloud fields. Cloud aspect ratio, as defined in this study, is found to best correlate this bias. While the IPA–2D bias is generally much smaller than the PP–IPA bias, it is not always negligible. Overall, it accounts for 20% rms of the bias. Negative IPA–2D bias occurs due to reflection from the cloud sides in the 2D solution. The dependence of albedo bias on five cloud parameters (cloud fraction, mean cloud optical depth, cloud aspect ratio, standard deviation of optical depth, and gamma function shape parameter) is also examined. The most useful parameters for predicting bias due to cloud field horizontal inhomogeneities are found to be the standard deviation of the optical depth distribution, the gamma function parameter, and the cloud aspect ratio. Finally, the bias is examined at larger mesoscale and climatological scales, and found to be in good agreement with earlier work on overcast and mostly cloudy cases. New results are given here for broken and scattered cloud fields, which show the albedo bias can be even larger in these conditions, depending on the within-cloud variability of the cloud field.

## Acknowledgments

Paul Stackhouse, of Colorado State University, provided valuable assistance in learning to run the SHDOM model. Lindsay Parker, of Lockheed, processed the Landsat reflectances into the cloud optical depths and threshholds used in this study. Robert Arduini, also of Lockheed, provided the Legendre coefficients of the Mie phase function. Three anonymous reviewers provided many helpful comments.

## REFERENCES

## Footnotes

*Corresponding author address:* Dr. Lin H. Chambers, Atmospheric Sciences Division, NASA/Langley Research Center, Hampton, VA 23681-0001.

Email: L.H.Chambers@larc.nasa.gov