## Abstract

The authors examine 3D solar radiative heating rates within tropical convective–cirrus systems to identify the scales that contribute significantly to the spatial average over a climate model’s grid cell (i.e., its grid mean), and determine their relationship to the cloud field properties (e.g., cloud-top height variation). These results are used to understand the spatial resolution and subgrid-scale cloud property information needed in climate models to accurately simulate the grid-mean solar heating of these systems.

The 3D heating rates are computed by a broadband Monte Carlo model for several regional-scale cloud fields [(400 km)^{2}] whose properties are retrieved from satellite data over the tropical western Pacific. The analyses discussed in this paper have identified two key subgrid-scale features within these systems that largely govern the grid-mean heating rates: the variability in the cloud-top height, and the structure of the cloud edge. These features give rise to hot spots—regions of intense local heating that occupy a small area but dominate the grid-mean value. For example for the fields considered here, 5%–25% of the grid area can contribute 30%–60% of the total heating rate, respectively. Explicitly resolving the hot spots requires a model grid of about (20 km)^{2}–(30 km)^{2} which is smaller than that currently used in general circulation models (GCMs) for weather forecasting and about a factor of 20 smaller than that used for climate studies. It is shown that, unless a grid of ∼(20 km)^{2} is used, GCM-style heating rate calculations that employ a standard cloud overlap-type treatment can significantly overestimate the solar heating aloft and underestimate it below. This might enhance the vertical velocity within the cloud layer and suppress it at cloud base. Thus, over the long term, biases in the GCM treatments of the vertical heating rate might have consequences to cloud evolution and feedback, particularly for clouds in weak local dynamical regimes.

## 1. Introduction

Climate models must explicitly resolve or else parameterize the physical processes that govern the climate system. Toward this goal, we should guide our model development by an understanding of the spatial scales upon which these physical processes operate. This scale-dependent information enables distinguishing those processes that are resolvable by the horizontal resolution of the climate model’s grid (i.e., grid-resolvable) versus those that are smaller than the model grid (i.e., subgrid scale) and therefore must be parameterized. [For reference, a typical general circulation model (GCM) grid size is about 310 km × 310 km.] Such information may also be useful for designing related observational programs.

Clouds are a key component of the climate system and, toward the goal of understanding cloud energetics and its parameterization, we need to determine the scales of clouds and their atmospheric diabatic heating. Such an approach was initiated by Boer and Ramanathan (1997) who used satellite data and a Lagrangian algorithm to track convective-stratiform clouds in the tropical Pacific and observe their evolution and determine the scale dependence of their radiative properties (e.g., visible albedo and brightness temperature). They found, for example, that 95% of the radiatively important clouds are resolvable by a GCM with a grid resolution of about (50 km)^{2}, while only 50% could be resolved by a (250 km)^{2} grid. Further examples of scale dependence studies are Roca and Ramanathan (2000), Wilcox and Ramanathan (2001), and Zhang et al. (1999).

This paper takes the next step by considering the 3D solar radiative heating in tropical convective systems. Atmospheric solar radiative heating rate is strongly effected by clouds, which can have a significant effect on atmospheric dynamics. For example, cirrus solar radiative heating can significantly effect the upper-troposphere thermal structure (e.g., Ramaswamy and Ramanathan 1989). Further, the interplay between the solar radiative heating and the infrared warming/cooling may also affect cloud evolution (e.g., Starr and Cox 1985; Xu and Randall 1995; Petch et al. 1997; Donner et al. 1999).

GCMs compute the mean solar radiative heating rate per model grid cell (i.e., its grid mean) using its computed vertical distribution of fractional cloud cover and a cloud overlap assumption. However, there is significant subgrid-scale cloud variability that cannot be explicitly accounted for in these treatments (e.g., Barker et al. 1999) and could be important for atmospheric energetics. Several recent broadband Monte Carlo modeling studies have investigated the effects of subgrid-scale cloud variability on the radiative properties (e.g., O’Hirok and Gautier 1998), and specifically their effects on grid-mean heating rate calculations (e.g., Barker et al. 1998, 1999; Fu et al. 2000; Oreopoulos and Barker 1999). The emphasis of such studies has been on the resultant grid-mean heating rates, but no one has performed a careful examination of the properties of the 3D, subgrid-scale heating rates.

Our objective is to identify subgrid-scale features that contribute significantly to the grid-mean heating rate, examine their scale dependence, and determine their relationship to the cloud field properties in observed cloud systems. It is our hope that such an understanding would pave the way for objective parameterization of subgrid-scale effects. We will address this question in two steps. The first step is of a taxonomical nature where we will examine the 3D solar radiative heating rates within several convective-cirrus systems. We retrieve the regional-scale, cloud field structure from satellite imagery, and compute the pixel-resolved, 3D radiative fluxes using a broadband Monte Carlo model. This is the first time such 3D, pixel-scale analyses have been performed for multiple cloud scenes.

The second step compares the solar radiative heating rates computed by the Monte Carlo model to what a GCM would compute using a typical cloud representation (e.g., cloud overlap assumption). For this comparison, the resolution of the pixel-level data are systematically degraded to determine the scale at which a GCM must operate to match the high-resolution model calculation. We will show that there is a strong scale dependence in the model accuracy, and that the resolution needed is comparable to that found by the separate analysis given in the first step. In the final section, we discuss the implications of these findings to the parameterization of grid-mean solar radiative heating rates, and to the model-simulated dynamics and cloud fields.

## 2. Methods

### a. Broadband Monte Carlo model

The general approach to the 3D Monte Carlo model follows that outlined in Marchuk et al. (1980). The model treats the scattering and absorption at solar wavelengths by 3D distributions of gases, clouds, and aerosols. The model’s spectral band locations and widths are adjustable. For these solar broadband calculations, we use 25 band intervals between 0.25 and 5.0 *μ*m. The model saturates the water vapor where clouds are present, and computes the photon transport in the 3D media using the Maximum Cross-Section Method (Marchuk et al. 1980). The varying cloud top and bottom heights are treated explicitly (i.e., the Maximum Cross-Section Method permits the cloud boundaries to be different from the model layer boundaries).

Correlated-*k* distributions are used to incorporate gaseous absorption by water vapor, ozone, oxygen, and carbon dioxide (e.g., Goody et al. 1989; Lacis and Oinas 1991). The correlated-*k* distributions are computed offline with line-by-line accuracy using a modified version of a code kindly provided by W. Ridgway of Applied Research Corporation working at the Goddard Space Flight Center of the National Aeronautic and Space Administration [see Vogelmann et al. (1998) for more details]. Up to 21 correlated-*k* coefficients are used per gas per band. Considering the permutations of the correlated-*k* coefficients for overlapping gases, the broadband fluxes are computed using a total of 3132 pseudo-monochromatic calculations.

The computational time for the Monte Carlo calculations is reduced by optimizing the distribution of photons to place the greatest number in the terms that contribute most to the total solution. The sampling method we use is similar to that by Barker et al. (1998), who partition the total number of photons spectrally to place larger numbers in the bands that contain the greatest fraction of the top-of-atmosphere solar flux. Here, we optimize the partitioning of photons among the bands for the broadband heating rate calculation, and between the correlated-*k* coefficients within a band. For example, photons per band are partitioned to place the greatest number in the correlated-*k* permutations with the largest weights. The relative error in estimating a quantity *f* [0, 1] for a particular *k*-coefficient is proportional to [*f*(1 − *f*)]^{1/2}; so the statistical errors should be small when a correlated-*k* coefficient is either small or large (*f* is either close to 0 or 1). Thus, the number of photons applied to the largest and smallest correlated-*k* coefficients in a band is relatively small, which is consistent with their relative contribution to the error.

By this photon partitioning, the Monte Carlo solution time is essentially independent of the number of correlated-*k* coefficients used per band, and the accuracy enhancement afforded by increasing the number of correlated-*k* coefficients is obtained at essentially no extra computational cost. This provides a distinct computational savings relative to standard approaches (e.g., plane-parallel codes such as discrete-ordinates method codes) that spend the same amount of computational time for each correlated-*k* permutation. This can be particularly time consuming in bands that have multiple, overlapping gases each with large numbers of correlated-*k* coefficients (e.g., water vapor and carbon dioxide overlap in the near infrared). The method for approximating an optimal distribution of the photons between bands for the broadband calculation of tropospheric heating rates is given in appendix A.

The result of these optimizations is that the model can compute on a workstation, for example, the broadband column absorptivity in approximately 1 min, or the height-dependent, spectral atmospheric absorption for a column in approximately 10 mins. This efficiency enables multiple sensitivity tests and ensemble statistics of multiple cloud fields. We have validated the model using the Discrete-Ordinates Radiative Transfer Model (DISORT; Stamnes et al. 1988) for broadband, plane-parallel calculations, and Trautmann et al. (1999) used SHDOM to validate 2D monochromatic calculations. In addition, we conducted a battery of studies to determine the number of photons required to achieve a desired level of model accuracy.

### b. Model inputs

The model profile uses 50 atmospheric layers with vertical resolutions of 0.5 km from the surface to 16 km, 2 km from 16 to 40 km, and 10 km from 40 to 100 km. Average profiles of pressure, temperature, and water vapor concentration were obtained from clear-sky soundings taken during the Central Equatorial Pacific Experiment (CEPEX; Williams 1993), and the ozone profile was obtained from CEPEX ozone sondes. Profile information above the maximum sonde heights were taken from the McClatchey model tropical atmosphere. Since we are focused on studying cloud effects, no aerosols are included in the profile. The surface is treated as being open ocean, with its bidirectional reflectance computed according to Podgorny et al. (2000). The top-of-atmosphere solar flux is from Kurucz (1992). For illustrative purposes, a solar zenith angle of 60° is used for most calculations, and sensitivity studies are computed for 0°.

The model uses a representative, although simplified, treatment of liquid water and ice clouds. Clouds are treated as pure ice for profile temperatures less than 245 K, as pure water for temperatures greater than 273 K; within the intervening mixed phase region, the ice and liquid water concentrations are assumed to vary linearly with temperature. The cloud scattering properties (spectral variation in extinction coefficient, single-scattering albedo, and asymmetry parameter) are computed using Mie scattering for the water clouds, and the Fu (1996) scattering parameterization for ice clouds. Mie calculations use an observed drop-size distribution for fair weather cumulus (Telford et al. 1984) with an effective radius of 9 *μ*m, which is consistent with those reported for tropical fair weather cumuli (Stephens and Platt 1987). The cirrus-scattering parameterization uses an average effective ice crystal size (*D*_{ge}) of 65 *μ*m, based on CEPEX aircraft data (McFarquhar and Heymsfield 1996). For simplicity, the temperature dependence of the crystal sizes were ignored here (sensitivity studies indicate this would slightly underestimate the albedo calculation).

The cirrus and water cloud visible extinction coefficients (at 0.5 *μ*m) are assigned separate, constant values for simplicity. The cirrus visible extinction coefficient is 1 km^{−1}. This value is consistent with the CEPEX lidar measurements (for cirrus with optical depths less than four), and in approximately the median of the values reported by Martin and Platt (1997). The water cloud visible extinction coefficient is 60 km^{−1}, consistent with cumulus observations (Bradley 1981). Although these extinction coefficient values are within the ranges observed, we recognize that they do not consider the possible variation between and within clouds. We know of no definitive way to treat the variation in extinction coefficient, particularly in the unsampled cores of deep convective systems. Thus, we chose these two values, rather than a more elaborate but as uncertain set of values, since it is easier to interpret the results. The effects of using these assumed values is assessed by a sensitivity experiment, and we note that this treatment will not effect the qualitative points of the results that we emphasize.

### c. Cloud fields

The Monte Carlo calculations require the horizontal distribution of the clear and cloudy pixels within the region, and for each cloudy pixel we need the cloud optical depth, cloud-top altitude, and the geometrical thickness. These properties are retrieved over the tropical western Pacific (TWP) from the Japanese geostationary satellite GMS-4 data, which were collected between 7 March and 7 April 1993 during the field phase of CEPEX. These properties are retrieved at a 5-km pixel resolution using the procedure described in appendix B.

We note that retrieving these properties from satellite data has the benefit of enabling us to study cloud systems and their variability as they occur, but constrains us to treating single-layer clouds within each pixel. Since heating rate profiles depend upon the vertical distribution of the absorbing matter, this constraint may not be optimal. However, as will be shown, the heating is generally concentrated in the top cloud layers encountered by the radiation, which are likely detected by the satellite. The unavoidable exception is the case of thin clouds (e.g., cirrus) overlying lower clouds.

Five diverse cloud scenes were selected to contain cloud types typically found in the TWP, including tropical convective systems, altocumulus, midlevel stratus systems, and mixtures thereof. Figure 1 shows the cloud-top height and optical depth fields. Note that these cloud properties were retrieved when the cosine of the solar zenith angle was greater than 0.6, which avoids the retrieval problems associated with low sun angles (Loeb and Davies 1996). The four inner grids in each scene cover a (400 km)^{2} area, which is slightly greater than the T42 grid [∼(310 km)^{2}] commonly used by GCMs for climate calculations. (The “T” notation is commonly used within the GCM community to define the grid resolution. It refers to the number of spectral wavenumbers kept in the triangular truncation; the larger the number, the finer the grid). All figures illustrate the large degree of subgrid-scale variability characteristic of these cloud systems, which cannot be treated explicitly in climate models. Table 1 summarizes descriptive statistics of the retrieved cloud properties given in Fig. 1.

### d. Computations

To efficiently compute the heating rates on a 3D, pixel-by-pixel basis (instead of only grid means) the Monte Carlo code was run on a Cray T3E Massively Parallel Processing Computer. This approach enables the computation and pixel-scale analysis of multiple cloud scenes for a variety of observed conditions. Over one billion photons are used for each scene. Each calculation takes approximately 250 h of processor time, which is typically distributed over 40 processors and requires less than 7 wall-clock hours to execute. Calculations use a cyclic (“wrap around”) boundary condition. Photons are injected into pixels in the inner domain and the surrounding border (outer domain) indicated in Fig. 1. Calculations in the outer domain only serve as boundary conditions for the inner domain; they are otherwise disregarded since they are contaminated by photons exiting the opposite side of the scene via the wrap around condition.

## 3. Results

### a. Anatomy of a 3D heating rate field

We first conduct a brief survey of the subgrid-scale structure of the 3D heating rate fields, with special attention paid to those features that ultimately have a significant impact on the GCM grid-mean heating rates. These properties are introduced for the cirrus and stratus decks that comprise scene 5. While our emphasis is on their effect on the grid-mean value, it is also possible that variations in the subgrid-scale heating rates of the magnitudes shown may affect subgrid-scale dynamical and cloud processes.

For reference, Fig. 2a gives the grid-mean heating rate for the (400 km)^{2} region of scene 5, which is computed by horizontally averaging the results of the 3D Monte Carlo calculation. The primary maxima at 9 km results from the rapid enhancement in the extinction coefficient moving downward from the cirrus cloud into the mixed phase region, and the secondary maxima at 5 km occurs at a common cloud height in this region.

Figure 2b illustrates one plane of a vertical heating rate cross section taken through the 3D grid, and Fig. 3 illustrates horizontal cross sections for two different atmospheric levels. Naturally, a multitude of finer features exist in the 3D fields that are not seen directly in the mean. The regions of maximum heating closely follow the boundaries of the cloud. The dark regions in the center of the cloud areas (Figs. 3b,c) occur where there is little or no radiative convergence (i.e., weak absorption). This caused by the attenuation of radiation by clouds in the overlying layers. The vertical gradient of water vapor absorption causes the broad background heating rate feature above the clouds, and produces the light streaks in the lower atmosphere where the direct beam penetrates through cloud gaps and is absorbed.

Cloud topography produces localized patches of enhanced heating rates that reach almost 20 K day^{−1} and exceed the layer mean by more than 15 K day^{−1}. These regions of large, local heating values or “hot spots” are largely governed by the variability in the cloud-top height, and by the structure of the cloud edge. (Note that these hot spots are not related to the same term referring to the region of solar retro-reflection by plant canopies.) The hot spots tend to occur: at the tops of clouds, where there is little above cloud attenuation; on the sun-side of cloud edges, for similar reasons; and where the cloud extinction coefficients are greatest, for example, the tops of mixed phase cloud region that begins ∼9 km or below ∼5 km where water clouds reside.

### b. Scale dependence of heating rates

We now illustrate the nature of the scale dependence of the 3D heating rate field. For this illustrative example (Fig. 4), we choose a threshold heating rate of 6 K day^{−1} because it is approximately the median value for the middle-tropospheric, 3D heating rates (not shown). The heating rates for each atmospheric layer are screened for pixels with values greater than this threshold. Such pixels that are connected, but separated from other such pixels by a gap, are defined as a single heating rate“feature” (Fig. 4b). For added scale-dependent information, we repeat the analysis for a second threshold value, which we take as 50% greater than the first or 9 K day^{−1}. We note that the choice of thresholds is somewhat arbitrary, but the general features that will be stressed hold regardless of the threshold used.

The scale analyses is applied to each scene to determine the contribution of the features to the total scene area, to the total heating rate, and to the scale dependence of the heating rate features. The results for scene 5 (Fig. 5a) indicate that heating rate features greater than 6 K day^{−1} occupy a maximum of almost 25% of the total scene area, but contributes as much as 60% of the total heating. The same general pattern holds for the 9 K day^{−1} threshold (Fig. 5b), where the features’ maximum contribution to the area are almost 15%, but defines almost 40% of the total heating rate for the scene. Thus, the small area occupied by the hot spots dominates the grid-mean heating rate.

The scales of these heating rate features in scene 5 are considerably smaller than the 10^{5} km^{2} scales (T42) commonly used in GCMs for climate simulations. For example, the total area greater than 6 K day^{−1} is a maximum of ∼10^{4} km^{2}, which is less by an order of magnitude (Fig. 5c). This total area is on par with the 4 × 10^{3} km^{2} scales (T213) typical of GCMs used for weather forecasting. However, even this smaller grid is not able to resolve the smaller features that comprise this area (e.g., the 25% contour), which are on the order of 10^{2} km^{2}–10^{3} km^{2} [approximately (20 km)^{2} to (30 km)^{2} grid]. Yet even a smaller grid may be needed to explicitly represent these features, because Fig. 4b indicates that they are not square and are often crescent shaped (i.e., one must resolve their minimum length or width). The 9 K day^{−1} features are even smaller, with the scales of the total area generally being between 10^{2} km^{2} and 10^{3} km^{2} (Fig. 5d). Thus, even the high resolutions used by GCMs for weather forecasting cannot explicitly represent these features, suggesting that a parameterization would be needed to include them in the grid mean (this is explored further in the next section).

The scale analyses for the other scenes reveal the same basic patterns shown for scene 5; an example of the variations are illustrated using scene 2 (Fig. 6). The bimodality in the vertical structure of the total area and total heating rate plots (Figs. 5a,b and 6a,b) are closely correlated with the sharp gradients in the layer average cloud extinction coefficients, which occur at cloud top and at the onset of the mixed phase region (Fig. 7). Overall for scenes 2–5, heating rates greater than 6 K day^{−1} constitute only about 5%–25% of the scene area but contribute between 30% and 60% of the total heating, respectively. Most of their total area contributions fall about or below the T213 resolution. (Scene 1 exhibits similar patterns, but its cloud fractional coverage was too small to produce significant heating rate perturbations.)

Sensitivity studies were conducted on scene 5 to examine the effects of using a different solar zenith angle, and altering the assumed cloud extinction coefficient profile.

Using a solar zenith angle of 0° (not shown) reveals the same patterns discussed for 60°; however, because of the larger amount of incident solar radiation, comparable features occur at larger thresholds. For example, results for a 18 K day

^{−1}threshold at 0° are similar to the 9 K day^{−1}threshold at 60°; both occupy about 12%–15% of the total area and contribute almost 40% of the total heating.The cirrus and water cloud extinction coefficients are increased by factors of 4 and 2, respectively. The resulting average cloud thickness is 28% of that for the standard cloud profile. Results for a solar zenith angle of 60° (not shown) show the same general patterns given in Fig. 5, but the heating rates are generally larger and the heating is shifted higher in the troposphere. For example, heating rates greater than 6 K day

^{−1}generally constitute only 10%–15% of the area, but contribute up to 60%–70% of the scene heating. This occurs because the greater cirrus opacity prevents radiation from penetrating deeper into the cloud and spreading the heating over larger region. Thus, the standard cloud extinction profile used may underestimate the effects of the hot spots in the upper layers of denser clouds.

### c. Minimum GCM scales for stimulating grid-mean heating rates

We have shown that a large portion of the solar radiative heating rates are attributable to horizontal scales that are smaller than those explicitly treated in typical GCMs. The objective of this section is to determine the optimal horizontal resolution needed by GCMs to accurately simulate the vertical distribution of the grid-mean heating rate while using their method of calculation. To investigate this issue, we systematically degrade the high resolution (pixel-scale) cloud inputs within a specified grid and compute the grid-mean heating rate.

We use the internal 64 × 64 pixels of each of the five scenes in section 2c, which yields a total grid area of (320 km)^{2} (roughly a T42 grid commonly used in GCMs). The total GCM grid is divided into subgrids of progressively smaller sizes until the resolution of an individual pixel is reached: 64 × 64 [(320 km)^{2} resolution, i.e., GCM equivalence], 16 × 16 [(80 km)^{2}], 4 × 4 [(20 km)^{2}],1 × 1 pixels [(5 km)^{2}, the reference case]. For each subgrid, the fractional cloud coverages per model layer are determined from the satellite retrieved cloud fields. The single-layer cloud retrieved at each pixel can lead to fractional cloud cover at different layers within each subgrid; this is because the pixel-scale, single-layer clouds reside at different cloud-top heights and span multiple model layers.

One fundamental issue we have to deal with concerns how to overlap the clouds between different layers. Since our intent is to test the accuracy of treatments in climates models, we adopt an approach similar to that used by GCMs that is designed to approximate the overlap in tropical convective systems: the convective clouds are assumed to be completely overlapping in the vertical, while other clouds are assumed to be randomly overlapped (conserving the cloud fraction per layer). This overlap treatment is complementary to the types of overlap combinations that Barker et al. (1999) studied for a fixed horizontal resolution. Following the cloud classification definitions used by the International Satellite Cloud Climatology Program (ISCCP; Rossow et al. 1996), we define convective clouds as the ISCCP“cumulus” or “deep convection” clouds with cloud top pressures >680 mb and ⩽440 mb, respectively, and optical depths (at 0.5*μ*m) <3.55 and >22.63.

Results for scenes 2 and 4 illustrate the systematic bias found in all scenes (Fig. 8). (The cloud fractional coverage in scene 1 was too small to produce heating rate perturbations large enough to warrant study.) Relative to the reference calculation, the (320 km)^{2} GCM grid overestimates the heating in the upper part of the cloud and underestimates it below. The reasons for this systematic bias are as follows. The GCM assumption of random overlap has the effect of filling the top layers with clouds such that photons do not penetrate deeper into the atmosphere. As a result, more of the solar radiation is absorbed in the upper layers with a positive bias and a compensating negative bias below. For example, consider a scene with 50% cloud cover that is 3-km thick and is distributed over three 1 km-thick layers. When the GCM uses the random overlap of clouds between the three layers, the resulting cloud cover would be 1 − (1 − 0.5)^{3} ≈ 87%, instead of the true cloud fraction of 50%. The second reason for the bias is that the averaging smears out the variation in cloud altitude (Fig. 2), removing the hot spots among other features. The height where the differences crossover from positive to negative in Fig. 8 coincides with the rapid increase in cloud extinction coefficient associated with the downward transition into the mixed phase zone (see Fig. 7).

The differences in Fig. 8 generally increase systematically with lower resolution. The agreement tends to be better for cases with the lower cloud amounts, or when the scene contains a large convective cloud fraction (in which the maximum overlap condition maintains a vertical cloud structure that is similar to the input data). For all scenes, the (20 km)^{2} case agrees best with the reference calculation, generally to within ±20% for all heights and to within ±10% for most heights. This resolution is in general agreement with the heating rate scale found in the previous section. Sensitivity studies were performed for a solar zenith angle of 0° (not shown). They produce grid-mean heating rates that are about double those shown for 60°, but the percentage differences from the reference case are very similar.

### d. Implications to GCM dynamics

We further investigate the differences in layer heating within these simulations to determine their potential effects on the model dynamics. For simplicity, we will use potential vorticity, a widely used quantity in atmospheric dynamics, in the following discussion. The time-rate-change in the quasigeostrophic potential vorticity, *q,* is related to the layer heating rate, ∂*Q*/∂*t,* by

where *f*_{o} is the Coriolis parameter, *σ* is the static stability, *R* is the gas constant, *C*_{p} is the specific heat at constant pressure, and *p* is the pressure. For simplicity, we define a term *G*_{q} as

Thus, if *G*_{q} is positive, the heating in the layer generates potential vorticity, and if *G*_{q} is negative, the heating destroys potential vorticity. The results, given in Fig. 9, indicate a consistent mean pattern of negative *G*_{q} aloft, and positive *G*_{q} below. Thus, the localized layer heating has a tendency to generate potential vorticity in the lower part of the cloud layer, and destroy it aloft. Note that the change in vertical motion is 90° out of phase with the rate of potential vorticity generation, so a decrease of potential vorticity tendency with height tends to increase upward motion.

Again, the differences from the reference case indicate the best agreement for the (20 km)^{2} case. The difference shown for the (320 km)^{2} case can be as large as a third more to double the reference value. Relative to the reference case, the vertical gradients of *G*_{q} are steepest (negative) around 9 km, suggesting that the GCM may overestimate the vertical velocity in this cloud region, and underestimate it below. This systematic pattern may have implications to cloud evolution, particularly for clouds in weak dynamical regimes. For overhead sun (not shown), the differences are even larger, approximately doubled. A discussion of the potential implications and consequences of these findings are presented in the next section.

## 4. Discussion and conclusions

We examine the scale dependence of solar radiative heating rates by combining cloud fields retrieved from satellite data with a broadband, 3D Monte Carlo model. Monte Carlo models can efficiently compute GCM grid-mean radiative properties (e.g., average broadband heating rate); however, problems that require the 3D subgrid-scale structure (pixel-level detail) from such computations are much more demanding. In this study, we alleviate this demand by using the computational power available through a massively parallel processing machine. This is the first time that 3D, pixel-scale analyses were performed for multiple cloud scenes.

We have identified within tropical-convective systems two key subgrid-scale properties that largely govern the grid-mean heating rates: the variability in the cloud-top height, and the structure of the cloud edge. These features produce hot spots—regions of intense local heating that can be several times greater than the grid mean. Although they occupy only a small fraction of the grid area, they dominate the grid-mean heating rate. For example, for the scenes considered here, heating rates greater than 6 K day^{−1} constitute only about 5%–25% of the grid area but contribute between 30% and 60% of the total heating rate, respectively. These effects increase with increasing cloud extinction coefficient, indicating the importance of accurately computing cloud properties in climate models and, particularly, the highly varying properties of cirrus.

The hot spots have scales of approximately (20 km)^{2}–(30 km)^{2}, which is an order of magnitude smaller than the grids typically used by GCMs for climate model simulations, and is smaller than that used for weather forecasting. Thus, these features would need to be treated by subgrid-scale parameterization. We show that a GCM-style calculation, employing a cloud overlap treatment similar to those currently used, requires data at a similar scale of (20 km)^{2} to accurately compute the grid-mean heating rate (to within ±10% for most heights, and ±20% for all heights). This suggests that current subgrid-scale cloud overlap treatments employed in GCMs may not enable accurate calculations of this type for grids much larger than (20 km)^{2}. Suggested improvements to such GCM subgrid-scale treatments (e.g., Oreopoulos and Barker 1999; Cairus et al. 2000) would require subgrid-scale information currently unavailable from the models. Thus, climate modeling and their associated observational programs will have to address the task of providing such subgrid information that can encompass these needed scales. While the emphasis of this paper has been on solar radiative heating rates, such scales might be suggestive of those needed for other subgrid cloud and cloud-radiative processes.

We find that GCMs tend to overestimate the solar heating within the upper cloud layer and underestimate it below. This differential layer heating might enhance the vertical velocity in the cloud layer and suppress it at cloud base, particularly for clouds in weak local dynamical regimes. If large enough, this would affect the vertical moisture transport and the associated cloud formation and/or maintenance. Thus, over the long-term, biases in the GCM treatments of the vertical heating may effect cloud evolution and feedback.

Of course, evaluating the potential impact of such radiative changes would need to be assessed in a GCM. The subgrid-scale parameterization used would need to be formulated in terms of variables computed by the GCM and, to be energetically complete, should include the effects of longwave radiation. Based on the disproportionate contribution of small-scale hot spots to the grid-mean heating rate, it seems that the parameterization would need to address the complicated task of accounting for the lower probability events that occur in the tails of the heating rate probability distribution.

## Acknowledgments

We thank E. Boer for discussions and the use of his satellite algorithms, including his BDM model. We also thank G. Zhang for his discussions and insights, and H. Barker for his comments on this manuscript. W. Collins provided information regarding satellite calibration, and P. Minnis provided the GMS-4 calibration coefficients used. W. Ridgway provided the original versions of the correlated-*k* generation algorithms used. We acknowledge the T3e support provided at the San Diego Super Computing Center and the Texas Advanced Computing Center. This work was funded by the Department of Energy (DOE-ARM Grant DE-FG03-91ER61198) and the NSF Science and Technology Center for Clouds, Chemistry and Climate (C4) Grant ATM-9405024.

## REFERENCES

### APPENDIX A

#### Optimal Sampling for a Broadband Monte Carlo Integration

We discuss an optimal sampling technique for broadband Monte Carlo integration based on the concept of stratified sampling (Press et al. 1990). Stratified sampling is a method of reducing the error in estimating the integral *F* = Σ^{M}_{j=1} *F*_{j} over a region (solar spectrum, in our case) based on the knowledge of the variance in the *M* individual subregions (spectral bands). For the case of broadband Monte Carlo integration, the variance of *F* is given by

where *N*_{j} is the number of photons used for the *j*th band and ɛ_{0j} is a coefficient.

To determine the optimal set of *N*_{j}, we first minimize *η* using computational time *T* = Σ^{M}_{j=1} Δ*t*_{j}*N*_{j} as a constraint. Here Δ*t*_{j} denotes average time required for simulating one photon trajectory in the *j*th band. The solution to the optimization problem is given by

Alternatively, the set of *N*_{j} can be found by minimizing *η* for a fixed number of photons *N,* assuming that all values of Δ*t*_{j} are the same. In this case, (A2) becomes

and is equivalent to the stratified sampling described in Press et al. (1990). Both modifications of stratified sampling produce nearly the same results for the set of spectral bands used in this study. We use a modification of (A3) that is computationally less expensive. We approximate ɛ_{0j} by

where *F*_{j} is (0, 1).

Barker et al. (1998) partition the photons among the spectral bands based on the fraction of the top-of-atmosphere solar irradiance contained within each band. However, we do not have a priori knowledge of how much each band contributes to the broadband heating rate. Thus, to optimize tropospheric heating rate calculations, we estimate the contribution per band by injecting the same number of trial photons into each band. For this process, only a rough estimate of the heating rate is needed, which is achieved by using only 1% of the total number of photons to be used in the broadband integration. An integrated heating term h_{j} is then computed for each band by vertically integrating the heating rate, with the values in the troposphere values weighted by one, and all other values weighted by 0.5. (This weighting factor could be refined further if, for example, the heating rate was desired in a region of the atmosphere more specific than the troposphere.) *F*_{j} is then the fractional contribution of *h*_{j} to its sum over all of the *j*-bands. The total number of photons used in the broadband integration are then distributed to according to the weighting computed for ɛ_{0j}.

### APPENDIX B

#### Cloud Property Retrievals from Satellite Data

This appendix describes how the model inputs of cloud properties are retrieved from satellite data. We use the Japanese geostationary satellite GMS-4 visible (0.5–0.75 *μ*m) and infrared (10.5–12.5 *μ*m) radiances. Their nadir resolutions are, respectively, 1.56 km^{2} and 25 km^{2}, which are transformed to an equal 0.044° lat–long grid with a 24 km^{2} resolution (∼5 km pixel). The radiances are converted to visible reflectances and brightness temperature, *T*_{B}, using the SeaSpace TerraScan software, and are calibrated using coefficients provided by P. Minnis (1999, personal communication). The four pieces of information that we need from the satellite data are: classification of each pixel as being clear or overcast (100% cloud cover), and the cloud optical depth, cloud-top altitude, and geometrical thickness for each overcast pixel.

*T*_{B} is the effective temperature at which the surface-atmosphere system radiates to space. Clouds radiate at their colder atmospheric temperatures and reduce *T*_{B} from its warm, tropical, clear-sky value. Based on the justification in Boer and Ramanathan (1997), we define pixels with *T*_{B} < 285 K as overcast; all other pixels are assumed in this study to be clear. Pixels may be partially cloud filled if their values are greater than 285 K and less than the value for clear skies, *T*_{clear}, estimated to be 292 K (determined from frequency histograms of *T*_{B} for the CEPEX region as being the warmest, common value). Due to the problem of separating the radiative contributions of the cloudy-sky and clear-sky portions in such partially cloudy pixels, their retrieved cloud properties are highly uncertain. Thus, to avoid contaminating our results with their uncertainty, we do not treat these pixels. Further, this *T*_{B} criterion also has the benefit of removing high, optically thin clouds for which the retrievals of cloud-top height are problematic for these satellite channels. We acknowledge that, by ignoring these pixels, we will underestimate the cloud cover. However, this will not effect the retrievals of the large, optically thick tropical-convective clouds systems that are the emphasis of this study.

The visible reflectances are converted to narrowband albedos using the observationally determined bidirectional models (BDMs) described in appendix A of Boer and Ramanathan (1997). The visible cloud optical depth (0.5–0.75 *μ*m), *τ*_{υ}, is determined from the albedo by lookup tables. (In principle, this procedure is essentially the same as converting radiances to optical depth directly by lookup tables.) The lookup tables tabulate *τ*_{υ} as a function of the visible albedo, solar zenith angle, and *T*_{B} (this follows from Boer and Ramanathan who base their cloud classification on *T*_{B} for their BDM constructions). The lookup tables are computed using the multiple-scattering radiative transfer code known as Discrete-Ordinates Radiative Transfer Model (DISORT;Stamnes et al. 1988). The calculations use the model atmospheric profile, ocean albedo, and the cloud models described in section 2b, with two modifications: a boundary layer sea salt aerosol was included with a 0.1 optical depth, and the cirrus cloud-scattering calculations use temperature-dependent crystal sizes, parameterized from the CEPEX aircraft measurements of McFarquhar and Heymsfield (1996). The maximum *τ*_{υ} in the lookup tables is 500, which is approximately where the albedo is constant with increasing *τ*_{υ} and retrievals of larger values are not reliable. The *τ*_{υ} determination from the DISORT lookup tables ignore any possible 3D effects that may be contained in the observationally determined BDMs.

The method for retrieving the cloud-top height, *Z*_{T}, is as follows. For a cloud whose optical depth is sufficiently large that its emissivity is approximately one, *T*_{B} is approximately the cloud-top temperature, and *Z*_{T} can be determined from the corresponding altitude in the temperature profile. In particular, if *Z*_{T} is 5 km, above which the water vapor absorption is small in the infrared window, *T*_{B} will be equal to the temperature at *Z*_{T} to within a few degrees. For smaller optical thicknesses, *T*_{B} includes the emission from the cloud and from the underlying atmosphere. In principle, we can use the infrared radiative transfer equation to solve for the cloud temperature from *T*_{B}, given the observed atmospheric profiles of temperature and water vapor concentration. In practice, however, we developed a simpler algorithm to relate *T*_{B} to cloud-top and cloud-base altitudes, and we tested it with the solution of the infrared radiative transfer equation. The starting point of the algorithm is the following equation,

where *T*_{cloud} is the effective brightness temperature for the cloud, and *B* is the Planck function at for a body radiating in the GMS-4 infrared channel at the indicated temperature. We estimate the cloud’s emissivity, ɛ, from ɛ = 1 − exp(−*τ*_{IR}/*μ*), where *μ* is the cosine of the satellite viewing angle, and *τ*_{IR} is the cloud’s infrared absorption optical depth taken as 0.49 of *τ*_{υ} [from cirrus single-scattering calculations based on Fu (1996) and Fu et al. (1998) for *D*_{ge} = 65 *μ*m, as in section 2b]. Thus, Eq. (B1) can be solved for *T*_{Cloud}, and the corresponding emission altitude is then determined from the observed CEPEX sounding. For optically thick clouds this emission altitude is approximately *Z*_{T}, and for thin clouds this is the altitude of the middle of the cloud. Specifically, we take *Z*_{T} for optically thick clouds (*τ*_{IR} > 2) as the emission altitude adjusted upward by a factor Δ*Z/τ*_{IR}, where Δ*Z* is the cloud’s geometrical thickness (discussed later); for *τ*_{IR} < 2, this factor is held fixed at Δ*Z*/2. We tested this method for retrieving *Z*_{T} using solutions of the infrared radiative transfer equation computed by the SBDART atmospheric radiative transfer model (Ricchiazzi et al. 1998). For the 28 000 overcast pixels involved in this study, *Z*_{T} is generally (about 98% of the cases) accurate to within 1.0 km, and almost always to within 1.3 km. We are not, however, claiming that this approximate procedure will work for cases not considered here and it should be used with caution, particularly for optically thin clouds.

The last parameter we need is the cloud-bottom altitude, *Z*_{B}. A cloud between *Z*_{T} and *Z*_{B} has a *τ*_{υ} given by

where *k*(*z*) is the cloud extinction coefficient; values used for the cirrus, water cloud, and mixed phase zone are described in section 2b. Based on the retrieved information and the *k*(*z*) assumed, Eq. (B2) can solved for *Z*_{B} and the geometrical thickness may be determined, Δ*Z* = (*Z*_{T} − *Z*_{B}). Because *Z*_{T} depends on Δ*Z,* and Δ*Z* depends on *Z*_{T}, their solutions are computed in three steps. First, *Z*_{T} is approximated as being the emission altitude and a trial Δ*Z* is found. Second, *Z*_{T} is recomputed using the trial Δ*Z* in the altitude adjustment factor discussed earlier. Third, the final Δ*Z* is recomputed for this *Z*_{T} (generally resulting in a minor revision in Δ*Z*). Obviously, *Z*_{B} and Δ*Z* will be affected by the uncertainties in the cloud extinction coefficients chosen for this study (discussed in section 2b). These effects are assessed by a sensitivity experiment.

As with any derived or modeled 3D cloud fields, we note that it is difficult to validate our retrieved cloud fields in an absolute sense. However, by limiting our retrieval to overcast pixels (with *T*_{B} < 285 K), we have reduced the possible errors in the fields. In addition, we note that the retrieved fields appear to have physically reasonable magnitudes (e.g., Table 1) and horizontal structure (e.g., Fig. 1). Also, the retrieved cloud-top heights are not inconsistent with CEPEX airborne lidar measurements of tropical–convective systems. Finally, the reconstructed cloud bases are relatively flat (since we reconstructed the cloud from the top down, sufficient errors in the cloud properties used could have produced unphysical bowl-shaped or umbrella-shaped cloud volumes). Although not a complete validation, these physically appealing constraints at least suggest that there are no dramatically unrealistic problems.

## Footnotes

*Corresponding author address:* Dr. A. M. Vogelmann, Center for Clouds, Chemistry and Climate, Scripps Institution of Oceanography, University of California, San Diego, Mail Code 0221, La Jolla, CA 92093-0221.

Email: avogelmann@ucsd.edu