On the Radiative Properties of Ice Clouds: Light Scattering, Remote Sensing, and Radiation Parameterization

Ice clouds, consisting of various kinds of nonspherical ice crystals and frequently occurring in the upper troposphere and into the lower stratosphere, play an important role in various physical processes in the Earth's atmosphere (Liou, 1986; Stephens et al., 1990; Liou et al., 2000, 2001; Lynch et al., 2002; Baran, 2012). For example, by reflecting solar radiation and absorbing thermal infrared emission from the lower atmosphere and the surface, ice clouds regulate the radiation budget of the Earth-atmosphere system, and thus affect the atmospheric cooling/heating profile. In comparison with water clouds, which consist of spherical liquid water droplets (Slingo, 1989), the current level of understanding of the microphysical and optical properties of ice clouds is relatively low; in particular, substantial uncertainties remain in the representation of ice clouds in radiative transfer models used for a variety of applications. To improve our understanding of the bulk optical properties of ice clouds, we review three lines of ice cloud studies including: (a) the simulation of the single-scattering properties of individual ice crystals; (b) the use of the ice cloud bulk scattering properties to infer ice cloud parameters from observations made by an array of passive and active satellite sensors that measure different parts of the electromagnetic spectrum, particularly the solar and thermal infrared bands; and (c) the ice cloud radiation parameterizations for use in weather and climate models on the basis of microphysical data obtained from numerous field campaigns.

With respect to the representation of ice clouds in broadband radiative transfer models, the transition from the use of ice spheres (e.g., Houghton and Hunt, 1971; Ramaswamy and Detwiler, 1986) and infinitely-long circular cylinders (Liou, 1972a, 1972b; Stephens, 1980a, 1980b) to more realistic hexagonal habits, which occur in the Earth's atmosphere, began after the seminal study presented by Takano and Liou (1989a, 1989b). The availability of the scattering and absorption properties of hexagonal columns and plates led to attempts to incorporate these properties into broadband radiation parameterization models (e.g., Ebert and Curry, 1992, 1993; Fu and Liou, 1993; Fu, 1996, 2007; Mitchell et al., 2006; Baran, 2012). Moreover, the phase function for hexagonal ice particles has been extensively utilized to support satellite remote sensing of cirrus clouds (e.g., Minnis et al., 1993a, 1993b; Sun and Shine, 1995). The optical properties of complex ice habits (Macke, 1993; Macke et al., 1996, 1998; Takano and Liou, 1995; Iaquinta et al., 1995; Yang and Liou, 1998) have also been reported for remote sensing application and radiation parameterization (Key et al., 2002; McFarquhar et al., 2002; Platnick et al., 2003; King et al., 2004; Hong et al., 2009a; Zhang et al., 2015).

Further insight into the microphysical and optical properties of ice particles has been obtained based on data collected during a number of field campaigns focused on cirrus clouds (e.g., Heymsfield et al., 1990; Francis, 1995; Heymsfield et al., 2002; Gayet et al., 2004, 2006). Initially, field campaigns were designed to facilitate ice cloud in situ measurements to support satellite retrievals of the spatial distribution and optical depth of cirrus clouds from the International Satellite Cloud Climatology Project (ISCCP) (Starr and Wylie, 1990; Rossow and Schiffer, 1999), but subsequent field campaigns have been designed to study ice clouds worldwide (Heymsfield et al., 2013). A number of issues have been discovered with regard to aircraft-based in situ measurements. For example, a high percentage of the small ice particle concentration measured from aircraft could be artificially produced by the shattering of large ice particles that occur at the inlets of the 2D-C probes (e.g., Field et al., 2006; Lawson, 2011; Korolev et al., 2011, 2013a, 2013b). (Jackson and McFarquhar, 2014) investigated the impact of potential ice particle shattering on the bulk scattering properties inferred from particle size distributions and showed that shattering increases the number concentration, but has much less impact on the properties of higher order moments such as the asymmetry parameter or the single-scattering albedo. Other important issues for the single-scattering properties of ice particles are related to their intricate and complex morphologies such as aggregation, fractures, embedded aerosol particles (e.g., soot particles) or air bubbles, as well as surface roughness. These imperfections alter the scattering and absorption properties of individual ice particles as well as their bulk radiative properties (Macke et al., 1996; Yang and Liou, 1998; C.-Labonnote et al., 2001; Shcherbakov et al., 2006; Ulanowski et al., 2006; Nousiainen and Muinonen, 2007; Yang et al., 2008; Baum et al., 2010; Neshyba et al., 2013; Yang et al., 2013; Ulanowski et al., 2014).

The strengths and weaknesses of global satellite-based ice cloud property retrievals need to be well characterized, including the uncertainties resulting from radiometric noise, use of ancillary data, algorithms, and choices made to filter and aggregate the data into global datasets. The resulting global ice cloud properties in turn provide crucial constraints on the parameterization of various cloud microphysical processes and cloud-radiation interactions. To this end, another line of investigation has been conducted to improve the inference of ice cloud properties from satellite sensors, including both polar-orbiting and geostationary imagers. To work with satellite sensors in an operational environment, radiative transfer (RT) models have been used to generate look-up tables (LUTs) of reflectivity and transmissivity (and other properties) as functions of optical depth, effective particle size, and viewing geometry (solar zenith, viewing zenith, and relative azimuth angles). The datasets of the optical properties of individual ice crystals and models of spectral or narrowband bulk single-scattering properties developed for various remote sensing applications have been reported for solar, infrared, far-infrared, and microwave portions of the electromagnetic spectrum (Hess and Wiegner, 1994; Hess et al., 1998; Yang et al., 2000; Yang et al., 2003; Yang et al., 2005; Baum et al., 2005a, 2005b; Kim, 2006; Baum et al., 2007; Liu, 2008; Hong et al., 2009b; Petty and Huang, 2010; Baum et al., 2012; Yang et al., 2013; Baum et al., 2014; and references cited therein). To improve the representation of ice clouds in RT models of varying complexity across a wide range of wavelengths, the incorporation of an adequate representation of ice cloud radiative properties is critical. (Baran, 2009) provided a review of issues similar to those discussed in this paper, and articulated the use of polarization measurements to provide further insight into the remote sensing of ice clouds. In this paper, we review the latest advances in the preceding research areas.

This review paper is organized into five sections. Section 2 reviews several widely used light-scattering computational techniques, including the improved geometric optics, finite-difference time domain (FDTD), pseudo-spectral time domain (PSTD), discrete dipole approximation (DDA), and the T-matrix methods, and presents pertinent numerical examples of the single-scattering properties of ice crystals. Section 3 gives a brief overview of the satellite remote sensing of ice cloud properties and discusses several operational retrieval algorithms to infer cloud optical depth, effective ice crystal size and habit, and cloud top pressure. Two case studies are presented to illustrate both the capabilities and shortcomings of the current retrieval algorithms. With advanced sensors being launched or under development, new approaches and algorithms are needed to interpret the available data. Furthermore, by combining data from multiple instruments that encompass a broad range of spectral channels as well as polarization capabilities, a more complete and accurate quantification of ice clouds can be determined. Section 4 reviews the parameterization of ice cloud bulk single-scattering properties for broadband models that can be used in global and regional numerical models, and includes the presentation of a number of relevant results. Section 5 summarizes the review.

Simulation of electromagnetic scattering is an interdisciplinary area between applied physics and computational mathematics. The most famous problem on this subject may be the scattering of light by a sphere, which is formulated as the Lorenz-Mie theory and is explained in detail in a number of texts (e.g., van de Hulst, 1957; Bohren and Huffman, 1983; Liou, 2002; Wendisch and Yang, 2012). Numerically robust computer codes based on the Lorenz-Mie theory have been developed and widely used (e.g., Wiscombe, 1980). The rich history of the study of the optical properties of spheres was thoroughly reviewed by (Logan, 1965). Although the Lorenz-Mie theory facilitates efficient computation of the optical properties of a sphere, many particles (e.g., ice crystals within cirrus clouds) are highly nonspherical. In the past three decades, the electromagnetic scattering research community has devoted significant effort to the study of the single-scattering properties (namely, the extinction efficiency, single-scattering albedo, and complete phase matrix) of nonspherical/inhomogeneous particles, from which various analytical formalisms and numerical techniques have been developed. Interested readers are referred to several excellent monographs and review articles (Mishchenko et al., 2000, 2002; Kahnert, 2003; Doicu et al., 2006; Borghese et al., 2007; Wriedt, 2009). Here, we review several methods that have often been utilized to simulate the scattering of light by ice crystals.

Geometric optics principles can be used to calculate the single-scattering properties of particles much larger than the wavelength of incident light. When an electromagnetic plane wave impinges on a particle, the wave front blocked by the particle's cross section is considered to consist of idealized narrow beams referred to as rays. Each ray propagates along a straight-line path until it impinges on the particle-medium (i.e., the ambient medium surrounding the particle) interface where reflected and refracted rays emerge in new directions. Over the course of successive reflections and refractions, a hierarchy of rays is generated and an angular distribution of the electric field associated with the scattered light is obtained by tracing the electric field vectors associated with the rays. In addition to the angular scattering pattern associated with the ray-tracing procedure, according to Huygens's principle, the blocking of the incident wave front also causes the scattering of the incident light. In other words, the incomplete wave front (the original wave front minus the blocked wave front) contributes to the scattering of the incident light, a phenomenon known as diffraction. Hence, the total scattered light can be obtained by a direct superposition of the scattered fields produced by diffraction and reflected/refracted rays.

The angular distribution of the scattered light associated with diffraction is dependent on a particle's geometric cross section and, once the particle shape is defined, can be computed by means of the Fraunhofer diffraction theory or the Kirchhoff formula (e.g., van de Hulst, 1957; Cai and Liou, 1982; Macke, 1993; Macke et al., 1996; Yang and Liou, 1998; Bi et al., 2011a; Takano et al., 2012). The scattered field generated by ray tracing is dependent on the particle geometry as well as the refractive index. The directions of the reflected and refracted rays can be determined by Snell's law, while the amplitude change and phase of the electric field vector associated with each ray can then be computed in a ray-tracing process according to the Fresnel formulae and the accumulated ray path-length. With modifications, the ray-tracing technique can be applied to various nonspherical particles. As an example, Fig. 1a shows the application of a Monte Carlo approach to initializing random ray-tracing processes for a hexagonal ice crystal. Each point on the particle surface facing the direction of the incident light represents an incoming ray localized on the incident wave front. Figure 1b is a schematic diagram showing the ray-tracing process involving a hexagonal particle. The scattered field consists of the superposition of the fields from all externally reflected rays and refracted rays from the particle to the medium.

Figure 1.  (a) Incident points for initial rays. (b) Schematic diagram illustrating the ray-tracing process involving a hexagonal ice crystal.

Application of the ray-tracing procedure to quantitatively compute the phase function of ice crystals can be traced back to (Jacobowitz, 1971) in whose calculations a hexagonal prism with an infinite length was assumed. Jacobowitz's results quantitatively produced the 22^° halo and illustrated the inadequacy of the Lorenz-Mie theory for simulating the optical properties of large ice crystals. (Wendling et al., 1979) applied the same principle to a finite hexagonal prism and quantitatively depicted the 46^° halo pattern. Both (Jacobowitz, 1971) and (Wendling et al., 1979) focused on the intensity distributions of the scattered light and neglected the potential interference among scattered rays and the polarization state. In a later development, (Cai and Liou, 1982) included the polarization configuration into the ray-tracing process and determined the scattering properties of the full phase matrix elements for randomly orientated columns and plates. (Takano and Liou, 1989a) incorporated ice crystal size distribution and the ice birefringence effect into the ray-tracing process by splitting refracted rays into ordinary and extraordinary components. Moreover, for large absorbing particles, Yang and Liou (1995, 2009a, 2009b) addressed the nature of inhomogeneous waves (Born and Wolf, 1959; Bohren and Huffman, 1983) in the ray-tracing process, and thus derived an effective refractive index and generalized Fresnel's formulae.

The ray-tracing approach has been widely used in modeling light scattering by various nonspherical ice crystals (e.g., Macke, 1993; Macke et al., 1996; Muinonen et al., 1996; Nousiainen and McFarquhar, 2004; Um and McFarquhar, 2007). However, the range of validity for this approach has not been precisely defined, particularly for particles with moderate values of the size parameter (defined as π times the ratio of the particle's circumference dimension to the incident wavelength). As the wavelength increases in relation to the particle size, the size parameter decreases. For this reason, much work has focused on the accuracy of ray tracing for small-to-moderate size parameters. For example, comparison of the ray-tracing solution with an analytic solution based on the exact Lorenz-Mie theory (Liou and Hansen, 1971) revealed that the geometric-optics solution is reasonably accurate for poly-disperse spheres when the size parameter is larger than 100. Similar findings have been reported in the cases of spheroids and cylinders, for which the solutions can be rigorously obtained from the T-matrix method (Macke et al., 1995; Yang et al., 2007; Bi et al., 2013a, 2013b). Note that spheres, spheroids, and the side surfaces of circular cylinders have both smooth and rounded surfaces, but realistic ice crystals are faceted. For ice particles that have flat local surfaces, the cross section of a ray does not diverge after transmitting through or being reflected by two parallel surfaces, leading to caustics in the far-field region. For this reason, pronounced singularities in the phase function may be caused for ice particles with a fixed orientation. Furthermore, even for randomly oriented ice crystals, the ray-tracing approach produces unrealistic delta-transmission peaks in the 0^° forward and 180^° backward directions related to the interactions between rays and the parallel facets of ice crystals (Takano and Liou, 1989a; Mishchenko and Macke, 1998).

To circumvent the aforementioned shortcomings, improvements have been suggested to incorporate fundamental electromagnetic theories into ray propagation from the near-field zone to the far-field zone. In these improvements, ray tracing is only applied to compute the near field within the particle's interior or on the particle's external surface. The single-scattering properties of an ice crystal can then be computed according to exact electromagnetic integral equations. Specifically, the scattered electric field E _sca at the position r when \(\bmr\to \infty\) can be related to the internal electric field E as (Saxon, 1973) \begin{eqnarray} \bm{E}_{\rm sca}(\bm{r})&=&\dfrac{k^2\exp(ik\bm{r})}{4\pi\bm{r}}\iiint_v[\varepsilon(\bm{r}')-1]\\}\bm{E}(\bm{r}')-\quad\nonumber\\ &&\hat{\bm{r}}[\hat{\bm{r}}\cdot \bm{E}(\bm{r}')]\}\exp(-ik\hat{\bm{r}}\cdot \bm{r}')\bm{dr}' , \end{eqnarray} where the integration domain v is the volume occupied by the particle, \(\hat\bmr\) is a unit vector pointing along r, k is the wavenumber and ε is the permittivity, and the extinction and absorption cross sections can be obtained from the following equations (Hage et al., 1991): \begin{eqnarray} \label{eq2} \sigma_{\rm ext}&=&{\rm Im}\left[\dfrac{k}{|\bm{E}_{\rm inc}|^2}\iiint_V(\varepsilon-1)\bm{E}_{\rm inc}^\ast(\bm{r}')\cdot \bm{E}(\bm{r}')d^3\bm{r}'\right] ,\ \ \\ \label{eq3} \sigma_{\rm abs}&=&\dfrac{k}{|\bm{E}_{\rm inc}|^2}\iiint_V{\rm Im}(\varepsilon)\bm{E}^\ast(\bm{r}')\cdot \bm{E}(\bm{r}')d^3\bm{r}' , \end{eqnarray} where E _inc is the incident electric field. An alternative method to compute the scattered field is to integrate the tangent components of the electric field E and magnetic field H on the particle surface as follows (Jackson, 1975): \begin{eqnarray} \label{eq4} \bm{E}_{\rm sca}(\bm{r})&=&\dfrac{k^2\exp(ik\bm{r})}{-ik\bm{r}4\pi}\hat{\bm{r}}\times \iint\\}\hat{\bm{n}}_{\rm s}\times \bm{E}(\bm{r}')-\hat{\bm{r}}\times\nonumber\quad\\ &&[\hat{\bm{n}}_{\rm s}\times \bm{H}(\bm{r}')]\}\exp(-ik\hat{\bm{r}}\cdot \bm{r}')\bm{d}^2\bm{r} , \end{eqnarray} where \(\hat\bm n_\rm s\) is a unit vector that is locally normal to the particle surface. (Muinonen, 1989) pioneered the application of Eq. (5) to ice crystal optics by considering the spreading effect of scattered rays. Yang and Liou (1996b, 1997) systematically developed approaches by mapping the near field computed from geometric optics to the far field through Eqs. (1) to (5). Another conceptual advance employs a broad-beam tracing technique to compute the near field for non-absorbing ice crystals (Popov, 1996; Borovoi and Grishin, 2003) instead of mapping a group of individual rays. As a follow up to Yang and Liou (1997), Bi et al. (2011b) proposed a new computational approach to simulate the scattering and absorptive properties of ice crystals by accounting for the inhomogeneity of electric field amplitude over the beam cross-section for absorptive particles on the basis of Eq. (1).

The method of combining the Fraunhofer diffraction and the angular distribution computed by ray tracing is termed as the conventional geometric optics method (CGOM). Figure 2 shows a comparison of the asymmetry factor, the extinction efficiency, and three phase matrix elements (P₁₁, P₁₂/P₁₁, and P₂₂/P₁₁) computed from the CGOM and the improved geometric-optics method (IGOM) based on an electromagnetic volume-integral relation, i.e., Eq. (1) (Yang and Liou, 1997; Bi et al., 2011b), which incorporated the use of fundamental electromagnetic theories for the mapping of the electric vector within a particle's interior to the far field described previously. Note, the IGOM can also be implemented on the basis of an electromagnetic surface-integral relation, i.e., Eq. (5). The IGOM results in Fig. 2 are computed with the algorithm reported in (Bi et al., 2011b). The CGOM asymmetry factor differs substantially from that obtained from the IGOM when the particle size parameters are on the order of 20-50. Furthermore, a significant difference between the two approaches is revealed in computing the extinction efficiency. As articulated by (van de Hulst, 1957), the CGOM extinction efficiency is given by a value of 2 as a result of the energy equivalence between the contributions of diffraction and ray reflections/refractions. However, when complementary electromagnetic theories are incorporated in the ray-tracing analysis, interferences between diffraction and optical rays occur, leading to oscillations of the scattered field and producing fluctuations in the results of the asymmetry factor and extinction efficiency.

Figure 2.  Comparison of the asymmetry factor, the extinction efficiency, P₁₁ (phase function), and P₁₂/P₁₁ and P₂₂/P₁₁ computed from the CGOM and the IGOM with the exact solutions (Bi and Yang, 2014a). Note that the IGOM used in the simulation is implemented on the basis of Eqs. (1) to (4) for mapping the near field computed with geometric optics to the corresponding far field.

The conventional CGOM halos were found to be independent of ice crystal size based on the phase function results. Results from IGOM, however, illustrate that halo patterns (22^° and 46^°) become broad and eventually disappear as the size parameter decreases in association with the ray-spreading effect for the light beams with finite cross sections. Subtle oscillations occurring in P₁₂/P₁₁ and P₂₂/P₁₁ are produced by the interference among scattered rays, and large differences in P₂₂/P₁₁ are shown in the backscattering directions.

The geometric ray-tracing approach is based on the concept of the localization principle, which does not account for rays tangential to the particle's surfaces and those traveling around the particle, whose optical effect is referred to as the edge effect. This complicated effect and the related surface waves are beyond the scope of geometric-optics applicability. For example, the geometric ray tracing can be used to explain rainbows and rainbow-type features produced by the scattering of light by spherical raindrops; however, the approach fails to generate the glory feature that occurs in the 180^° backscattering direction (van de Hulst, 1957).

Nussenzveig (1979, 1992) and (Nussenzveig and Wiscombe, 1980) developed the complex angular momentum theory, a mathematical analysis based on the asymptotic behavior of the exact Lorenz-Mie solution for the physical explanation of glory and rainbow features. (Liou et al., 2010) followed the edge effect formulations given in (Nussenzveig and Wiscombe, 1980) and incorporated the contribution into a volume-integral-based IGOM (Yang et al., 1997), the geometric-optics surface-wave (GOS) approach, the single-scattering results of which illustrate close agreement with those computed from the exact Lorenz-Mie theory. Moreover, (Liou et al., 2011) applied the GOS approach to nonspherical particles by developing a correction factor using a nondimensional volume parameter of 1 for spheres and 0 for elongated particles. By using this principle, the extinction coefficient, single-scattering albedo, and asymmetry factor results for randomly oriented columns and plates compare reasonably well with those computed from the numerical methods that solve Maxwell's equations (Liou et al., 2011). (Bi et al., 2010) and (Bi and Yang, 2014b) have discussed the importance of the edge effect in the cases of circular cylinders and spheroids. For spheroids, (Bi and Yang, 2014b) extended the extinction formula obtained from the complex angular momentum theory (Nussenzveig and Wiscombe, 1980) from spheres to spheroids by modifying the divergence factor of the central ray and the radius of curvature in the computation of the edge effect. In comparison with the T-matrix results, the extended formula is reasonably accurate when the size parameter is larger than ∼20 and can even be accurate up to three or four decimal points when the size parameter is larger than ∼ 200. Moreover, (Takano et al., 2013) demonstrated that by employing the Rayleigh-Gans-Debye (RGD) adjustment, which includes a form factor developed by (Debye, 1915), the GOS approach could be efficiently applied to the computation of the optical properties of black carbon aggregates.

In summary, within the last 30 years a spectrum of geometric-optics approximations has been developed with respect to the size parameter range, including the conventional geometric-optics method (CGOM, geometric ray tracing plus diffraction), the improved geometric-optics method (IGOM, the mapping of the ray-tracing results for the near field within the particle's interior or on the particle's surface to the far field based on the exact volume or surface electromagnetic integrals), and a form of the IGOM with the inclusion of the edge effects. In the developments of the aforementioned approximate methods based on physical principles and rationale, the results computed by other more rigorous and accurate methods were utilized, which are applicable for certain size parameter ranges discussed in the following sub-section.

The part reflections and refractions play within the context of geometric optics is based on the principle of localization, i.e., optical rays can be identified on the incident wavelength. However, it should be noted that the localization principle is valid only in the case of particles much larger than the incident wavelength. (van de Hulst, 1957) articulated "A pencil of light of length l can exist only if its width at its base is large compared to \(\sqrt\lambda l\)" where Λ is the incident wavelength. Thus, based on the limited applicability of the localization principle postulated by (van de Hulst, 1957), the conventional geometric-optics method may not be applicable to the scattering of light by a particle if the size parameter is less than approximately 20. For particles of moderate size parameters in the range of \(\kappa=20\)-50, the accuracy of geometric-optics solutions was checked with the following three numerical approaches, which provide "numerically exact" solutions: the finite-difference time domain (FDTD) method (Yee, 1966; Yang and Liou, 1995; Yang and Liou, 1996a; Sun et al., 1999), pseudo-spectral time-domain (PSTD) method (Liu, 1997; Chen et al., 2008; Liu et al., 2012a), and discrete dipole-approximation (DDA) method (Purcell and Pennypacker, 1973; Draine and Flatau, 1994; Yurkin and Hoekstra, 2007).

Pioneered by(Yee, 1966),\,the FDTD first\,solves\,Maxwell's equations in the time domain as follows: \begin{eqnarray} \label{eq5} \nabla\times \bm{H}(\bm{r},t)&=&\dfrac{\varepsilon(\bm{r})}{c}\dfrac{\partial\bm{E}(\bm{r},t)}{\partial t} ,\\ \label{eq6} \nabla\times \bm{E}(\bm{r},t)&=&-\dfrac{\mu(\bm{r})}{c}\dfrac{\partial \bm{H}(\bm{r},t)}{\partial t} , \end{eqnarray} where E and H denote electric and magnetic vectors, respectively; μ is the permeability; and c is the speed of light in a vacuum. The total field can be obtained by the superposition of incident and scattered fields given by \begin{eqnarray} \label{eq7} \nabla\times \bm{H}_{\rm sca}(\bm{r},t)&=&\dfrac{{\rm Re}[\varepsilon(\bm{r})]}{c}\dfrac{\partial \bm{E}_{\rm sca}(\bm{r},t)}{\partial t}+\nonumber\\ &&k{\rm Im}[\varepsilon(\bm{r})]\bm{E}_{\rm sca}(\bm{r},t)+\bm{S}(\bm{r},t) ,\quad\\ \label{eq8} \nabla\times \bm{E}_{\rm sca}(\bm{r},t)&=&-\dfrac{1}{c}\dfrac{\partial \bm{H}_{\rm sca}(\bm{r},t)}{\partial t} ,\\ \label{eq9} \bm{S}(\bm{r},t)&=&\dfrac{[1-{\rm Re}(\varepsilon)]}{c}\dfrac{\partial \bm{E}_{\rm inc}(\bm{r},t)}{\partial t}-\nonumber\\ &&k{\rm Im}(\varepsilon)\bm{E}_{\rm inc}(\bm{r},t) , \end{eqnarray} where S represents the source function. Note, the permeability is assumed to be unity in Eq. (9), a reasonable assumption for ice particles. The solution of these equations requires that the computational space domain be finite so that simulations can be implemented numerically. Figure 3a depicts a conceptual diagram to illustrate the principle of the FDTD technique developed by (Yee, 1966).

First, an absorbing boundary box is designed to enclose the particle such that any wave incident on the boundary can be absorbed without reflection. A number of absorbing boundary conditions have been reviewed in (Yang and Liou, 2000). The space between the particle and absorbing boundary layer is called "white space". The electric and magnetic fields within the particle and white space are the same as those in an unbounded scattering medium.

Second, Eqs. (8) to (10) must be discretized for numerical implementation by assigning the electromagnetic field on Cartesian grid cells in the computational domain. Each grid cell is called a Yee cell, as shown in Fig. 3b, where the electric field is assigned at the cube edges, while the magnetic field is defined at the center of the faces. In order to compute the derivative of electric and magnetic fields with respect to the space, the FDTD employs a finite-difference technique, while the PSTD (Liu, 1997; Panetta et al., 2013) adopts the spectrum method. For example, \begin{eqnarray} \label{eq10} &&{\rm FDTD}: [\nabla\times\bm{E}_{\rm sca}(\bm{r}_{i+1/2,j,k},t)]_x\nonumber\\[1mm] &\!\!=\!\!&\dfrac{1}{\Delta z}\\}\bm{E}_{{\rm sca},y}(\bm{r}_{i+1/2,j,k+1/2},t)- \bm{E}_{{\rm sca},y}(\bm{r}_{i+1/2,j,k-1/2},t)\}+\nonumber\\[1mm] &&\dfrac{1}{\Delta y}\\}\bm{E}_{{\rm sca},z}(\bm{r}_{i+1/2,j+1/2,k},t)\!-\! \bm{E}_{{\rm sca},z}(\bm{r}_{i+1/2,j-1/2,k},t)\} ,\qquad \end{eqnarray} \begin{eqnarray} \label{eq11} &&{\rm PSTD}: [\nabla\times\bm{E}_{\rm sca}(\bm{r}_{i+1/2,j,k},t)]_x\nonumber\\[1mm] &\!\!=\!\!&\bm{F}^{-1}\\}-ik_z\bm{F}[\bm{E}_{{\rm sca},y}(\bm{r},t)]\}\!-\! \bm{F}^{-1}\\}-ik_y\bm{F}[\bm{E}_{{\rm sca},z}(\bm{r},t)]\} ,\qquad \end{eqnarray} where F and F^-1 indicate forward and inverse Fourier transformations, respectively.

Figure 3.  (a) Conceptual FDTD diagram for the computation of the near field for the scattering of an electromagnetic wave by a nonspherical bullet rosette ice crystal. (b) Configuration for the electric and magnetic fields on a grid cell used in the Yee algorithm (Yee, 1966).

Third, the FDTD/PSTD derives direct solutions of the near field in the time domain. The electric and magnetic field in the frequency domain is obtained by applying the discrete Fourier transform. Having the near field defined, the optical properties are subsequently computed from Eqs. (1) to (5) by means of numerical integrations. The computational capabilities of the PSTD to model the optical properties of ice crystals have been described in Liu et al. (2012a, 2012b) and (Panetta et al., 2013).

Different from the FDTD/PSTD, the DDA directly solves field equations in the frequency domain such that the computational domain is the particle volume. Similar to the FDTD, the DDA, in principle, is applicable to arbitrarily shaped particles, including those without rotational symmetry or having inhomogeneous chemical composition. (Purcell and Pennypacker, 1973) introduced the DDA to study the optical properties of interstellar dust particles. However, the basic idea of the DDA can be traced to the work of (DeVoe, 1964), who studied the optical properties of molecular aggregates. As shown in Fig. 4a, the DDA assumes a particle to be a finite array of dipoles. With an incident plane wave, each dipole responds to the incident wave, then immediately emits radiation that interacts with the other dipoles (an example of a three-dipole system is displayed in Fig. 4b) until the local field for each dipole, which includes the incident field plus the fields produced by other dipoles, is stable. Consider the ith dipole having an electric dipole moment P_i, and its contribution to the radiation field at the position of the jth dipole is given by (Jackson, 1975) \begin{eqnarray} \label{eq12} \bm{E}(\bm{R}=\bm{r}_j-\bm{r}_i)&=&\bigg\\}k^2(\hat{\bm{R}}\times \bm{P}_i)\times\hat{\bm{R}}+[3\hat{\bm{R}}(\hat{\bm{R}}\cdot\bm{P}_i)-\bm{P}_i]\times\quad\nonumber\\ &&\left(\dfrac{1}{\bm{R}^2}-\dfrac{ik}{\bm{R}}\right)\bigg\}\dfrac{e^{ik\bm{R}}}{\bm{R}} , \end{eqnarray} where \(\hat\bmR\) is a unit vector pointing along position vector R. The vector P_i is proportional to the local field, leading to a self-consistent equation as follows: \begin{equation} \label{eq13} \bm{P}_i=\alpha_i\bm{E}_i=\alpha_i\left(\bm{E}_{{\rm inc},i}+\sum_{i\ne j}\bm{A}_{ij}\cdot \bm{P}_j\right) , \end{equation} where α_i is polarizability and A_ij is implicitly defined by comparing Eqs. (12) and (13) with an explicit form given by \begin{eqnarray} \label{eq14} \overline{\bm{A}}(\bm{R}=\bm{r}-\bm{r}')&=&\dfrac{k^2\exp(ikR)}{R}\bigg[\left(\overline{\bm{I}}-\dfrac{\bm{RR}}{R^2}\right)-\quad\nonumber\\[1mm] &&\dfrac{1-ikR}{(kR)^2}\left(\overline{\bm{I}}-\dfrac{3\bm{RR}}{(kR)^2}\right)\bigg] , \end{eqnarray} where \(\overline\bmI\) is the identity dyad. The polarizability is related to the refractive index \(\tildem\) through the Clausius-Mossotti relation (Jackson, 1975) in the form \begin{equation} \label{eq15} \alpha=d^3\dfrac{3}{4\pi}\dfrac{\tilde{m}^2-1}{\tilde{m}^2+2} , \end{equation} where \(\tildem\) is the complex refractive index, and d is the dipole length. Once the amplitudes of dipole moments denoted in Eq. (14) are solved, the extinction and absorption cross sections can be obtained from volume integrals given by Eqs. (3) and (4) as follows (Draine and Flatau, 1994): \begin{eqnarray} \label{eq16} \sigma_{\rm ext}&=&\dfrac{4\pi k}{|\bm{E}_{\rm inc}|^2}\sum_{j=1}^N {\rm Im}(\bm{E}_{{\rm inc},j}^\ast\cdot\bm{P}_j) ,\\ \label{eq17} \sigma_{\rm abs}&=&\dfrac{4\pi k}{|\bm{E}_{\rm inc}|^2}\sum_{j=1}^N \\}{\rm Im}(\bm{E}_j^\ast\cdot\bm{P}_j)\} . \end{eqnarray}

Figure 4.  (a) Schematic geometry to illustrate the approximation of a bullet rosette in terms of discretized dipoles. (b) An example illustrating the interactions of three dipoles.

According to Eq. (1), the scattered far field is given by \begin{equation} \label{eq18} \bm{E}_{\rm sca}(\bm{r})|_{r\to\infty}=\dfrac{\exp(ik\bm{r})}{-ik\bm{r}}(-ik^3)(\overline{\bm{I}}-\hat{\bm{r}}\cdot\hat{\bm{r}})\cdot \sum_j \bm{P}_j\exp(ik\bm{r}_j\cdot\hat{\bm{r}}) . \end{equation}

As shown in (Purcell and Pennypacker, 1973), Eq. (14) is derived on the basis of physical rationale rather than rigorous mathematical analysis. In later developments, however, the DDA was more rigorously refined by starting from an electric volume-integral equation in the form \begin{equation} \label{eq19} \bm{E}(\bm{r})=\bm{E}_0(\bm{r})+\int_V k^2(\varepsilon-1)\overline{\bm{G}}(\bm{r}-\bm{r}')\cdot\bm{E}(\bm{r})dv , \end{equation} where \(\overline\bmG\) is the dyadic Green's function satisfying the following differential equation: \begin{equation} \label{eq20} \nabla\times\nabla\times\bm{G}(\bm{r}-\bm{r}')-k^2\bm{G}(\bm{r}-\bm{r}')=\overline{\bm{I}}\delta(\bm{r}-\bm{r}') , \end{equation} where δ(r-r') is the Dirac delta function. The solution of the dyadic Green's function is given by \begin{equation} \label{eq21} \overline{\bm{G}}(\bm{R}=\bm{r}-\bm{r}')=\left[\overline{\bm{I}}+\dfrac{\nabla\nabla}{k^2}\right]\dfrac{\exp(ik\bm{R})}{4\pi\bm{R}}= \dfrac{\overline{\bm{A}}(\bm{R}=\bm{r}-\bm{r}')}{4\pi k^2} . \end{equation}

The dyadic Green's function in the source region has a singularity problem associated with the 1/R³ dependence term defined by \begin{equation} \label{eq22} \overline{\bm{G}}(\bm{R}=\bm{r}-\bm{r}')=\dfrac{1}{4\pi k^2}\nabla\nabla\dfrac{1}{R}=\dfrac{1}{4\pi R}\left[-\dfrac{1}{k^2R^2}- \hat{\bm{R}}\hat{\bm{R}}\dfrac{3}{k^2R^2}\right] . \end{equation}

To facilitate numerical calculations, the electric field is assumed to be constant within a small volume so that Eq. (20) can be written as (Yurkin and Hoekstra, 2007) \begin{eqnarray} \bm{E}(\bm{r})&=&\bm{E}_0(\bm{r})+\int_{V-V_\varepsilon}k^2(\varepsilon-1)\overline{\bm{G}}(\bm{r}-\bm{r}')\cdot \bm{E}(\bm{r}')dv+\nonumber\\ &&(\varepsilon-1)[\overline{\bm{M}}(\bm{r})+\overline{\bm{L}}(\bm{r})]\cdot\bm{E}(\bm{r}), \end{eqnarray} where \begin{eqnarray} \label{eq23} \overline{\bm{M}}(\bm{r})&=&k^2\left\\}\int_{V_\varepsilon}[\overline{\bm{G}}(\bm{r}-\bm{r}')-\overline{\bm{G}}^{s}(\bm{r}-\bm{r}')]d^3\bm{r}'\right\} ,\\ \label{eq24} \overline{\bm{L}}(\bm{r})&=&k^2\left\\}\int_{V_\varepsilon}\overline{\bm{G}}^s(\bm{r}-\bm{r}')d^3\bm{r}'\right\} . \end{eqnarray}

Based on the divergence theorem (Morse and Feshbach, 1953), Eq. (25) is equivalent to \begin{equation} \label{eq25} \overline{\bm{L}}=-\dfrac{1}{4\pi}\int_s\dfrac{\hat{\bm{n}}\hat{\bm{R}}}{R^2}ds , \end{equation} where s is the surface of an elementary volume V_ε. If the elementary volume is a cube or a sphere, it can be shown that \(\overline\bmL=-\overline\bmI/3\). According to Eq. (23), the source of the field (or the exciting field) at the position r (i.e., the contribution from incident and scattered fields from the remaining volume) is given by \begin{equation} \label{eq26} \bm{E}_{\rm exc}(\bm{r})=[\overline{\bm{I}}-(\varepsilon-1)(\overline{\bm{M}}+\overline{\bm{L}})]\cdot \bm{E}(\bm{r}) . \end{equation}

The amplitude of the dipole moment associated with the volume V_ε is governed by \begin{equation} \label{eq27} \bm{P}=\dfrac{V_\varepsilon(\varepsilon-1)}{4\pi[\overline{\bm{I}}-(\varepsilon-1)(\overline{\bm{M}}+\overline{\bm{L}})]}\bm{E}_{\rm exc}(\bm{r}) . \end{equation}

The polarizability can now be defined as: \begin{equation} \label{eq28} \alpha=\dfrac{V(\varepsilon-1)}{4\pi[\overline{\bm{I}}-(\varepsilon -1)(\overline{\bm{M}}+\overline{\bm{L}})]} . \end{equation}

(Lakhtakia, 1992) proposed strong and weak forms of the DDA to distinguish the treatment of a singularity problem in the dyadic Green function. The weak form corresponds to the Clausius-Mossotti relation (i.e., \(\overline\bmM=\overline\textbf0\)), whereas the strong form considers nonzero \(\overline\bmM\). For the strong form of DDA, several polarizability relations (Yurkin et al., 2007) are formulated replacing the simple Clausius-Mossotti relationship, leading to different approximations for \(\overline\bmM\). Polarizability relations based on different physical postulations are also discussed in (Yurkin et al., 2007).

The DDA can be used to efficiently determine the scattering and absorption properties when \(\kappa<\sim 20\). The computational resources (such as memory and the number of processors) increase dramatically with increasing the size parameter. However, with the development of state-of-the-art computing clusters, the optical properties for moderate-to-large size parameters can be obtained by means of parallel computations. For example, the optical properties of a sphere with a size parameter of 320 and a refractive index of 1.05 have been obtained by using the DDA with 512 processors and 698 Gb memory (Yurkin and Hoekstra, 2011). Due to the existing shape errors and the approximate nature of polarizability relations, the DDA is often viewed as an approximate or semi-exact method. The accuracy of the DDA method has been assessed (Penttila et al., 2007) through the comparison of DDA results with those computed from analytical or semi-analytical methods; e.g., the Lorenz-Mie theory for spheres, and the extended boundary condition method (EBCM) for spheroids and circular cylinders (Mishchenko et al., 2002).

The FDTD/PSTD and DDA have comparable accuracy. As an example, Fig. 5 shows a comparison of the six phase matrix elements of randomly oriented bullet rosette ice crystals computed with the PSTD and DDA. The refractive index used in the simulation is 1.3078+i10^-6 and the size parameter, defined with respect to the dimension of the circumscribed sphere, is 30. The two different numerical methods show quite similar numerical results. The slight differences may result from different numerical implementations and shape errors associated with particle geometry discretization.

Figure 5.  Comparison of the phase matrix of randomly oriented bullet rosette ice crystals computed with the PSTD and the DDA.

However, the computational efficiencies of the FDTD/ PSTD and DDA are quite different, particularly for large size parameters with large refractive indices. (Yurkin et al., 2007) reported a comparison of the DDA and FDTD in the case of spheres with the refractive index values ranging from 1.02 to 2 and the size parameters up to 80. The DDA was found to be faster than the FDTD when the refractive index is smaller than 1.4. The PSTD technique has the advantage over the FDTD because of coarse spatial resolution in particle discretization. (Liu et al., 2012b) and (Podowitz et al., 2014) reported a comparison between the numerical performance of the DDA and PSTD techniques. Their studies show that, for a fixed size parameter, a critical value of the refractive index exists above which the PSTD is faster than the DDA. Figure 6 shows a comparison of the computational efficiency of the PSTD and DDA based on the data presented in (Liu et al., 2012b). The DDA fails to produce convergent results in the case of very large refractive indices. For a fixed refractive index the ratio of computational time of the DDA over that of the PSTD was found to increase with respect to size parameter. From the aforementioned comparison studies, the time domain methods are more powerful than the DDA with reference to the handling of optically dense particles. Because of slow convergence or failure to converge, the DDA had not been used in the past when the refractive index was larger than 2. An attempt to extend the application of DDA to large refractive indices (>2) can be found in (Yurkin et al., 2010).

Figure 6.  Ratio of the computational time required by the PSTD to that by the DDA for various refractive indices.

The T-matrix method was first introduced by Waterman (1965, 1971) for solving acoustic and electromagnetic scattering problems. The method has been further improved and numerically implemented by a number of researchers (e.g., Wiscombe and Mugnai, 1986; Barber and Hill, 1990; Lakhtakia et al., 1984; Mackowski and Mishchenko, 1996; Mishchenko et al., 1996; Mishchenko and Macke, 1998; Doicu et al., 2006; Baran et al., 2001a, 2001b). In particular, (Mishchenko, 1991) proved that the T-matrix is independent of the orientation of the scattering particle with respect to the direction of the incident light, and this feature provides an exact and efficient approach to averaging the single-scattering properties of nonspherical particles over particle orientations. In the following, we briefly summarize the formalism of the T-matrix method.

Within the framework of the T-matrix method, the incident and scattered fields are expanded in terms of suitable vector spherical wave functions in the forms (Mishchenko et al., 2002) \begin{eqnarray} \label{eq29} \bm{E}_{\rm inc}&=&\sum_{l=1}^\infty a_l{\rm Rg}\bm{M}_l(k\bm{r})+b_l{\rm Rg}\bm{N}_l(k\bm{r}) ,\\ \label{eq30} \bm{E}_{\rm sca}&=&\sum_{l=1}^\infty p_l\bm{M}_l(k\bm{r})+q_l\bm{N}_l(k\bm{r}) , \end{eqnarray} where RgM_l and RgN_l are regular vector spherical functions, and M_l and N_l are irregular vector spherical functions. In Eqs. (30) and (31), l is an index defined to combine two indices m₁ and m₂ via m₁(m₁+1)+m₂, where m₁ is the total angular momentum (\(m_1=1,\infty\)) and m₂ is the projected angular momentum (m₂≤ m₁). Because of the linearity of Maxwell's equations and electromagnetic boundary conditions, the expansion coefficients of the scattered and incident fields are related by means of a transition matrix, called the T-matrix (symbolically) in the form \begin{equation} \label{eq31} \left[ \begin{array}{c} \bm{p}\\ \bm{q} \end{array} \right]=\bm{T}\left[ \begin{array}{c} \bm{a}\\ \bm{b} \end{array} \right] . \end{equation}

If the T-matrix is known for an arbitrary incident field, the corresponding scattered field can be directly obtained through a matrix product. Thus, the T-matrix contains all required information associated with scattering and absorption, from which the commonly defined single-scattering properties, including phase matrix, extinction efficiency and single-scattering albedo, can be derived in a relatively straightforward manner.

A critical component of the T-matrix algorithm for computing the optical properties of particles is the numerical capability to obtain the T-matrix from Maxwell's equations or their mathematical equivalents. Although the extended boundary condition method (EBCM) is a pioneering method for computing the T-matrix of a single homogeneous nonspherical particle, a number of alternative techniques have been suggested for different applications. In the simplest case of a sphere, the T-matrix is a diagonal matrix whose elements can be obtained from the method of separation of variables (i.e., the Lorenz-Mie theory). The T-matrix for a cluster of spheres can be obtained from the T-matrices of individual spheres by employing the addition theorem of vector spherical wave functions to account for mutual interactions among the spheres (Peterson and Str\"om, 1973; Mackowski and Mishchenko, 1996). An alternative approach for computing the T-matrix is the invariant imbedding T-matrix (II-TM) method (Johnson, 1988), which is stable and efficient for arbitrarily shaped homogeneous or inhomogeneous particles (Bi et al., 2013a, 2013b; Bi and Yang, 2014a). In particular, the II-TM can be combined with the separation of variables method, the EBCM, or the superposition principle to improve its computational efficiency, and the choice of a specific combination scheme depends on the particle geometry. Other methods for T-matrix calculations include the null-field method with discrete sources (e.g., Doicu et al., 2006), the EBCM assisted with group theory (Kahnert, 2013), the generalized point-matching method (Nieminen et al., 2003), the superposition method (Mackowski and Mishchenko, 1996, 2011), the discrete-dipole moment method (Mackowski, 2002), the iterative EBCM (Lakhtakia et al., 1984), and the many-body T-matrix for a single particle with a large aspect ratio (Yan et al., 2008; Sun et al., 2013). In the following, we briefly highlight the EBCM and the II-TM and focus on demonstrating the II-TM's applicability to the simulation of the optical properties of ice crystals.

The EBCM computes the T-matrix through integral equations over the particle's surface in the forms (Mishchenko et al., 2002) \begin{eqnarray} \label{eq32} \bm{E}_{\rm sca}(\bm{r})|_{\bm{r}\in {\rm medium}}&=&\oint ds(\bm{r}')\bigg[\dfrac{i\omega\mu}{c}(\hat{n}_+\times \bm{H}(\bm{r}'))\cdot \overline{\bm{G}}(\bm{r}',\bm{r})+\nonumber\\ &&(\hat{n}_+\times\bm{E}(\bm{r}'))\cdot(\nabla\times\overline{\bm{G}}(\bm{r}',\bm{r}))\bigg] ,\\ \label{eq33} \bm{E}_{\rm inc}(\bm{r})|_{\bm{r}\in {\rm particle}}&=&\!-\!\oint ds(\bm{r}')\bigg[\dfrac{i\omega\mu}{c}(\hat{\bm{n}}_+\!\times \bm{H}(\bm{r}'))\cdot \overline{\bm{G}}(\bm{r}',\bm{r})+\nonumber\\ &&(\hat{\bm{n}}_+\times\bm{E}(\bm{r}'))\cdot(\nabla\times\overline{\bm{G}}(\bm{r}',\bm{r}))\bigg] , \end{eqnarray} where \(\hat\bmn_+\) is the outward normal directions of the particle surfaces. The unknown surface fields are defined through boundary conditions to provide a connection to the internal field, which is expanded with respect to vector spherical functions as follows: \begin{equation} \label{eq34} \bm{E}_{\rm int}=\sum_{l=1}^\infty c_l{\rm Rg}\bm{M}_l(\tilde{m}k\bm{r})+d_l{\rm Rg}\textbf{\emph{N}}_l(\tilde{m}k\bm{r}) , \end{equation} where \(\tildem\) is the refractive index. The dyadic Green's function can be formulated with respect to the vector spherical functions. According to Eqs. (33) to (35), the relation between (a, b) and (p, q) can be established through the T-matrix equation in the form \begin{equation} \label{eq35} \bm{T}=-{\rm Rg}\bm{Q}\cdot \bm{Q}^{-1} , \end{equation} where RgQ and Q are defined by \begin{eqnarray} \label{eq36} \left[ \begin{array}{c} \bm{a}\\ \bm{b} \end{array} \right]={\rm Rg}\bm{Q}\left[ \begin{array}{c} \bm{c}\\ \bm{d} \end{array} \right] ,\\ \label{eq37} \left[ \begin{array}{c} \bm{p}\\ \bm{q} \end{array} \right]=-\bm{Q}\left[ \begin{array}{c} \bm{c}\\ \bm{d} \end{array} \right] . \end{eqnarray}

The EBCM has been shown to be an efficient method for the computation of randomly oriented axially symmetric particles, such as spheroids, circular cylinders, bispheres and Chebyshev particles (Mishchenko et al., 2002). Unlike the purely numerical methods (e.g. FDTD and DDA), the EBCM is semi-analytical and provides an efficient and stable algorithm for axially symmetric particles. Due to its high accuracy, the specific T-matrix method is often employed as an accurate reference to validate numerical methods and geometric-optics approximations and to examine the numerical performance. For example, the T-matrix method is employed to study the accuracy of the DDA with different numerical implementations (Penttila et al., 2007). (Yang et al., 2007) employed the T-matrix program developed by (Mishchenko and Macke, 1998) to study the accuracy of the IGOM for the single-scattering properties of randomly oriented spheroids with size parameters from the Rayleigh regime to the lower end (\(\kappa=50\)) of the geometric optics regime. However, employing the EBCM to compute the optical properties of complex ice particles is challenging because of the ill-conditioned problem inherent in numerical algorithms involving complicated particle geometries.

Different from the EBCM, the II-TM computes the T-matrix through Eq. (20) representing an electromagnetic volume-integral equation; see the pioneering work presented in (Johnson, 1988). In principle, the computational procedure of determining the T-matrix within the II-TM framework can be summarized as: (1) the nonspherical particle is assumed to be an inhomogeneous sphere circumscribing the nonspherical particle (note, the portion of the inhomogeneous sphere, exterior to the nonspherical particle, is empty, i.e., with a refractive index of 1——see the schematic diagram shown in Fig. 7a); (3) the inhomogeneous sphere in Fig. 7a is discretized into a number of thin layers (see Fig. 7b); and (4) the T-matrix of the particle can be obtained from inner spherical layers to outer spherical layers through an iterative procedure based on the following equation (Johnson, 1988; Bi et al., 2013a, 2013b): \begin{eqnarray} \label{eq38} \bm{T}(r_p)&=&\bm{Q}_{11}(r_p)+[\bm{I}+\bm{Q}_{12}(r_p)][1-\nonumber\\ &&\bm{T}(r_{p-1})\bm{Q}_{22}(r_p)]^{-1}\bm{T}(r_{p-1})[\bm{I}+\bm{Q}_{21}(r_p)] ,\quad \end{eqnarray} where I indicates the identity matrix, T(r_p-1) represents the T-matrix for the layer with a radius of r_p-1, and the Q_ij matrices are computed from surface integrals over the r_p spherical layer. Note, for a given spherical layer, only the portion associated with a non-unity refractive index contributes to the aforementioned integrals. Figures 7c and d schematically illustrate the II-TM iterative procedure for computing the T-matrix for a hexagonal particle. Figure 7c shows the intersection of a spherical layer and the hexagonal particle in the process of iteratively computing the T-matrix, starting from a sphere inscribed by the hexagonal particle to the nonspherical particle geometry of interest (Fig. 7d).

Figure 7.  (a) A nonspherical particle circumscribed by an imaginary inhomogeneous sphere. The portion of the sphere, exterior to the nonspherical particle, is empty, i.e., the corresponding refractive index is 1. (b) Discretization of the inhomogeneous sphere into multi-layers. (c) Intersection of one spherical layer with a hexagonal ice crystal. (d) Hexagonal ice crystal represented as an inhomogeneous sphere.

To illustrate the accuracy of the II-TM method, Fig. 8 shows a comparison of six non-zero phase matrix elements of a randomly oriented circular cylinder computed from the EBCM and the II-TM. The EBCM T-matrix code developed by (Mishchenko and Macke, 1998) and the II-TM code developed by (Bi et al., 2013a) are used in this study. For the results in Fig. 8, the size parameter defined with respect to the particle length or the diameter (note, L/2a=1 for the particle) is 60. The refractive index is assumed to be 1.3078+i1.66× 10^-8. The results from the two methods are virtually the same.

Figure 8.  Comparison of the phase matrix elements of randomly oriented circular cylinders computed with the\quad EBCM and the II-TM.

To illustrate the numerical applicability of the II-TM to the simulation of the optical properties of complex ice crystals, Fig. 9 shows the phase matrix elements of an aggregate ice crystal consisting of eight columns under the random-orientation condition. The II-TM code developed by Bi and Yang (2014a) is used in this study. The wavelength is assumed to be 0.66 μm and the particle maximum dimension is 6 μm. The \(\kappa\) defined in terms of the maximum dimension is approximately 57. The results from the two methods agree with slight differences due to numerical errors. In the DDA simulation, 128α angles, 33β angles, and 32γ angles are used to average over the random orientations of the particle, where α,β, and γ are Euler angles defined according to the z-y-z convention. The ADDA code developed by (Yurkin et al., 2007) is used in this study. Because the T-matrix method is more efficient than the FDTD/PSTD and DDA in performing the average of the optical properties over particle orientations, the T-matrix method is more computationally efficient in modeling randomly oriented ice particles.

Figure 9.  Comparison of six phase matrix elements of an aggregate consisting of eight solid columns computed with the DDA and the II-TM under randomly oriented conditions.

The term "remote sensing" refers to the use of the measurement of electromagnetic radiation based on passive (e.g., satellite imager) or active (e.g., radar) sensors to infer the physical properties of a remote object of interest. The sensors can take satellite-, aircraft-, or ground-based measurements. The inference of cloud parameters, such as cloud height or ice water path, is made based on the measurements. In general, retrieval approaches make use of radiative transfer simulations based on the atmospheric state and the bulk single-scattering properties of liquid or ice clouds. In other words, remote sensing of ice cloud characteristics can be viewed as an inverse problem to convert radiometric or polarimetric measurements into ice cloud properties. To accomplish this, one needs to incorporate fundamental light scattering and radiative transfer simulations into the inversion process (Minnis et al., 1993a, 1993b; Liou et al., 2000; Rolland et al., 2000; Rolland and Liou, 2001; Platnick et al., 2003; King et al., 2004; Li et al., 2005; Baum et al., 2005b; Minnis et al., 2011; Baum et al., 2014). Ice cloud properties that yield the best agreement between the numerical simulations and measurements are generally accepted as the retrieved properties.

To determine and understand ice cloud properties over large spatial and temporal scales, remote sensing based on satellite observations is the only practical choice. The parameters necessary to describe ice clouds include cloud top height/pressure/temperature and ice water path, which is obtained from the inferred optical thickness and particle effective size. These quantities play an essential role in determining the role of ice clouds in the hydrological cycle and on the Earth's energy budget. The analysis of satellite data generally begins with identifying those pixels (or fields of view) that are cloudy, and in particular clouds composed primarily of ice, which means that one must infer the cloud thermodynamic phase. Once an ice cloud is identified, the height of the cloud is determined, followed by estimating the ice water path (i.e., the cloud optical thickness and effective particle size). (King et al., 2004) presented some concise, yet insightful, flow-charts for retrieving cloud properties from high-altitude airborne observations, which are applicable to satellite-based retrievals as well. In this article, we concentrate on discussing the process to infer the ice water path that most closely uses the scattering models described.

A variety of methods have been reported in the literature to infer cloud optical thickness and particle size from satellite measurements. The most comprehensive satellite record is provided by the series of Advanced Very High Resolution Radiometer (AVHRR) imagers, which have been in continuous orbit on NOAA polar-orbiting platforms since 1978. The AVHRR sensors provide measurements at solar wavelengths and in infrared windows. Based on the AVHRR sensors, decadal records of the ice cloud properties are now available (Walther and Heidinger, 2012; Foster and Heidinger, 2013; Karlsson et al., 2013; Stengel et al., 2014). The most familiar polar-orbiting imagers are the Moderate Resolution Imaging Spectroradiometer (MODIS) (King et al., 1992; Platnick et al., 2003; Minnis et al., 2011; Baum et al., 2012) and the Visible Infrared Imaging Radiometer Suite (VIIRS), although others exist. Decades of measurements from geostationary imagers are also available (e.g., Rossow and Schiffer, 1999; Roebeling et al., 2006). Currently several space-borne IR spectrometers provide hyper-spectral information useful for ice cloud analyses (Huang et al., 2004; Li et al., 2005; Yue et al., 2007; Wang et al., 2011; Kahn et al., 2014). In addition, narrow-band infrared measurements have been used to retrieve ice cloud properties (Chiriaco et al., 2004; Garnier et al., 2012, 2013; Iwabuchi et al., 2014). Furthermore, satellite observations with polarization properties also show great potential for inferring ice cloud properties compared with those only measuring the intensity (Chepfer et al., 1998, 2001; C.-Labonnote et al., 2001; van Diedenhoven et al., 2012, 2013; Cole et al., 2013, 2014). What each of these retrieval studies has in common is the necessity for accurate light scattering calculations in support of the remote sensing of ice cloud properties. The critical need for consistent ice habit scattering properties over wide wavelength and large particle size ranges has provided impetus for much of the recent development in light scattering methodology. The solar-band bulk scattering properties used in remote sensing applications (Baum et al., 2011) are given by: \begin{eqnarray} \langle Q_{\rm ext}\rangle&\!\!=\!\!&\dfrac{\int_{\lambda_{\min}}^{\lambda_{\max}}\int_{r_{\min}}^{r_{\max}}[Q_{\rm ext}(\lambda,r)A(r)F(\lambda)S(\lambda)n(r)]drd\lambda} {\int_{\lambda_{\min}}^{\lambda_{\max}}\int_{r_{\min}}^{r_{\max}}[A(r)F(\lambda)S(\lambda)n(r)]drd\lambda}\,,\\ \langle Q_{\rm sca}\rangle&\!\!=\!\!&\dfrac{\int_{\lambda_{\min}}^{\lambda_{\max}}\int_{r_{\min}}^{r_{\max}}[Q_{\rm sca}(\lambda,r)A(r)F(\lambda)S(\lambda)n(r)]drd\lambda} {\int_{\lambda_{\min}}^{\lambda_{\max}}\int_{r_{\min}}^{r_{\max}}[A(r)F(\lambda)S(\lambda)n(r)]drd\lambda}\,,\\ \langle \varpi\rangle&\!\!=\!\!&\dfrac{\langle Q_{\rm sca}\rangle}{Q_{\rm ext}} ,\\ \langle g\rangle&\!\!=\!\!&\dfrac{\int_{\lambda_{\min}}^{\lambda_{\max}}\!\int_{r_{\min}}^{r_{\max}} \![g(\!\lambda,\!r\!)AQ_{\rm sca}(\!\lambda,\!r\!)(\!r\!)F(\!\lambda\!)S(\!\lambda\!)n(\!r\!)]drd\lambda} {\int_{\lambda_{\min}}^{\lambda_{\max}}\!\!\int_{r_{\min}}^{r_{\max}}[Q_{\rm sca}(\!\lambda,r\!)A(\!r\!)F(\!\lambda\!)S(\!\lambda\!)n(\!r\!)]drd\lambda}\,,\nonumber\\ \\ \langle P(\!\theta\!)\rangle&\!\!=\!\!&\dfrac{\int_{\lambda_{\min}}^{\lambda_{\max}}\!\!\int_{r_{\min}}^{r_{\max}} \![P(\!\theta,\!\lambda,\!r\!)Q_{\rm sca}(\!\lambda,\!r\!)A(\!r\!)F\!(\!\lambda\!)S\!(\!\lambda\!)n(\!r\!)]drd\lambda} {\int_{\lambda_{\min}}^{\lambda_{\max}}\!\!\int_{r_{\min}}^{r_{\max}}[Q_{\rm sca}(\!\lambda,r\!)A(\!r\!)F(\!\lambda\!)S(\!\lambda\!)n(\!r\!)]drd\lambda}\,,\nonumber\\ \end{eqnarray} where Λ_max and \(\lambda_\min\) are the upper and lower wavelength limits of a given channel of the satellite instrument; Q _ext and Q _sca are the extinction and scattering efficiencies; \(\varpi\) is the single scattering albedo; g is the asymmetry factor; P is the phase matrix; n(r) is the particle size distribution; F(Λ) is the response function of the instrument; S(Λ) is the solar spectrum; and, A is the projected area. For IR-band remote sensing applications, S(Λ) in the preceding equations needs to be replaced by the Planck function. Furthermore, for a given particle size distribution, the effective particle radius r _eff (note, the effective particle diameter is given by D _eff=2r _eff) is defined as (Foot, 1988): \begin{equation} r_{\rm eff}=\dfrac{3}{2}\dfrac{\Sigma_{h=1}^M[\int_{r_{\min}}^{r_{\max}}V_h(r)f_h(r)n(r)dr]} {\Sigma_{h=1}^M[\int_{r_{\min}}^{r_{\max}}A_h(r)f_h(r)n(r)dr]} , \end{equation} where M is the number of ice crystal habits; V_h and A_h are the volume and area of an ice crystal with the habit h; and f_h is the habit fraction of habit h.

In the following, we review the development of satellite-based remote sensing on ice clouds and introduce some well-established algorithms to infer cloud properties such as cloud top height/pressure, optical thickness, effective particle size, and particle shape. The retrievals based on solar-band and IR-band observations as well as those considering the polarization properties are discussed in detail.

The inference of cloud optical thickness and/or effective particle size (effective diameter, or D _eff) from solar wavelengths was suggested by (Hansen and Pollack, 1970), and further developed for application to aircraft measurements (e.g., Twomey and Cocks, 1982; King, 1987). The use of reflected intensity measurements at absorbing and non-absorbing wavelengths was reported by (Twomey and Cocks, 1989) and implemented as a bi-spectral algorithm by (Nakajima and King, 1990). The approach has been widely adopted for operational retrievals based on satellite imager measurements to infer ice cloud optical thickness and effective particle size. The approach uses a weakly absorbing, visible or near-infrared window band (VIS/NIR) and a moderately absorbing shortwave IR band (SWIR). The non-absorbing channel is mainly a function of optical thickness, but the latter is sensitive to both the effective particle size and the optical thickness. Operational retrievals make use of look-up tables (LUTs) that contain the transmission and absorption characteristics associated with a wide range of optical thickness and effective particle size values. With the LUTs, the cloud optical thickness and effective size are determined in the solution space under the constraint that the simulations most closely agree with the measurements. This assumes, of course, that the solution space does not permit multiple solutions for a given set of measurements.

Figure 10 shows an example of reflectivity LUTs for a pair of bands at 0.87 and 2.25 μm from the VIIRS (Lee et al., 2010; Hillger et al., 2013). The ice cloud bulk single-scattering properties used are similar to those used to derive the MODIS Collection 6 cloud product, i.e., the ice habit model assumes solely the use of ice crystal aggregates, each of which consists of eight solid columns (Baum et al., 2014). The solid curves represent the reflectivity values for specified values of effective particle size, and the dashed curves illustrate those of given values of cloud optical thickness. The reflectivity at 0.86 μm (the non-absorptive band) is clearly a strong function of optical thickness with little dependence on the effective particle size, whereas the reflectivity at 2.25 μm, in contrast, is sensitive to both the parameters at low values of the optical thickness, but primarily a function of the particle size for larger values of the optical thickness.

Because solar bands are available on many satellite imagers, the solar reflectivity algorithms have been widely employed for ice cloud property retrievals. Examples include MODIS (Platnick et al., 2003; Minnis et al., 2011), AVHRR (Heidinger et al., 2005), Spinning Enhanced Visible and Infrared Imager (SEVIRI) (Roebeling et al., 2006), and Polarization and Directionality of the Earth's Reflectances (POLDER) (Buriez et al., 2005).

However, there are limitations to the solar reflectivity approach. As illustrated in Fig. 10, the isolines converge as clouds become optically thin (optical thickness less than 1), and the reflectivity becomes sensitive to the surface conditions (i.e., surface reflection characteristics), which indicates that the method becomes less robust for optically thin clouds. The single-scattering properties of ice clouds in the solar bands are sensitive to particle microphysical properties, e.g., particle size, shape, and surface texture (i.e., the degree of particle surface roughness). The technique is limited to daytime measurements and uncertainties increase at high solar zenith angles and over bright surfaces (non-vegetated surfaces such as desert, snow, and ice). Further uncertainties are introduced by shadows resulting from 3D cloud effects where a cloud may have bright (i.e., illuminated) and shadowed sides if the sun is not directly overhead. Furthermore, the vertical variations of ice crystal habit and size may significantly complicate the retrievals of ice cloud properties (Yang et al., 2001).

Figure 10.  Reflectivity look-up tables for ice clouds at the VIIRS 0.87- and 2.25-μm bands. Surface albedo is 0.05 for both bands, and solar and viewing zenith angles are 60^° and 0^°. Solid curves are isolines of specified particle effective size (D _eff), and dashed curves are isolines with constant optical thickness (τ).

A case study is performed to illustrate the performance and sensitivity of the solar reflectivity approach. MODIS measurements, i.e., reflectivities for bands 2 and 7 (with central wavelengths at 0.86 and 2.13 μm, respectively), are used to retrieve ice cloud optical thickness and effective diameter. An Aqua MODIS granule (5 min of data) recorded at 0820 UTC 10 December 2013 is used for the case study. Figure 11 shows the true color image for the scene over the Indian Ocean, which is mostly covered by ice clouds. The retrievals are carried out for pixels over the ocean and identified as ice clouds in the MODIS Level 2 Collection 5.1 cloud product (MYD06).

Figure 11.  Aqua MODIS granule over the Indian Ocean at 0820 UTC 10 December 2013.

Figure 12 illustrates the retrieved cloud optical thickness and particle effective diameter based on the single-scattering properties assuming the sole use of either ice spheres or hexagonal solid columns. The hexagonal column model scattering properties used in this study are taken from (Bi and Yang, 2014a). The upper panels are the results based on the ice sphere model, and the lower panels are those from the hexagonal column model. The use of ice spheres results in the retrieval of larger cloud optical thickness values than those from hexagonal columns, but the opposite is true for the effective diameters. The results based on the two models are quite different.

Figure 12.  Retrieved cloud optical thickness (COT) and particle effective diameter (D _eff) based on the scattering properties calculated using the sphere (upper panels) and hexagon (lower panels) models.

To better interpret the results in Fig. 12, Fig. 13 shows the histograms of occurrence, on a logarithmic scale, for retrievals based on the two different ice habit (sphere and hexagon) models, with red indicating the highest frequencies of occurrence. The black 1:1 lines are included to ease interpretation of the results. The sphere model is shown to systematically give larger optical thickness values compared with those from the hexagonal solid column model. The comparison of the effective diameters is less straightforward, and provides both under- and overestimates compared with the results from the solid columns. The results clearly show the importance of the assumed ice habits and their scattering properties on inferring cloud properties using the solar band retrievals, and thus indicates the importance of developing more accurate and practical ice crystal habit models to determine the corresponding ice cloud optical properties.

Figure 13.  Comparison of retrieved ice cloud optical thicknesses and particle effective diameters based on the scattering properties calculated using the sphere and hexagon models.

In addition to the aforementioned retrieval algorithms and related bands, we wish to briefly elaborate on the usefulness of a near-infrared band. Specifically, the 1.38-μm water-vapor-absorbing band is unique and quite useful for detecting clouds located high in the atmosphere (Gao et al., 1993; Gao and Kaufman, 1995; Gao et al., 1998; Roskovensky and Liou, 2003a, 2003b). An operational algorithm has been developed to infer cirrus cloud reflectivity with the MODIS 1.38-μm channel (Gao et al., 2002). The MODIS 1.38-μm channel measurements have been used to estimate the spatial distribution and optical/radiative properties of thin cirrus clouds (Dessler and Yang, 2003; Meyer et al., 2007a, 2007b; Lee et al., 2009; Sun et al., 2011). Because of strong water vapor absorption at this wavelength, little outgoing radiance from either low clouds or the surface reaches the satellite, whereas radiance reflected by the high-level clouds, above which there is little water vapor, experiences little attenuation. Furthermore, the wavelength is long enough that Rayleigh scattering from the atmosphere is negligible. (Gao et al., 2004) used a combination of the 1.38- and 1.88-μ m bands to infer the optical thickness and particle effective size. (Roskovensky et al., 2004) and (Roskovensky and Liou, 2005) used the 1.38-μm channel to study co-existing aerosol (e.g., airborne dust) and cirrus cloud. (Wang et al., 2012) used the 1.24- and 1.38-μm bands to study cloud top height using MODIS observations, and more recently, (Wang et al., 2014) used the 1.38-μm channel to infer the phase function of ice clouds.

Considering the temporal and spatial limitation of solar-band-based retrievals, the use of IR bands can provide complementary information on ice clouds that is independent of solar illumination. In the IR, the influence of any individual habit is reduced over that in the shortwave because of increased absorption within the particles (e.g., Wendisch et al., 2007; Baum et al., 2014), which increases the benefit of using IR channels. The split-window method (Inoue, 1985; Parol et al., 1991; Giraud et al., 1997; Chiriaco et al., 2004) exploits the difference in ice absorption in the IR window bands (i.e., 8.5, 11, and 12 μm), and has been applied to Infrared Interferometer Spectrometer (IRIS) (Prabhakara et al., 1988), AVHRR (Heidinger and Pavolonis, 2009), and MODIS observations (Wang et al., 2011). Note, the 11-μm band is located within one of the Christian bands of the refractive index of ice, where the extinction of light by ice crystals is predominately due to absorption rather than scattering (Arnott et al., 1995; Yang et al., 1997). The primary limitation is that the IR radiances become saturated as the ice cloud optical thickness becomes larger than about 6. Nevertheless, the method can provide a very useful complement to solar-band retrievals.

To illustrate how the split-window method works, Fig. 14 shows the relationship between the simulated brightness temperature (BT) of three VIIRS IR-window bands with central wavelengths of approximately 8.5, 11, and 12 μm, and the associated brightness temperature differences (BTD) between two of the channels. The simulations are based on the US standard atmosphere containing ice clouds based on different optical thickness and effective particle diameter values, and the clear-sky and cloud top temperatures are 288 K and 223 K, respectively. A viewing zenith angle of 20^° is assumed. The solid curves are isolines of specified effective particle size, and the dashed curves are isolines of constant optical thickness. The BTs are mainly dependent on optical thickness, whereas the BTDs are sensitive to both the optical thickness and the effective particle size. The figure indicates that the split-window method is more sensitive to optically thin clouds (optical thickness between 0.5 and ∼6) with relatively small particle sizes (effective diameter less than ∼50 μm). As the optical thickness becomes larger than 10, the isolines converge, and the method is not applicable in practice.

Figure 14.  Relationships between simulated brightness temperature (BT) and brightness temperature difference (BTD) with respect to the VIIRS 8.5-, 11-, and 12-μm bands.

A case study of the IR-band-based retrieval is performed on the same MODIS granule as that used for the solar-band retrieval, and Fig. 15 is similar to Fig. 12 but for retrieval results based on three MODIS IR bands (i.e., bands 29, 31, and 32 with central wavelengths at 8.5, 11, and 12 μm). For the IR retrieval, only optically thin pixels (with optical thickness inferred by the MOD06 product of less than 10) are selected. The same hexagonal column and sphere models are used for the retrievals. Note that the IR bands infer most pixels to be optically thin, i.e., optical thicknesses less than 5, and the particle effective diameters are smaller than those given by the solar-band-based retrievals. Furthermore, different from the solar-band-based retrievals, the IR retrievals given by the two models show not only similar overall patterns but also amplitude, which indicate that the IR-band-based retrieval is less sensitive to the optical properties of the ice models.

Figure 15.  As in Fig. 12 but for retrieval results based on an IR-band-based algorithm.

For a better interpretation of the IR-band-based retrieval, Fig. 16 compares the retrieved optical thickness and particle effective diameter values based on the hexagon and sphere models, and the figure is organized in the same form as Fig. 13. The optical thicknesses given by the two models are clearly shown to agree much better than those derived with the solar-band-based retrievals, and the spherical model gives optical thicknesses slightly smaller than those associated with the hexagonal model. The retrieved effective particle size is distributed over a much wider range, which is similar to the solar-band-based retrievals.

Figure 16.  As in Fig. 13 but for retrieval results based on an IR-band-based algorithm.

Previous comparisons focus on retrievals based on either solar or IR bands but different ice cloud models, and Fig. 17 compares the occurrence of thin ice cloud pixels retrieved based on the solar and IR bands. The x-axis is for the optical thickness given by the solar-band retrieval, and the y-axis is for the corresponding IR-band retrieval with the same ice model. The red color corresponds to the high frequency of occurrence. Although the same ice cloud model is used for the solar- or IR-band-bases retrievals, the solar-band retrieved optical thicknesses are higher than those from the IR-band retrievals based on either the hexagon model or the sphere model. The difference between the optical thicknesses retrieved from different wavelengths further demonstrates a pressing need to develop more realistic and robust ice cloud models and, subsequently, to retrieve consistent corresponding single-scattering property results derived from different spectral methods. To this end, (Liu et al., 2014) developed a novel two-habit ice crystal habit model for the optical properties of ice clouds, which leads to spectrally consistent retrievals between the visible/near-infrared method and the thermal infrared method.

Figure 17.  Comparisons of ice cloud optical thicknesses retrieved from the solar-band-based and IR-band-based algorithms. The results are based on the single hexagonal column model (left panel) and the sphere model (right panel).

Furthermore, in the cloud retrieval community, implementing an optimal estimation approach (e.g., Rodgers, 2000) to infer different parameters is becoming increasingly common. A primary benefit of the optimal estimation approach is providing an estimate of the retrieval uncertainties as part of the retrieval process. The process provides a measure of the simulations to fit the measurements through minimization of a cost function, and uses error estimates of the various algorithm input parameters. The optimal estimation approach has been applied to different sensors, such as the AVHRR (Heidinger and Pavolonis, 2009), the Along Track Scanning Radiometer (ATSR) (Poulsen et al., 2012), the SEVIRI (Watts et al., 2011), and MODIS (Iwabuchi et al., 2014).

This sub-section briefly discusses the benefits of lidar and polarimetry, both of which measure ice cloud polarization properties. One limitation of passive satellite sensors is their limited ability to infer the vertical structure in an atmospheric column that contains multiple cloud layers. Active remote sensing by lidar offers a unique capability to observe a vertical profile of the atmosphere. The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) is on board the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) (Winker et al., 2009) platform. CALIOP is a dual-wavelength backscatter lidar (0.532 and 1.064 μm) and is the first polarization (perpendicular and parallel at 0.532 μm) sensor to provide global high-resolution vertical profiles of aerosols and clouds. Energy backscattered from the atmosphere is received for the 0.532 μm perpendicular and parallel bands and the 1.064 μm band, and the attenuated backscatter and depolarization ratio are calibrated for aerosol/cloud detection, scene classification, and various cloud and aerosol parameters. Of pertinence to this article, (Hu et al., 2009) found a relationship between the depolarization ratio and the backscatter coefficient to discriminate between randomly oriented and horizontally oriented planar ice crystals. Ice clouds containing even a small fraction of horizontally oriented ice crystals tend to have small values of the depolarization ratio but large values of backscatter. Following (Hu et al., 2009), several studies were conducted to use the relationship between layer-integrated depolarization ratio and backscatter to investigate cloud thermodynamic phase and the percentage of horizontally oriented ice plates in cirrus clouds (Cho et al., 2008; Noel and Chepfer, 2010; Zhou et al., 2012, 2013). Furthermore, the observations by active sensors provide an effective way to validate the retrievals based on passive sensor observations (e.g., Cho et al., 2009).

As another unique instrument in the NASA A-Train satellite constellation, the POLDER on board the Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar (PARASOL) is a passive optical imaging radiometer and polarimeter that measures I, Q, and U components of the Stokes parameter at three wavelengths with up to 16 viewing angles for each pixel (Deschamps et al., 1994). The normalized polarized reflectivity, which is specified by observed linearly polarized radiance and satellite-viewing geometry, is used to represent the polarization properties of the target. Chepfer et al. (1998, 2001) showed that polarized reflectivities from an ice cloud have some sensitivity to the ice habit, and simulations based on various ice habits were tested to fit the measured polarized reflectivities from the POLDER observations. Based on POLDER, the polarized reflectivity measurements have become an important tool to guide ice cloud scattering models (C.-Labonnote et al., 2001; Baran and C.-Labonnote, 2006; Cole et al., 2013, 2014).

Polarized reflectivities become saturated for ice clouds with optical thickness values larger than 5, at which point the scattering properties are only sensitive to the ice habit. (Cole et al., 2014) used the polarization measurements to infer the ice cloud habit and surface texture (i.e., surface roughness) by minimizing the difference between the observed polarized reflectivities from POLDER and the simulations. Based on the methodology used in Labonnote et al. (2001), Fig. 18 shows the behavior of different ice habits on modeling the polarization properties of ice clouds. The color contours in the figure represent the frequency of occurrence of observed polarized reflectivities from the POLDER 0.86 μm band for ice cloud pixels on 1 August 2007. Superimposed on the color contours are 3000 black dots in each panel that indicate the modeled values for a randomly chosen set of solar-viewing geometries provided in the POLDER data. The scattering properties of four ice habits (i.e. sphere, smooth-faced hexagonal solid column, severely roughened hexagonal solid column, and severely roughened aggregate of solid columns) are used, and the same effective diameter of 50 μm is assumed for the four models. For the two smooth (i.e., non-roughened) models, the simulations show significant oscillations over the entire range of the scattering angles (especially for the spherical model), which are quite different from the observations. When the severely roughened ice habits are used, the modeled polarized reflectivity becomes a relatively smooth function of the scattering angle, and tends to agree more closely with the observations. However, neither the roughened column nor the aggregate model can match the observations over the entire scattering angle range from 60^° to 180^°. The results in Fig. 18 demonstrate a need for improvement with regards to the ability of the models to represent the polarization properties of ice clouds from global satellite observations. For example, the two-habit ice cloud microphysical property (habit/shape) model developed by (Liu et al., 2014) overcomes the aforementioned problem in the polarized radiative transfer simulations involving ice clouds. Furthermore, accounting for the polarization state is important in the forward light scattering and radiative transfer simulations (i.e., the generation of LUTs) involved in cloud property retrievals using passive radiometric measurements (Yi et al., 2014).

Figure 18.  Comparison of normalized polarized reflectivities obtained from one day of PARASOL data over ocean in the color contour to numerical results based on different ice crystal habit models: (a) sphere; (b) smooth hexagonal column; (c) roughened hexagonal column; and (d) roughened aggregate.

Fast, efficient, and accurate parameterizations for ice cloud microphysical and optical properties are employed in broadband radiative transfer models, generally for use in numerical weather prediction (NWP) and general circulation models (GCMs). Advances in the development of efficient ice cloud parameterizations are possible because of the improvement in the light scattering calculations that account for complexities in particle size and shape in a wide spectral region (e.g., from ultraviolet to far-IR), as well as the increasing knowledge about ice cloud microphysical properties (i.e., ice particle size and habit distributions) derived from in situ observations. As noted in the introduction, broadband parameterization progress tends to occur sporadically rather than from a concerted effort.

To derive a parameterization for the ice cloud bulk optical properties, some prognostic variables in the model must be related. For example, (Platt and Harshvardhan, 1988) developed a parameterization in which the cirrus optical properties depend solely on ice cloud temperature, based on some observational evidence that ice water content and ice particle size increase with temperature. But such a parameterization apparently fails to reflect the ice particle size dependence of the optical properties on ice particle size.

At one time, GCMs treated ice clouds as being composed of spherical ice particles. For example, HadAM3 (Hadley Centre Coupled Model, version 3) (Pope et al., 2000) assumed a spherical ice particle model with an effective radius of 30 μm to represent ice clouds. However, the use of ice spheres introduces significant errors in both the flux and heating/cooling rate simulations, thus causing errors in assessing the ice cloud radiative forcing.

After (Slingo, 1989) established a new GCM parameterization for water cloud radiative properties, (Ebert and Curry, 1992) adopted the same idea and first proposed the use of ice water content and the ice particle effective radius in parameterizing the shortwave ice cloud mass extinction coefficient, single-scattering albedo, asymmetry factor, and longwave ice cloud mass absorption coefficients. The parameterization was used in the National Center for Atmospheric Research (NCAR) Community Climate Model, Version 3 (CCM3) (Kiehl et al., 1998). (Ebert and Curry, 1992) used the ice scattering properties of randomly oriented hexagonal particles (Takano and Liou, 1989a). The parameterization scheme developed by (Fu, 1996) and the subsequent updates (Fu, 2007; Fu et al., 1999) successfully incorporated the scattering properties of hexagonal ice crystals and the delta-four-stream approximation into the Fu-Liou radiative transfer model (Fu, 1996). With a similar methodology, an improved parameterization scheme for the ice cloud bulk radiative properties was implemented in the UCLA (University of California, Los Angles) general circulation model (Gu et al., 2003) and the Weather and Research Forecast (WRF) model (Gu et al., 2011). (Key et al., 2002) and (Chou et al., 2002) parameterized the ice cloud optical properties assuming either different individual ice particle habit or habit mixtures, and found large differences in the simulated ice cloud fluxes. (Edwards et al., 2007) adopted an ice aggregate model (Baran and Francis, 2004) that additionally considered the roughened surface textures of ice crystals for the parameterization in the Met Office's GCM. Parameterizations using other particle shape assumptions (i.e., Chebyshev particles) also exhibit their applicability in different modeling studies (McFarquhar et al., 2002; 2003). Note, substantial discrepancies may exist in radiative flux calculations based on different parameterization schemes for the radiative properties of ice clouds (e.g., Zhang et al., 2015).

Furthermore, it is worth noting that Mitchell et al. (1996, 2006) developed a parameterization for the NCAR CESM (Community Earth System Model) based, in part, on the modified anomalous diffraction approximation (MADA) and the database of (Yang et al., 2000). The MADA method provides approximate solutions for the extinction and absorption efficiencies of ice crystals, but does not provide solutions for the scattering phase matrix or the asymmetry factor. Thus, the MADA method must be combined with other more rigorous methods to provide complete single-scattering properties. Note, inconsistencies in the optical properties calculated with different methods cause additional uncertainties in the parameterization.

But how much complexity in the ice habits should be built into a parameterization that will be used for global radiance calculations? Some research teams prefer the use of a single habit, with the observation that the actual microphysical properties of the habit are not as important for this application as the scattering properties. However, a variety of ice habits are found in field campaign observations, and different habits have been shown to have quite distinctive scattering properties. Additionally, finding any individual habit whose polarization properties match closely with POLDER/PARASOL measurements is difficult. The adoption of a single habit or a mixture of habits to develop a broadband property parameterization is currently a significant topic of discussion in the research community.

(Yang et al., 2005) studied and parameterized the longwave scattering properties of ice clouds using the ice habit mixture model employed in the MODIS Collection 4 satellite cloud retrievals. (Hong et al., 2009a) revisited and improved the ice cloud optical property parameterization for the shortwave and longwave spectrum using an updated ice single-scattering property database with 6 ice particle habits, but also incorporated new information regarding ice particle habit and size distributions. In this study, ice cloud radiative forcing simulations displayed large differences in the shortwave for optically thick clouds and in the longwave for small effective diameters.

Except for the parameterization scheme developed by (Edwards et al., 2007), ice particles were assumed to have smooth surfaces in the studies described previously, although increasing observational evidence indicates roughened particles (e.g., Shcherbakov et al., 2006; Ulanowski et al., 2006; Nousiainen and Muinonen, 2007; Neshyba et al., 2013; Ulanowski et al., 2014). With the update to a new ice scattering property database (Yang et al., 2013), which contains nine ice habits and three degrees of surface roughening with the particle size range of 2-10 000 μm and the spectral wavelength range of 0.2-100 μm, (Yi et al., 2013) investigated the influence of ice particle surface roughening on the ice cloud radiative effects simulated by the Rapid Radiative Transfer Model GCM version (RRTMG) (Iacono et al., 2008) and the Fu-Liou single-column radiative transfer models. These broadband simulations were performed assuming either completely smooth or severely roughened ice cloud models with a general habit mixture (Baum et al., 2014), and used ice cloud particle size distributions (over 14 000) from 11 field campaigns.

Based on the same principles as the Fu scheme (Fu and Liou, 1993; Fu, 1996; Fu et al., 1999; Fu, 2007), (van Diedenhoven et al., 2014) presented a flexible parameterization to accommodate the derivation of ice cloud optical properties with any combination of volume, projected area, aspect ratio, and distortion parameter at any shortwave wavelength, unlike the conventional schemes which were restricted to applications in certain preselected spectral bands. However, the parameterization used the CGOM rather than the IGOM, and thus the calculations did not account for some advances in light scattering reported earlier in this article. A separate study concerning the problem was reported in (Bi and Yang, 2014a), and focused on the assessment of the CGOM's applicability to remote sensing and radiation parameterization involving ice clouds.

(Baran, 2012) discussed and pointed out the deficiencies of the traditional parameterization techniques based on a relationship with the effective particle size. (Baran et al., 2014) developed a new, coupled parameterization based on a weighted ice habit mixture model, with the bulk optical properties parameterized as functions of wavelength and ice water content (IWC). Basically, the approach couples the ice cloud bulk optical properties to the IWC, and has been implemented in the Met Office Unified Model Global Atmosphere 5.0 (GA5). The parameterization was tested in a 20-year simulation, with the result that a tropical troposphere cold temperature bias was reduced, while the southern midlatitudes were warmed by approximately 0.5 K. Furthermore, we wish to point out that ice clouds induced by human activities, such as contrails, have a non-negligible impact on the radiation budget, particularly in air-traffic corridors. Contrails are generally treated as natural cirrus clouds containing primarily small particles. (Xie et al., 2012) employed the depolarization ratio derived from space-borne lidar measurements to develop an appropriate contrail particle model consisting of a variety of ice habits. A contrail optical property parameterization was also developed and used in the RRTMG and tested using the CAM5 atmospheric general circulation model (AGCM) to estimate the regional and global contrail radiative forcing (Yi et al., 2012).

The differences in ice cloud radiative effects based on the spherical and nonspherical ice cloud particle habit models are quantified with the NCAR CAM5 AGCM. Two sets of parameterization schemes are employed: one developed using the Lorenz-Mie theory assuming spherical ice cloud particles, and the other developed using II-TM+IGOM (Invariant Imbedding T-matrix method and the Improved Geometric-Optics Method) (Yang and Liou, 1996b; Bi and Yang, 2014a) assuming ice clouds consisting of solid hexagonal column ice aggregates with severely roughened surfaces. Both the cases use particle size distributions that follow the modified Gamma size distribution with a mean effective variance of 0.1. Following the formulae given in (Yi et al., 2013), the ice cloud optical properties, i.e., mass extinction and absorption coefficients, single-scattering albedo, and asymmetry factor, are parameterized as functions of the effective diameter for each of the shortwave and longwave radiation bands of the RRTMG. The fitting equations for the parameterized bulk radiative properties are: \begin{eqnarray} \left\langle\dfrac{\beta_{\rm e}}{{\rm IWC}}\right\rangle&=&\dfrac{a_0+a_1D_{\rm eff}^{-1}}{1+a_2 D_{\rm eff}^{-1}+a_3D_{\rm eff}^{-2}} ,\\ \left\langle\dfrac{\beta_{\rm a}}{{\rm IWC}}\right\rangle&=&\dfrac{b_0+b_1D_{\rm eff}^{-1}}{1+b_2 D_{\rm eff}^{-1}+b_3D_{\rm eff}^{-2}} ,\\ 1-\langle\varpi\rangle&=&\sum_{n=0}^{n=4}c_{\rm n}D_{\rm eff}^n ,\\ \langle g\rangle&=&\sum_{n=0}^{n=4}d_{\rm n}D_{\rm eff}^n , \end{eqnarray} where β _e/ IWC and β _a/ IWC are the mass extinction and mass absorption coefficients; \(\varpi\) is the single-scattering albedo; g is the asymmetry parameter; and, a, b, c, and d are the fitting coefficients.

The two parameterization schemes are respectively implemented in the NCAR Community Atmosphere Model version 5 (which uses RRTMG) for a 10-year modeled climatology of the cloud radiative effect. The AGCM model horizontal resolution is 1.9^° (lat) × 2.5^° (lon), and has 30 vertical layers. Figure 19 shows the comparison of bulk ice cloud optical properties calculated with the II-TM and the Lorenz-Mie theory for selected bands. Little impact is found in the mass extinction coefficient and single-scattering albedo, but apparent differences exist in the asymmetry factor. The differences in the asymmetry factor are the most significant for shorter wavelength bands and smaller effective particle diameters, and can range from 0.05-0.1.

Figure 19.  Ice cloud optical properties parameterized as functions of effective diameter based on the II-TM+IGOM (upper panels) and the Lorenz-Mie theory (lower panels) calculations for selected bands of RRTMG.

Figure 20 shows the differences (Lorenz-Mie theory case minus II-TM+IGOM case) in the shortwave, longwave, and net (shortwave + longwave) cloud radiative effects (CREs) at the top of the atmosphere (TOA) between the cases with II-TM+IGOM and Lorenz-Mie theory optical property parameterizations. The shortwave CRE differences are mostly positive and are larger than the longwave CRE differences. The results indicate that the shortwave cloud radiative effect is underestimated, whereas the longwave cloud radiative effect is overestimated. The global averaged CRE difference is 2.41 W m^-2 for the shortwave and 1.13 W m^-2 for the longwave. Thus, the net CRE difference can be up to 3.54 W m^-2 globally, and ∼ 18 W m^-2 regionally (i.e., in the tropical warm pool).

Figure 20.  The differences in simulated cloud radiative effect between the cases with the Lorenz-Mie theory and the II-TM+IGOM ice optical properties: (a) shortwave cloud radiative effect (SWCRE) difference; (b) longwave cloud radiative effect (LWCRE) difference; and (c) net cloud radiative effect (NTCRE) difference. Units: W m^-2.

Figure 21 shows the simulated difference in the surface temperature using the ice cloud optical properties from the Lorenz-Mie theory and II-TM+IGOM. Given the facts that the shortwave TOA CRE of the Lorenz-Mie theory case is weaker and the longwave CRE is stronger, it is not surprising that the Lorenz-Mie theory case simulates a warmer climate. The warmer regions mostly exist in the middle and high latitudes over land.

Figure 21.  (a) 10-year averaged surface temperature simulated with the II-TM+IGOM ice optical properties. (b) Difference in simulated surface temperature between the cases with the Lorenz-Mie theory and the II-TM+IGOM ice optical properties. Units: K.

The simulated shortwave and longwave heating rates in the vertical direction are also affected. Figure 22 shows the zonal average of shortwave and longwave heating rates during the 10-year simulation, as well as the corresponding differences between the Lorenz-Mie theory case and II-TM+IGOM case. The Lorenz-Mie theory case contributes to both a weaker heating rate for the shortwave and a weaker cooling rate for the longwave than the II-TM+IGOM case. This study illustrates well the importance of the improvements to the ice cloud optical properties for better simulating the global cloud radiative forcing and atmospheric radiation budget.

Figure 22.  10-year zonally averaged vertical sections of (a) shortwave heating rate, (b) shortwave heating rate difference between the Lorenz-Mie theory and the II-TM+IGOM cases, (c) longwave heating rate, and (d) longwave heating rate difference between the Lorenz-Mie theory and the II-TM+IGOM cases. Units: K d^-1.

Electromagnetic scattering by dielectric particles is a mature branch of physics. Although the scattering of light by homogeneous spheres was formulated (commonly known as the Lorenz-Mie or Mie theory) more than 100 years ago, the single-scattering properties of more complex nonspherical and inhomogeneous particles remain a significant challenge in applied and computational physics. Within the past 30 years, incremental progress has been made to develop numerical methods and approximations to derive the single-scattering properties of nonspherical particles. In this paper, we review the application of the geometric-optics, FDTD/PSTD, DDA, and T-matrix methods to the computation of the extinction efficiency, single-scattering albedo, and full phase matrix of various nonspherical ice crystals. With the current state of computational capability, the numerically rigorous methods, in practice, are still confined to a limited range of size parameters. The current approach for calculating ice particle optical properties, by combining numerically exact methods with geometric-optics approximations, is unlikely to change in the near future. To bridge the gap between rigorous solutions and geometric-optics approximations, there is a pressing need to develop computationally efficient capabilities for the semi-classical region, where the diffraction effects, including edge effects and surface waves, cannot be taken into account in the geometric-optics approximations from first principles.

We also review various ice cloud retrieval algorithms developed for satellite sensor measurements of ice clouds at different wavelengths and present some case studies using the MODIS and POLDER sensors to illustrate the current capabilities of inferring ice cloud properties. The ice cloud macrophysical and microphysical properties inferred from satellite remote sensing include cloud top pressure/height, optical thickness, particle effective diameter, and particle habit. Many research teams use a bi-spectral method that employs two solar bands (one non-absorbing and one moderately absorbing band) to infer cloud optical thickness and particle effective diameter. The retrievals are limited to only daytime measurements, and the related uncertainties increase at high solar zenith angles, at low optical thicknesses, and over bright surfaces. However, in simulating the single-scattering properties, uncertainties tend to decrease if the ice crystal surfaces are assumed to be rough. Furthermore, the comparison between the retrieved results based on spherical and hexagonal models indicates that the algorithm is very sensitive to the scattering properties of the assumed ice particle habit model.

In comparison, the IR-band-based retrievals, which are based on the IR emissions directly from the atmosphere and the cloud, can be applied to both daytime and nighttime observations. IR-based retrievals are best for ice clouds of lower optical thickness, making this approach somewhat complementary to solar-based methods. However, there is some sensitivity to surface emissivity, while there is little sensitivity to the assumed ice habit or particle roughness. A case study based on three MODIS IR bands confirms that the IR-band-based retrieval is insensitive to the assumption of ice cloud habits. Note that the IR-band-based retrievals are insensitive to optically thick clouds, and are only applicable to thin ice clouds.

A comparison between the solar-band and IR-band retrievals indicates some inconsistencies in the inferred optical thickness between the two approaches, suggesting the necessity for further improvement of the ice cloud models. We also illustrate the unique capabilities of using polarization observations to infer ice cloud properties. Most current retrieval techniques have limitations that reduce the consistency of the results retrieved from different algorithms and sensors. With the increasing number of sensors flown on aircraft and in space, and considering the sensors under development, more accurate and robust retrieval algorithms and more representative ice cloud models are greatly needed to better infer the ice cloud microphysical properties. Consistency may be achieved through combinations of multiple spectral channels, multiple instruments, and polarimetric observations.

From the review of the parameterizations of ice cloud optical properties proposed over the past few decades, we can see that improvements are driven by the continuous advances in the capabilities of calculating light-scattering properties of ice particles with a wide range of sizes, complex habits, and various surface roughness features, as well as the increasing observational evidence made possible by aircraft in-situ measurements. Early parameterizations were based on rather limited theoretical calculations, and show the limitations and drawbacks of assuming spherical ice particles or solely hexagonal column models. The early parameterizations assumed dependence of ice cloud optical properties on a single variable (particularly temperature). At this time, the parameterization scheme assuming a randomly oriented single hexagonal ice particle with optical properties parameterized as functions of ice water content and effective particle radius is the most widely used approach. There is, at the time of writing, some controversy over whether to adopt a single habit or a habit mixture to represent ice clouds for retrievals and parameterizations. Other complexities, such as the ice cloud particle surface roughness features, have been found to be critical for radiative transfer simulations of ice clouds. The importance of better representing ice cloud optical properties in modeling studies is illustrated through a comparison between two 10-year GCM simulations using spherical and nonspherical ice particle assumptions.

Although new ice cloud optical property parameterizations have been proposed, some open questions are still to be addressed. First, a more suitable variable is needed to establish the parameterization. The effective radius/diameter is commonly used, but is not a direct output of the cloud microphysical modules of many GCMs, while ice water content may be a more appropriate surrogate. Second, ice habits that are physically constrained by observations should be used to facilitate the validation and comparison of modeling results with observations. We suggest that a worthwhile goal is to work towards a uniform ice cloud model that would be consistent for both satellite remote sensing implementations and weather/climate modeling applications.