The authors propose a new cloud property retrieval technique that accounts for cloud side illumination and shadowing effects present at high solar zenith angles. The technique uses the normalized difference of nadir reflectivities (NDNR) at a conservative and an absorbing (with respect to liquid water) wavelength. It can be further combined with the inverse nonlocal independent pixel approximation (NIPA) of Marshak et al that corrects for radiative smoothing, thus providing a retrieval framework where all 3D cloud effects can potentially be accounted for. The effectiveness of the new technique is demonstrated using Monte Carlo simulations. Real-world application is shown to be feasible using Thematic Mapper (TM) radiance observations from Landsat-5 over the Southern Great Plains (SGP) site of the Atmospheric Radiation Measurement (ARM) Program. For the moderately oblique (45°) solar zenith angle of the available Landsat scene, NDNR gives similar regional statistics and histograms when compared with standard independent pixel approximation (IPA), but significant differences at the pixel level. Inverse NIPA is also applied for the first time on observed high-resolution radiances of overcast Landsat subscenes. The dependence of the NIPA-retrieved cloud fields on the parameters of the method is illustrated and practical issues related to the optimal choice of these parameters are discussed.
It is natural to compare novel cloud retrieval techniques with standard IPA retrievals. IPA is useful in revealing the inadequacy of plane parallel theory in certain situations and in demonstrating sensitivities to parameter choices, parameterizations, and assumptions. For example, it is found that IPA has problems in matching modeled and observed band-7 (2.2 μm) reflectance values for ∼6% of the pixels, most of which are at cloud edges. For simultaneous cloud optical depth–droplet effective radius retrievals (where a conservative and an absorptive TM band are needed), it is found that the band-4 (0.83 μm)–band-7 pair was the most well behaved, having less saturation, smaller changes in nominal calibration, and better overall consistency with modeled values than other bands. Mean values of optical depth, effective radius, and liquid water path (LWP) for typical IPA retrievals using this pair are τ = 22, re = 11 μm, and LWP = 157 g m−2, respectively. Inclusion of aerosol scattering above clouds results in ∼8% decrease in mean cloud optical depth for an aerosol optical depth of 0.2. Degradation of instrument resolution up to ∼2 km has small effects on the optical property means and histograms, suggesting small actual cloud variability at these scales and/or radiative smoothing. Comparisons with surface instruments (microwave radiometer, pyranometer, and pyrgeometer) verify the statisitical adequacy of the IPA retrievals. Last, cloud fractions derived with a simple threshold method are compared with those from an automated procedure called Automatic Cloud Cover Assessment now in operational use for Landsat-7. For the northernmost 2000 scanlines of the scene, the cloud fraction Ac is 0.585 from thresholding, as compared with Ac = 0.563 for the automated procedure, and the full scene values are Ac = 0.870 and Ac = 0.865, respectively. This suggests that the Landsat-7 automated procedure will likely give reliable scene-averaged cloud fractions for moderately thick clouds over continental U.S. scenes similar to SGP.
Cloud optical property retrieval is one of the main goals of the Earth Observing System (EOS) (Wielicki et al. 1995), and improving the algorithms developed to accomplish this goal is an area of active research. Traditionally, to infer cloud optical depth and particle size, plane parallel theory is used. The standard approach is to construct look-up tables for at least two bands, one absorbing and one conservative with respect to water, using a radiative transfer code such as Stamnes et al.’s (1988) discrete ordinate radiative transfer (DISORT) (Baum et al. 2000) or adding–doubling (Minnis et al. 1998) and compare observed and modeled reflectances for each pixel as a function of optical depth, particle size, and the geometry of observation. Various variants of the above methodology have been applied to both aircraft (Twomey and Cocks 1989; Rawlins and Foot 1990; Nakajima et al. 1991) and satellite (Arking and Childs 1985; Minnis et al. 1992; Platnick and Valero 1995) radiometric observations. None of the above studies used any corrections for 3D radiative effects.
With increasing recognition in recent years that 3D radiative transfer can have a significant impact on cloud retrievals, interest in the measurements of the Thematic Mapper (TM) radiometer aboard the Landsat satellites has been renewed. This instrument provides a unique high resolution dataset where the effects of 3D cloud structure are expected to appear more prominently than in Advanced Very High Resolution Radiometer (AVHRR) or other low-resolution imagery. Attempts to account for 3D effects in Landsat optical property retrievals have been scarce and on limited datasets (e.g., Barker and Liu 1995). Instead, plane parallel theory still remains the standard in all Landsat studies (Wielicki et al. 1990; Nakajima et al. 1991; Harshvardhan et al. 1994). What this means in practice, is that the independent pixel approximation (IPA), which neglects horizontal photon transport (Cahalan et al. 1994a,b) is assumed to be a good approximation. Using simulations, Marshak et al. (1998) have shown that as a result of this simplified approach, at high sun, optical depth variability is underestimated, and suggest using inverse nonlocal IPA (NIPA) to bring variability closer to truth. For low sun, on the other hand, with radiative smoothing still present at small scales, “roughening” (or “sharpening”) takes place at intermediate scales, resulting in overestimates of optical depth variability (Zuidema and Evans 1998; Várnai 2000). The roughening has also been reproduced in Monte Carlo simulations (Marshak et al. 1999; Oreopoulos et al. 1999), and its signature has been traced in power spectra of Landsat radiances (Oreopoulos et al. 2000).
The purpose of this paper is to suggest a new type of IPA that uses normalized differences of nadir reflectivities (hence named “NDNR” method) for an absorbing and a nonabsorbing (for liquid water) wavelength. NDNR has been designed to remove radiative roughening and improve retrievals at oblique solar illuminations. Before NDNR and NIPA are applied on real world observations they are tested on Monte Carlo simulations with bounded cascade clouds (Cahalan et al. 1994a) as input. Since both methods require prior development of an IPA-type retrieval procedure (which can also be used for comparisons), we also include in this paper IPA Landsat retrievals that incorporate (albeit crudely) atmospheric effects, and discuss the limitations and challenges when validation is attempted with surface observations. Our satellite dataset consists of a full Landsat-5 TM scene over the Atmospheric Radiation Measurement (ARM) Program Southern Great Plains (SGP) site in Oklahoma.
The paper starts with a description of the satellite dataset and the methodology for the standard IPA retrievals in section 2, continues with cloud detection and cloud fraction comparisons between two methods in section 3, applies NIPA and NDNR on simulated data in section 4, and presents Landsat cloud retrievals with the various techniques in section 5. We conclude with summary remarks and discussion in section 6.
Satellite dataset and retrieval approach
The satellite image used in this paper (Fig. 1) is from the TM instrument aboard Landsat-5, and covers an area of ∼196 × 185 km2 at a resolution of 28.5-m (6888 × 6489 pixels from which ∼7.5% are blank border “filler” pixels) around the ARM program SGP site near Lamont, Oklahoma. The time of the satellite overpass was 1625 UTC on 24 September 1996. The solar zenith angle (SZA) at that time was 45° and the azimuth angle (AZA) was 138°. The image digital numbers (DNs) for each band Bi (for a summary of TM bands and their wavelength span, refer to Table 1) are converted to (nadir) radiances using the calibration coefficients in the header file accompanying the image files, as provided by the National Landsat Archive Production System (NLAPS), and the exoatmospheric solar irradiances given by Markham and Barker (1987) (Table 1). We use both solar and infrared TM bands to distinguish between clear and cloudy pixels (details are given in the following section) and several combinations of a conservative and an absorbing band pair to retrieve cloud optical depth (τ), effective radius (re), and their by-product, liquid water path (LWP). We look for differences in standard IPA retrievals of τ and re and evaluate their performance with the pairs (B2, B7), (B4, B7), (B2, B5), and (B4, B5). All sensitivity calculations and extensions beyond standard IPA are conducted using the pair (B4, B7). The values B2 and B4 are suitable for retrieving optical depth since cloud droplets are virtually nonabsorbing at the wavelength span of these bands, and radiances depend weakly on cloud particle size. While B1 and B3 also have the same attributes, they are not used here because they show a significantly greater percentage of saturated pixels (i.e., pixels with radiance values above those that can be resolved by the 8-bit TM discretization) than the other two. At B5 and B7 wavelengths there is size-dependent droplet absorption that makes these bands suitable for inferring droplet effective radius. The use of an absorbing and a nonabsorbing pair of radiometer bands for simultaneous τ and re retrievals is, as mentioned in the introduction, a classic procedure. The reader can find the theory and applications for various instruments and conditions in the papers (among others) of Arking and Childs (1985), Twomey and Cocks (1989), Foot (1988), Nakajima and King (1990), Nakajima et al. (1991), Platnick and Valero (1995), and Minnis et al. (1998).
The essence of our implementation of the procedure is a search in (τ, re) space for values at which modeled and observed reflectivities combine to minimize the quantity (Nakajima and King 1990):
Ri = πIi/μ0Fi is the (bidirectional) reflectivity, Ii is the nadir radiance for band i, Fi is the extraterrestrial solar irradiance for that band, and μ0 is the cosine of the SZA; the superscripts “calc” and “obs” refer to modeled and observed reflectivities, respectively. We restrict ourselves to N = 2 (the band pairs mentioned above). As will be shown later in the description of NDNR, one of the pair components can be a quantity other than a single band reflectivity. The calculated values Rcalci are from a modified version of the radiative transfer code described by Tsay et al. (1990), assuming that the top of the atmosphere reflectivity can be approximated as in Platnick and Valero (1995):
where Rcld is the reflectivity of the cloud-surface system, t(μ0) is the atmospheric transmissivity to the incoming beam, t(μ) is the atmospheric transmissivity to the reflected beam, μ (=1) is the cosine of the viewing angle, and Ratm is the reflectivity of the atmosphere above the cloud. As explained in Platnick and Valero (1995), Eq. (2) is a crude approximation suitable only for nonrigorous retrievals, but has the advantage of being flexible for sensitivity analysis. We have also adopted Eq. (2), because it can easily handle many of the tests we performed, and because a more sophisticated method would not be justified given the uncertainties in atmospheric constituents (particularly aerosols), and cloud properties (e.g., cloud top and base height), within the Landsat scene (surface measurements of clouds at SGP are hardly representative of the entire Landsat scene because they are point measurements). The filter function of the radiometer for each band is taken into account in our calculations as in Platnick and Valero (1995). For the calculations the cloud was placed between 1 and 2 km and was assumed to be composed of water droplets (with an exception mentioned later) forming a lognormal drop size distribution with effective radii ranging from 4 to 25 μm and effective variance (defined in Hansen and Travis 1974) of 0.13. Because of the existence of multiple solutions (Nakajima and King 1990), we avoided retrievals below re < 4 μm. Thus, the first bin in effective radius histograms should be understood to contain all droplets with re ⩽ 4 μm. Similarly, the 25-μm bin potentially contains droplets with re ≥ 25 μm. Table 2 summarizes these and other modeling assumptions and inputs for each band.
Pixel identification and cloud fraction retrievals
A step that precedes the retrieval of cloud optical properties is the identification of the cloudy pixels. Thus, we first check whether our simple threshold algorithm is capable of separating clear from cloudy skies, by comparing it with the operational algorithm Automatic Cloud Cover Assessment (ACCA) developed for operational clear/cloud discrimination for Landsat-7 (whose sensor Enhanced TM+ has bands similar to TM) by R. Irish (1999, personal communication).
The properties of clear-sky pixels (for vegetated surfaces as in SGP) that are used to distinguish them from cloudy pixels are: 1) they are usually less reflective than cloudy pixels in bands 1–4 (0.45–0.90 μm); 2) their B2 (0.55 μm) and B4 (0.83 μm) reflectivities (R2 and R4) differ more than for cloudy skies; and 3) their B6 (11 μm) brightness temperatures (BT6) are higher than those of cloudy pixels. Based on these characteristics we select thresholds from R4, R4/R2, BT6 one- and two-dimensional histograms. Based on this threshold selection, we classify pixels as clear when they satisfy all of the following criteria: R4 < 0.3, R4/R2 > 1.6, and BT6 > 285K. A pixel that has not been classified as clear because of failure to satisfy one or more of the above criteria, can be classified as clear at a later stage if its reflectivity value R2 or R4 is lower than the minimum value of the look-up table used for the retrieval of optical depth (calculated for some assumed surface albedo, which can be found in Table 2).
ACCA consists of a “first-pass” series of threshold tests involving bands 1–6 and some quantities formed by their combinations, a “second pass” using a BT6 threshold determined from the results of the first pass, and a final filtering that fills cloud holes whenever more than 50% of a clear pixel’s immediate neighbors are classified as cloudy by the first two passes. For the northernmost segment of the Landsat scene (first 2000 scanlines), which is covered mostly by fair-weather cumulus clouds, ACCA gives a cloud fraction of Ac = 0.563, whereas our (simpler) algorithm gives Ac = 0.585 for the B4 surface albedo, and Ac = 0.576 for the B2 surface albedo. For the entire scene, the numbers are Ac = 0.865 for ACCA, Ac = 0.870 for the B4 surface albedo, and Ac = 0.868 for the B2 surface albedo. The closer agreement in cloud fraction over the full scene is due to the dominance in the bottom ⅔ of the scene of extensive thick clouds that are easily detected by both algorithms.
Despite the agreement of the cloud fraction values from the two algorithms, a detailed examination of the cloud masks reveals that for the fair weather cumulus portion of the scene, discrepancies in the classification of individual pixels are not uncommon, that is, some of the cloud fraction agreement is due to compensating differences. This is not surprising since ACCA was designed only to give good statistical estimates of cloud fraction for large portions (i.e., quadrants) of a Landsat scene with little concern about errors for individual pixels. In contrast, our algorithm is scene specific, and aims at detecting clouds at the pixel level. An overestimation of cloud fraction by our algorithm may be due to the conservative selection of thresholds for which the condition of simultaneous fulfillment cannot sufficiently compensate. Another indication that our algorithm may be overestimating cloud fraction is pixel misclassification at cloud edges. We suspect that this may be taking place because we are occasionally unable to match observed and modeled B7 reflectivities in these regions. On the other hand, cloud fraction underestimation by ACCA could be due to the strictness of requiring that seven thresholds are satisfied during the first pass for a pixel to be classified as cloudy (for the present scene, we found that the second pass does not significantly change the cloudy pixel population). This rather severe demand has its roots in the desire for ACCA to be a global algorithm, and capable of identifying clouds in difficult situations, such as over highly reflective surfaces (e.g., desert, snow, ice). Last, we note that ACCA would underestimate cloud fraction even more without the cloud hole “filling” procedure; when filling is turned off, ACCA values drop to 0.467 for the topmost segment and 0.822 for the entire scene.
Cloud property retrievals on model clouds
Before carrying out optical depth retrievals from the Landsat image, we compare the performance of standard IPA versus the novel retrievals techniques with simulated data. The data consists of nadir reflectivities from 3D Monte Carlo (MC) simulations on 1D bounded cascade clouds generated following Cahalan et al. (1994a). The specifics of the experiments (each consisting of 10 realizations) are given in the caption of the relevant figures. The advantage of using simulations is, of course, that the truth (the optical depth field), is known, so that the various retrieval approaches can be evaluated in a straightforward manner. Note that these simulations do not address the problem of dual (optical depth–effective radius) retrievals, since the appropriate cloud model is not yet available. We plan to study this problem in the future.
Three-dimensional radiative transfer effects are expected to cause errors in IPA retrievals of optical depth. For simplicity, we classify 3D effects into two categories: those causing smoothing and those causing roughening relative to IPA. The relative amounts of smoothing and roughening and their scale dependence are a function of illumination conditions, wavelength, and the properties of the cloud field itself. Smoothing has been discussed in Stephens (1988), Cahalan (1989), Cahalan et al. (1994b), Marshak et al. (1995), Marshak et al. (1998), Várnai (2000), and others, and is due to net horizontal transfer of photons between pixels in conditions where multiple scattering is dominant. As explained in Marshak et al. (1998) smoothing causes IPA retrievals to underestimate the variability of the true optical depth field.
When SZA exceeds ∼40°, roughening starts to occur (Zuidema and Evans 1998; Marshak et al. 1999; Várnai 2000). Roughening is more prominent at absorbing wavelengths and can be seen in power spectra of Landsat radiance observations (Oreopoulos et al. 1999, 2000). Roughening is a manifestation of side illumination and shadowing due to both geometrical and optical thickness gradients. It results in overestimation by the IPA of both the mean and standard deviation of the optical depth field (Zuidema and Evans 1998; Várnai 2000). IPA tends to overestimate optical depth for cloudy columns that receive large amounts of side illumination (i.e., when optically or geometrically thinner clouds are placed between the column in question and the solar beam), and to underestimate optical depth for cloudy columns that exist in the geometrical or optical shadow of thicker cloudy columns. Marshak et al. (2000) have noted that these 3D radiative effects can also be observed from the ground under broken cloud conditions. They indicate that the 3D radiative effects are similar for conservative and moderately absorbing wavelengths, and can be removed using a normalized transmissivity difference ratio.
Based on the above, one can conclude that the following actions need to be taken in order to improve cloud properties retrieved from inversions of high resolution satellite radiances: (a) apply a transformation that roughens the observed field at small scales; (b) apply a transformation that smooths the field at intermediate scales, and (c) account for wavelength (absorption) dependencies. Here, the inverse NIPA of Marshak et al. (1998) is used to tackle (a), NDNR deals with (b), while thorough investigation of (c) (related to effective radius retrievals) will be left as the subject of a future study. Still, we present some initial results of applying NIPA at absorbing wavelengths and demonstrate the potential of a NIPA–NDNR combination for retrieving optical depth.
IPA versus NIPA
This subsection presents results similar to the ones reported by Marshak et al. (1998), with the addition of MC simulations on a cascade cloud exhibiting cloud top variability. Cloud top variability is generated following Marshak et al. (1999) with parameter values given in the caption of Fig. 2.
The core of the inverse NIPA transformation lies in the following equation (for details see Marshak et al. 1998):
where G(k) is the Fourier transform of the cloud Green’s function, approximated by a gamma distribution:
is a stabilization function to deal with the ill-posed nature of the problem. In the above equations RNIPA is the reflectivity obtained after applying the inverse NIPA transformation on the MC reflectivities, RMC; x is the location coordinate, α, η are parameters defining the shape of the gamma distribution (c is its normalization constant), and γ is a free parameter modulating the stabilization function. The transformation is carried out in Fourier space (hence the appearance of wavenumber k). The premise of inverse NIPA is that the RNIPA field is a good approximation to the IPA field corresponding to the actual optical depth distribution, and is more appropriate for lookup table inversion than RMC, since it is devoid of the smoothing due to net photon horizontal transport.
Figures 2a,b compare IPA- and NIPA-retrieved optical depths with the actual optical depths for a 128-column segment of one of the (10) realizations simulated for each case [“flat” and “bumpy,” the latter obtained from the former by superimposing cloud top variations generated from a fractional Brownian motion that are uncorrelated with optical depth as in Marshak et al. (1999)]. Table 3 (second and third row) shows the first two moments of the retrieved optical depth for this realization, which can be compared with the moments of the true optical depth field (first row). Similar results are obtained for the other realizations. Figure 2 shows that indeed, the NIPA field captures better the variations of the true field than the IPA field, which is too smooth. IPA underestimates both the mean and variance of optical depth while NIPA brings them closer to the actual values. The impact, however, of NIPA is far greater on individual pixels: taking into account all 10 realizations the rms error of IPA is 1.63 versus 1.13 for NIPA (flat cloud). Note that the performance of NIPA for a particular realization can be improved by selecting its own (optimal) set of α, η, and γ parameters. The parameters are determined from iterative trials until the NIPA rms errors are minimized. The iterative procedure can be accelerated if the initial values of α and η are the best pair from forward NIPA trials (in which the IPA field is smoothed to approach the MC field). The parameters used in the realization of Fig. 2 are given in the caption.
IPA versus NDNR
In this subsection it is shown that in contrast to the IPA underestimates of optical depth at high sun elevation (previous subsection), inversion of 3D nadir reflectivities at low sun elevation yield overestimates of optical depth mean and variability. This is due to shadowing and side illumination (roughening) arising from low-order scattering. Further, it is shown that retrievals using the normalized index,
formed by the reflectivity at a conservative wavelength Rc and the reflectivity at an absorbing wavelength Ra, behave similarly to IPA retrievals at high sun, and are an improvement over retrievals using only Rc. A visual confirmation of this claim is given in (c) and (d) of Fig. 2 comparing the true optical depth field with the optical depth field retrieved from IPA and NDNR for a segment of one of the 10 cloud realizations (the same shown in a and b). Means and standard deviations are given in the last two rows of Table 3. NDNR retrievals do indeed smooth the optical depth field in a manner similar to that of the IPA retrievals at high SZAs (see Figs. 2a,b), but give better statistics than IPA, which retrieves a much rougher optical depth field. The rms errors of IPA and NDNR for all 10 realizations are contrasted in Figs. 3a,b (curves labeled “IPA” and “NDNR”). It can be seen that NDNR does better overall, especially for the bumpy cloud, though there is evidence of occasional oversmoothing. Note that rms errors are greater than in the SZA = 0° case, suggesting that retrievals are more difficult at low solar illuminations. In the above figures Ra is calculated for ω0 = 0.98, but the results are very similar for ω0 = 0.99 (the volume extinction coefficient and asymmetry factor are the same as in the conservative case). Obviously, the success of NDNR owes to the fact that low order 3D scattering effects that cause roughening have similar patterns for the two wavelengths, yet differ enough in magnitude to be unambiguously mapped as a function of optical depth (of which NDNR is a monotonic function). These requirements are not met if the selected Ra corresponds to either very weakly or very strongly absorbing wavelengths.
IPA versus NIPA–NDNR
Because MC simulations provide evidence that NDNR retrievals at low sun behave similarly to IPA retrievals at high sun (Fig. 2), it would seem appropriate to apply inverse NIPA to the NDNR index to compensate for the smoothing. The rms error of this NIPA–NDNR combination for all 10 realizations of the flat and bumpy bounded cascade is shown in Fig. 3 (curve labeled “NIPA–NDNR”). It can be seen that NIPA on NDNR improves relative to plain NDNR. Note, however, that there is no well-founded physical interpretation in such a transformation since the smoothing we attempt to remove does not originate from a physical process. Therefore, there is no experimental justification for approximating the kernel G(k) with a gamma distribution, or choosing its parameters based on the physical characteristics of the cloud. Despite this, the NIPA–NDNR transformation works quite well.
Cloud property retrievals from Landsat
We start the Landsat analysis by performing the IPA retrievals according to the model described in section 2. We examine various aspects of the IPA retrievals, such as sensitivity to assumptions, the choice of band pair, effect of resolution degradation, and consistency with ground-based measurements. NDNR and NIPA retrievals are discussed in subsections 5b,c.
Cloud property histograms
Figure 4 shows histograms of retrieved τ, re, and LWP = (⅔)τre for various band pairs. The τ distribution (Fig. 4a) looks significantly different for the two bands and is generally not as skewed as in previous Landsat studies (Barker et al. 1996; Harshvardhan et al. 1994). The distributions derived from B2 end abruptly at τ ≈ 35 because of saturation, while the distributions derived from B4 are much better behaved, with a longer tail. As expected, the choice of absorbing band for droplet size retrievals does not significantly affect the τ distribution. When a constant effective radius re = 10 μm is used throughout, the inferred mean value of τ, namely 22.3 for (B4, B7), is reduced only by about 10%. This is another indicator of the stability of the τ retrieval. The means and standard deviations for the different cases are summarized in Table 4.
The re distributions (Fig. 4b) are wider than those for oceanic boundary layer clouds shown by Nakajima et al. (1991) and Platnick and Valero (1995); this is not surprising given the larger spatial extent of the Landsat scene and the greater atmospheric variability over land. As explained in Nakajima et al., the retrievals should be viewed as representing droplet sizes near cloud top. It is obvious that the retrievals using B5 are plagued with problems. First, the re retrievals depend more (as compared with B7 retrievals) on the choice of the conservative band used for τ retrieval. Second, there is a pronounced peak at the uppermost 25-μm bin [0.130 for the (B4, B5) pair] indicating too much apparent absorption. We note that this bin can contain pixels with re > 25 μm, but not pixels with R5 < min(Rcalc5). This peak remains even when the atmospheric absorption of the model is doubled (value changes to 0.079), or when the surface is set to be black (change to 0.081). It can only be eliminated if the calibration coefficients suggested by Thome et al. (1997) are used (change to 0.022), but this results in such overreflective clouds, that ∼50% of the retrievals assign re to the first (4 μm) bin! The frequency histogram of BT6 for the pixels with retrieved re ≥ 25 μm, does not show any peaks at low values of brightness temperature suggesting that there is no systematic tendency for this value to be retrieved for ice-contaminated pixels. Neither is there a preference for these retrievals to occur for very thin clouds where the procedure is less accurate (Nakajima et al. 1991). We must therefore conclude that there is either poor calibration knowledge for B5, or substantial errors have been committed in modeling the B5 reflectivities. The first possibility is reinforced by the fact that among all bands, it is B5’s vicariously determined gain coefficients (Thome et al. 1997) that differ the most from nominal values—about 30%. Given the problems with B2 and B5, we suggest that the (B4, B7) pair is used for cloud retrievals over land whenever a substantial fraction of clouds consists of optically thick (τ > 30) clouds.
The LWP maxima (Fig. 4c) are very similar to the maxima shown by Nakajima et al. (1991) for marine stratocumulus for two of their days (10 and 16 July), but our histograms are enhanced at small values of LWP due to the numerous fair-weather cumulus in the northern part of the scene. The discrepancies between the values derived from the different pairs are consistent with the discrepancies in optical depths (B2 retrieving smaller optical depths than B4), and effective radius (B5 retrieving larger sizes than B7). Our retrieved values of LWP are occasionally much higher than the maximum values inferred by the surface Microwave Radiometer (MWR) [see 5a(3)], but the high values are usually found in the bottom half of the scene and belong to clouds that probably did not advect over the MWR.
Excluded from the histograms is the fraction of clear sky (slightly different for retrievals using B2 and B4—as explained in section 3), the fraction of saturated pixels for B2 (0.043) and B4 (0.013), and the fraction of pixels, which, while classified cloudy, have B5 or B7 reflectivities smaller than the minimum value of the lookup tables (0.01 and 0.06, respectively). The latter pixels occur mostly at cloud edges, thus suggesting a direct influence of cloud geometry on the transfer of solar radiation. The fact that retrievals using B5 have a smaller fraction of these dubious pixels does not necessarily mean that B5 is less affected by radiative transfer at cloud edges; B5 reflectivities are within lookup table bounds, but result in re = 25 μm, which is suspect, given the magnitude of the anomalous peak for that bin (Fig. 4b). Some of the discrepancy is probably also the result of erroneous pixel classification as discussed in section 3, that is, shadow pixels classified as cloudy with apparent surface nadir reflectivity smaller than our universal B5 and B7 lookup table values.
The benchmark (B4, B7) retrievals are repeated with the addition of a continental aerosol layer of optical thickness 0.2 above the cloud, which was assumed to affect only B4 reflectivities. The presence of aerosols triggers a drop in the mean-retrieved optical depth to 20.5 (from the original 22.3) and a drop of the maximum value to 65 (from the original 73). The drop is, of course, due to the greater atmospheric reflectance (∼0.023 at the center of B4) in comparison with the pure Rayleigh scattering case, and the resulting decrease in the apparent cloud reflectivity [see Eq. (2)]. The histogram of optical depth retrieved in the presence of aerosol is compared to the original histogram in Fig. 4d; the effective radius retrievals remain practically unaffected.
As already mentioned, when re is set to 10 μm, the mean optical depth drops to 20.3, and when cloudy pixels with BT6 less than 266 K (∼10% of the cloudy pixels) are considered to be made exclusively from ice crystals of single effective size with constant scattering properties within the B4 wavelength range, the mean optical depth further decreases to 19.3. The latter decrease is due to our use of the International Satellite Cloud Climatology Project (ISCCP) ice phase function (Mishchenko et al. 1996), which has smaller asymmetry factor (more backscattering).
Figure 5 shows the effects of coarsening the measurement resolution on the retrieved (B4, B7) means. The coarsening is accomplished by successive averaging of 2 × 2, 4 × 4, and so on, arrays of raw count values, followed by implementation of the IPA retrieval procedure (including cloud masking) as before. Coarsening does not have a substantial effect on the histograms (not shown), while the changes in the mean values (∼4% for optical depth and ∼5% for effective radius at 64 × 64 pixel degradation) suggest that the effect of assuming homogeneous pixels and thus introducing plane parallel homogeneous (PPH) bias (Cahalan et al. 1994a; Davis et al. 1997; Oreopoulos and Davies 1998) is small for the scene as a whole for degradation up to ∼2 km. Larger differences do, of course, occur locally. The apparent homogeneity at these scales is probably a combination of weak actual cloud heterogeneity and radiative smoothing (Marshak et al. 1995). Careful examination of Fig. 5 reveals that, contrary to the direction of the PPH bias, there is a slight increase of the retrieved mean optical depths up to a certain degradation of resolution (4 × 4 pixels), which is due to saturated values (not used in the retrieval process at the original resolution) entering the averaging process and forming nonsaturated degraded pixels. This effect is still present at coarser resolutions, but its impact is smaller than that stemming from the convex dependence of reflectance on optical depth, and the inclusion of clear pixel contributions in the degraded reflectivities. The mean effective radius dependence on resolution on the other hand, follows the averaging laws of concave functions:R7 drops with re for a constant optical depth, so degrading leads to retrievals of higher effective radii.
LWP comparisons with MWR
Around the time of the satellite overpass, the surface MWR at the SGP site [described in the ARM Web site http://www.arm.gov/docs/instruments/static/mwr.htm—also check Liljegren (1995)] reports LWP values that are in general smaller than the values shown in the histograms of Fig. 4c. The LWPs inferred from MWR are actually more consistent with the Landsat retrievals for the northernmost ∼2000 scanlines (which contain the instrument site) of the Landsat scene (Fig. 1). Obviously, even after assuming that the clouds that have passed over the MWR instrument remain unaltered until the time of satellite overpass (as in the “frozen turbulence” hypothesis), it is very difficult to pinpoint the location to which they have been advected within the Landsat scene. Thus, we use the Landsat scene’s topmost 2000 scanlines and the 1-min MWR measurements centered in an 8-h window around the satellite overpass, for LWP histogram comparison. We discard values< 20 gm−2 for both satellite retrievals and surface measurements, and renormalize the histograms with respect to the remaining values (Fig. 6). The justification for not taking into account the small values is that the residual error of the multiple linear regression used to convert measured brightness temperatures to LWP inherently sets a threshold for unambiguous cloud detection (see Cahalan et al. 1995 and the Web address given above). The mean values are 55.2 g m−2 for MWR, 60.9 g m−2 for the (B2, B7) pair, and 72.2 g m−2 for the (B4, B7) pair. The clear sky fraction (<20 g m−2) is 0.47 for MWR, 0.53 for the (B2, B7) pair, and 0.44 for the (B4, B7) pair. Bear in mind that part of the discrepancy in these comparisons is also due to the different resolution of the instruments: TM resolution is ∼30 m while MWR resolution depends on cloud-base height and the wind speed that determines the rate of cloud advection. Since the MWR data are of coarser resolution (∼500 m), one would expect the difference in the mean values between satellite and surface retrievals to be reduced slightly upon coarsening the satellite resolution. There is indeed a small improvement; for 16 × 16 averaging for example, the (B2, B7) pair has a mean of 57.8 g m−2, and the (B4, B7) pair has a mean of 70.4 g m−2.
Pyranometer comparisons with model calculations
The temperature and water vapor profile measured by the radiosonde launched at the Central Facility at 1429 UTC and the retrieved cloud optical properties from the (B4, B7) pair are used as input in the National Center for Atmospheric Research (NCAR) Column Radiation Model (CRM) version 2 (Briegleb 1992), in order to calculate the broadband shortwave and longwave fluxes reaching the surface. Comparisons are made with the surface observations at the Central Facility location in SGP where pyranometers and pyrgeometers operate continuously as part of the ARM Program. The 1-min time series of shortwave and longwave fluxes at 2 m above the surface are shown in Fig. 7a for the period 1500–1800 UTC of 24 September. The CRM calculations use pixel arrays of various sizes centered around the Central Facility. Results from two array sizes are shown as triangle and diamond symbols in Fig. 7a. These are averages of independent column calculations. There is very good agreement in the longwave (only one of the two points can be seen because the results for the two array sizes are virtually identical), while shortwave values do not agree as well because of the very pronounced shortwave flux variability at the time of Landsat overpass. This variability is clearly due to broken clouds. For pixel arrays ranging from 32 × 32 to 160 × 160, the calculated values change very little, partly because the fraction of clear sky (∼70%) in the arrays is significantly higher than the cloudy fraction. Thus, the comparison of shortwave flux averages cannot be considered as a robust validation of our cloud retrievals. However, the range of calculated SW and LW single-pixel flux values is realistic, as shown by the transect results of Fig. 7b. For SW, the minimum single-pixel value for the 128 × 128 array centered around the Central Facility is ∼225 W m−2 (compare with pyranometer local minimum value around the time of Landsat overpass of ∼215 W m−2) and the maximum value (corresponding to pixels identified as clear) is 729 W m−2 (compare with pyranometer value of 795 W m−2). Note that for broken cloud conditions it is not unusual for the solar irradiance beneath clouds to exceed that of the clear sky (Ackerman and Cox 1981; Leontieva and Stamnes 1996) because of the downward streaming of solar flux from cloud sides, something the CRM cannot simulate. The calculated LW surface fluxes range from ∼341 to 399 W m−2, consistent with the range of pyrgeometer values shown in Fig. 7a. Some of the remaining discrepancies owe to the fact that contributions to surface radiometer measurements come from sky areas much larger than that of a Landsat pixel.
For Landsat retrievals NDNR is defined as
This definition differs from that used in the MC analysis in that it includes atmospheric, cloud drop size, and surface effects, but is still suitable for inversion since NDNR remains a monotonic function of visible optical depth (Fig. 8). For (τ, re) retrievals, NDNR takes the place of R4 or R2 in Eq. (1).
Fig. 4d shows that the retrievals from the (NDNR, B7) pair result in an optical depth histogram that is very similar to that of the (B4, B7) pair. While this is expected for effective radius retrievals, it is somewhat surprising for optical depth retrievals. Similarly, NDNR does not significantly change the (B4, B7) LWP histogram of Fig. 6. However, further analysis reveals that despite the similarity of the histograms, pixel-level optical property differences can be substantial. Figure 9 shows rms differences for 25 subscenes [of size 1024 × 1024 pixels, i.e, ∼(30 km)2] into which the original scene was divided. For some subscenes, rms differences exceed 3 for optical depth; effective radius rms differences are small (as expected) with the exception of the first few subscenes with numerous broken clouds. The overall rms for the entire scene is 2.30 for optical depth and 1.39 μm for effective radius. At this point, there is no clear way to assess the quality of NDNR retrievals for Landsat. Analysis of subscene results reveals that NDNR-retrieved optical depth fields are not always smoother than their corresponding IPA-retrieved fields. This deviation from MC behavior may be due to surface and atmospheric effects, R7 dependence on re, and the sun angle not being low enough for substantial roughening. Indeed, at SZA = 45°, MC simulations showed that radiative roughening due to low-order scattering is not as strong as at 60° (which itself depends on the specifics of the cloud field), and that NDNR retrievals do not always improve upon IPA (not shown).
Sensitivity to parameters
We have applied inverse NIPA on the observed B4 and B7 radiances of several 256 × 256 [i.e., ∼(7.6 km)2] overcast segments of our Landsat scene, while retrievals of optical depth and effective radius were performed as before. NIPA fields are obtained from a straightforward extension of Eq. (3) to 2D, with observed reflectivities R4 replacing RMC, followed by look-up table inversion. Figure 10 compares the (B4, B7) optical depth field retrieved from IPA with that from (NDNR, B7) and several incarnations of NIPA for one of these subscenes. The numbers next to each panel (except IPA, which is used to compare all others) are the optical depth rms differences from the original IPA retrievals. Each panel is a different incarnation resulting from different sets of the parameters α, η of the NIPA kernel G [see Eq. (4)], and γ, the parameter appearing in the regularization function f [Eq. (5)]. At present, because of insufficient understanding of how NIPA should be applied at absorbing wavelengths, we use for B7 the same parameters for B4. Since smoothing is weaker in absorbing bands, a reasonable assumption perhaps would be to use parameters α and η that would make G peak closer to the photon’s entry point [in the searchlight problem of Marshak et al. (1995)]. However, since our emphasis in this work is on τ rather than re retrievals, we ignore this difference for now. It can be seen that the NIPA transformation does indeed lead to optical depth fields more variable than IPA fields, but only for a certain range of the regularization parameter γ [a fact stressed by Marshak et al. (1998)]. If, for example, γ is changed to 10 times its Fig. 10 value, NIPA will give smoother instead of rougher fields (field not shown, see histogram of Fig. 11). For parameter values η = 0.025, α = 0.5, the standard deviation στ of the optical depth field is 7.36 compared to 6.87 for IPA. Here, στ increases to 8.18 for η = 0.035, α = 0.5 and is 8.04 for η = 0.025, α = 1.0. From the above it is obvious that increasing either η or α results in rougher reflectivity and optical depth fields. Note that none of the above transformations affects significantly the mean value of optical depth (∼24) as explained in Marshak et al. (1998), and that NIPA is stable for small changes in the parameters α and η. Figure 11 shows histograms of the τ and re differences between the original IPA and four NIPA retrievals (including two cases from Fig. 10). In the histogram representation an increase (decrease) in the value of α and η manifests as an increase (decrease) in the width of the histograms, indicating more pronounced deviations from IPA. Using for the NIPA transformation of B7 the same parameters as for B4, leads to changes in re, which are mostly confined to ±2 μm. The re difference histograms would have been even narrower had the α and η parameters been chosen to reflect the fact that smoothing is weaker at absorbing wavelengths.
Contrary to the controlled experiments with model clouds where the true optical depth field is known and can be used to guide the choice of parameters and to evaluate the performance of NIPA (as in Figs. 2a,b), the parameter values used in Figs. 10 and 11 were chosen only to illustrate the magnitude of roughness change and its direction. The process of choosing optimal parameters for the NIPA transformation will be aided by experience gained from MC simulations. Another possibility is to use scale-by-scale (spectral) analysis. For example, Fig. 12 shows the power spectra of optical depth for IPA and two of the NIPAs of Fig. 10. Smoothing (steepening of the slope) is evident in the IPA retrievals, but not in the NIPAs, although the NIPA with α = 1, seems to be overroughening. At this stage, we cannot state with certainty which of the NIPA transformations of Figs. 10 and 11 is closer to the truth. An important difference between NIPA experiments with simulated and with Landsat data is that in the latter case roughening is applied on reflectivities that include atmospheric and surface effects, which can either smooth or roughen cloud reflectivities, depending on scene specifics (Oreopoulos et al. 2000). A possible remedy to this problem would be to apply the NIPA transformation directly on the retrieved fields instead of the reflectivities.
NIPA on NDNR retrievals
For Landsat observations, the roughening can be applied on either the NDNR values alone or to both NDNR and R7. Note, that even in the first case, the effective radius retrievals are affected, albeit slightly, inasmuch as the retrievals require minimization of the corresponding χ2 [see Eq. (1)]. Application of NIPA to absorbing bands (or indices involving absorbing bands such as NDNR) is, as already mentioned, only an empirical extrapolation of our understanding of NIPA at conservative bands. Here, we examine behavioral differences when the same NIPA parameters are used on both the (B4, B7) and (NDNR, B7) pair.
The NDNR optical depth field is compared with a realization of the corresponding NIPA–NDNR field in Fig. 10. The effect of NIPA is similar on NDNR retrievals as on the standard (B4, B7) retrievals: στ increases from 7.28 for IPA–NDNR to 8.31 for NIPA–NDNR. Figure 11 also shows histograms of τ and re differences between the original IPA–NDNR retrievals and several of their NIPA counterparts (right column) for the subscene of Fig. 10. As discussed already, increases in the values of the η and α parameters leads to wider histograms. The effective radius histograms are not much affected by the method chosen for the retrieval of optical depths. The histograms are only slightly asymmetric around zero, indicating small effects on the regional mean of the retrieved optical properties. The most interesting result divulged by the histograms are greater deviations from IPA by NIPA–NDNR compared to NIPA on B4 when identical NIPA parameters are used. Investigating whether one of the two NIPA transformations is consistently better calls for further MC studies.
Summary and conclusions
We have introduced a dimensionless normalized difference index (NDNR) formed by conservative and nonconservative reflectivities and shown that it can be used to improve optical depth retrievals. Monte Carlo simulations show that NDNR retrievals at oblique sun illuminations perform better than IPA retrievals, eliminating optical depth extrema due to intermediate scale shadowing and side illumination, and giving more accurate statistics. NDNR retrievals at low sun have similar characteristics as standard IPA retrievals at high sun. For Landsat we defined NDNR = (R4 − R7)/(R4 + R7), where R4 and R7 are reflectivities of a conservative and an absorbing Thematic Mapper band. Retrievals with NDNR for a September 1996 Landsat scene over the Oklahoma ARM site do not have a significant effect on overall statistics and histograms of optical properties when compared with IPA, but introduce important differences at the pixel level, especially for optical depth. The inverse NIPA of Marshak et al. (1998) was applied for the first time on satellite observations (in this case overcast 256 × 256 pixel portions of the Landsat scene) in order to remove small-scale radiative smoothing and to examine the dependence of the degree of roughening on the parameters of the NIPA kernel. MC experiments with bounded cascade clouds illuminated at SZA = 60° show that hybrid NIPA–NDNR retrievals are better than plain NDNR retrievals. It was also shown that retrievals using a combination of NIPA and NDNR are feasible with Landsat observations.
Standard IPA retrievals of optical depth, effective radius, and liquid water path for clouds were also performed for the Landsat scene. The successes and difficulties of the IPA method, which is the standard assumption used in all operational cloud retrievals, were demonstrated and discussed. Liquid water path comparisons with the surface microwave radiometer at the ARM site gave statistical consistency when the northernmost 2000 Landsat scanlines were used and retrievals smaller than 20 g m−2 were eliminated from both datasets. Comparisons of modeled surface solar fluxes with pyranometer measurements at the Central Facility were of somewhat dubious value because of the presence of broken clouds at the time of the satellite overpass and the ensuing large fluctuations in the measured irradiance. IPA encountered problems matching modeled and observed B7 (2.2 μm) reflectivity values for pixels near cloud edges.
The retrievals were also sensitive to the choice of conservative-absorbing band pair: pairs that include B5 (1.6 μm) gave less reliable effective radii than pairs with B7, while B4 (0.83 μm) retrieved larger optical depths than (the more sensitive to saturation) B2 (0.55 μm). The means and histograms of retrieved optical properties for the entire scene were rather insensitive to the degradation of resolution indicating that either the cloud fields did not have inherent variability up to ∼2-km scales or that radiative smoothing gave them apparent smoothness. Last, cloud amounts determined from a global operational algorithm to be used with Landsat-7 compared well with amounts retrieved from a simple (scene-specific) threshold method.
Clearly, NIPA, NDNR, and a potential hybrid scheme that would remove both high-order (smoothing) and low-order (roughening) 3D scattering signatures are concepts still in their formative stages. The dependencies of the NIPA parameters on cloud properties and conditions are not yet thoroughly known, and the procedure for tracking down their optimal values is still under development. The suitability of NIPA-type extensions for broken clouds and NDNR difficulties for optically thin clouds remain subjects of further study. The selection of the most appropriate absorbing channel for NDNR (for instruments allowing different choices) or even the search for other normalized indices capable of removing 3D effects are also topics that merit further investigation. We are currently performing extensive MC simulations to answer these questions. While the simulations have the advantage of providing the true answer, they are computationally expensive, and force some compromises/approximations in the simulated cloud fields that may prolong the search for definitive answers. Nevertheless, NIPA and NDNR provide a solid framework on which to continue building techniques that account for 3D effects in optical property retrievals.
Data were obtained from the Atmospheric Radiation Measurement (ARM) Program sponsored by the U.S. Department of Energy, Office of Energy Research, Office of Health and Environmental Research, Environmental Sciences Division. We thank Richard Irish for providing us information about the ACCA algorithm, and W. Wiscombe, R. Pincus, A. Davis for helpful discussions. Thanks also go to M. Mischenko of NASA GISS for making the ISCCP ice crystal phase function available to us. This research was supported by funding provided under Landsat Science Team activities, which are part of NASA-MTPE, under Proposal 1996-MTPE-00116, the Department of Energy’s ARM program under Grants DE-A102-97ER62369, and DE-A105-90ER61069, and NASA Grant NAG5-6675 as part of the EOS validation program.
Corresponding author address: Lazaros Oreopoulos, NASA GSFC, Code 913, Greenbelt, MD 20771.