1 Introduction

The ultimate boundary conditions of Earth’s climate system are defined by the flow of solar radiation in and longwave emission out of the top-of-atmosphere (TOA). All else constituting climate and life falls between these boundaries. Understandably then, radiative transfer through, and optical properties of, the Earth-atmosphere system lie at the core of our concepts of climatic forcing and feedbacks. Moreover, many remote sensing techniques, both passive and active, involve measurement and subsequent model interpretation of (largely atmospheric) radiative transfer. Mathematical treatments of radiative transfer were developed initially by astrophysicists (e.g., Chandrasekhar 1950), but the latter half of the 20th century saw substantial advancements from Earth scientists in their response to the myriad conditions and scales of Earth’s climate system. Arguably, these advancements are exemplified best by work centred on radiative transfer for cloudy atmospheres.

Over the past three decades there has been a string of diagnostic studies aimed at demonstrating errors associated with neglect of 3D radiative transfer within cloud, weather and climate models (e.g., McKee and Cox 1974; Davies 1978; Ellingson 1982; Welch and Wielicki 1984; Stephens 1988; Davis et al. 1990; Barker and Davies 1992; Cahalan et al. 1994; Marshak et al. 1995; O’Hirok and Gautier 1998; Barker et al. 1999, 2003; Cole et al. 2005a). Generally, these studies used small, and often selective, samples of 3D cloud fields derived either from idealized models, such as arrays of simple forms (Welch and Wielicki 1985; Kobayashi 1988) and scaling algorithms (Marshak and Davis 2005a), passive satellite data (Cahalan et al. 2005), aircraft data (Barker 1992; Räisänen et al. 2004) or limited-area domains simulated by cloud system-resolving models (Barker et al. 1999; Hinkelman et al. 2005). Ideally, a large number of 3D cloud fields based on observational data should be used to demonstrate errors due to neglect of 3D solar transfer. This would enable proper hypothesis tests to be carried out involving clouds whose geometric properties resemble most faithfully those that actually occur.

Cloud properties derived from synergistic retrievals using data from lidars, cloud-profiling radars and passive radiances provide, thus far, the most complete pictures of true cloud structure (ESA 2001; Stephens et al. 2002; Hogan et al. 2006; Winker et al. 2007; Dupont et al. 2011). However, because current active-passive systems point fixedly at either zenith (ground-based) or nadir (satellite-based), they provide 2D vertical cross-sections, not 3D domains. Those derived from satellite data have the advantage of providing a global perspective on a regular grid whose horizontal spacings are not subject to height-dependent advection rates. Their disadvantage is that they are less resolved than, and not as sensitive as, their surface counterparts.

To get around the issue of satellite retrievals being essentially 2D cross-sections, 3D fields can be constructed using the algorithm of Barker et al. (2011) that was developed for use by the EarthCARE satellite mission (ESA 2001). EarthCARE is slated for launch in 2015 and will carry a high-spectral resolution lidar, a 94-GHz cloud Doppler radar, a 7-channel imager and a broadband radiometer. Their algorithm combines conventional passive satellite imagery with nearby coeval information from active-passive profiles; recipient columns associated with off-nadir passive pixels get assigned, by proxy, an active-passive donor column associated with a nadir pixel thereby leading to production of a 3D domain of cloud.

This study represents a first step towards comparing 1D and 3D solutions for large samples of cloud fields derived from satellite data. A-train data and retrieval products (Kato et al. 2010) construction algorithm developed by Barker et al. (2011) to construct 3D cloudy atmospheres. A broadband solar transport solver (capable of affecting 1D and 3D solutions) was then applied to them. Reflected fluxes from the 3D model were compared to corresponding (20 km)2 CERES values (Wielicki et al. 1996). For (50 km)2 domains that resemble the size of weather and climate model grid-cells, flux differences arising from 1D and 3D solutions were compared.

In the following section the A-train dataset is described briefly. This is followed by descriptions of the 3D construction algorithm, the radiative transfer model and the nature of the simulations. Results are presented in Sect. 6, and conclusions and recommendations are in the final section.

2 Merged A-train Data

Data gathered by several A-train satellites were used to assess the importance of 3D solar radiative transfer. A merged A-train dataset has been compiled at NASA-Langley in which CloudSat (Stephens et al. 2002) and CALIPSO (Winker et al. 2007) profiles were mapped to 1-km resolution thereby associating each column with a 1-km MODIS pixel (Kato et al. 2010). It also includes the line of (20 km)2 CERES broadband radiances (and fluxes) along the CloudSat-CALIPSO track as well as MODIS pixels either side of it out to ∼40 km.

A merged cloud mask was created by interpolating CloudSat’s and CALIPSO’s cloud masks onto CERES’s arbitrary vertical grid which includes temperature, humidity and ozone (from the Goddard Earth Observing System—Data Assimilation System—GEOS-4; see Bloom et al. 2005; Kato et al. 2005). If a CloudSat or CALIPSO cloudy cell overlaps a CERES layer, the cell is designated to be filled with cloud.

Because a large fraction of columns have just cloud masks from CloudSat and CALIPSO, vertical profiles of cloud properties were set using MODIS quantities and the following algorithm. For each MODIS pixel along the CloudSat-CALIPSO track, effective visible optical depth τeff, effective particle size r eff and particle phase were retrieved based on inversion of a 1D radiative transfer model (Minnis et al. 2008, 2010a, b). If there are n ice and n liq cells in a column, the former having temperatures <273 K and the latter >273 K, it is assumed that

$$ \frac{3}{2}\left(\sum_{i=1}^{n_{\rm ice}}\frac{{\mathcal{I}}}{\rho r_{\rm ice}}\Updelta z_{i}+\sum_{i=1}^{n_{\rm liq}}\frac{{\mathcal{L}}}{r_{\rm liq}}\Updelta z_{i}\right) =\tau_{\rm eff} $$
(1)

where \(\mathcal{I}\) and \(\mathcal{L}\) are ice and liquid water contents, respectively, ρ ≈ 0.917 g cm−3, and \(\Updelta z_{i}\) is layer thickness. If CloudSat’s columnar classification (Sassen and Wang 2008) is cirrus, altostratus, altocumulus or deep convection, all ice cells in the column use r ice = r eff, and any liquid cells use r liq = 10 μm. When the columnar classification is stratus, stratocumulus, cumulus or nimbostratus, all liquid cells use r liq = r eff, and any ice cells use r ice = 50 μm (e.g., Yang et al. 2007). Equation 1 is then solved for water contents assuming that \(\mathcal{I}=\mathcal{L}\) but restricting \(\mathcal{I}\) to be no larger than 0.3 g m−3 (e.g., Protat et al. 2007).

This simple algorithm produces clouds that are, for the most part, locally homogeneous in the vertical. Obviously this can influence both 1D and 3D radiative transfer solutions, as well as their differences which are of concern here. Nevertheless, when clouds produced this way are operated on by 1D shortwave and longwave radiative transfer algorithms, local TOA flux estimates agree well with CERES data, particularly longwave values (Barker et al. 2011).

Without defending this vertical allocation scheme too much, it is worth noting that for geometrically thin clouds that get resolved into a small number of layers, such as many boundary layer clouds, the impact of vertical structure on 1D–3D transfer differences is a minor issue. The same goes for very cloudy scenes, which as will be seen constitute the majority of the cloudy (50 km2) domains. Finally, preliminary studies, not reported here, using data from cloud system-resolving models suggest that vertical homogenization of water for towering convective clouds generally exaggerates 1D–3D differences, implying that results presented below are approximately an upper-bound for 1D–3D differences.

3 Constructing 3D Cloud Domains from A-train Data

For the convenience of readers, Barker et al. (2011) construction algorithm is reiterated in this section. A-train data make-up the 2D retrieved cross-sections (RXS). Let all pixels/columns along the RXS at positions (i, 0), where for this study i = 1 was in the southern hemisphere of an ascending portion of an orbit, form the set of potential donor columns that can act as proxies for off-nadir recipient columns associated with pixels located at \(\left\{(i,\,j)\mid j\in \left[-J,-1\right] \cup \left[1,J\right] \right\}\) . For each recipient pixel, begin by computing

$$ F\left(i,\,j;\,m\right) =\sum_{k=1}^{K}\left[\frac{r_{k}\left(i,\,j\right) -r_{k}\left(m,0\right)}{r_{k}\left(i,\,j\right) }\right]^{2}{:}\, m\in \left[i-m_{1},i+m_{2}\right] $$
(2)

where r k are MODIS radiances, and k denotes spectral interval. For this study m 1 = m 2 = 200 pixels which corresponds to searching for ∼±200 km along the RXS beginning at the pixel on the RXS that is in the recipient’s cross-track line. Next, F(i, j; m) are ordered from smallest to largest to form the set

$$ \begin{aligned} \left\{\min \left[F\left(i,\,j;m\right) \right] ,\ldots ,\max \left[F\left(i,\,j;m\right) \right] \right\} & = \left\{{\mathcal{F}}\left(i,\,j;1\right) ,\ldots ,{\mathcal{F}}\left(i,\,j;m_{1}+m_{2}+1\right) \right\}{:} \\ & \qquad {\mathcal{F}}\left(i,\,j;n\right) \leqslant {\mathcal{F}} \left(i,\,j;n+1\right). \end{aligned} $$
(3)

Defining Euclidean distance between a potential donor at (m,0) and recipient at (i, j) as

$$ D\left(i,\,j;m\right) =\Updelta L\sqrt{\left(i-m\right)^{2}+j^{2}}, $$
(4)

where \(\Updelta L\) is imager resolution, let D(i, j;m) and m go passively along with the ordering of F(i, j; m) so that \(\left\{\mathcal{F}\left(i,\,j;n\right) \right\} \) has an associated co-ordered set of distances \(\left\{\mathcal{D}\left(i,j;n\right) \right\} . \) One then solves for the index m * as

$$ {\mathop{\arg \min}\limits_{{m^{\ast}\in \left[1,\left(m_{1}+m_{2}+1\right) f\right] }}}\left\{{\mathcal{D}}\left(i,\,j;m^{\ast}\right) \right\} : f\in \left(0,1\right) , $$
(5)

which means: find the index m * that corresponds to the smallest distance between recipient and those pixels that constitute the smallest 100f% of F(i, j; m) . For this study f = 0.03 so that the smallest 3% of the (m 1 + m 2 + 1) values of F(i, j; m) are considered (provided they have the same surface type as the recipient). Knowing m * leads directly to m and thus the donor column. Finally, all the properties associated with the column at (m, 0) get replicated at (i, j). In addition to cloud properties, this includes profiles of temperature, moisture and aerosol, as well as surface conditions. The vast majority of cross-sections used here were completely over ocean and so there were 133 layers from the surface up to 65 km.

Figure 1 shows a schematic of this process which gets applied until the desired 3D scene is constructed. Clearly, RXS columns (j = 0) identify themselves so for them there is no need to apply the algorithm; the RXS forms the centre of the constructed domain.

Fig. 1
figure 1

The thin RXS is shown along with the sequence of MODIS visible pixels associated with it. The objective is to fill the volume marked by the wider MODIS swath with cloud properties drawn from the RXS. For example, the column associated with the pixel at (m, 0) has been designated as the proxy for the pixel at (i, j) so the cloud-radiation attributes associated with (m, 0) get donated to (i, j). The algorithm is applied until all desired off-RXS pixels are filled by donor RXS columns

Full size image

Verification of the construction algorithm using A-train data is not straightforward. One can, however, go a certain distance by attempting to reconstruct the RXS itself. In so doing, one attempts to fill an RXS column at (i, 0) by searching the RXS and applying (5) to potential donor pixels in [i − m 1, i − n] ∪ [i + n, i + m 2] which bars the first ±n pixels next to i; hence defining a dead-zone in the search process. This test is meant to mimic filling off-RXS columns that are ±n pixels away from the RXS. For example, when n = 5, searching for a proxy column begins five pixels away, just as for off-RXS pixels at (i, ±5).

Figure 2 shows attempts to reconstruct a 400-km-long stretch of RXS. The upper image is the actual RXS merged cloud mask. This is an especially demanding case as it involves fairly dense multi-layer clouds. Over much of this domain passive-only retrievals would yield very little, if any, useful information about cloud vertical structure. Lower images are reconstructions for discrete values of n = 1, 5, 10 and 20. These results stem from using four spectral channels (0.62–0.67; 2.105–2.155; 8.4–8.7; and 11.77–12.27 μm) in (2). By n = 5, which corresponds to the outer edge of an 11-km-wide domain the likes of those to be computed for EarthCARE, it is clear that some error was creeping in; nevertheless, a significant amount of detail was captured. Even out at n = 20, multi-layers of clouds were replicated well. The region where the greatest difficulty was encountered between 100 and 200 km along the horizontal where clouds were transitioning or dying out entirely.

Fig. 2
figure 2

Topmost image is merged cloud mask (1-km horizontal resolution) for a stretch of tropical RXS along an orbit on 19 April 2007. Black and grey indicate ice and liquid, respectively. Sequences of lower images represent corresponding masks produced by the construction algorithm for various dead-zone lengths as listed

Full size image

Figure 3 shows mean profiles of several variables accumulated out to n as functions of n for the field shown in Fig. 2. Results for accumulated fields, of widths 2n + 1, are shown because the algorithm is intended to produce full 3D domains not single rows. For these accumulations averaging included the original RXS as it is included in constructed fields, yet gave double weight to the reconstructed lines so as to represent scene construction on both sides of the RXS. For clouds higher than 10 km, layer cloud fraction A c , mean cloud water contents 〈w〉, and ν = (〈w〉/σ w ) 2, where σ w is standard deviation of w, were reproduced extremely well for all n. The largest errors on A c occurred for n = 20 at altitudes near 5 km. 〈w〉 contents were reconstructed very well except for n = 5 and 10 for clouds between 2 and 5 km high where associated values of ν were overly small, indicating that these errors were the result of having selected a small number of columns with anomalously large w. On the other hand, for n = 20, almost all reconstructed clouds below 10 km lacked sufficient horizontal variability as indicated by values of ν being twice as large as they should be, despite corresponding 〈w〉 being fine. This was due to too many occurrences of a small number of RXS columns being used multiple times.

Fig. 3
figure 3

Domain-average profiles of layer cloud fraction, cloud water content \(\left\langle w\right\rangle\), and \(\nu =\left(\left\langle w\right\rangle /\sigma_{w}\right)^{2}\), where σ w is standard deviation of w, for the RXS (actual) shown in Fig. 2 as well as for reconstructed cross-sections corresponding to various dead-zone lengths as listed

Full size image

As mentioned, the exemplary field used here was a difficult case of multi-layer tropical cloud. Most 400-km sections of cloud do not exhibit as much intricacy as this one. Numerous other examples were examined and almost all reconstructions performed equally well or better than those shown here. In general, the more extensive and planar the clouds and the fewer the number of definite layers, the better the reconstruction.

4 Radiative Transfer Model

Broadband shortwave (SW) fluxes and radiances were computed by a 3D Monte Carlo algorithm (Barker et al. 1998, 2003). For this study, full 3D results pertain to horizontal grid-spacings \(\Updelta L\simeq 1\) km. Cyclic horizontal boundary conditions were employed. Independent Column Approximation (ICA) results were produced by the same model using \(\Updelta L=10^{8}\) km. Gaseous transmittances (H2O, CO2, O3) were computed using the correlated k-distribution method with 31 quadrature points in cumulative probability space (Scinocca et al. 2008). Optical properties for liquid droplets (Dobbie et al. 1999; Lindner and Li 2000) and ice crystals (Fu 1996; Fu et al. 1998) were resolved into four bands. Aerosol effects were not included but Rayleigh scattering was. For each grid-cell a cumulative distribution of extinction was computed including Rayleigh scattering, gases and cloud particle sizes segregated into 23 radius bins progressing in width following a power law from 0–1 μm to 128.0–161.3 μm. Size-integrated Mie scattering functions were used for each bin, and phase functions for ice crystals were treated as though they were spherical droplets. To reduce variance in radiance estimates, Mie functions were used for a photon packet’s first 5 scattering events with the Henyey-Greenstein function, with proper asymmetry parameter, thereafter if needed (Barker et al. 2003).

The majority of data used here were collected over the Pacific and Indian Oceans. Surface reflectance for oceans was represented by the Cox and Munk (1956) ergodic wave-facet model. Probability of reflection and direction of reflected photons were derived from zenith angle of incident photon and surface wind speed as reported in the A-train dataset (Bloom et al. 2005).

Land was treated as Lambertian with spectral albedos set to values inferred from MODIS data as reported in the A-train dataset (Rutan et al. 2009). For the model’s 0.2–0.7 μm band, values obtained from 0.469 μm radiances were used. For the three near-IR bands, values for 1.24 μm were used.

5 Radiative Transfer Simulations

Radiative transfer simulations were performed on constructed domains that measured ∼61 km wide by ∼14,000 km long and consisted of 85.4 × 106 columns and ∼113.6 × 106 cells. Going much wider with the current version of the construction algorithm is not recommended (see Figs. 2, 3 as well as Barker et al. 2011). Less memory was required than one might think as optical properties were unique only for RXS profiles and each off-nadir profile was indexed to its RXS donor. While photons showered entire domains uniformly, the outer 5 km in the across-track direction and 500 km at the along-track ends were omitted from the analysis. The omitted portions did, however, serve in the 3D simulations as buffer-zones that at once shielded the inner-domain from potentially adverse effects set-up by cyclic horizontal boundary conditions and facilitated horizontal transport of photons through the sides of the inner-domain.

Simulations were performed using solar geometry at satellite time-of-crossing. As such, solar azimuth angle relative to A-train ground-track φ0(i, j) and solar zenith angle θ0(i, j) were set to values reported in the merged A-train dataset. Photons were injected uniformly across constructed domains and given an initial weight of μ0(i, j) = cos θ0(i, j).

6 Results

Two lines of investigation were explored and these are discussed here. First, TOA reflected fluxes inferred from CERES radiances were compared to their simulated counterparts which were obtained by applying the 3D radiative transfer model to constructed A-train domains. Second, TOA reflected solar fluxes obtained by running the radiative transfer model in 3D mode were compared to results obtained by running it as a 1D model.

6.1 Reflected TOA Fluxes: CERES Versus Modelled

Comparison of CERES quantities to modelled values constitutes an attempt to perform a radiative closure experiment as CERES data were not used in the retrieval process. Figure 4 shows the CloudSat quick-look image for granule 999. While the constructed 3D domain was 14,000 km long and 61 km wide, the dots on the track run from ∼44°S to ∼72°N and demarcate the portion of the orbit considered for comparison.

Fig. 4
figure 4

Upper panel shows a CloudSat quick-look image for one of the orbits used here. The line indicates the satellite’s ground-track and the dots on it mark the start (∼44°S) to the end (∼72°N) of the portion used for radiation calculations. Solar geometry is shown in the lower left while solar zenith θ0 and azimuth \(\varphi_{0}\) angles along with surface wind speed for the marked area on the quick-look image are shown on the graph

Full size image

Figure 5 shows TOA fluxes inferred from CERES radiances (Loeb et al. 2005, 2006) as well as model estimates for the constructed domain. Each value corresponds to a (20 km2) CERES pixel. 256 × 106 photons were injected over the full domain making for ∼120 × 103 per CERES footprint with maximum Monte Carlo noise errors for TOA reflectance of ∼0.002. Visually, the agreement is strikingly good for the most part, and encouraging too given the lack of synergistic retrieval and relatively little attention paid to ice scattering properties and surface albedo and no attention given to aerosols.

Fig. 5
figure 5

CERES and 3D Monte Carlo estimates of TOA reflected solar fluxes for the marked portion of the orbit shown in Fig. 4. The 60-km-wide MODIS 0.645 μm image is shown along the top

Full size image

Figure 6 shows the cumulative histogram of differences between modelled–CERES fluxes (neglecting values between 50°N and 60°N where some cloud-related data were missing). Only about 30% of the differences fall within the ±10 W m−2 goal of EarthCARE (ESA 2001; Domenech et al. 2011) and the bias for this particular orbit is −15 W m−2. This large bias, however, arises primarily from the near-cloudless, high-Sun cases between 5°S and 25°N. Underestimation of albedo in this region could result from undetected small clouds (that were not included in the model), systematically too small ocean albedo (stemming from underestimated surface wind speeds) or neglect of aerosol. The plot also shows the cumulative histogram of differences between modelled–CERES nadir radiances scaled by π in order to be of similar magnitude to fluxes. The fact that histograms for fluxes and radiances are almost identical suggests that, for this orbit, there is no clear advantage to using either nadir radiances or fluxes to assess the retrievals. It is worth pointing out, however, that while radiances are a direct measurement CERES fluxes carry the additional error of radiance-to-flux conversion. It is expected that inclusion of off-nadir radiances, as in EarthCARE, will make for a more challenging assessment of retrievals.

Fig. 6
figure 6

Cumulative frequency distributions of (3D model–CERES) TOA reflected solar fluxes for the data shown in Fig. 5. Also shown is the distribution of differences between CERES nadir radiance and 3D Monte Carlo nadir radiances (multiplied, for plotting purposes, by π so as to put them on roughly the same scale as fluxes)

Full size image

Note also that measured radiances, and hence 1D MODIS retrievals, are subject to radiative smoothing and roughening (e.g., Marshak and Davis 2005b). Thus, when they are used at the 1 km scale to initialize a radiative transfer model, even a 3D model as here, secondary smoothing and roughening takes place. This undoubtedly introduces errors into the comparison that are far from straightforward to estimate. The hope is that column properties retrieved by a robust synergistic algorithm will rely less on passive retrievals and thereby reduce smoothing/roughening errors. There is also the possibility that an expanded rendition of the construction algorithm might partially address the effects of horizontal transport (cf., Barker et al. 2011).

Before leaving this section it is interesting to give an example that goes well beyond the initial steps of radiative closure using domain averages and consider the much more demanding case of imagery reconstruction. Figure 7 shows a 50 × 150 km domain of MODIS cloud optical depth τ and its construction algorithm counterpart. Next to this is the corresponding 0.645 μm MODIS image and nadir reflectances produced by the 3D model using constructed τ. While the reconstructed radiances are not perfect, the main features were captured well and will only be that much better once full synergistic retrieval comes online.

Fig. 7
figure 7

Image on the extreme left (actual) shows cloud visible optical depths inferred by applying an inverse 1D RT model to MODIS radiances for a 50 × 150 km scene. Next to it is its 3D constructed version. To their right are MODIS radiances at 0.645 μm, and next to it corresponding radiances from the 3D RT model using the constructed field. Solar zenith angle was 40° and the Sun was coming in from the direction shown

Full size image

6.2 Reflected TOA Flux Differences: 1D Versus 3D

So, while a serious attempt at radiative closure using the current A-train dataset would be premature, on account of known difficiencies in retrieved properties it is likely that cloud structural features are defined well enough to conduct a comparison between 1D and 3D mean solar fluxes for mesoscale-size domains. In so far as the RXSs and their 3D constructed counterparts capture realistic cloud structure down to the 1 km scale with sufficient accuracy, results presented here constitute a fair assessment of radiometric errors due to neglect of 3D solar transport.

When comparing domain-average characteristics of 1D and 3D radiative transfer it is common to employ standalone domains that are assuredly isolated in space via cyclic horizontal boundary conditions. This is because one knows from the start that all photons either get absorbed by the domain’s atmosphere or underlying surface, or exit through the top. An arbitrarily delineated domain in the real atmosphere, however, is not isolated, cyclic boundary conditions do not apply and a clean comparison between 1D and 3D transfer becomes complicated—the more so as the domain’s horizontal extent diminishes. For this reason, particular attention was paid in these experiments to fluxes through the vertical sides of domains.

Results for each simulation were averaged into 50 × 50 km sub-domains, a nominal area representing weather and climate model grid-cells. As each simulation used 256 × 106 photons, each sub-domain received ∼750 × 103 photons. This translates into very small Monte Carlo errors for domain-average fluxes and heating rates, ∼2.5 times smaller than those for results presented in the previous sub-section.

Results were partitioned into sub-domains with total cloud fractions A c : (i) 0.05 < A c  < 0.3; (ii) 0.3 ≤ A c  ≤ 0.7; and (iii) A c  > 0.7. This is because previous aforementioned studies have suggested that cloud structural effects should maximize at intermediate cloud amounts and taper off for larger and smaller A c . Data from five orbits were used and the number of 50 × 50 km sub-domains for the three ranges of A c were: (i) 205; (ii) 246; and (iii) 852, for a total of 1,303 sub-domains.

Figure 8 shows flux differences (3D–1D) as functions for μ0 for the three ranges of A c . For small cloud fractions the only clear signal is greater reflected flux at TOA when 3D transfer was considered. The same trend is visible for A c  ∈ [0.3,0.7]. This is the well-documented effect of cloud-side illumination for \(\mu_{0}\,\lesssim\,0.8\) (Welch and Wielicki 1985). Conversely, at large μ0 there is a slight tendency for 3D transfer to produce smaller reflectances due to photons emerging from cloud-sides in predominantly downwelling directions. Additionally for A c  ∈ [0.3,0.7] there are trends for larger atmospheric absorption and weaker surface absorption (irradiance) by 3D transfer as μ0 decreases. Again, this is due to increased cross-sectional area of cloud presented to direct-beam as μ0 decreases. The vast majority of differences for TOA reflected and surface absorption are smaller than ±20 W m−2; for atmospheric absorption, differences are largely confined to ±5 W m−2.

Fig. 8
figure 8

Differences between fluxes predicted by 3D and 1D transfer for (50 km2) scenes constructed from A-train data. Left column is for reflected flux at TOA, centre is for atmospheric absorption, while the rightmost is for surface absorption. Each row corresponds to scenes with A c as listed on the plots. Grey lines are least squares fits, with the correlation coefficient r > 0.1

Full size image

When cloud fractions become greater than 0.7, the dominant feature is a large scatter of positive and negative differences that cluster around zero with a very slight tendency for 3D transfer to yield smaller reflectances by about 3 W m−2 and larger atmospheric absorption by ∼2 W m−2. These values are not too dissimilar to those obtained by Cole et al. (2005b) using 2D model atmospheres, and possibly stem from photons undergoing internal reflections between multiple cloud layers thereby increasing photon pathlengths, number of scattering events, and hence absorption. As μ0 decreases, the direct beam can undercut layers of clouds when 3D transfer is admitted; for 1D transfer, however, photons must diffuse through upper clouds before undergoing internal reflections.

Figure 9 shows cumulative frequency distributions of differences in reflected flux at TOA predicted by 3D and 1D radiative transfer for all sub-domains partitioned into the three cloud fraction ranges. Median values for A c  < 0.7 are very close to zero but for cloudier scenes it is −8 W m−2. It can be expected that errors in TOA reflectance will exceed ±10 W m−2 for 10–30% of scenes having A c  < 0.7 if 1D rather than 3D transfer is used. This rate jumps to ∼50% of scenes with A c  > 0.7. Preliminary results suggest that these percentages increase by 15–20% when smaller domains the size of CERES footprints are considered. The implication therefore for EarthCARE is that 30–65% of its constructed scenes can expect errors in computed reflected short wave (SW) flux at TOA to exceed its ±10 W m−2 target if 1D rather than 3D transfer is adhered to. More dramatic errors in excess of ±50 W m−2 appear as though they will be very rare regardless of cloud conditions.

Fig. 9
figure 9

Cumulative frequency distributions for differences in reflected solar flux at TOA as predicted by a 3D transport model and its 1D counterpart for 50 × 50 km constructed scenes (i.e. the values shown on the left column of Fig. 8)

Full size image

Figure 10 shows surface and atmospheric absorption histograms that correspond to those in Fig. 9. Regarding surface absorption, the most interesting point is that errors due to use of 1D transfer are often large (i.e. greater than 20 W m−2) for A c  > 0.3 and are dominated by overestimation relative to 3D values. This stems from no illumination of cloud-sides in 1D calculations, especially apparent for A c  ∈ [0.3, 0.7]. Biases in total atmospheric absorption due to use of 1D models appear to be of secondary importance for (50 km2) domains as they rarely exceed ±10 W m−2. Naturally more extreme localized differences exist at smaller scales but they are largely averaged out at 50 km.

Fig. 10
figure 10

As in Fig. 9 except these distributions are for atmospheric and surface absorption (i.e. data shown in the centre and rightmost columns of Fig. 8)

Full size image

6.3 Horizontal Transport Through Domain Sides

For a column of atmosphere associated with a 50 × 50 km sub-domain, let R, T and A be net fluxes out the top, at the surface and absorbed in the column. Since the sub-domains considered here were effectively standalone and not subject to cyclic boundary conditions,

$$ R+T+A\neq \mu_{0}F_{0} $$
(6)

holds in general where F 0 is incident flux on a normal surface at the TOA. Horizontal transport of radiation through the sides of a sub-domain can be defined as

$$ H=\mu_{0}F_{0}-\left(R+T+A\right). $$
(7)

Marshak and Davis (2005b) provided an in-depth discussion of H for solar transfer through cloudy atmospheres with (7) corresponding to their (12.13). Figure 11 illustrates how to interpret (7). For 1D transfer (6) is an equality and so H = 0. For 3D transfer, however, if the shadow of a cloud outside the sub-domain is cast into the sub-domain, but the majority of reflected photons come from, for example, a large planar cloud aloft (so that R 1D  ≃ R 3D ), H > 0 which might seem odd given that transport through the column’s sides is less than it would have been had the external cloud been absent. If the external cloud is now moved so that direct-beam radiation undercuts the layer cloud and some radiation gets scattered by the external cloud back into the sub-domain, H < 0. Since non-zero H can arise through differences between 1D and 3D values of R, T or A as well as actual horizontal flux through vertical boundaries, perhaps it is best to think of H as implied horizontal flux.

Fig. 11
figure 11

Schematic diagram illustrating horizontal transport H of solar radiation. The column of concern is in the centre of each figure and consists of an overcast cloud. R, A and T denote reflectance, atmospheric and surface absorptance, respectively. Left figure shows 1D results where there is no transfer through the column’s vertical sides and so H 1D  = 0. The other figures show the impact on H 3D due to clouds outside the column of concern

Full size image

Figure 12 shows histograms of H for the three classes of sub-domains using data from the five orbits considered here. Values of \(\left\vert H\right\vert \gtrsim 50\) W m−2 occur for ∼10% of the scenes with A c  > 0.7 and point directly to the potential importance of 3D effects on the radiation budget of atmospheric volumes resembling those of global model grid-cells. Although the scale at which H was evaluated here (50 × 50 km domains) and the nature of the domains differ much from those assessed by Marshak and Davis (2005b), values are similar to those they arrive at (cf. their Fig. 12.4). When sub-domain size is reduced to that of a CERES footprint, distributions of H for cases with A c  ∈ [0.3,0.7] and A c  > 0.7 resemble one another closely but the fractions of occurrences of \(\left\vert H\right\vert \gtrsim 50\) W m−2 increase to roughly 20%. Again then, it appears to be crucial that 3D solar transfer be utilized by the EarthCARE mission if it is to make a serious go at both estimating surface-atmosphere radiation budgets and assessing the quality of active-passive retrievals by demanding that simulated TOA fluxes be within ±10 W m−2 of values inferred from its broadband radiometer measurements for (10 km2) scenes (Domenech et al. 2011), especially for broken cloud scenarios.

Fig. 12
figure 12

As in Fig. 9 except these distributions are for implied horizontal transport of radiation as defined in (7)

Full size image

Finally, Fig. 13 shows sub-domain average heating rates computed by the 3D model for a 13,000 km stretch of the orbit on 5 July 2006. It also shows percentage differences between 3D and 1D transfer. The thing to point out is that differences rarely exceed ±20%. When they do, however, the bias persists from the surface to the base of those clouds that are responsible for maximum heating. Once a larger number of orbits get processed the intention is to examine heating rates for scenes sorted according to various cloud properties including cloud-type (e.g., Xu et al. 2005; Bacmeister and Stephens 2010).

Fig. 13
figure 13

Beneath the 60-km-wide visible MODIS image are solar heating rate profiles computed by the 3D Monte Carlo model using 50 × 50 km sections of the 13,000-km-long constructed domain for the orbit shown in Fig. 4. The lower panel shows percentage differences between solar heating rates predicted by the 3D and 1D transport models

Full size image

7 Conclusions and Recommendations

The primary purpose of this study was to begin to get a global picture of errors in estimates of solar fluxes, for cloudy domains roughly the size of those found in global models, that are incurred by applying the 1D independent column approximation (ICA). Both the ICA and 3D benchmark values were obtained by applying a Monte Carlo photon transport model to cloud fields constructed from A-train satellite data. The results are of relevance to both global models that use ICA-based methods to compute mean radiative fluxes (e.g., Randall et al. 2003; Morcrette et al. 2008; Shonk and Hogan 2008) and to radiative closure studies.

3D cloud scenes were created from 2D A-train satellite data via application of the Barker et al. (2011) 3D scene construction algorithm. Their algorithm was developed for the EarthCARE mission which intends to perform continuous radiative closure experiments to assess the quality of its active-passive retrievals. As many A-train cloud profiles consisted of a cloud mask only, profiles of cloud properties had to be allocated subject to constraint by MODIS retrievals. Based on preliminary results using cloud system-resolving model (CSRM) data, it appears that the simple technique used here slightly overestimates differences between 1D and 3D transfer models. Hence, one way view of the results is that they represent an approximate upper-bound for differences between 1D and 3D solutions.

Other limitations of this study were its: inability to address impacts of horizontal fluctuations in cloud extinction at scales smaller than ∼1 km; treatment of scattering by ice crystals as though they were liquid spheres; and neglect of aerosols. It did, however, include important fluctuations at scales larger than 1 km (e.g., Zuidema and Evans 1998), and according to simulations by Cole et al. (2005b), cloud albedo depends weakly on horizontal resolution for grid-spacings less than ∼2 km. The general impression, therefore, was that 3D cloud structural features, plus some information pertaining to microphysical variations, were inferred well enough from A-train data to begin setting the stage for a global statistical analysis of differences between 1D and 3D solar transfer. The quality of this comparison is expected to improve once greater attention is given to ice cloud, aerosol and surface optical properties, and synergistic active-passive retrieval methods mature (e.g., Hogan et al. 2006).

A Monte Carlo transport algorithm was applied to constructed scenes using A-train time-of-crossing values of cosine of solar zenith angle μ0. Calculations were performed on domains measuring 61 × 14,000 km and results were partitioned into 1303 50 × 50 km sub-domains which nominally represent the size of cells found in conventional global models. In corroboration with previous studies, 3D–1D differences are largest for sub-domains with intermediate values of total cloud fractions A c , 3D transfer leading to larger reflectances at TOA and atmospheric absorption, and smaller surface irradiances. Results were, however, somewhat dependent on μ0 with mean differences being almost zero at large μ0 and about ∼10 W m−2 for TOA reflectance, −10 W m−2 for surface irradiance, and ∼5 W m−2 for atmospheric absorption for μ0 ≃ 0.2. Similarly, as shown for a 13,000 km stretch of a single orbit, CERES-inferred and model-predicted solar TOA fluxes agreed to with ±10 W m−2 about 30% of the time. Similar values were found for other orbits.

More so than biases the most defining characteristic of 3D–1D flux differences, regardless of both μ0 and A c , was large random variations. For instance, for A c  > 0.3 it was not unusual to find surface irradiance differences exceeding ±30 W m−2. This scatter of results even occurs for near-overcast conditions. The reason for this may be related to multiple layers of clouds that lead to larger photon pathlengths, and thus slightly larger atmospheric absorption, but can sway reflectance and surface irradiance deviations one way or another.

It is difficult to say whether deviations of these magnitudes will have significant ramifications on simulations of clouds, weather and climate. The few experiments with CSRMs that have attempted to explore deviations like these suggest not (Mechem et al. 2008; Cole et al. 2005a; Pincus and Stevens 2009). But those tests were far from exhaustive. More telling, perhaps, are the Monte Carlo Independent Column Approximation (McICA) experiments with large-scale models (Pincus et al. 2003; Barker et al. 2008; Morcrette et al. 2008). In those studies the amount of radiative noise associated with random sampling of clouds far exceeded that shown in Fig. 8. Nevertheless, there was very little impact on the evolution of simulations.

What these results do suggest is that in order to limit errors in radiative closure studies, regardless of whether TOA fluxes or radiances get used, one should at least conditionally apply a 3D solar radiative transfer algorithm to retrieved/constructed scenes that exhibit highly variable cloud. This is especially so if one seeks to assess retrieval performance for small individual scenes using solar measurements. This is indeed the plan for the EarthCARE mission which, as far as the authors are aware, will represent the first attempt in the atmospheric sciences to employ a 3D radiative transfer model operationally.