## Abstract

Different approaches for parameterizing the effects of vertical variability of cloudiness on radiative transfer are assessed using a database constructed from observations derived from lidar and millimeter cloud radar data collected from three different locations. Five different methods for dealing with the vertical overlap of clouds were incorporated into a single radiation model that was applied to the lidar/radar data averaged in time. The calculated fluxes and heating rates derived with this model are compared to broadband fluxes and heating rates calculated with the independent column approximation using the time-resolved cloud data. These comparisons provide a way of evaluating the effects of different overlap assumptions on the calculation of domain-mean fluxes. It was demonstrated how two of the most commonly used overlap schemes, the random and maximum-random methods, suffer a severe problem in that the total cloud amount defined by these methods depends on the vertical resolution of the host model thus creating a vertical-resolution-dependent bias in model total cloudiness and radiative fluxes. A new method is introduced to overcome this problem by preserving the total column cloud amount.

Despite these problems, the comparisons presented show that most methods introduce a relatively small bias with respect to the single-column data. This is largely a consequence of the nature of the cloud cover statistics associated with the lidar/radar observations used in this study and might not apply in general. Among the three best-performing methods (random, overcast random, and maximum random), the more commonly used maximum-random method does not perform significantly better than the other two methods with regard to both bias and rms error despite its relative high computational cost. The comparisons also reveal the nature and magnitude of the random errors that are introduced by the subgrid-scale parameterizations. These random errors are large and an inevitable consequence of the parameterization process that treats cloud structure statistically. These errors may be thought of as a source of noise to the general circulation model in which the parameterization is embedded.

## 1. Introduction

The climate of Earth is largely governed by its reservoirs of water, the exchanges of water between these reservoirs, and modifications of heat associated with these exchanges. A quantitative understanding of the way water cycles between these reservoirs, the so-called hydrological cycle, and understanding the influence of this cycle on the energy budget of Earth are fundamental for understanding climate and ultimately climate change.

The mass of water in clouds constitutes a minute fraction of the total mass of water in the atmosphere, let alone on Earth. However, this tiny amount of water exerts an enormous influence on the energy balance of the planet (e.g., Stephens 1999). For example, condensed water in the form of precipitation-sized particles is an essential source of energy that fuels the circulation systems of the atmosphere and oceans and is the fundamental supply of fresh water to life on Earth. Absorption of radiation by condensed water in the form of liquid water droplets and ice crystals is a further critical source of energy sustaining these circulations.

Uncertainties in the way clouds are specified in climate models and, in particular, how the radiative processes associated with clouds are dealt with is widely accepted as one of the principal obstacles preventing accurate climate change prediction (e.g., Webster and Stephens 1983; Cess et al. 1989; Senior and Mitchell 1993; Houghton et al. 2001; see the Web site http://www.ipcc.ch). Over the past decade, significant resources have been broadly directed toward addressing these uncertainties although relatively little effort per se has focused specifically on the radiation parameterization schemes used in climate models.

These parameterization schemes involve three main levels of simplifying approximations.

The first and most basic level of approximation concerns the method for solving the radiative transfer equation. Methods used are simple by design, set on one-dimensional geometry, and in some instance omit certain processes (e.g., Stephens 1984). It is a simple matter to assess the accuracy of these solution methods and a number of papers have been devoted to this topic, such as the two-stream model papers of Meador and Weaver (1980) and King and Harshvardhan (1986), among others. These studies, however, do not map the solution errors to broadband flux errors as needed for judging the merits of these solution methods for climate applications. Stephens et al. (2001) attempt such an assessment of the two-stream solutions and the intercomparison activities of Barker et al. (2002) as well as Ellingson and Fouquart (1991) offer further relevant climate model tests.

The next level of approximation concerns the specification of the optical properties of the atmosphere. These are the input variables required by the radiative transfer solution module and include the layer optical depth and associated scattering properties. The difficulty confronting the evaluation of these different parameterization methods varies according to the processes being considered. For example, it is relatively straight forward to determine how well gaseous absorption is handled by comparison to existing spectroscopic line-by-line data and appropriate models (e.g., Barker et al. 2002). The evaluation of the parameterizations of cloud optical properties is far more problematic and less developed. Aerosol property parameterization is also an area of active, ongoing research.

The third level of parameterization involves the treatment of subgrid variability of the optical properties and the effects of this variability on the radiative transfer. Subgrid variability is of primary concern in the parameterization of cloud–radiation interactions. This problem is generally approached by factoring effects of horizontal variability and vertical variability. Treatment of the effects of horizontal variability is mentioned by Harshvardhan and Randall (1985) and addressed theoretically by Stephens (1988). Other than the studies of Rossow and Cairns (1995), Tiedtke (1996), and Li and Barker (2002), the parameterization of these effects in GCMs has barely been addressed. Effects of vertical variability of cloudiness are parameterized via statistical treatments of the vertical overlap of clouds. A number of empirical schemes have been introduced to address this problem over the past 20 years (Geleyn and Hollingsworth 1979; Liang and Wang 1997; Collins 2001, and others). As we show below, the use of these methods raises serious concerns for at least two reasons. The first major area of concern is that all methods are empirical with no physical basis. It is particularly troublesome, given their empirical basis, that the use of these methods has continued for so long with practically no evaluation. The second concern emerges from this study as well as the recent study of Barker et al. (2002). We will show that certain overlap parameterizations can introduce intolerably large bias errors in the estimation of areal averaged fluxes. These biases are typically and practically tuned by adjusting the model so that the statistics of simulated top-of-atmosphere fluxes at large space and time scales match those derived from satellite measurements, thus burying the real problem deeper into the cloud–radiation parameterization problem.

Despite the potentially large errors introduced by the parameterization of subgrid variability, very little attention has been placed on systematically evaluating existing schemes used in models today. The study of Barker et al. (2002) provides an initial glimpse at the magnitude of the problems associated with these parameterization methods and points to the need for a more systematic and critical assessment of the uncertainties introduced by them. The purpose of the present paper, and one to follow (N. B. Wood et al. 2003, unpublished manuscript, hereafter Part II), is to provide such a systematic assessment. The present paper focuses on the topic of cloud overlap and seeks to quantify both the systematic errors and noise introduced by different overlap methods. This assessment is based on a large number of different cases derived using real cloud data. The assessment introduced is constructed around the radiation scheme introduced by Stephens et al. (2001) (hereafter BugsRad) and differs from other recent studies, such as those of Hogan and Illingworth (2000) and Mace and Benson-Troth (2002), that report on cloud overlap statistics. The present study seeks to quantify the effects of different overlap parameterizations in terms of quantities that matter for climate models, namely broadband radiative fluxes and heating rates.

The following section outlines the salient features of the BugsRad model used to simulate the radiative fluxes necessary for this assessment. The implementation of five different overlap schemes in BugsRad is described in section 3. The performance of these schemes is presented in section 4 for two reference cases introduced by Barker et al. (2002) as a continuation of the Intercomparison of Radiation Codes for Climate Models (hereafter ICRCCM-III). Section 5 presents the approach for developing a far more extensive test database for the issues described. Five 3-month periods of Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) program cloud radar data taken from three different ARM sites are used to construct cloud distributions that when inputted into the radiation scheme provide independent column fluxes and heating rates. These column values are averaged over fixed domains and compared with the GCM version of the code that uses matching domain-averaged cloud properties and cloud amount. The results of these comparisons are presented in section 6 and the final section of the paper provides a discussion of the results, key conclusions, and the broader implications of the study.

## 2. The BugsRad radiation model

A radiation model has been developed for use in both global-scale and cloud-resolving models. This radiation model, BugsRad, has been constructed from the two-stream solution of the radiative transfer equation set on the geometry of a plane-parallel atmosphere (Stephens et al. 2001). The basic difference between the version of the model used in general circulation climate models (GCCMs) over that used in cloud-resolving models are the subgrid-scale parameterizations introduced in the former. The version of the model implemented in cloud-resolving models (we will refer to this as the single-column model), unlike the GCM counterpart, has no fractional cloudiness (thus any given layer is either fully clear or overcast) making cloud overlap parameterization moot. This single-column version of the model has been used in other studies including Miller and Stephens (2001) and L'Ecuyer and Stephens (2003) among others and serves as the basis for the radiation budget analysis method developed for the National Aeronautics and Space Administration (NASA) CloudSat program (Stephens et al. 2002).

BugsRad exploits the two-stream approximation to solve the radiative transfer equation for both solar and IR fluxes and uses the same overlap assumption for both. This solution accommodates the processes of gaseous absorption, molecular scattering, and particle absorption and scattering. One advantage of the solution of BugsRad is the demonstrated linear dependence of the model computation performance with model layer count (Stephens et al. 2001). The correlated-*k* distribution method is used to treat the gaseous absorption process and is based on the *k*-distribution data of Fu and Liou (1992). The parameterization of the scattering (and absorption) properties of clouds is based on the use of anomalous diffraction theory first introduced for such an application by Stephens (1984) and developed since by many others (e.g., Mitchell et al. 1996).

As described in Stephens et al. (2001), the two-stream solution for diffuse fluxes can be written in its interaction form as

for the *level* diffuse fluxes *F*^{±}_{n,n+1} given the respective *layer* reflection *R*_{n} and transmission *T*_{n}. The source Σ^{±}_{n} coefficients have the form

where 𝒫&⊚ _{,IR} are the particular solutions of the two-stream equations representing contributions by solar sources (indicated by the ⊚ subscript) and longwave (IR) sources. The specific forms of these particular solutions depend on the details of the two-stream assumptions. Examples of their form are given, for example, in Stackhouse and Stephens (1991) and Stephens et al. (2001). We also adopt the nomenclature such that *n* or *n* + 1 refers to layer quantities when subscripted and level quantities when in parentheses. The ± superscript indicates upward (+) and downward (−) quantities. In dealing with cloud overlap, two sets of quantities must be specified: one set defines the clear-sky properties (hereafter identified by the superscript *f*) and the second set defines cloudy-sky properties (hereafter identified by the superscript *c*). As described below, a given overlap scheme determines how these properties are mixed in any given layer.

For later purposes, we write the particular solution associated with solar sources as

where the layer coefficients *p*&⊚ _{,n} are specific to the given two-stream model, *τ*(*n*) is the optical depth of the entire atmosphere above level *n*, and *μ*_{o} is the solar zenith angle. With this form of this particular solution, the source function becomes

where Δ*τ*_{n} = *τ*(*n* + 1) − *τ*(*n*) is the optical depth of the *n*th layer.

Equation (1) can be written in the multilayer form as described by Stephens et al. (2001) as

yielding a tridiagonal system of 2*n* + 2 equations. The solutions of (4) produce *diffuse* fluxes *F*(1, 2, … , *N* + 1). For solar fluxes, the flux associated with the direct beam must be added to (1). In this case, the downward *global* flux becomes

where *F*_{⊚} is the extraterrestrial solar flux *F*_{⊚} incident at a solar zenith angle *μ*_{⊚}. Hereafter we refer to the direct-beam flux as

The reader is referred to Stephens et al. (2001) for specific details on the form of the coefficients.

## 3. Cloud overlap parameterizations

Four overlap methods currently in use in different climate models have been implemented in BugsRad. These methods are described below as well as a new method that attempts to overcome problems inherent in two methods. We also incorporated the method of Collins (2001), which in some ways is similar to that of Liang and Wang (1997). The approach of Collins (2001) is to divide each GCM column into independent subcolumns, set the cloud fields in the subcolumns so that a particular overlap configuration is produced for the column as a whole, perform independent column calculations for each subcolumn, then derive area-weighted averages of the fluxes and heating rates for all subcolumns. Since this method essentially reproduced the results reported below for the maximum-random method, the specific results are not shown.

As we show below, a severe problem with two of the more common methods in use today is the dependence of the total cloud amount on the vertical resolution of the host model. This tends to create a vertical-resolution-dependent bias in the implied total cloud amount and a subsequent bias in derived radiative fluxes.

The specific implementation of these five methods in BugsRad is now described. The most complicated implementation of all methods is that of the maximum random overlap (MRO) method and the specific details of this implementation are given in the appendix. The complexity of the MRO arises from the fact that the cloud properties of any given layer are not independent of the properties of adjacent layers. This coupling of adjacent layers effectively increases the dimensionality of the radiative transfer problem producing a system of 4*n* + 2 equations compared to the 2*n* + 2 system used to produce the single-column model solutions. This increased dimensionality increases the computational cost of the method relative to other overlap schemes. The other methods introduced below in one way or another effectively decouple the properties of one layer from the properties of adjacent layers thus facilitating the solution via the simpler and more computationally efficient single-column solution and thus preserving the 2n + 2 system of equations noted above.

### a. The random overlap method

The random overlap method (ROM) assumes the cloudiness in any given layer is independent of the cloudiness of other layers. Although unrealistic, this assumption offers the computational benefits noted. The actual implementation of this scheme is simple, involving the weighting of the reflection, transmission, and source coefficients for a given layer of index *n* by the cloud fraction *C*_{n}. In the case of the reflection coefficient, the layer reflection becomes

where the superscripts *c* and *f*, respectively, refer to cloud and cloud-free portions of the layer. As a consequence, it is unnecessary to mix clear- and cloudy-sky solutions.

The total cloud amount of the column follows as

for *N*_{c} layers of cloudiness. Equation (5) implies that the total cloudiness depends on the number of cloudy layers that define the column and thus the vertical resolution of the model. As a consequence, the total cloud fraction, *C*_{tot}, implied in these schemes typically exceeds the actual total cloud fraction by an amount that systematically increases as *N*_{c} increases.

### b. The maximum random overlap

The maximum random overlap method was introduced by Geleyn and Hollingsworth (1979) as an attempt in part to do away with the unrealistic assumptions underlying the random overlap method. The maximum random overlap method assumes that cloudiness in the adjacent layers above and below a given layer in question are maximally overlapped with the cloudiness of that layer. Elsewhere clouds in layers separated by layers of clear air above and below are randomly overlapped. The method, however, is computationally cumbersome for reasons mentioned and it does not alleviate the problem of a lack of conservation of total cloud cover.

For the maximum random overlap method, we combine the interaction principle (1) for the clear and cloudy sky to produce a modified version of (4),

where 𝗔 is a (4*n* + 2) × (4*n* + 2) matrix of the form

which is a band-diagonal matrix constructed by the 4 × 2 block matrices:

where *T*^{c/}_{n}^{f} and *R*^{c/}_{n}^{f} are the layer reflection and transmission functions defined for cloudy/cloud-free portions of the *n*th layer. The cloud overlap coefficients are given by (e.g., Geleyn and Hollingsworth 1979)

and these coefficients are calculated recursively beginning with *O*_{11} = 1 − *C*_{1} and *O*_{31} = 1. An interpretation of these overlap coefficients is given in Geleyn and Hollingsworth (1979). The block elements 𝗯_{1} and 𝗯_{2} of matrix 𝗔 accomodate boundary conditions. The first is a 4 × 2 matrix such that 𝗯_{1} = 𝗮_{2} whereas the second is a 2 × 4 matrix of the form

where *α* is the prescribed surface albedo.

The solution of (6) requires the inversion of the 𝗔 matrix to determine the vector of required clear- and cloudy-sky fluxes of the form

given the specified source function vector:

where *δ*_{IR} = 1 for longwave fluxes and zero otherwise and *B*_{n+1} is the blackbody emission from the surface and *α* is the IR surface albedo. The source elements *S*^{f/c}_{n} are weighted combinations of cloudy and cloud-free sources. The thermal sources are trivially defined in terms of the cloud fraction of layer *n* according to

The solar sources, however, are more complex, being influenced by the cloudiness above the layer in question, and require the direct-beam component incident on the top of the layer. The global solar fluxes also require the addition of this direct-beam component to the diffuse fluxes obtained by inversion of (6). This direct beam is computed recursively beginning with a weighting of the particular solutions

### c. The threshold random method

The threshold random method (TRM) was introduced by Hansen et al. (1983) and is the procedure currently used in the Goddard Institute for Space Studies (GISS) climate model. The method reduces layers containing partial cloudiness to layers that are either purely overcast or clear thus reverting the model to the single column version described above. The procedure generates a random number between 0 and 1 at each grid point and compares this number to the cloud fraction at each layer. Full cloud cover is assigned for cloud fractions that exceed the random number otherwise the layer is assumed to be clear. There is no a priori justification for this method. The method clearly alters the specified cloud cover, increasing coverage to overcast on some occasions and removing clouds entirely on others. The expectation is that when averaged over many realizations, the resultant mean fluxes contain insignificant biases. The comparisons below support this assumption but the method is inferior on an single case-by-case basis.

### d. The scaling method

The scaling method was introduced by Briegleb (1992) in an early version of the radiation scheme used in the National Center for Atmospheric Research (NCAR) Community Climate Model. The intent of the method is to approximate the random overlap approach described above. As with the TRM, the scaling method effectively converts layers containing partial cloud cover to overcast layers but in this case by scaling the optical depth of the cloudy portion of the layer by the cloud amount. This scaling has the empirical form

where, as before, the fractional area of cloud in layer *n* is *C*_{n} and the optical depth of the cloudy portion of the layer is *τ*_{c}. Although this method was originally used to calculate shortwave fluxes, we also apply this method in BugsRad to calculate longwave fluxes for these tests.

### e. A new hybrid overcast-random overlap scheme

The new scheme introduced here divides the model grid into three regions: a purely clear-sky region of fractional area 1 − *C*_{tot}, a purely overcast region defined by the fractional area *C*_{over}, and a region *C*_{ran} wherein the clouds are purely randomly overlapped such that *C*_{ran} = *C*_{tot} − *C*_{over}. The overcast region thus corresponds to a region in which any layer containing cloud is filled completely by that cloud.

This scheme requires the specification of the total cloudiness *C*_{tot}, which, for the current study, is derived from the cloud radar data as discussed below. Given this cloud amount, the portion of overcast sky follows from the solution of

where the product is taken over the layers for which *C*_{n} > *C*_{over}. Although the cloudy-sky portion of the solution requires the combination of two separate solutions, both are of the single-column-mode form wherein the properties of each layer remain independent from the properties of the other layers. One solution applies to purely overcast layers for that fraction of the domain of fractional area *C*_{over} and the second solution uses the purely random approximation for the fractional area *C*_{ran}.

The domain-averaged flux is derived by combining the three solutions according to

This scheme has some similarity to that of Morcrette and Jakob (2000) except that an overcast portion is now singled out. The proportion of the grid box divided into overcast, random, and clear is determined by the total cloudiness, which, for the applications reported in this paper, is given from the data used in this study. However, the implementation of this method in a climate model that does not provide *C*_{tot} but only layer cloudiness is to be described in a forthcoming paper.

A relevant yardstick for assessing the different overlap schemes is their relative computational cost. Figure 1 presents a comparison of the computational time it takes to execute a complete calculation of broadband solar and IR fluxes. This computational time is normalized to the SM method for reference. As expected, the time to execute BugsRad with the MRO exceeds all other methods, being almost 2.5 times slower than any of the other four methods. Also noteworthy is the computational cost of the new overcast-random scheme, which is comparable to the other three schemes.

## 4. The ICRCCM-III comparisons

In the ICRCCM-III study, Barker et al. (2002) compare the performance of several GCM solar radiative transfer codes applied to a limited selection of cases developed using different cloud configurations. These prescribed cloud conditions were derived from a selection of outputs of three-dimensional cloud-resolving models. We apply BugsRad with the overlap schemes introduced above to two of the ICRCCM-III cases merely for reference purposes. The results obtain with these different methods are similar to the results reported in Barker et al. (2002) obtained for the different models and implementation of the overlap methods considered in this paper. The ICRCCM-III cases are from the 3D simulations of tropical convection by Grabowski et al. (1998). One case (case A) consists of nonsquall clusters of organized convection composed of liquid-only clouds in as low-shear environment and the clouds are distributed in a quasi-organized array reminiscent of trade wind cumulus. The second case (case B) is of a more complex squall line consisting of both water and ice phase clouds.

The data provided to BugsRad are of the form of domain-averaged profiles of cloud liquid water and ice water contents and cloud fractional cover. In the original study of Barker et al. (2002), the atmosphere is divided into 42 layers. For case A, *N*_{c} = 25 and the total cloudiness is 0.462. For case B, *N*_{c} = 25, and the total cloud cover is 0.251 when liquid only is assumed and *N*_{c} = 32 and *C*_{tot} = 0.999 when liquid plus ice are assumed.

Figures 2a–d present the top-of-the-atmosphere (TOA) albedo and the profile of solar heating derived from BugsRad with the inputs as described. As in Barker et al., (2002), the TOA albedos are presented as a function of the cosine of the solar zenith angle (*μ*_{0}) and the solar heating rates are those calculated for *μ*_{0} = 0.5. Also shown for reference are the benchmark values derived from a 3D Monte Carlo model as well as the independent column values. From these comparisons, we note the following: (i) As shown elsewhere, the independent column fluxes closely approximate the benchmark fluxes for the two cases. (ii) The performance of the different overlap schemes applied in one single model produces a spread in albedo results, relative to the benchmark, that is similar to model-to-model differences reported by Barker et al. (2002). (iii) The overlap schemes introduce significant albedo biases relative to the benchmark. The sign of these biases differs from case to case. The exception is the random overlap scheme, which introduces a positive albedo bias for both. (iv) With the exception of the TRM, the effect of overlap assumption on solar heating rates is small.

Figure 3 presents a similar comparison of albedo but now only for the random, MRO, and the overcast-random methods applied to case A (upper two panels) and a version of case B where the ice cloud layers are removed (lower two panels). These results are presented as a function of the changing vertical resolution of the model as indicated by the given total number of cloudy layers. The resolution was altered by binning the original data into progressively thicker layers while preserving total cloud water and optical thickness. These comparisons highlight the problems already noted, namely that the performance of some schemes is dependent on the vertical resolution the model. For example, the performance of the random and maximum random methods deteriorates as the number of cloudy layers is increased. The finer is the resolution of vertical cloudiness, the less relevant are assumptions about independence from one layer to the next. The increasing cloudiness with layer count produces a corresponding increasing albedo bias. The performance of the overcast random method, now constrained to the given total cloudiness, is improved over the random method but even in this case a bias with respect to the benchmark results remains. One source of this bias is the effect of horizontal variability, which is not included in our simulations by design. Barker et al. (2002) use these same cloud fields and show that even when overlap is modeled exactly, the bias that results from treating horizontally inhomogeneous clouds as homogenous can be significant. The issue of horizontal variability is addressed in Part II.

## 5. Test datasets

The ICRCCM-III test comparisons described above reveal the nature of the overlap errors on solar fluxes and heating. Since these errors depend on the cloud configuration under study, more extensive assessment is needed. We now describe how this very limited type of assessment is extended through the creation of a more extensive test dataset based on observed clouds rather than a few modeled cloud fields.

### a. The ARM cloud data

Data extracted from the Active Remotely Sensed Clouds Locations (ARSCL; Clothiaux et al. 2000) developed under the auspices of the U.S. Department of Energy Atmospheric Radiation Measurement (ARM) program provide the basis for the test datasets developed in this study. This dataset is a composite of lidar, millimeter-wave cloud radar (MMCR), and ceilometer data. The specific data extracted from ARCSL for this research include the cloud locations and cloud radar reflectivities. The cloud radar reflectivities correspond to a frequency of 35 GHz (W band) and to the radar pointed in a fixed zenith direction. In this mode of operation, the radars have a vertical range resolution of 90 m (up to about 20-km altitude) and a sensitivity of −48 dB*Z* at 5 km. These radars have operated quasi-continuously at three sites for more than a year, Nauru (deployed in 1998), at the Southern Great Plains (SGP) site in Oklahoma (deployed in 1996), and at the North Slope of Alaska (NSA) site at Barrow, Alaska (deployed in 1998). A fourth MMCR was recently deployed at the Darwin site in late 2001 but none of these data are used in this study.

Figure 4a provides an example of profiles of W-band reflectivity (*Z*) collected over a single 3-h period from the Nauru radar. Individual profiles are composited together to produce this time–height cross-section segment. All radar data from the three ARM sites were composited into similar 3-hourly segments, hereafter referred to as domains. A total of 2005 domains were created from five three-month periods of radar data.

If we take 10 m s^{−1} as a coarsely representative value of the advecting wind across the 3-h domain, then the equivalent horizontal dimension of this domain would be about 100 km. This domain is slightly smaller than current climate model resolutions but slightly larger than current global NWP models and thus falls in between the two types of global models. Each 3-h domain is constructed from 1080 radar reflectivity profiles taken at 10-s intervals and the variability from profile to profile within these domains can be thought of as representing the subgrid-scale distribution of cloudiness that is represented by this approximate 100-km model grid box. The domain-mean cloud amount and domain-mean cloud properties derived from these profiles are the necessary input into the overlap approximations and are meant to resemble the GCM gridpoint data.

Figures 5a and 5b present a composite of the cloud amount and the number of cloudy layers of the 2005 cases. The probability density function (pdf) shown in Fig. 5a has the form of a bimodal distribution analogous to other cloud cover datasets (e.g., Rossow and Schiffer 1999) indicating a prevalence for the domains to be either clear or overcast. The pdf of the number of cloud layers is broadly distributed about a maximum that corresponds to two layers of clouds.

### b. Independent column fluxes

The MMCR radar reflectivity data were converted to cloud liquid and ice water contents following the procedure outlined in Miller and Stephens (2001). The reflectivities measured at temperatures below 253 K were converted to cloud ice water content (IWC) using the *Z*–IWC relationship of Brown et al. (1995) that corresponds to the case of particle effective radius *r*_{e} < 50 *μ*m. For temperatures exceeding 253 K, reflectivities were converted to cloud liquid water content (LWC). For simplicity, the mixed-phase region, nominally occurring between the freezing level and 253 K, was assumed to be liquid. The LWC was determined as in Miller and Stephens (2001) except where spuriously large LWC values are created in regions of precipitation most notably occurring below the melting layer where reflectivities typically exceed 0 dB*Z.* We eliminate these spurious values in a simple way by screening for precipitation below the melting level and only converting those reflectivities that have values below 0 dB*Z.*

The cloud LWC and IWC data derived using the 45-m vertically resolved reflectivity data were then averaged in the vertical to produce averaged properties that matched the 19 vertical layers of the GCM radiation code. This averaging procedure resulted in the data presented in Fig. 4b.

Since the focus of this study is directed toward assessing the effect of different overlap assumptions on the radiative transfer, we also remove the complicating effects of in-cloud horizontal variability. These effects are removed by replacing the individual profiles of LWCs and IWCs at the specified levels with the in-cloud average values derived for that level. Comparison of Fig. 4c with Fig. 4b shows how this horizontal homogenization of the data affects the LWC and IWC distributions.

The profiles of Fig. 4c, together with temperature and relative humidity derived from matched radiosonde data, are the required inputs of the single-column version of BugsRad. The cloud optical properties used in this version of the model are derived from the IWC and LWC as described elsewhere (Stephens et al. 1990) with *r*_{e} = 10 *μ*m for water clouds and *r*_{e} = 30 *μ*m for ice clouds. For each domain, we obtain a distribution of independent column fluxes and radiative heating rates averaged to produce domain-mean values. We hereafter consider these domain-mean values to be a form of truth for assessing the different overlap methods.

Before presenting the results, it is necessary to comment on the evaluation dataset created in this way. First of all, the cloud distributions are two-dimensional. It is neither practical nor possible to construct extensive datasets on 3D cloud distributions from existing observational data. Second, the horizontal homogenization of cloud properties removes other effects discussed in Part II.

## 6. Comparison results

The GCM version of BugsRad requires domain-averaged profiles of in-cloud liquid and ice water contents and the profiles of cloud amount, temperature, and humidity. The overcast-random method also requires the additional total cloud fraction information, *C*_{tot}, which is determined from the cloud radar data. Solar zenith angle is also required and the values of 20° and 60° were used to match the independent column fluxes. The output from BugsRad are profiles of broadband fluxes and related heating rates that supposedly represent the domain averages of these quantities. The key results of this study are then presented in the form of a comparison of these fluxes and heating rates to the equivalent domain-averaged values derived by integrating the respective independent column quantities.

### a. Exitant boundary fluxes

An example of the comparison of the exitant fluxes are shown in Fig. 6 in the form of a scatter diagram. The fluxes shown in this case correspond to those derived using the Nauru MMCR data, one point for each 3-h domain for each parameterization method. Figures 6a and 6b correspond to the longwave fluxes leaving the atmosphere whereas Figs. 6c and 6d are the respective solar fluxes for the given 20° solar zenith angle.

Domains with more cloudiness generally show increased scatter relative to the independent column fluxes. The threshold random method shows the largest errors, at times on the order of several hundred watts per square meter for solar fluxes. These domains appear to be partly cloudy, perhaps with one or two fractionally cloudy layers, and the threshold random method converts them to purely clear or purely overcast.

The scaling method appears to underestimate the albedo of the cloudy sky. For example, for downwelling solar fluxes at the surface (Fig. 6d) of 800 W m^{−2}, the scaling method estimates the flux to be about 950 W m^{−2}. Under more heavily clouded skies, with downwelling solar fluxes at the surface of 200–400 W m^{−2}, the scaling method overestimates the albedo, giving fluxes that are often 50% or less of the independent column values. A similar pattern of scatter, but of smaller magnitude and with sign reversed, occurs in the longwave fluxes. As can be seen in the plots, the scatter patterns for the scaling method are completely different from those of the random method, which the scaling method is intended to approximate.

For the three remaining methods, the longwave fluxes are in good agreement with the independent column fluxes for the full range of sky conditions.

The three methods however differ for shortwave fluxes. The maximum random method exhibits moderately large errors for moderate to high amounts of cloudiness. The albedo of the cloudy sky appears to be underestimated by this method under these conditions. This result could occur if there were a tendency for the method to maximally overlap cloudy layers to a greater degree than actually occurs in the radar data. This issue is not investigated in this study, but Hogan and Illingworth (2000) show evidence that in vertically continuous cloud layers, the overlap configuration transitions from maximum to near random as the vertical separation increases to about 4 km. This transition is not implemented in the maximum random overlap method.

For near-clear skies, the random method overestimates the albedo of the cloudy sky. Downwelling solar fluxes at the surface are thus underestimated and upwelling TOA solar fluxes are overestimated compared to independent column results. These cases probably consist of domains with few cloudy layers that in reality overlap more so than that predicted by simple random overlap.

The overcast-random method shows low scatter and little bias over a large range of cloudy conditions. For solar fluxes, the method exhibits a slight tendency to underestimate the albedo of the cloudy sky for high amounts of cloudiness. For example, when downwelling solar fluxes at the surface (Fig. 6d) are between 40 and 200 W m^{−2}, this method estimates the fluxes to be between 140 and 250 W m^{−2}. These cases are roughly 1% of the total number of cloudy domains examined at Nauru.

Results obtained using data from the other two locations are similar to those presented above for Nauru. Tables 1a and 1b, respectively, and Fig. 7 compare the exitant TOA and surface fluxes composited for all sites for each of the five different overlap methods. The flux information contained in Tables 1a and 1b is presented in terms of a bias error, an rms error, and the maximum and minimum outlier values derived from the 2005 domains. The comparisons shown in Fig. 7 similarly contrast the bias errors and random errors.

The results presented in this way can be summarized as follows:

Contrary to the ICRCCM-III cases, the random and overcast-random methods introduce an almost negligible bias. While at first sight this is unexpected, it is a consequence, in part, of the nature of the cloud cover statistics noted in reference to Fig. 5a. The small bias of the random overlap method also reflects the smaller number of layers occupied by clouds in the ARM data compared to the ICRCCM-III cases as well as the removal of the effects of horizontal variability.

It is also noteworthy that the MRO method does not perform substantially better than either the overcast random or random methods with regard to both the bias and the rms error but that it produces a tighter error distribution overall as is evident by comparison with the maximum and minimum outliers.

When averaged over a large number of cases, the TRM introduces a relatively small bias as originally hypothesized by Hansen et al. (1983). The method, however, is inferior in all other respects, most notably in terms of the outlier values of fluxes and the rms error introduced.

### b. Radiative heating rates

The effects of overlap parameterization on atmospheric radiative heating are presented in the form of a scatter diagram of differences plotted as a function of altitude and the rms heating rate error is described in the form of a histogram. These scatter diagrams of heating rate differences are shown in Fig. 8 for all 2005 cases and the rms heating rate errors of each method are presented in Fig. 9 after being binned into four broad atmospheric layers.

The effects on heating rate occur predominately in the middle and lower troposphere, as would be expected. These regions are either influenced by cloud directly or are influenced by clouds at higher altitudes.

Overall, the heating rate bias errors are small for all five methods (less than 0.05 K day^{−1}) but the rms errors are considerably larger. In terms of the rms errors, the methods fall into two distinct groups. Both the TRM and SM introduce rms errors on longwave heating that exceed 0.5 K day^{−1} and similarly large solar heating rate errors for the case of *θ*_{o} = 20°. By contrast, the random, MRO, and overcast-random methods introduce modest rms heating errors of less than 0.1 K day^{−1}. These heating rate rms errors seem to be strongly influenced by the filling of cloudy layers (as is done by the TRM and SM) or the emptying of cloudy layers (as is done by the TRM).

## 7. Discussion and conclusions

This paper attempts to assess different approaches for parameterizing the effects of vertical variability of cloudiness on radiative transfer. To this end, a database was constructed from observations derived from lidar and millimeter cloud radar data collected from three ARM sites.

Five different treatments of the vertical overlap of clouds were incorporated into a single radiation model that was applied to the lidar/radar data averaged in time. The calculated fluxes and heating rates derived with this model are compared to broadband fluxes and heating rates calculated with the independent column approximation using the time-resolved cloud data. These comparisons provide a way of evaluating the effects of different overlap assumptions on the calculation of domain-mean fluxes.

Two of the most commonly used overlap schemes, the random and maximum-random methods, suffer a severe problem in that the total cloud amount defined by these methods depends on the vertical resolution of the host model. This tends to create a vertical-resolution-dependent bias in the implied total-cloud amount and a subsequent bias in derived radiative fluxes. A new method, referred to as the overcast-random method, is introduced to overcome this problem by preserving the total column cloud amount.

Based on the comparisons using a 19-layer atmospheric model, both the random and overcast-random methods introduce an almost negligible bias in fluxes and heating rates. This is expected, in part, given the nature of the cloud cover statistics associated with the lidar/radar observations used in this study. The small bias error of the random overlap method also reflects both the relatively small number of layers occupied by clouds in the dataset and the removal of the effects of horizontal variability. The more commonly used MRO method does not perform significantly better than either the overcast random or the random method with regard to both the bias and the rms error despite the fact that the computational cost of the method is almost 2.5 times that of the other methods.

The comparisons also reveal the nature and magnitude of the random errors that are introduced by the subgrid-scale parameterizations. These random errors are large and an inevitable consequence of the parameterization process. These errors may be thought of as a source of noise to the general circulation model in which the parameterization is embedded. A study is currently under way to assess the extent of the influence of this noise on a GCM climate model.

## Acknowledgments

This research was supported by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research, Environmental Sciences Division under Grant DE-FG03-94ER61748 as part of the Atmospheric Radiation Measurement (ARM) program. Data used were obtained from ARM, which is also sponsored by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research, Environmental Sciences Division. We also acknowledge Dr. H. Barker for a number of helpful conversations and supply of the ICRCCM-III data.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Implementation of the MRO Method in a Two-Stream Model

Here we provide more detail on the derivation of Eqs. (6)–(10) and describe methods for solving (6). Figure A1 is a useful reference point for discussion including the layer and level conventions used. In the MRO approach, cloud layers above and below some reference layer are maximally overlapped with respect to clouds in that layer and randomly overlapped when separated by a clear layer. There are six fluxes whose solution must be constructed: the upwelling diffuse fluxes for the clear and cloudy sky and the corresponding downwelling diffuse fluxes and for solar fluxes the additional downwelling direct-beam components in clear and cloudy skies.

The solutions outlined in the body of this paper derive from the interaction principle as expressed by (1a) and (1b). For a given layer, we simply apply this principle to (1a) and (1b) each for the cloudy and cloud-free portions of the layer and then empirically combine these by weighting by cloud fraction. Consider the diffuse fluxes at level *n* and *n* + 1. The two downwelling diffuse fluxes emerging from the *n*th layer at level *n* + 1 follow from (1b), (8a), and (8b) as

and

and the upwelling fluxes at level *n* similarly follow as

When the flux vector is constructed according to (9), then Eqs. (A1)–(A4) map directly into Eq. (6) producing the matrix (7) and source vector (11).

We offer a relatively simple example to facilitate interpretation of (A1) and (A2). Consider the three-layer example of Fig. A1 with cloud fractions of 0.3, 0.5, and 0.8 from top to bottom. Consider the clear and cloudy fluxes at level 2. The diffuse cloudy-sky downwelling component follows from (A1) as

The factor of [3/5] represents the ratio of the cloudy-sky fractions between layers 1 and 2 and weights the upwelling cloudy flux issuing from layer 2, which is then multiplied by the cloudy-sky reflection coefficient of layer 1. The factor [3/10] accounts for the diffuse clear-sky flux entering layer one, which is transmitted by the cloudy region in the first layer. To understand why the upwelling diffuse clear-sky flux is multiplied by a weight of 0, we appeal to the definition of the MRO, embodied by the definitions of the *O*_{i,n}. For the given geometry, the upwelling clear-sky flux can only enter the clear-sky layer immediately above it. In fact as long as the cloud fraction is monotonically increasing in the manner of Fig. A1, the weighting of the clear-sky diffuse fluxes entering into a cloudy portion will be zero.

##### Numerical solution

The 𝗔 matrix in (7) is a band-diagonal matrix of bandwidth 11. This suggests two techniques for its solution: one exact, the other, approximate. The exact technique uses a band-diagonal solver to perform an LU decomposition as described in the book *Numerical Recipes* (Press et al. 1989). This is vastly more efficient than inverting the full 𝗔 matrix, but is still suboptimal because the solver does not fully exploit the sparsity of the matrix: there are still some zero elements in the bands. We replace the 𝗔 matrix in (6) with a transformed matrix 𝗔′ having the following fill pattern:

where each element of the matrix correspondingly relates to each element of the 𝗔 matrix of (7). For example, the *a*_{1,5} element noted in (A6) is the element of the first row, fifth column, of the 𝗔 matrix of (7). The two routines used to solve this new system are “bandec” and “banbks” found in *Numerical Recipes.*

We have also explored other more approximate iterative methods of solution. These require an initial guess and then the system is iterated to within some prescribed level of accuracy. The iterative algorithms used were 1) Jacobi iteration, 2) Gauss–Seidel iteration, and 3) successive overrelaxation. These will now be discussed briefly. We have not implemented these approximate methods in the operational version of BugsRad. To check the correctness of the solution we employed three different tests: the clear-sky limit characterized by the absence of clouds in all layers, the cloudy-sky limit wherein all layers contain cloud (cloud fraction of unity for all layers), and the stacked cloud limit where a fixed, single cloud fraction is applied to all layers. In this case all of the overlap coefficients become unity.

###### Jacobi iteration

We split the attenuation matrix as follows:

Then Eq. (6) may be written as

The Jacobi iterative scheme becomes

The method converges *unconditionally* when 𝗔 is diagonally dominant. If it is not, convergence may be possible but at the expense of more iterations. We have exploited fully the sparsity of the 𝗔 matrix to maximize computational speed.

###### Gauss–Seidel iteration

In the Jacobi iteration, every *k* + 1 iteration is evaluated using only components of the *k*th iteration. In the Gauss–Seidel method values of all previous components are used to perform updates within the same iteration. The iterative scheme is

This is equivalent to iteration with back-substitution. This improves the rate of convergence by about a factor of 2 over the Jacobi iteration. As with the Jacobi iteration, the method is unconditionally stable if the 𝗔 matrix is diagonally dominant. Otherwise, more iterations will be required for convergence, if convergence is indeed attainable.

###### Successive overrelaxation

Both the Jacobi and Gauss–Seidel schemes can be expressed as

We can multiply the right-hand side by a constant factor *w* to enhance the convergence rate. This factor has the effect of reducing the size of the eigenvalues in the iteration matrix. In the case of the Gauss–Seidel scheme, we have

Therefore,

The optimum choice of *w* is problem dependent and is usually found by trial and error. In the IR, we found *w* = 1 to be optimal. In the shortwave, depending on the band, we found that *w* > 1 (e.g., *w* = 1.08), to decrease the number of iterations.

## Footnotes

*Corresponding author address:* Dr. Graeme L. Stephens, Dept. of Atmospheric Science, Colorado State University, Room 408 ATS, Fort Collins, CO 80523-1371. Email: stephens@atmos.colostate.edu