A two-habit model for the microphysical and optical properties of ice clouds

To provide a better representation of natural ice clouds, a novel ice cloud model is developed by assuming an ice cloud to consist of an ensemble of hexagonal columns and 20-element aggregates with specific habit fractions at each particle size bin. The microphysical and optical properties of this two-habit model (THM) are compared with both laboratory and in situ measurements, and its performance in downstream satellite remote sensing applications is assessed. The ice water contents and median mass diameters calculated based on the THM closely agree with in situ measurements made during 11 field campaigns. In this study, the scattering, absorption, and polarization properties of ice crystals are calculated with a combination of the invariant imbedding T matrix, pseudo-spectral time domain, and improved geometricoptics methods over an entire practical range of particle sizes. The phase functions, calculated based on the THM, show close agreement with counterparts from laboratory and in situ measurements and from satellite-based retrievals. When the THM is applied to the retrievals of cloud microphysical and optical properties from MODIS (the Moderate Resolution Imaging Spectroradiometer) observations, excellent spectral consistency is achieved; specifically, the retrieved cloud optical thicknesses based on the visible/near infrared bands and the thermal infrared bands agree quite well. Furthermore, a comparison between the polarized reflectivities observed by the PARASOL satellite and from theoretical simulations illustrates that the THM can be used to represent ice cloud polarization properties.


Introduction
Ice clouds, i.e., high clouds containing ice crystals with various sizes and shapes, on average cover over 20 % of the Earth and up to 60-70 % of the tropical areas (Lynch et al., 2002;Nazaryan et al., 2008;Baran, 2009).Not surprisingly, ice clouds significantly influence both the climate system radiation budget and large-scale circulations in the atmosphere (Herman et al., 1980;Liou, 1986;Minnis et al., 1993a, b;Sassen and Comstock, 2001;Stephens et al., 1990;Stephens, 2005;Ramanathan et al., 2007;Loeb et al., 2009).However, owing to uncertainties in the ice microphysical properties (particle habit and size distribution) and optical properties (extinction coefficient, single-scattering albedo, and scattering phase matrix), ice clouds are still one of the least understood atmospheric components from the perspective of remote sensing and radiative transfer simulations involved in general circulation models (GCMs).Thus, a realistic and robust ice cloud model is being sought and of vital importance to atmospheric research.
A numerical model of ice clouds normally assumes either a single particle habit (i.e., particle shape) or an ensemble of habits (Baran and Labonnote, 2007;Baran et al., 2009;Baran, 2009Baran, , 2012;;Yang et al., 2013;Baum et al., 2014).The optical properties determined based on the particle habit/habits are fundamental to the downstream applications in remote sensing, radiative transfer and GCMs (Ebert and Curry, 1992;Fu and Liou, 1993;Fu, 1996Fu, , 2007;;Minnis et al., 1998;Katagiri et al., 2010;Heymsfield and Miloshevich, 2003;Baum et al., 2011;Edwards et al., 2007;Yi et al., 2013;Baran et al., 2014a).Thus, in order to reduce the uncertainties in downstream applications, an improved repre-C.Liu et al.: A two-habit model for the microphysical and optical properties of ice clouds sentation of ice cloud particle habits and optical properties is needed for the construction of a robust ice cloud model.Numerous laboratory and in situ measurements have been made to improve our knowledge of ice clouds, and various satellite observations have also played important roles in determining their microphysical and optical properties (Minnis et al., 1993b;Heymsfield and Miloshevich, 1995;Gayet et al., 2006Gayet et al., , 2012;;Lawson et al., 2008;Febvre et al., 2009;Heymsfield et al., 2013).The observations from different perspectives serve as the most practical and insightful standards from which to develop an ice cloud model.This study considers the currently available data in an attempt to improve the representation of ice clouds with a theoretical model based on two particle geometries.
As one of nature's greatest artworks, ice crystals show a myriad of variations in shape/habit for different meteorological conditions.Ice cloud habit study begins with an understanding of the microphysical processes necessary for nucleation, diffusion growth, collision and aggregation within the atmosphere.Both laboratory and in situ observations have contributed meaningful information about ice crystal shapes (Magono and Lee, 1966;Heymsfield et al., 2002Heymsfield et al., , 2005;;Lawson et al., 2006).The studies indicate that ice crystals occur with geometries with various degrees of complexity, e.g., pristine hexagonal columns, plates, and bullets, rosettes of different forms, and complicated and irregular aggregates.Furthermore, some detailed structures, such as surface roughness, hollow structure, and inhomogeneity (with air bubbles or ice nuclei inside), have been widely noted in the observations and considered for numerical studies (Ulanowski et al., 2012(Ulanowski et al., , 2014;;Schmitt and Heymsfield, 2007;Labonnote et al., 2001).Both the ice particle overall geometry and detailed structure have a significant effect on the optical properties.Thus, constructing an idealized model both geometrically and optically representative of natural particles is quite challenging.
In addition to observations of particle geometries, various measurements have been attempted to obtain reliable information on the microphysical and optical properties of ice clouds (Curry et al., 2000;Heymsfield et al., 2013).A series of field campaigns were conducted at a variety of midlatitude and tropical locations in both hemispheres over the period between 1999 and 2006 to study the microphysical properties of ice clouds (Heymsfield et al., 2013).The microphysical data collected include the particle size distribution (PSD), ice water content (IWC), and median mass diameter (D mm ).The D mm is defined as the diameter at which the mass in the PSD with smaller particles equals that with larger ones.Moreover, the ice cloud optical properties were obtained in numerous laboratories and field campaigns.The polar nephelometer (PN) was widely used to measure the scattering phase function of an ensemble of cloud particles (water droplets, ice crystals, or a mixture of both) (Sassen and Liou, 1979;Gayet et al., 1998Gayet et al., , 2004;;Barkey and Liou, 2001;Auriol et al., 2001;Febvre et al., 2009).Although lim-ited spatially and temporally, the measurements have played an essential role in the numerical studies of ice clouds, and will be fully considered in this study.
Satellite observations are used to infer cloud properties by comparing sensor measurements and radiative transfer simulations for a set of known cloud and atmospheric conditions (Wielicki et al., 1998;Chepfer et al., 2002;Winker et al., 2003;Platnick et al., 2003;Knap et al., 2005;McFarlane and Marchand, 2008;Minnis et al., 2011;Baran et al., 2012b).The satellite measurements may be at either visible/nearinfrared solar bands or thermal infrared (IR) bands, and may also include polarization.Sensors on board satellites flying as part of the NASA Earth Observing System A-Train constellation simultaneously provide measurements encompassing all of these characteristics (L'Ecuyer and Jiang, 2010).The satellite-based retrieval of ice cloud properties, e.g., the effective particle diameter (D eff ) and optical thickness (τ ), relies on the use of accurate and efficient radiative transfer models to simulate the cloud radiances at the top of atmosphere, and the optical properties of a given ice cloud model are required for those simulations (Minnis et al., 1993a(Minnis et al., , b, 2011;;Platnick et al., 2003).However, when applied to satellite remote sensing (e.g., based on the Moderate Resolution Imaging Spectroradiometer (MODIS) observations), most ice cloud models encounter a challenge known as spectral inconsistency, i.e., significant differences occur in cloud optical thicknesses retrieved with different spectral bands (e.g., solar or thermal IR bands) for the same cloud model (Baum et al., 2014).Thus, another goal of this study is to construct an ice cloud model that can infer consistent optical properties in solar and thermal IR-band retrievals.
Using the ice cloud polarization properties (e.g., polarized reflectivity) has increasingly gained attention as a means to infer the microphysical and optical properties (van Diedenhoven et al., 2012Diedenhoven et al., , 2013;;Labonnote et al., 2001), and such applications can be widely found with available observations from the PARASOL (Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar) satellite (Labonnote et al., 2001;Baran and Labonnote, 2007;Baran, 2009;Cole et al., 2013).The radiometer/polarimeter on board POLDER (POLarization and Directionality of the Earth's Reflectances) measures the I , Q, and U components of the Stokes vector at three wavelengths with up to 16 viewing angles for each pixel (Deschamps et al., 1994).Previous studies have indicated that the polarized reflectivities simulated based on several ice cloud models, using either an individual habit or a mixture of multiple habits, can approximately match those of PARASOL observations (Doutriaux-Boucher et al., 2000;Baran, 2009;Cole et al., 2013).The polarized reflectivity data from the satellite have also been used to retrieve the ice particle habit (Sun et al., 2006) and degree of surface roughness (Cole et al., 2014).Thus, as another unique perspective of ice cloud, the polarization properties of a numerical model must be carefully checked.
With the variety of laboratory experiments, field campaigns, and satellite sensors to measure the microphysical and optical properties, constructing an ice cloud model that can consistently represent a wide range of perspectives is extremely challenging.This study strives to develop a robust ice cloud model based on two particle geometries, the twohabit model (THM), and to verify its performance in modeling the microphysical and optical properties of natural ice clouds.Section 2 reviews some of the previous ice cloud models and introduces the novel THM.Section 3 discusses the THM microphysical properties.Section 4 shows the optical properties of the THM, and comparisons with measurements.The THM performance in satellite remote sensing applications is presented in Sect. 5. Section 6 contains the conclusions.

Two-habit model
Ice cloud particle geometries show significant variations for different meteorological conditions, especially temperature and relative humidity (Magono and Lee, 1966;Heymsfield and Miloshevich, 1995).A tremendous amount of effort has been devoted to the detailed study of ice crystal geometries from both laboratory and in situ observations, and to classify ice crystals into multiple categories based on the general geometries.The most widely observed ice crystal types include hexagonal columns, hexagonal plates, bullet rosettes, and aggregates of various pristine particles (Magono and Lee, 1966;Korolev et al., 1999;Heymsfield et al., 2002Heymsfield et al., , 2005;;Evans et al., 2005;Lawson et al., 2006).A number of studies have been reported on building numerical models for ice clouds with idealized geometries, and in using the corresponding microphysical and optical properties to represent natural clouds.For different applications, either the microphysical or the optical properties at certain wavelengths are normally considered in the models, and seldom does an ice cloud model consistently represent all ice cloud properties for different applications (Baran et al., 2014b).
Due to the limitation in numerical simulations of light scattering by non-spherical particles, the single hexagonal column model for ice clouds was introduced to atmospheric applications in the 1970s and 1980s (Wendling et al., 1979;Cai and Liou, 1982;Takano and Liou, 1989a, b), but additional particle habits have since been developed and applied (Macke, 1993;Takano and Liou, 1995;Macke et al., 1996;Yang and Liou, 1998;Um andMcFarquhar, 2007, 2009).Various commonly occurring ice cloud habits are now widely used in radiative transfer and remote sensing, and popular examples include hexagonal columns and plates of various aspect ratios, droxtals, polycrystals, solid or hollow bullet rosettes, and aggregates of columns, plates or rosettes (e.g., Baran, 2009;Yang et al., 2013;and references cited therein).The models use either an individual particle habit or an ensemble of habits.When multiple habits are used, the habit fractions normally vary for different particle sizes.The optical database and parameterization based on the numerical models are normally made for further applications in remote sensing, radiative transfer and GCMs (Fu, 1996(Fu, , 2007;;Minnis et al., 1998;Edwards et al., 2007;Letu et al., 2012;van Diedenhoven et al., 2014).Ice cloud models with a single habit have been found useful for some special applications as the computational burden for single-scattering simulations is minimized.However, the single habit models are limited in the aspect of consistently representing multiple ice cloud properties (Baum et al., 2014).Natural ice clouds show unclear variation of particle habits and have different habit preferences at different size ranges.Multiple ice habits definitely introduce much more freedom to accurately represent the microphysical and optical properties.However, without explicit theoretical or observational data, the choice of the habits and habit fractions is arbitrary.Furthermore, the accurate calculation of the scattering properties of different non-spherical habits is very time consuming.Thus, we will attempt to construct an ice cloud model that can capture and represent all major properties by using as few particle habits as possible, which will simplify the model and minimize the computational burden for computing the single-scattering properties.
A study by Schmitt and Heymsfield (2014) suggests that atmospheric ice particles can be separated into two categories in terms of particle complexity (i.e., simple and complex) by using particle imagery data from high-resolution aircraft particle imaging probes.A dimensionless parameter representing the particle "complexity" is defined based on particle projected area, area ratio and perimeter, and a cutoff value is chosen to identify pristine and complex ice particles imaged by the Cloud Particle Imager (CPI) probe with a resolution of a few microns.The results of two example data sets from in situ measurements indicate that complex particle habit fraction increases as the particle maximum dimension increases.The idea of separating ice crystals into two general categories and the conclusions obtained from the study of Schmitt and Heymsfield (2014) are of great importance in the simplification of ice cloud models.Moreover, numerical studies indicate that the optical properties of particles of the same kind are not strongly affected by the number and orientations of monomers of, e.g., complex aggregates of hexagonal plates (Xie et al., 2011) or bullet rosettes (Um and McFarquhar, 2007), particularly if the monomers are sufficiently separated, so that the multiple scattering among the monomers is negligible.Therefore, as a quite accurate approximation, it is possible to use a relatively simple particle morphology to represent a group of more complicated counterparts in the computation of the particle optical properties.
Based on the preceding physical rationale and the observations and classifications given by Schmitt and Heymsfield (2014), in addition to the consideration of the computational burden, this study explores the feasibility of using a simple habit and a complex habit to represent ice clouds.To construct the model, we need to determine the most rep- L is equal to the maximum dimension D for the hexagonal column.resentative particle habits.A hexagonal column, considered as the simple/pristine particle, is the primary candidate, and a type of complex aggregate is the second most widely observed and studied particle habit.We use hexagonal columns as the aggregate monomers.By changing the geometric parameters related to the aspect ratio, number of the monomers, and aggregation configuration, the optical properties of the particles are optimized to match those of natural ice clouds.The aspect ratio of a hexagonal column is defined as 2a / L, where a is the semi-width of the hexagonal cross section and L is the column length.A hexagonal column with an aspect ratio equal to one, and an aggregate with 20 hexagonal column monomers are used in this study.The relative sizes and aspect ratios of the 20 monomers are randomly generated, and they are point attached to form the aggregate.The particle maximum dimension D is used to specify the particle size (i.e., L for the hexagonal column, and the maximum distance between two points on the particle for the aggregate), and the size parameter, which is an important parameter for light scattering simulations, is defined as π D / λ.
In addition to the overall geometries of the two habits, the detailed structures of natural crystals are considered.In situ measurements have indicated that ice crystals have predominantly hollow structures (Walden et al., 2003;Schmitt and Heymsfield, 2007) and irregular geometries, and, for consideration of these facts in the THM, the hexagonal monomer of the aggregate is assumed to have hollow structures similar to those used by Yang et al. (2013).Figure 1 illustrates the geometry of a hollow hexagonal element.The depth of the hollow structure is specified by d and d/L = 0.25 in this study.From observations, particle surface roughness is widely noted as an important ice crystal feature (Cross, 1969;Ulanowski et al., 2006Ulanowski et al., , 2012Ulanowski et al., , 2014;;Neshyba et al., 2013), and numerical studies indicate that the surface roughness has significant influence on the particle optical properties and cloud radiative effect (Yi et al., 2013), especially the angulardependent scattering phase matrix elements (Peltoniemi et al., 1989;Macke et al., 1996;Shcherbakov et al., 2006;Yang et al., 2008).The hexagonal columns and the aggregates over the entire size range considered will be treated as roughened particles with the same degree of surface roughness.The technical details of the roughened surface definition can be found in Liu et al. (2013).In this study, severely roughened particles (Yang et al., 2013) are used.
Figure 2 shows the two particle geometries used for the THM, and both the hollow structure and surface roughness are illustrated in the figure.The column is clearly a single but "compact" particle, whereas the aggregate is very complex and loose in the space.The two habits represent the simple and complex ice crystals classified by Schmitt and Heymsfield (2014).The processes necessary to form the hexagonal aggregate and its geometric parameters are detailed in the Appendix A.

Microphysical properties
With the explicit geometries of the two particle habits defined, the habit fraction, as a function of particle maximum dimension, becomes a key parameter to determine the microphysical properties of the THM.This section introduces the habit fraction used for the THM and compares the simulated microphysical properties, i.e., IWC and D mm , with those from in situ measurements.
As discussed in Sect.2, Schmitt and Heymsfield (2014) separate ice crystals into simple and complex categories by analyzing CPI images, and show that the complex habit fraction increases as the maximum dimension increases.In the THM, the hexagonal column and aggregate correspond, respectively, to the simple and complex particles.Although that the exact percentages of the simple and complex particles may differ from case to case, the first role in qualitatively determining the habit fractions is to increase the aggregate fraction as the particle maximum dimension increases, which ensures that the geometric model represents natural crystals.
A more quantitative way to determine the habit fraction is to consider the microphysical data sets from the in situ measurements.The IWC and D mm of ice clouds are highly related to the ice particle volume, which is a strong function of the particle maximum dimension.With a given PSD and the two particle habits in the THM, the IWC and D mm are determined by and where ρ ice is the density of solid ice (a value of 0.917 g cm −3 is used in this study), D min and D max are the minimum and maximum particle sizes in the distribution, and n(D) is the number concentration of particles with a maximum dimension of D. f c (D) and f a (D) are the habit fractions of the column and aggregate in the THM, and, for any size, and V a (D) indicate particle volume for the two particle habits.We use the microphysical data collected from 11 field campaigns, and a detailed summary of these data can be found in Heymsfield et al. (2013) and Baum et al. (2014).Coefficients for the gamma size distribution are fitted to the data sets of the particle number concentration versus size and provided for each individual PSD, and a total of over 14 000 PSDs with cloud temperatures colder than −40 • C are used to ensure the measured clouds are indeed ice clouds.With certain habit fractions, the IWC and D mm based on the THM can be computed for each PSD by the integrals given by Eqs. ( 1) and ( 2), and may then be compared with the observations.Fitted gamma size distributions from the data are used for the aforementioned integrals.After we tested different habit fractions to minimize the differences between the simulated and observed IWC and D mm , we chose a continuous habit fraction for the hexagonal column that leads to close agreement of the microphysical properties, which is given by the following: and the fraction of aggregate is given by f Figure 3 shows the THM habit fractions obeying Eq. ( 3).The fraction of the aggregate, i.e., the complex particle, smoothly increases with increasing particle diameter, and the trend is the same as that obtained from ice crystal image analysis.For small particles with maximum dimensions less than 100 µm, we assume over 80 % of the ice crystals to be hexagonal columns, but the fraction drops to only 1.7 % for particles larger than 1500 µm.
Note that, considering the uncertainties in the observations, the final habit faction we use for the THM does not necessarily give the best fit to all in situ data, but we find that all habit fractions with similar trends lead to similar agreement in the microphysical properties.Furthermore, considering the significant variation of ice clouds under different meteorological conditions, no single "best" exists for all ice clouds, because the best for one condition may not represent the ice cloud properties under another condition.For applications of ice clouds having very different microphysical properties, the fractions of the two habits can be easily modified to match the specific properties.Thus, the continuous habit fraction given by Eq. ( 3) is used in the THM and the following simulations.
With the habit fractions given, the upper panels of Fig. 4 compare the measured and calculated IWC and D mm values for each of the PSDs from the 11 field campaigns.The names of the field campaigns are listed in the figure and differentiated by both colors and symbols.The values, for both IWC and D mm , calculated with the THM are in close agreement with the observations.Slight differences are noticed for D mm at values larger than 500 µm.The largest differences in the IWC are shown for data from the CRYSTAL-FACE campaign.Overall, Fig. 4 indicates that, with two particles and the corresponding habit fractions given by Eq. ( 3), the THM can reasonably represent the microphysical properties of ice clouds.The lower panels of Fig. 4 show the histograms of the distributions of the measured and calculated IWC and D mm .As expected, the THM-based distributions are essentially the same as the measured counterparts.
The relative differences (RDs) between the theoretical microphysical properties based on the THM and the in situ measurements for each of the 11 field campaigns are listed in Table 1.The names of the field campaigns and the related numbers of PSDs with temperatures colder than −40 • C are given.For both D mm and IWC, the means and standard deviations (SDs) of the RDs are also listed.We can see that the mean  RDs for D mm are generally less than 5 %.The only exception is the case of the Stratospheric-Climate Links with Emphasis On the Upper Troposphere and Lower Stratosphere (SCOUT) field campaign with an average RD of 17 % because relatively small D mm (less than 50 µm) values were observed during this field campaign and the measurements are less reliable with such small particle sizes.Averaged for all the field campaign data, the model shows a mean RD of −0.27 % with a standard deviation of 5.2 % for D mm .The RDs for IWC are almost one order larger in magnitude compared with the case of D mm , because the values of measured IWC span more than 6 orders of magnitude.Furthermore, the mean RDs can be as large as 148 % with a standard deviation of 50 % for the CRYSTAL-FACE campaign, although the model works well in the case of data obtained during other campaigns such as the TRMM, ARM-IOP and MPACE.Overall, the present model overestimates the IWC by approximately 13 % with a standard deviation of 24 %.
To further quantify the performance of the THM for modeling the microphysical properties of ice clouds, Fig. 5  els of Fig. 5).We call special attention to the fact that there are significant uncertainties related to the measurements of small particles, and the RDs at small D mm bins show quite large standard deviations.The performance of the model is also sensitive to the IWC.Both the RDs and the SDs for D mm and IWC tend to decrease as IWC increases; particularly when IWC is larger than 10 −2 g m −3 .
Furthermore, the relationship between particle volume (V ) and the particle maximum dimension (D) determines IWC and D mm for a given PSD.The V -D relationship based on the THM is given by 100 µm ≤ D < 1500 µm 0.036D 3 D ≥ 1500 µm. (4) In the above expression, D is specified in units of µm, and V is in units of µm 3 .Figure 6 illustrates the V -D relationship given by Eq. (3).Given the large amount of in situ measurements we used in this study and the quite reasonable agreement between the model results and measurements, the preceding V -D relationship can be used as a reasonably accurate expression to estimate the variation of ice crystal volume as a function of particle maximum dimension for other relevant applications.

Optical properties
With the geometrical and microphysical model (i.e., the two particle habits as well as their habit fractions) discussed above, we turn to the optical properties of the THM.First, we give a brief introduction of the numerical algorithms used to obtain the optical properties, and illustrate the singlescattering properties of the THM.The second subsection compares the modeled phase functions with the results from both measurements and satellite retrievals.
To better illustrate the advantages of the THM, here we also consider a single hexagonal column model for comparison.This single column model (SCM) is based on a smooth surface, and the aspect ratio decreases as the particle size increases.The details of the single column model as well as its microphysical and optical properties can be found in Yang et al. (2013) and Bi et al. (2014).It should be noticed that the SCM we used for this study is based on pristine particles with smooth surfaces, and the conclusions obtained with the present SCM should not be generalized to other single column models.Furthermore, models based on single columns or plates are still widely used for radiative flux calculation and remote sensing implementations (e.g., Fu, 2007;van Deidenhoven et al., 2014), which are articulated to be rational with demonstrated success for some specific applications.

Single-scattering simulations
Numerical simulations of light scattering by randomly oriented non-spherical particles are a major challenge limiting the development of ice cloud models.The conventional geometric-optics method (CGOM), which is relatively simple and computationally efficient, is one of the most popular methods for the solution of light scattering by ice crystals (Cai and Liou, 1982;Takano and Liou, 1989a;Macke et al., 1996;Baran, 2009), although its accuracy in the cases of small and moderate size parameters is questionable due to the inherent shortcomings of the ray-tracing technique.Bi et al. (2014) elaborate on the uncertainties with the CGOM in remote sensing applications and radiative transfer simulations by comparing with results from a benchmark scattering data set obtained with a combination of the Invariant Imbedding T matrix method (II-TM) (Bi and Yang, 2014) and the improved geometric-optics method (IGOM) (Yang and Liou, 1996;Bi et al., 2009).The study indicates that the CGOM errors in inferring the optical thickness and effective diameters from the MODIS observations can be up to 20 %, and on the order of 10 Wm −2 in ice cloud radiative forcing calculations.
The single-scattering properties of ice crystals given by the II-TM (Bi and Yang, 2014) can be considered as a benchmark, because the II-TM solves Maxwell's equations from first principles.Note, the II-TM is applicable to moderately large size parameters for which the IGOM has reasonable accuracy.However, due to the loose structure of the hexagonal aggregate considered in the THM, the computational memory used by the II-TM simulations increases significantly as the particle size increases.To optimize numerical computations, in this study the pseudo-spectral time domain method (PSTD) that is a numerically accurate technique (Liu, 1997;Liu et al., 2012a, b) is employed for the size parameters in the regime between the II-TM and IGOM simulations.The applicability and accuracy of these three methods have been extensively studied in previous studies; they are thus not repeated here.Without discussing technical details, we use a synergic combination of those three numerical models to minimize the bias introduced by light scattering simulations, and the single-scattering properties of the two-particle habits with maximum dimentions from 2 to 10 000 µm at interested wavelengths are simulated.Furthermore, the scattering properties involved in this study are associated with particles with random orientations.
Figure 7 shows the THM and the SCM extinction efficiencies, single-scattering albedos and asymmetry factors as functions of the particle maximum dimension.The singlescattering properties at three wavelengths, 0.67, 2.13 and 12.0 µm, are illustrated.The SCM data are obtained from a combination of the II-TM and IGOM as shown by Bi et al. (2014).The II-TM, PSTD and IGOM are used to cover the entire practical size range that we consider for the THM.With the edge effect included in the results following the approach used in Yang et al. (2013), we can see that smooth curves are obtained for the extinction efficiency and the single-scattering albedo.The SCM and the THM show quite similar patterns for both the extinction efficiency and singlescattering albedo, whereas differences are evident in their asymmetry factors.At visible wavelengths, i.e., 0.67 µm, the THM exhibits an almost constant asymmetry factor with a value of approximately 0.76, whereas the SCM values increase to almost 0.9 as the particle maximum dimension increases.Based on a climatic feedback sensitivity study, Stephens et al. (1990) suggest that a reduction of the Mietheory-based asymmetry factor (∼0.87) to a lower value of 0.7 may be necessary to achieve broad agreement between theory and observation.Thus, reduction of the asymmetry factor from its SCM value (as large as ∼0.9) to the THM (∼0.76) is in alignment with the previous speculation.
For remote sensing applications, the bulk scattering properties of an ensemble of ice particles with specified size distributions are normally used.We assume a Gamma size distribution (Hansen and Travis, 1974) to integrate the bulk scat-   tering properties of the THM.The dimensionless effective variance is assumed to be 0.1, and the effective diameter values increase from 10 to 180 µm in steps of 10 µm (McFarquhar and Heymsfield, 1998).The effective diameter of the particle is defined to be 1.5 × V / A following Foot (1988), where V and A are the volume and projected area of the particles.
Figure 8 shows the bulk non-zero phase matrix elements of the THM and SCM at the wavelengths of 0.67, 2.13 and 12.0 µm in the case of an effective particle diameter of 30 µm for both models.With the surface roughness considered in the THM, the halo peaks observed in the case of pristine hexagonal columns are smoothed out, and featureless phase matrix elements are obtained with THM.In Fig. 8, the bulk extinction efficiency (Q ext ), single-scattering albedo (SSA), and asymmetry factor (g) are also given for both the THM and SCM.We call special attention to the fact that the SCM has larger asymmetry factors at all three wavelengths.Although the oscillations of the phase matrix elements of the SCM consisting of pristine ice crystals can be smoothed out by surface roughness, the effects of surface structure on the values of the integral scattering properties (e.g., the extinction efficiency and the asymmetry factor) are relatively small.
The bulk extinction efficiency, single-scattering albedo and asymmetry factor of the THM and SCM are shown in Fig. 9 as functions of the effective particle diameter.As expected, the extinction efficiency of the THM converges to 2 for larger particles, and the single-scattering albedo given by the two models are quite similar.It should be noticed that, at the 0.67 µm visible wavelength, the THM results give an almost constant asymmetry factor with a value of approximately 0.76, whereas the values for the SCM increase as the effective particle diameter increases (from 0.78 to almost 0.84).Larger asymmetry factors are also obtained with SCM at the other two wavelengths.

Comparison with observations
Compared with the large number of data sets on the microphysical properties of ice clouds and images of their geometries, our observations and understanding of the optical properties are relatively limited.The measured ice cloud phase functions have been widely used to construct and verify the numerical models (Baran et al., 2001(Baran et al., , 2012a)), and results from the THM are compared with those from laboratory and in situ measurements as well as satellite retrievals.
The PN probe has been used in various laboratories and field campaigns to measure the scattering phase function of ice clouds simultaneously with the size distribution.The PN measurements suggest that ice clouds show featureless phase functions with a relatively flat trend at backscattering angles (Barkey and Liou, 2001;Gayet et al. 1998Gayet et al. , 2004;;Febvre et al., 2009).It should be noticed that unusual scattering phase functions with certain features were also observed from some in situ measurements (Gayet et al., 2012;Baran et al., 2012), and we will not consider these special cases when building our THM.However, we compare the phase functions simulated based on the THM with measurements from laboratory and in situ measurements at a visible and a near infrared wavelength.Barkey and Liou (2001) reported the light scattering measurements of small ice crystals gen- erated in a cloud chamber at a wavelength of 0.67 µm.In situ measurements of light scattering and microphysical characteristics presented by Febvre et al. (2009) show the phase function of ice clouds at 0.804 µm.In addition, both studies measured the ice crystal number concentrations.For the case we use, Febvre et al. (2009) articulated that the effects of ice crystal shattering on the in situ measurement are probably not very important; thus, they will not be considered in our study.Figure 10 shows comparisons of the bulk phase functions between the THM and the observations (left panels), and the corresponding number concentrations are given in the right panels.The effective diameters of the two cases are approximately 5 and 35 µm, respectively.The THM exhibits a reasonable agreement in both cases.Note that the phase function of the in situ observations is normalized to the values at 30 • .Both the modeled and measured phase functions show similar and relatively smooth overall trends.The absence of halo phenomena, i.e., scattering peaks commonly seen at 22 and 46 • , is an indication of the irregularity or surface roughness of ice crystals, and demonstrates the necessity and importance of including the hollow structure and surface roughness in the THM.For the laboratory results (upper panel), the THM slightly overestimates the phase function values with scattering angles between 60 and 90  Figure 11.Comparison between the phase functions from the numerical models (both single-column and two-habit models) and MODIS retrieval at a wavelength of 1.38 µm (Wang et al., 2014).
The effective diameter used for the THM is 50 µm.Wang et al. (2014) retrieve the scattering phase function of ice clouds from satellite observations.To reduce the impact from surface reflection and highlight thin ice clouds in the upper troposphere, the reflectance at MODIS 1.38 µm channel is used to statistically derive the scattering phase function, and the phase function values at 30 scattering angles between 90 and 180 • are obtained for ice clouds over ocean and land.Figure 11 illustrates the modeled (both single-column and two-habit models) and retrieved phase functions of ice clouds, and the upper and lower panels are for the retrieved results respectively over ocean and land.The red circles in the figure represent the averaged phase functions, and the error bars indicate the standard deviations.Because the variation of a phase function with a change of effective diameter for the THM can be ignored compared with the standard deviations of the retrieved phase functions, especially for the backward scattering, the THM bulk phase functions with D eff of 50 and 100 µm are used for comparison.The phase functions of the THM at the two sizes are almost indistinguishable except for the forward peaks, illustrating that the THM-based phase function in the side and backward scattering directions are not sensitive to particle effective sizes at visible and near-infrared wavelengths.The phase function of the SCM at a single effective particle diameter of 50 µm is used.The phase functions given by the THM at two sizes almost perfectly match the retrieved values over ocean, which was also achieved by Wang et al. (2014) using three particle habits, whereas the one based on the SCM shows significant oscillation in the region.For ice cloud over land, the agreement between the satellite retrieval and numerical re-sult is relatively limited, although the modeled results are within the standard deviations of the retrieval over the entire backward direction.The THM underestimates the phase function values for scattering angles larger than 120 • .Note that, even considering the phase functions of ice habits with three different degrees of surface roughness given by Yang et al. (2013), Wang et al. (2014) cannot accurately match the inferred results over land, and this may be due to the larger uncertainties associated with the inferred phase function over land.
Overall, considering the comparisons between the phase functions calculated based on the THM and those from measurements retrievals, the THM does show excellent representation of the optical properties of ice clouds at visible and near infrared wavelengths.However, we are far from claiming the THM to be an optimal ice cloud model, and the optical properties of the THM over longer wavelengths are not verified because of a lack of observations with which to compare (Cox et al., 2010).

Satellite remote sensing applications
Both the microphysical and optical properties of the THM match the measurements closely, and another important goal in the development of the THM is to improve the consistency in the downstream remote sensing of ice cloud properties.One issue is the significant difference between ice cloud optical thicknesses retrieved from solar and infrared bands (Wang et al., 2013b;Baum et al., 2014).The polarization properties observed from the PARASOL satellite are an important aspect widely used to test ice cloud models (Baran, 2009;Cole et al., 2013).Note, the plane-parallel radiative transfer model with single cloud layer is assumed in this study, and the vertical inhomogeneity and 3-dimensional effects of clouds (Yang et al., 2001;Fauchez et al., 2014) are not considered in this study.

Comparison between the solar-and IR-band retrieved optical thicknesses
Two popular methods are normally used to retrieve ice cloud properties from satellite observations: the first is a bi-spectral method employing solar reflectance bands (the solar-band retrieval) (Nakajima and King, 1990); and the second is based on the IR bands (the IR-band retrieval) (Inoue, 1985;Heidinger and Pavolonis, 2009).To infer the optical thickness and effective particle diameter of ice clouds, an ice cloud model, i.e., the optical properties obtained from the given particle habit or habits, is fundamental for both the solarband and IR-band retrievals.Thus, identical cloud properties are expected to result from using the solar-band and IR-band retrievals for the same target based on the same ice model; however, this does not hold true for most ice cloud models.The optical thicknesses retrieved from IR-band observations Atmos.Chem.Phys., 14, 13719-13737, 2014 www.atmos-chem-phys.net/14/13719/2014/are generally smaller than those from solar-band retrievals (Baum et al., 2014).Specifically, the solar-band retrieval is based on two solar reflectance bands, i.e., a weakly absorbing, visible or nearinfrared window band (VIS / NIR) mainly sensitive to the cloud optical thickness τ , and an ice absorbing shortwave infrared (SWIR) band sensitive to both τ and D eff .The approach based on a VIS/NIR and a SWIR band is used by the MODIS operational cloud-property retrieval (Platnick et al., 2003).Another method to obtain τ and D eff is the splitwindow technique (Inoue, 1985) based on multiple IR window channels (e.g., 8.5, 11.0 and 12.0 µm for the MODIS observations), and the application of the algorithm can be found in the Advanced Very High Resolution Radiometer (AVHRR) (Heidinger and Pavolonis, 2009), as well as some studies based on MODIS observations (Minnis et al., 2011;Wang et al., 2013b).Note that the IR-band retrieval is not strongly sensitive to the explicit scattering properties, because of the strong absorption within ice crystals; whereas, the scattering properties are essential for the solar-band retrievals.Thus, the optical properties of the ice cloud model at the solar bands become the key parameters to determine the spectral consistency of the models.
A case study, using MODIS observations, is conducted to assess the spectral consistency of optical thickness retrievals based on both solar-band and IR-band observations.The solar-band retrieval uses MODIS reflectances at 0.86 and 2.13 µm bands and the fast radiative transfer model (FRTM) developed by Wang et al. (2013a).By using pre-computed bidirectional reflectance/transmittance distribution functions and a numerical integral over a twisted icosahedral mesh, the FRTM is approximately 2 orders of magnitude faster than that of the standard 128-stream discrete-ordinate radiative transfer code.The IR-band retrieval is based on the three MODIS IR bands at 8.5, 11, and 12 µm, and a fast highspectral resolution radiative transfer model (HRTM) developed by Wang et al. (2013b) is used to simulate radiances and resulting brightness temperatures at the three bands.The HRTM accounts for the gas absorption using a pre-computed transmittance database, and the optical properties of the ice cloud model are used to calculate the look-up-tables for cloud reflectance, transmittance, effective emissivity, and effective temperature functions.
Data sets from the Aqua/MODIS and the Modern Era Retrospective-Analysis for Research and Applications (MERRA) are used for the retrievals.The MODIS level-1B calibrated radiances (MYD021KM) product provides top of the atmosphere radiance/reflectance and brightness temperatures for the solar and IR bands.The 1 km-resolution geolocation and solar-satellite geometry are obtained from the MOD03 data sets.The MODIS level-2 cloud product (MYD06) is used to give cloud phase, cloud optical thickness and cloud top height.The over-ocean pixels identified as ice cloud by MYD06 are retrieved, and because the IRband retrieval is inherently less sensitive to optically thick clouds, the cases with MODIS optical thicknesses larger than 5 are ignored.The atmospheric profile used for radiative transfer simulations and gas absorption is collocated from the MERRA data.The MERRA 3-hourly instantaneous atmospheric profile data (Int3_3d_ams_Cp) provides temperatures, water vapor densities, and ozone densities at 42 pressure levels with a spatial resolution of 1.25 • × 1.25 • .Figure 12  the relatively small asymmetry factors of the THM at visible wavelengths (as shown in Fig. 9), which yield larger optical thickness from the solar-band retrieval compared with those based on the SCM.

Comparison between the simulated and observed polarized reflectivities
The polarization property of ice clouds obtained from the PARASOL satellite is important and useful perspective for evaluating the performance of numerical models, because the measured polarized reflectivity is very sensitive to the P 12 element of the phase matrix.We use PARASOL observations over ocean at 0.865 µm from 1 August 2007, and the data set details can be found in Cole et al. (2013).Data from only 1 day of observations is used, because a previous study (Baum et al., 2014) indicates that the occurrence frequency of the PARASOL polarized reflectivities exhibits a very similar pattern over time.A vector adding-doubling radiative transfer model is used (Huang et al., 2015), and the simulation assumes a single-layer ice cloud with an optical thickness of 5 at a height of 9 km over an ocean surface.Cole et al. (2013) demonstrated that an ice cloud with an optical thickness of 5 is sufficient for the polarized reflectivity to be saturated.The simulated polarized reflectivities based on the THM are not strongly sensitive to the particle effective diameter, and the optical properties of ice clouds with effective diameters of 50 µm are used.Figure 13 illustrates ice cloud polarized reflectivities from the PARASOL measurements over ocean (color contours) and the simulations (black dots) based on the bulk scattering properties developed using the single-column (left panel) and two-habit model (right panel).The same D eff of 50 µm is used for both models.The color contours in Fig. 13 are the occurrence frequency of the PARASOL polarized reflectivities of ice clouds over ocean, and the red color indicates the region of high occurrence for the measurements.The black dots in the figure correspond to the model calculations of a given set of solar-satellite geometries (i.e., solar zenith, view-ing zenith and relative azimuth angles), and 3000 different geometries from the PARASOL data are used for the simulations.Due to the scattering peaks in both the phase function and the other phase matrix elements for the smooth hexagonal column, the single-column-based results show significant oscillations as well, and exhibit very different variations than the PARASOL observations.However, the numerical results based on the THM accurately match the satellite observations over the entire range of scattering angles.Considering similar patterns in the occurrence frequency of the PARASOL polarized reflectivities over time and the similar scattering phase matrices of the THM at different D eff , the THM is expected to perform consistently in matching the observed polarized reflectivities of ice clouds.
Again, the THM not only infers similar optical thickness from the solar-band and IR-band retrievals, but also provides polarization properties similar to satellite observations.The excellent performance of the THM indicates a great potential for the remote sensing applications.More research is needed to further confirm whether the THM is a robust model for referring ice cloud properties based on observations from different wavelengths and sensors.

Conclusion
This study constructs an ice cloud model with two particle habits, and the performance and consistency of the THM in representing the microphysical and optical properties of ice clouds are investigated in detail.The THM includes a hexagonal column with an aspect ratio of unity and an aggregate containing 20 hexagonal columns, and both hollow structure and surface roughness are considered.The habit fractions of the two particle habits are determined to match the in situ measurements of ice cloud microphysical properties and the general trends, from analyses of particle imagery data sets, in the percentages of simple and complex crystals (i.e., more complex particles as particle maximum dimension increases).The simulated IWC and D mm values based on the THM agree closely with the in situ data sets.Furthermore, an expression for ice crystal volume as a function of particle maximum dimension is also presented, which leads to the aforementioned agreements in the cases of IWC and D mm .
The optical properties of the THM are calculated with a combination of the II-TM, PSTD and IGOM models for particle sizes from 2 to 10 000 µm at wavelengths of interest, and the data library contains the extinction coefficient, singlescattering albedo, asymmetry parameter, and six independent nonzero phase matrix elements.The simulated phase functions based on the THM show excellent agreement with both the laboratory and in situ measurements at 0.67 and 0.804 µm as well as satellite retrievals at 1.38 µm.
In addition to the excellent performances in representing the microphysical and optical properties of natural ice clouds, an initial retrieval analysis demonstrates that the THM significantly improves the spectral consistency in the remote sensing of ice cloud properties from different wavelengths.The optical thicknesses retrieved based on the two MODIS solar bands show close agreement with those inferred from the MODIS IR window measurements.Furthermore, a comparison between the simulated polarized reflectivities based on the THM and those measured from the PARASOL satellite indicates that the THM can closely represent the polarization properties of ice clouds.
We focused on the development and performance of the THM in representing ice cloud properties, but their effect on radiative forcing is not tested.Developing the THM optical property database over the whole spectral domain, obtaining the parameterized optical properties, performing retrievals over all ranges of viewing and illumination conditions, and investigating the radiative effects in the radiative transfer models (RTMs) and GCMs are straightforward and will be discussed in further studies.
Furthermore, we would like to emphasize that the SCM used for comparison is based on pristine ice crystals with smooth surfaces and particular aspect ratio values, and the findings based on the assessment of the performance of SCM in remote sensing applications may not necessarily be applicable to a different single column/plate model; particularly when particle surface roughness is considered.

C. Liu et al.: A two-habit model for the microphysical and optical properties of ice clouds Appendix A: Geometry of the hexagonal aggregate
The THM uses an aggregate of hexagonal columns as the complex particle, and 20 hexagonal columns with different sizes and aspect ratios are used to build the aggregate.Four steps are necessary to build the final aggregate as shown in Fig. 2b.
First, we randomly generate each of the column elements by giving its length and aspect ratio: where A 1 and A 2 are constants related to the geometries of the hexagonal columns, and ξ 1 and ξ 2 are independent random numbers distributed uniformly in [0, 1].A 1 determines the range of the column sizes, and A 2 limits the minimum aspect ratio.We use values of 0.2 and 0.8 for A 1 and A 2 , respectively, to generate the 20 hexagonal columns used to build the aggregate.Here, L o is the column average length.Once the aggregate is generated, the dimensions can be scaled to fit the ice crystal size in the single-scattering computations.
Secondly, the 20 hexagonal columns are attached to form an aggregate without overlapping.An improved particlecluster aggregation algorithm, normally used for a fractal aggregate (Filippov et al., 2000;Liu et al., 2012c) with spherical monomers, is adapted, and the only difference is that each hexagonal column is randomly rotated to attach to another.For simplification, the rotation is managed to make only a vertex and a surface point attached (without overlapping or surface-surface attachment).The criteria used by Xie et al. (2011) to detect overlapping between two hexagonal particles are used to avoid intersecting particle faces.
The third step is to introduce a hollow structure into the hexagonal columns by replacing one hexagonal surface of each column using the hollow structure as shown in Fig. 1.The depth of the hollow is fixed at d / L = 0.25.To ensure the attachment of the aggregate, the hollow structure is only added to a surface without any attached particles, and a hexagonal column with both hexagonal surfaces connected with other monomers is kept solid.The aggregate has only one solid column.
As a final step, surface roughness is added to the particle by replacing each of the smooth surfaces with roughened ones.In the II-TM and PSTD simulations, explicit particle geometries are achieved by following the rough surfaces defined by Liu et al. (2013).The titled-facet approximation (Yang et al., 2008) is applied for the IGOM simulations because of the efficiency without significant loss of accuracy (Liu et al., 2013).
The completed roughened aggregate is shown in panel b of Fig. 2. Numerically, the aggregate is defined by the explicit vertices and surfaces of the columns.Thus, the volume and averaged projected area of the aggregate can be rigorously and numerically calculated, and will be used for the microphysical and optical properties of the THM.Note that the surface roughness has little effect on the microphysical properties of the THM, but its influence is shown in the optical properties.
The geometrical parameters used to determine the aggregate geometry are listed in Table A1, including the lengths, aspect ratios, coordinates of three points and hollow depths of 20 monomers.It should be noticed that all length parameters are normalized to L o .The coordinates of three points for each monomer, i.e., the center of a column (point O in Fig. 1), the center of a particle face (point A in Fig. 1), and a vertex (point B in Fig. 1) are listed.The last column in Table A1 indicates whether a monomer has a hollow structure as shown in Fig. 1.The maximum dimension of the aggregate is numerically calculated, which is 7.137 L o .In addition, the volume and projected area of the aggregate are numerically found to be 0.0255 D 3 and 0.260 D 2 , respectively.
Although Baran (2009) demonstrated that adding hexagonal monomers with the element number beyond 3 does not significantly alter the asymmetry factor, in this study we select 20 monomers for three reasons: (1) as an appropriate particle geometry is sought to mimic the complicated morphologies of realistic aggregates within ice clouds and the use of only a few monomers seems to be an oversimplification; (2) an aggregate geometry corresponding to a potentially lowest value of the asymmetry factor is desired, and it is found that the asymmetry factor slightly decreases as the number of monomers increases; (3) with the trial and error method, the use of 20 monomers is optimal in terms of the balance between the computational effort in light scattering simulation and the performance of the particle habit model in fitting the measured microphysical properties (specifically, IWC and D mm ).

Figure 1 .
Figure 1.Geometry of a hexagonal column with a hollow structure.L is equal to the maximum dimension D for the hexagonal column.

Figure 2 .
Figure 2. Particle geometries for the two-habit model (THM): (a) single hexagonal column with an aspect ratio of unity, and (b) hexagonal aggregate with 20 solid or hollow columns.

Figure 3 .
Figure 3. Ice crystal habit fraction as a function of particle maximum dimension for the two-habit model.

Figure 4 .
Figure 4. Upper panels: comparison between the measured and calculated microphysical properties (D mm and IWC) for each of the PSDs from 11 field campaigns.Lower panels: Histograms of the distributions of the measured and calculated D mm and IWC.

Figure 5 .
Figure 5. Relative differences (RD) of the calculated microphysical properties at different bins of median mass diameter (left panels) and ice water content (right panels).Error bars indicate the corresponding standard deviations.

Figure 6 .
Figure 6.Relationship between ice crystal volume and maximum dimension.

Figure 7 .
Figure7.Extinction efficiency (upper), single-scattering albedo (middle), and asymmetry factor (lower) of the single-column and two-habit model as functions of particle maximum dimension at wavelengths of 0.67, 2.13 and 12.0 µm.

Figure 8 .
Figure 8.Comparison of bulk non-zero phase matrix elements of the two-habit model and the single-column model with an effective diameter of 30 µm at wavelengths 0.67, 2.13 and 12.0 µm.

Figure 10 .
Figure 10.Comparison between the phase functions (left panels) from the two-habit model and the polar nephelometer (PN) measurements from the following: (a) laboratory at a wavelength of 0.67 µm and (b) in situ at a wavelength of 0.804 µm.The right panels are observed particle number concentration of the corresponding measurements.

Figure 12 .
Figure 12.(a) RGB image of an Aqua/MODIS granule from 24 February 2014 at 03:50 (b) Retrieved optical thickness of thin ice clouds.(c) and (d) are comparisons of ice cloud optical thicknesses retrieved from the MODIS solar bands and IR bands, and the results are based on the single-column model and two-habit model, respectively.The histograms illustrate occurrences of thin ice cloud pixels, and the red color corresponds to the high frequency of occurrence.

Figure 13 .
Figure 13.Comparisons between the normalized polarized reflectivities obtained from 1 day of PARASOL data over ocean (color contours) and calculations (black dots) based on the single-column model (left) and the two-habit model (THM, right).

Table 1 .
Baum et al. (2014)(2013)HM in representing the microphysical properties obtained during 11 field campaigns.The mean and standard deviation (SD) of the relative errors of the theoretical median mass diameter (D mm ) and ice water content (IWC) are listed.The details of the 11 field campaigns can be found inHeymsfield et al. (2013)andBaum et al. (2014).

2014 13728 C. Liu et al.: A two-habit model for the microphysical and optical properties of ice clouds
• compared with the in situ observations (lower panel).The asymmetry factors of the laboratory and in situ measurement is approximately 0.76 and 0.79, respectively, and the corresponding modeled values of the THM are 0.77 and 0.78.www.atmos-chem-phys.net/14/13719/2014/Atmos.Chem.Phys., 14, 13719-13737,

Table A1 .
Geometric parameters of the hexagonal aggregate with 20 monomers.All length parameters are normalized by L o , the average monomer length.The last column indicates whether the monomer has a hollow structure as shown by Fig.1("Y" indicates "yes", and "N" indicates "no").Note, only the third monomer does not have a hollow structure.