Retrieval of Cirrus Cloud Optical Depth under Day and Night Conditions from MODIS Collection 6 Cloud Property Data

This paper presents a technique to generate cirrus optical depth and particle effective size estimates from the cloud emissivities at 8.5, 11 and 12 μm contained in the Collection-6 (C6) MYD06 cloud product. This technique employs the latest scattering models and scattering radiative transfer approximations to estimate cloud optical depth and particle effective size using efficient analytical formulae. Two scattering models are tested. The first is the same scattering model as that used in the C6 MYD06 solar reflectance products. The second model is an empirical model derived from radiometric consistency. Both models are shown to generate optical depths that compare well to those from constrained CALIPSO retrievals and MYD06. In terms of effective radius retrievals, the results from the radiometric empirical model agree more closely with MYD06 than those from the C6 model. This analysis is applied to AQUA/MODIS data collocated with CALIPSO/CALIOP during January 2010.


Introduction
The recently released MYD06 Collection 6 (C6) cloud products have been expanded in several ways compared to the Collection 5 (C5) data [1]. One of the new features of C6 MYD06 is the inclusion of cloud emissivity values in three InfraRed (IR) channels located in the longwave IR-window (8.5, 11 and 12 μm). The emissivities are a natural byproduct of the MYD06 cloud top height approach [2]. Band radiance thresholds in C6 MODIS are reduced, making high thin clouds detection more effective [1]. As shown in this paper, the addition of the emissivity values provides the basis for an accurate and computationally efficient estimation of the optical and microphysical properties of semitransparent cirrus clouds. These IR products can be generated both day and night and complement the suite of optical and microphysical properties in MYD06 derived during the day.
This paper describes the methodology and provides an evaluation through comparisons with the V3 Cloud-Aerosol LIdar with Orthogonal Polarization (CALIOP), the Infrared Imaging Radiometer (IIR) flown on CALIPSO, and the C6 solar reflectance based products in MYD06 [3]. The goal here is to make estimates of optical depth and particle size for cirrus clouds that are accurate and spectrally consistent, meaning that the properties are consistent from both the IR Longwave Window (8-12 μm) and the visible/shortwave-infrared retrievals from MYD06.
The technique described here is physically similar to that described by Garnier et al. [4,5], which uses data from CALIOP and the IIR. The Garnier technique uses information from CALIOP during the retrieval process and therefore the results are limited spatially to the narrow (~80 m) CALIOP/IIR swath. The technique described here does not require any information beyond that included in the C6 MYD06 data and that will likely be included in future MODIS Collections. This work builds on past studies that have explored the use of the channels for cirrus property estimation (e.g., [6][7][8]).

Motivation
As stated above, the methods developed here apply only to semitransparent ice or cirrus clouds. To motivate the importance of this subset of cloudiness, an analysis is done to determine the fraction of observations that would potentially yield results from this technique. To estimate this number, four filters are applied to the data. First the observations have to be cloudy. Second, the phase of the cloud has to be ice. Third, the cloud has to be semi-transparent, and fourth, the cloud's IR spectral emissivities vary within the expected range. Observations that survive these four filters represent the subset of cloudiness that could potentially yield results from the methods described later. Figure 1 shows the impact of these four filters on the MODIS/AQUA C6 observations for data (day + night) from January, 2010. Figure 1a shows the total cloud amount derived from the MODIS cloud mask. Figure 1b shows the distribution of ice clouds determined from the MODIS IR cloud phase product. The distribution of semi-transparent ice cloud is displayed in Figure 1c. Semi-transparency is determined using the 11 μm emissivity in the MODIS C6 product. This is defined here as ice clouds with the 11 μm emissivity ≤ 0.95, which corresponds to an absorption optical depth of 3. The technique discussed in this study is applicable to the clouds shown in Figure 1c. Finally, Figure 1d shows the distribution of semi-transparent ice clouds where the 12 μm emissivity is greater than the 11 μm emissivity. As described later, the microphysical models will yield physical results only when this condition holds. The mean global cloud fractions in order in Figure 1 are 69%, 28%, 19% and 14%. The final number, 14%, represents the potentially maximum yield of data for this month. The actual yield may be less, as consistency with other channels is required. These numbers indicate that roughly half of all ice clouds can be used in the retrieval described later. Another obvious factor is that MYD06 provides no cloud microphysical products at night and any information at night is valuable. Figure 1. Global distribution of mean cloud fractions for January, 2010 (day + night) from MODIS/AQUA MYD06 C6. Panel (a) shows the total cloud fraction; Panel (b) shows the ice cloud fraction; Panel (c) shows the semitransparent ice cloud fraction; Panel (d) shows the cloud fraction for clouds where the method described here can potentially retrieve an optical depth and a particle size. In terms of all pixels, the mean values are 69%, 28%, 19% and 14% for a-d. From these numbers, 50% of ice clouds produce viable optical depths and particle sizes.

Data Used in the Study
The data used in this analysis are from January 2010 and were provided by the University of Wisconsin NPP Atmosphere Product Evaluation and Test Element (PEATE). The MODIS level 2 products were generated at the PEATE using the most recent version of Collection 6 software, which includes the 8.5, 11 and 12 μm cloud emissivity values and an updated IR cloud phase [1]. This MODIS cloud phase is used later to isolate ice clouds. The updated IR phase approach employs cloud emissivity and cloud emissivity ratios and is described further in Section 4. The benefit to the use of cloud emissivity ratios is that the radiance contribution from the surface is accounted for through the use of a radiative transfer model. The signal from the cloud is enhanced by minimizing the influence of the surface and permits the algorithm to take advantage of the surface emissivity, which can depart significantly from unity over non-vegetated regions. The NASA Langley CALIPSO/CALIOP overpasses were collocated with MODIS using the procedure described in Nagle et al. [9]. From CALIPSO/CALIOP data, the optical depth profiles, their quality flags and the number of cloud layers in a column were extracted and used in this analysis. Also included in the files were the C6 MYD06 optical depths and particle sizes.
Included in the PEATE CALIPSO/MODIS matchup files are the optical depth and effective radius from the CALIPSO/IIR official product described by Garnier et al. [4,5], which have been validated by field campaigns [10]. For comparison, the IR optical depths were converted to visible optical depths using the suggested factor of 2.25 and the IIR-based diameters were converted to radii.
Another key component of this study is the ice crystal scattering library [11]. This library provides the single scattering properties for 445 wavelengths between 0.2 and 100 µm and 189 particle sizes ranging from 2 to 10,000 µm. The data are provided for three levels of particle roughening (smooth, moderate and severe). The data are given for 11 habits including three quasi-spherical ice crystal habits (droxtals, prolate and oblate spheroids), solid bullet rosettes, hollow bullet rosettes, solid columns, hollow columns, plates, aggregate of solid columns (henceforth the aggregate column), small aggregate plates and large aggregate plates. For the quasi-spherical habits, we only consider the droxtal in this study, that is, we only consider nine habits. For inclusion in the retrieval scheme, the single scattering albedo (ωo), the asymmetry factor (g) of the phase function, and the extinction efficiency (Qe) are integrated for each habit over the MODIS/AQUA spectral response functions. The values are fit as functions of effective radius using fourth order polynomials. The cirrus particle size distribution was assumed to be a modified Γ distribution [12] with an effective variance of 0.1. This size distribution is the same one used for the daytime MODIS cloud products. For the remainder of this paper, particle size will refer to the effective particle size based on this assumed distribution.

Retrieval Method
This section provides the basis of the retrieval, first by explaining the general principles and then the two variants of the technique used later in the study.

General Method for Effective Radius Retrieval
As stated above, C6 MYD06 contains cloud emissivity (ec) for the 8.5, 11 and 12 μm MODIS channels (bands 29, 31 and 32). ec is computed in the MYD06 cloud height algorithm using the following equation (1) where I is the observed radiance, Iclr is the computed clear-sky radiance, Iac is the radiance component coming from the atmosphere above the cloud, tac is the transmission along the path above the cloud and B(Tc) is the Planck emission computed at the cloud temperature (Tc). The clear-sky radiative transfer model used in the MYD06 processing is the PFAAST model [13]. The MYD06 CO2 slicing cloud-top pressure algorithm assumes the clouds to be isothermal, non-scattering layers and does not account for multiple cloud layers [2]. Ratios between different pairs of CO2 absorption channels are first used to retrieve cloud-top pressure; cloud temperature is derived subsequently by comparing cloud pressure with the National Center for Environmental Prediction Global Data Assimilation System [14] model profile for upper level clouds.
As described by Parol et al. [15], cloud emissivity at two different spectral bands (x, y) can be used to compute ratios of absorption optical depths (β) as follows: As further described in Parol et al. [15], β values can be approximated very accurately based solely on the single scattering properties of ωo, g, and Qe using the following relationship: The spectral variation of the single scattering properties for solid bullet rosettes from 8 to 13 μm is shown in Figure 2. The β values shown in Figure 2 are referenced to 11 μm. Note the significant spectral variation of the properties through the 8-12 μm IR window region. For the ice cloud bulk scattering properties, the values of the single scattering properties in the ice particle database are integrated over C6 MYD06 size distribution and over the spectral response functions shown in Figure 2. Though roughening has little impact at IR wavelengths, only the severely roughened crystals were used to be consistent with the MODIS C6 optical property retrievals. Figure 3 shows the variation in the 11-12 μm β value as a function of effective radius for each of the nine habits. Figure 4 shows the same plot for the 11-8.5 μm β values. Based on the β values, the retrieval of particle size is straightforward. The three values of ec generate two values of β. As Figures 3 and 4 show, the relationship between each β value and effective radius (re) is monotonic and two values of β generate two values of re. In this technique, the final effective radius value is the average of the individual effective radii deduced from the 11-12 μm β and the 11-8.5 μm β.

General Method for Optical Depth Retrieval
The next step in this analysis is the estimation of the cirrus cloud optical depth (τvis), which can be derived after effective particle size is retrieved. Based solely on the cloud emissivity (ec), one can estimate the absorption optical depth (τabs) using the following relationship: where μ is the cosine of the viewing zenith angle. Evoking the scaling approximations of van de Hulst [16] allows us to estimate the full optical depth at the IR wavelength by the following relationship: 1 If we assume the extinction efficiency at visible wavelengths is approximately 2, we can derive the full visible wavelength using this relationship 2 At this point it is worth exploring the veracity of the scattering approximations using explicit radiative transfer calculations. The radiative transfer calculations were performed using an adding-doubling model to simulate a single layer isothermal cloud above a black-body surface. All atmospheric effects were ignored. The input to the model was the optical depth and effective particle radius. The model computed top-of-atmosphere radiances at 8.5, 11 and 12 μm. Cloud emissivities were then computed using Equation (1) and the values of β were computed using Equation (2). Figure 5a shows the relative errors in optical depth caused by using Equation (3) as a function of the 11 μm ec and the 11-12 μm β value. This figure shows the approximation for computing the full optical depth from the absorption depth and the single scattering properties is accurate to within 1% for most of the β and ec values. Only for small values of ec and large values of β do the errors exceed 3%. Higher 11-12 μm β values correspond to smaller particle sizes where scattering plays a larger role in the total extinction process. In Figure 5b a density plot using MODIS data is plotted, which demonstrates that data mostly concentrate when emissivity is large and 11-12 μm β is small. Therefore, retrieval errors are expected to be small using the approximations in Equation (3).

C6 Habit
As stated earlier, we employed two variants of this technique. The first was to use a microphysical model that is consistent with that used in C6 MYD06 solar channel retrievals. The ice cloud optical property model employed for MODIS C6 is different from that assumed for earlier C5 products. The model for C6 assumes a distribution of severely roughened aggregate columns. The model for C5 assumed a habit mixture, but the crystals had smooth surfaces, with the exception of the aggregate columns. Several investigations raised issues with the C5 models that led to this change for C6. The C5 MYD06 optical depth values for cirrus were significantly higher than those reported by CALIPSO/CALIOP (e.g., [17]). While both MODIS and CALIPSO algorithms are in active changes and new versions are being released, the latter one retrieves cloud properties from lidar signal returns and are generally considered more accurate than passive satellite retrievals. The difference between the CALIPSO and MODIS C5 optical depth values was attributed to the C5 microphysical values having a large asymmetry parameter (g) in the solar channels that resulted in an overestimation of the optical depth (τ). The same issue was found in a comparison of MODIS C5 cloud properties to those from POLDER (POLarization and Directionality of the Earth's Reflectances) [18], which assumed the use of a roughened ice particle. The adoption of severe particle roughening resulted in lower values of g, which in turn brought the optical depth/particle size values more into line with those from CALIOP [19]. In defining the models of C6, consistency with CALIPSO/CALIOP was paramount. The nine habits of Yang's database were analyzed and the severely roughened aggregate columns were chosen due to their agreement with CALIPSO/CALIOP. Therefore the first variant of the IR technique is to use the β curves in Figures 3 and 4 corresponding to the aggregate columns (solid grey line).

Empirical Method
As described above, the final re value is calculated as the mean of the two spectral values (11-12 μm and 11-8.5 μm). It is planned to run this method on the MODIS C6 data using the same habit as that used in the C6 solar reflectance retrievals. However, it was noted that the aggregate column habit did not often produce spectral re values that agreed with each other. A choice of other habits often reduced this spectral inconsistency. Therefore, we have implemented a second variant of the technique where the results are based on the re value from the habit that produced the highest level of IR spectral consistency. Spectral consistency is defined as the spectral re value whose difference was either less than 20% of the mean re value or 1 μm. For the data analyzed here, roughly 70% of the retrievals met the criteria for spectral consistency. For simplicity, all of the results that achieved spectral consistency were used to construct empirical "radiometric" models. Figures 6 and 7 show the resulting empirical models for the two β curves. Similarly, the resulting models for the other single scattering parameters (Qe, g, ωo) are constructed by selecting data from the habit when the spectral consistency is achieved. Figures 6 and 7 show the density of the retrievals with each β and re bin. The β bin size is 0.01 and the re bin size is 2 μm. The colors show the relative density of the position of spectral consistent solutions in the β and re space. The white solid curves plot the smoothed mean β values within each re bin. This smoothing has been done many times for the purpose of a smoother transition between nearby re bins. The values shown by white solid lines therefore constitute the radiometric empirical model. In the radiometric empirical model, no additional habit assumption is made. Because the empirical model is based on the spectrally consistent retrievals, the spectral consistency of its re retrievals is much superior to that with the aggregate column habit results. This method does not attempt to infer meaning in the choice of habits that result in spectral consistency, nor do we imply that the empirical model is better than the aggregate column model, but one could view this approach as suggesting that the variation in the scattering properties in Yang's database supplies the needed range and the empirically derived curves in Figures 6 and 7 represent the optimal radiometric model. The inference of the habit mixtures and size distributions that could produce this model is left to future studies.

Results
To demonstrate this retrieval, the two variants (C6 and Empirical) were applied to all of the CALIOP/MODIS matchups from January 2010. The data were filtered as follows: in all results showing CALIOP optical depth comparisons, only single layer constrained CALIOP retrievals between 60°N and 60°S were shown (CALIOP QC flag is equal to 1), since they are generally considered having the highest quality. Here, constrained retrievals are retrieved from a measured two-way transmittance, whereas unconstrained retrievals use an a priori value of the lidar ratio [20]. Lidar ratio uncertainties are one of the dominant factors for optical depth estimation. Constrained CALIOP results are rare during the day and therefore most of the CALIOP optical comparisons occur at night. The CALIPSO/CALIOP QC filtering was not applied in the comparisons to C6 MYD06 or CALIPSO/IIR. In all comparisons, all results with optical depths greater than 2 were excluded. For the effective radii comparisons, IR results with values greater than 50 μm were excluded. The following results show comparisons of the IR retrievals against CALIPSO/CALIOP, MYD06 and CALIPSO/IIR over both land and ocean using both the C6 aggregate column habit and the Radiometric Empirical Model. Figure 8 shows the comparisons of visible optical depths based on IR retrieval to those from CALIPSO/CALIOP, C6 MYD06 and CALIPSO/IIR, which all fall on CALIPSO track. Other than C6 MYD06, the other two products are available during both day and night time. The results are shown as density scatterplots with the IR-derived results on the vertical axis. The left column show the IR results using the aggregate column habit and the right column shows the results using the Radiometric Empirical Model. The top row shows the comparisons to the CALIPSO/CALIOP product. The middle rows show the comparisons to the C6 MYD06 product. The bottom row shows the comparison to the CALIPSO/IIR product. The magnitudes of the biases are lower and the correlations are higher for the Empirical Model results than Aggregate Column Habit results for each of the three comparisons. The worst results in terms of bias and correlation are for the aggregate column comparisons to the C6 MYD06 results, and the best results are for the Empirical Model compared to the CALIPSO/IIR. The range of variability in the biases and correlations is small. Therefore, if optical depth is the desired parameter, there is little justification for using the Empirical Model over the aggregate columns. The agreement of the IR-derived optical depth with these three sources is the major accomplishment of this technique. This empirical technique, which requires no more additional information than that contained within C6 MYD06 and computationally efficient radiative transfer approximations, yields accurate optical depth retrievals over the range of 0.1 to 2. The consistency of the solar reflectance and IR-derived optical depth values is a major improvement over C5 MYD06. This consistency will enable progress to be made on inclusion of more spectral information in future version of MODIS retrievals. This consistency also adds confidence to use of the IR-derived cloud properties for diurnal studies.  Figure 9 shows the comparison of the effective radii retrievals. CALIPSO/CALIOP makes no effective radius product. The top row shows the comparison to the C6 MYD06 results using the 2.1 μm channel. The bottom row shows the comparison to the CALIPSO/IIR results. An obvious feature found from comparison of Figures 8 and 9 is the decrease in the level of agreement in the effective radii results compared to the optical depth comparisons. With use of C6 MYD06 as a reference, the Empirical Model results give the smallest bias magnitude and highest correlation. The opposite is true when using the CALIOP/IIR results as the reference. As described above, the Radiometric Empirical Model is constructed using results that maximize spectral consistency across the 8-12 μm region. The conclusions from the top row in Figure 9 are encouraging, since they imply that the radiometric empirical model may produce results that are spectrally consistent across the 8-12 μm region, and also in the VIS/SWIR and IR regions. However, these comparisons among different satellite-based retrievals are not conclusive and a more objective validation is needed. One omission from this study is the impact of the vertical variation in particle size within the cloud and its influence on the retrievals. While work is ongoing in this area [21], a method to include these effects in the technique described here is unknown.

Conclusions
This paper presents a method to derive the optical depth and particle size for cirrus clouds from the 8.5, 11 and 12 μm cloud emissivities included within the C6 MYD06 series of MODIS cloud products. The method employs scattering approximations and a state-of-the-art scattering library that provides single-scattering properties for a set of nine different ice particle habits. Two models are employed here. One model is based on an aggregate of solid column particles and is the same model used in C6 MYD06 products. The other model is a radiometric empirical model that emphasizes spectral consistency in re retrievals. Both models produce optical depth values that agree well with the constrained CALIPSO/CALIOP, C6 MYD06 and CALIPSO/IIR products. The effective radii results are less consistent. While the radiometric empirical model generated effective radii results that are in better agreement with C6 MYD06, the use of the aggregate of solid column model results in better agreement with CALIOP/IIR. Future work includes applying this technique to C6 MODIS data to generate IR-based global maps of these cirrus properties and their diurnal variation and more research into extracting information on crystal habit and size distribution. of many figures. Bryan Baum provided his expertise in the scattering models and expected behavior of non-spherical particles. He also edited the manuscript extensively. Robert Holz provided the CALIPSO matchup files and gave guidance on the use and interpretation of the LIDAR results. Steven Platnick gave guidance as co-investigator on all aspects and specifically guided the analysis of the MODIS solar-reflectance results. Ping Yang provided the scattering property library and provided guidance on its use.

Conflicts of Interest
The views, opinions, and findings contained in this report are those of the author(s) and should not be construed as an official National Oceanic and Atmospheric Administration or U.S. Government position, policy, or decision.