Inference of an Optimal Ice Particle Model through Latitudinal Analysis of MISR and MODIS Data

The inference of ice cloud properties from remote sensing data depends on the assumed forward ice particle model, as they are used in the radiative transfer simulations that are part of the retrieval process. The Moderate Resolution Imaging Spectroradiometer (MODIS) Collection 6 (MC6) ice cloud property retrievals are produced in conjunction with a single-habit ice particle model with a fixed degree of ice particle surface roughness (the MC6 model). In this study, we examine the MC6 model and five other ice models with either smoother or rougher surface textures to determine an optimal model to reproduce the angular variation of the radiation field sampled by the Multi-angle Imaging Spectroradiometer (MISR) as a function of latitude. The spherical albedo difference (SAD) method is used to infer an optimal ice particle model. The method is applied to collocated MISR and MODIS data over ocean for clouds with temperatures ≤233 K during December solstice from 2012–2015. The range of solar zenith angles covered by the MISR cameras is broader at the solstices than at other times of the year, with fewer scattering angles associated with sun glint during the December solstice than the June solstice. The results suggest a latitudinal dependence in an optimal ice particle model, and an additional dependence on the solar zenith angle (SZA) at the time of the observations. The MC6 model is one of the most optimal models on the global scale. In further analysis, the results are filtered by a cloud heterogeneity index to investigate cloudy scenarios that are less susceptible to potential 3D effects. Compared to results for global data, the consistency between measurements and a given model can be distinguished in both the tropics and extra-tropics. The SAD analysis suggests that the optimal model for thick homogeneous clouds corresponds to more roughened ice particles in the tropics than in the extra-tropics. While the MC6 model is one of the models most consistent with the global data, it may not be the most optimal model for the tropics.


Introduction
The inference of ice cloud optical thickness τ and effective particle size r eff from passive spaceborne radiometric measurements requires an assumed ice particle model that provides the bulk scattering and absorption properties.Based on the ice particle model, look-up tables (LUTs) for cloud property retrieval are generated by using radiative transfer simulations.The LUTs provide the transmission, scattering, and emission characteristics as functions of, for example, optical thickness and the sun-satellite geometric configuration (e.g., solar zenith angle, viewing zenith angle, and relative azimuth angle).In practice, the LUTs are applied to global satellite measurements so the assumed ice particle model should be applicable to a large spatial domain.
There are numerous constraints on choosing an ice particle model for global data processing.For operational retrievals, the Clouds and the Earth's Radiant Energy System (CERES) and Moderate Resolution Imaging Spectroradiometer (MODIS) adopted single-habit models.Specifically, CERES Edition 4 adopted a severely roughened hexagonal column model, whereas MODIS Collection 6 (MC6) adopted a model with aggregates consisting of severely roughened columns [1,2].In addition, a Voronoi particle model was suggested for use with geostationary satellite data [3] and the inhomogeneous hexagon model (IHM) was defined for the Polarization and Directionality of the Earth's Reflectances (POLDER) on the Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar (PARASOL) satellite ice cloud property retrievals.Multiple models are not generally adopted by an individual team to avoid potential discontinuities in retrievals resulting from the transition between models.
Recent research suggests that the ice particle models should have a substantial degree of surface roughness [4,5], or at least some amount of inhomogeneity.The IHM model employs a different approach than surface roughening to increase photon dispersion by including air bubbles.Generally speaking, the surface roughness or inhomogeneity of ice particles tends to smooth the phase function and results in a relatively low asymmetry factor at solar wavelengths.Because the phase function is fundamental to the remote sensing of global cloud properties, a better understanding of the appropriate degree of ice particle surface roughness for a given ice particle habit is important for improving the consistency of the retrievals based on observations by different sensors.
The overarching goal of this study is to identify an appropriate degree of surface roughness adopted for ice particle models through the use of collocated Multi-angle Imaging Spectro Radiometer (MISR) and MODIS data.Compared to a satellite sensor with only a nadir-viewing camera, multi-angle cameras measure the reflectance of a cloud over a wide range of scattering angles as the satellite passes over a location.Because different physical ice particle habits lead to significantly different angular distributions of reflectance simulated at the top of the atmosphere (TOA), comparing theoretical radiative transfer simulations with multi-angle camera measurements provides valuable constraints on the particle morphology.Several previous studies [6][7][8][9][10] attempted to infer the predominant atmospheric ice particle habits based on multi-angle satellite measurements.In particular, the spherical albedo difference (SAD) method [6] was developed to quantify the comparison between spherical albedo values computed with an assumed ice particle model and their counterparts derived from multi-angle satellite measurements.Furthermore, this algorithm was used to validate cloud particle models and their single-scattering properties [11][12][13].
This study assesses the latitudinal dependence of six ice cloud models based on the application of the SAD method to the fused MISR and MODIS 0.86-µm channel data.Section 2 describes the satellite measurements and the SAD method.Section 3 presents the results, including the latitudinal consistency of ice particle models with MISR observations via the SAD method, including analyses under different cloud heterogeneity conditions.Section 4 discusses potential uncertainties of the results.A summary and conclusions are given in Section 5.

Data and Methods
In the SAD method, the spherical albedo difference value (A diff ) is computed for every pixel in the spatial domain of interest by comparing theoretical radiances and the measured counterparts.To compute A diff values for ice clouds based on MISR observations, we apply a look-up table approach in this study.Section 2.1 summarizes available satellite data and the procedure for selecting only ice cloud pixels, Section 2.2 describes the ice cloud particle models used for comparisons with observations, Section 2.3 describes LUTs, and Section 2.4 describes the application of the SAD method to define an optimal ice cloud particle model.

Satellite Data
In this study, we use collocated datasets obtained by the MISR and MODIS instruments onboard the Terra satellite during December solstices from 2012-2015.The details of collocation are described by Liang et al. [14] and Liang and Di Girolamo [15].The MODIS product at 1-km resolution and all MISR camera views are co-registered to the MISR nadir camera (AN) pixel positions to generate a 1.1-km resolution fused dataset.Measured radiances used in this study are from the MISR observations, and the MC6 products are used to classify cloud pixels.
MISR uses nine cameras to measure radiance along the satellite track [16,17].In addition to one nadir-viewing camera (AN), four cameras (AF, BF, CF, DF) point forward and four cameras (AA, BA, CA, DA) point aft along the orbital track.The viewing zenith angles (VZAs) for the AA/AF, BA/BF, CA/CF, DA/DF cameras are 26.1 • , 45.6 • , 60.0 • , and 75.0 • , respectively.The three near-nadir-viewing cameras (AA, AN, AF) are used to avoid potential 3D effects due to large VZAs.Each camera measures radiances in four narrow spectral bands.The band centered at 0.86 µm is selected because this channel is less affected by ozone and ice absorption.The measured radiances are converted to reflectances ( R) as follows: where I is the measured radiance, P is the phase function, µ 0 is the cosine of solar zenith angle (SZA), µ is the cosine of VZA, φ 0 is solar azimuthal angle, φ is viewing azimuthal angle, d is Sun-Earth distance in astronomical units (AU), and E 0 is the solar irradiance at 1 AU.
To avoid cloud pixels containing liquid phase particles, pixels are selected by applying two criteria based on the MODIS products: (1) cloud phase identified with infrared channels as ice; and (2) cloud top temperature lower than 233 K. To avoid potential effects of land reflectance and associated complexity for radiative transfer computations, observations are limited to an ocean surface.All MISR camera measurements are removed that have sun glint angles smaller than 35 • , or specifically satellite-viewing directions within a 35 • cone around the sunlight direction.Avoiding sun glint reduces the number of selected camera views (n c ) for some pixels.Moreover, pixel locations are restricted to 60 • N-60 • S to avoid the effects of sea ice since an ice-free ocean surface is assumed.
The cloud optical thickness (τ) and cloud heterogeneity index (H σ ) from MODIS products are used to stratify clouds for statistical analysis.The MODIS cloud optical thickness product is not used elsewhere in this study, such as in Sections 2.3 and 2.4.H σ estimates the degree of cloud horizontal heterogeneity, which is computed with the 4 × 4 sub-pixel array composing a MODIS 1-km pixel (MODIS has two channels at 0.25-km resolution).The H σ is defined as the standard deviation divided by the mean measured reflectance for such a pixel [18].
SZA generally increases with increasing latitude.To reduce differences in the SZAs at the same latitude, data are collected and analyzed on nearly the same date for four years.Due to the solar-earth-satellite viewing geometry, the range of SZA values covered by the MISR cameras is broader at the solstices than at other times of the year, and fewer scattering angles are associated with sun glint on the December solstice in comparison to the June solstice.In other words, using data on the December solstice provides information over the broadest scattering angle range with fewer sun glint issues.For these reasons, this study uses data during December solstices from 2012-2015, except 2013.In 2013, we selected 25 December 2013, not the solstice, to avoid complexities caused by a cyclonic storm over the Indian Ocean.The specific dates chosen for detailed analyses are December 21, 2012, December 25, 2013, December 21, 2014, and December 22, 2015.

Ice Particle Models
Before applying the SAD method to the MISR datasets, LUTs are prepared by assuming specific ice particle models.Six ice particle models are employed in this study to examine the consistency between the theoretical radiances and the MISR observations at various latitudes.
The MC6 roughened hexagonal ice aggregate model is used in the MODIS Collection 6 operational products [19], which assume ice cloud particles to be an ensemble of randomly oriented aggregates.Each aggregate is composed of eight hexagonal columns with a fixed degree of surface roughness σ 2 = 0.5.The roughness parameter σ 2 in the light scattering calculations is defined as the standard deviation of a two-dimensional Gaussian distribution for the tilting of a particle facet [20].The standard deviation of this roughness parameter σ is approximately equivalent in value to the roughness parameter δ defined by Macke et al. [21][22][23].A higher degree of roughness means a higher probability density of strongly distorted surfaces.Detailed descriptions of the ice particle habit and associated degree of surface roughness are discussed by Yang and Liou [20] and Yang et al. [24].In general, the ice particle phase function at scattering angles between 50 • and 175 • becomes more featureless with an increasing degree of surface roughness (Figure 1) for six models assuming a range of surface roughness of σ 2 = 0.001, 0.03, 0.14, 0.5, 1.0, and 3.5 (hereafter, R0001, R003, R014, R05, R10, and R35, respectively).Note that the MC6 model corresponds to the R05 model.Both the ice particle habit and ice particle surface roughness can modify the single-scattering properties [25], although other factors are also important, such as particle impurities, internal fractures, and air bubbles.The asymmetry factor (g) for each model is also included in Figure 1.Among these selected models, the asymmetry factor decreases with an increasing degree of roughness until σ 2 = 0.14, beyond which the asymmetry factor increases in even rougher models.Perhaps some of this behavior in the most severely roughened models could be a modeling artifact associated with the ray-tracing technique when it is applied to a very rough particle surface.Although the MC6 ice particle habit model has a better match to global satellite measurements than other ice particle habits [1,4,26], the "optimal" degree of ice particle surface roughness for the MC6 ice model has not been studied rigorously yet.To test the degree of roughness used in MC6, we developed five other ice particle models using the MC6 ice particle habit, but with different degrees of ice particle surface roughness.The additional five degrees of roughness have different phase functions over the MISR observational range of scattering angles in comparison with the case of σ 2 = 0.5.
Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 17 The MC6 roughened hexagonal ice aggregate model is used in the MODIS Collection 6 operational products [19], which assume ice cloud particles to be an ensemble of randomly oriented aggregates.Each aggregate is composed of eight hexagonal columns with a fixed degree of surface roughness σ 2 = 0.5.The roughness parameter σ 2 in the light scattering calculations is defined as the standard deviation of a two-dimensional Gaussian distribution for the tilting of a particle facet [20].The standard deviation of this roughness parameter σ is approximately equivalent in value to the roughness parameter δ defined by Macke et al. [21][22][23].A higher degree of roughness means a higher probability density of strongly distorted surfaces.Detailed descriptions of the ice particle habit and associated degree of surface roughness are discussed by Yang and Liou [20] and Yang et al. [24].In general, the ice particle phase function at scattering angles between 50° and 175° becomes more featureless with an increasing degree of surface roughness (Figure 1) for six models assuming a range of surface roughness of σ 2 = 0.001, 0.03, 0.14, 0.5, 1.0, and 3.5 (hereafter, R0001, R003, R014, R05, R10, and R35, respectively).Note that the MC6 model corresponds to the R05 model.Both the ice particle habit and ice particle surface roughness can modify the single-scattering properties [25], although other factors are also important, such as particle impurities, internal fractures, and air bubbles.The asymmetry factor (g) for each model is also included in Figure 1.Among these selected models, the asymmetry factor decreases with an increasing degree of roughness until σ 2 = 0.14, beyond which the asymmetry factor increases in even rougher models.Perhaps some of this behavior in the most severely roughened models could be a modeling artifact associated with the ray-tracing technique when it is applied to a very rough particle surface.Although the MC6 ice particle habit model has a better match to global satellite measurements than other ice particle habits [1,4,26], the "optimal" degree of ice particle surface roughness for the MC6 ice model has not been studied rigorously yet.To test the degree of roughness used in MC6, we developed five other ice particle models using the MC6 ice particle habit, but with different degrees of ice particle surface roughness.The additional five degrees of roughness have different phase functions over the MISR observational range of scattering angles in comparison with the case of σ 2 = 0.5.

Look-Up Table Approach
LUTs are developed separately for each ice particle model, i.e., R0001, R003, R014, R05, R10, and R35.For each ice cloud particle model, the LUTs are calculated using an adding-doubling radiative transfer model [27].The LUTs include model reflectances as a function of solar geometry, viewing geometry, and τ; the model reflectance R is equivalent to the reflectance R from MISR measurements in Equation ( 1).The bulk scattering properties for each ice particle model are the ensemble-mean single-scattering properties integrated over a Gamma distribution with an effective variance of 0.1 and an effective radius of 30 µm.The reflectance at wavelength 0.86 µm is insensitive to the particle effective size due to weak ice absorption.Atmospheric molecular scattering is considered in the model, but aerosols are neglected.The ocean surface reflection is based on the Cox-Munk model [28] with a wind speed of 10 m s −1 .The reflectance for each MISR camera is calculated by assuming a homogeneous cloud layer with a cloud-top pressure of 200 hPa.A sensitivity study shows that the effect of cloud-top pressure on retrievals is not substantial.
The R values can be integrated over µ and φ (viewing directions) to obtain the planetary albedo (A p ): The cloud spherical albedo, A s , is the integration of A p over all incident directions (µ 0 , φ 0 ): By carrying out the integrals in Equations ( 2) and ( 3), R is converted to A s for a given τ in conjunction with a specific ice particle model.LUTs are developed for each ice particle model and contain both R and the corresponding A s as functions of τ.

The SAD Method
The SAD method for examining angular variations of ice particle phase functions is described by Doutriaux-Boucher et al. [6] and C.-Labonnote et al. [29], who applied this method to POLDER data.These studies found that A s has a one-to-one relationship with τ for a given ice particle model [6], implying that there is a one-to-one theoretical relationship between A s and R. Because R and A s are provided in the LUTs, an R value from each MISR camera measurement is used to search for the corresponding value of R in the LUTs.The associated A s in the LUTs is identified as the retrieved cloud spherical albedo ( A s ) for this camera.This procedure is iteratively performed for every camera measurement to compute A s .Because only three near-nadir cameras are selected from a set of MISR measurements, one cloudy pixel may correspond to 3 A s values.
Each MISR camera retrieves R for a pixel at a different scattering angle (Θ).If an assumed ice particle model matches observations accurately over the range of three observed Θ values, the model R values would be equal to the observed R values for every Θ.When τ is identical for all cameras, the same A s could be obtained.Given that idealized ice models may not fit a real cloud perfectly for every available Θ, there may be differences in the computed A s values.A di f f is defined as the difference between the A s of a selected camera and the mean of A s averaged for the selected cameras for a given pixel: where i is the camera index in the pixel and n c is the total number of selected cameras in the given pixel.Because of the sun glint screening process, n c is not necessarily 3.As defined in Equations ( 2) and (3), A s does not depend on either solar or viewing geometries but solely on P for a given τ.Therefore, smaller absolute values of A di f f indicate better consistency between simulations and observations.To identify an optimal ice particle model, the same procedure is applied to each pixel using all LUTs.The smallest A diff value from the six LUTs is chosen as the optimal model for that pixel.Finally, to quantify the consistency between ice particle models with observations, the standard deviation of A diff (χ 2 ) for each ice particle model in a given latitude bin is defined as follows: where n R is the total number of selected reflectances in a latitude-scattering angle bin and n b is the total number of selected latitude-scattering angle bins in a particular latitude band.VZAs of each MISR camera, the cameras provide reflectances at different Θ for a given pixel, and the Θ changes with latitude.The minimum available Θ is 76 • and the maximum available Θ is 172 • in these selected pixels.Because the daytime orbit of the Terra satellite is from north to south with an equator crossing time at 10:30 am, the forward camera measures at smaller Θ than the aft camera in the northern high latitudes.The three cameras primarily take measurements in the side and backward scattering directions.Note that the presence of arc-shaped strips in Figure 2a is a result of filtering out the pixels with sun glint.For the same reason, some bins in Figure 2b have no observations.and (3),  does not depend on either solar or viewing geometries but solely on P for a given .Therefore, smaller absolute values of  indicate better consistency between simulations and observations.

Sampling Scattering Geometry Characteristics
To identify an optimal ice particle model, the same procedure is applied to each pixel using all LUTs.The smallest  value from the six LUTs is chosen as the optimal model for that pixel.Finally, to quantify the consistency between ice particle models with observations, the standard deviation of  (χ ) for each ice particle model in a given latitude bin is defined as follows: where  is the total number of selected reflectances in a latitude-scattering angle bin and  is the total number of selected latitude-scattering angle bins in a particular latitude band.VZAs of each MISR camera, the cameras provide reflectances at different Θ for a given pixel, and the Θ changes with latitude.The minimum available Θ is 76° and the maximum available Θ is 172° in these selected pixels.Because the daytime orbit of the Terra satellite is from north to south with an equator crossing time at 10:30 am, the forward camera measures at smaller Θ than the aft camera in the northern high latitudes.The three cameras primarily take measurements in the side and backward scattering directions.Note that the presence of arc-shaped strips in Figure 2a is a result of filtering out the pixels with sun glint.For the same reason, some bins in Figure 2b have no observations.The detected range of Θ changes with latitude, in large part due to changes in SZA, as shown in Figure 2c.The camera geometry discussed above causes the detected range of Θ variation with SZA.The narrowest range of scattering angle measurements occurs at the lowest SZA (~20 • S), where no reflectances are measured at Θ < 140 • .However, all reflectances are measured at Θ < 140 • at the highest SZA (~50 • -60 • N).The latitudes where measurements are made at the largest Θ are different for each camera, but roughly speaking, all three cameras measure in side scattering to backscattering directions with increasing SZA.The frequencies of Θ observations also reveal latitudinal variations.When the SZA is low, the same scattering angles can be recorded by two cameras viewing the same pixel.

Consistency of Models and Measurements with Latitude
The quantity χ 2 in Equation ( 5) as a function of latitude for each particle model is shown in Figure 3a, and quantitative comparisons between the R05 and other models based on Equation ( 5) are displayed in Figure 3b.The χ 2 values for all models in the Northern Hemisphere are larger than in the Southern Hemisphere and decrease with latitude from north to south.The six χ 2 values are more similar in value to each other in the low latitudes than in the high latitudes of both hemispheres.Figure 3 shows that the R05 and R014 models have similar χ 2 values over the latitude range, and both models have lower χ 2 values than the other models in all latitudes.Note that these two models have lower g than the other models (Figure 1).However, the relationships between χ 2 and g are not simple.The R35 model has the highest g but has a lower χ 2 value than the R0001 model in most latitudes.Also, the R0001 and R10 models have similar g values, but the consistency of A diff results for the R0001 model is much less than in the case of the R10 model.The detected range of Θ changes with latitude, in large part due to changes in SZA, as shown in Figure 2c.The camera geometry discussed above causes the detected range of Θ variation with SZA.The narrowest range of scattering angle measurements occurs at the lowest SZA (~20° S), where no reflectances are measured at Θ < 140°.However, all reflectances are measured at Θ < 140° at the highest SZA (~50°-60° N).The latitudes where measurements are made at the largest Θ are different for each camera, but roughly speaking, all three cameras measure in side scattering to backscattering directions with increasing SZA.The frequencies of Θ observations also reveal latitudinal variations.When the SZA is low, the same scattering angles can be recorded by two cameras viewing the same pixel.

Consistency of Models and Measurements with Latitude
The quantity χ in Equation ( 5) as a function of latitude for each particle model is shown in Figure 3a, and quantitative comparisons between the R05 and other models based on Equation ( 5) are displayed in Figure 3b.The χ values for all models in the Northern Hemisphere are larger than in the Southern Hemisphere and decrease with latitude from north to south.The six χ values are more similar in value to each other in the low latitudes than in the high latitudes of both hemispheres.Figure 3 shows that the R05 and R014 models have similar χ values over the latitude range, and both models have lower χ values than the other models in all latitudes.Note that these two models have lower g than the other models (Figure 1).However, the relationships between χ and g are not simple.The R35 model has the highest g but has a lower χ value than the R0001 model in most latitudes.Also, the R0001 and R10 models have similar g values, but the consistency of  results for the R0001 model is much less than in the case of the R10 model.To better understand the contributions of  to χ , Figure 4 shows the median value of  (hereafter,  indicates the median value in 5° × 5° latitude-Θ bins) on December solstices during 2012-2015.The variations of  values more closely follow the changes in SZA.The axis of symmetry of  is located at about 20° S, where there is a minimum in the Θ range.For SZA < 30°, the  values computed with all six models display nearly the same pattern.The  values are close to zero for most Θ, but slightly negative at the highest Θ (Θ~170°).As SZA increases,  values are still close to zero when 30° < SZA < 50°.However,  values become positive for the To better understand the contributions of A diff to χ 2 , Figure 4 shows the median value of A diff (hereafter, A diff indicates the median value in 5 • × 5 • latitude-Θ bins) on December solstices during 2012-2015.The variations of A diff values more closely follow the changes in SZA.The axis of symmetry of A diff is located at about 20 • S, where there is a minimum in the Θ range.For SZA < 30 • , the A diff values computed with all six models display nearly the same pattern.The A diff values are close to zero for most Θ, but slightly negative at the highest Θ (Θ~170 • ).As SZA increases, A diff values are still close to zero when 30 • < SZA < 50 • .However, A diff values become positive for the largest measurable Θ with an increasing degree of roughness and become negative for the smallest measurable Θ with a decreasing degree of roughness.At high latitudes (SZA > 50 • ), the A diff values broadly become highly positive or negative as shown in Figure 4a,b,e,f (especially in Figure 4a,f), but not in Figure 4c,d.The absence of large extreme values at high latitudes in Figure 4c,d indicate that the R014 and R05 models have lower χ 2 values than the other models in high latitudes in Figure 3.Some of the latitudinal A diff value differences could be a result of the change of measurable Θ as shown in Figure 2a.Another issue that might influence the results is incomplete filtering of sun glint regions.Note that Figure 4a has significantly negative A diff values in most latitudes at Θ~150 • where the phase function of the R0001 model has an obvious peak in Figure 1.In general, the consistency of A di f f results for the R05 and R014 models is better than for the other models in high latitudes, and similar in the tropics.
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 17 largest measurable Θ with an increasing degree of roughness and become negative for the smallest measurable Θ with a decreasing degree of roughness.At high latitudes (SZA > 50°), the  values broadly become highly positive or negative as shown in Figure 4a,b,e,f (especially in Figure 4a,f), but not in Figure 4c,d.The absence of large extreme values at high latitudes in Figure 4c,d indicate that the R014 and R05 models have lower χ values than the other models in high latitudes in Figure 3.Some of the latitudinal  value differences could be a result of the change of measurable Θ as shown in Figure 2a.Another issue that might influence the results is incomplete filtering of sun glint regions.Note that Figure 4a has significantly negative  values in most latitudes at Θ~150° where the phase function of the R0001 model has an obvious peak in Figure 1.In general, the consistency of  results for the R05 and R014 models is better than for the other models in high latitudes, and similar in the tropics.

Sensitivity Test with Synthetic Dataset
To investigate which ice particle model produced  values at high SZAs such as the case shown in Figure 4, a sensitivity study is performed with simulated data (Figure 5).The simulated radiances are generated using the same December solstice observations and the same pixel filtering

Sensitivity Test with Synthetic Dataset
To investigate which ice particle model produced A diff values at high SZAs such as the case shown in Figure 4, a sensitivity study is performed with simulated data (Figure 5).The simulated radiances are generated using the same December solstice observations and the same pixel filtering process as in Section 2.1, the viewing geometry from each MISR camera observation, and the τ from MODIS with the R05, R10, and R35 models.In other words, the simulated radiances are constructed using the same satellite geometry as the MISR data but the reflectances are generated with the R05, R10, and R35 model LUTs using the τ from the collocated MODIS MC6 product.Subsequently, A diff values are computed using the LUTs from the other two ice particle models.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 17 process as in Section 2.1, the viewing geometry from each MISR camera observation, and the τ from MODIS with the R05, R10, and R35 models.In other words, the simulated radiances are constructed using the same satellite geometry as the MISR data but the reflectances are generated with the R05, R10, and R35 model LUTs using the  from the collocated MODIS MC6 product.Subsequently,  values are computed using the LUTs from the other two ice particle models.Compared to the simulated dataset produced by the R10 model (Figure 5c,d) and the R35 model (Figure 5e,f), A diff values computed from the simulated dataset produced by the R05 model (Figure 5a,b) have better consistency with A diff values computed with real data (Figure 4) at high SZAs.Specifically, Figures 5a and 4e illustrate A diff values based on the same R10 LUTs.The simulated A diff pattern in Figure 5a reproduces the negative-positive patterns of A diff values (i.e., positive values of A diff become negative with increasing Θ in a latitude) at high latitudes in Figure 4e.Similar negative-positive patterns in high latitudes are seen in both Figures 4f and 5b, which use the same R35 LUTs.However, neither Figure 5d nor Figure 5f are similar to Figure 4f even using the same R35 LUTs.The similarity of using the same LUTs between Figures 4 and 5 indicate that the R05 model is closer to the actual ice particle in the measurable range of scattering angles at high latitudes.With use of the same approach to compare Figures 4 and 5 under SZA < 30 • , it is found that both the R05 and R10 models do not match the actual ice cloud particle as well as the R35 model when SZA < 30 • , but the difference with the R35 results are rather small.To summarize, the similarity between Figures 4 and 5 implies that the R05 model fits observations better than the R10 and R35 models for SZA > 30 • .However, it is not straightforward to select the most optimal model when SZA < 30 • .

Classification by Heterogeneity of Clouds
To better understand the effect of cloud heterogeneity on variations of A diff values in latitude, we further classify the A diff values by H σ and τ. Figure 6 shows A diff computed with the R05 model in a 4 × 4 matrix with four τ bins: 0-3, 3-8, 8-16, and 16-64, and four H σ bins: <0.4,0.4-1.6,1.6-3.2, and 3.2-15.A previous study [30] selected 0.3 as the threshold for a homogenous cloud, but very few selected pixels in this study meet this threshold.Therefore, a higher threshold (i.e., H σ = 0.4) was used here to select homogeneous clouds (i.e., the first row in Figure 6).The same pixel filtering process as in Section 2.1 was applied here.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 17 Similar negative-positive patterns in high latitudes are seen in both Figures 4f and 5b, which use the same R35 LUTs.However, neither Figure 5d nor Figure 5f are similar to Figure 4f even using the same R35 LUTs.The similarity of using the same LUTs between Figures 4 and 5 indicate that the R05 model is closer to the actual ice particle in the measurable range of scattering angles at high latitudes.With use of the same approach to compare Figures 4 and 5 under SZA < 30°, it is found that both the R05 and R10 models do not match the actual ice cloud particle as well as the R35 model when SZA < 30°, but the difference with the R35 results are rather small.To summarize, the similarity between Figures 4 and 5 implies that the R05 model fits observations better than the R10 and R35 models for SZA > 30°.However, it is not straightforward to select the most optimal model when SZA < 30°.

Classification by Heterogeneity of Clouds
To better understand the effect of cloud heterogeneity on variations of  values in latitude, we further classify the  values by  and  .Figure 6 shows  computed with the R05 model in a 4 × 4 matrix with four  bins: 0-3, 3-8, 8-16, and 16-64, and four  bins: <0.4,0.4-1.6,1.6-3.2, and 3.2-15.A previous study [30] selected 0.3 as the threshold for a homogenous cloud, but very few selected pixels in this study meet this threshold.Therefore, a higher threshold (i.e.,  = 0.4) was used here to select homogeneous clouds (i.e., the first row in Figure 6).The same pixel filtering process as in Section 2.1 was applied here.The  values become more negative or more positive with increasing  in Figure 6.The larger absolute value of  for larger  show the R05 results are more consistent with measurements for optically thin clouds than in optically thick clouds.Unlike the impact of increasing τ, the sign of  changes with increasing  .In a broad sense,  decreases with increasing Θ in low  bins but increases with increasing Θ in high  bins.Compared with all  bins, 1.6 <  < 3.2 bins have a trend of  with increasing Θ that is similar to the trend in Figure 4a.When  increases, the decreasing trend of  values with increasing Θ reverses to an increasing trend, indicating that the scattering angular distribution of reflectance depends on cloud homogeneity.The A diff values become more negative or more positive with increasing τ in Figure 6.The larger absolute value of A diff for larger τ show the R05 results are more consistent with measurements for optically thin clouds than in optically thick clouds.Unlike the impact of increasing τ, the sign of A diff changes with increasing H σ .In a broad sense, A diff decreases with increasing Θ in low H σ bins but increases with increasing Θ in high H σ bins.Compared with all H σ bins, 1.6 < H σ < 3.2 bins have a trend of A diff with increasing Θ that is similar to the trend in Figure 4a.When H σ increases, the decreasing trend of A diff values with increasing Θ reverses to an increasing trend, indicating that the scattering angular distribution of reflectance depends on cloud homogeneity.

Latitudinal Variations of Consistency in Measurements of Thick Homogeneous Clouds
Since low H σ clouds are less likely to be influenced by 3D effects, the results are recomputed for thick homogeneous clouds (H σ < 0.4 and τ > 16).Figures 7 and 8 are the same as Figures 3 and 4, but only for thick homogeneous clouds.Compared to the results for global clouds (Figure 3), the results for thick homogeneous clouds show that: (1) the χ 2 of all selected models are significantly lower in thick homogeneous clouds; (2) the absolute differences of χ 2 among models are much larger in both the tropics and extra-tropics; and (3) the χ 2 of each model in the Northern Hemisphere is not significantly larger than χ 2 in the Southern Hemisphere.While the R05 model generally has the lowest χ 2 in all latitudes in Figure 3, the χ 2 R10 and χ 2 R35 results are smaller than the χ 2 R05 counterparts at latitudes between 20 • S and 30 • N in Figure 7.In latitudes between 20 • S-60 • S and between 30 • N-60 • N, the R05 model and the less roughened ice models produce better fits among the six.In other words, the optimal model suggests the need for rougher ice particles in the tropics than in the extra-tropics.Furthermore, Figure 7 does not show a simple relationship between χ 2 and g.For instance, at latitudes between 20 • S and 30 • N, χ 2 of the R14 model, which has the lowest g, is within the range of the other models, and χ 2 of the R0001 and R35 models are not similar even though these two models have similar g values.

Latitudinal Variations of Consistency in Measurements of Thick Homogeneous Clouds
Since low  clouds are less likely to be influenced by 3D effects, the results are recomputed for thick homogeneous clouds ( < 0.4 and τ > 16).Figures 7 and 8 are the same as Figures 3 and  4, but only for thick homogeneous clouds.Compared to the results for global clouds (Figure 3), the results for thick homogeneous clouds show that: (1) the χ of all selected models are significantly lower in thick homogeneous clouds; (2) the absolute differences of χ among models are much larger in both the tropics and extra-tropics; and (3) the χ of each model in the Northern Hemisphere is not significantly larger than χ in the Southern Hemisphere.While the R05 model generally has the lowest χ in all latitudes in Figure 3, the χ and χ results are smaller than the χ counterparts at latitudes between 20° S and 30° N in Figure 7.In latitudes between 20° S-60° S and between 30° N-60° N, the R05 model and the less roughened ice models produce better fits among the six.In other words, the optimal model suggests the need for rougher ice particles in the tropics than in the extra-tropics.Furthermore, Figure 7 does not show a simple relationship between χ and g.For instance, at latitudes between 20° S and 30° N, χ of the R14 model, which has the lowest g, is within the range of the other models, and χ of the R0001 and R35 models are not similar even though these two models have similar g values.Figure 8 is the same as Figure 4 but computed for thick homogeneous clouds ( < 0.4 and  > 16).For SZA < 30°, the  values computed with all six models display nearly the same pattern as shown in Figure 4, and values are close to zero at most measurable scattering angles.Every model in Figure 8 shows positive  values for the lowest measurable Θ in each latitude when 30° < SZA < 50°, particularly with less roughened models.At high latitudes (SZA > 50°), the  values show a relationship with g.Specifically, the  values at the smallest measurable Θ turn negative with a higher g model and at the largest measurable Θ turn positive.For SZA < 30 • , the A diff values computed with all six models display nearly the same pattern as shown in Figure 4, and values are close to zero at most measurable scattering angles.Every model in Figure 8 shows positive A diff values for the lowest measurable Θ in each latitude when 30 • < SZA < 50 • , particularly with less roughened models.At high latitudes (SZA > 50 • ), the A diff values show a relationship with g.Specifically, the A diff values at the smallest measurable Θ turn negative with a higher g model and at the largest measurable Θ turn positive.

Discussion
To examine the latitudinal variations of six ice particle models with a range of surface roughness, we compute spherical albedo difference ( ) values using different LUTs for each model.We then compute the χ (the residual sum of squares of mean  ) in latitude.The  is defined as the difference between the cloud spherical albedo ( ) of a selected MISR camera and the mean of  averaged across all available cameras for a given pixel.The χ values in the extra-tropics are significantly larger than the counterparts in the tropics.Hexagonal column ice particles occur more frequently at low latitudes than high latitudes [31], and this may be a possible reason that the aggregated hexagonal column habit better explains the ice particle in the tropics than in the extratropics.The differences of χ values among ice models in the global data are not significant in the tropics, and MISR data are unable to discriminate between the models.However, in the thick homogeneous cloud regime, the differences of χ values among models were significantly heightened in both the tropics and extra-tropics.The optimal particles are more roughened in the

Discussion
To examine the latitudinal variations of six ice particle models with a range of surface roughness, we compute spherical albedo difference (A diff ) values using different LUTs for each model.We then compute the χ 2 (the residual sum of squares of mean A diff ) in latitude.The A diff is defined as the difference between the cloud spherical albedo ( A s ) of a selected MISR camera and the mean of A s averaged across all available cameras for a given pixel.The χ 2 values in the extra-tropics are significantly larger than the counterparts in the tropics.Hexagonal column ice particles occur more frequently at low latitudes than high latitudes [31], and this may be a possible reason that the aggregated hexagonal column habit better explains the ice particle in the tropics than in the extra-tropics.The differences of χ 2 values among ice models in the global data are not significant in the tropics, and MISR data are unable to discriminate between the models.However, in the thick homogeneous cloud regime, the differences of χ 2 values among models were significantly heightened in both the tropics and extra-tropics.The optimal particles are more roughened in the tropics than in the extra-tropics in thick homogeneous clouds, which is in agreement with previous results [4].The high frequency of strong convection in the tropics likely has an influence on the MISR observations, leading to an indication of heightened surface roughening.
In addition to the χ 2 variation with latitude, the variation of A diff as a function of increasing scattering angle Θ also varies with latitude but more precisely changes with SZA when using any of the LUTs.These differences of A diff also appear in the clouds stratified by cloud homogeneity index H σ and optical thickness τ.If a given camera-retrieved τ is higher than the counterparts associated with the other cameras, a higher A s value would be computed for this one camera than for the other cameras, and this higher A s value leads to a positive A diff value.Therefore, the magnitudes of A diff values at different Θ imply differences between the retrieved τ and the mean τ of all cameras in each Θ.A positive (negative) A diff value broadly means that the retrieved τ is larger (smaller) than mean τ of all cameras.Furthermore, the inference of τ is determined primarily by phase function, and the impact of ice particle phase function on reflectance is important [32,33], particularly for optically thin clouds.Thus, the sign of the A diff value is related to the variation of the phase function with a scattering angle.Note that the A diff method is not sensitive to the absolute value of the phase function, but only to the relative variation of the phase function with scattering angle.This approach, also described in C.-Labonnote et al. [32], is able to explain why the A diff value at Θ. ~150 • in Figure 4a is negative at almost every latitude, but this difference does not appear in other panels in Figure 4 because the phase function of the R0001 model has a significant peak at Θ~150 • .Similarly, this is the reason that the R05 model fits the measurements better than the other models at most latitudes.The phase function of the R05 model fits the measurements better than the other selected models over the MISR available Θ range.
High latitudes tend to have a high SZA, and the 3D radiative effect is one potential explanation for how the A diff appears to change with SZA.Previous studies [15,34,35] investigate the impact of solar geometry on the 3D effect, particularly variations of the retrieved τ with SZA.Since τ can easily be converted to A s , there may be a similar potential impact of the 3D effect on the results.In general, low H σ clouds have uniform cloud tops, which approximately satisfy the plane-parallel approximation.Thus, the impact of the 3D effect on reflectance of this type of cloud is not as strong as on a high H σ cloud.Therefore, low H σ clouds are usually less susceptible to the 3D effect than high H σ clouds.The stratified clouds in Figure 6 have better agreement with this theory.The trend of A diff values with increasing Θ as a function of SZA is significantly larger for high H σ clouds than for low H σ clouds.Note that we consider only one definition of cloud heterogeneity, and a different definition may lead to different results [36].
The microphysical ice particle habit distribution changes with latitude, dynamical environment, and temperature [37][38][39] could impact the trend of A diff values with increasing Θ at different latitudes as well.In the thick homogeneous cloud regime (Figure 8), not all models show different trends of A diff values with Θ in different latitudes compared to global clouds (Figure 4).In addition, high latitudes frequently have ice clouds with a lower cloud top height [40].If we remove clouds with cloud-top temperatures >233K in the data filtering process, some ice clouds at high latitudes are probably removed because the tropopause is lower there.It is interesting to note that for the thick homogeneous cloud regime, the larger differences of χ 2 values between the models enable some discrimination between values of surface roughness.The best-fit particles are more roughened in the tropics than extra-tropics in thick homogeneous clouds.
The limited and variable range of available Θ measurements might be one source of uncertainty for computing the A diff value.Only limited Θ ranges can be sampled by the MISR sensor, but the whole Θ range must be used in retrievals to compute reflectance.It means that there is insufficient data to know how the phase function in the undetectable Θ range changes in each latitude, and to evaluate how large the influence of the phase function in the undetectable Θ are to our results.The range of available Θ values from MISR varies systematically with SZA (Figure 2).McFarlane and Marchand [26] also suggested that the retrieved ice particle habit depends on the minimum observed Θ.Not only the Θ range but also the pixel frequencies in each Θ-latitude bin changes with latitude.Because of varying sun-earth-satellite viewing geometry, it is not straightforward to evaluate these factors on the computed A diff value.
The consistency of six ice models (R0001, R003, R014, R05, R10, and R35 model) is compared to MISR measurements at different latitudes.Our results, based on the fused MODIS-MISR dataset, can be used to estimate the latitudinal (or SZA) bias in retrievals of MC6 operational products because the MC6 operational ice cloud property retrievals assume that the same R05 model applies globally.Our analysis shows that the R05 model works very consistently with global MISR measurements but there may need to be more roughened particles in the tropics.
In the bi-spectral shortwave technique [41], a large asymmetry factor ice particle model causes a higher retrieved τ and lower r eff [42].The asymmetry factors of the R05 and R014 models are the lowest of the six models, and the R35 model has the largest asymmetry factor.As noted earlier, the extreme randomly distorted particle surface facets of the R35 model may increase the amount of forward scattering as an artifact of the numerical modeling.Since the best-fit model is hard to define in the tropics because of the rather small differences in χ 2 values, the ice cloud τ and effective radius could have large uncertainties by assuming the best-fit model from MISR measurements in the retrieval.

Summary and Conclusions
This study examines the latitudinal dependence of ice particle models, based on the MISR-MODIS fused dataset on December solstices from 2012 to 2015.To evaluate the latitudinal consistency of ice particle models with MISR measurements, we apply the SAD method to the fused dataset by using the ice particle habit of the MODIS Collection 6 (MC6) ice model (R05 model) with six different assumptions of ice particle surface roughness: σ 2 = 0.001, 0.03, 0.14, 0.5, 1.0, and 3.5 (referred to as R0001, R003, R014, R05, R10, and R35 models, respectively).
Of the six models used in this study, the R05 model is one of the most consistent in comparison with the MISR measurements globally but has a similar consistency with other models at lower latitudes.The results of computing the residual sum of squares of mean A diff (χ 2 ) show that all six χ 2 are much higher at >30 • N than all other latitudes.For thick homogeneous clouds, the consistency among models is significantly greater in both the tropics and extra-tropics, and the larger differences of χ 2 values between the models enable some discrimination between values of surface roughness.The optimal particle model should have rougher ice particles in the tropics than in the extra-tropics.Specifically, the MC6 model and less roughened models fit thick homogeneous clouds better than more roughened models in the zonal means for latitudes between 20 • S-60 • S and 30 • N-60 • N; for latitudes between 20 • S and 30 • N, more roughened models produce better fits among the six.
Comparisons of computed A di f f values in latitude bands show that the trend of A di f f values with increasing Θ for each ice model changes with latitude, or more precisely, these variations are primarily a function of SZA.This result demonstrates the latitudinal variation of the best-fit model because of the relationship between SZA and latitude (e.g., the SZA in high latitudes is large).The variation of A di f f values with Θ also appears to change with SZA.However, homogenous clouds (low H σ values) do not show the change of A di f f trend with SZA.In addition, for the classified results, the sign of A di f f changes with increasing H σ , but not with increasing τ, as only the magnitudes of the A di f f values increase with increasing τ.
Because MISR provides reflectances over a wider Θ range than MODIS, the collocated MODIS and MISR data/products provide an unprecedented opportunity to identify the optimal ice particle model for retrievals.The variations of the χ 2 value with latitude found in this study suggest that the invariant ice particle shape assumption for retrieving ice particle properties needs to be reconsidered.

Figure 1 .
Figure 1.The phase functions of an aggregate of eight randomly attached hexagonal particles with an effective radius of 30 µm for six different degrees of surface roughness (σ 2 ): 0.001 (R0001), 0.03 (R003), 0.14 (R014), 0.5 (R05), 1.0 (R10), and 3.5 (R35).The asymmetry factor (g) of each ice model is listed in the legend.The Moderate Resolution Imaging Spectroradiometer (MODIS) Collection 6 ice particle model (MC6 model) is the same as the R05 model here.The phase functions are at the wavelength of 0.86 µm.

Figure 1 .
Figure 1.The phase functions of an aggregate of eight randomly attached hexagonal particles with an effective radius of 30 µm for six different degrees of surface roughness (σ 2 ): 0.001 (R0001), 0.03 (R003), 0.14 (R014), 0.5 (R05), 1.0 (R10), and 3.5 (R35).The asymmetry factor (g) of each ice model is listed in the legend.The Moderate Resolution Imaging Spectroradiometer (MODIS) Collection 6 ice particle model (MC6 model) is the same as the R05 model here.The phase functions are at the wavelength of 0.86 µm.

Figure
Figure 2a,b show the MISR camera names that are used and the normalized frequency of occurrences of Θ and latitude based on the four chosen December solstice days during 2012-2015 from the MISR-MODIS fused dataset over oceans.Due to the varying SZA with latitude and differentVZAs of each MISR camera, the cameras provide reflectances at different Θ for a given pixel, and the Θ changes with latitude.The minimum available Θ is 76 • and the maximum available Θ is 172 • in these selected pixels.Because the daytime orbit of the Terra satellite is from north to south with an equator crossing time at 10:30 am, the forward camera measures at smaller Θ than the aft camera in the northern high latitudes.The three cameras primarily take measurements in the side and backward scattering directions.Note that the presence of arc-shaped strips in Figure2ais a result of filtering out the pixels with sun glint.For the same reason, some bins in Figure2bhave no observations.

Figure
Figure 2a,b show the MISR camera names that are used and the normalized frequency of occurrences of Θ and latitude based on the four chosen December solstice days during 2012-2015 from the MISR-MODIS fused dataset over oceans.Due to the varying SZA with latitude and differentVZAs of each MISR camera, the cameras provide reflectances at different Θ for a given pixel, and the Θ changes with latitude.The minimum available Θ is 76° and the maximum available Θ is 172° in these selected pixels.Because the daytime orbit of the Terra satellite is from north to south with an equator crossing time at 10:30 am, the forward camera measures at smaller Θ than the aft camera in the northern high latitudes.The three cameras primarily take measurements in the side and backward scattering directions.Note that the presence of arc-shaped strips in Figure2ais a result of filtering out the pixels with sun glint.For the same reason, some bins in Figure2bhave no observations.

Figure 2 .
Figure 2. Multi-angle Imaging SpectroRadiometer (MISR) camera (a) names and (b) normalized frequency of occurrences as a function of scattering angle and latitude on the four December solstice days (2012-2015) from the MISR-MODIS fused datasets over ocean.(c) The median value of the solar zenith angle (SZA) as a function of latitude on the same dates.

Figure 2 .
Figure 2. Multi-angle Imaging SpectroRadiometer (MISR) camera (a) names and (b) normalized frequency of occurrences as a function of scattering angle and latitude on the four December solstice days (2012-2015) from the MISR-MODIS fused datasets over ocean.(c) The median value of the solar zenith angle (SZA) as a function of latitude on the same dates.

Figure 3 .
Figure 3. (a) The residual sum of squares of mean  value ( ; see Equation (5)) using the R0001, R003, R014, R05, R10, and R35 ice particle models on December solstices from 2012 to 2015.The dotted curve (and top scale) is the median SZA as a function of latitude.(b) The residual sum of squares of mean  value using each corresponding model minus the residual sum of squares of mean  value using the R05 ice particle model as a function of latitude.

Figure 3 .
Figure 3. (a) The residual sum of squares of mean A diff value (χ 2 ; see Equation (5)) using the R0001, R003, R014, R05, R10, and R35 ice particle models on December solstices from 2012 to 2015.The dotted curve (and top scale) is the median SZA as a function of latitude.(b) The residual sum of squares of mean A diff value using each corresponding model minus the residual sum of squares of mean A diff value using the R05 ice particle model as a function of latitude.

Figure 5 .
Figure 5.The variations of median  values with latitude and scattering angle computed from a synthetic dataset generated with the R05 model (a, b), the R10 model (c, d), and the R35 model (e, f), and then retrieved using the LUTs with the other two models.Compared to the simulated dataset produced by the R10 model (Figure 5c,d) and the R35 model (Figure 5e,f),  values computed from the simulated dataset produced by the R05 model (Figure 5a,b) have better consistency with  values computed real data (Figure 4) at high SZAs.Specifically, Figures 5a and 4e illustrate  values based on the same R10 LUTs.The simulated  pattern in Figure 5a reproduces the negative-positive patterns of  values (i.e., positive values of  become negative with increasing Θ in a latitude) at high latitudes in Figure 4e.

Figure 5 .
Figure 5.The variations of median A diff values with latitude and scattering angle computed from a synthetic dataset generated with the R05 model (a, b), the R10 model (c, d), and the R35 model (e, f), and then retrieved using the LUTs with the other two models.

Figure 7 .
Figure 7.The same as Figure 3 but computed for thick homogeneous clouds only.(a) The using the R0001, R003, R014, R05, R10, and R35 ice particle models on December solstices from 2012 to 2015.The dotted curve (and top scale) is the median SZA as a function of latitude.(b) The  value using each corresponding model minus the  value using the R05 ice particle model as a function of latitude.

Figure 7 .
Figure 7.The same as Figure 3 but computed for thick homogeneous clouds only.(a) The χ 2 using the R0001, R003, R014, R05, R10, and R35 ice particle models on December solstices from 2012 to 2015.The dotted curve (and top scale) is the median SZA as a function of latitude.(b) The χ 2 value using each corresponding model minus the χ 2 value using the R05 ice particle model as a function of latitude.

Figure 8
Figure8is the same as Figure4but computed for thick homogeneous clouds (H σ < 0.4 and τ > 16).For SZA < 30 • , the A diff values computed with all six models display nearly the same pattern as shown in Figure4, and values are close to zero at most measurable scattering angles.Every model in Figure8shows positive A diff values for the lowest measurable Θ in each latitude when 30 • < SZA < 50 • , particularly with less roughened models.At high latitudes (SZA > 50 • ), the A diff values show a relationship with g.Specifically, the A diff values at the smallest measurable Θ turn negative with a higher g model and at the largest measurable Θ turn positive.

Figure 8 .
Figure 8.The same as Figure 4 but computed for thick homogeneous clouds only.The A di f f value in latitude-scattering angle bins on December solstices from 2012 to 2015 using (a) R0001, (b) R003, (c) R014, (d) R05, (e) R10, and (f) R35 ice particle models.The dashed and solid lines correspond to SZA 50 • and 30 • , respectively.