CERES MODIS Cloud Product Retrievals for Edition 4—Part I: Algorithm Changes

The Edition 2 (Ed2) cloud property retrieval algorithm system was upgraded and applied to the MODerate-resolution Imaging Spectroradiometer (MODIS) data for the Clouds and the Earth’s Radiant Energy System (CERES) Edition 4 (Ed4) products. New calibrations for solar channels and the use of the 1.24- $\mu \text{m}$ channel for cloud optical depth (COD) over snow improve the daytime consistency between Terra and Aqua MODIS retrievals. Use of additional spectral channels and revised logic enhanced the cloud-top phase retrieval accuracy. A new ice crystal reflectance model and a CO2-channel algorithm retrieved higher ice clouds, while a new regional lapse rate technique produced more accurate water cloud heights than in Ed2. Ice cloud base heights are more accurate due to a new cloud thickness parameterization. Overall, CODs increased, especially over the polar (PO) regions. The mean particle sizes increased slightly for water clouds, but more so for ice clouds in the PO areas. New experimental parameters introduced in Ed4 are limited in utility, but will be revised for the next CERES edition. As part of the Ed4 retrieval evaluation, the average properties are compared with those from other algorithms and the differences between individual reference data and matched Ed4 retrievals are explored. Part II of this article provides a comprehensive, objective evaluation of selected parameters. More accurate interpretation of the CERES radiation measurements has resulted from the use of the Ed4 cloud properties.

MCAT-ML radiance and first guess radiance in channel k. r U , r L MCAT-ML upper and lower layer cloud particle effective radius. r e , D e Cloud particle effective radius and diameter. r ea Cloud particle effective radius from in-situ data. T Temperature. Cloud optical depth at Z maxRH . τ U , τ L MCAT-ML upper and lower layer cloud optical depth. τ 4 , τ 12 Cloud optical depth for channels 4 and 12.

I. INTRODUCTION
T HE CERES Project [1] has been making global measurements of the earth's radiative energy budget and cloud properties from Terra (1030 LT equatorial crossing time) and Aqua (1330 LT equatorial crossing time) since March 2000 and July 2002, respectively. The primary measurements used by CERES are taken by two broadband scanners [2] and a MODIS (see [3]) on each satellite. These measurements and their conversion to physical parameters are designed to create a long-term climate data set for monitoring the radiation budget and for exploring the relationships among clouds, aerosols, radiation, and various surface and atmospheric variables (see [4]- [9]). As they are valuable for tracking climate parameters (see [10]- [12]), for improving our understanding of climate (see [13]- [19]), assessing climate model output (see [20]- [23]), and other applications (see [24]), the CERES data are occasionally reanalyzed with the latest calibrations and methods to increase their accuracy and consistency.
The CERES data processing system first calibrates and geolocates the scanner data. To compute fluxes, it then branches into two sequences of subsystems: 1) a mimic of an older simpler analysis system relying only on the scanner LW and SW radiance observations and 2) the CERES standard approach, a more complex methodology that incorporates radiances and cloud and aerosol properties from a collocated conterminous satellite imager (MODIS for this article) and geostationary satellite (GEOsat) data that cover the entire diurnal cycle [25], [26]. The former methodology determines only the TOA fluxes using older ADMs to convert radiances to fluxes (see [27]) and to assist interpolation of the fluxes over the diurnal cycle to estimate the entire radiation budget [28]. After classifying each imager pixel as cloudy or clear [29], the CERES standard approach follows a sequence of steps that are contingent on each previous step. The first of these is the retrieval of the cloud properties from MODIS that are merged with the aerosol properties from MODIS [30], [31] and the CERES SW and LW radiances to compute the TOA and surface fluxes. This is followed by gridding and averaging of the instantaneous fluxes and cloud properties; estimation of the hourly fluxes from the instantaneous values with and without inclusion of GEOsat data; merging of the Terra and Aqua data; and, finally, computation of the parameter averages at several space and time scales. Each set of calibrations, methods, algorithms, and auxiliary data used in various subsystems comprises a CERES Edition. The Edition number is incremented if a major change occurs in the calibrations, processing method, or auxiliary data, in more than one subsystem. Notable, but less comprehensive changes having minor effects on particular subsystems are identified by appending a letter to the Edition number.
In the first step of the CERES standard methodology, the cloud property retrievals are merged with the aerosol and scanner flux data to create the SSF (see [32]) product for each satellite. The SSF parameters form the basis for the products computed in the subsequent steps of the processing system. In generating the SSF, the cloud properties retrieved from the sampled 1-km MODIS data are convolved into the nominal 20-km scanner footprint. Their averages and the accompanying aerosol information define the scene identification, which determines the ADM for converting the footprint's broadband radiances to the TOA fluxes. The ADMs were developed using those same parameters merged with scanner footprints from the CERES rotating the azimuth plane scanner [33]. The cloud properties are also used in the parameterized methods for estimating the surface fluxes that are included in the SSF [34]. Downstream in the CERES processing, the SSF cloud and aerosol properties are crucial for calculating the CERES Cloud and Radiation Swath product and CERES System Synoptic product, which together provide instantaneous footprint and global 3-or 1-hourly 1 • -gridded TOA, surface, and inatmosphere fluxes [35]- [38]. These multiple uses make the accuracy and consistency of the cloud properties derived from MODIS a critical component of the CERES system.
The Edition 2 (Ed2) cloud detection and retrieval algorithms described in [39] and [40] were used to provide the MODIS-based cloud properties for the CERES Ed2 and Ed3 SSF products. Through comparisons with other observations and retrievals, a variety of issues were identified that if resolved could improve the accuracy and consistency of the CERES-MODIS (CM) cloud properties. Among others, these issues include calibration, coding errors, model and algorithm shortcomings, and biases in certain parameters. Technological and theoretical advances since 2002, when Ed2 processing began, have the potential for improving the characterization of clouds from passive imager data. To address the Ed2 issues and take advantage of progress during the following decade, the CERES Ed2 cloud property retrieval system has been updated to provide both revised and new cloud properties that are included in the CERES Edition 4 (Ed4) products. The CERES Ed4 mask [29] and retrieval algorithms also comprise part of the Satellite ClOud and Radiative Property Retrieval System (SatCORPS). SatCORPS processes the global constellation of GEOsats [41] along with the regional AVHRR and MODIS data sets to provide near-real-time historical cloud properties for weather (see [42]- [44]) and climate (see [45], [46]) applications. The SatCORPS also supports CERES, producing the hourly GEOsat cloud property information used in the greater CERES Edition 4 processing system (SatCORPS-Ed4-GEO). A later paper will describe the SatCORPS-Ed4-GEO in detail. The Ed4 retrieval algorithms were implemented in 2012 and, at the time of this writing, are still being used to analyze the MODIS data as they become available.
This article provides an overview of alterations to the retrieval process and their impact. It is the first of two parts. Significant changes include, among others, the use of a new ice crystal model for retrieving the ice cloud properties, a different spectral channel for retrievals over snow, a revised cloud phase determination scheme, a new multilayer detection and retrieval method, and new approaches to determining CTH and thickness. These changes have resulted overall in more accurate cloud properties compared with Ed2. A brief review of the data used for Ed4 and the Ed2 methodologies is provided in Section II. Section III summarizes the major changes in the analysis, followed by Section IV, which presents some results. After Section V discusses the Ed4 results, the concluding remarks are given in Section VI. An extensive analysis comparing the CERES Aqua cloud products with those from the active sensor data is presented in Part II [47].

II. DATA
The input data consist of imager radiances and several ancillary data sets and models used to estimate the expected cloud-free spectral radiances and to simulate the cloudy sky radiances for different heights, optical depths, and particle sizes for both ice and liquid water clouds. Other data are used for evaluating the results.

A. MODIS Radiances
For Ed4, the CERES ingests a 19-channel subset of the 36-channel MODIS Collection 5 Level 1B geo-located and calibrated radiance data. The 1-km MODIS data are sampled every fourth pixel and every other scan line to yield an effective nominal resolution of 4 km × 2 km or 8 km 2 . The CERES Ed4 cloud retrievals use eight MODIS channels having central wavelengths at 0. 65, 1.24, 2.13, 3.75, 8.55, 11.0, 12.0, and 13.3 μm. These channels along with four others are used in the Ed4 cloud mask. In referring to each channel in algorithm descriptions, the CERES uses a different channel numbering system than the original MODIS numbers because similar versions of the algorithms are applied to data from other satellites having different numbers for comparable channels. For the solar and thermal channels, the radiance parameters are given as reflectance ρ k and brightness temperature T k , where the subscript k denotes the CERES channel number. The CERES numbering system for these retrieval channels is given in Table I along with the acronyms used as the reference names and the radiance parameter variable names.
The nominal MODIS Collection-5 (C5) calibrations, based on [48] and references therein, produced significant differences between Terra and Aqua for some channels. These biases caused some discrepancies in the retrieved parameters [49], [50]. To account for calibration differences between Terra and Aqua MODIS, the Ed4 Terra reflectances were normalized to their Aqua counterparts for several channels using the calibration adjustments reported by Sun-Mack et al. [51].
The impacts of the calibration normalization process on the Terra Ed4 CODs, effective particle size, and NP cloud phase selection are also reported in [51]. In addition to the Terra calibration adjustments, the nocturnal 3.8-μm channel data from both the satellites were processed through a destriping algorithm because the response characteristics were not uniform for all the MODIS detectors [29].
The MODIS C5 generation ended in February 2016. Thus, continued Ed4 processing from March 2016 onward required the use of MODIS Collection 6.1 radiances, which, in some cases, have different calibrations relative to the C5 data [51], [52]. To minimize calibration discontinuities between the Ed4 C5 and C6.1 cloud retrievals, all C6.1 radiances were renormalized to their respective C5 radiances using matched MODIS granules for particular days. The data used in the normalization process and the linear correlation coefficients are given for Terra and Aqua at https://satcorps.larc.nasa.gov/ SCATTER-TERRA5 and https://satcorps.larc.nasa.gov/ SCATTER-AQUA5, respectively. The official CERES Ed4 record uses C5 radiances through February 2016. Thereafter, it uses the C5-normalized C6.1 data sets. The results presented here reflect that record.

B. Ancillary Input
As in Ed2, the Ed4 MODIS data are processed on a tile basis to minimize the number of computations for the cloud mask and retrieval algorithms. Each tile nominally represents a 32 km × 32 km area consisting of 8 pixels across the track and 16 along the track. The ancillary data used in the retrievals are kept in their native resolution and the data from the grid box center closest to the center pixel of the tile are used in the processing. For time-dependent variables, the data are interpolated to the time of the observation. Thus, only one set of viewing and illumination angles, soundings, and other data listed below are used for analyzing all the pixels in a given tile. As in the Ed2 processing, the global surface skin temperature; the surface wind speed; and the atmospheric temperature, ozone, and humidity profiles; as well as the total PW vapor are taken from the CERES MOA data set for Ed4. Reanalyses from version 5.4 of the Global Modeling Assimilation Office Global Earth Observing System GEOS-5, an update of the versions described in [45], provide the MOA with the algorithm-consistent surface skin temperature and vertical profiles of temperature, humidity, and ozone throughout the Ed4 record. The ancillary and clear-sky radiance data used for Ed4 are described in [29]. The ozone profiles from the GEOS-5.4 products replace those used in Ed2. GEOS-4 analyses [54] used through December 2007 in the Ed2 MOA data set were replaced with the GEOS-5.4 reanalyses, so that, unlike Ed2, Ed4 processing uses the same reanalysis product throughout its record. The GEOS-5.4 vertical profiles are available at a nominal horizontal resolution of 0.5 • × 0.625 • every 6 h, while the surface skin temperature T s is provided every 3 h. The total column water vapor values are taken over ocean from the Special Sensor Microwave Imager product at a 25-km resolution [55]. The Ed4 MOA interpolates all the data to time and space resolutions of 3 h and 0.5 • × 0.5 • , respectively, while Ed2 used a 1 • latitude-longitude grid.

C. Other Data
In 2006, CloudSat [56] and the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO, [57]) were launched into orbits nearly identical to that of Aqua as part of the Afternoon Constellation, known as the A-Train. The CALIPSO, CloudSat, CERES, and MODIS (C3M) merged product [58] was used to develop a new ice cloud thickness parameterization (Section III-B.2).

III. ALGORITHM CHANGES FOR Ed4
Only four cloud properties are retrieved directly for each pixel from the MODIS radiances. The direct retrievals are cloud emissivity ε c , CET T c , COD τ , and CER r e . All the remaining cloud properties are estimated based on those retrievals, radiances, and auxiliary data using various parameterizations and tests. In Ed4, the surface skin temperature, T s , was added to the products retrieved by the cloud algorithms. It is retrieved for each tile whenever the clear fraction exceeds 15% or T s is greater than the MOA skin temperature T s (MOA) and the clear fraction exceeds 0.0. In addition, whenever the retrieved value is more than 10 K greater than T s (MOA), the retrieval was discarded. That coding error mainly affected high temperatures over arid regions. Selected MODIS radiances and cloud properties are matched with the corresponding CERES scanner footprint and classified as belonging to one of the two cloud layers, as multilayered, or clear (see [40]). The properties of the pixels are convolved with the point spread function of the CERES scanner footprint to compute the averages and standard deviations for each of the four categories. The results are then included with various CERES radiation parameters, navigation data, viewing and illumination angles, auxiliary data, aerosol data from external sources, and various flags to comprise an SSF. A full listing of the SSF parameters is provided in [32]. This article deals only with the pixel-level properties before they are convolved into the SSF product.

A. Primary Cloud Property Retrievals
For a given pixel, the parameters, CET, COD, and CER, are retrieved using one of the three methods that were also used for Ed2 [40], albeit with some changes. The VISST, applied during the daytime (SZA < 82 • ) over snow-free surfaces, uses the VIS (0.65 μm) channel to estimate τ , the IRW (11.0 μm) channel for T c , and the SIR (3.75 μm) channel for the primary particle size. The SPW (12.0 μm) channel is used to aid the phase selection. If the underlying surface is determined to be snow-or ice-covered either from the snow-ice maps or from identification of the nearby pixels as clear snow, then the SINT is applied during daytime [40]. The SINT is the same as the VISST, except that the SNI (1.24 μm) channel replaces the VIS channel for optical depth retrievals. This differs from the Ed2 SINT, which used the 1.6-μm channel on Terra and the 2.1-μm channel on Aqua to replace the VIS channel. The SIST, used for nighttime retrievals over all surfaces, is the same as that in Ed2. The maximum COD allowed in the Ed4 retrievals is 150 compared with 128 for Ed2. While larger values of COD occur in reality, the limit is imposed to minimize severe overestimates at higher SZAs due to 3-D cloud effects.
All three of the techniques seek to minimize the differences in the observed radiances and computed TOA radiances. The latter are calculated using liquid water and ice reflectance and emission models as functions of CET, COD, and CER (or T c , τ , and r e ), as described in [40]. There are four possible outcomes of this process: 1) a solution is found for a water cloud model (i.e., CETW, CODW, and CERW); 2) a solution is found for an ice cloud model (CETI, CODI, and CERI); 3) a solution is found for both a water and an ice cloud model; or 4) no solution is found using any of the models. The phase for the variable symbol, T c , τ , and r e , is denoted by adding "W" or "I" to the subscript. A decision to select one of the results or to fill with default values is made in the phase selection process described in Section III-A.5.
The water droplet reflectance LUTs used in Ed4 are the same as those used for Ed2, except that new LUTs were computed for the 1.24-and 2.13-μm channels using the same approach of [59] and the indices of refraction were taken from [60]. Corrections for atmospheric absorption use the radiative transfer calculations described in [40], using the correlated K -distribution method with the absorption coefficients computed for the spectral response functions of various channels.
A number of coding errors and other programming bugs discovered after implementing Ed2 processing were corrected for Ed4. While too numerous to cover here entirely, the two major ones dealt with the atmospheric absorption calculations. Other small constraint changes were also implemented. For example, the valid maximum brightness temperatures were assumed to be 330 K. Any temperatures exceeding that value are reset to 330 K. Ed2 had no such restriction. Overall, these changes and corrections tend to decrease the number of no-retrieval pixels, reduce the blockiness in the output, and produce fewer CER extremes.
1) Ice Crystal Model: A significant change for all the standard methods is the use of a new ice crystal reflectance model. To represent the composition of ice clouds, the Ed2 used smooth hexagonal ice columns that have asymmetry factors g between 0.77 and 0.85 at 0.65 μm depending on the crystal size [59]. These relatively high values of g, on average, lead to larger values of τ than observed with other techniques (see [61], [62]), as the number of samples taken at scattering angles greater than 90 • far exceeds those taken at smaller angles [63]. For the VISST, an overestimate of τ for thin clouds will cause an overestimate of T c and, subsequently, an underestimate of the CEH and sometimes an incorrect phase selection. The ice crystal reflectances used in Ed4 are based on radiative transfer computations using distributions of hexagonal ice columns with roughened surfaces having the normalized roughness parameter set equal to 1.0 [63]. These models, with g ranging from 0.77 to 0.81, should reduce the COD and increase the CER for a given ice cloud retrieval [64]. The roughened crystal model also produces, relative to obser-vations, better angular consistency in reflectance than the smooth crystal model and five other crystal formulations [65]. Thus, the roughened crystal model was adopted for Ed4 and used to compute reflectance and albedo LUTs for the 0.65-, 1.24-, 1.64-, 2.13-, and 3.75-μm channels.
Instead of using the ice crystal effective diameter D e definition used in Ed2, which is based on the formula of [66], Ed4 uses the ice crystal effective radius, CERI or r eI . CERI is defined in the same manner as the water droplet effective radius, CERW or r eW , which is the ratio of the third moment to the second moment of the particle size distribution. For CERI, the ice crystal dimensions used in a given distribution are defined as the radii of equivalent spheres. The relationship between the Ed4 and Ed2 dimensions can be approximated as The Ed4 ice crystal size distributions are the same as those used for Ed2, but the crystal facets are roughened. Because it is consistent with other retrievals of ice crystal effective radius, this revision in the ice crystal size definition will facilitate its comparison with those from other retrievals.
2) Clear-Sky Values: Clear-sky reflectances ρ csk and brightness temperatures T csk used for each channel k in the VISST, SINT, and SIST are described, for the most part, in [29] and references therein. Additional information about the clear albedos and reflectances for the 1.24-and 2.13-μm channels are provided in [67] and [68] and [69], respectively. A change to the procedures for specifying the clear-sky reflectance in the retrieval was made for Ed4. In Ed2, the clear-sky reflectances were changed from the predicted values to the average observed value in clear portions of the scene whenever 10% or more of the tile was classified as clear, but in Ed4 this practice was discontinued. Ed4 always uses the predicted clear-sky reflectances for all solar channels, and the clear-sky temperatures are revised in the same manner used in Ed2.
It should be noted that the use of the 1.24-μm channel instead of the 2.13-μm reflectances over snow surfaces comes with increased uncertainty because of the greater variability of the snow surface albedo and bidirectional reflectance at 1.24 μm compared with those of the more absorptive 2.13-μm channel. In addition, that uncertainty is exacerbated by the 1.24-μm snow albedos (see [70]), which are much greater than their 2.13-μm counterparts that have values typically between 0.02 and 0.06 [29]. The greater 1.24-μm surface albedo reduces the contrast between the cloud and surface reflectances necessitating more accurate characterization of the clear-sky radiance field. Unfortunately, the approach taken in Ed4 to estimate clear-sky 1.24-μm reflectance, on average, overestimated the true clear-sky value. Following the approach of [29] for comparing the observed and estimated clear-sky reflectances, it was found that the mean difference between the predicted and observed reflectances is 0.030 ± 0.045 over permanent snow-ice areas. For other surface types, this difference ranges from 0.021 ± 0.083 over bog/tundra to 0.030 ± 0.066 over pastureland. Because of these overestimates, the predicted value often exceeded the observed values for both clear and cloudy pixels. Furthermore, at some viewing-illumination angle combinations, the reflectance of the clear areas exceeds that of the cloudy pixels. To obtain a solution in those instances, the 1.24-μm clear-sky reflectance was estimated as the lowest observed reflectance minus 0.01 of the observed cloudy 1.24-μm reflectances in the tile. The impact of this re-estimate procedure is likely to cause an increase in the retrieved COD for those cases having optically thin clouds.
3) Multispectral Retrievals: Additions to the cloud property complement for Ed4 include retrievals of CER over nonsnow surfaces using the 1.24-and 2.13-μm channels, CER 7 and CER 2 , respectively. The standard retrieval at 3.75 μm is denoted simply as CER, since it is the CERES default value. Fig. 1 plots the diffuse albedos for the water droplet distributions and smooth and roughened ice crystal models at 1.24, 1.62, and 2.13 μm. The approximate limits of the retrievable optical depth for each wavelength are indicated with the green arrows. The green arrows roughly correspond to the lowest optical depth associated with the maximum albedo for the middle particle size. The actual reflectance limit depends on the particular particle size and the viewing and illumination conditions. The 1.24-μm albedos vary over a greater COD range for both liquid [ Fig. 1 Aqua, but the 1.62-μm smooth crystal models were used in the Terra Ed2 retrievals. It is noteworthy that for a given CODI, the roughened ice crystal model albedos exceed their smooth crystal counterparts. The greater COD limits (60-80 for liquid and 30-50 for ice) for 1.24 μm makes it more valuable for COD retrievals compared with 1.62 μm (limits: 20-30 for liquid and 8-16 for ice) and 2.13 μm (limits: 8-25 for liquid and 4-15 for ice), but the greater sensitivity of the longer wavelengths to changes in particle size renders them more valuable for particle size retrieval. Yet, the greater τ -range for 1.24 μm, accompanied by some sensitivity to particle size, indicates it could possibly be valuable for providing a value of CER that may be more representative of the entire cloud, particularly for thicker clouds. To retrieve these new experimental parameters, the VISST-retrieved values of T c and phase are used to force a solution for CER 7 and CER 2 given the observed and predicted clear-sky reflectances for their respective channels. Unfortunately, an undetected error in the 1.24-μm ice cloud LUTs affected the retrievals of CERI 2 and CERI 7 . Therefore, those values are unreliable and should not be used. This error did not affect the retrievals of CODW over snow in Ed4.

4) Modified CO 2 -Absorption Technique:
To provide an alternative effective CTH estimate and a seed for an ML cloud detection algorithm, the two-channel MCAT of [71] and [72] is used to detect upper level clouds and retrieve the CO 2 cloud emissivity ε M , cloud top temperature T t M , height Z t M , and pressure P t M . The MCAT is a modified version of the original CAT. It has been modified to improve the calculation of the below-cloud upwelling radiance, which relies on the accuracy of the surface skin temperature and temperature and relative humidity profiles from the ancillary data. The MCAT processing first uses the original CAT to derive an effective cirrus cloud emissivity based on the surface skin temperature. The effective cirrus cloud emissivity is then used to inversely retrieve a background effective temperature through comparisons between the model-computed and observed upwelling IRW radiances. The background effective temperature can be the surface or a low-cloud temperature. It is derived with the presence of the cirrus cloud when the computed upwelling IRW radiance matches the observed upwelling IRW radiance. Then, the MCAT uses the effective background temperature to retrieve the cirrus cloud top properties. The MCAT reduces detection of false thin cirrus caused by biased surface skin temperatures or, to some extent, the presence of a lower cloud. When the MCAT and CAT retrieve different CTH results, a weighted average based on the relative humidity values corresponding to the two CTHs is calculated to serve as the MCAT CTH.
The MCAT also does not rely on the conventional assumption that cloud emissivities are equal at channels 4 and 12, ε 4 = ε 12 ; rather, it uses a parameterized scheme that accounts for the cloud emissivity changes through and where σ 4 e and σ 12 e denote the corresponding spectral extinction coefficients, and μ denotes the cosine of the satellite VZA [71]. The MCAT further uses the corresponding spectral extinction coefficients at channels 1 and 4 to derive the associated VIS COD τ M from the retrieved cloud emissivity ε M = ε 4 . If the MCAT finds a solution for P t M < 600 hPa, then the flag, MCAT_flg, is set equal to 1 and the results are used in determining the final retrievals of T c , r e , and τ for ice clouds; otherwise, the results are not used and MCAT_flg = 0.
There are several exceptions to using the VISST to estimate COD, CER, and CET during the day that are reserved for certain conditions. These include those clouds identified as thin cirrus for which a solution was not found with either VISST or MCAT. In these cases, it is assumed that the CET is the same as the temperature T maxRH at the atmospheric level Z maxRH having the highest relative humidity in the MOA vertical profile above the height corresponding to a temperature of 233 K. The cloud emissivity ε maxRH , computed using CET and the observed and clear-sky temperatures at 11 μm, is used to estimate CODI = τ maxRH using [73, eq. (28)] with their C20 ice crystal model and the assumption that CERI = 10 μm. All positive retrievals having τ I < 0.02 are reset, so that τ I = 0.02. The application of the above approach or the use of the MCAT result by itself depends on some very specific conditions that are tested during the phase selection process.

5) Cloud Thermodynamic Phase:
There is a dramatic difference in the Ed2 cloud thermodynamic phase partitioning between daytime and nighttime retrievals, primarily due to the reduced information content resulting from the lack of solar reflectance at night. The Ed2 nighttime algorithms used only the 3.8-, 10.8-, and 12.0-μm channels in the retrievals. For Ed4, the algorithms exploit information in other MODIS channels. As in Ed2, the Ed4 cloud phase is determined using a combination of direct model matching during the retrieval process for a majority of the pixels and a series of logical steps for the remaining pixels. The former process remains the same in Ed4. That is, each of the three retrieval methods attempts to match the observed radiances with the model calculations using both the ice crystal and water droplet models. Changes were made to the second part of the algorithm, which is used when the model matching does not produce a unique phase solution. Both the day and night Ed4 phase determinations use results of the BSM of [74]. It uses the channels 4 and 6 brightness temperatures and is applied to all the cloudy pixels outside of the main phase selection process. The BSM results serve as the input to the phase selection. a) Daytime phase selection: Fig. 2 shows an overview flowchart of the daytime phase selection. Detailed flowcharts of the entire phase selection process are provided at https://satcorps.larc.nasa.gov/CERES_algorithms. Much of the diagram in Fig. 2 is very similar to the Ed2 version [40] with the following exceptions.

1) Thin Cirrus Test:
This routine, represented by the blue box on the left side of the diagram, is applied in every case after a solution is reached in the base algorithm.
If the original solution is ice or it was determined for a cloud over a snow/ice surface, it remains unchanged. Otherwise, a series of tests are applied to decide whether the pixel should remain liquid or whether a phase for a no-solution pixel should be assigned. These tests use reflectances from channels 1, 2, 7, and 9 and brightness temperatures from channels 4 and 5 as well as T t M from the MCAT retrieval. Further tests use the parameters associated with the maximum relative humidity level approach for thin cirrus clouds: T maxRH , ε maxRH , and τ maxRH . 2) Dual-Solution Override: This subalgorithm is represented by the blue box at the bottom of Fig. 2, where it follows the dual-solution phase discrimination process. This routine replaces the "1.6/0.6" ratio and "ice cloud likely" subalgorithms used in Ed2 (see [40]), but is located downstream of the "Check LBTM" subalgorithm instead of being an alternative to those two Ed2 algorithms. It applies tests that use CET, T t M , ρ 9 , and the BSM phase result, as well as particle size constraints to select a phase. The LBTM is the layer bispectral threshold method, an older algorithm that uses only the VIS and IRW channels [75].

b) Night phase selection:
The nighttime algorithm discriminates between thin and thick clouds because less information about optically thick clouds is available in the infrared radiances. The test for distinguishing thick from thin clouds is as follows: then the cloud is considered to be "thick." The main part of the nocturnal phase algorithm in Fig. 3 shows that when the cloud is thick, a few simple tests are applied to select a final phase and assign default values for τ and r e or to continue on to the algorithm used for all "thin" clouds. That latter process applies several tests based on which phase possibilities came from the SIST solutions. The tests are primarily temperature or BTD 34 comparisons. In one case, the SIST errors, err w and err i , for the water and ice solutions, respectively, are used to reach a decision. Ultimately, the process goes directly through to "Postretrieval" or indirectly through the "Check BSM" subalgorithm as indicated in most of the blue boxes. The "Postretrieval" procedure, indicated in the separate diagram at the bottom of Fig. 3, checks for bad channel 3 data and for clouds warmer than the surface to either select a final phase based on T 4 or keep the incoming phase as the final selection. In the "Check BSM" subalgorithm shown in Fig. 4, the phase estimated with the BSM is used to determine whether the incoming phase remains the same or is changed prior to the "Postretrieval" procedure. It uses results from the MCAT and the BSM results to arrive at a final selection. The first step determines whether the BSM was able to return a solution, as indicated by the flag, BSM_flg, being equal to 1. If not, or if the no-phase solution has been made, it goes directly to Postretrieval. The procedure will then change the phase if certain criteria are met. These criteria are compared with T cW , T cI , Z t M , BTD 64 , MCAT_flg, and CERW in the tests. The comparison using CER is shown in the lower right portion of Fig. 4. If CERW > R fit, the phase is changed to ice, based on the assumption that the CERW is likely to be smaller if the cloud is liquid. The formula for R fit is shown in the gray box in Fig. 4. c) Clear-sky restoration: For cloudy pixels meeting certain criteria during the daytime, the original cloud mask is changed from cloudy to clear. The tests used to decide finally whether the pixel is clear or cloudy are given in a flowchart at https://satcorps.larc.nasa.gov/CERES_algorithms.

B. Secondary Cloud Properties and Vertical Structure Estimation 1) Cloud Top Height, Temperature, and Pressure:
The main retrieved parameters for estimating cloud vertical structure are T c and τ from the standard retrievals (i.e., retrievals assuming a single-layer cloud). Using T c , which corresponds to the radiating center of the cloud, it is possible to determine the CEH Z c and CEP P c . They are estimated as the lowest (highest) altitude Z (pressure P) occurring in a vertical sounding, where T (Z or P) = T c . Because the goal of CERES is to provide a data set for relating clouds and radiation and validating models, it is important to provide cloud properties that are both radiatively and physically correct. For example, a single-layer cloud may span several vertical layers in a climate model, but radiate effectively from a specific level between the top and the bottom of the cloud. For ML clouds, the vertical column could have several low layers filled with liquid water droplets and a number of upper layers filled with ice cloud particles. Because the satellite retrieval interprets the cloud in terms of its effective radiance or brightness temperature, the corresponding level could be located anywhere between the top of the ice cloud and the base of the lowest water cloud. While providing the effective radiating temperature or the height of cloud is important, particularly for computing TOA fluxes, the vertical structure of the clouds is also important because it affects the surface and in-atmosphere radiation calculations. Furthermore, accurate characterization of cloud height, thickness, and layering is important for relating clouds to atmospheric dynamics and hydrology. To that end, the vertical boundaries of the cloud are also estimated from CET, COD, phase, CEH, and CEP using various parameterizations. But first CEH and CEP must be determined from the sounding.
Nominally, the MOA sounding is used for T (z), except in the boundary layer. From previous analyses (see [68]), it was found that the detail necessary to fully resolve the temperature profile around the boundary layer inversion, where the tops of many low-level clouds occur, is absent in most model profiles. Furthermore, the clouds in these inversion cases often do not radiate at the air temperature, but at a lower value that is not even observed in the rawinsonde profiles (see [69]) or in aircraft in-situ measurements (see [70]) in the vicinity of the cloud top or the inversion. The result of using model soundings typically causes significant overestimates of CEH for boundary-layer clouds. To account for these and other sources of bias, modifications of the MOA profiles are necessary. These usually take the form of splicing lapserate () corrected temperatures into the MOA sounding.
For Ed2, the boundary-layer temperature profile was calculated with a constant lapse rate where Z o is the surface elevation and T o is the sea surface temperature over ocean and, over land, the 24-h running mean surface (2 m) air temperature from the MOA sounding [38].
In the standard definition of a lapse rate, T o would be the actual air temperature above the surface, and T c would be the temperature of the air at the top of the cloud. In addition, when using (5) to solve for Z c , CET is used instead of the air temperature at cloud top, so b is referred to as the apparent lapse rate. This lapse rate is blended with the MOA profile to construct a new sounding for each tile. Although the mean global regional low CTHs were found to be within ±0.5 km of their CALIPSO counterparts, some significant regional biases over land and ocean were found, particularly over land [41] due to the use of the constant value of b [38]. To minimize the regional biases, Sun-Mack et al. [79] developed seasonal regionally dependent values of b using matched CALIPSO and Aqua MODIS data. These apparent lapse rates are specific for ice-free water, snow-free land, and snow-covered surfaces, respectively, for both daytime and nighttime. They consist of monthly averages at the seasoncenter months: January, April, July, and October. Monthly lapse rates used in Ed4 result from interpolation between two adjacent seasonal centers. In addition to replacing the constant lapse rates with variable ones, Ed4 uses a latitudinally dependent blending approach to produce the modified soundings in the retrievals. A detailed description of this approach and how it is applied in Ed4 processing is provided in [79].
During the daytime, an adjustment is applied to Z c for certain clouds identified as nonopaque cirrus to minimize underestimation of Z c . If MCAT retrievals are available, Z t M is used to alter the standard value of Z c , whenever τ M < 4.5 and τ M > 1.5 * τ over snow-or ice-covered surfaces, or when τ < 4 and P c − P t M > 100 hPa over all other surfaces. The MCAT retrieval is also used whenever the standard methods return a no-retrieval condition and MCAT returns a valid value of τ M . The value of Z t M corresponds more closely to the cloud top rather than the effective radiating height, so it is not used directly to replace Z c . Rather, a new value of Z c is computed using a third-order polynomial.
If τ M < 2, then Otherwise, The coefficients were determined from matched MCAT and VISST CTHs for ice phase clouds. Then, Z c is set equal to Z cM , and a new value of T c is determined from the sounding temperature corresponding to Z cM . The SIST is used to compute r e and τ , while forcing the solution with this new value of T c . This MCAT adjustment to Z c was applied to 7% of the pixels and affected daytime results the most. For NP ocean and land areas, Z cM exceeded the standard CEH by 2.55 and 3.26 km, respectively; whereas over the PO regions, the difference is 1.72 km.

2) Cloud Thickness and Base Height:
To estimate the physical extent of the clouds from the standard retrievals, it is necessary to compute the CTH and CBH, Z t and Z b , respectively. The CERES approach first estimates CTH from CEH and then an estimate of cloud thickness H is subtracted from CTH to yield CBH. Once the top and base heights are found, the CTP and CBP and temperatures, P t and P b and T t and T b , respectively, are taken from the soundings. For Ed2, Z t was determined using a set of parameterizations based on ε c , T c , and τ [40]. Those same parameterizations are used in Ed4, except for the following differences that pertain only to ice clouds.
For most ice clouds, T c corresponds to the temperature at a level Z c that is significantly lower than CTH, Z t . In Ed2, the difference between Z t and Z c was estimated for all ice clouds using parameterizations based on nonopaque cirrus clouds. For clouds having τ > 6, T t was limited to values less than or equal to T c − 2 K and constrained to be no less than the tropopause temperature as described in [40]. This constraint limited the difference, Z t − Z c , to values generally smaller than 0.8 km. It was found that for deep convective clouds the difference typically exceeds 1.0 km [80]. Thus, Minnis et al. [81] developed a new parameterization, based on matched CALIPSO and CERES MODIS data that greatly increase the Z t − Z c differences for ice clouds having τ > 8. A few adjustments were made to that parameterization yielding where μ is the cosine of the VZA. This parameterization was included in the Ed4 code, but its results were mistakenly overwritten by the Ed2 correction later in the code so that the archived results do not reflect the contribution of this new parameterization. It is recommended that users seeking a better estimate of Z t apply (7) to Z c for any CERES MODIS ice cloud data having τ > 6. An exception to that recommendation is for the case of cloud tops detected above the tropopause. The algorithm of [82], introduced in Ed4 to permit cloud tops in the stratosphere, uses infrared data to detect overshooting convective cloud tops. In Ed4, the pixels identified as being overshooting tops are flagged and Z c is estimated by assuming a lapse rate of −8 K·km −1 above the tropopause. CTH is then assumed to be equal to Z c .
The Ed2 methods for estimating ice cloud thickness are based on a few limited studies and tend to grossly underestimate the thickness for optically thick clouds. To more accurately estimate the cloud thickness, new parameterizations were developed using the C3M product [58], which comprises near-nadir, 1-km matched cloud property data from those four A-Train sensors. The matched data from April 2007 were divided into tropical (20 • S-20 • N) and PO (poleward of 50 • ) land and ocean sets. Only those data having an SL with cloud top phase identified as ice by CALIPSO were selected for analysis. Thus, the thickness is obtained by subtracting CBH from CTH. The data were fitted with least squares multiple regression to the following formula to make the initial thickness estimate, where IWP is the ice water path in g·m −2 and the coefficients a i are found in the fitting process. This formula was determined through trial and error. Table II, which lists the coefficients for (8), shows that the thickness is negatively dependent on T c and positively correlated with τ and IWP as might be expected because ice cloud extinction decreases as the temperature drops. For the tropical regions, it was found that a more accurate estimate of the cloud depth can be obtained using the following correction. A third-order polynomial was fitted to the results from (8)  For tropical land, the ice cloud thickness is estimated as follows: For midlatitude regions, the solutions to (8) Finally, H is constrained to be between 0.1 and 19.0 km for all regions. Thickness fits were performed for the midlatitude zones, but it was found that the interpolation process described above produced more accurate results. Fig. 5 shows the results of the parameterization as applied to the 94 678 daytime training samples for April 2007 over tropical ocean. The Ed2 ice cloud thickness parameterization [ Fig. 5(a)] yields a mean thickness of H(Ed2) = 2.41 ± 1.58 km, a significant underestimate compared with H(obs) = 4.77 ± 3.56 km, the mean thickness observed using the C3M data. The difference is especially large for thicknesses exceeding 5 km. This average 2.36-km underestimate is reduced to 0.13 km with the new parameterization [ Fig. 5(b)], although the Ed4 thickness remains too low compared to observed thicknesses above 10 km. In addition, the correlation improves with Ed4 and the RMSD is reduced from 3.44 to 2.02 km for Ed4.
Having Z t and the cloud thickness, the CBH is constrained to be a minimum of 0.1 km above the surface.

3) Multilayer Detection and Retrieval:
To more fully relate the observed radiances to the cloud field, CERES had planned to include overlapping ML clouds in the SSF complement; however, no complete ML retrieval technique was available in 2002, so no ML products are included in the Ed2 SSF. By 2011, when the Ed4 cloud property components were finalized, only a few algorithms were available. One of these, the MCAT-ML method, was developed for CERES and uses the MCAT and the results from the single-layer retrievals [83]. It is used on an experimental basis for Ed4. The ML results are explicitly included within the Ed4 SSF, but they are not used in the SSF 4-category summary, which allows for clear, lower cloud, upper cloud, and overlapped cloud categories for each CERES footprint. The Ed4 overlapped cloud category remains empty and all the retrievals in the upper and lower layer categories are based on the single-layer methods. While this experimental product is not used in the CERES cloud category parameters, it is mentioned here to inform users of its inclusion in the SSF. An in-depth description of the technique and its results, uncertainties, and potential improvements are beyond the scope of this article and will be addressed in a future manuscript.

4) Cloud Water Path:
The Ed4 LWP is computed in the same manner as in Ed2, assuming that CER represents the mean for the entire cloud column. That is, LWP = 4r e τρ w /3Q (11) where Q is the cloud particle extinction coefficient, and ρ w is the density of liquid water. Various studies have shown that a formulation using the adiabatic assumption often yields a more accurate estimate of LWP than (11) [84]- [86]. Although the LWP recorded in the CERES data set is based on (4), the adiabatic LWP can be computed directly from the CERES result by multiplying the LWP by 0.83. The Ed4 ice particle effective radius is based on the spherical equivalent of the ice crystal volume-area ratio, and therefore, the ice water path in Ed4 is computed in the same manner as (11), except that the density of ice replaces ρ w . The retrievals of CER and COD assume that the cloud is a SL composed of only one phase. Thus, the retrievals in the ML ice-over-water cloud systems that are identified as ice cloud have interpreted the radiative contributions of the water phase part of the cloud as being those of an ice cloud. The converse is true for ML clouds classified as liquid.
The resulting LWP or IWP in all cases is thus assumed to represent the total cloud water path TWP. As demonstrated in previous analyses [87], [88], the retrieved TWP in these cases is likely to be significantly biased if the radiative contributions of the phases are not treated explicitly.

5) Other Changes:
In Ed2, the IRW clear-sky brightness temperature, T cs4 , provided to the SSF was used by another CERES subsystem to retrieve the surface skin temperature T s from clear pixels. The retrieved value was included in the SSF. For Ed4, the same retrieval method was applied within the cloud retrieval framework. In (12), ε s4 is the IRW surface emissivity, and B 4 and B −1 4 are the IRW Planck function and its inverse, respectively. The atmospheric attenuation and contribution factor Tr is computed as in [89, eq. (1)]. The values for ε s4 are the same as those used in [29] and [40]. Omission of the surface reflectance term (1 − ε s4 ) in (12) will cause some overestimation of T s , particularly for low-emissivity surfaces.

IV. RESULTS
The differences between the Ed2 and Ed4 cloud amounts were presented and discussed in [29]. Here, the cloud properties from Ed4 are summarized and their differences relative to those from Ed2 are delineated. The averages computed for 2008 results are used as examples in the comparisons. Except for cloud phase fraction, all the reported parameter averages are first computed regionally for each month and year. To obtain global, PO (PO, latitudes poleward of 60 • ), and NP (NP, 60 • S-60 • N) means, the averages are first computed for each 5 • latitude zone by weighting each regional parameter mean by its corresponding cloud fraction. Then the zonal averages, weighted by the cloud fraction mean and the cosine of latitude, are used to compute the means for the specific latitudinal boundaries. Cloud fraction weighting was not used for the Ed2 averages reported in [50], so the results are different in many cases.  Fig. 6(a)]. The only exception appears to be in the heart of the Sahara Desert where few clouds occur. Overall, in NP areas, CFW increased from 30.5% to 39.8%, while in the PO regions it jumped from 31.9% to 41.5%. This contrasts with generally less uniform changes in CFI. The Ed4 marine ice cloud cover [ Fig. 6(d)] is noticeably smaller than its Ed2 counterpart [ Fig. 6(b)] south of 30 • S, while it is larger over Antarctica. The drop over the southern oceans mostly compensates the rise in CFW within the same zones, while the change over the southern PO region represents more cloud detection overall. The same can be said for the increase in CFW over the marine stratus and trade wind zones, where CFI is much the same in the two editions. Over the NP regions,  the mean Ed4 ice cloud amount is 3.2% less than the 26.7% Ed2 average. Overall, the PO averages see a drop of only 1.1% from Ed2 to Ed4. The increases over Antarctica are canceled by the decrease in ice cloud cover in the Ed4 retrievals over the Arctic. These results are summarized in Table III. At night, the 2008 Aqua Ed4 retrievals produced similar increases in water cloud amount [ Fig. 6(e) and (g)], averaging 37.0% compared with 27.6% in the NP regions (Table III) Fig. 6(h)], while the ice cloud amount dropped from 46.8% to 43.1% over the PO regions. The rise in the liquid water clouds is not nearly compensated by decreased ice cloud amounts, indicating that much of the Ed4 gain at night is due mainly to detecting more low clouds than to the changes in the phase algorithm. The differences between the Terra Ed2 and Ed4 cloud phase amounts are similar to those for Aqua, except that the discrepancies between the Aqua and Terra results in Ed4 are generally smaller during the day than their Ed2 counterparts.

A. Cloud Phase
This improved consistency is evident in the daytime trends of ice and water cloud amounts in Fig. 7, except for ice clouds in the PO regions [ Fig. 7(d)]. For Ed2, NP CFW [ Fig. 7(a)] and CFI [ Fig. 7(b)] parallel each other, but CFW(Terra) exceeds CFW(Aqua) by nearly 4% while the converse is true for CFI in the NP areas. Though less parallel, the Ed4 Terra CFW(Terra) tops CFW(Aqua) by 0.0%-0.7% and the Aqua Ed4 ice fraction is only 1%-2% higher than its Terra counterpart. Some discrepancy is expected between the Terra and Aqua mean ice cloud fractions because of somewhat systematic diurnal variations in high cloud cover. High clouds tend to occur more frequently during the afternoon than during midmorning in the NP areas (see [90]). The changes in the Terra calibration and the phase selection algorithm from Ed2 to Ed4 reduced the diurnal difference, perhaps, to a more reasonable value. In the PO zones, CFW(Terra) and CFW(Aqua) are the same until 2007 [ Fig. 7(c)], when the former drops gradually by 0.04 before rapidly rising again in 2015. The PO Ed2 CFW fractions diverge mainly after 2010 reaching a separation of 6% after 2014. This 2010 change in Ed2 CFW(Terra) is mirrored in the sudden rise of 4% in CFI(Terra) such that it matches the Aqua values and remains with them after 2011 [ Fig. 7(d)]. Similarly, the daytime Ed4 CFI(Terra) average [ Fig. 7(d)] mirrors its liquid counterpart by trending upward from a 2.5% difference with Aqua in 2008 up to a 4% discrepancy in 2015. It then drops dramatically after 2015. In general, the Ed2 and Ed4 Aqua daytime means show little trending in the NP regions, while in the extreme latitudes, both Ed4 Aqua CFW and CFI, respectively, tend to drop and rise slightly over the period. The trends in those same parameters are more obvious in the Ed2 results.
Nocturnal cloud phase trends for the NP regions were presented and discussed in [51] and, therefore, are omitted from Fig. 7. Over the PO regions, the nighttime Terra and Aqua Ed2 CFW means [ Fig. 7(e)] are steady and nearly identical over the period. The Ed2 Terra and Aqua PO night ice cloud fractions mostly parallel each other [ Fig. 7(f)], but CFI(Terra) is 2% higher until the latter years, when it converges toward its Aqua counterpart. While Ed4 CFW(Aqua) and CFI(Aqua) are relatively constant, Ed4 CFW(Terra) steadily decreases after 2004 and Ed4 CFI(Terra) rises with more variability until the end of the time period, when the trends sharply reverse. These trends are similar to the NP trends reported in [51], but the variability in cloud phase fractions at all times of day is much greater in the PO than in NP regions. The stronger seasonal and interannual variations are primarily due to the changing geography of the areas in dark and sunlight over the course of the year and the smaller size of the PO averaging area.
The consistent trends in Terra Ed4, decreasing liquid cloud cover and increasing ice cloud fraction, much like those for NP night clouds, are mostly due to the trend in the Terra MODIS C5 8.55-μm channel radiances, which are used at all times of day to help the phase selection for Ed4. Some effects of the Terra C5 6.7-μm calibration are also included in the PO trends. The impact of the phase fraction trend is more dramatic at night because the 8.55-μm channel is used less during the daytime than at night. The sudden changes in the Terra Ed4 parameters after 2015 are due to the replacement of the input C5 radiances by their C6.1 counterparts beginning in February 2016. The running means obscure the discontinuities that occur at the Collection 5-to-6.1 switch. The revised Terra calibrations of the C6.1 data set tend to bring the Terra Ed4 phase fractions back to their 2002 levels. Thus, much of the apparent trending in the Terra Ed4 record is simply the effect of the degrading C5 calibrations. The Terra Ed2 parameters were also impacted since they used the 6.7-μm channel to select phase. As shown in [51] and below, the phase fraction trend impacts the other Terra Ed4 variables retrieved from the Terra MODIS data.
To determine the stability of the phase selection with VZA, 2008 Ed4 CFI and CFW averages were computed as a function of VZA and are plotted in Fig. 8. As expected, there is a noticeable increase in both types of cloudiness during both day and night. Yet, the fraction of ice relative to the total cloud amount (ice plus liquid) hovers around 0.365 during the daytime with a range of only 0.02. At night, the range in the ice fraction relative to the total is 0.012 centered at 0.445. The values show no monotonic behavior at either time of the day. Thus, it can be concluded that the phase selection is, on average, unaffected by VZA.
The revised daytime retrieval algorithms dramatically decreased the fraction of no-retrieval pixels. For example, for 2008, the Ed4 Aqua mean no-retrieval fractions for the globe, NP region, and PO region are 0.7%, 0.6%, and 1.1%, respectively. These values can be compared with 3.0%, 2.1%, and 7.9% for the corresponding Ed2 averages. Similar results were found for Terra. There was little change in the fraction of no-retrievals, ∼0.2%, at night.

B. Cloud Heights
As shown in Fig. 9, the Aqua Ed4 CEH means for 2008 are generally higher than those from Ed2. Fewer areas with liquid CEH, CEHW < 2 km occur in Ed4 [ Fig. 9(c)] than in Ed2 [ Fig. 9(a)] and more of the regional means exceed 4 km in Ed4. Table IV indicates that during daytime, Aqua CEHW(Ed4) averages ∼0.5 km higher than CEHW(Ed2) in the NP regions and 0.34 km higher in the PO areas. At night, the Ed4-Ed2 CEHW differences are much smaller except over the poles, where the mean difference is nearly 0.3 km. The results from Terra (not shown) are very similar. The Ed4 CEHI averages [ Fig. 9(d)] are significantly greater than their Ed2 counterparts [ Fig. 9(b)] and a bit more zonal in distribution. More areas have CEHI > 9 km in Ed4 and regions with CEHI > 11 km, widespread in the Ed4 results, are not seen in the Ed2 data. Table IV shows that during the day and night, the CEHI(Ed4) global averages are almost 1 km greater than CEHI(Ed2) for both Terra and Aqua. The CEHI means are generally greater at night, except for the PO regions where the night is accompanied by a more compressed atmosphere.
The CEH results from 2008 are generally representative of those for all years, as shown in Fig. 10. During the daytime, the Terra and Aqua Ed4 CEHW NP averages are nearly identical between 2002 and 2017, while their Ed2 Terra and Aqua counterparts differ by 0.3 km [ Fig. 10(a)]. The NP daytime Ed4 Aqua CEHI means are roughly 0.2 km higher than the Terra values, but that difference remains fairly steady over the period of record [ Fig. 10(b)]. The same difference is seen in the early years of the Ed2 results, but the Terra CEHI means converge to the Aqua values in the latter half of the period. Over the PO regions during the day [ Fig. 10(c)], the Ed4 CEHW means begin diverging from equivalence  Fig. 10(h)] means, as well as in the NP nocturnal CEHI averages [ Fig. 10(f)]. In the latter case, the Ed4 Terra values exceed the Aqua heights by ∼0.2 km before dropping below the Aqua means by 2011. It is clear that at night, unlike the daytime results, the Terra and Aqua effective heights are more consistent in Ed2 than in Ed4 because the latter relies more on the degrading channels. Like the daytime, the Terra Ed4 nocturnal heights return to their 2002 levels after 2016.  The CEH and CET averages are plotted as a function VZA in Fig. 11 for NP ice and water phase clouds from 2008 Aqua data. During the day, CEHW increases by ∼100 m from 0 • to 65 • , a rise of about 4% [ Fig. 11(a)]. This contrasts with a nocturnal increase of 12%. Similarly, the changes in CEHI with VZA are only 5% (0.46 km) during the day with CEHI, while at night the rise is 0.75 km, or 8%. The CEH is expected to increase with VZA because it is based on the effective radiating temperature of the cloud. CET is a function of the cloud particle concentration in the cloud, the temperature profile, and the optical path. That last item generally increases with VZA, so the level corresponding to the radiating center from the cloud top also rises, yielding a lower value of CET as shown in Fig. 11(b) for all conditions, except for daytime water clouds. For those clouds, CET and CEH are relatively flat and begin rising and dropping, respectively, around VZA = 40 • . Since cloud particle concentrations in ice cloud tops are typically lower than those in water clouds and the radiating center can be 2 km or more below the physical cloud top [81], the change in CEHI with VZA greatly exceeds that for CEHW. The greater rise of CEH at night could be due to the different interpretation of COD using the VIS and IRW channels during daytime and nighttime, respectively, particularly for multilevel thin-ice-over-low stratus cases. The divergence in day-night CETW at 40 • is also probably linked to the differential phase classification of multilevel clouds. Fig. 12 maps out the mean daytime 2008 Ed4 ice CTHs, CTHI, and bases, CBHI, from Ed4. Because CTHI [ Fig. 12(a)] is based on CEHI, its patterns are very similar to, but higher than those in Fig. 9(d). Means exceeding 11 km dominate in the tropics with some regions having CTHI > 12 km. The differences between CTHI and CEHI over the PO regions are much larger, particularly at night (not shown). Overall, the 2008 CTHI means are ∼0.7 km greater than their CEHI counterparts (Table IV), less so over the NP regions and more so over the poles. For NP regions, the CTHW values are ∼0.03 and ∼0.07 km higher than CEHW for NP day and night, respectively. The CTH differences between Terra and Aqua are similar to the CEH differences discussed above.
In the tropics, the ice cloud bases [ Fig. 12(b)] are highest (>9 km) in areas corresponding to the descending branches of the Hadley cell and in many cases where CTHI is highest. The lowest tropical CBHI values (<6 km) are seen in convergence zones and areas where deep convection is common. In the NP regions, the estimated average ice cloud  (Table IV). This difference arises from the limited COD information at night and the greater CTHI values at night. Similarly, CBHW is higher at night, mainly due to the greater CTHW values. The Terra and Aqua CBHW global averages are quite close for both day and night, while the mean ice cloud bases differ by ∼0.2 km.
The CTH and CBH results were not retained in Ed2 so they cannot be compared directly with the Ed4 values. CTP and CBP were saved, so it is possible to measure how much the cloud vertical structure in Ed4 differs from its Ed2 counterpart in pressure units. In 2008 (Table IV), Ed2 CTPI averaged ∼60 hPa more than CTPI from Ed4 during the day globally and in the NP zones. This difference translates to a height difference of ∼1.2 km. Over the PO regions, the CTHI difference is roughly 1.0 km. These values exceed the Ed4-Ed2 CEHI differences, indicating that the Ed4 changes in the CTH algorithms yielded additional gains in CTHI above those attributable to increased CEHI means. At night, the global CTPI differences are, on average, nearly the same as daytime values. Therefore, the CTHI Ed4-Ed2 differences should be around 1.2 km, a value that is still larger than the corresponding CEHI differences.
As shown in Table IV, the Ed4 ice cloud base pressures are all greater than their Ed2 counterparts. For Ed4, the Terra and Aqua CBPI means differ by less than 10 hPa, but together their global averages are 59 and 91 hPa greater than those from Ed2 during the day and night, respectively. The differences are greater in the PO than in the NP areas during the day, while the opposite is true at night. The general day-night differences reflect the thickness differences that arise from the limited COD information at night. Based on having lower CTPI means and greater CBPI averages, it is clear that the ice cloud thickness parameterizations used for Ed4 produced significant increases in the cloud physical depth as might be expected from Fig. 5. Since the liquid cloud thickness parameterization did not change, no comparisons with Ed2 are presented.

C. Cloud Optical Depth, Particle Size, and Water Path
The CERES Aqua mean 2008 daytime CODs are mapped in Fig. 13. The relative patterns in liquid optical depths, CODW, changed little from Ed2 [ Fig. 13(a)] to Ed4 [ Fig. 13(c)], although Ed4 means are higher in the PO regions and lower in some tropical areas. As seen Table V, Ed4 CODW is greater by ∼6% in the NP areas, but nearly a third higher than its Ed2 counterpart in the PO regions. Regionally, the Ed4 ice cloud COD means, CODI, [ Fig. 13(d)] are both larger and smaller than the Ed2 values [ Fig. 13(b)] in the NP areas, but overall are about the same as the Ed2 means (Table V). In the PO regions, CODI(Ed4) is larger everywhere and, on average, is nearly twice its Ed2 counterparts for Aqua and 50% greater than Ed2 for Terra. Over Antarctica and its ice shelf [ Fig. 13(d)], CODI is as much as 16 times greater than it is in Ed2. This dramatic increase in CODI over snow is essentially due to the use of the 1.24-μm channel instead of 2.13 μm for Aqua Ed2 and 1.61 for Terra Ed2, and to the defective absorption LUTs noted earlier. The CODI values over snow may be highly uncertain and biased as result of the input and clear-sky reflectance errors. The limitations of the 1.6-and 2.13-μm channels for retrieving COD, discussed in Section III, are apparent in the Terra-Aqua differences in Ed2 COD means for both ice and water. Trends in the CODs and the differences between Terra and Aqua for both Editions were discussed in [51].
Table V also lists the mean 2008 logarithmic averages of COD, which are two to three times smaller than the linear averages listed in the upper left quadrant. This approach to averaging is used for radiative equivalence because of the exponential relationship between reflected radiance and COD. It is used by the ISCCP [91] and is included for each tile in the CERES SSF.
The 2008 probability distributions of daytime COD are plotted in Fig. 14 in powers-of-2 bins to emphasize the lower ends of the distributions. For water clouds in the northern PO regions [ Fig. 14(a)], the Ed4 COD distribution is shifted to higher values relative to that from Ed2, as expected from the previous discussion. The north PO Ed4 COD mode is in the 8-16 bin compared with the 4-8 bin for Ed2. In addition, 5% of the Ed4 pixels are in the 128-150 bin, which was not available in Ed2. These pixels were originally confined to COD < 128 in Ed2. This very large COD feature is common to all the histograms in Fig. 14. In the northern midlatitudes [ Fig. 14(b)] and the tropics [Fig. 14(c)], the Ed2 and Ed4  histograms are not dramatically different, except that Ed4 detects more clouds having COD <0.25 than found in Ed2. This may be due to using the revised clear-sky reflectance model over ocean [29] and the infrared approach to retrieving τ for cirrus clouds detected only by the thin cirrus test (see [29]). Ice CODs in the Arctic [ Fig. 14(d)] are distributed much like their water counterparts with Ed4 shifted to greater values than found with Ed2. The Ed2 mode is in the 1-2 interval, while it is between 2 and 4 for Ed4. More low-τ pixels are retrieved in the Arctic with Ed2. Ed4 retrieves significantly more CODI < 0.25 than Ed2 in the northern midlatitudes [ Fig. 14(e)] and, especially, the tropics [ Fig. 14(f)], where nearly 12% of the pixels are in the lowest bin. This is due to several factors including the new ocean reflectance model and the use of the 1.38-μm channel and the MCAT in Ed4. Fig. 15 maps the 2008 mean cloud particle effective radii for Ed2 and Ed4. In Fig. 15(a), the mean droplet radius, CERW(Ed2), appears to be slightly larger in many tropical marine areas compared with its Ed4 counterparts [ Fig. 15(c)], but is smaller in the midlatitudes. On average, CERW(Ed4) is 0.2 μm greater than CERW(Ed2) over NP oceans. Over land, CERW(Ed4) exceeds the Ed2 averages in most areas, such that it also averages 0.9 and 0.2 μm greater than Ed2 over NP and PO land, respectively. While the Ed2 CERW averages [ Fig. 15(a)] appear to be slightly larger than their Ed4 counterparts in many PO regions, the PO means (Table V) are the same at the 0.1-μm level. This is mainly because the PO means are calculated using cloud fraction weighting and the Ed2 averages are smaller (greater) than the Ed4 means in the Arctic (Antarctic). Globally, the CERW means are 0.3 and 0.6 μm greater in Ed4 than in Ed2 (Table V) for Aqua and Terra, respectively. As reported in [51], the Terra and Aqua Ed4 means agree as a result of a calibration change in the Terra Ed4 3.7-μm radiances.  Thus, part of the Terra Ed4-Ed2 CERW difference is due to the calibration change. The overall increase in CERW for Aqua Ed4 and part of the Terra Ed4 increase can be explained by the additional cloud cover and liquid phase selection in Ed4.
The Ed2-to-Ed4 changes in ice crystal effective radius, CERI, are more dramatic. The range in regional mean CERI(Ed2) is relatively small [ Fig. 15(b)], ∼15-40 μm with few values on the extremes, and its geographical distribution is relatively flat except for some variations in the tropics. In Fig. 15(d), CERI(Ed4) has a larger range and a quasi-zonal distribution. Despite the distribution differences, CERI(Ed4) for Terra and Aqua combined is, on average, only about 1 μm greater than CERI(Ed2) in the NP regions, but it is ∼9 μm larger in the PO areas (Table V). Globally, CERI(Ed4) is 28.7 μm compared with 26.6 μm for Ed2.
Ice crystal particle size typically decreases with cloud temperature (see [92]) or height [93], and the value of CERI based on the 3.8-μm channel roughly corresponds to the portion of the cloud that contributes to CET. Thus, the retrievals of CERI should be related to CETI, on average. To examine this possibility, the zonal mean values of CERI and CETI were computed for land and ocean separately using July 2013 Aqua retrievals from Ed2 and Ed4. The NP averages are plotted in Fig. 16. The scatterplots for Ed2 [ Fig. 16(a)] show little correlation between CERI and CETI, while those for Ed4 [ Fig. 16(b)] show a significant relationship with squared linear correlation coefficients R 2 of 0.75 and 0.65 over ocean and land, respectively. Using land and ocean together yields R 2 = 0.76. CERI decreases with dropping CETI for Ed4.
These results suggest that the changes in cloud detection, phase selection, ice cloud model, and algorithms made a significant improvement in the retrieval of CERI. An important contributor to the lower CERI values in the tropics is the increased application of the CERI default value, 10 μm, as a result of the additional very thin (τ < 0.25) cirrus detection in Ed4. Thus, some of the correlation may have been inadvertently forced by assuming the default value for thin cirrus. The histograms of 2008 Aqua Ed4 CERI in Fig. 17 show that the values of CERI between 5 and 15 μm account for 37% of the tropical retrievals. This percentage drops to 19% and 12%, respectively, over the northern midlatitude and PO zones. For Ed2, the fractions of CERI in that range are 26%, 22%, and 32% for the tropics, northern midlatitudes, and northern PO regions, respectively, and the histograms are relatively flat (not shown). More use of the default CERI values in Ed4, however, is not the only reason for the correlation observed in Fig. 16(b). It is clear in Fig. 17 that the fraction of CERI between 15 and 30 μm (no default values included) in the tropics is also greater than its counterparts in the northern zones. Similarly, CERI in the midlatitudes occurs more frequently Fig. 17.
Histograms of ice cloud effective radius from the 2008 Aqua Ed4 retrievals. Blue: northern polar zone, red: northern midlatitudes, and gold: tropics. between 15 and 40 μm than over the PO regions indicating that the zonal (cloud temperature) dependence of CERI is a real phenomenon and not due simply to the enhanced number of default values in Ed4.
Since it is proportional to the product of COD and CER, the TWP averages generally follow those for COD and the Ed2-Ed4 comparisons are not plotted here. The 2008 daytime means, though, are listed in the lower right quadrant of Table V. Globally, the mean retrieved LWP increased by 19 g·m −2 or 25% from Ed2 to Ed4, while smaller (∼2%) relative gains in IWP occurred in Ed4. In the NP regions, the mean IWP dropped by ∼5% in Ed4, while LWP increased by 12%. As expected from the Ed4 COD increases, the Ed4 PO LWP and IWP means are roughly double their Ed2 counterparts. If the adiabatic assumption were used, the LWP means would be 17% lower. Because these averages, like those for CER and COD, are based entirely on cloudy pixels alone, they represent only the cloudy portions of the Earth. Lower values would be realized if the zeros corresponding to the clear portions of the Earth were included in the averaging. Fig. 18 plots the distributions of the mean 2008 LWP and IWP, as well as their cloud-fraction-weighted counterparts to show how much condensate overall is estimated to have occurred within each region. In addition to reducing the average water path for liquid water clouds, the weighting [ Fig. 18(b)] results in more regional variability compared with the relatively zonal structure of the unweighted LWP [ Fig. 18(a)]. The minima in IWP [ Fig. 18(c)] within the tropical subsidence areas are enhanced when weighted by the cloud fraction [ Fig. 18(d)]. Similarly, the anomalously high maximum in IWP over eastern Antarctica is highlighted in the product of IWP and ice cloud fraction. Overall, the weighting reduces the mean LWP to 47 and 89 g·m −2 for the NP and PO regions, respectively, decreases of ∼50% relative to the values in Table V. The corresponding IWP averages in Table V drop by ∼63% when weighted by the ice cloud fraction.
The viewing angle variations in the daytime mean optical depth, particle effective radius, and water path are plotted in Fig. 19. Over ocean, liquid cloud τ decreases by 18% between 5 • and 65 • [Fig. 19(a)], while over land it drops   Fig. 19(b)]. Over ocean, droplet effective radius increases by only 6%, while over land it rises by 4% [ Fig. 19(b)]. For ice cloud particles, r e increases by 6% and 9% over ocean and land, respectively. The changes partially complement each other such that the LWP over ocean decreases by 12% from nadir to VZA = 65 • and drops by only 12% over land [ Fig. 19(c)]. The mean IWP decreases by 13% over ocean and land. These changes are likely due to vertical structure and horizontal inhomogeneity in the clouds that is not accounted for in the plane parallel model used in the retrievals [94], [95].
In general, CER retrieved using the 2.1-μm channel, CER 2 , is significantly greater than its 3.8-μm counterpart for both phases. The results for CER 2 are discussed in the following section. They are representative of this parameter for Ed4.

V. DISCUSSION
To better understand the accuracies of the Ed2 cloud retrievals, they were compared with a variety of other data sets including some that could be considered ground truth. Those studies are summarized in [50]. Additional comparisons have been performed for the Ed4 retrievals to assess the uncertainties for several parameters. Some of those comparisons are summarized in [97]. Other studies are reviewed below along with new comparisons to field observations and to similar results from the mean October 2008 MAST MYOD08 Collection 6 products [96], [98]. The Ed4-MAST comparisons are summarized in Table VI for all the considered cloud parameters. In Part II [47], cloud phase, CTH, CBH, and thin ice cloud microphysical properties are examined in detail using CloudSat and CALIPSO data. The reader is referred to that article for a more comprehensive assessment of those parameters.

A. Cloud Phase
The cloud top phase is a difficult parameter to verify in a statistically and globally meaningful way without the aid of a satellite-based reference data set. To that end, [47] used phase determinations from CALIPSO measurements to assess the cloud-top phase assigned to the MODIS pixels in the Ed4 process. Even with an objective reference like CALIPSO, such comparisons are fraught with complications due to the nature of real clouds, particularly those having temperatures in the supercooled droplet range of −40 • C to 0 • C. Overlapped ice and water clouds, mixed phase tops, and mixtures of cloud types in a given MODIS pixel all exacerbate the evaluation of a given phase retrieval. Part II covers some of those aspects of validation. The impact of errors in phase or any of the parameters ultimately depends on the use of the data.
Nevertheless, it is informative to understand how the Ed4 results differ from other long-term satellite data sets. The water and ice phase fractions from four different cloud data sets, including CERES Ed4 Aqua, were compared in [45] for daytime October 2008 data. The Ed4 liquid cloud fraction was found to be higher than that MAST Aqua retrievals by 0.05-0.10 at nearly all latitudes. Similarly, it is greater than its counterpart from the PATMOS-X Advanced Very-High Resolution Radiometer (AVHRR) data set [99] by an even larger amount poleward of 45 • latitude, but less than the SatCORPS AVHRR means. Conversely, the Ed4 ice fraction is less than the PATMOS-X values by 0.00-0.06 in the tropics, but by up to 0.30 or more in the midlatitudes and PO regions. However, the Ed4 ice fractions exceed both the MAST Aqua and SatCORPS AVHRR amounts by 0.02-0.10 or more at all latitudes.
Some of the discrepancies in the MAST-Ed4 phases may lie in sampling. To examine this, the Ed4 and MAST cloud phase fractions for October 2008 Aqua data are summarized in Table VI and plotted in Fig. 20, which shows the average water and ice cloud fractions from CERES and computed from the MYD08 products [100] using the sum of the overcast, partly cloudy fractions, and undetermined fractions. The undetermined pixels are assumed to be liquid water [98].  (Table VI) over the globe. Since the MAST total daytime cloud fraction is typically ∼0.01 greater than that for Ed4 [29], the MAST phase retrievals are not provided for ∼19% of the total identified cloud population. This is not unexpected as the MAST retrieval cloud fraction should be noticeably less than the cloud mask amount as some of the cloud properties were not estimated for various reasons, for example, some pixels were reclassified as clear in the retrieval process [101]. The Ed4 cloud properties and phase are also retrieved using the 3.8-μm channel instead of the 2.1-μm channel, which is used for the MAST standard product. The fraction of MAST clouds having a retrieval using 3.8 μm, shown in parentheses in Table VI, increases by 0.024 in the NP regions but is roughly the same as the standard retrieval fraction in the PO areas.
Using 8 months of matched MODIS and CALIPSO data, Yost et al. [47] found that during the day Ed4 overestimates the topmost liquid cloud phase by an average of 0.03 and underestimates the ice fraction by 0.04-0.06, but by less than 0.10 at all latitudes. The phase discrepancies are primarily due to the interpretation of optically thin ice clouds over a layer of water clouds as being liquid phase [47]. Thus, much of the phase difference between Ed4 and other techniques likely arises from the differences in interpreting radiances from ML clouds. At night, the Ed4 ice fraction agrees with the CALIPSO retrievals, but liquid cloud fraction is underestimated, as the phase selection is more responsive to the upper layer cloud signal.

B. Cloud Vertical Structure
To validate the regional lapse rate method for low clouds described earlier, Sun-Mack et al. [79] compared retrieved low cloud heights with CALIPSO for 2 months using fully sampled MODIS data. Using only single-layer, water-phase clouds having a CTH below 4 km according to CALIPSO, the study also showed how the Ed4 algorithm performed relative to several other techniques. The agreement between CALIPSO and Ed4 is better and the distribution of the differences more random than they are for the Ed2 results.
Xi et al. [102] compared SL stratus CTH and CBH determined from a cloud radar at Graciosa Island, Azores with matched Terra and Aqua Ed4 heights. The cloud top and base height differences are 0.06 +0.25 and −0.07 +0.28 km, respectively, during the daytime. At night, the corresponding differences are 0.14 +0.34 and −0.04 +0.37 km. The mean CTH biases are similar to those seen in the CALIPSO comparisons in [47, Table II], while the mean base height differences are similar to those from CALIPSO during the day and much better than those at night [47]. The standard deviations of the differences are quite small, less than half those seen in the global CALIPSO comparisons.
Minnis et al. [45] compared the averaged October 2008 cloud top pressures determined from several satellites using various algorithms. These include CERES Ed4 Aqua MODIS, SatCORPS-A1 NOAA-18 AVHRR, MAST C6 Aqua MYD08, PATMOS-X NOAA-18 AVHRR [99], ISCCP-D2 [91], and CLARA-A2 [103]. Between 15 • S and 30 • N, the Ed4 mean zonal CTP is lowest among the various results, but is generally close to the mean for all algorithms at other latitudes. It is significantly less than its MYD08 counterpart at all latitudes indicating that the CERES cloud heights, on average, are greater than the MAST cloud heights. This is confirmed in Fig. 21, which shows the zonal mean MAST and Ed4 Aqua CTH s for October 2008. The MAST heights (gray) are nearly the same for both day (dashed) and night (solid), while the Ed4 daytime heights (black dotted) are significantly less than those at night (black solid), as expected from Table IV and the results in [47]. The differences between the Ed4 and MAST daytime heights range from near-zero at 50 • S to 1.5 km at 30 • N. At night, the zonal differences vary from 0.5 to 3.0 km. In the NP areas, the mean day and night differences are 0.65 and 1.77 km, respectively. The corresponding PO differences are 1.05 and 1.35 km. The MAST algorithm [96] uses the same channels for both day and night retrievals and, hence, produces more consistency between the day and night CTHs. Since Ed4 uses a combination of VIS and several infrared channels during the day and only infrared channels at night, it yields the kind of day-night differences seen in Fig. 21. The sources of the differences are discussed in Part II. Combining the day and night results in Fig. 21 and similar results found for April 2008 (not shown) corroborates the Ed4-MAST CTP differences found in [45].
A global comparison of Ed4 CBHs and single-layer cloud thicknesses with the CALIPSO data was conducted by Yost et al. [47]. Details of the comparisons and results are provided therein.

C. Cloud Optical Depth, Particle Size, and Water Path
The Ed4 changes in COD retrievals over snow, the use of the rough ice crystal model, and the additional cloudiness detected in Ed4 likely had minimal effect on the properties of the Ed2 liquid water clouds that have already been compared with various surface and airborne instruments as summarized in [50]. Thus, it is expected that similar comparisons of the Ed4 results with the same surface observations used for Ed2 would yield comparable agreement for liquid water clouds over snow-free areas. A number of studies performed to validate the Ed4 retrievals against surface and aircraft measurements can be added to or contrasted with the Ed2 comparisons. These include both MODIS and GEOsat retrievals (see [95], [104], [105]) that used the Ed4 algorithms or their equivalents. Initial Ed4 comparisons performed over other regions provide entirely new information.
1) Liquid Clouds: Xi et al. [102] compared Terra and Aqua Ed4 retrievals to 19 months of observations taken at the Azores AMF. The Ed4 CERW means over a 30-km 2 area are 1.3 μm greater than the ARM retrievals (12.8 μm), while the Ed4 mean LWP is 13.5 g·m −2 less than its ARM counterpart (114.2 g·m −2 ) due to its smaller optical depth (9.6 versus 13.7). The differences are reduced by 50% when the Ed4 averages are computed using only the MODIS pixel nearest the AMF site. The 10% differences between the ARM and CM LWP and CERW retrievals are within the uncertainties of the ARM LWP (∼20 g·m −2 ) and CERW (∼10%) retrievals; however, the 30% difference in optical depth is significant. Possible reasons contributing to this discrepancy are increased sensitivities in optical depth from surface retrievals when τ ∼ 10 and an island effect. Similar results were found using the MAST product [106].
Those results contrast with Painemal et al. [86] who found that the 4-month Ed4 LWP averages computed using the adiabatic assumption exceed the AMSR2 means by 1.3% and 9.3% when using CERW and CERW 2 for computing LWP, respectively, over the northeastern Pacific for overcast AMSR2 pixels. The Ed4 CERW-based LWP means also matched ship-based ARM-AMF MWR retrievals for all types of marine stratus sky conditions, nearly clear, partly cloudy, mostly cloudy, and overcast. The average LWP computed with CERW 2 for all the cloud conditions is between the AMSR2 means and LWP based on CERW. Thus, CERW provides a reliable estimate even for broken cloudy scenes, while CERW 2 is likely overestimated in the same conditions.
To further illustrate the value of the adiabatic assumption for estimating LWP, the Ed4 LWP data from overcast scenes used in [86] were averaged within a 10-km radius of the ship location at each Terra and Aqua overpass time and matched with the mean ship LWP for an hour centered on the overpass time (Fig. 22). A scatter plot of the matched ship and Ed4 LWP shows that the standard LWP, based on the homogeneous cloud assumption (6) exceeds the ARM-AMF 3-channel MWR LWP by 18.2 g·m −2 [ Fig. 22(a)]. Given the mean MWR LWP, 63.8 g·m −2 , the bias is 29%, in the same direction as the Ed2 comparisons over land [49]. However, using the adiabatic approximation, the satellite LWP bias reduces to only 4.5 g·m −2 or 7% [ Fig. 22(b)]. The standard deviation of the differences is also reduced when using the adiabatic approach. Absent an island effect, these differences and those reported in [86] are likely to be more representative of the retrieval accuracy over ocean than the Azores comparisons. These results confirm that the adiabatic approach offers better performance, at least, for marine stratus LWP estimation.
Dong et al. [107] compared Ed4 MODIS overcast stratus properties with retrievals based on surface measurements at the ARM surface site at Barrow, Alaska, for both snow-free and snow-covered conditions. For snow-free cases, they found that the Ed4 CERW is 1.8 ± 3.5 μm (14%) larger, on average, than that found from the surface site. However, the average COD is 0.4 (5%) smaller than that from the surface retrievals. The corresponding mean LWP bias relative to the ARM retrieval is −0.6 ± 35.0 g·m −2 (−1%). With the adiabatic assumption, the bias would be −12 g·m −2 (−15%). For snow-covered cases, the Ed4 r e and τ biases are 1.6 ± 4.0 (14%) and 0.6 ± 4.3 μm (8%), respectively. For these cases, the mean LWP is 5.5 ± 47.2 g·m −2 (9%) too high. Here, the adiabatic assumption yields a similar bias, −5.8 g·m −2 , but with a different sign.
Wood et al. [108] compared the in-situ measurements of r e using a wide-spectrum cloud probe for a flight through drizzling marine stratocumulus with matched Ed4-like retrievals from a GOES imager. They and Glienke et al. [109] found that although the mean GOES CERW of 23.6 μm was only 1-2 μm larger than the in-situ average, the probability distributions of the two data sets are markedly different with the satellite having a narrow range (8-32 μm) compared with the in-situ data that ranged from 6 to 60 μm. The study confirms the existence of clouds having large (>20 μm) values of CERW found in satellite retrievals over some open ocean areas. It also suggests that the satellite retrievals tend to overestimate r e for many drizzling clouds, while at the same time underestimating r e for heavily drizzling clouds because the retrieval LUT is limited to CERW < 32 μm. These results must be taken into consideration in designing future retrievals. This impact of drizzle is discussed further in a later subsection.
Given the larger values of r e retrieved using the 2.1-μm channel (see Section V-C.3), it can be concluded that LWP computed using those values will yield significant overestimates for most of the stratus cases considered. The one exception is for the Azores matchup, where using r e2 to calculate LWP would reduce the difference between the Ed4 and ARM MWR retrievals from 13.7 to 2.1 g·m −2 .
2) Ice Clouds: Validation of the Ed4 ice cloud properties has been conducted using CALIPSO, surface radar, and aircraft in-situ measurements. Part II of this article discusses the CALIPSO comparisons, which only consider thin cirrus. Other comparisons are recapped below.
Ice cloud particle size and habit tend to vary with altitude or temperature, while the value of CERI generally corresponds to the upper 1-3 optical depths of the cloud [110]. Thus, for thin ice clouds, CERI should generally be representative of the entire cloud. Mace et al. [111] found that CERI from CERES Ed2 underestimated the radar retrieval in cirrus with τ < 2.2 by 2.9 μm. Given the average increase of 1.9 μm in CERI relative to Ed2 for NP land (Table V), it is expected that the Ed4 CERI would be closer than its Ed2 predecessor for those cases.
For thick ice clouds, it is expected that CERI will be less than r e for the entire cloud because crystal size tends to increase with temperature. The temperature (height) dependence of CERI is illustrated in Fig. 23, which shows r e averages from in-situ CERI averages (solid symbols) plotted as a function of normalized distance below cloud top using measurements taken in ice clouds during the 2013 Studies of Emissions and Atmospheric Composition, Clouds and Climate Coupling by Regional Surveys field program [112]. The data were taken during flights on August 8,12,14,16,21,23, and 30 and on September 2,4,11,13,18,and 21. The normalized altitude, or depth below cloud top, is Z n = (Z t − Z a ) / H , where Z a is the aircraft altitude, and Z b and H are taken from the Ed4 retrievals when cloud top and base height information is not available from the aircraft. Most of the sampled clouds are from deep convective anvils with cloud thickness ranging from 2 to 9 km and averaging ∼5.5 km. The in-situ effective radii r ea were computed using measurements from the 2-D-S probe [113] onboard the NASA DC-8 and SPEC IC9 Lear Jet. Each solid point represents the average r ea in one of the seven equal Z n intervals for flights in ice clouds at times within 20 min of the Terra or Aqua overpass.
The profile in Fig. 23 shows that r ea is smaller at the top and increases to a point corresponding to a depth of ∼0.65 in the cloud before decreasing to cloud base (Z n = 1). Retrievals from Terra and Aqua MODIS pixels were matched to the aircraft flight track and used to compute averages corresponding to r ea for each normalized flight level. The differences between the Ed4 and in-situ means are shown as open symbols in Fig. 23. Those differences are positive above the 0.2 level and drop to as low as −12 μm before increasing again toward cloud base. The absolute values of the differences may be somewhat uncertain as the 2-D-S probe does not measure particle radius smaller than 5 μm and likely underestimates concentrations for sizes smaller than 10 μm [113]. The vertical changes, however, are probably representative of the relative vertical profiles.
Waliser et al. [20] found very good correspondence between average IWP distributions from CloudSat and Ed2 retrievals. The level of agreement may be fortuitous given the shortcomings of the two different approaches to estimating IWP, particularly in light of the findings in [88], as discussed below. Nevertheless, similar results are expected for Ed4 since the mean IWP for NP areas is nearly identical to that of Ed2 (Table V). Over the PO regions, Ed4 IWP is about double the Ed2 value, but in [20], Ed2 underestimates IWP relative to CloudSat in the PO regions. The Ed2 and Ed4 IWP and LWP values are actually the estimates of the TWP since it is assumed that the cloud is single phase and single-layered. If liquid water is in the column below the ice cloud, its contributions to the reflectance will be interpreted as if it were ice.
Smith [88] compared TWP derived for thick ice-topped clouds from radar and radiometer data with Ed4 IWP values over the ARM Oklahoma site. He found excellent agreement between the surface-based TWP and a surrogate Ed4 IWP from the GOES data for IWP < 0.5 kg·m −2 . For larger values, IWP underestimated TWP by 15% for TWP = 1.0 kg·m −2 up to 25% for TWP = 4.0 kg·m −2 . Thus, for IWP > 0.5 kg·m −2 , it is expected that the true IWP is greater than the retrieved value.
Ice water paths for deep convective clouds determined from the Next-Generation Radar (NEXRAD) data over the south central United States were compared in [114] with Ed4 Terra and Aqua retrievals and with Ed4-like retrievals from GOES imager data. It was determined that the GOES and Ed4 retrievals underestimated the NEXRAD IWP by 4 ± 61% and 9 ± 59%, respectively, for the thick anvil portions of the convective systems, which had a mean IWP of 2.24 kg·m −2 . For the stratiform rain parts of those systems (mean IWP = 3.0 kg·m −2 ), the corresponding underestimations increased to 25 ± 67% and 36 ± 71%. Part of the stratiform IWP differences is due to the Ed4 cap of 150 on the possible retrieved COD. Some portion of it is likely due to using CERI from the 3.8-μm channel, which is not likely representative of the total column effective radius size (e.g., Fig. 23). The comparisons also do not account for any liquid cloud water that may have been in the bottom layers of the clouds. Short of increasing the maximum COD retrieval and using a shorter wavelength, for example, 1.6 μm, for the CERI retrieval or applying a correction such as that suggested in [88], it will not be possible to obtain much improvement in the IWP retrieval for very deep convective cloud systems. Nevertheless, the current approach appears to be reasonable for all but the thickest cloud systems.
3) Comparisons With Other Satellite Retrievals: CERES Ed4 CODs and particle sizes from Aqua MODIS were compared with the results from several retrieval methods in [45] using the October 2008 data. The Ed4 optical depths for liquid and ice clouds exceed those from PATMOS-X AVHRR retrievals by 1-2 and 1-3, respectively, in the NP regions, but are considerably smaller than the PATMOS-X means in the PO regions. Zonally, the Ed4 mean CERW is 0-2 μm less than its PATMOS-X counterpart in the NP areas and 0-5 μm less over Antarctica. For CERI, the Ed4 retrievals are within ±2 μm of those from PATMOS-X in the tropics, but elsewhere CERI from Ed4 is typically 15 μm greater than its PATMOS-X counterpart. The latter decreases toward the poles from a latitude of 30 • , while the former increases.
The geographical distributions of Ed4 and MAST COD regional averages are very similar, except over the PO regions (not shown). Both the Ed4 and MAST algorithms retrieve an optical depth when retrieving r e at a given wavelength; however, τ values returned from the Ed4 2.1-μm r e retrieval are not saved and cannot be compared with their MAST counterparts. Because the 2.1-μm retrieval, COD 2 , is the standard for MAST, it is shown in Table VI along with the 3.8-μm retrieval. The latter is listed in parentheses. For NP liquid clouds, CODW from Ed4 falls between the MAST CODW 2 and CODW from MAST differing by an average of 3%. The PO Ed4 CODW averages ∼40% less than its MAST counterparts. The mean Ed4 CODI in the NP areas is 26% greater that the corresponding MAST averages for either channel, while it is between the two MAST means over the PO regions, differing by ±7%. In addition to the sampling differences, the CODI differences may be due to the smaller MAST asymmetry factors, which are nearly constant at 0.75. Fig. 24 plots the global distribution of the mean CERW retrieved by CERES and MAST for the October 2008 Aqua MODIS data as determined using channels 2 and 5. The global patterns are very similar for all the four panels, but there are some notable differences. CERW from MAST [ Fig. 24(b)] generally tends to be slightly smaller than Ed4 CERW [ Fig. 24(a)] over land, while differences between Ed4 and MAST are both positive and negative over ocean. Overall, the Ed4 retrieval is 0.3 μm greater than its MAST counterpart over the NP zones and 0.4 μm smaller over the PO regions, mainly due to the large droplets seen by MAST over Antarctica. Nevertheless, the two 3.8-μm results are very close despite the differences in the cloud fractions that were sampled.
The MAST CERW 2 [ Fig. 24(d)] is noticeably smaller than CERW 2 from CERES [ Fig. 24(c)] over ocean areas, but is greater over land areas. Overall, the NP Ed4 averages are 1.6 μm greater and 0.6 μm smaller than the MAST means (Table VI) over ocean and land, respectively, but larger by 1.2 μm overall. Over the PO regions, they are nearly identical. The PO agreement may be coincidental in that most of the CERES pixels are sampled over snow-free land and ocean areas, while the MAST averages also include extensive snow-covered areas. The larger 2.1-μm droplets from CERES could be due to various factors including the different reflectance parameterizations and the model reflectance LUT interpolation.
Zhang et al. [106] found that the Ed4-MAST CERW 2 differences over the Azores grow with increasing droplet size with near-zero differences around r e2 = 11 μm. Since clouds having r e > 15 μm often contain significant drizzle droplets, the drizzle droplets may cause a overestimate of CERW 2 and less so for CERW because the drizzle droplets increase the absorption of the cloud at the retrieval wavelengths depending on the fraction of drizzle-sized droplets in the cloud (see [115]). Zhang and Platnick [94] found that the drizzle-generated overestimate for one set of model runs is only 2 μm at the largest values of r e2 . However, the droplet distribution widens as the amount of drizzle increases [109], which could produce additional overestimates. Thus, if the variance used in the Mie scattering calculations to create the reflectance and absorption LUTs for the r e retrievals does not represent the spread of the actual droplet distribution, the retrieved value will likely be biased. Arduini et al. [116] demonstrated that increasing the effective variance beyond the value of 0.1, used for Ed4, for calculations of 3.7-μm reflectances tends to decrease the retrieved CERW. It is likely that the differences will be even larger for 2.1-μm reflectances. Table VI indicates that on average, Ed4 CERW and CERW 2 exceed their MAST counterparts by 18% and 24%, respectively. This difference is also reflected in the ratios of CERI 2 /CERI, which are 1.36 and 1.42, respectively, for MAST and Ed4. Over land, the two retrievals at 2.1 μm are closer with CERI 2 means of 29.6 and 33.2 μm for MAST  and Ed4, respectively, a difference of only 12%. There are minimal ocean-land differences between the Ed4 CERI means and for the MAST CERI and CERI 2 means. The differences are likely found in the sampled populations, the retrieval parameterizations, and differences in the ice cloud models used in the retrievals. Given the ocean-land differences for both Ed4 CERI 2 and CERW 2 , potential errors in the treatment of the Ed4 clear ocean reflectance must also be considered as a possible source of the Ed4-MAST biases and should be investigated further.
The MAST and Ed4 mean cloud water path means are determined from the product of τ and r e and therefore should reflect the Ed4-MAST differences. This is generally true in Table VI, where the MAST NP averages of LWP and LWP 2 bracket the Ed4 mean LWP of 0.88 kg·m −2 and MAST NP means for IWP and IWP 2 are both less than the Ed4 mean IWP. Over the PO areas, the MAST LWP averages exceed those from Ed4. Conversely, the MAST IWP averages are both less than the Ed4 mean because the MAST ice particle sizes are both less than theirEd4 counterparts.
The VZA variations in mean τ , r e , and cloud water path seen in Fig. 19 have been observed in previous analyses of both MAST and CERES Ed2 retrievals. Heck et al. [117] found that Ed2 CER tended to increase with VZA for both liquid and ice clouds, especially over land, while COD was relatively constant for water clouds and decreased with VZA for ice clouds. The Ed4 results in Fig. 19 show COD decreasing with VZA for both phases. The differences may be due to the reduced sampling for the earlier study, which used only 1 month of data. While not providing actual means, Maddux et al. [118] determined that CERW and CERI from MAST C5 retrievals, on average, typically increased from nadir to the end of the scan (VZA ∼ 66 • ) by ∼2 and ∼4 μm, respectively, in the NP regions. In that same analysis, CODW and CODI dropped by ∼3 and ∼5, respectively, over the same VZA range. In Fig. 19(b), the increases in Ed4 CERW and CERI are not as great, but the drop in optical depths is comparable. A more focused study in [119] used only warm liquid clouds over ocean to show that the MAST C5 retrievals of r eW rise by 2%-4% out to VZA = 50 • , similar to the results in Fig. 19(b). The corresponding decrease in τ W for that data set is less than that in Fig. 10(a) over ocean with the result that the CERES LWP diminishes by ∼12% compared with ∼5% for the warm clouds analyzed in [119]. It is not clear that the discrepancies in the VZA changes are due to sampling or algorithmic differences because the Ed4 averages include all clouds identified as liquid, which includes warm, supercooled, and some cirrus over liquid clouds. Horvath et al. [119] also found that although they decrease overall with VZA, τ W and LWP increase with VZA for backward scattering angles. It is expected that the Ed4 values behave in a similar fashion. Horvath et al. [119] and Liang and Di Girolamo [120] provide explanations for the observed variations.

D. Surface Skin Temperature
The new surface skin temperature variable added to the Ed4 SSF is based on a retrieval that does not account for the surface-reflected IRW radiance and uses surface emissivities that differ from those retrieved or used by other groups (see [121], [122]). To explore how these and other error sources impact the skin temperature retrieval, the results are compared with other T s estimates.  Table VII, show that T s averages over ocean from CERES are 0.2 and 0.5 K less than those from GEOS-5 during day and night, respectively, for the NP area. The corresponding differences over PO oceans are The differences are slightly negative at night for both the regions. The negative nocturnal differences over land tend to balance the positive differences during the day, resulting in differences of less than 1.0 K overall. The larger magnitude of the nocturnal differences over ocean compared with the daytime values is probably due to more cloud contamination in the retrievals at night [47]. While the sea surface temperatures appear quite accurate relative to GEOS-5, the apparent overestimation of T s over many land areas in some regions could arise for a variety of reasons. The GEOS-5 averages include cloudy scenes, which would typically reduce the mean during the day and increase it during the night over land or snow. In addition, the GEOS-5 is a reanalysis and, although it assimilates sea surface temperature observations, the land skin temperature seems to be a prognostic variable and is unaffected by the assimilated satellite radiances [53]. While the differences in Fig. 25 could be due to model sampling and some GEOS-5 uncertainties in formulating the land skin temperatures, other error sources are also in play. Kato et al. [123], in a more direct comparison of T s (Ed4) using only matched clear-sky GEOS-5 and other reanalysis data sets, found that T s (Ed4) is generally about 5 K warmer than all the reanalyses over dry areas such as North and South Africa. Furthermore, compared with CERES observations, the outgoing LW radiation computed with T s (Ed4) tends to be significantly overestimated in many dry areas, due, in large part, to the surface emissivities that were used in the retrieval. In dry areas, the CERES IRW surface emissivities have been found to be lower than their counterparts from [124] and others and, hence, produce noticeably higher surface temperatures.

VI. CONCLUSION
A significant number of changes have been made to various components of the CERES cloud retrieval algorithms in Ed4. These have been accompanied by the addition of new (to CERES) retrieval parameters that though not essential to the downstream CERES processing have been included for experimental purposes. The notable changes to the algorithms for the essential properties are summarized as follows.
1) As reported here and in [51], CERES calibration adjustments to several Terra MODIS channels resulted in better long-term Terra-Aqua consistency in cloud phase, effective radius, effective height, and optical depth, except for those parameters affected by calibration drifts in one or more infrared channels. 2) The revised cloud mask for Ed4 [29] allowed for the retrieval of more clouds with small optical depths than found in Ed2. 3) An updated set of algorithms produced improved cloud phase detection as indicated in [47], particularly for single-phase clouds. 4) The use of the 1.24-μm channel for retrieving COD over snow produces better Terra-Aqua consistency than resulted from Ed2. Thick water CODs appear to be reasonable, but COD is likely overestimated for clouds having τ < 8 due to a parameterization error and surface albedo input. Correction of the reflectance parameterization, better surface albedos, and, possibly, use of a longer wavelength, such as 1.6 μm, should lead to more accurate CODs over snow at all COD values. 5) Use of a roughened ice crystal reflectance model and the CO 2 channel in the MCAT increased the effective heights for ice clouds by ∼1 km relative to the Ed2 heights, a marked improvement. Nevertheless, the daytime cloud heights are lower than those at night and, by implication, inconsistent with the results that would be obtained with an infrared technique as confirmed in Part II [47]. Future editions of the CERES algorithms should consider using a model, such as the two-habit model of 125] and [126], to achieve consistency between the VIS and infrared COD retrievals. 6) Correction for the radiating depth of ice clouds was applied and mistakenly overwritten in Ed4. It is a simple correction and can be applied by users to obtain a more accurate height for ice clouds having τ > 6. Its impact is quantified in [47]. 7) A new ice cloud thickness parameterization resulted in better characterization of the vertical extent of single-layer ice clouds and is verified in [47]. 8) Use of regionally dependent apparent boundary-layer lapse rates has improved the liquid water CTHs over both land and water surfaces compared with Ed2 [51]. Comparisons of averages with those of other retrieval methods and of point data with field measurements indicate that the Ed4 retrievals are quite reasonable. In addition to the differences in algorithms, the discrepancies among the data sets can be attributed to sampling differences, interpretation of multilayered clouds, and treatment of no-retrieval pixels.
The experimental cloud parameters added to the CERES SSF product include an ML retrieval that is evaluated elsewhere, multispectral particle size retrievals, and surface skin temperature retrievals. The multispectral size retrievals are less reliable than expected, mainly as a result of uncertainties in the input parameters and coding and parameterizations used in their retrievals. Additional research to fully examine the sources for those errors and to correct them is needed to facilitate their upgrade to standard parameters in the next edition of CERES cloud properties. Surface skin temperatures in arid regions appear to be overestimated due to underestimates of the surface emissivities in those areas and the lack of a surface reflectance term in the retrieval parameterization. These sources of bias can be easily remedied in future editions using a different set of emissivities (see [124]) and the parameterization of [89]. For Ed4, these new parameters can be used confidently for certain conditions, as outlined in this article.
The changes to the essential CERES cloud parameters represent a marked improvement relative to the Ed2 retrievals and to "ground truth" reference data sets. Those alterations to the retrieval algorithms have contributed significantly to higher accuracies in the fluxes determined from the CERES Terra and Aqua radiances [26], [33], [36] and lend greater confidence to studies of cloud and radiation interactions using the CERES data. Remaining shortcomings in the retrieved parameters, both new and old, have been identified and will be addressed in future editions of the CERES cloud algorithms or, in some cases, can be rectified to some degree by interested users of the data.

ACKNOWLEDGMENT
The authors would like to thank N. Loeb, S. Kato, D. Doelling, and D. Kratz for discussions leading to the development and assessment of the Ed4 cloud products. They would like to thank J. K. Ayers, R. Brown, E. Heckert, and T. Chee for computing, Internet, and graphical support of these algorithms. They would like to thank K. Bedka for providing his overshooting top algorithm. The CERES and C3M data are available at the NASA LaRC Atmospheric Sciences Data Center (https://eosweb.larc.nasa.gov). The MODIS Monthly Global Product data were obtained at http:// modis-atmos.gsfc.nasa.gov/MOD08_M3/. The authors also thank P. Lawson for the use of the 2-D-S data. The DC-8 2-D-S data were taken from the NASA Airborne Science Data for Atmospheric Composition webpage at https:// www-air.larc.nasa.gov/cgi-bin/ArcView/seac4rs.