Data Collection 5 processing for the Moderate Resolution Imaging Spectroradiometer (MODIS) on board the NASA Earth Observing System (EOS) Terra and Aqua spacecraft includes an algorithm for detecting multilayered clouds in daytime. The main objective of this algorithm is to detect multilayered cloud scenes, specifically optically thin ice cloud overlying a lower-level water cloud, that present difficulties for retrieving cloud effective radius using single-layer plane-parallel cloud models. The algorithm uses the MODIS 0.94-μm water vapor band along with CO2 bands to obtain two above-cloud precipitable water retrievals, the difference of which, in conjunction with additional tests, provides a map of where multilayered clouds might potentially exist. The presence of a multilayered cloud results in a large difference in retrievals of above-cloud properties between the CO2 and the 0.94-μm methods. In this paper the MODIS multilayered cloud algorithm is described, results of using the algorithm over example scenes are shown, and global statistics for multilayered clouds as observed by MODIS are discussed. A theoretical study of the algorithm behavior for simulated multilayered clouds is also given. Results are compared to two other comparable passive imager methods. A set of standard cloudy atmospheric profiles developed during the course of this investigation is also presented. The results lead to the conclusion that the MODIS multilayer cloud detection algorithm has some skill in identifying multilayered clouds with different thermodynamic phases.
Plane-parallel single-layered cloud radiative transfer (RT) models are used by global passive imager algorithms like Moderate Resolution Imaging Spectroradiometer (MODIS) (Barnes et al. 1998) for cloud thermodynamic phase, cloud-top pressure–temperature, and optical and microphysical properties retrievals (King et al. 2003; Platnick et al. 2003). The use of such an RT model works reasonably well, as confirmed by many field campaigns and theoretical calculations (King et al. 2004; Mace et al. 2005; Chiriaco et al. 2007; Bedka et al. 2007; Otkin and Greenwald 2008). The model can work for some retrievals if there are multilayered clouds in a vertical column (e.g., an ice cloud overlapping a liquid water cloud) and the uppermost layer is optically thick. In particular, use of the RT model can result in biases with cloud effective radius retrievals when liquid water clouds are overlaid by relatively thin cirrus clouds (Davis et al. 2009). The retrieved effective radius of what is thought to be single-layer ice clouds decreases significantly in areas overlying the water clouds. When the cirrus is too optically thin to dominate the upwelling radiance and the cloud is identified as being liquid water phase, the retrieval tends toward abnormally large water droplets. There is not a large detrimental effect on cloud optical thickness to the extent that the combined optical thickness of all layers is retrieved with little dependence on the assumed phase.
It is important to flag areas where there are problematic effective radius retrievals due to multilayer clouds of differing thermodynamic phases since those retrievals can adversely affect cloud statistics and should be excluded from further analysis.
There have been other algorithms designed to identify multilayer clouds with passive imagers. The algorithm developed by Pavolonis and Heidinger (2004) is a pixel-level algorithm that uses ratios and differences of reflectances and brightness temperatures in various bands. This approach can be applied to historical and current multispectral imager data such as the Advanced Very High Resolution Radiometer (AVHRR) on the National Oceanic and Atmospheric Administration (NOAA) spacecraft and MODIS. Such an approach may also be continued with future measurements from the Visible/Infrared Imager Radiometer Suite (VIIRS) that will be flown on the National Polar-Orbiting Operational Environmental Satellite System (NPOESS) platforms. This algorithm is used in NOAA’s Extended Clouds from AVHRR (CLAVR-x) processing system (Heidinger and Pavolonis 2005). We had not chosen to use this method for MODIS Collection 5 as its skill specifically in detection of multilayer clouds that may compromise cloud effective radius retrievals had not been examined. We are currently performing such examination for MODIS Collection 6 and will show a comparison of this AVHRR–VIIRS algorithm with the MODIS algorithm in section 5.
The algorithm developed by Baum and Nasiri (Baum et al. 2000; Nasiri and Baum 2004) is a statistically based algorithm that is executed in shifting steps over a box area of user-defined size, typically 200 × 200 pixels, with a restriction that some clear sky is available in the area; the algorithm retrieves a probability that the cloud is multilayered. This algorithm was developed for the MODIS instrument, but has not been used extensively outside of case studies. The need to use a large area to work on and a requirement for the presence of clear-sky pixels within each work area reduces the effective algorithm resolution and usefulness as many multilayered cloud retrievals occur within synoptic systems that span a wide area with extensive cloud cover. We will show a comparison of our algorithm with the Nasiri–Baum algorithm in section 5.
Another approach for multilayer cloud detection has been presented by Chang and Li (Chang and Li 2005a,b). The method of Chang and Li uses an estimation of cirrus cloud emissivity based on the difference of cloud-top temperature retrieved by using the CO2 slicing result (not assuming an opaque cloud) and the 11-μm band (i.e., assuming an opaque cloud). The algorithm relies on being able to identify single-layer liquid water clouds and clear-sky pixels in an area of 250 km × 250 km centered on the point of interest. The cloud effective emissivity is then computed, from which the infrared (IR) cloud optical thickness is derived. If that cloud optical thickness is significantly different from the cloud optical thickness retrieved using a visible or shortwave infrared (SWIR) band, the cirrus cloud likely has a liquid water cloud underneath it. We did not use this method in Collection 5 because it is not a pixel-level method and thus would not function well within the infrastructure of our operational code. It also requires knowledge of IR cirrus emissivity. That quantity is produced by the MODIS cloud-top properties product (Menzel et al. 2008), but the uncertainty is not specified, especially in the presence of multilayer clouds of the type we study, making it difficult to assess algorithm skill. Depending on improvements made by the MODIS cloud-top properties algorithm team it may be possible to include this algorithm in final analysis in the future.
The MODIS operational multilayer cloud detection algorithm relies on a difference in above-cloud precipitable water retrievals obtained from using the 0.94-μm band versus above-cloud precipitable water computed from the CO2 slicing–derived cloud-top altitude. The 0.94-μm band is relatively insensitive to optically thin cirrus and so the column moisture is integrated from the top of the atmosphere (TOA) to the lower-level cloud, if such is present. The CO2 slicing retrieval of cloud-top height, and subsequent calculation of the above-cloud precipitable water from a forecast model profile, occurs from the TOA to the level of the higher cloud. From that difference, and several other tests such as the difference between retrieved IR and SWIR cloud thermodynamic phases and reflectance ratios to screen for single-layer clouds over bright surfaces, a determination is made as to whether the cloud is multilayered in a way that affects the applicability of the plane-parallel single-layer cloud models used in retrievals of cloud effective radius.
Initial CloudSat evaluations of the MODIS multilayer cloud detection algorithm have been done by Joiner et al. (2010) as part of a study that developed a global multilayer cloud detection algorithm via cloud-top pressure derived from the Aura Ozone Monitoring Instrument (OMI). Joiner et al. found good agreement between MODIS multilayer algorithm and CloudSat with 83.4% agreement at 5 km × 5 km area comparison and improving to 91% for 12-km OMI footprint area comparison. We present a detailed discussion of evaluation results in section 5b.
In the following discussion we present the MODIS operational multilayer cloud detection algorithm, describe how the multilayer cloud information is stored in the MODIS cloud optical and microphysical properties product (MOD06/MYD06) level-2 hierarchical data format (HDF) files, present results of executing the algorithm on data produced by forward simulations of multilayered clouds, and compare the algorithm to other methods.
To summarize briefly the discussion that will follow, section 2 describes the algorithm and the data format in which the results are stored. Section 3 presents the details of the RT simulations and describes in detail the method used to create the simulation dataset. Section 4 provides results from applying the MODIS operational multilayer cloud detection algorithm over the simulated scenes as well as selected MODIS data granules. Section 5 discusses a direct comparison of our results with other passive remote sensor methods for detecting multilayer clouds and presents initial validation with CloudSat/Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation platform (CALIPSO). We show global statistical aggregation of multilayer cloud data in section 6. Conclusions, ongoing work, and future directions are discussed in section 7.
2. Algorithm description
The operational MODIS multilayer cloud retrieval uses a number of bands in addition to individual retrievals of physical quantities, such as above-cloud precipitable water and cloud optical thickness, to arrive at a decision. The main component of the retrieval is a test for the difference of above-cloud precipitable water retrievals obtained by two different methods.
a. MODIS CO2 slicing
The first method is based on the cloud-top pressure retrieval obtained from CO2 slicing using ratios of MODIS channels centered at 13.3, 13.6, 13.8, and 14.2 μm (Menzel et al. 2008). The retrieved cloud-top pressure is then used to obtain above-cloud water vapor amount (PWCO2) by adding up the layer averaged water vapor amounts from the National Centers for Environmental Prediction (NCEP) global 6-h atmospheric profile product, produced at 1° resolution. Because of the nature of CO2 absorption, the algorithm is sensitive to high clouds of optical thickness (τc) greater than 0.5 (Menzel et al. 2008) when multilayer clouds are present and will return a low value of above-cloud precipitable water.
b. The 0.94-μm above-cloud precipitable water
The second method uses water vapor absorption in the MODIS 0.94-μm band. Above-cloud precipitable water is retrieved using an iterative approach. That is possible because cloud reflectance is flat in the spectral range between 0.86 and 0.94 μm and the difference in measured cloud reflectance is due to the water vapor amount between the cloud and the sensor. If the visible optical thickness of thin cirrus layer is less than 6, the 0.94-μm band is sensitive to the low clouds when multilayer clouds are present and will return a higher value of above-cloud precipitable water than the CO2 slicing method would. The discrepancy in retrieved amounts of above-cloud precipitable water can be attributed to the presence of multilayered clouds.
The MODIS operational multilayer algorithm first assumes that a single layered cloud exists, with a cloud-top temperature based on the 11-μm brightness temperature. Cloud-top pressure is then inferred by mapping the temperature into the NCEP pressure profile. The mapping is done from the top downward so as to avoid the high likelihood of temperature inversions nearer the surface.
This cloud-top pressure together with the view geometry is used to index a MODIS atmospheric transmittance table for 0.86 and 0.94 μm, which is generated by using the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) atmospheric profile database as input to moderate resolution atmospheric transmission (MODTRAN) version 4.2r1 (Berk et al. 1998). This lookup results in a vector of two-way atmospheric transmittance as a function of above-cloud precipitable water for each band. These transmittance vectors are then applied to the measured reflectances. The dominant contributor to absorption in the 0.94-μm band is water vapor. If there were no water vapor between observer and cloud, measured reflectances can be assumed to be identical. Using that assumption we look for a point where the two vectors intersect. The closest table index value of above-cloud precipitable water at the intersection point is our retrieval of above-cloud precipitable water (PW0.94). We choose to neglect a very small amount of ozone absorption in the 0.86-μm band (<0.001 additional absorption amount) as it has no discernible impact on the location of the intersection point because of lookup table resolution.
We then use the retrieved water vapor amount to perform a crude atmospheric emission correction on the 11-μm radiance. Measured 11-μm radiance consists of three components: emission from ground, emission from cloud and emission from atmosphere above cloud. We assume that cloud emissivity is unity; therefore, we do not deal with emission from ground. This is the exact assumption made by MODIS CO2 slicing–based cloud-top properties retrieval method. Now we must subtract the atmospheric emission from measurement and also correct the result for water vapor absorption in the 11-μm channel. So the final corrected radiance takes on the following form:
where Imeas is the measured radiance, Tmean_above_cloud is the integrated layer mean temperature from given atmospheric profile, and trans is the 1-way atmospheric transmittance at 11 μm.
The entire process is repeated using the corrected 11-μm radiance as a source of cloud-top temperature. We have found that one additional iteration is enough for the retrieval to converge to within 0.25 K, which we consider to be a sufficient degree of accuracy for our purpose. This same type of retrieval, but without iteration, is also performed one additional time with the assumption that the cloud in question is located at 900 hPa (PW0.94@900). If a high, cold cloud (Tc < 265 K) with little water vapor above it, is moved vertically in the atmosphere, its retrieved temperature and pressure stay nearly constant because atmospheric transmittance for amounts of water vapor less than 0.5 cm shows very little dependence on pressure. Moving such a cloud from 200 hPa down to 900 hPa changes 11-μm transmittance by only 0.8%, 0.94-μm transmittance by 1.05%, and 0.86 μm by 0.01%. However, this is not so for a warm, low cloud with a significant amount (>1 cm) of water vapor above it, which is fairly typical for boundary layer clouds. For such cloud, 11- and 0.86-μm transmittances change by about the same amount as for a high cloud, but the 0.94-μm transmittance changes by 8% if such cloud with 1 cm of above-cloud precipitable water above it is moved between 600 and 900 hPa. The error in retrieved above-cloud precipitable water amount for the lower-level cloud will increase as the optical thickness of the overlying ice cloud increases. The result is similar regardless of where the lower-level cloud lies between 800 and 1000 hPa. A low-level cloud pressure of 900 hPa is chosen as the default value. We mitigate the effect of ground elevation because the NCEP profiles extrapolate every profile down to 1000-hPa level, regardless of terrain. An above-cloud precipitable water retrieval based on this assumption acts to mitigate the “cooling” effect of an upper ice cloud and results in the inference of a more realistic high above-cloud precipitable water amount above the lower-level cloud. This process does not affect the results for single-layered ice clouds or multilayered clouds where the upper ice cloud layer is optically thick, and permits the tracking of more multilayered clouds.
As both 0.94- and 0.86-μm channels are much more sensitive to the presence of lower-level clouds in multilayer situations, the retrieved above-cloud precipitable water value is quite different from the same retrieval performed based on the inference of high clouds from CO2 slicing. That difference, weighted by the total column precipitable water (TPW), is a key determinant of whether or not there may be multilayered clouds present. A value of
is used as the threshold for marking the pixel as potentially containing multilayered clouds based on case studies and estimates regarding the occurrence of effective radius biases (see following example). Forward radiative transfer simulations, discussed in sections 3 and 4, confirm that this is an appropriate choice.
c. Example retrieval
Figure 1 illustrates the effect of this retrieval on a portion of a MODIS data granule. These data were collected from Terra MODIS at 0015 UTC 25 October 2008 in the western Pacific Ocean just east of Japan. The panels show the process of obtaining a multilayer result using the above-cloud precipitable water method. Figure 1a shows a false-color image of MODIS bands 6, 2, 26 (1.64, 0.86, and 1.38 μm, respectively). Thin cirrus is advecting over a field of cumuliform clouds. In this false-color composite, liquid water clouds appear gold, ice clouds appear blue and white, and the ocean surface appears black. A number of areas where thin cirrus overlaps the liquid water clouds are visible in the image and take on a greenish hue. Figure 1b is an image of above-cloud precipitable water from the MODIS cloud-top properties algorithm that uses CO2 slicing (PWCO2). The figure shows a low above-cloud water vapor amount retrieved throughout, indicating that the IR bands are seeing the high clouds and not the low ones. There is barely a trace of the low-level clouds in the image. Very low values of above-cloud precipitable water seen for the high clouds are as expected. Figure 1c is an image of the standard 0.86–0.94-μm retrieval of above-cloud precipitable water (PW0.94), which is more sensitive to low clouds and so gives higher above-cloud precipitable water values that are more typical for those clouds. Figure 1d is the difference image between the above-cloud precipitable water from the CO2 slicing and the 0.94-μm algorithm. Outlines of low-level clouds are becoming clearly visible in the difference image. Small differences in above-cloud precipitable water correspond to either thicker cirrus, which is not sensitive to multilayer clouds, or breaks in the low-level cumulus clouds. But more cloud could be flagged as multilayer as the cirrus becomes thicker to the west and affects the vertical placement of the cumulus by the 0.94-μm-based cloud-top properties retrieval. Figure 1e shows the above-cloud precipitable water retrieval in which the low-level clouds are assumed to be at the 900-hPa level. It is not that different from the main 0.86–0.94-μm result with the exception that it captures some of the cloud features covered by somewhat thicker cirrus to the west. Even though the clouds are thicker, they still contain some contribution from the underlying low-level cloud. Figure 1f shows the difference image resulting from the 900-hPa retrieval versus that from the CO2 slicing.
The final two images are the retrieved cloud optical thickness and effective radius for the scene. The warm colors indicate liquid water clouds with cold colors for ice cloud retrievals. The optical thickness image indicates that the cirrus is quite thin and fairly uniform over the overlap area. There is no significant impact of multilayered clouds on optical thickness as the overlying cirrus is thin and its contribution to the combined visible optical thickness is very small. In contrast, the impact on the cloud effective radius retrieval is much greater. The outlines of low-level clouds are clearly seen in the effective radius image as areas of small ice effective radii. The breaks of open water in the cumulus cloud fields return effective radius values of around 25 μm, so it is unlikely that the actual cloud microphysics is changing in the overlap area.
Figure 2 shows the net statistical effect of multilayer clouds on cloud optical thickness and cloud effective radius. While there is not a large effect on cloud optical thickness, there is a significant shift in effective radius distribution toward smaller radii when multilayered clouds are not removed from the scene. In this particular case, 19.2% of ice clouds in the scene were multilayer and ∼54 000 pixels were removed from the distribution.
The MODIS CO2 slicing algorithm is applied with the most confidence for clouds at pressures lower than about 700 hPa (Menzel et al. 2008). In a typical MODIS scene, however, the CO2 slicing algorithm is rarely applied for clouds at pressures larger than 600 hPa. If the CO2 slicing algorithm is unable to converge on a solution, the 11-μm band is used under the assumption that there is a low-level opaque cloud present. The choice was made to ignore CO2 slicing results at pressures larger than 550 hPa to minimize the potential for false positive retrievals. In light of improvement in vertical resolution to 101 levels used in MODIS CO2 slicing algorithm beginning with Collection 6, this 550-hPa restriction may be eased in the future, although uncertainties due to resolution of the NCEP profiles will remain.
d. Possible limitations
Because of uncertainties in inferring cloud emissivity from passive sensors, it is possible to obtain a false positive multilayer retrieval for the case when an optically thin cirrus cloud is present with τc < 4. If the cloud is very optically thin, upwelling radiance from surface will cause that cloud to be placed at pressure much higher than truth. That means the 0.94-μm cloud-top properties method will retrieve much higher above-cloud precipitable water amount than CO2 slicing would because of surface contamination and not because of multilayer situation. We assume that if the total column optical thickness is <4, the likelihood is that there is not a lower cloud underneath it. The liquid water cloud layer underneath would most likely push the total optical thickness above 4. If a liquid water cloud is so thin that threshold of 4 is not reached, then we would have difficulty with retrieving effective radius because of the shape of the forward library space (Platnick et al. 2003), any multilayer situation aside. False negatives do arise from use of this threshold, but with overall effective radius retrieval uncertainty being well above 20% for thin clouds, the weight of such retrievals should be greatly reduced in any statistical studies anyhow.
We also must consider cases of single-layer clouds over bright surfaces. It is possible for the algorithm to mistake a thin cirrus cloud over a bright surface for a cloud that is multilayer. The 0.65- and 1.24-μm reflectances are used to check for vegetation and snow–ice, respectively. Cloud reflectance is reasonably flat in that spectral region, while surface albedo changes significantly. So for a true multilayer cloud situation, the reflectance ratio would be close to 1.0, but not so for a single layer of thin cirrus over a bright surface. It is useful to use ratios of 0.86-μm reflectance to 0.65- and 1.24-μm reflectance to check for bright surfaces, with thresholds set as follows:
These thresholds were empirically derived on the basis of case studies; however, our forward simulations indicate that a parameterization based on ecosystem type may be more appropriate in the future. We will investigate such parameterization in MODIS data for Collection 6.
e. Cloud thermodynamic phase test
In addition to the above-cloud precipitable water difference, another test is based on retrievals of cloud thermodynamic phase from two different methods. The first method is the MODIS SWIR thermodynamic phase (SP) algorithm (Platnick et al. 2003) that uses a number of cloud mask tests and reflectance ratios in visible, near-infrared (NIR), and SWIR bands to arrive at cloud thermodynamic phase. The second method is the IR bispectral cloud phase (IP) algorithm based on brightness temperature differences between the 8.5- and 11-μm bands, which is a modification of the Baum IR trispectral algorithm (Baum et al. 2000). When these two methods infer different thermodynamic phases, that can be an indication of a multilayered cloud situation. This particular test tends to be sensitive to cirrus over liquid water clouds in which thin cirrus is too thin to result in an ice phase retrieval but still biases the liquid water cloud retrievals as the cloud effective radius retrieval is larger than expected.
The main uncertainty associated with using the thermodynamic phase test tends to arise in polar regions. At latitudes above 60°, the IR method results in quite a few undetermined phase answers because of inherent difficulties of an IR method over very cold surfaces, so we assign a lower degree of confidence to multilayered clouds that are flagged only by the cloud phase test and no other test.
f. Stored product
The 0.94-μm above-cloud precipitable water retrieval performed at both pressure at cloud top and at 900 hPa, together with a test on retrieved cloud thermodynamic phase, combine to create a final integer answer that tells the user whether the multilayer detection algorithm arrived at a positive result and what method(s) were positive as shown in Table 1. We store the final value in the MOD06/MYD06 level 2 HDF file in two places. The values from Table 1 are stored in a scientific dataset (SDS) named Cloud_Multi_Layer_Flag. The multilayer cloud information is also stored in the fifth byte of the Quality_Assurance_1km SDS as information about the thermodynamic phase of the cloud and its multilayer status. The full description of the Quality_Assurance_1km SDS is given in Hubanks (2006) and a brief listing of relevant values is given in Table 2.
The discussion in this section is summarized in Fig. 3. The algorithm flowchart shows the overall logical flow of the algorithm.
3. Radiative transfer models
We have conducted an extensive set of forward RT modeling studies of multilayer clouds under varying atmospheric conditions, layer separations, surface types, and layer thicknesses to thoroughly test the sensitivities and skill of the MODIS multilayer cloud detection algorithm.
Zonal and temporal average profiles are calculated from the ECMWF sampled 60-level global atmospheric profile database aggregated from ERA-40 data over 48 days for two years using 1 and 15 of each month between January 1992 and December of 1993 (Chevallier 2001). The database profiles were separated to represent a typical midlatitude summer (MLS), midlatitude winter (MLW), tropical atmosphere (TRP), and polar oceanic (POL) profile. Profiles over polar landmasses, dominated by profiles from the Antarctic continent, were not included as they would contain strong inversions and would be likely used disproportionately for pressures lower than 700 hPa. The polar oceanic profile consists of daytime profiles only. Nighttime profiles are not used since for our purposes, cloud optical and microphysical property retrievals are performed in daytime only. We define the tropical region as 30°S < latitude < 30°N, midlatitudes as 30° < |latitude| < 60° and the polar regions as above 60° latitude. For midlatitudes, winter profiles occur between 1 November and 30 April; summer profiles are the remainder of the year. Within each latitude belt, profiles are chosen from regions that had cloud fraction (CF) > 0.85 to match the conditions of interest. Profiles were separated further by land and ocean using the ECMWF land fraction flag with threshold set at 0.5.
Starting with these averaged profiles, we set relative humidity to 100% for the specific levels that contained clouds in our simulations. The profiles were interpolated from the native 60-level resolution to 36 levels spaced at 1 km vertically between 0 and 25 km with sparser resolution in the upper atmosphere. These profiles are available at http://modis-atmos.gsfc.nasa.gov/MOD06_L2/validation.html.
Figure 4 shows a combined plot of the temperature and moisture profiles used in the simulations. These particular plots show the liquid water cloud layer at 2 km. Simulations were run for a variety of solar and view zenith angles with the solar zenith angles appropriate for the time of year in question. We sampled the solar zenith angle from the MODIS level 3 global monthly product (Hubanks et al. 2008). The cosine of the view zenith angle corresponded to μ = 1.0, 0.8, and 0.6. For detailed examination, simulations were run for ice cloud effective radii of 10, 30, and 50 μm and water radii of 6, 10, and 20 μm. An ice cloud layer of 2-km physical thickness was fixed at the base of the tropopause as indicated by temperature in each of the different profiles shown in Fig. 4: 8 km (MLW and POL), 12 km (MLS), and 14 km (TRP). The ice cloud optical thickness varied between 0 and 20, with increments as follows: 0.0, 0.1, 0.25, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0, and 20.0. These increments were chosen specifically to examine the thin cirrus region and also to appropriately capture the point where the ice cloud becomes too thick to be affected by the underlying liquid water cloud. Water cloud layers were assumed to be 1 km thick and were placed at two different altitudes: 2 and 4 km. For liquid water clouds, optical thicknesses ranged as follows: 0.0, 2.0, 5.0, 10.0, and 20.0.
Radiances were simulated for 16 MODIS bands, which was necessary to perform the relevant cloud mask tests (Ackerman et al. 2006; Frey et al. 2008), in particular the 3.7–11-μm brightness temperature test, the CO2 slicing cloud-top properties retrieval, and the full MODIS cloud optical and microphysical property retrievals. The set included the 0.65-, 0.86-, 0.94-, 1.24-, 1.38-, 1.6-, 2.1-, 3.7-, 3.9-, 8.5-, 11.0-, 12.0-, 13.3-, 13.6-, 13.8-, and 14.2-μm MODIS channels (Ackerman et al. 2006).
Each simulation was repeated over a wide variety of surfaces. The oceanic profiles only had one option (dark ocean with surface albedo of 0.05) with the exception of polar ocean that also included a sea ice surface. The land surface profiles presented options of vegetated, desert, or snow cover. Midlatitude land included mixed forest and desert with or without snow, appropriately, while tropical land included desert and evergreen broadleaf forest. All classifications were based on definitions of the International Geosphere–Biosphere Programme (IGBP) ecosystem map and the surface albedo values taken from MODIS albedo product (MOD43)-based 1-km-resolution surface albedo product (Moody et al. 2005, 2007, 2008). Figure 5 shows a plot of the white-sky (diffuse) surface albedo as a function of wavelength for the various surfaces considered in this investigation. MODIS bands that contain no solar component were given a zero surface albedo.
The RT simulations were performed using Discrete Ordinate Radiative Transfer (DISORT) code (Stamnes et al. 1988) using liquid water cloud phase function results from Mie calculations based on the water droplet size distributions using a gamma distribution with an effective variance of 0.1 (Platnick et al. 2003) and bulk ice cloud phase functions developed by Baum et al. (2005a,b). The same phase functions for both ice and liquid water are used in the lookup tables (LUTs) employed in the MODIS cloud optical and microphysical properties algorithm for Collection 5. The correlated-k method (Kratz 1995) was used to account for water vapor and other gaseous absorbers. The DISORT code, in conjunction with the correlated-k method, then produced the simulated MODIS band radiances. We used 32 streams in our radiative transfer calculations, which, together with truncation of strong forward peaks and use of delta-fit method by Hu et al. (2000), can be considered sufficient computational accuracy as described by Ding et al. (2009).
With the parameter ranges described above, the forward RT calculations resulted in 26 files corresponding to combinations of atmospheric profiles and surface types. Each file contained 7560 individual data points for each geometry, optical thickness, and effective radius tested. Results are provided in the following section for application of the MODIS multilayer cloud detection algorithm to a cross section of this database of simulated MODIS radiances. In section 5 similar results are provided for the Pavolonis–Heidinger and Nasiri–Baum algorithms to this same dataset with comparison to the results from the MODIS operational algorithm.
4. Model retrieval results
In this section we show results of applying the MODIS multilayer cloud retrieval simulated MODIS data.
We could not run the Nasiri–Baum algorithm on DISORT simulations because it is a statistical aggregate algorithm that depends on natural variability of the data within a given area.
We show the results from a cross section of our forward RT simulations. Figure 6 shows a set of combined results from the DISORT forward simulations. To facilitate the interpretation of results, we group individual runs having all but one identical parameter to illustrate the effect of the differing parameter on the multilayer cloud retrieval result. Figure 6a combines the results of simulations conducted with a nadir view, solar zenith at 32°, dark ocean surface, and liquid water cloud located at an altitude of 2 km. The atmospheric profile is varied in terms of the overall column moisture content. The plot in Fig. 6a effectively shows multilayer cloud detection as a function of the total column water vapor. The “bits” in the effective binary numbers that result from this data combination indicate whether or not a multilayer cloud was detected. The bit significance was arranged as a function of the column moisture with the least significant bit for the most moisture. For example, a value of 011, which is light green in the plot, means that a multilayer cloud was detected under the conditions specified above using TRP and MLS profiles, but no multilayer cloud was detected for the MLW profile. The algorithm is more likely to detect a multilayer situation when the ice cloud is optically thin if the atmospheric moisture content is higher.
Figure 6b shows the same basic situation as Fig. 6a with the exception that the altitude of the lower-layer liquid water cloud was placed at 4 km and thus decreases the cloud-layer separation. When the cloud-layer separation is smaller, the amount of atmospheric water vapor between the two cloud layers is also lower and so the absorption in the 0.94-μm channel is decreased over the previous case. The sensitivity of the algorithm decreases as the ice optical thickness increases compared to the case where the liquid water cloud is at 2-km altitude. Some false positives occur in which multilayer cloud is detected for thicker liquid water clouds where there is no ice cloud above. These false positives come from growing uncertainties in retrieving IR cloud phase and CO2 cloud-top properties as the cloud gets colder. The detection results can be inspected further by looking at individual tests, some of which have lower confidence than others as mentioned in section 2. The detection status is reported as a binary answer and may result in a false positive result.
Figure 6c illustrates the multilayer detection result as a function of underlying surface type, assuming a single MLS profile and a liquid water cloud placed at 2-km altitude. The surface types are arranged such that the least significant bit corresponds to the lowest overall surface albedo with no snow on the ground. The plot shows that multilayered clouds are not detected for a desert ecosystem with thin liquid water clouds below, since the liquid cloud emissivity is likely somewhat less than 1.0, thereby indicating that we may need a separate detection threshold for deserts since the surface albedo of deserts is significantly different in spectral shape from vegetation and snow–ice surfaces. The desert spectral albedo tends to be somewhat flatter than vegetation, as Fig. 5 shows, and so may require a somewhat different approach. The effect of this on our global statistics is not very significant as the actual cloud fraction over deserts is rather low (cf. Fig. 9). Figure 6d shows multilayered cloud detection as a function of cosine of the viewing zenith angle (μ) for a MLS profile with a dark ocean surface and a lower-layer liquid water cloud placed at 2 km. The points are ordered in μ space such that a more oblique angle, that is, lower μ, is the least significant bit in the binary number displayed. The relative azimuth angle for this comparison was set to 0°. The figure indicates that the algorithm is more likely to detect a thinner ice cloud over a liquid water cloud at more oblique angles. On the other hand, it is possible to flag cases with higher ice cloud optical thicknesses at more nadir view angles.
The Pavolonis–Heidinger method, originally developed for the AVHRR and adapted for the upcoming VIIRS instrument, uses a series of reflectance and brightness temperature difference thresholds described in detail in (Pavolonis and Heidinger 2004). For the algorithm comparison purposes we have been provided with their most recent development of the method, with improvements and modifications made since the publication of their paper. Similar to the MODIS operational algorithm, it is a single-pixel method that works on samples individually without using any spatial aggregation. Because of this similarity we were able to execute the Pavolonis–Heidinger algorithm on the results of our DISORT simulations of multilayer clouds. Figure 7 shows the comparison of these results. The figures on the left are the MODIS results from Figs. 6a,b,d, and on the right are corresponding Pavolonis–Heidinger results. We compared the algorithms for three out of four database cross sections shown in Fig. 6. It was not possible to perform the exact comparison for the surface-type section, since the Pavolonis–Heidinger algorithm uses a 0.41-μm band over desert regions that we did not include in the original DISORT band set. The Pavolonis–Heidinger algorithm uses a LUT derived from simulations of multilayered clouds over a various surfaces. The LUT includes the difference in brightness temperatures (BTD) between the 11- and 12-μm bands. A threshold function is defined since the multilayered clouds (i.e., ice over water cloud) display a BTD as a function of visible reflectance that is quite different from single-layered liquid water and ice clouds. In addition to that threshold, a number of constraints are placed on reflectances at 0.65 and 1.38 μm to help with the identification of single-layer clouds over a variety of surfaces. The 1.65-μm band is used by the algorithm to aid in identifying the thermodynamic phase of clouds since ice clouds have greater absorption than liquid water clouds at 1.65 μm.
There are similarities in the results as well as some differences, but overall the comparison is favorable. The MODIS algorithm has a somewhat wider section where multilayer clouds are detected for the entire range of the varied conditions, be it atmospheric moisture content or view angle. However, the detection rate generally drops off as the ice cloud thickens, with only the thickest simulated liquid water cloud showing at ice cloud optical thickness of 10. Both algorithms show that once the ice cloud optical thickness reaches 20, no detection of multilayer clouds is possible. Both algorithms also show that detection is a function of layer separation with detection rate being lower when the liquid water cloud is placed at 4-km as opposed to 2-km cloud-top altitude. The Pavolonis–Heidinger algorithm shows more detection when both cloud layers thicken, but not as much when the cloud layers are thin.
Our overall conclusion from examining all these results is that the MODIS multilayer cloud detection algorithm is robust and performs as intended under a wide variety of conditions.
5. Comparisons with other retrieval methods
In this section we show an example case study from MODIS and comparisons of our method against two other multilayer cloud detection algorithms, which we mentioned in section 1.
a. Case study comparison with different retrieval methods
Figure 8 shows an example of multilayer cloud detection for a Terra MODIS granule acquired at 0015 UTC 25 October 2008 off the coast of Japan. This is a full granule, a portion of which was shown in Fig. 1. Figure 8a shows an atmospherically corrected true-color image formed as a composite of MODIS bands 1, 4, and 3 (0.65, 0.55, and 0.47 μm, respectively). While the true-color image indicates where the clouds are, it provides very little information about the various cloud layers or the thermodynamic phase of the clouds. Figure 8b shows a false-color image formed as a composite of bands 6, 2, and 26 (1.64, 0.86, and 1.38 μm, respectively), which more readily separates clouds of different thermodynamic phase by color. There is a significant amount of multilayer cloud in this scene, indicated by areas where the yellow liquid water clouds show through the more blue and white ice clouds. Figure 8c shows the results from applying the multilayer cloud detection algorithm. Different values on the color scale correspond to tests flagging the cloud as clear-sky (0), single-layer cloud (1), and multilayer (2–8) cloud, as described in Table 1. These results are not an absolute measure of multilayer cloud amount, but rather provide a map of areas where the presence of multilayer clouds adversely affects cloud effective radius retrievals.
The three multilayer cloud detection algorithms previously discussed are now applied to the MODIS granule shown in Fig. 8, with the results shown in Fig. 9. Figure 9a shows the true-color composite constructed from bands at 0.65, 0.55, and 0.47 μm; Fig. 9b shows the false-color composite constructed from bands at 1.64, 0.86, and 1.38 μm; and Fig. 9c shows the false-color composite constructed from bands at 0.55, 1.64, and 2.13 μm. Figures 9d–f show the results of applying the multilayer cloud detection using the MODIS operational algorithm (Fig. 9d), Pavolonis–Heidinger algorithm (Fig. 9e), and Nasiri–Baum algorithm (Fig. 9f).
As there is a wide range of options described in the code documentation that the Nasiri–Baum algorithm can be executed under, we took for the purposes of this comparison, the suggested default values. The Nasiri–Baum algorithm can only be executed under conditions that some clear sky, liquid water cloud, and ice cloud exists within the box being currently analyzed, so the algorithm does not attempt retrievals over a portion of this granule. The Nasiri–Baum algorithm also outputs its result as a probability of multilayer cloud being present. For clarity we display only a nonzero overlap probability as a positive answer. We performed a similar procedure with the results from the MODIS operational multilayer cloud algorithm, combining the multilayer values 2–8 into a single positive identification value. The Pavolonis–Heidinger algorithm returns its result as a single value so no additional data conversion was necessary to visualize the results.
Overall many of the same areas flagged as multilayer, even though the results may not look exactly the same, as the different multilayer algorithms were developed with different applications in mind. The Nasiri–Baum algorithm gives the fewest multilayer occurrences, but that can be attributed to a limited area over which the algorithm attempts retrievals. The main disagreement between Pavolonis–Heidinger and our algorithm arises in the flagging of thicker high clouds as being part of multilayer scenes (e.g., left side portion of the granule). One might argue that most, if not all, ice phase clouds in that part of the granule are multilayered clouds because of the apparent wide presence of low clouds in that region as well as there being some indication in the 1.38-μm false-color composite. The result given by the Pavolonis–Heidinger algorithm is consistent with detection achieved for simulated DISORT data, where clouds with combined extinction optical thickness as large as 30, with the upper-layer thickness of 10, can be flagged as multilayer, as shown in Fig. 7. However, that result is not directly useful for our purpose as when upper-layer thicknesses are so high, there is no effect on cloud effective radius retrievals and so the Pavolonis–Heidinger result is not optimal for detecting multilayer clouds that compromise effective radius retrievals.
The decision whether to flag a cloud as multilayer depends on the issue being addressed. In our case, we are looking for multilayer clouds that challenge the applicability of our single-layer plane-parallel cloud models used in cloud optical and microphysical property retrievals. Our goal is to create a map of areas where the model application is problematic. From our RT simulations we have found that the effect of ice cloud overlapping a liquid water cloud on cloud effective radius retrieval diminishes quite rapidly with increasing ice cloud optical thickness and is barely detectable when ice cloud optical thickness becomes greater than about 6. While it may be the case that the thicker upper-level clouds in this granule are also multilayer, having those clouds flagged as such does not impact our primary objective regarding microphysical biases.
b. Comparison with CloudSat/CALIPSO
Joiner et al. (2010) have provided initial comparisons of the MODIS multilayer cloud detection algorithm to CloudSat/CALIPSO as part of a study which developed a multilayer detection algorithm via OMI-derived cloud-top pressure (UV rotational-Raman scattering) compared with MODIS thermal emission retrievals. The OMI approach is philosophically similar to the MODIS approach reported here in that both take advantage of solar reflectance pathlengths being affected by gas species between cloud layers (O2 for OMI; water vapor in this study) in a fundamentally different way than CO2 slicing spectral bands; an O2 method, of course, has the advantage of using a well-mixed gas (see Figs. 6a,b). With multilayer detection from CloudSat defined as layer separations greater than 200 hPa, global analyses from a single day in Joiner et al. found that the MODIS algorithm correctly identified single- and multilayer cloud layers 83.4% of the time, with false positives and negatives, 9.8% and 6.8% of the time, respectively when used on a 5 km × 5 km sample. However, the MODIS algorithm is sensitive to upper-layer cirrus optical thicknesses on the order of 0.2 (Figs. 6, 7). False positives from missed thin cirrus are not considered in the Joiner et al. study. As already noted, the MODIS algorithm was intended to flag cases where upper-layer clouds have optical thicknesses too small to significantly screen lower-level clouds in the shortwave infrared bands used for particle size retrievals; such cases were not considered in assessing false negatives. For both types of false detection, information from CALIPSO cirrus optical thickness retrievals are needed, as well as monthly and seasonal comparisons (versus single day). The agreement with CloudSat improves to 91% detection with 3.7% and 5.3% false positives and negatives when applied over OMI footprint of 12 km. The sample was considered multilayered if any pixels within the sample being compared to CloudSat had the multilayer flag set. This result arises from the use of larger effective area to compare with CloudSat and statistics help smooth over the fact that the CloudSat footprint does not align with and is not the same size as MODIS. Multilayer cloud algorithm gives a binary answer and it is difficult thus to directly compare when resolutions and locations do not exactly align as many multilayered situations tend to occur near cloud feature edges. Thus much more work is needed to have an agreeable collocation scheme and to create a reasonable comparison basis between MODIS and CloudSat/CALIPSO.
Joiner et al. conclude that in addition to compromising the cloud effective radius retrievals, multilayer clouds must be identified for calculations of cloud radiative forcing as well. As the difference in cloud radiative forcing between a high and low cloud can be on the order of tens of W m−2, it is very important to use the cloud-top pressure that would be appropriate for shortwave view, that is, cloud-top pressure that does not include the thin cirrus that may be on top of the liquid water cloud layer.
6. Global results
In Collection 5, MODIS multilayer cloud retrievals are aggregated to the global level 3 daily, eight-day, and monthly products as an average of data down-sampled to 5 km and combined into a 1° grid. The multilayer cloud fraction is computed as fraction of cloudy pixels within a single grid box for which the multilayer cloud flag is set. This fraction is calculated using either all clouds or clouds of a particular thermodynamic phase as a basis. So there is a combined multilayer cloud fraction that divides the count of all multilayer cloud pixels by the number of all cloudy pixels. There is an ice multilayer cloud fraction that divides the count of pixels of ice thermodynamic phase for which multilayer cloud flag had been set over the count of all ice pixels in the grid box. Same operation is performed for liquid water and undetermined phase clouds.
The multilayer cloud information is also used to produce mean values of cloud optical and microphysical properties retrievals, both with and without multilayer clouds (Hubanks et al. 2008). Figure 10 shows an example of such an aggregation for the month of October 2008 derived from Terra MODIS data. Figure 10a shows the fraction of all cloudy pixels that have the multilayer flag set, and Fig. 10b shows the mean monthly cloud fraction. The small black area on the very top of the images corresponds to polar darkness or low sun where no retrievals are attempted (cosine of the solar zenith angle μ0 < 0.15).
A monthly global map like this is useful for providing the spatial distribution of multilayered clouds. Based on observational evidence, one might expect a higher frequency of multilayered clouds to occur in the vicinity of low pressure systems and their frontal boundaries. Higher frequencies of multilayer clouds tend to occur in the Southern Ocean and in the storm tracks (higher latitude zones) of both hemispheres. Tropical anvil cirrus is also a likely candidate to create multilayer cloud situations. A good portion of the intertropical convergence zone (ITCZ) is flagged as multilayer in the eastern Pacific Ocean. Strong convective zones over rain forest areas also tend to generate anvil cirrus, resulting in high frequencies of multilayered clouds in the Congo basin, Borneo, and New Guinea. One can also note the effect of advection of anvil cirrus over the marine stratocumulus zones in the Southern Hemisphere off the coasts of Peru and Ecuador, and in the Gulf of Guinea.
These global maps are also a useful measure of where there is a potential for higher uncertainty in retrieved cloud effective radius because of multilayer clouds.
7. Conclusions and future directions
In this paper we present the MODIS operational multilayer cloud detection algorithm used in the MODIS Collection 5 cloud optical and microphysical properties product (MOD06 for Terra; MYD06 for Aqua). The multilayer cloud detection method was developed to address a need to indicate areas of cloud where an assumption of single-layer plane-parallel cloud models was challenged because of the presence of two distinct cloud layers with differing thermodynamic phases, with the upper cloud layer being optically thin. Such situations manifest themselves as areas of abnormal cloud effective radius retrievals. Our method uses the difference between retrieved above-cloud precipitable water amounts from the 0.94-μm band and from the CO2 slicing cloud-top height together with a number of other tests. The physical basis of the multilayered cloud detection algorithm is provided, with examples of results from forward simulations as well as case studies involving MODIS data and global aggregations of results. Results from this approach are compared to two other methods of multilayer cloud detection. We also present a set of standard cloudy atmospheres that we developed to perform our studies. Wherever possible we perform all comparisons using a single source dataset, so the differences in retrieved results are solely due to differences in methodology.
Our results and analysis indicate that the multilayer cloud detection algorithm presents a reliable means of identifying situations that would create difficulties for retrievals of cloud effective radius. The forward simulations indicate that there are very few false positive results and that they arise under conditions that would result in high retrieval uncertainty due to one of the cloud layers being extremely thin. Forward radiative transfer simulations, performed under a wide variety of surface and atmospheric conditions, are used in our analyses to provide further insight as to the robustness of the algorithm.
We are currently investigating a number of improvements for the MODIS operational multilayer cloud detection algorithm that may be implemented for MODIS Collection 6. Those improvements involve bringing in additional retrievals of physical quantities performed using different methods, which in our experience has shown to contain multilayer cloud information. We also intend to continue our ongoing comparisons by performing more extensive MODIS to CALIPSO comparisons. However, for that work we require more data products than what is currently available from CALIPSO. We have begun such studies (R. Holz 2009, personal communication) and are awaiting the next release of CALIPSO products (version 3).
The authors thank Brad Wind for developing the groundwork for simplifying the modifications to the operational MODIS code that made most of these studies possible. This work was funded by the MODIS Science Team and NASA’s Radiation Sciences Program.
Corresponding author address: Galina Wind, Code 613.2, NASA Goddard Space Flight Center, Greenbelt, MD 20771. Email: firstname.lastname@example.org