A Review of Ice Cloud Optical Property Models for Passive Satellite Remote Sensing

The current wealth of spaceborne passive and active measurements from ultraviolet to the infrared wavelengths provides an unprecedented opportunity to construct ice cloud bulk optical property models that lead to consistent ice cloud property retrievals across multiple sensors and platforms. To infer the microphysical and radiative properties of ice clouds from these satellite measurements, the general approach is to assume an ice cloud optical property model that implicitly assumes the habit (shape) and size distributions of the ice particles in these clouds. The assumption is that this ice optical property model will be adequate for global retrievals. In this review paper, we first summarize the key optical properties of individual particles and then the bulk radiative properties of their ensemble, followed by a review of the ice cloud models developed for application to satellite remote sensing. We illustrate that the random orientation condition assumed for ice particles is arguably justified for passive remote sensing applications based on radiometric measurements. The focus of the present discussion is on the ice models used by the Moderate Resolution Imaging Spectroradiometer (MODIS) and the Clouds and Earth’s Radiant Energy System (CERES) science teams. In addition, we briefly review the ice cloud models adopted by the Polarization and Directionality of the Earth’s Reflectance (POLDER) and the Himawari-8 Advanced Himawari Imager (AHI) for ice cloud retrievals. We find that both the MODIS Collection 6 ice model and the CERES two-habit model result in spectrally consistent retrievals.


Introduction
Although it has been recognized that ice clouds significantly influence the Earth's climate system through their regulation of the atmospheric radiative energy budget [1][2][3][4], definitive knowledge and fundamental understanding about the microphysical and optical properties of ice clouds are still limited, particularly from the perspective of climate studies based on modeling and remote sensing capabilities.For example, Zhang et al. [5] showed substantial discrepancies of longwave, shortwave, and net cloud radiative forcings based on climate model runs in comparison with three major satellite products.While the discrepancies are caused in part by unrealistic cloud occurrence frequencies and heights, there are substantial uncertainties related to unrealistic assumptions of cloud microphysical and radiative properties.In cloud classification, numerous types of clouds are defined [6]; however, in the study of the transfer of radiation in a cloudy atmosphere, we usually classify clouds into three major categories in terms of cloud-base height: low clouds with cloud-base heights below 3 km; mid-level clouds with cloud-base heights between 3 and 6 km; and high clouds with cloud-base heights above 6 km.An alternative classification makes use of cloud thermodynamic phase as follows: water clouds, consisting of liquid phase droplets; ice clouds, consisting of exclusively non-spherical ice particles; and mixed-phase clouds that consist of various mixtures of ice and liquid phase particles.The total global cloud coverage lies within 65%-70%, of which the ice cloud fraction is approximately 20%-30% globally, but with higher fractional coverage (30%-40%) in the tropics [7][8][9][10][11].
The current level of understanding of ice cloud optical and microphysical properties is much lower than that of liquid water clouds.For example, Waliser et al. [12] found that there are pronounced discrepancies in the global average of ice water path (IWP) obtained from simulations based on 19 general circulation models (GCMs).The best way to achieve an improved representation of ice cloud optical and microphysical characteristics over the entire globe is to use satellite retrievals of cloud microphysical properties to diagnose model outputs and validate ice cloud schemes [12].One such parameterization was given by Elsaesser et al. [13].Furthermore, it should be pointed out that significant improvements have been made to infer IWP by using active remote sensing techniques (e.g., [14,15]).
Satellite measurements are used to infer global cloud microphysical and radiative properties, particularly cloud effective particle size and optical thickness as well as other properties such as cloud top height/pressure/temperature and cloud thermodynamic phase.This review paper focuses on optical thickness and the effective particle size; the product of these two parameters in conjunction with a constant scaling factor provides an estimate of IWP.
As noted above, spaceborne sensors are currently providing both passive and active measurements of global clouds.The passive imager measurements span the ultraviolet (UV) through the infrared (IR).Current polar-orbiting sensors include the Advanced Very High Resolution Radiometer (AVHRR) [16], the Moderate Resolution Imaging Spectroradiometer (MODIS) [17,18], and the Visible Infrared Imaging Radiometer Suite (VIIRS) [19].Furthermore, microwave imagers provide unprecedented measurements of precipitation (rain and snow); in particular, the Global Precipitation Measurement (GPM) [20] Microwave Imager (GMI) covers a latitudinal measurement area up to 65 • based on a specially designed orbit [21].MODIS provides cloud products on both the NASA Earth Observing System (EOS) Terra and Aqua platforms.AVHRR spectral measurements are more limited than newer imagers, ranging from 0.65 to 12 µm but have been used to infer cloud microphysical and optical properties [16,22,23].MODIS spectral measurements range from 0.4 to 14.2 µm, and both the Terra and Aqua MODIS sensors are still providing measurements long past their intended lifetime.The current MODIS cloud products were described by Platnick et al. [24].The VIIRS takes measurements in 22 imaging and radiometric channels ranging from 0.41 to 12.5 µm.Passive microwave measurements are useful to estimate not only liquid hydrometeor properties (e.g., [25]) but also ice counterparts (e.g., [26]).Additionally, the sub-millimeter spectral region is known to be sensitive to IWP and the bulk ice particle size, and Baran et al. [27] described the bulk single-scattering properties for ice clouds at frequencies between 89 and 874 GHz.
The NASA A-Train constellation has provided years of data from a variety of sensors.We also note that the French satellite Polarization and Directionality of the Earth's Reflectances (POLDER) [28] takes polarimetric measurements that augment the information content provided by other imagers.Additionally, active lidar depolarization measurements of ice clouds from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) aboard the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) [29] platform have proven to be invaluable for comparisons with products from passive sensors (e.g., [30]).CALIOP has been used to infer ice water content [31].Use of multiple sensors [32,33] combined CloudSat, CALIPSO, and MODIS data to infer ice cloud properties.Further evaluations of different A-Train ice cloud retrievals may be found in the studies of Deng et al. [14] and Matrosov [34].
In addition to these sensors that take narrowband measurements, there are several broadband and hyperspectral IR sensors to note.On the EOS Aqua platform is the Atmospheric Infrared Sensor (AIRS), which takes measurements at 2378 IR channels.AIRS data have been used to infer cloud properties [35].The Suomi NPP and JPSS platforms have the Cross-track Interferometer Sounder (CrIS) [36].Historically, the AVHRR and High-resolution Infrared Radiometer Sounder (HIRS) sensors have flown in pairs on NOAA operational polar-orbiting platforms since about 1978.The HIRS measurements continue at the time of this writing, although hyperspectral IR sensors are supplanting these IR measurements.The European Metop A/B satellite platforms include the AVHRR, HIRS, and Infrared Atmospheric Sounding Interferometer (IASI) sensors [37].The IASI is a hyperspectral sounder that provides a continuous spectrum in the same wavelength range as CrIS and AIRS but without any spectral gaps.
In addition to this variety of sensors on polar-orbiting platforms, the new generation of geostationary satellites is providing spectral data at higher spectral and spatial resolution.We call special attention to cloud retrieval algorithm comparisons (e.g., [38]) involving the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) onboard the Meteosat Second Generation (MSG) platform.In addition to working towards consistency between various polar-orbiting platforms, there is considerable work to be done to improve consistency between polar-orbiting and geostationary satellite products (e.g., [39]).A new generation of advanced geostationary imagers including Himawari-8 Advanced Himawari Imager (AHI) [40] and the GOES-16 Advanced Baseline Imager (ABI) [41] offers unprecedented data.In particular, both Himawari-8 AHI and GOES-16 ABI have 16 channels with resolutions of 0.5-2 km, depending on specific channels.
While several methods have been developed to infer cloud optical thickness and effective particle size, in the following discussion we focus on two approaches.The first approach, referred to as the Nakajima-King algorithm [42], is based on a pair of visible (VIS) and shortwave-infrared (SWIR) channels.The background of inferring cloud properties from spaceborne observations in the solar spectrum can be traced to earlier studies [43][44][45][46].The second approach is based on the satellite observations made in the 8-12 µm IR window region, known as the split window technique.In the earliest application of thermal IR radiation, Inoue [47] presented the inference of temperature and effective emissivity of semitransparent cirrus clouds using a pair of two IR window channels.The split-window technique has been extended to the retrieval of ice cloud properties, the details of which have been discussed in [48][49][50][51], and references cited therein.For both methods, it is critical to use an appropriate ice cloud bulk optical property model to generate look-up tables (LUTs) for the shortwave channel-based retrievals or in-line radiative transfer simulations for the thermal IR based retrievals.
In this paper, we review the evolution of ice cloud optical property models with a focus on the applications to retrievals conducted by the MODIS [14,15] and the Clouds and Earth's Radiant Energy System (CERES) [52] science teams.In addition, we also briefly summarize the ice cloud optical property models used for ice cloud property retrievals based on observations by the POLDER [28] and the Japanese geostationary satellite Himawari-8 AHI [53].The remainder of this paper is organized as follows.In Section 2, we present the theoretical basis pertaining to ice cloud optical property models for satellite remote sensing.In Section 3, we first review the early development of ice cloud optical property models, followed by a discussion of the ice cloud models used by several satellite science teams as well as the assessment of the performance of the aforementioned models.Finally, a summary is given in Section 4.

Theoretical Basis for the Optical Properties of Ice Clouds
Before we discuss the detailed optical properties of ice clouds, it is necessary to introduce several key optical and radiative properties involved in the transfer of radiation in a cloudy atmosphere.For atmospheric remote sensing applications based on the spectral or narrowband observations between 0.5 and 15 µm, the specification of a radiative quantity is essentially based on optical analysis (i.e., the relevant measurements are made by optical instruments).For this purpose, it is sufficient to use the Stokes vector to fully specify the radiation field, given the fact that two radiation beams with the same Stokes vector cannot be distinguished through optical analysis [54].To define the Stokes vector, let us consider a plane electromagnetic wave.The electric vector of the wave can be decomposed into perpendicular and parallel components (E // and E ⊥ ) with respect to a given reference plane, as illustrated in Figure 1.In conjunction with the decomposition of the electric field vector, the Stokes vector, given by I = (I, Q, U, V) T where the superscript T indicates matrix transpose, are defined as follows: where the asterisk indicates the complex conjugation and i = √ −1.The four elements of the Stokes vector are usually referred to as the Stokes parameters that have the dimension of irradiance.Among the Stokes parameters, I indicates the intensity of the radiation, Q and U indicate the state of linear polarization, and V indicates the state of circular polarization.An extensively used quantity is the degree of linear polarization (DLP) defined in the form [55,56] In the case of U = 0, DLP is also defined in the form [57] wave can be decomposed into perpendicular and parallel components ( E // and E ^) with respect to a given reference plane, as illustrated in Figure 1.In conjunction with the decomposition of the electric field vector, the Stokes vector, given by where the superscript T indicates matrix transpose, are defined as follows: , where the asterisk indicates the complex conjugation and . The four elements of the Stokes vector are usually referred to as the Stokes parameters that have the dimension of irradiance.Among the Stokes parameters, I indicates the intensity of the radiation, Q and U indicate the state of linear polarization, and V indicates the state of circular polarization.An extensively used quantity is the degree of linear polarization (DLP) defined in the form [55,56] . (5) In the case of , DLP is also defined in the form [57] .Figure 2 shows the geometric configuration for the scattering of a plane electromagnetic wave propagating along the y-axis by a non-spherical particle.Both the incident and scattered electric vectors can be decomposed into parallel and perpendicular components with respect to the scattering plane, a plane containing the incident and scattering directions.The decomposition is described as follows: Figure 2 shows the geometric configuration for the scattering of a plane electromagnetic wave propagating along the y-axis by a non-spherical particle.Both the incident and scattered electric vectors can be decomposed into parallel and perpendicular components with respect to the scattering plane, a plane containing the incident and scattering directions.The decomposition is described as follows: where the superscripts s and i indicate scattered and incident waves, respectively.k is referred to as the modified wavenumber, given by k = 2π/λ and λ is the wavelength.r is the distance between the scattering particle and the observer in the radiation or far-field zone (r >>λ).S i, i=1−4 are the elements of the amplitude scattering matrix.In the case of a non-spherical particle, without axial symmetry, the amplitude scattering matrix depends on not only the scattering zenith angle (θ) but also the azimuthal angle (φ) of the scattering plane illustrated in Figure 2.
where the superscripts s and i indicate scattered and incident waves, respectively.k is referred to as the modified wavenumber, given by and  is the wavelength.r is the distance between the scattering particle and the observer in the radiation or far-field zone (r >>  ). are the elements of the amplitude scattering matrix.In the case of a non-spherical particle, without axial symmetry, the amplitude scattering matrix depends on not only the scattering zenith angle ( q ) but also the azimuthal angle ( f ) of the scattering plane illustrated in Figure 2. In conjunction with the transformation of the incident electric field to the scattered electric field in Equation ( 7), the corresponding transformation for the Stokes parameters is in the form , (8) where matrix F is often referred to as the scattering Mueller matrix, which is a 4 × 4 matrix with elements that can be derived from the amplitude scattering matrix in Equation (7).For example, F 11 is given as follows: The mathematical expressions for the other elements of matrix F are described in van de Hulst [58].For many applications, e.g., in atmospheric remote sensing implementations and plane-parallel radiative transfer simulations, the scattering medium is assumed to be an ensemble of particles that are randomly oriented with equal numbers of mirror-imaging positions.The complexity greatly increases if the particles are given non-random orientations [59,60].In this paper, we assume random orientations for the particles with equal numbers of mirror-imaging positions unless specifically defined otherwise.In this case, matrix F can be simplified in the following form [58,59,61]: In conjunction with the transformation of the incident electric field to the scattered electric field in Equation ( 7), the corresponding transformation for the Stokes parameters is in the form where matrix F is often referred to as the scattering Mueller matrix, which is a 4 × 4 matrix with elements that can be derived from the amplitude scattering matrix in Equation (7).For example, F 11 is given as follows: The mathematical expressions for the other elements of matrix F are described in van de Hulst [58].For many applications, e.g., in atmospheric remote sensing implementations and plane-parallel radiative transfer simulations, the scattering medium is assumed to be an ensemble of particles that are randomly oriented with equal numbers of mirror-imaging positions.The complexity greatly increases if the particles are given non-random orientations [59,60].In this paper, we assume random orientations for the particles with equal numbers of mirror-imaging positions unless specifically defined otherwise.In this case, matrix F can be simplified in the following form [58,59,61]: Furthermore, matrix F in Equation ( 10) can be expressed in a normalized form as follows: P 11 P 12 0 0 P 12 P 22 0 0 0 0 P 33 −P 43 0 0 where matrix P is generally referred to as the normalized scattering phase matrix, and C sca is the scattering cross section given by In addition to the scattering cross section, one must also know the absorption cross section C abs to quantify the scattering particle's capability of absorbing the incident radiation.The sum of the absorption and scattering cross sections provides the extinction cross section: For practical light-scattering computation, the extinction cross section is usually computed by using the amplitude scattering elements in the exact forward direction via the optical theorem [58].The ratio of the scattering cross section to the extinction cross section is referred to as the single-scattering albedo, given by The term P 11 (θ) (usually called the phase function) is normalized by: In many practical applications, we do not consider the detailed angular variation of the phase function.Instead, we introduce a parameter called the asymmetry factor g to quantify the degree of the deviation of the scattered light from the incident direction, which is given by The above optical properties are referred to as the single-scattering properties for an individual particle.In practice, we need to consider an ensemble of scattering particles, which is referred to as a polydisperse system.The particle size distribution, n(D) where D is the characteristic dimension specifying the size of the corresponding particle, has units of 1/volume/length (e.g., 1/cm 3 /um).The physical meaning of the size distribution is that n(D)dD is the number of particles per unit volume with particle characteristic dimensions between D − dD/2 and D + dD/2.For a given size distribution, the bulk extinction coefficient with units of per length is given by To describe wavelength-dependent attenuation in a cloud given a dispersion of particles, a useful dimensionless quantity is the optical thickness, defined by where z = ∞ indicates the top of the atmosphere (TOA).
The above single-scattering properties are fundamental to the equation governing the transfer of radiation in a cloudy atmosphere.For example, in the case of plane-parallel atmosphere, the governing equation for radiance is given by where I indicates the total (direct plus diffuse) radiance, B indicates Planck's function associated with thermal emission, and cos Θ is the cosine of the scattering angle, given by Forward radiative transfer simulations based on Equation ( 19) are essential to the implementation of a remote sensing retrieval algorithm based on either LUTs or an inline radiative transfer solver.As the single-scattering properties are fundamental quantities in the radiative transfer equation, the use of the most representative bulk optical properties of cloud particles is crucial for simulating the radiation field accurately.Note that the choice of a representative bulk optical property model depends on the application.The model chosen for inference of global cloud properties based on a satellite imager may be quite different than the model chosen for a ground-based or aircraft-based sensor, or for a polarimeter where other elements of the phase matrix must be considered.
Early studies (e.g., [62]) assumed ice particles to be spherical, which is no longer a valid assumption given the numerous measurements of actual ice particles (e.g., [63][64][65][66][67]).Over the years, the research community has made significant progress towards solving for the optical properties of non-spherical particles from very small to very large size parameters, and from the UV through the far-infrared spectrum, from either theoretical or practical computational perspectives.Numerous light-scattering computational methods have been developed, including the separation of variables method for spheroids [68], the finite-difference time domain (FDTD) technique [69][70][71], the discrete dipole approximation (DDA) [72][73][74], the T-matrix method [75][76][77], and the geometric optical method and its improved counterparts [78][79][80][81][82][83][84][85].Because light scattering by a non-spherical particle is a complicated subject, we do not review detailed light-scattering computational methods.An interested reader is referred to Mishchenko et al. [59,86] for non-spherical particles in general and Liou and Yang [4] for ice particles more specifically.In particular, we note that light scattering in the entire size parameter region may be effectively simulated by a synergetic combination of efficient implementations of the invariant imbedding T-matrix method [87] and the physical-optical method [83,84] in conjunction with recent improvements on its computational efficiency [88,89].

Ice Particle Shape and Orientation
Ice particles in the atmosphere are exclusively non-spherical particles (e.g., [63,64,90,91]).Beyond the spherical approximation for non-spherical ice particles, one of the earliest non-spherical models considered horizontally oriented cylinders [92], as illustrated in Figure 3. Specifically, ice particles are assumed to be long circular cylinders that are oriented horizontally in relation to their rotational axes.In Figure 3a, the cylinder axes point along the same direction, whereas, in Figure 3b, the cylinders are horizontally oriented in a random manner.While analytical light-scattering analysis for an infinitely long circular cylinder [93,94] is applied to approximate the single-scattering properties of a finite circular cylinder through an appropriate scaling procedure, Liou's ice cloud models shown in Figure 3 were used by others (e.g., [95,96]) who showed that the significant differences in the single-scattering properties of spheres and cylinders resulted in pronounced difference in radiance simulations.
Ice particles in the atmosphere are exclusively non-spherical particles (e.g., [63,64,90,91]).Beyond the spherical approximation for non-spherical ice particles, one of the earliest non-spherical models considered horizontally oriented cylinders [92], as illustrated in Figure 3. Specifically, ice particles are assumed to be long circular cylinders that are oriented horizontally in relation to their rotational axes.In Figure 3a, the cylinder axes point along the same direction, whereas, in Figure 3b, the cylinders are horizontally oriented in a random manner.While analytical light-scattering analysis for an infinitely long circular cylinder [93,94] is applied to approximate the single-scattering properties of a finite circular cylinder through an appropriate scaling procedure, Liou's ice cloud models shown in Figure 3 were used by others (e.g., [95,96]) who showed that the significant differences in the single-scattering properties of spheres and cylinders resulted in pronounced difference in radiance simulations.The scattering of light by oriented particles is quite complicated.As an example, Figure 4 shows the two-dimensional phase functions of an ensemble of hexagonal plates (Figure 4, upper panel) and hexagonal columns (Figure 4, lower panel) assuming a horizontally random orientation condition.The plates are randomly oriented in terms of their rotations around the asymmetry axes that are aligned vertically and columns are horizontally oriented with two parallel faces perpendicular to the vertical direction (known as Parry oriented columns).The physical-geometric optical computational program developed by Sun et al. [89] is used for this simulation.
As shown in Figure 4, the angular distribution of scattered energy associated with oriented ice particles is intriguing.For example, note that in the 2-D contour of the phase function of a horizontally oriented plate, there are two bright spots that correspond to a scattering polar angle of approximately 22°; these are the basic optical mechanism of sundogs shown in Figure 5a [4,97].The scattering of light by oriented particles is quite complicated.As an example, Figure 4 shows the two-dimensional phase functions of an ensemble of hexagonal plates (Figure 4 As shown in Figure 4, the angular distribution of scattered energy associated with oriented ice particles is intriguing.For example, note that in the 2-D contour of the phase function of a horizontally oriented plate, there are two bright spots that correspond to a scattering polar angle of approximately 22 • ; these are the basic optical mechanism of sundogs shown in Figure 5a [4,97]. Atmosphere 2018, 9, x FOR PEER REVIEW 9 of 32       The false color image in Figure 5a is produced by combining the radiance simulations at three wavelengths in a false color composite (Red:Green:Blue (RGB) form of Red: I 1.064um × 3, Green: I 0.532um × 0.7, and Blue: I 0.355um × 2).Only 7.9% of plates are assumed to be horizontally oriented, whereas 92.1% of the particles are 3-D randomly oriented.The simulated radiances are mapped to the coordinates of fish-eye camera image.In Figure 5a, note the presence of an optical phenomenon called a sun pillar in addition to the sundogs [98].The 2-D phase function of the Parry oriented columns shown in Figure 4 has complexities that may be captured with appropriate sensor measurements under single scattering conditions.The false color image in Figure 5b is produced by combining the radiance simulations at three wavelengths (Red: 610 nm, Green: 532 nm, and Blue: 460 nm) by assuming an optically thin ice cloud consisting of 100% Parry ice columns.Multiple scattering effects are excluded in the simulations.Note that rare optical phenomena such as helic, subhelic, and Parry infralateral arcs ( [4], and the references cited therein) are reproduced in Figure 5b.
The preceding discussion leads to a question: Is it a valid simplification to assume that ice particles are randomly oriented in the atmosphere for passive remote sensing applications?Fortunately, measurements made by CALIOP aboard CALIPSO provide unprecedented observations to quantify atmospheric ice particle orientations.In particular, Hu et al. [99] showed that the relationship between the layer-averaged attenuated backscatter (γ) and layer-averaged depolarization ratio (δ) from the CALIPSO measurements can be used to effectively discriminate between horizontally oriented and 3-D randomly oriented ice particles.Here, γ and δ are defined as [99] where β b,⊥ and β b,// are the perpendicular and parallel components of the attenuated backscatter observed at CALIOP 0.532-µm channel, respectively.Hu et al. [99] explained that, in the γ-δ relation shown schematically in Figure 6a, the segment corresponding to strong backscatter and low depolarization ratio is associated with horizontally oriented ice particles, whereas the segment corresponding to weak backscatter and large depolarization ratio is associated with 3-D randomly oriented ice particles.In following Hu et al. [99] to investigate the aforementioned theoretical speculation, Figure 6b shows the CALIOP data for the γ-δ relation as observed during April 2007.At that time, the CALIOP lidar was in a nadir-viewing configuration and the specular reflection by horizontally oriented ice particles was so substantial that the lidar pointing configuration was soon adjusted to point away from nadir. Figure 6c is similar to Figure 6b except the data were obtained during April 2009 when CALIOP lidar pointed 3º off the nadir.In the case of Figure 6c, the specular reflection is not observed by the lidar receiver and the segment of large γ and small δ is missing in Figure 6c.The comparison of panels 6b and c of Figure 6 unequivocally confirms the existence of oriented ice particles.With observations made by CALIOP and the Infrared Imaging Radiometer (IIR), both of which are on the CALIPSO platform, Saito et al. [100] derived the percentage of oriented ice particles.The results are shown in Figure 7.The fraction of oriented ice particles is less than 1% at cloud temperatures warmer than −30 °C but further decreases with temperature.These results are essentially consistent with those by Zhou et al. [101] and Chepfer et al. [102].Based on the current evidence, it is quite reasonable to assume ice particles to be randomly oriented for passive remote sensing applications that do not involve polarimetric measurements.Note that polarimetric microwave measurements are sensitive to hydrometeor orientations [26], and the assumption may introduce errors in ice cloud retrievals based on the measurements.With the assumption of the random orientation condition for the scattering particles, the single-scattering properties and the equation of radiative transfer are significantly simplified.Specifically, the phase matrix depends With observations made by CALIOP and the Infrared Imaging Radiometer (IIR), both of which are on the CALIPSO platform, Saito et al. [100] derived the percentage of oriented ice particles.The results are shown in Figure 7.The fraction of oriented ice particles is less than 1% at cloud temperatures warmer than −30 • C but further decreases with temperature.These results are essentially consistent with those by Zhou et al. [101] and Chepfer et al. [102].Based on the current evidence, it is quite reasonable to assume ice particles to be randomly oriented for passive remote sensing applications that do not involve polarimetric measurements.Note that polarimetric microwave measurements are sensitive to hydrometeor orientations [26], and the assumption may introduce errors in ice cloud retrievals based on the measurements.With the assumption of the random orientation condition for the scattering particles, the single-scattering properties and the equation of radiative transfer are significantly simplified.Specifically, the phase matrix depends only on the scattering polar angle and is independent of the azimuthal angle of the scattering plane.The phase matrix is in a diagonal-block form as in Equation ( 10) (note: the diagonal-block form implies that the particles have an equal number of mirror-imaging orientations in addition to the random orientation condition).Furthermore, a scalar extinction coefficient can be used instead of the extinction matrix for radiative transfer simulations.
Atmosphere 2018, 9, x FOR PEER REVIEW 12 of 32 only on the scattering polar angle and is independent of the azimuthal angle of the scattering plane.The phase matrix is in a diagonal-block form as in Equation ( 10) (note: the diagonal-block form implies that the particles have an equal number of mirror-imaging orientations in addition to the random orientation condition).Furthermore, a scalar extinction coefficient can be used instead of the extinction matrix for radiative transfer simulations.

Bulk Optical Properties for Passive Remote Sensing Applications
For passive remote sensing applications, the first group of ice particle models discussed here includes three models developed by the MODIS science team.The MODIS science team applies the Nakajima-King bi-spectral method [42] to simultaneously retrieve cloud optical thickness and effective particle size.As illustrated in Figure 8, the reflectivity correlation between a pair of MODIS channels contain information about cloud optical thickness and effective particle size.The information content stems from the fact that, in a conservative scattering channel (e.g., a channel centered at 0.85 μm), ice cloud reflectivity is insensitive to the effective particle size while the reflectivity at weakly absorbing shortwave-infrared channel (e.g., a channel centered at 2.13 μm) depends on both optical thickness and effective particle size.In particular, the optical thickness isolines are orthogonal to the effective particle size isolines at large optical thicknesses, allowing straightforward retrievals of the two parameters given a sufficiently high optical thickness.In the literature, there are several forms to define ice cloud effective particle size, but the following definition gives consistent definitions in both the cases of water and ice clouds [103,104]: where V total and A proj,total are the total volume and projected area, respectively, given by A proj,total = A proj (r)n(r)dr

Bulk Optical Properties for Passive Remote Sensing Applications
For passive remote sensing applications, the first group of ice particle models discussed here includes three models developed by the MODIS science team.The MODIS science team applies the Nakajima-King bi-spectral method [42] to simultaneously retrieve cloud optical thickness and effective particle size.As illustrated in Figure 8, the reflectivity correlation between a pair of MODIS channels contain information about cloud optical thickness and effective particle size.The information content stems from the fact that, in a conservative scattering channel (e.g., a channel centered at 0.85 µm), ice cloud reflectivity is insensitive to the effective particle size while the reflectivity at weakly absorbing shortwave-infrared channel (e.g., a channel centered at 2.13 µm) depends on both optical thickness and effective particle size.In particular, the optical thickness isolines are orthogonal to the effective particle size isolines at large optical thicknesses, allowing straightforward retrievals of the two parameters given a sufficiently high optical thickness.In the literature, there are several forms to define ice cloud effective particle size, but the following definition gives consistent definitions in both the cases of water and ice clouds [103,104]: where V total and A proj,total are the total volume and projected area, respectively, given by A proj,total = r max r min A proj (r)n(r)dr.
In Equations ( 24) and ( 25), r indicates half of the maximum characteristic dimension (i.e., the radius in the case of a sphere) of an individual particle.In the case of a water cloud consisting of spherical droplets, Equation ( 23 The definition given in Equation ( 26) is exactly the same as that in Hansen and Travis [103] for water cloud particles assumed to be spherical.Thus, Equation ( 23) is a rational extension of the definition introduced by Hansen and Travis to non-spherical ice cloud particles.Another advantage of the definition by Equation ( 23) is that the ice water path (IWP) can be expressed in the form where ρ ice = 0.916g/cm 3 is the mass density of bulk ice, and <Q ext > is the mean extinction efficiency given by Atmosphere 2018, 9, x FOR PEER REVIEW 13 of 32 In Equations ( 24) and ( 25), r indicates half of the maximum characteristic dimension (i.e., the radius in the case of a sphere) of an individual particle.In the case of a water cloud consisting of spherical droplets, Equation ( 23 The definition given in Equation ( 26) is exactly the same as that in Hansen and Travis [103] for water cloud particles assumed to be spherical.Thus, Equation ( 23) is a rational extension of the definition introduced by Hansen and Travis to non-spherical ice cloud particles.Another advantage of the definition by Equation ( 23) is that the ice water path (IWP) can be expressed in the form where is the mass density of bulk ice, and < Qext > is the mean extinction efficiency given by Figure 8. Schematic diagram illustrating the principle of the Nakajima-King bi-spectral method [42] for simultaneous retrieval of cloud optical thickness and effective radius.The look-up table shown in the diagram is generated with a solar zenith of 60°, a viewing zenith angle of 45°, and a relative azimuthal angle of 105° for a severely roughened ice particle model.Courtesy of Yi Wang (Texas A&M University).
A complementary approach for inferring the optical properties when the optical thickness is lower than the optimal range for a reflectance-based method in the VIS-SWIR bands is to use A complementary approach for inferring the optical properties when the optical thickness is lower than the optimal range for a reflectance-based method in the VIS-SWIR bands is to use split-window measurements in infrared channels at 11 µm and 12 µm (e.g., [105]).These measurements are provided by most weather satellite imagers.Figure 9 shows the look-up table constructed by split-window measurements from MODIS bands 31 and 32 (11 µm and 12 µm, respectively).
To estimate an optimal solution of a state vector (x) in the look-up table generated by the Nakajima-King bi-spectral method or split-window IR measurements, we use a maximum likelihood method [106] by minimizing a cost function (J) given by where y, F(x), and S y are the measurement vector, the forward model, and the measurement error covariance matrix, respectively.The errors of the retrieved state and measurement vectors are assumed to follow normal distributions.Furthermore, we assume that there is no correlation between measurement errors.The forward model employs the correlated-K method [107] to account for gas absorption and the discrete ordinate radiative transfer (DISORT) model [108] or the adding-doubling method [109,110] for radiance simulations.The retrieval process is schematically shown in Figure 9. Based on an initial state vector, x = (τ 1 , D e1 ) T , we obtain an updated state vector, x = (τ 2 , D e2 ) T from a Jacobian matrix in conjunction with the Newton method.Then, we proceed this process iteratively to finally obtain the optimal state vector, x = (τ 3 , D e3 ) T , which meets a criterion of the cost function (J ~m, where m is the number of elements of the measurement vector) [106].In this study, retrievals of ice clouds are conducted over ocean, and the atmospheric profiles and ocean surface temperature are provided by ECMWF reanalysis data.
In the computation of MODIS particle models (see Figure 10), the definition of maximum dimension (D max ) depends on particle habits.D max is defined as the particle height for a hexagonal (solid and hollow) column, diameter of basal facets for a hexagonal plate, and the longest distance between two vertices for a column aggregate, bullet rosette, or droxtal.The mixed scattering properties averaged over the particle size distribution are called "bulk" scattering properties.The precise definitions of the bulk properties were summarized by Baum et al. [111].The other important information necessary for the derivation of the bulk scattering properties is knowledge of the particle size distribution (PSD).Historically, the PSDs have been defined for this effort in two ways for ice cloud remote sensing by passive imagers: using a well-defined gamma distribution or in-situ measurements.The obvious benefit of using a gamma distribution is that one has total control over the size limits and width of the PSD.The use of actual in-situ data is much more involved than the use of analytical PSDs, but the relevant bulk single-scattering models can be assessed for different dynamical regimes (e.g., low updraft velocity such as synoptic cirrus versus tropical anvils formed in a convective dynamical environment).Different approaches may be adopted for active sensors, which can sense much larger particles and infer much higher values of IWP, and millimeter to sub-millimeter measurements [27].
Among the aforementioned three ice particle models employed for the calculation of ice cloud bulk scattering properties, the first model (Figure 10a) was employed in MODIS Collection 4 [112,113].The particle habit fractions are fixed in either a small or large particle regime as follows: (1) 50% four-element bullet rosettes, 25% hexagonal plates, and 25% hollow columns for D max < 70 µm; and (2) 30% aggregates, 30% bullet rosettes, 20% hexagonal plates, and 20% hollow columns for D max > 70 µm.The bulk properties are calculated by integrating the single scattering properties over only twelve PSDs, each of which is discretized in terms of five coarse size bins (10-30, 3-70, 7-170, 170-430, and 430-1070 µm).For the 12 PSDs, the effective particle radii are 6.7, 9.6, 14.5, 19.7, 25.0, 28.0, 32.0, 40.0, 47.4, 48.6, 50.7, and 58.9 µm.The effect of the particle-size bin resolution was studied by Nasiri et al. [114], who basically found that a finer discretization was necessary.Liu et al. [115] found that a high bin resolution in terms of using the kernel technique can uncover more nuanced information of the bulk microphysical and optical properties of ice particles.Given the coarseness of the MODIS Collection 4 model, effort was made to increase the library of ice particle single scattering properties (e.g., [116,117]).Based on improvements in both the library of single scattering properties and the use of microphysical PSD data archived from various field campaigns, a second model [118,119], illustrated in Figure 10b, was developed for MODIS Collection 5.The model was more complex because single scattering properties were available in 45 size bins, with habit fractions as follows: (1) 100% droxtals for D max < 60 µm; (2) 50% solid columns, 35% hexagonal plates, and 15% six-element bullet rosettes for 60 < D max < 1000 µm; (3) 45% solid columns, 45% hollow columns, and 10% column aggregates for 1000 < D max < 2000 µm; and (4) 97% bullet rosettes and 3% column aggregates for D max > 2000 µm [118].These particle habit mixtures were determined to best explain 1117 samples of airborne in-situ particle size distributions, ice water contents, median mass diameters, and pictures of ice particles captured during the field campaigns.While bulk single-scattering properties were computed for each of the 1117 particle size distributions, the single-scattering properties used in retrievals were based on an average of 10-20 size distributions distributed around each of 18 target effective diameters [119].The MODIS Collection 5 ice model was also used for AIRS version 6 retrievals [35].
While the MODIS Collection 5 models were based on a set of 1117 different PSDs, the measurements detailed in [90] greatly increased the number of available PSDs by an order of magnitude, with resulting ice cloud bulk single-scattering models detailed in [120] for visible through far-infrared wavelengths.
The most current model used by the MODIS science team is that for the operational MODIS Collection 6 [18].As shown in Figure 10c, the MODIS Collection 6 model is rather straightforward in terms of particle habit: it consists solely of aggregates of severely roughened hexagonal column particles, regardless of particle maximum dimension, and gamma particle size distributions are assumed with a fixed effective variance parameter (v e f f = 0.1).The development of the MODIS collection 6 ice cloud model was largely motivated by the consistency between shortwave-based and IR-based retrievals [18,121].
In addition to the MODIS science team, the CERES science team retrieves cloud properties as part of its data reduction effort.Both the MODIS and CERES sensors are on the EOS-Aqua and EOS-Terra platforms.CERES has employed two different ice particle models to date, with a new two-habit model in development.CERES Edition 2 and Edition 3 products used smooth hexagonal column particle models based on eight discrete size bins and ten particle size distributions from measurements by airborne field campaigns [122][123][124][125]. Edition 4 keeps the particle shape and particle size distributions the same, while applying surface roughness to the particle model.The planned Edition 5 may use a two-habit model [126].Figure 11 schematically depicts the CERES Edition 2 through five particle models.The two-habit model consists of a roughened hexagonal column and the ensemble average of 20 random aggregates of 20 distorted hexagonal columns.The proportion of the single column is defined as follows to satisfy the in-situ airborne measurements: and the proportion of the aggregate column is obtained by f aggregate = 1 − f column .In MODIS Collection 4 and 5 models, the ice particles are assumed to have smooth surfaces, with the exception of the aggregate particle that is assumed to be moderately roughened.Because smooth particles dominate in these ice models, the corresponding phase functions have angularly dependent features that are evident in the phase function comparison between MODIS Collections 5 and 6 in Figure 12.The phase function of the MODIS Collection 6 is featureless because particle surface roughness smooths out the scattering peaks and maxima.Additionally, the asymmetry parameter increases with particle size for solar wavelengths in both the Collection 4 and 5 bulk models, while the asymmetry parameter is essentially independent of particle size for the Collection 6 model (Figure 13).
For two clouds with optical properties (specifically, optical thickness and the asymmetry factor) in a non-absorbing band (t , g) and (t ', g '), a quasi-invariant relation known as the similarity relation [127,128] holds: With the CERES Edition 2 model and the two-habit model as an example for comparison, it is straightforward to speculate that the optical thickness values retrieved with the two-habit model are smaller than their counterparts based on the CERES Edition 2 model.This is because Edition2 (1-g Edition2 ) / (1-g THM ) from Equation ( 31) and (1-g Edition2 ) / (1-g THM ) < 1, as evident in the upper middle panel of Figure 13.In MODIS Collection 4 and 5 models, the ice particles are assumed to have smooth surfaces, with the exception of the aggregate particle that is assumed to be moderately roughened.Because smooth particles dominate in these ice models, the corresponding phase functions have angularly dependent features that are evident in the phase function comparison between MODIS Collections 5 and 6 in Figure 12.The phase function of the MODIS Collection 6 is featureless because particle surface roughness smooths out the scattering peaks and maxima.Additionally, the asymmetry parameter increases with particle size for solar wavelengths in both the Collection 4 and 5 bulk models, while the asymmetry parameter is essentially independent of particle size for the Collection 6 model (Figure 13).
For two clouds with optical properties (specifically, optical thickness and the asymmetry factor) in a non-absorbing band (τ, g) and (τ , g ), a quasi-invariant relation known as the similarity relation [127,128] holds: With the CERES Edition 2 model and the two-habit model as an example for comparison, it is straightforward to speculate that the optical thickness values retrieved with the two-habit model are smaller than their counterparts based on the CERES Edition 2 model.This is because τ THM ≈ τ Edition2 (1 − g Edition2 )/(1 − g THM ) from Equation ( 31) and (1 − g Edition2 )/(1 − g THM ) < 1, as evident in the upper middle panel of Figure 13.Recent studies based on measurements from sensors in the NASA A-Train [129] indicate that, to obtain ice cloud properties that are more independent of the sensor, it is important to include a mechanism that acts to randomize the scattering within an ice particle, so that the phase function has reduced halo features.This finding has been the result of various sensor comparisons.Comparisons  Recent studies based on measurements from sensors in the NASA A-Train [129] indicate that, to obtain ice cloud properties that are more independent of the sensor, it is important to include a mechanism that acts to randomize the scattering within an ice particle, so that the phase function has reduced halo features.This finding has been the result of various sensor comparisons.Comparisons Recent studies based on measurements from sensors in the NASA A-Train [129] indicate that, to obtain ice cloud properties that are more independent of the sensor, it is important to include a mechanism that acts to randomize the scattering within an ice particle, so that the phase function has reduced halo features.This finding has been the result of various sensor comparisons.Comparisons between results using the MODIS ice models and CALIOP measurements have shown that a roughened particle is preferred (e.g., [111,121]).Comparisons between MODIS Collection 5 and POLDER indicate that certain ice habits better match the polarization measurements and also indicate particles should not have smooth facets [120,130,131].Given the variety and scope of current sensors, our goal is to work towards construction of an ice model that minimizes differences in ice cloud properties based on observations made by various sensors.

Spectral Consistency of Bulk Ice Particle Models
Polarization measurements by the POLDER sensor illustrate an aspect of the intrinsic consistency of ice cloud models for solar channels.The modified polarized reflectance product (or equivalently, modified polarized normalized reflectance) is not sensitive to cloud optical thickness when clouds are optically thick (τ > 4) but is sensitive to the selection of ice cloud models.In this study, we analyzed the modified polarized reflectance product in comparison with model simulations based on the CERES Edition 2 model and the two-habit model.Here, the modified polarized reflectance product is defined as [132] where E 0 µ 0 is the flux density crossing through a horizontal plane, E 0 is the beam irradiance, and η is the sign (±1) defined by C.-Labonnote et al. [132].
Figure 14 shows the distribution of cloud pixels of POLDER observations that are used for the following analyses shown in Figure 15.We selected pixels in September 2006 that meet the following criteria: (1) ice phase clouds; (2) MODIS Band 31 (11 µm) brightness temperatures < 233 K; (3) the cloud optical thickness > 4; and (4) POLDER measurements over ocean.We sampled these pixels on a global scale to test the two-habit model because the model is intended to be used globally for retrieval of ice cloud properties.The first two criteria ensure that ice particles are present at cloud top.Optically thin cirrus clouds and ice cloud over land were removed from the analysis to reduce uncertainty associated with surface reflectance.Many cloud pixels are located in the midlatitude storm tracks and in the tropics, especially over Maritime continent and tropical Indian Ocean.Because there is more land cover in the northern hemisphere, the Southern Ocean contributes significantly to the midlatitude data collection.As September is in the tropical cyclone season for the northern hemisphere, the temporal sampling window is primarily responsible for the contrast of the data population in the northern and southern subtropics.In total, 283,594 pixels met our criteria and were used for the analysis.
Atmosphere 2018, 9, x FOR PEER REVIEW 19 of 32 between results using the MODIS ice models and CALIOP measurements have shown that a roughened particle is preferred (e.g., [111,121]).Comparisons between MODIS Collection 5 and POLDER indicate that certain ice habits better match the polarization measurements and also indicate particles should not have smooth facets [120,130,131].Given the variety and scope of current sensors, our goal is to work towards construction of an ice model that minimizes differences in ice cloud properties based on observations made by various sensors.

Spectral Consistency of Bulk Ice Particle Models
Polarization measurements by the POLDER sensor illustrate an aspect of the intrinsic consistency of ice cloud models for solar channels.The modified polarized reflectance product (or equivalently, modified polarized normalized reflectance) is not sensitive to cloud optical thickness when clouds are optically thick (τ > 4) but is sensitive to the selection of ice cloud models.In this study, we analyzed the modified polarized reflectance product in comparison with model simulations based on the CERES Edition 2 model and the two-habit model.Here, the modified polarized reflectance product is defined as [132] ( , , , ) = ( + ), where E 0 m 0 is the flux density crossing through a horizontal plane, E 0 is the beam irradiance, and η is the sign (±1) defined by C.-Labonnote et al. [132].
Figure 14 shows the distribution of cloud pixels of POLDER observations that are used for the following analyses shown in Figure 15.We selected pixels in September 2006 that meet the following criteria: (1) ice phase clouds; (2) MODIS Band 31 (11 μm) brightness temperatures < 233 K; (3) the cloud optical thickness > 4; and (4) POLDER measurements over ocean.We sampled these pixels on a global scale to test the two-habit model because the model is intended to be used globally for retrieval of ice cloud properties.The first two criteria ensure that ice particles are present at cloud top.Optically thin cirrus clouds and ice cloud over land were removed from the analysis to reduce uncertainty associated with surface reflectance.Many cloud pixels are located in the midlatitude storm tracks and in the tropics, especially over Maritime continent and tropical Indian Ocean.Because there is more land cover in the northern hemisphere, the Southern Ocean contributes significantly to the midlatitude data collection.As September is in the tropical cyclone season for the northern hemisphere, the temporal sampling window is primarily responsible for the contrast of the data population in the northern and southern subtropics.In total, 283,594 pixels met our criteria and were used for the analysis.Figure 17 shows the comparison of the retrieved optical thickness from SW and LW techniques.In the upper left panel of Figure 17, it can be seen that, if ice particles are assumed to be spheres, the optical thicknesses retrieved with the IR split window technique are much smaller than their counterparts based on the solar bands using the Nakajima-King algorithm.Thus, spectral consistency in the retrieval is not achieved for the spherical model.This is also true for the CERES Edition 2 model (upper right panel of Figure 17).However, it can be seen in the two lower panels of Figure 17 that both the MODIS C6 model and THM lead to spectrally consistent retrievals for optical thickness up to 3. For larger optical thicknesses, the IR signals at TOA lose sensitivity to the lower portion of the cloud and thus the IR band retrievals are systematically lower than the solar band counterparts.The solar and IR based methods are complementary in this respect, with the IR-based method having greater sensitivity to lower optical thicknesses than the solar method.Based on the comparisons shown in Figure 17, it can be said that both the MODIS C6 model and THM are optimal from the perspective of spectral consistency in retrieving ice cloud properties.Figure 17 shows the comparison of the retrieved optical thickness from SW and LW techniques.In the upper left panel of Figure 17, it can be seen that, if ice particles are assumed to be spheres, the optical thicknesses retrieved with the IR split window technique are much smaller than their counterparts based on the solar bands using the Nakajima-King algorithm.Thus, spectral consistency in the retrieval is not achieved for the spherical model.This is also true for the CERES Edition 2 model (upper right panel of Figure 17).However, it can be seen in the two lower panels of Figure 17 that both the MODIS C6 model and THM lead to spectrally consistent retrievals for optical thickness up to 3. For larger optical thicknesses, the IR signals at TOA lose sensitivity to the lower portion of the cloud and thus the IR band retrievals are systematically lower than the solar band counterparts.The solar and IR based methods are complementary in this respect, with the IR-based method having greater sensitivity to lower optical thicknesses than the solar method.Based on the comparisons shown in Figure 17, it can be said that both the MODIS C6 model and THM are optimal from the perspective of spectral consistency in retrieving ice cloud properties.Figures 18 and 19 show the optical thickness retrieved using different ice particle models by the VIS-SWIR bi-spectral technique.The retrievals were performed from Terra MODIS data acquired over ocean between 55° S and 50° N over a week from 12 to 18 September 2013.Figure 18 indicates that optical thickness values are consistent when retrieved with the MODIS Collection 6 model and two-habit model.The scatter plots for a spherical ice particle model (Figure 18a) and CERES Edition 2 (Figure 18d) have significantly broader distributions than (Figure 18b) CERES Edition 4 and (Figure 18c) MODIS Collection 6 models.
Figure 19 shows the ratios of optical thicknesses between four different ice models and the two-habit model to show consistency/inconsistency between models.In this figure, there are approximately 30,000 data points for every 10° scattering angle bins centered at about 90° to 170°.As indicated in these results, the best consistency is shown between the CERES Edition 4 model to two-habit model in Figure 19b, and between the MODIS Collection 6 model to two-habit model in Figure 19c.The comparisons involving the CERES Edition 2 model and the ice spheres show the least consistency.If the goal is to be able to compare products from different teams, this sort of information is critical for understanding the strengths and weaknesses of the products, and for modelers to make sense of them.
In addition to the ice cloud models discussed above, an inhomogeneous ice particle model (Figure 20a) was developed for ice cloud property retrieval from POLDER observations [132] and an ice cloud model based on the Voronoi shape (Figure 20b) was developed for retrievals from the Himawari-8 AHI observations [133,134].Moreover, van Diedenhoven et al. [135,136] assumed hexagonal particles with incorporation of surface roughness for retrievals based on airborne scanning polarimetric measurements.A thorough comparison is made between the POLDER and MODIS Collection 5 ice cloud optical property models in [137].A case study in [137] shows Figures 18 and 19 show the optical thickness retrieved using different ice particle models by the VIS-SWIR bi-spectral technique.The retrievals were performed from Terra MODIS data acquired over ocean between 55 • S and 50 • N over a week from 12 to 18 September 2013.Figure 18 indicates that optical thickness values are consistent when retrieved with the MODIS Collection 6 model and two-habit model.The scatter plots for a spherical ice particle model (Figure 18a) and CERES Edition 2 (Figure 18d) have significantly broader distributions than (Figure 18b) CERES Edition 4 and (Figure 18c) MODIS Collection 6 models.
Figure 19 shows the ratios of optical thicknesses between four different ice models and the two-habit model to show consistency/inconsistency between models.In this figure, there are approximately 30,000 data points for every 10 • scattering angle bins centered at about 90 • to 170 • .As indicated in these results, the best consistency is shown between the CERES Edition 4 model to two-habit model in Figure 19b, and between the MODIS Collection 6 model to two-habit model in Figure 19c.The comparisons involving the CERES Edition 2 model and the ice spheres show the least consistency.If the goal is to be able to compare products from different teams, this sort of information is critical for understanding the strengths and weaknesses of the products, and for modelers to make sense of them.
, illustrating the sensitivity of the retrieved ice cloud optical thickness to the forward optical property model used in implementing the retrieval algorithm.In addition to the ice cloud models discussed above, an inhomogeneous ice particle model (Figure 20a) was developed for ice cloud property retrieval from POLDER observations [132] and an ice cloud model based on the Voronoi shape (Figure 20b) was developed for retrievals from the Himawari-8 AHI observations [133,134].Moreover, van Diedenhoven et al. [135,136] assumed hexagonal particles with incorporation of surface roughness for retrievals based on airborne scanning polarimetric measurements.A thorough comparison is made between the POLDER and MODIS Collection 5 ice cloud optical property models in [137].A case study in [137] shows τ POLDER ≈ 0.7τ MODIS , illustrating the sensitivity of the retrieved ice cloud optical thickness to the forward optical property model used in implementing the retrieval algorithm.

Summary
In this review, we survey recent studies related to the construction of ice cloud bulk scattering models and some of the choices that teams made for their derivation.We review the history of the development of ice cloud models, beginning with an early model assuming ice particles to be horizontally oriented circular cylinders, and moving to the models being used in current satellite retrieval products.We discuss whether horizontally-oriented particles exist in ice clouds, and to what degree.This analysis is based on CALIOP data analysis.This is important because ice particle orientation increases the complexity of ice particle single-scattering properties.We illustrate that it is valid to assume that ice particles are randomly oriented in the atmosphere for passive remote sensing purposes without considering particle orientation.
While a team may construct these ice models for a single sensor, there is increasing work being done to enhance the consistency of cloud property retrievals across a variety of spaceborne sensors.That is, there is more work being done to merge data from multiple sensors to improve the ice cloud

Summary
In this review, we survey recent studies related to the construction of ice cloud bulk scattering models and some of the choices that teams made for their derivation.We review the history of the development of ice cloud models, beginning with an early model assuming ice particles to be horizontally oriented circular cylinders, and moving to the models being used in current satellite retrieval products.We discuss whether horizontally-oriented particles exist in ice clouds, and to what degree.This analysis is based on CALIOP data analysis.This is important because ice particle orientation increases the complexity of ice particle single-scattering properties.We illustrate that it is valid to assume that ice particles are randomly oriented in the atmosphere for passive remote sensing purposes without considering particle orientation.
While a team may construct these ice models for a single sensor, there is increasing work being done to enhance the consistency of cloud property retrievals across a variety of spaceborne sensors.That is, there is more work being done to merge data from multiple sensors to improve the ice cloud retrievals, and subsequently to compare the satellite-based ice cloud properties with those from models.One would hope that the ice cloud properties derived from two sensors on the same platform, e.g., VIIRS (an imager) and CrIS (a hyperspectral IR sounder) would result in consistent global products.As noted in this review, there is a tendency for each team to develop their own ice cloud bulk single-scattering model in isolation.For example, there are teams that use a gamma particle size distribution and a single ice habit, and other teams that make more extended use of in-situ microphysical measurements of PSDs and different habits.Since various retrieval teams may develop an ice bulk scattering model that fits their needs, there must be studies that compare the results between sensor teams.Our goal is to reach better consistency in retrievals from imagers, hyperspectral infrared sounders, polarimeters, and active lidar measurements.This will extend in the future to encompass millimeter and sub-millimeter measurements and perhaps polarimetric measurements over a broader spectral range.
As examples shown in this study, comparisons of optical thickness are made between a visible-shortwave infrared technique and an infrared channel technique.We demonstrate that the MODIS Collection 6 model and the CERES two-habit model give spectrally consistent retrievals of ice cloud microphysical and optical properties.Additionally, the same bulk scattering models are used to compare results from a space-based polarimeter (POLDER).These two models incorporate severely roughened particle surfaces that lead to featureless bulk phase functions.At present, the MODIS Collection 6 and the CERES two-habit models are optimal ice particle models for passive remote sensing of global ice clouds.Polarimetric-based observations of ice clouds tend to focus on ice particles near the top of the cloud, while shortwave-infrared channels may penetrate the ice cloud to different depths.Infrared-based ice cloud retrievals have a more limited optical thickness range due to the increased particle absorption.Active sensors provide a wealth of information about the columnar ice cloud optical/radiative and microphysical properties but the measurements are limited by their narrow beams.With the effort to assemble multiple datasets, there is then an opportunity to gain a wider perspective of the range of ice cloud optical/radiative and microphysical properties.The eventual goal is to find a framework for inferring these properties, which leads to more consistent results across the sensors where new measurements can add to the body of knowledge of ice clouds.

Figure 1 .
Figure 1.The decomposition of the electric field vector associated with a plane electromagnetic wave into parallel and perpendicular components ( E // and E ^) with respect to a reference plane.

Figure 1 .
Figure 1.The decomposition of the electric field vector associated with a plane electromagnetic wave into parallel and perpendicular components (E // and E ⊥ ) with respect to a reference plane.

Figure 2 .
Figure 2. A schematic diagram illustrating the incident and scattering configuration for the scattering of light by a non-spherical particle.

Figure 2 .
Figure 2. A schematic diagram illustrating the incident and scattering configuration for the scattering of light by a non-spherical particle.

Figure 3 .
Figure 3. Schematics showing horizontally oriented cylinder models following Liou [92]: (a) uniform orientation in the horizontal plane; and (b) random orientation in the horizontal plane.Note that particles in (b) are rotated about dotted lines from particles in (a).

Figure 3 .
Figure 3. Schematics showing horizontally oriented cylinder models following Liou [92]: (a) uniform orientation in the horizontal plane; and (b) random orientation in the horizontal plane.Note that particles in (b) are rotated about dotted lines from particles in (a).
, upper panel) and hexagonal columns (Figure 4, lower panel) assuming a horizontally random orientation condition.The plates are randomly oriented in terms of their rotations around the asymmetry axes that are aligned vertically and columns are horizontally oriented with two parallel faces perpendicular to the vertical direction (known as Parry oriented columns).The physical-geometric optical computational program developed by Sun et al. [89] is used for this simulation.

Figure 4 .
Figure 4.The two-dimensional phase function (P11) of: (upper panel) horizontally oriented ice plates (i.e., the asymmetry axes of hexagonal plates are aligned vertically, but the particles are oriented in a rotationally random manner); and (lower panel) Parry oriented columns (i.e., the symmetry axes of hexagonal columns are aligned horizontally with two parallel prism side facets perpendicular to the vertical direction).The particle aspect ratio is a/L = 1000 μm/250 μm for a plate and a/L = 250 μm/1000 μm for a column (a and L are the semi-width and thickness, respectively) at wavelength =0.532 μm.The incident configuration is illustrated in the left diagram.

Figure 5 .
Figure 5. Simulated optical phenomena associated with light scattering by: (a) a mixture of horizontally oriented hexagonal ice plates (7.9%) and randomly oriented counterparts (92.1%;) and (b) the Parry oriented columns.For the simulated sky image, multiple scattering is not included.See also Liou and Yang [4].

Figure 4 .
Figure 4.The two-dimensional phase function (P 11 ) of: (upper panel) horizontally oriented ice plates (i.e., the asymmetry axes of hexagonal plates are aligned vertically, but the particles are oriented in a rotationally random manner); and (lower panel) Parry oriented columns (i.e., the symmetry axes of hexagonal columns are aligned horizontally with two parallel prism side facets perpendicular to the vertical direction).The particle aspect ratio is a/L = 1000 µm/250 µm for a plate and a/L = 250 µm/1000 µm for a column (a and L are the semi-width and thickness, respectively) at wavelength λ = 0.532 µm.The incident configuration is illustrated in the left diagram.

Figure 4 .
Figure 4.The two-dimensional phase function (P11) of: (upper panel) horizontally oriented ice plates (i.e., the asymmetry axes of hexagonal plates are aligned vertically, but the particles are oriented in a rotationally random manner); and (lower panel) Parry oriented columns (i.e., the symmetry axes of hexagonal columns are aligned horizontally with two parallel prism side facets perpendicular to the vertical direction).The particle aspect ratio is a/L = 1000 μm/250 μm for a plate and a/L = 250 μm/1000 μm for a column (a and L are the semi-width and thickness, respectively) at wavelength =0.532 μm.The incident configuration is illustrated in the left diagram.

Figure 5 .
Figure 5. Simulated optical phenomena associated with light scattering by: (a) a mixture of horizontally oriented hexagonal ice plates (7.9%) and randomly oriented counterparts (92.1%;) and (b) the Parry oriented columns.For the simulated sky image, multiple scattering is not included.See also Liou and Yang [4].

Figure 5 .
Figure 5. Simulated optical phenomena associated with light scattering by: (a) a mixture of horizontally oriented hexagonal ice plates (7.9%) and randomly oriented counterparts (92.1%;) and (b) the Parry oriented columns.For the simulated sky image, multiple scattering is not included.See also Liou and Yang [4].

Figure 6 .
Figure 6.(a) Schematic diagram illustrating the relation between the layer-integrated backscatter (γ) and the layer-integrated depolarization (δ) associated with ice clouds; (b) ice cloud γ-δ relation observations made by CALIOP during April 2007 when the CALIOP lidar was in nadir-viewing configuration; and (c) ice cloud γ-δ relation observations made by CALIOP during April 2009 when the CALIOP lidar pointed 3° off nadir.

Figure 6 .
Figure 6.(a) Schematic diagram illustrating the relation between the layer-integrated backscatter (γ) and the layer-integrated depolarization (δ) associated with ice clouds; (b) ice cloud γ-δ relation observations made by CALIOP during April 2007 when the CALIOP lidar was in nadir-viewing configuration; and (c) ice cloud γ-δ relation observations made by CALIOP during April 2009 when the CALIOP lidar pointed 3 • off nadir.

Figure 7 .
Figure 7.The fraction of oriented ice particles derived from one month (April 2007) of observations made by CALIPSO CALIOP and IIR based on the method reported by Saito et al. [100].

Figure 7 .
Figure 7.The fraction of oriented ice particles derived from one month (April 2007) of observations made by CALIPSO CALIOP and IIR based on the method reported by Saito et al. [100].

Figure 8 .
Figure 8. Schematic diagram illustrating the principle of the Nakajima-King bi-spectral method [42] for simultaneous retrieval of cloud optical thickness and effective radius.The look-up table shown in the diagram is generated with a solar zenith of 60 • , a viewing zenith angle of 45 • , and a relative azimuthal angle of 105 • for a severely roughened ice particle model.Courtesy of Yi Wang (Texas A&M University).

Figure 9 .Figure 10 .
Figure 9. Infrared retrieval lookup table generated by using radiances from MODIS bands 31 and 32 and applying the two-habit ice cloud model (see Figure 11).Black lines indicate different visible optical thickness from 0.01 to 10. Red, blue green, purple, and yellow lines indicate different effective diameter from 6 to 156 μm.Arrowed dashed lines show retrieval iteration processes.

Figure 9 .
Figure 9. Infrared retrieval lookup table generated by using radiances from MODIS bands 31 and 32 and applying the two-habit ice cloud model (see Figure 11).Black lines indicate different visible optical thickness from 0.01 to 10. Red, blue green, purple, and yellow lines indicate different effective diameter from 6 to 156 µm.Arrowed dashed lines show retrieval iteration processes.

Figure 9 .Figure 10 .
Figure 9. Infrared retrieval lookup table generated by using radiances from MODIS bands 31 and 32 and applying the two-habit ice cloud model (see Figure 11).Black lines indicate different visible optical thickness from 0.01 to 10. Red, blue green, purple, and yellow lines indicate different effective diameter from 6 to 156 μm.Arrowed dashed lines show retrieval iteration processes.

Figure 11 .
Figure 11.(a) Ice cloud models used for CERES Editions 2-4; and (b) a two-habit model as a candidate ice model for future CERES Edition 5.

Figure 11 .
Figure 11.(a) Ice cloud models used for CERES Editions 2-4; and (b) a two-habit model as a candidate ice model for future CERES Edition 5.

Figure 12 .
Figure 12.Comparison of the phase functions at wavelength 0.86 μm based on: (a) MODIS Collection 4, 5, and 6; (b) CERES Edition 2 and 4, and the two-habit model; and (c) MODIS Collection 6 and the two-habit model.(d-f) Counterparts of (a-c) at wavelength 2.13 μm.Effective radius is fixed at 30 μm.

Figure 13 .
Figure 13.Comparison of the asymmetry parameter at wavelength 0.86 μm based on: (a) MODIS Collection 4, 5, and 6; (b) CERES Edition 2 and 4, and the two-habit model; and (c) MODIS Collection 6 and the two-habit model.(d-f) Counterpart of (a-c) at wavelength 2.13 μm.

Figure 12 . 32 Figure 12 .
Figure 12.Comparison of the phase functions at wavelength 0.86 µm based on: (a) MODIS Collection 4, 5, and 6; (b) CERES Edition 2 and 4, and the two-habit model; and (c) MODIS Collection 6 and the two-habit model.(d-f) Counterparts of (a-c) at wavelength 2.13 µm.Effective radius is fixed at 30 µm.

Figure 13 .
Figure 13.Comparison of the asymmetry parameter at wavelength 0.86 μm based on: (a) MODIS Collection 4, 5, and 6; (b) CERES Edition 2 and 4, and the two-habit model; and (c) MODIS Collection 6 and the two-habit model.(d-f) Counterpart of (a-c) at wavelength 2.13 μm.

Figure 13 .
Figure 13.Comparison of the asymmetry parameter at wavelength 0.86 µm based on: (a) MODIS Collection 4, 5, and 6; (b) CERES Edition 2 and 4, and the two-habit model; and (c) MODIS Collection 6 and the two-habit model.(d-f) Counterpart of (a-c) at wavelength 2.13 µm.

Figure 14 .
Figure 14.Distribution of cloud pixels used to produce the observation density distribution in Figure 15.

Figure 14 .
Figure 14.Distribution of cloud pixels used to produce the observation density distribution in Figure 15.

Figure 15 32 Figure 15
Figure 15 shows the probability density distribution of the modified polarized reflectance from the POLDER observations and the counterparts from the simulations assuming the two-habit model and the CERES Edition 2 model.Overall, the observed modified polarized reflectance decreases with increasing the scattering angle over a range of 100-170 • .The simulations based on the two-habit model closely match with those from the observations except at the scattering angles with 100-120 • because of the weak polarized side scattering by the Two-habit model, whereas the simulations based on the CERES Edition 2 substantially deviate from the observations.Therefore, the two-habit model is much more appropriate than the CERES Edition 2 model in terms of the consistency between the model-based modified polarized reflectance and the measured counterparts at the scattering angles within 100-170 • .Atmosphere 2018, 9, x FOR PEER REVIEW 20 of 32

Figure 15 .
Figure 15.Comparison between observed modified polarized reflectance product values and the simulations based on the two-habit and CERES Edition 2 ice cloud models.

Figure 15 .
Figure 15.Comparison between observed modified polarized reflectance product values and the simulations based on the two-habit and CERES Edition 2 ice cloud models.Furthermore, to demonstrate the dependence of retrieved cloud properties on particle habit, we present results from shortwave (SW) and longwave (LW) retrievals.The data were acquired by Aqua MODIS over the ocean at latitudes between 55 • S and 50 • N for two orbits on 10 September 2006.The visible composite image and the MODIS cloud phase are presented in Figure 16.The analysis was limited to pixels that are identified as ice by the MODIS Collection 6 cloud flag and with optical thickness less than 4. The results indicate that SW and LW retrievals are consistent with use of the two-habit or MODIS Collection 6 ice particle model.Figure17shows the comparison of the retrieved optical thickness from SW and LW techniques.In the upper left panel of Figure17, it can be seen that, if ice particles are assumed to be spheres, the optical thicknesses retrieved with the IR split window technique are much smaller than their counterparts based on the solar bands using the Nakajima-King algorithm.Thus, spectral consistency in the retrieval is not achieved for the spherical model.This is also true for the CERES Edition 2 model (upper right panel of Figure17).However, it can be seen in the two lower panels of Figure17that both the MODIS C6 model and THM lead to spectrally consistent retrievals for optical thickness up to 3. For larger optical thicknesses, the IR signals at TOA lose sensitivity to the lower portion of the cloud and thus the IR band retrievals are systematically lower than the solar band counterparts.The solar and IR based methods are complementary in this respect, with the IR-based method having greater sensitivity to lower optical thicknesses than the solar method.Based on the comparisons shown in Figure17, it can be said that both the MODIS C6 model and THM are optimal from the perspective of spectral consistency in retrieving ice cloud properties.

Figure 16 .
Figure 16.(a) Visible composite image of two orbits used in the SW and LW comparison; and (b) reflectivity image colored by the MODIS cloud phase.The orange color indicates the clouds are identified as ice clouds.

Figure 16 .
Figure 16.(a) Visible composite image of two orbits used in the SW and LW comparison; and (b) reflectivity image colored by the MODIS cloud phase.The orange color indicates the clouds are identified as ice clouds.

Figure 17 .
Figure 17.Comparison of retrieved optical thickness values from the shortwave technique (shortwave bi-spectral method) and the longwave technique (split-window method): (a) ice sphere; (b) CERES Edition 4 model; (c) MODIS Collection 6 model; and (d) two-habit model (CERES Edition 5 model).

Figure 17 .
Figure 17.Comparison of retrieved optical thickness values from the shortwave technique (shortwave bi-spectral method) and the longwave technique (split-window method): (a) ice sphere; (b) CERES Edition 4 model; (c) MODIS Collection 6 model; and (d) two-habit model (CERES Edition 5 model).

Figure 18 .
Figure 18.Comparison of retrieved optical thickness values from the shortwave technique (shortwave bi-spectral method): (a) ice sphere and two-habit model (THM; potential CERES Edition 5 model); (b) CERES Edition 4 model and two-habit model; (c) MODIS Collection 6 model and two-habit model; and (d) CERES Edition 2 model and two-habit model.

Figure 18 . 32 Figure 19 .
Figure 18.Comparison of retrieved optical thickness values from the shortwave technique (shortwave bi-spectral method): (a) ice sphere and two-habit model (THM; potential CERES Edition 5 model); (b) CERES Edition 4 model and two-habit model; (c) MODIS Collection 6 model and two-habit model; and (d) CERES Edition 2 model and two-habit model.Atmosphere 2018, 9, x FOR PEER REVIEW 24 of 32

Figure 19 .
Figure 19.Ratio of retrieved optical thickness values to that of the two-habit model with retrieval results from the shortwave bi-spectral method: (a) ice sphere to two-habit model (THM; potential CERES Edition 5 model); (b) CERES Edition 4 model to two-habit model; (c) MODIS Collection 6 model to two-habit model; and (d) CERES Edition 2 model to two-habit model.

Figure 19 .
Figure 19.Ratio of retrieved optical thickness values to that of the two-habit model with retrieval results from the shortwave bi-spectral method: (a) ice sphere to two-habit model (THM; potential CERES Edition 5 model); (b) CERES Edition 4 model to two-habit model; (c) MODIS Collection 6 model to two-habit model; and (d) CERES Edition 2 model to two-habit model.