To relate the error associated with 1D radiative calculations to the geometrical scales of cloud organization and/or in-cloud optical inhomogeneities, a new idealized methodology, based on a Fourier statistical technique, has been developed. Three-dimensional cloud fields with variability over a selected range of horizontal spatial scales and consistent vertical structure can be obtained and controlled by a small number of parameters, which relate directly to the dynamical and thermodynamical meteorology of the situation to be examined. This initial study deals with marine stratocumulus. Two experiments are conducted: an overcast situation and a broken cloud case with maximum cloud cover of 80%. For each experiment, five cloud fields are generated with the dominant organizational scale changing from 1.4 to 22 km, while all other quantities—such as cloud cover, cloud liquid water, and total water variance—remain constant. For each scene, three radiative calculations are performed for solar zenith angles of 0° and 60°: a plane parallel (PP) calculation, similar to that commonly implemented in general circulation models; an independent pixel approximation (IPA); and a full 3D calculation. The “PP bias” (IPA-PP) is used to assess cloud optical homogeneity approximation, while the “IPA bias” (3D-IPA) measures the impact of horizontal photon transport.
For the overcast scenes, the neglect of horizontal photon transport was found to be unimportant, and the IPA calculation gives accurate results. For the broken cloud case, this was only true for clouds with dominant horizontal spatial scales exceeding 10 km. With a scale of 2 km or less, the IPA bias in reflectivity, transmissivity, and absorption could exceed 5%. This indicates that even for shallow cloud systems, cloud geometry can play an important role. The sign of the bias depends critically on the solar angle, with IPA over- (under) estimating reflectivity for high (low) sun angles.
The PP bias in reflectivity was also found to be around 5% for both cases, comparable to the IPA bias, and smaller than previous estimates for this cloud type. Additional sensitivity tests prove this to be due to the vertical cloud structure. Vertically resolving the subcloud adiabatic liquid profile leads to a more opaque cloud upper boundary, reducing photon penetration into the cloud layer, and thus also PP biases. Additionally, for the broken cloud case it was found that changes in cloud fraction with height are translated by cloud overlap rules into effective horizontal variability in the liquid water path, further reducing biases. Taking a vertically uniform slab with identical integrated properties led to much larger PP biases comparable to previous estimates. Thus, models with sufficient vertical resolution are likely to suffer from smaller PP biases than previously estimated.
Three-dimensional (3D) solar radiative effects could be relevant for dynamical and climate issues; nevertheless, a full 3D radiative calculation is computationally unaffordable for many practical applications. When a 1D radiative calculation is performed, such as in the two-stream approximation used in most general circulation models (GCMs), two effects are neglected: the internal optical inhomogeneity of the clouds themselves, and the interaction between clouds. The radiative transport equation is solved for two independent horizontally homogeneous columns for the clear and cloudy regions. However, 3D radiative effects may need to be parameterized if proven to be important.
Many previous studies have concentrated their attention on in-cloud variability (e.g., Chambers et al. 1997a; Cahalan et al. 1994b; Marshak et al. 1998), examining the “plane parallel (PP) bias.” This is defined as the difference between a standard GCM-type calculation that neglects horizontal inhomogeneity,1 and the independent pixel approximation (IPA) method [sometimes referred to as the independent column approximation (ICA)], where the average of many PP column calculations for a much higher horizontal resolution grid is taken, in order to resolve the horizontal in-cloud variability. Cahalan et al. (1994b) and Barker and Davies (1992b) claimed that the nonlinearity of albedo with increasing liquid water content (LWC) results in substantial PP biases of around 15%. Cahalan et al. (1994b) even reported biases approaching 30% for very skewed liquid water distributions.
An analog to the PP bias can be defined, the “IPA bias,” that accounts for the inaccuracy of neglecting the spatial organization of the cloud inside the domain, and can be assessed from the differences between IPA and full 3D radiative transfer calculations. The obvious relationship between the radiation field and the photon pathlength would suggest that the importance of 3D processes could be conditioned by the different horizontal and vertical scales of organization of the clouds within the system in question. In broken cloud fields, horizontal radiative transport can induce enhancements of absorption and reflection between neighboring clouds (e.g., McKee and Cox 1974; Aida 1977; Welch and Wielicki 1984; O'Hirok and Gautier 1998a; Marshak et al. 1998; Barker et al. 1999). For example, Cahalan et al. (1994a) and Chambers et al. (1997b) found IPA albedo biases of 1% in overcast cloud systems, one order of magnitude smaller than the PP bias, while investigations of broken cloud systems (Welch and Wielicki 1984, 1989; Coakley and Kobayashi 1989; Breon 1992; O'Hirok and Gautier 1998a) found that the influence of cloud arrangement gave IPA biases of roughly 5%–15%, sometimes as much as 30%, due to side illumination, intercloud interaction and shadowing effects. These studies in general appear to show that broken cloud fields tend to decrease scene albedo with respect to the stratified homogeneous approximation, but with the aspect ratio of the clouds, the intercloud distances and sun position all play important roles such that even the sign of the IPA bias can be uncertain.
Part of the difficulty in assessing IPA biases arises from the method or source used to provide the cloud field investigated, which can be categorized as observations, numerical models, or idealized techniques. Using observations (e.g., Zuidema and Evans 1998) has the advantage that the cloud field generated is realistic in structure. However, the difficulties of retrieving cloud properties are not restricted to the condensate amount, but also to the dimensionality, with satellites or radar usually only providing a horizontal or vertical two-dimensional view. An additional problem of reconstructing cloud fields from satellite imagery such as Landsat data (e.g., Chambers et al. 1997b) is that 3D radiative effects smooth the retrieved cloud field variability at the smallest scales (Davis et al. 1997). The consequential application of a 3D solar transfer code to such a derived field is by construction obliged to underestimate the effect of 3D transport.
Like observations, numerical models are often also restricted to providing 2D cloud images (e.g., Fu et al. 2000), due to the computational expense of running both the cloud model and the 3D radiative code in three dimensions over large spatial scales. Some 3D studies of deep convective situations have been conducted (Barker et al. 1998, 1999; Tompkins and Di Giuseppe 2003) but for a limited set of cloud scenes.
A less obvious drawback of using both observations and complicated cloud models to provide the input field is that it is very difficult to conduct systematic sensitivity studies in which parameters such as the spatial organization can be changed in isolation to gain an understanding of the respective importance for 3D radiative transfer. This is where the strength of idealized studies lies.
Previous idealized investigations have taken a number of approaches. Welch and Wielicki (1984, 1989), Kite (1987), and Kobayashi (1993) have placed idealized cloud shapes such as cubes, with certain specified radiative properties, in a variety of highly geometrical arrangements. While instructive, it is not immediately clear if all the conclusions drawn from such highly idealized cases can be applied to natural systems. Additionally, the number of possible arrangements of an array of arbitrary cloud elements leads to the system having an indeterminant number of degrees of freedom, rendering practically unattainable the task of exploring the full parameter space.
Alternatively, several authors have adopted the fractal cascade model, which produces a linear log–log power spectrum of the liquid water cloud in the horizontal, usually for stratocumulus regimes in overcast conditions (Barker and Davies 1992a; Cahalan et al. 1994b; Davis et al. 1997; Marshak et al. 1998). A great deal has been gained from the use of such fractal models, but the cloud structures are disjunct from the meteorology that produces them in nature and it is difficult to introduce realistic vertical structure into such models. These investigations adopt a single pixel, vertically homogeneous cloud. It should also be noted that fractal models express the most variability on the largest length scales, equivalent to the domain length for which the radiation calculation is performed, since no long wavelength cutoff is applied. Introducing dominant scales of organization related to turbulent eddy or vertical cloud structure in idealized studies has remained a challenge.
Building on the experience of these previous idealized investigations, the aim here is to study the 3D radiative transport in cloud fields generated by a new idealized cloud model. The system for generating the artificial cloud fields is based on a Fourier transform technique such that cloud fields with variability on many spatial scales can be obtained and controlled by a limited number of parameters, which relate directly to the dynamical and thermodynamical meteorology of the situation to be examined. The model incorporates realistic vertical structure in the cloud fields using a simple model of a turbulent boundary layer. This flexible idealized model is used to conduct a systematic investigation into the effect of the spatial scale of organization of cloud fields on the IPA and PP biases.
After introducing the radiative model in section 2 of this paper, section 3 outlines in detail the new cloud generation technique. Section 4 summarizes the model parameters adopted for this study and presents the array of cloud fields to be studied. This initial investigation deals with a single-layer stratocumulus case, since there exists a wide bibliography for these clouds for comparison of the results and also to avoid the complications due to the mixed-phase and vertical cloud overlap in complex multilayer cloud scenes. The extension of the model to such multilayer clouds will be the subject of future research. Sections 5, 6, and 7 discuss the results of the full 3D radiative transfer investigation and compare the results to both IPA and PP calculations, and section 8 draws the conclusions of this study.
2. Radiative models and assumptions
a. Radiative assumptions
In this initial investigation, only the cloud mass varies in the horizontal; water vapor, temperature, ozone, and CO2 are functions of height. A multispectral band calculation is performed, using the k-distribution model of Fu and Liou (1992). Six bands cover the solar part of the spectrum (0.2–4 μm). Water droplet interaction is described by Mie theory and Rayleigh scattering is also included. A gamma function is used to describe the water droplet size distribution. A fixed drop concentration number of 65 cm−3 is assumed, implying that the effective radius takes on a range of values dependent on the cloud mass mixing ratio, which will be given later. Above the troposphere, seven additional atmospheric levels are placed between 20 and 100 km, which are interpolated using tropical standard profiles (McClatchey et al. 1972). The surface albedo is set to 0.07 (ocean).
b. Numerical aspects
The domain used for the calculations consists of 128 by 128 grid points in the horizontal using a resolution of 350 m. For the 3D transfer calculation, periodic boundary conditions are assumed. The vertical resolution is set to 50 m in the lowest 2 km of the atmosphere, but is coarser above, with a total of 50 levels describing the troposphere. Each radiative code makes modifications to this grid resolution to improve solution accuracy, which are described below. In all calculations, the accuracy is driven by the angular resolution and the base grid. This study adopts a medium angular resolution (Nμ = 12, Nϕ = 16, where Nμ and Nϕ define the number of the integration zenith and azimuth angles).
c. Radiative models
Most existing three-dimensional radiation models have been based on Monte Carlo or spherical harmonic approaches. Monte Carlo methods explicitly model the passage of individual photons through the medium. The advantage is that energy conservation is guaranteed and that error estimates are derivable. The drawback is that more complex statistics, such as radiances, 3D fluxes in domain subregions, or fluxes integrated in any direction are cumbersome and time consuming to obtain. This is where the strength of the spherical harmonics approach lies, which explicitly solves the radiative transfer equation and offers the possibility of simultaneous calculation of these complex diagnostics, but with an inferior achievable accuracy.
Considering these points, we draw on the strengths of both techniques in this study. For all domain average statistics where accuracy is important, the results of Monte Carlo integrations are shown. When more complex statistics are required, the spherical harmonics method is employed. For both the Monte Carlo and spherical harmonics calculations, care is taken to ensure that the assumptions concerning radiative properties are identical, and a comparison of the two methods reveals that the generic results and trends are very similar, if not exactly equal, for the reasons given above.
1) Monte Carlo code
The radiative transfer calculations have been performed using a modified version of the forward Monte Carlo algorithm (GRIMALDI). Details of the main code can be found in Scheirer and Macke (2001). For gaseous interaction, Chandrasekhar (1960) and Bucholtz (1995) give the Rayleigh phase function and the extinction cross-section analytic forms, respectively. In each band, 106 photons are traced. The accuracy in area-averaged quantities is better than 0.5%.
Monte Carlo techniques usually perform the photon transport calculation on a finer grid than the base grid to ensure solution accuracy. Here, the photon transport uses a grid that is 10 times finer with a nominal horizontal resolution of 35 m. The cloud properties are interpolated onto this finer mesh.
2) Spherical harmonics code
The spherical harmonics approach was implemented using the spherical harmonic discrete ordinate method (SHDOM; Evans 1998). At each grid point the radiative source function is represented in spherical harmonics, which is converted to discrete ordinates and integrated to determine the radiances. Most radiative quantities can be consequently calculated from the spherical harmonic representation of the source function.
Similar to the Monte Carlo calculation, the spherical harmonics approach adopts a higher-resolution grid for accuracy, particularly where strong gradients exist in the source function, for instance at cloud boundaries. Rather than implementing a fixed higher-resolution grid, SHDOM employs a more flexible approach, where only cells with source function gradients exceeding a specified threshold are split into a number of smaller cells, with properties derived by interpolation. This makes efficient use of memory and ensures that excessive gradients are absent in the final field. Thus, in addition to the base resolution, the cell splitting threshold is important. To handle the high variability of the cloud fields to be investigated, a relatively sensitive threshold that only allows a 10% variation in cloud properties between cells is used. Investigations with higher angular resolution and a cell splitting threshold of just 1% estimated the resolution-dependent error in fluxes at approximately 2% using this setup.
3. Idealized Cloud Field
We present a generalized method for the generation of cloud fields with specified horizontal spatial variability and realistic vertical structure, while still being defined in terms of a small number of parameters.
Fractional cloudiness over a certain horizontal area is a direct result of horizontal variability of total water (cloud plus vapor) and temperature (affecting the saturation vapor pressure qs). Observations and theoretical considerations indicate that variability in total water (qt) is more important than that of temperature (Price and Wood 2002) and therefore only total water variations are considered here. Including temperature variability in this model would be straightforward.
If the probability density function (PDF) of total water G(qt) is known, then the cloud fraction C and mean cloud water qc are given by
where qs is the saturation mixing ratio.
These equations assume that no supersaturation can exist, a good assumption for warm clouds, but not necessarily for the ice phase (e.g., Heymsfield and Miloshevich 1995). Previous statistical cloud schemes (e.g., LeTreut and Li 1991; Bony and Emanuel 2001; Tompkins 2002) have used such an equation to calculate C by specifying G(qt) and estimating its moments.
However, specifying G(qt) does not provide any information concerning the arrangement of the clouds within the GCM grid space, or indeed how the clouds overlap in the vertical. This information is required in order to make the radiation calculation. A method is therefore used in which an idealized function for the horizontal total water anomaly is instead defined in frequency space. The cloud positions and liquid water content are determined using the Fourier transform of this idealized function, giving the total water anomaly q′t in real space with fluctuations over a range of spatial scales, on a suitably defined grid. The cloud water mixing ratio at each grid point is then simply
where qt is the background mean total water profile, which like qs, will be a function of height. In this way, information concerning the spatial arrangement of the clouds can be summarized in a limited number of parameters.
The generation of cloud fields using this Spectral, Idealized, Thermodynamically Consistent Model (SITCOM) can thus be divided into four sequential tasks, namely, the specification of
the basic state of the atmosphere through the definition of the mean vertical profiles for total water and temperature (giving qs) appropriate for the cloud regime to be simulated,
an analytical function in the frequency space that defines the spatial scale and organization of the total water anomalies about the mean atmospheric state,
a vertical variance profile that will determine the magnitude of the total water perturbations and therefore the maximum optical depth of the clouds, and
how the anomalies overlap in the vertical.
The remainder of this section explains the approach used to fulfill each of these four tasks of SITCOM. In the next section the assigned values of the parameters described below for the situation in study here are summarized.
b. Mean profiles
A cloud-topped mixed layer is generated starting from the well-known atmospheric anatomy for these clouds (Lilly 1968). The stratocumulus develops in a radiatively active turbulent cloud layer over a sea surface under a strong subsidence inversion. Thus, an appropriate lower-atmospheric boundary condition is a 1.5-K temperature (θ) difference to the imposed sea surface temperature, and a relative humidity (RH) of 80%, giving the boundary layer total water qBLt. The well-mixed boundary layer, which continues to the lifting condensation level (LCL), is assumed to have a constant dry static energy (s = cpT + gz, where cp is the specific heat capacity of dry air) and constant qt equal to qBLt.
The mean cloud-base height is directly determined by the LCL (a strong function of the imposed surface RH), while the cloud-top height remains a free parameter to be specified. Its position is approximately determined by placing a temperature and total water inversion. If the boundary layer total water mixing ratio were to be maintained through this cloud layer to the inversion, (i.e., nonprecipitating pseudoadiabatic parcel ascent), the mean adiabatic liquid water content qadiac would simply be the difference between qBLt and the saturation mixing ratio at each height, qadiac = qBLt − qs. In reality, precipitation and entrainment processes would prevent this adiabatic value from being attained. Observations (e.g., Paluch and Lenschow 1991; Wood and Field 2000) indicate that this can be simplistically incorporated by reducing the mean adiabatic cloud content by the adiabatic fraction ξ, such that the mean total water between the LCL and the inversion is defined as
Above the inversion, the temperature follows a moist adiabatic profile of constant moist static energy (h = cpT + gz + Lqυ, where L is the latent heat of vaporization) while the total water mixing ratio exponentially decreases with a scale height of 1.5 km.
c. Frequency space
The next task is to define the total water mixing ratio anomaly [q′t(νx, νy)] function in a Fourier frequency (ν) space. This will be transformed to the real space on a square domain of horizontal dimension L with equal grid resolution in the x and y directions of Δx(=Δy). The frequency ranges between L−1, corresponding to one wavelength across the domain, and the Nyquist frequency 0.5Δ−1x. The Nyquist frequency essentially defines the highest frequency resolved by the selected resolution Δx.
The “synthetic” power spectrum is constructed with both the desired cloud spatial organization in mind, but also the mathematical constraints imposed by the Fourier analysis technique. Regarding the latter, in order to obtain a real function q′t(x, y) in the spatial space with a zero imaginary component, it is necessary to construct a complex Hermitian function in the frequency space, consisting of a real (even) and imaginary (odd) functions. The correct parity of the real and imaginary parts is ensured by incorporating the use of cosine and sine functions, respectively.
Since the aim of this study is to identify the effect that the cloud spatial scale of organization has on radiative properties, a function is desired that can isolate a specific range of frequencies, tailing off to zero power at long and short cutoff wavelengths. Mathematical tractability and simplicity are practical additional qualities. A gamma function fulfills these requirements. The use of a gamma function will not provide the power-law dependency often observed in stratocumulus liquid water path (LWP) power spectra (Cahalan and Joseph 1989; Davis et al. 1997). Due to its flexibility, it is straightforward to substitute an alternative function into SITCOM to render more realistic 3D cloud fields that reproduce these power-law features, and such a function will be used in a companion paper investigating sensitivity to other stratocumulus properties. However, the present aim of investigating sensitivity to cloud organization is not achievable with a wideband power-law function, hence the choice of the gamma function.
A random white noise function (constructed to preserve the correct respective parity) is added to break the symmetry of the real space fields. Taking all these features into account, the real and imaginary functions are therefore defined as follows:
where rν is the frequency radius
ℱℛ(rν) and ℱℐ(rν) are two random functions that add 10% variance preserving the parity of the Hermitian complex function, and σqt is the standard deviation of qt. The constants α and β determine the distribution shape and are discussed further below. The factor β(α+1)/2[Γ(α + 1)]−1/2 is a normalization coefficient. The resolution in frequency space is defined by the constant b, directly related to the Nyquist frequency, and can be defined in terms of the desired real space resolution as b = π/(2Δx).
From Parseval's theorem (e.g., Gabel and Roberts 1987, 263–265), imposing the conservation of energy between the spatial and frequency spaces, the power spectrum integral of the full frequency domain represents the total variance associated with the qt(x, y) field. The power spectrum is also called the covariance spectrum and expresses the variance of qt(x, y) at each particular frequency. Thus, the presence of σqt and the normalization factor in Eqs. (5) and (6) above becomes clear when the power spectrum function 𝒫𝒲(rν) is constructed, neglecting the random functions ℱℛ(rν) and ℱℐ(rν):
Thus, it is seen that 𝒫𝒲 is the product of a normalized gamma function and the variance of qt, the vertical profile of which shall be defined in the following section.
The length scale of organization of q′t(x, y) is determined by the configuration of the parameters α and β in Eq. (8), which set the position of the frequency peak (Fν) and the slope of the curve in the frequency space. The peak position of the 𝒫𝒲 function is given by the ratio Fν = α/β, and it expresses the dominant spatial scale over which the maximum power of the total water anomalies is concentrated. Thus, Fν can be considered the parameter that synthetically defines the scale of cloud organization for a generated cloud field.
For a given Fν, selecting small values of α and β renders a sharp power spectrum, with one dominant frequency, while large values of α and β render more equity across a wide range of spatial scales. A useful metric of this is the variance of the power spectra function. For the gamma function expression used in Eq. (8), this is σ2𝒫𝒲 = αβ2.
The previous two sections have described how the mean profiles are defined for T and qt, and have given the idealized frequency function for q′t used for this particular investigation. The power spectrum in Eq. (8) is scaled by the variance, which is defined in this section. The approach is to solve a diagnostic form of the turbulent variance equation for total water.
The first term represents the production of variance by qt fluxes in the presence of a vertical qt gradient, the second is a transport term, and the third term represents dissipation. Turbulence is the unique mechanism considered for these terms and is represented by simple K-diffusion theory, thus
where K is the eddy diffusivity factor.
The dissipation term ε is treated as a Newtonian relaxation to isotropy, as is usual, with a timescale τ:
This timescale is the ratio of a turbulent eddy velocity scale w (assumed to be 1 m s−1) and H, a turbulent height scale, which is taken to equal the inversion height; thus, τ = w/H. The equation for the variance consequently reduces to
The solution of the equation with zero boundary conditions at the surface and Z = ∞ is given by
The advantage of this approach, rather than imposing a typically observed variance profile, is that the variance is derived entirely from the imposed mean profile for qt and is therefore also completely consistent with the latter. Maximum variance values will be located at the inversion, where the vertical gradient of qt is greatest, and variance will also be nonzero in the homogeneous boundary layer due to the inclusion of the transport term. Within the limits imposed by the simplicity of the approach (e.g., the necessity of providing a constant value for the eddy diffusivity), no extra zero-order degrees of freedom are introduced.
e. Vertical overlap
The final task is to specify the vertical overlap of the clouds. Since we are dealing here with a single stratocumulus layer, the simplest assumption of maximum overlap is adopted. This premise, while valid for the stratocumulus regime, cannot be considered widely applicable for other cloud types, and will be modified in future studies with this model for multilayer clouds. To implement maximum overlap the same random functions ℱℛ(rν) and ℱℐ(rν) are applied at each height in the vertical. In other words, only one transform of the idealized frequency function is required, which is simply rescaled by σ2qt at each height level according to Eq. (8).
4. Experimental setup
a. Cloud fields
Two sets of cloud fields are generated by SITCOM, for which 3D, IPA, and PP radiative fluxes are then calculated. These illustrate the power of the approach that permits the selection of the scale of cloud system organization in isolation, while other system characteristics such as cloud cover remain invariant.
Two stratocumulus-type regimes, indicated as experiments A and B, have been generated from the same type of atmospheric profile. The parameters that describe the mean atmospheric profiles and the variance equation are given in Table 1, along with a summary of the domain size and resolution used. The inversion height is set at 0.8 km, with a temperature inversion of 4 K and a total water inversion of 2 g kg−1, as seen in the tephigram showing the mean humidity and temperature profiles (Fig. 1).
The only difference between experiments A and B is the setting of the adiabatic fraction ξ, which equals 0.8 in experiment A, resulting in a maximum cloud cover of 100%. For experiment B, ξ is given a much lower value of 0.2, for which the maximum cloud cover is roughly 80%. Using the basic parameter values given in Table 1, each experiment consists of five scenes that progressively change the scale of cloud organization, by varying the frequency at which the power spectral function peaks (Fν), which we recall is controlled by the ratio α/β. In order to retain the spread in frequency, σ2𝒫𝒲(=αβ2) is kept constant. The value for σ2𝒫𝒲 is taken from cloud-resolving model investigations (Tompkins 2001). In Table 2, a summary of the controlling parameters for 𝒫𝒲 as function of Fν for the whole set of experiments is reported.
Changing the wavenumber for the maximum of the power q′t changes the scale of organization of the cloud field without actually altering the mean properties of the cloud, such as the variance profiles, the PDF for the LWP, or the cloud cover (to the limit of the stochastic fluctuations due to the random functions), as can been observed in Fig. 2. The peak in variance at the inversion level is clear. The mean LWP for the A and B experiments is 126 and 57 g m−2, respectively. A fixed drop concentration number of 65 cm−3 guarantees that the effective radii reff will range between 7 and 23 μm with a mean of 11.2 μm for experiment A and between 5 and 20 μm with a mean of 8.6 μm for experiment B. A nonuniform effective radius across the cloud agrees with observations of the water droplet size dependency on the LWP (Stephens and Platt 1987). Approximately, the mean optical thickness is 16 for the A experiment and 10 for the B experiment, calculated from 1.5 LWP(ρreff)−1 as in Stephens (1976), where ρ is the density.
Note that the maximum value of the cloud cover does not occur at the inversion height as one might expect, due to the qadiac being largest there, but occurs lower in the cloud layer for both experiments. The reason is that the variance is much larger at the inversion. In regimes where qt > qs, increasing variance reduces cloud cover, while the opposite is true when qt < qs, as is apparent in the simulations of Tompkins (2002).
Figures 3 and 4 show the LWP fields for each scene investigated for experiments A and B, and reveal the difference in terms of internal inhomogeneity and how the scale of organization is effectively modulated by the frequency function 𝒫𝒲. Note that there is overlap between the power spectra of the sequential cases, especially between cases 1 and 2, which explains their resemblance. To illustrate the vertical coherence of the clouds and the appearance of the eddies, Fig. 5 shows a transect through the cloud field of experiments A3 and B3. The function of the maximum overlap and the effect of dry entrainment at the inversion are obvious.
In summary the method has been shown to produce three-dimensional cloud fields effectively with the specification of a limited number of parameters, for which, importantly, reasonable values can be obtained without excessive ad hoc assumptions.
b. Radiation calculation
The fields generated for experiments A and B have been used to calculate fluxes using a full 3D radiative transfer calculation in addition to simplified calculations based on the IPA and PP approximations. To avoid uncertainty, the same radiative algorithm is used for each mode of calculation, with simplifications applied to mimic the IPA and PP methods.
The IPA calculation is obtained by constraining each single column of the 3D atmospheric field to have periodic boundary conditions. In this way the horizontal transport of radiation is eliminated while conserving the horizontal optical inhomogeneity. The PP calculation is performed by first horizontally averaging the cloud condensate mass at each height and rearranging the clouds according to the maximum overlap rule usually applied in GCMs with vertically adjacent cloud layers. Examining the cloud overlap at each height, the atmosphere is then divided up into the number of columns required to capture the vertical cloud structure. The radiative transfer is then performed separately for each of these columns as for the IPA mode. Since the experiments have seven vertical cloud layers, this results in the PP calculation having exactly 8 separate columns for case B (one additional clear-sky column is required) and less for case A where several cloud layers are overcast. The clear and cloudy fluxes are then averaged according to the cloud fraction overlap following F↑↓ = Σ wiF↑↓i, where wi are the overlap weights and F↑↓ the upwelling/downwelling fluxes for each column.
Although the PP method employed in this paper has recently been implemented in a GCM by Collins (2001), it is not identical to that performed by most other GCMs, since, after a GCM calculates how much flux is transferred from cloudy to clear regions (or vice versa) at the interface between two vertical layers, the flux entering each cloud/clear region is averaged to obtain just two columns. This involves an implicit and artificial horizontal radiative flux at each layer, rendering the GCM calculation less accurate than the PP calculation performed here. We maintain that the PP method employed in this paper permits a more accurate assessment of the PP bias, since the differences between the IPA and PP calculation can only result from the neglect of horizontal inhomogeneity within cloud, and not from inadequacies in the method of calculation itself. In any case, the differences are minimal for the cloud scene studied here, with maximal overlap and cloud cover monotonically changing with height, when tested with the GCM radiation code of Morcrette (1991).
In Fig. 6 the A and B fields reduced to the “PP mode” are shown. The horizontal axis represents the weight associated with the column. Note how no clear-sky column is present in the overcast experiment A, and that it appears that almost a two-column calculation is performed since the cloud fraction is close to unity in a number of the cloud layers. In experiment B the weights of the columns are more similar.
In summary, the IPA-PP differences demonstrate the influence of the optical inhomogeneity, while the 3D-IPA comparison renders information about the importance of the horizontal photon transport. Calculations are made for solar zenith angles of 0° and 60°, which for brevity will be referred to as the SZA0 and SZA60 cases, respectively.
a. Reflection, transmission, and absorption
The domain mean statistics are examined in terms of top of atmosphere (TOA) reflection (R), surface transmission (T), and total atmospheric absorption (A). For each quantity x the relative IPA and PP bias is calculated as 100(R3D − RIPA)/RIPA and 100(RIPA − RPP)/RPP, respectively.
1) Case A
Examining first case A (Fig. 7a) for the sun overhead, there is very good agreement between the 3D and IPA calculations for all quantities, and for all scales of organization, with the bias limited to less than 1%. This is also the case for low sun angles (Fig. 7b). Thus, it seems that for overcast clouds horizontal photon transport perhaps plays a minor role, in agreement with Cahalan et al. (1994a) and others.
The PP biases are larger than the IPA bias. The bias in reflectivity is 4% and 5% for the two sun angles, while the transmissivity bias reduces from 5% to 3% as the solar zenith angle increases. The absorption bias is bigger for SZA60. Note that the IPA and PP calculations give identical results for all scales of organization, since by construction they cannot appreciate the reorganization of the columns as the scale of organization changes, as noted by Cahalan et al. (1994a). The insensitivity of PP bias to solar angle is expected from simple energy rescaling arguments presented by Varnai and Davies (1999). The magnitude of the PP reflectivity bias at 4% and 5% for the two sun angles is substantially lower than many previous investigations have reported. The reasons for this will be investigated.
2) Case B
For the partially cloudy B scene (Figs. 8a,b), the PP bias results are similar to the overcast scene, restricted to less than 5% for both the reflectivity and transmissivity for both sun angles.
With gaps now present in the cloud the 3D behavior is markedly different from the overcast scene. Taking first the sun overhead case, horizontal photon transport produces a general decrease of the reflectivity and an increase in transmissivity. At the smallest cloud scale investigated here, the IPA biases in reflectivity and transmissivity are comparable to the PP biases. As the scale of organization increases, the IPA biases decrease asymptotically toward a much smaller (but nonzero) value, which is attained for the largest cloud scale of 22 km.
The decrease in reflectivity and increase in transmissivity is due to the “spilling” of radiation from cloudy to clear regions. Any photon incident on a cloud may be scattered into an adjacent clear (or optically thinner) region and can continue, relatively unimpeded, perhaps to the surface. The effect, which is not present in an IPA calculation, obviously becomes more efficient as the scale of cloud organization decreases, increasing cloud boundaries and the vicinity of clear-sky gaps.
With a lower sun angle, the IPA biases in reflectivity and transmissivity are slightly larger. It is striking that the sensitivity to cloud-scale organization is opposite to the sun overhead case. With the cloud cells having a dominant horizontal spatial scale of 1.4 km, taking into account 3D radiation transfer actually increases reflectivity and leads to a corresponding reduction in transmissivity. When the cloud scale increases, the biases change sign at a dominant scale of around 5 km and tend toward the values found in the sun overhead case.
The bias sign reversal is due to the cloud side illumination effect where the effective cloud cover increases when the sun is lower in the sky. As the solar zenith angle increases, the chance of a photon passing relatively “unimpeded” through cloud gaps decreases, an effect not appreciated in an IPA calculation. This geometrical effect is obviously dependent on the aspect ratio between the cloud depth and the horizontal size of the clear-sky gaps in addition to the sun angle. For the shallow cloud layer investigated here, the effect is already minimal when the cloud scale reaches 5 km for an SZA of 60°. As the cloud exceeds this horizontal scale the opposing effect of photon spilling dominates. For a SZA of 75° this 5-km threshold would likely approximately double. Barker et al. (1998), Zuidema and Evans (1998), and O'Hirok and Gautier (1998a) have previously discussed some aspects of these opposing spilling and cloud side illumination effects in the context of a variety of cloud types.
Note that for a cloud scale of 1.4 km the magnitude of the IPA bias actually exceeds that of the PP bias and would probably be even larger for lower sun angles and smaller horizontal cloud scales.
b. Heating rates
For case B, which produced larger IPA biases in the main radiative statistics, the heating rates are also examined (Fig. 9). For the sun overhead case, the heating rates at cloud top are substantial, at around 30 K day−1, enhanced by the choice of high vertical resolution within the cloud layer (Petch 1995). The heating rate differences amount to around 5 K day−1 for the PP bias, while they are restricted to a maximum of 1.8 K day−1 for the IPA bias. Since horizontal transport is disallowed in both the IPA and PP methods, any heating rate PP biases in the cloud layer will be compensated by a bias of opposite sign in the subcloud layer. This is not the case for the IPA bias, which can be of one sign throughout the troposphere. Thus, while the peak heating rate difference is larger for the PP bias, the column mean bias statistics are of a similar order of magnitude in each case.
For both solar angles, the IPA bias is scale dependent as expected, but it is worth highlighting that for low solar angles, the IPA bias in heating rates reduces almost to zero when the cloud organization is concentrated on the larger scales (Fig. 9b).
6. Discussion: The IPA bias
In a distinct contrast to the PP bias, the IPA bias was shown to be highly sensitive to the interaction between the scale of cloud organization and the sun's position. It was suggested that this is due to a combination of effects involving side illumination, shading, and the scattering from cloud to clear regions.
In Fig. 10, the horizontal flux divergence/convergence normalized by the vertical flux convergence as function of LWC is shown for cases A and B. There is a clear signal of flux divergence in the optically thick regions and convergence in the clear-sky regions, as expected. As outlined above, this is due to the scattering of radiation from the cloudy to the clear sky. This effect increases dramatically in importance as the scale of organization decreases, due to the increase in the number of cloud boundaries and also when clear-sky regions are present since this provides clear-sky windows for photon transmission. This explains the reduction of total absorption for smaller organizational scales since photons are more easily able to “escape” through scattering to neighboring optically thin or clear-sky regions. The lower part of Fig. 10 shows the same statistic for case B. The same signal as in experiment A is observed both for the sun overhead position and when the sun is low while the magnitude of the convergence/divergence increases, showing that the sun position is of zero-order importance with respect to cloud cover.
Regarding previous investigations of IPA bias, Barker et al. (1998) and Chambers et al. (1997a) appeared to reach the conclusion that internal inhomogeneity of clouds is the dominant effect, while cloud geometry was a secondary issue, when considering the prevalent cause of inaccuracy in a PP calculation. This would mean that the IPA approximation, while not derived for broken clouds, could always be applicable. However, in both these studies, the success of the IPA method could be due to peculiarities of the isolated cases examined. For example, in Barker et al. (1998) clouds with very extended horizontal dimensions may minimize the effect of side geometry and aspect ratio. In Chambers et al. (1997a), the variability of cloud over smaller scales could be less than reality since the cloud structure is derived from radiation fields, which suffer from radiative smoothing effects (Davis et al. 1997), saturation effects, and also have no internal vertical structure.
It is therefore possible that investigations using cloud fields that place the emphasis on the large-scale component of organization (possibly generated through cascade models with no long-wavelength cutoff, or from observations of limited horizontal resolution) may underestimate the gravity of the 3D radiative transfer effects. Other studies such as O'Hirok and Gautier (1998a,b) and Welch and Wielicki (1984) seem to show that the cloud geometry may play a more important role than optical inhomogeneities. On the other hand, some of these studies have exceptionally high optical thicknesses or use the diffusion approximation that could lead to an overestimation of the geometrical effects.
In summary, IPA has been previously tested and validated mainly for scenes that have most of the horizontal variations over large spatial scales. At smaller scales, also important for remote sensing, the independent pixel concept is not able to improve the radiative transfer solution. The issue of subgrid-scale inhomogeneity loses relevance at these scales where horizontal fluxes remove the correlation between the radiative and thermodynamic properties of a column. This will be also clear from Fig. 11, which, in addition to showing the albedo for PP and IPA columns, also displays a scatterplot for the reflection of each point in the full 3D calculation. Markedly, no significant correlation is found between the LWP and the optical properties of the column.
7. Discussion: The PP bias
a. Contrast to existing studies
The second aspect of this investigation possibly unexpected from existing literature is that the PP biases appear to be limited for transmission and reflection with respect to previous studies of this cloud type, always less than 5%. The seminal study of stratocumulus clouds by Cahalan et al. (1994b) documented substantial differences in the albedo between IPA and PP calculations due to optical inhomogeneity. For solar zenith angle of 60° and a mean LWP of 126 g m−2 (approximately as in case A here), a minimum PP relative bias of 10% was reported, while an LWP of 57 g m−2 (as in case B, but overcast) gave approximately 15% biases, depending on the variance of the LWP distribution. The question naturally arises as to why the biases documented here might be less than those previously reported.
The study of Cahalan et al. (1994b) neglects vertical structure, considering only horizontal inhomogeneities in a single-layer cloud slab. The PP bias is calculated as the difference between the reflection for a uniform slab (with mean X and zero variance) and the reflection for a distribution of LWP (with same mean X and variance σ2).
There are two reasons why this approach may increase the magnitude of the PP bias estimate. The first concerns the neglected vertical structure in cloud liquid water. For a given optical depth, a vertically well-resolved cloud that has an adiabatic liquid water profile will have a higher liquid content near the upper boundaries, relative to a vertically homogeneous slab cloud layer. This would reduce the number of photons penetrating deep into the cloud layer, and could thus also act to decrease PP biases with respect to previous estimates based on slab cloud models.
The second reason why PP biases could be smaller with high vertical resolution is more subtle, and relates to the cloud fraction. To illustrate this, Fig. 11 shows the reflection (we have chosen scene B1, SZA60 as an example, but this does not change the generality of the argument presented) for each single IPA column as a function of the LWP. The reflection at the TOA has a strongly nonlinear regime for LWP smaller than 50 g m−2 and saturates at a value of around 0.45 for thicker clouds that is dependent on the assumptions concerning the effective radius distribution, the asymmetry parameter value, and the single scattering albedo. The square symbol represent the calculations performed for each of the seven PP columns, with the symbol size scaled according to the logarithm of the column fraction (weight). The mean values for the PP and IPA calculation are also plotted showing the larger PP reflection. Comparing the PP and IPA curves, it is clear that the PP is resolving some of the horizontal inhomogeneity of the LWP. This is due to the cloud slab being divided vertically into several (in this case seven) layers. Even though the cloud is horizontally homogeneous in each vertical layer, the fact that each layer is associated with a different cloud cover implies that the PP effectively appreciates horizontal inhomogeneity in LWP through the application of the cloud overlap rules (cf. Fig. 6 above). To emphasize this, the right panel contrasts the LWP “PDF” for the PP and IPA calculations, revealing that the PP constitutes an under sampling of the IPA calculation. Moreover, the largest weight of 80% is constrained to be the thickest part of the clouds and therefore in the saturation regime, where albedo increases slowly. The other columns in the PP calculation, while only responsible for 20% of the domain, are situated in the strongly nonlinear part of the LWC curve, thus effectively reducing the PP bias.
One might argue that this ability of a GCM is lost if cloud fraction is identical throughout a number of vertical layers. However, unless the distribution is strongly skewed, the cloud is in this case likely to be homogeneous, or of substantial optical thickness such that that PP biases will be reduced to a minimum (due to the saturation regime). The reasoning is that the addition of any variability to an optically thin (i.e., in the strongly nonlinear reflectance regime) cloud layer will necessarily lead to partial subsaturation of the grid box and thus cloud fractions less than unity, consequently partially “resolvable” by a PP approach, provided sufficient resolution is used.
To summarize, cases A and B have relatively small, almost identical, PP reflectance biases. The first mechanism proposed for restraining the PP bias magnitude in both cases is simply that, by resolving the stratocumulus adiabatic vertical structure in liquid water, the cloud tops are more opaque, preventing deep photon penetration. However, it was additionally suggested that the bias could be limited for the overcast cloud since it is optically thicker and thus resides in the flat, saturated reflectance regime, while the PP bias of the optically thinner case B is small because the LWP variability is resolved partially by the PP approach.
c. Sensitivity tests
Since it is not obvious which, if any, of these two hypotheses dominate, two additional experiments have been conducted. First, we remove the ability to resolve the vertical liquid water structure by replacing each cloudy grid point with the vertical mean in-cloud value. In this way the cloud region becomes vertically homogeneous, while the cloud mask is retained.
The second test goes further by vertically averaging the liquid water for all clear-sky and cloudy points throughout the layer between cloud base and cloud top. This results in a cloud layer that is vertically homogeneous and with constant cloud fraction. The rearrangement of the clouds is summarized schematically in Fig. 12. Note that both experiments do not alter the PDF of LWP and that the radiation calculation is still conducted with a 50-m vertical resolution in all cases in order to preserve the numerical accuracy.
Figure 13 shows how the PP relative bias increases drastically for experiment B when the cloud is vertically homogenized but the cloud mask is retained. The increase is a doubling from 4% to 9%. When the cloud fraction structure is lost by the second averaging technique, rendering a slab cloud as in Cahalan et al. (1994b), the bias further increases to 13%. These figures are more in line with those of Cahalan et al. (1994b).
As expected, the change is more limited for the optically thicker experiment A, with the bias increasing to 6% as the liquid water is homogenized. Averaging the cloud mask has little additional impact, since nearly all cloud layers were overcast. This differential growth in PP bias between cases A and B is exactly as predicted. The sensitivity tests would indicate that the vertical subcloud variability is equally important as horizontal variability. Cloud fraction plays a secondary, but nethertheless nonnegligible role; the greater the number of vertical layers, the more horizontal LWP variability is resolvable by cloud fraction changes.
Cahalan et al. (1994b) also stated that skewed LWP distributions exacerbate the PP bias for a given variance. The case A overcast simulations reported here have a slightly asymmetrical PDF for LWP (skewness = 0.4), while the PDF for lower fractions is strongly positively skewed (skewness = 1), since only the tail of the total water distribution is reflected in the LWP. These characteristics of the LWP skewness versus cloud cover reproduce the observations of Wielicki and Parker (1994). However, the fact that the PP bias remains limited in case B despite the skewed distribution would indicate that the influence of the third moment of the distribution is relatively minor with respect to the vertical resolution. Again this is easy to appreciate since skewness will not produce a difference in the asymptotic saturated regime.
Cahalan et al. (1994b), in their seminal study of single-layer fractal stratocumulus clouds, stated that future research should include cloud models that were also capable of generating realistic vertical structure. In this paper we present a new model, the Spectral Idealized Thermodynamically Consistent Model, or SITCOM, which generates three-dimensional idealized cloud fields with horizontal structure over a defined range of spatial scales and a self-consistent vertical structure, with the definition of a very limited number of parameters.
The approach of the model is to define an idealized power spectrum in Fourier space, the inverse Fourier transform of which produces a two-dimensional anomaly field for the total water. The perturbations are about a specified mean state appropriate for the meteorological situation to be modeled. The anomaly field is scaled to give the correct variance obtained by the solution of a diagnostic turbulent variance equation. Comparing the final total water field to the saturation mixing ratio allows the cloud mass to be diagnosed. The maximum overlap assumption was adopted for the vertical correlation of cloud layers in this initial simplified single-layer cloud investigation, which will be generalized in future work.
The strength of the model is that properties such as the domain average cloud fraction and diagnosed liquid water content can be held constant, while other cloud field features are altered in isolation. This study chose to investigate the effect of the horizontal spatial scale and thus a power spectra function was selected that allowed a range of organizational scales to be isolated in turn. The model was set up to produce an idealized stratocumulus layer with two sets of five scenes with peak spatial scales ranging from 22 to just 1.4 km, in a domain size of 44 km by 44 km. The first set of experiments generated overcast conditions, while the second set had a mean cloud cover that peaked at 80%. Experiments were conducted for solar zenith angles of 60° and for overhead sun conditions. Finally, three methods of the radiative transfer calculation were employed consisting of a full 3D calculation, an independent pixel approximation (IPA) approach, and a plane parallel (PP) method. The difference between the IPA and PP calculation determined the influence of horizontal in-cloud inhomogeneity and was termed the PP bias, while the 3D-IPA contrast assesses the effect of horizontal photon transport and was referred to as the IPA bias. From the analysis of the two different biases as a function of the scale of organization, it was possible to quantify in which metric the mean optical properties of the atmosphere are affected by the cloud geometrical arrangement and/or their optical internal inhomogeneities.
Assuming that the full 3D calculation represents “truth,” in the overcast scene, the IPA approximation was found to give very accurate assessments of domain mean quantities of reflection, transmission, and absorption for all scales of cloud organization, agreeing with previous studies. This was also the case for the broken cloud scene with 80% cover, but only when the cloud horizontal variability was chiefly organized on scales exceeding 10 km. As the scale of organization decreased, so the IPA bias increased. With the sun overhead, the IPA approach overestimates the reflection, since photons are unable to scatter and escape through the adjacent clear-sky gaps. For low sun angles the standard IPA calculation underestimates the reflection, since it cannot appreciate the apparent increase in cloud cover. For cloud scales of less than 2 km, the bias rivals that from the neglect of horizontal subgrid-scale variability at around 5%, indicating that even for shallow stratocumulus type systems, the cloud geometry (neglected by IPA and PP approaches) can be just as important as internal optical variability.
This assessment of PP reflection bias at around 5% is, on the other hand, smaller than previous investigations of stratocumulus clouds. Earlier studies using fractal models focused on clouds occupying a single vertical layer. If a GCM employs a vertical resolution that resolves the cloud slab, two factors can reduce PP bias. First, resolving the adiabatic liquid water structure increases the opacity of the cloud upper boundary, reducing deep photon penetration of the cloud layer, and subsequently reducing PP bias. Additionally, variations in cloud fraction with height are effectively converted into horizontal variability of LWP by the application of cloud overlap rules. Sensitivity tests proved that both these mechanisms could be important, but that the structure of liquid water predominates and can exceed the bias due to the neglect of horizontal water variability.
Thus GCMs can, in theory, improve on previous estimates of their PP biases, provided adequate vertical resolution is used. The necessary vertical resolution required to reduce the bias obviously depends on the structure, location, and opacity of the cloud system, and it is not claimed that the 50-m resolution used in this study is generally appropriate. Moreover, current resources prevent present GCMs approaching such high vertical resolutions, even in the boundary layer (Tompkins and Emanuel 2000). The main point to emphasize is the importance of vertically resolving clouds and cloud processes and to show that if this is achieved, the method employed by GCMs is potentially capable of producing smaller biases with respect to an IPA calculation than previously estimated.
In summary, this study suggests that GCMs with sufficient vertical resolution could perform better than previously assessed with their PP calculations. On the other hand, if clouds are organized on small horizontal scales, biases due to the neglect of horizontal photon transport could exceed previous assessments if the scene is not overcast. These biases became significant as the dominant horizontal cloud scale fell below 5–10 km. These conclusions could only be revealed by the use of a controllable idealized model that can generate fully three-dimensional cloud scenes with horizontal, and perhaps more pertinently, consistent vertical structure.
That said, this study is only a first step, investigating just two realizations of one cloud type. Although the results are likely to be robust to changes in the mean temperature and humidity profiles, other parameters such as the cloud-top height are likely to have a greater influence. Moreover, multilayer cloud systems that occupy the entire tropospheric depth involve a still greater level of complexity that has to be addressed in order to generalize the conclusions drawn here.
We are grateful to Franklin Evans, Ronald Scheirer, and Andreas Macke for making their codes publicly available and giving advice on their use. The aid of Jeff Settle was also appreciated, as were the views and support of Keith Shine. An early version of this manuscript benefited from the critical attention of Anthony Slingo and three anonymous reviewers and Paquita Zuidema provided further helpful comments. The first author was supported by NERC Contract GT/4/00/217.
Corresponding author address: F. Di Giuseppe, Environmental Systems Science Centre, University of Reading, 3 Earley Gate, Reading RG6 6AL, United Kingdom. Email: firstname.lastname@example.org