Degree of ice particle surface roughness inferred from polarimetric observations

The degree of surface roughness of ice particles within thick, cold ice clouds is inferred from multidirectional, multi-spectral satellite polarimetric observations over oceans, assuming a column-aggregate particle habit. An improved roughness inference scheme is employed that provides a more noise-resilient roughness estimate than the conventional best-fit approach. The improvements include the introduction of a quantitative roughness parameter based on empirical orthogonal function analysis and proper treatment of polarization due to atmospheric scattering above clouds. A global 1-month data sample supports the use of a severely roughened ice habit to simulate the polarized reflectivity associated with ice clouds over ocean. The density distribution of the roughness parameter inferred from the global 1month data sample and further analyses of a few case studies demonstrate the significant variability of ice cloud singlescattering properties. However, the present theoretical results do not agree with observations in the tropics. In the extratropics, the roughness parameter is inferred but 74 % of the sample is out of the expected parameter range. Potential improvements are discussed to enhance the depiction of the natural variability on a global scale.


Introduction
Satellite observations at visible and infrared wavelengths can characterize global cloud microphysical parameters and radiative properties.Numerous techniques have been developed to retrieve ice cloud optical and microphysical properties from radiometric measurements (e.g., Inoue, 1987;Nakajima and King, 1990;Minnis et al., 1993) and have been adopted in operational retrieval efforts (Rolland et al., 2000;Platnick et al., 2003;Minnis et al., 2011).A synergetic combination of satellite and in situ observations (e.g., Heymsfield et al., 2002Heymsfield et al., , 2013) ) serves as a constraint for the parameterization of bulk ice cloud optical properties for remotesensing implementations as well as general circulation models (GCMs).
The accuracy of these retrieval techniques and the validity of downstream applications including GCM radiation parameterization hinges on steady improvements in the singlescattering calculations involving ice crystals.As ice clouds consist of nonspherical particles with characteristic sizes much larger than the wavelengths of interest, the singlescattering properties depend on the size, shape, and microscopic morphology of the particles (Macke et al., 1996;Yang et al., 2008a;Xie et al., 2009;Baum et al., 2010;Um andMcFarquhar, 2007, 2009; Ulanowski et al., 2006Ulanowski et al., , 2014)).In the solar shortwave spectrum, particle shape, surface texture, and crystal imperfections have a substantial influence on the single-scattering properties.Recent improvements in scatter-ing calculation techniques are being incorporated into models that represent diverse ice particle populations in clouds.
However, it is challenging to quantify some of these influential microphysical parameters, given current satellite sensors.As a result, little information about their variability is available.The discrepancies among climate models (Waliser et al., 2009) in terms of ice water path (IWP) indicate that the GCM parameterizations need more reliable constraints on IWP.Recent work by Sourdeval et al. (2015), which includes direct retrieval of IWP, is a novel approach to this problem.
The ability of passive and lidar sensors to correctly infer IWP requires knowledge of ice cloud radiative properties.Application of an unrealistic ice model, e.g., with only smooth (unroughened) surfaces, results in an overall global bias (Macke and Mishchenko, 1996;Yang et al., 2007Yang et al., , 2008b;;Holz et al., 2016) as well as seasonal biases (Zhang et al., 2009) in cloud property retrievals.The overarching goal of this paper is to gain a better understanding of the constraints in the microphysical parameters of global ice clouds using angular polarimetric observations and state-of-the-art light-scattering computational capabilities.
Multidirectional polarimetric observations can constrain the representative particle shape and surface texture condition (specifically, the degree of surface roughness), owing to the sensitivity of the polarization state of reflected light to small-scale particle structures.These measurements have been used to infer both particle habit (Chepfer et al., 1998;C.-Labonnote et al., 2001;Masuda et al., 2002;Knap et al., 2005;Baran and C.-Labonnote, 2007) and surface roughness (Baran and C.-Labonnote, 2006;Cole et al., 2013Cole et al., , 2014)).Since polarized reflectivity saturates at relatively small optical thickness (generally about τ = 5, Masuda and Takashima, 1992), the conventional "best-fit" approach to this problem computes the residual sum of squares (RSS) from the multi-angle observations of polarized reflectivity and reflectivity simulations, and selects the ice particle model that minimizes the RSS when τ > 5.
The previous studies imply that the use of roughened particles is necessary to achieve maximum consistency between observations and numerical scattering calculations.Furthermore, the spectral consistency of visible/near-infrared and thermal infrared retrievals (Baran and Francis, 2004) was recently investigated by Liu et al. (2014) and Holz et al. (2016), who report that retrieved ice cloud optical thicknesses are more consistent when particles are roughened.
The treatment of particle surface roughness here is not a rigorous approach.Rather, it is an approximation of the effects of roughened surface texture (Neshyba et al., 2013) and other kinds of imperfections present in natural ice cloud particles.The scattering properties calculated by this approximate method are in reasonable agreement with those calculated by rigorous ray-tracing methods (Yang et al., 2008a).Although previous studies suggest that some degree of roughness is desirable, the issue remains as to the amount of roughness that should be adopted for global satellite-based retrievals, or used in numerical models.
Recent work by van Diedenhoven et al. (2012van Diedenhoven et al. ( , 2014) ) simultaneously infers both the aspect ratio and the degree of roughness from a combination of polarimetric and intensity observations over a virtually continuous parameter space, assuming that simple hexagonal ice particles can explain observations.The ability to infer a representative ice cloud particle aspect ratio adds yet another dimension to the problem.Such exploration into the variability of ice particle microphysical properties can lead to a more reliable satellite climatology of ice clouds.This study focuses on the quantitative inference of ice particle roughness parameter for a specific particle habit, and will not include a detailed investigation of aspect ratio.
While a conventional "best-fit" approach can constrain the range of the average roughness parameter at the global scale, it is not suitable for pixel-by-pixel inferences.This is because the signal-to-noise ratio for particle roughness is low, and in the conventional "best-fit" approach, even random observational errors can modify the inferred histogram significantly when it is applied to individual pixels.Figure 1 illustrates how such a modification takes place if the method is applied to a synthetic signal with random noise.To produce Fig. 1, viewing geometries are extracted from 1 month (September 2005) of cloud observations by the POLarization and Directionality of the Earth's Reflectance (POLDER) sensor (Deschamps et al., 1994) onboard the Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observation from Lidar (PARASOL) satellite (Fougnie et al., 2007).The "best-fit" inference is applied to synthetic multi-angle cloud polarized reflectivities (L p , defined in Sect.2.1) with and without random noise.In synthesizing the signal, a column aggregate particle shape (e.g.Yang et al., 2013) is assumed with a roughness parameter of σ 2 = 0.15 (variance of the slope of random facet tilts, see Yang et al., 2008b, for details), and the random error has a normal distribution with variance equivalent to the POLDER observational error, which is estimated in Sect.2.1.The hatched bar is the histogram with noise and the gray bar is that without noise.Note that the distinct peak at σ 2 = 0.15 is no longer apparent when instrumental noise is included, indicating the necessity of appropriate treatment of the error distribution in the analysis.
This paper demonstrates how a continuous parameter space for the roughness retrieval is constructed and how it can be used to infer the particle roughness of optically thick ice clouds on a pixel-by-pixel basis.Section 2 provides the details of the data and inversion method we employ, and the result of the application to 1 month of global data is described in Sect.3. Concluding remarks are given in Sect. 4.

Methodologies
To establish a method resilient to observational error, we first examine random errors in POLDER data and select pixels based on the MODIS Collection 6 cloud product, as given in Sect.2.1.Then, a continuous parameter space for inferring roughness is constructed by using an empirical orthogonal function analysis, and used in the retrieval scheme with the maximum likelihood method.The construction of the parameter space, and the design and the performance of the forward model, are discussed in Sect.2.2.

Reflectivity from POLDER
The POLDER sensor onboard the PARASOL satellite provides multispectral polarimetric observations at up to 16 viewing geometries for a single overpass (Fougnie et al., 2007).The PARASOL satellite was in the A-train satellite constellation from 2004 to 2009 and continued operation in a separate orbit until late 2013, providing a total of 9 years of global polarimetric observation data.The design of the instrument is inherited from previous POLDER sensors on the ADEOS (ADvanced Earth Observing Satellite) platforms.POLDER sensors provide the first three elements of the Stokes vector from three images taken successively with linear polarization filters (Deschamps et al., 1994).This study uses the single-pixel data set in the PARA-SOL Level 1B product.The approximate resolution is 6 km × 6 km.
PARASOL products report the intensity of reflection in terms of normalized radiance L n , which is equal to the reflectivity R of the surface-atmosphere system multiplied by the factor µ 0 = cos θ 0 (cosine of solar zenith angle).
In a similar manner, the polarized reflectivity is reported in terms of normalized radiance, so (L n QU ) becomes the first three Stokes parameters.In other words, the normalized polarized radiance L np = Q 2 + U 2 is equal to the polarized reflectivity R p multiplied by µ 0 .
where Q i and U i are defined to form the first three Stokes parameters in terms of radiance (I Q i U i ).It is worth noting the similarity between Eqs.
(1) and (3).We conduct the analysis in terms of L np = µ 0 R p defined in Eq. (3) to simplify the error estimate.
The distribution of random errors in L np observed with POLDER is estimated in the following procedure.A reflection property of an optically thick ice cloud is that the modified polarized reflectivity L nmp = η (µ + µ 0 ) L np /µ 0 (where η = ±1, C.- Labonnote et al., 2001) crosses zero at scattering angle ≈ 170 • as shown in Fig. 2.This implies that the polarization signal at ≈ 170 • is primarily due to the observational noise with additional contributions from the variation of cloud particle scattering properties.We utilize this reflection property to estimate the magnitude of observational noise from the POLDER data at scattering angles between 168 and 172 • , and further estimate the noise level at other angles with a typical polarization state of cloud reflection.
The POLDER observational noise consists of radiometric noise and misregistration noise.The misregistration noise is inherent in the POLDER sensor's design that extracts polarimetric information from three images successively taken with different polarizers.The co-registration process of these three images is an inevitable source of error.As the distribution of misregistration noise is unknown, our instrument model attempts to explain both noise components with a radiometric noise model in the following analysis.
We define a random variable L np that serves as a statistical model of observed L np as follows.
where random variables X 1 , X 2 , and X 3 represent the radiances of a pixel in the original three images with different polarizers (not available in a product).With the statistical model outlined in Eq. ( 4), we first assume that X 1 , X 2 , and X 3 follow the same normal distribution centered at 0.5 with variance s 2 (i.e., X i ∼ N (0.5, s 2 )) because the expectation of polarized radiance L np is assumed to be zero at scattering angles between 168 and 172 • .With this assumption, we apply the parametric bootstrap method (e.g., Evans and Rosenthal, 2010) to obtain the distribution of L np as a function of variance s 2 .The observational distribution of L np at 0.865 µm in the scattering angles between 168 and 172 • (within the rectangular box in Fig. 2) is shown in the bar chart of Fig. 3, and compared with the theoretical distribution with s = 0.00095 (solid line).Figure 4 justifies our selection of s = 0.00095 by showing that the sum of squared errors of the density in each bin of the histogram (Fig. 3) is minimized when s = 0.00095.Therefore, we take s = 0.00095 as the standard error for X 1 , X 2 , and X 3 .In Fig. 3, the distribution from observations is slightly more skewed than the distribution from bootstrapping, but their agreement justifies the use of the simple statistical model formulated in Eq. ( 4) to quantify the magnitude of measurement errors.
To obtain the approximate magnitude of the L np error at other scattering angles, the same parametric bootstrap method is applied with the degree of linear polarization fixed at 5 %, which is the upper limit for typical ice cloud reflection.This selection does not significantly affect the following analysis.When the signal is polarized, random variables X 1 , X 2 , and X 3 do not follow the same distribution, but we assume that the standard errors for X 1 , X 2 , and X 3 still stay the same.Figure 5 shows the estimated magnitude of error (variance) as a function of normalized radiance L n .The variance of L np asymptotes to a near-constant value once L n reaches L n = 0.2.As shown in insets, the distribution becomes closer to a normal distribution with increasing L n .Based on the discussion above, we conclude that the error distribution of L np approximately follows a normal distribution with variance var(L np ) = 1.35 × 10 −6 for a reflective target (L n ≥ 0.2).This estimate of error is about the same  magnitude as the value by Fougnie et al. (2007).Note that we assume that the error is purely from observational noise, neglecting any natural cloud variability.Therefore, the actual radiometric noise level should be somewhat smaller than our estimate.We estimate the magnitude of error using the 0.865 µm channel because the channel is likely to be the least contaminated by other sources of uncertainty such as ozone absorption (0.67 µm) and Rayleigh scattering (0.49, 0.67 µm).We apply the same variance to all three POLDER channels used in the analysis (0.865, 0.67, and 0.49 µm).

Ancillary data from MODIS and AIRS
The moderate resolution imaging spectroradiometer (MODIS) instruments onboard the Aqua and Terra satellites measure radiance at multiple visible and infrared wavelengths, providing various products (King et al., 2003) that are complementary to those from PARASOL.Of interest here is the Collection 6 Level 2 cloud product (MYD06) from Aqua MODIS, with which the PARASOL satellite was flying in formation until 2009.Cloud top temperature  and thermodynamic phase are extracted from MYD06 and are collocated to POLDER data to be used in the analysis described later in this section.In addition, Level 3 monthly mean ozone concentration from the Atmospheric Infrared Sounder (AIRS) on the Aqua satellite is also used, in particular to account for absorption by ozone that attenuates reflected radiation in the visible range.

Collocation and selection
The PARASOL Level 1 radiometric data are first collocated with the MODIS Level 2 cloud product (Platnick et al., 2015) to select pixels containing ice clouds.Only PARASOL pixels that have corresponding MODIS observations are selected, and filtered by the criteria summarized in Table 1.The intent of the filtering process is to avoid cloud edge contamination, to avoid supercooled water droplets, and to select pixels where clouds are optically thick.The selection criterion of 208 K is a threshold used to identify convective precipitation in the tropics (Mapes and Houze, 1993).The analysis is applied only over oceans so the influence of surface reflection is minimal.
A "pixel" in the PARASOL Level 1 product contains reflectivity data observed from up to 16 viewing angles.An in-dividual reflectivity value stored in a pixel is called a "view", and we select valid views using criteria relating to scattering angle and sunglint angle (see Table 1).When five or more valid views are contained in a pixel that satisfy all pixel criteria previously mentioned, the pixel is marked as valid, and the roughness inference is attempted.

Selection of retrieval parameters
To overcome the problem of the conventional "best-fit" approach that uses a discrete set of roughness parameters, we construct a continuous parameter space for the particle roughness with empirical orthogonal functions (EOF).The goal of the EOF analysis is to find the parameter space that describes the variation of the −P 12 element of the phase matrix when varying the particle roughness.An ideal approach would be to use a collection of −P 12 values from observations (Rodgers, 2000), but such a dataset is unavailable.For this reason, we apply EOF analysis to the −P 12 simulated with light-scattering calculations.The first and second EOFs together explain 99.3 % of the entire variation of −P 12 in the scattering angle range from 60 to 160 • .This implies that the following approximation is valid in the scattering angle range 60 where is scattering angle, Q 1 ( ) and Q 2 ( ) are the first and second EOFs, and x 1 σ 2 and x 2 σ 2 are weights for EOFs (EOF scores).Note that the set of EOFs and EOF scores obtained in this way depends on the selection of particle shapes and the degree of roughness.In our EOF analysis, 10 prescribed roughness parameter (σ 2 ) values are used: 0, 0.03, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, and 0.7.These roughness parameters are selected to outline the variation of −P 12 over the course of roughness changes, including the roughness parameter used in MODIS Collection 6 (σ 2 = 0.5).The EOF scores are shown in Fig. 6.The EOF 1 primarily describes the degree of roughness, and the EOF 2 score has sensitivity to large roughness parameters.Therefore, the EOF 1 and EOF 2 scores are selected as retrieval parameters related to particle roughness.Since the phase matrix follows a linear mixing rule, −P 12 of a mixture containing multiple degrees of roughness is also approximated by Eq. (5).For example, a mixture of MODIS Collection 6 particle (σ 2 = 0.5) and moderately roughened particle (σ 2 = 0.03) produces EOF scores (x 1 x 2 ) on a straight line between (x 1 (0.5), x 2 (0.5)) and (x 1 (0.03), x 2 (0.03)).Constructing a continuous parameter space using EOF scores (x 1 x 2 ) is a powerful approach because the method guarantees that the parameter space contains any mixture of prescribed shapes or degree of roughness.In Sect.2.2.1, we show that these EOF scores can be used to accurately parameterize the normalized polarized reflectance L np for a given direction.
The details of the particle model are as follows.The EOF analysis is applied to −P 12 elements of the phase matrices calculated by the method described by Yang et al. (2013), which is a combination of the improved geometric optics method (IGOM, Yang and Liou, 1996) and the Amsterdam discrete dipole approximation method (ADDA, Yurkin et al., 2007).Surface roughness is applied only in the IGOM computation (D max > 10 µm).The column aggregate shape is chosen because the most extensive previous study on a global scale (Cole et al., 2014) implies that this habit produces the most consistent agreement with observations.In addition, this habit is used in the retrieval scheme for the operational MODIS Collection 6 cloud products.This particle shape is an aggregate of eight column elements that are solid hexagonal particles with slightly different particle aspect ratios (originally defined in Yang and Liou, 1996, see Yang et al., 2013, for geometric parameters).A gamma particle size distribution with an effective size (diameter) of 60 µm and an effective variance of 0.1 is used in this study since we expect little impact on our analysis due to this size distribution selection (Cole et al., 2014).
As the EOF 1 score is a monotonic function of the roughness parameter and explains most of the −P 12 variation (85.6 %), it can be considered as an effective roughness pa- Figure The relationship between the particle roughness parameter and the EOF1 score.The natural logarithm of the particle roughness parameter is nearly linear to the EOF 1 score.This implies that the particle roughness can be directly inferred from the EOF 1 score.
rameter for this shape.The relation between EOF 1 scores and the natural logarithm of roughness parameters is nearly linear (Fig. 7), indicating that the roughness parameter can be subsequently inferred after the inference of the EOF 1 score.The straight line in Fig. 7 is the regression line defined in the form As the roughness parameter computed from Eq. ( 6) does not account for the variation of EOF 2 score, it is inaccurate for the mixture of particles containing multiple degrees of roughness.Equation ( 6) is introduced to compare our retrievals to the conventional discrete parameter space.
Another factor that can impact the roughness retrieval is atmospheric Rayleigh scattering above the cloud.Abovecloud Rayleigh scattering has previously been used to infer cloud top pressure from polarimetric measurements (e.g., Buriez et al., 1997), with results comparable to those from O 2 A band retrievals and ISCCP (Parol et al., 1999).With the POLDER instruments, Rayleigh scattering is primarily detected as a spectral and directional difference of polarized reflectivities.Figure 8 shows the change of L np at 0.865 µm, as a function of scattering angle, in response to a 300 hPa change in cloud top pressure (i.e., from 200 to 500 hPa, the red line) and a change of similar effect in roughness parameter (from σ 2 = 0.15 to 0.5, the dashed green line).The effects of cloud top pressure and roughness parameter changes on L np have different directional patterns but comparable magnitudes.The variation of the cloud top height must therefore be well constrained or retrieved simultaneously when attempting to infer the roughness parameter.

Construction of forward model
Once the inverse problem is formulated, the next step is to construct a forward model that is fast enough to be em- bedded in the inversion algorithm.From the discussion in Sects.2.1 and 2.2.1, the inverse problem is formalized as follows: (1) the parameters to be inverted are the EOF 1 and EOF 2 scores, and cloud top pressure; and (2) observations are MODIS-AIRS-collocated L np from POLDER at central channel wavelengths 0.865, 0.67, and 0.49 µm.To satisfy the requirements for numerical efficiency, the present forward model is based on a look up table.The adding-doubling radiative transfer program is used to compute L np for every phase matrix with seven atmospheric scattering optical thicknesses above the cloud: 0, 0.02, 0.05, 0.1, 0.15, 0.2, and 0.3.The result at a specific viewing geometry (denoted by subscript i) and an optical thickness (denoted by subscript j ) is parameterized by simple linear regression model defined as where x 1 and x 2 are EOF scores obtained in Sect.2.2.1, and a (i,j ) , b (i,j ) , and c (i,j ) are regression coefficients determined by the polarized reflectivities for multiple phase matrices.The viewing geometry is gridded as follows: solar zenith angles from 0 to 81 • , viewing zenith angles from 0 to 75 • , and relative azimuth angles from 0 to 180 • , with an interval of 3 • for each.The regression is repeated for seven atmospheric scattering optical thicknesses above the cloud and more than 40 000 viewing geometries.With this fast forward model, once cloud top height and EOF scores are given, L np can be obtained for each specific viewing geometry and wavelength.Note that when a set of EOF scores (x 1 x 2 ) is not exactly at the values corresponding to the 10 prescribed phase matrices, the forward model linearly interpolates the polarized reflectivity.We confirmed that the interpolation usually produces a reliable polarized reflectivity simulation for a phase matrix of intermediate roughness and a mixture of phase matrices.
The fast model constructed in this way is accurate enough to solve our inverse problem.A typical difference between an exact calculation and our forward model is shown in Fig. 9.
. Difference in L np between exact radiative transfer calculations and our simplified forward model.At almost all angles, the difference is less than 1 × 10 −4 .The polar plot shows the distribution of bias when the particle roughness parameter is σ 2 = 0.15.The bias is a function of scattering angle.However, the magnitude of error is acceptably small compared to the random observational error.
The overall accuracy is within 1 × 10 −4 in terms of L np and the peak-to-peak variation is 5 × 10 −4 even in the worst case (σ 2 = 0.03).The overall error of 1 × 10 −4 implies that the model bias is less than 10 % of the observation error given by ( var(L np ) = 1.35 × 10 −6 = 1.16 × 10 −3 ).The bias may be detected in the residual of the inversion, but the influence on the roughness inference is negligible.
In calculating cloud reflectivity, a single-layer homogeneous cloud is assumed, and the cloud optical thickness is set to 5 (roughly the saturation point of polarized reflectance).No aerosol is assumed to be present above and below clouds.As optically thick cold ice clouds occur in the upper troposphere, the radiometric contribution from lower tropospheric aerosols is neglected.For the same reason, the surface is assumed to be dark.There may be an influence from aerosols above the cloud layer, such as transported mineral dust and stratospheric sulfates, but we disregard them to be in line with previous studies.The influence of such aerosol layers on inferences of cloud properties is beyond the scope of this paper but should be investigated in the future.
The adding-doubling radiative transfer program formulated by de Haan et al. (1987) with significant improvements by Huang et al. (2015) is used in the calculation.The first-order scattering is calculated analytically and combined with the multiple scattering results from the addingdoubling model, following the TMS method (Nakajima and Tanaka, 1988).Further, the cloud reflectivity is multiplied by the transmissivity that changes due to ozone absorption; the transmissivity is calculated from the monthly mean AIRS ozone concentration.= 0.14.

Maximum likelihood estimation
Once the inverse problem is formulated and the forward model is built, the last step is to find the set of parameters for each pixel based on observations.The simple but powerful maximum likelihood method with a normal error distribution is appropriate for our problem because we have little knowledge about the distribution of parameters (EOF scores).As each pixel typically contains five to eight valid views (Sect.2.1.2) at three wavelengths, the number of observations in each pixel ranges from 15 to 24.
The standard deviation (SD) and correlation (Corr) of inferred parameters are calculated in the framework of maximum likelihood estimation, and used to avoid underconstrained inferences.The pixel is rejected if SD (EOF 1 Score) > 0.02, SD (EOF 2 Score) > 0.02, or Corr (EOF 1 Score, EOF 2 Score) > 0.3.The standard deviation and the correlation depend strongly on the observation geometry and particle model and are almost independent of the observed polarized reflectivity.Therefore, this rejection process can be interpreted as the refinement of pixels based on the information content to achieve a reliable inference.
The error distribution is confirmed to be normal (see Sect. 2.1.1),so if the problem is not strongly nonlinear, the parameters' error distributions are expected to be normal as well (Rodgers, 2000).As expected, Fig. 10 demonstrates that the application of the maximum likelihood method with synthetic L np data results in a symmetric distribution about the EOF 1 score corresponding to the true roughness parameter σ 2 = 0.15.The distribution is not strictly normal because the number of observations in each pixel varies, but the error distribution of each pixel is theoretically derivable, as well as the confidence interval.
For the synthetic retrieval in Fig. 10, the median of the inverted EOF 1 score is −0.00336 and the corresponding roughness parameter is σ 2 = 0.14.The interquartile range of the EOF 1 score distribution is [−0.01146: 0.00476], which corresponds to the roughness parameter range of [0.05 : 0.36].The result indicates that our approach has a practical skill in estimating the particle roughness parameter out of observations superimposed with noise.This is a remarkable contrast with the traditional "best-fit" approach (cf.Fig. 1).
The distribution of the χ 2 values for the synthetic retrieval is presented in Fig. 11.The χ 2 value is a variance-normalized residual squared sum that is defined for each pixel, and follows the χ 2 distribution with degrees of freedom N d if the inversion is successful, where N d is the observational degrees of freedom (approximately, the number of observations in a pixel).As the χ 2 distribution of N d degrees of freedom has a peak about N d , the distribution of the χ 2 value indicates whether the inversion is successful.If the location of the peak of a distribution of χ 2 values is smaller than N d , the observation error may be overestimated, and if the location of the peak is larger than N d , the observation error is underestimated, or the forward model does not represent reality (Rodgers, 2000).The distribution in Fig. 11 has a peak at about 12, and very few pixels have a χ 2 value over 40.This is a reasonable distribution because the number of observations (≈ N d ) is about 15 to 24 for most pixels.Because the 95th percentile for the χ 2 distribution with 24 • of freedom is 36.42, it is no surprise that very few pixels have a χ 2 value over 40.
Figures 10 and 11 demonstrate the validity of our inference framework under an idealized situation, where the error distribution and the true roughness parameter are constant.In application to actual satellite data, however, the true roughness parameter varies from pixel to pixel while the error distribution stays the same.Therefore, the distribution of the EOF 1 score must be more spread out as a result of convolution of the error distribution and the true roughness parameter distribution.In contrast, the χ 2 distribution is expected to be about the same.The result of the application to actual data is given in the next section.
3 Results and discussion

Roughness parameter of cold ice cloud over oceans
With the cloud selection criteria listed in Table 1, 79 192 pixels based on 1 month of collocated PARASOL/MODIS data over oceans during September 2005 were selected for inversion.The information content was sufficient for full analysis of 23 359 pixels, for which results are presented in this section.
The histogram of the inferred EOF 1 score is presented in Fig. 12 for the extratropical (latitude > 30 • ) oceans.The width of the histogram in Fig. 12 is broader than the monodispersive roughness case (Fig. 10), indicating significant variability in the microphysical properties of clouds.The median of the distribution is −0.0293, corresponding to a surface roughness parameter of 2.82.The interquartile range of the EOF 1 score is [−0.0429: −0.0165], implying 50 % of the data are within the roughness parameter (σ 2 ) range of [0.65 : 13.6].The result supports the use of the roughened particle model in extratropical ice cloud retrievals as suggested by previous studies.While our analysis is limited to very cold ice clouds over ocean, the validity of using roughened crystals in the MODIS Collection 6 ice model is supported, although further explorations into warmer and optically thinner clouds are desirable.In general, cloud particles become more complex as the cloud temperature increases (Heymsfield et al., 2002), thus we expect more roughened particles in warmer clouds that are not included in our analysis.The distributions of the χ 2 value in the tropics and extratropics are separately presented in Fig. 13.As discussed in the previous section, the distribution of χ 2 values indicates the validity of the inversion.While the distribution of the χ 2 values in the extratropics shows reasonable behavior (Fig. 13a), the distribution of the χ 2 values in the tropics has a very long tail with the mean χ 2 being 59.7, which is unacceptably large (Fig. 13b).This long tail implies that our forward model does not properly reproduce the observed L np field in the tropics, presumably because some underlying assumptions are not appropriate or the information content is not enough.Some possibilities that violate our underlying assumptions include sub-pixel scale cloud heterogeneity, the presence of ice particles with other habits or aspect ratios, their vertical heterogeneity, cloud 3-D effects, and the effect of aerosols.

Unexpectedly large roughness values in the extratropics
As the roughness parameter of 2.82 lies outside of our prescribed roughness parameter range (0 to 0.7), it is an estimate by extrapolation.Yet, this projection of roughness parameter implies that the conventional degree of roughness may not be sufficient to represent actual cloud particles with the aggregate column model.The proportion of pixels that contains inferred roughness parameter σ 2 > 0.7 is 74 %, which also indicates the limit of this particle shape.As the accuracy of roughness approximation for such a large roughness parameter is questionable, a particle shape that can fit observations with less intense roughening may be suitable for the representation of natural clouds.To attribute the cause of unphysically large roughness value in the extratropics, the same retrieval process is repeated assuming three additional particle shapes.Figure 14a shows the original inference with aggregate of columns shape, in which the observation density peaks away from the line connecting 10 points that corresponds to prescribed roughness values.The aggregate of plates (Fig. 14d) performs worst among the tested particles, and the solid bullet rosette shape (Fig. 14c) shows the largest overlap of parameter space and observation density.These results indicate that the roughness retrieval is sensitive to an assumed particle shape.
We also investigated the contamination by multi-layer clouds and aerosol above clouds by collocating the cloudaerosol lidar with orthogonal polarization (CALIOP) vertical feature mask and cloud layer products.As September 2005, which is analyzed in this study, is before the launch of the CALIPSO satellite, we analyzed the collocated POLDER-MODIS-CALIOP dataset in September 2006 in the extratropics.According to the CALIOP vertical feature mask, on the CALIOP track, about 20 % of pixels that are colder than the brightness temperature threshold of 233 K are possibly contaminated by either multi-layer cloud, aerosol above clouds, or a stratospheric feature.However, the distribution of the retrieved EOF scores is approximately the same even when assuring the absence of aerosol above cloud and limiting the analysis to single-layer clouds (Fig. 15).Therefore, we do not consider that aerosol contamination and multi-Figure 15.The distribution of retrieved EOF 1 and EOF 2 scores when using CALIOP data to filter out clouds with multiple layers or with aerosols above the cloud.The observation frequency is color shaded, and the EOF scores for column aggregate particles (circles) are connected by a line.This analysis is conducted on a different EOF space from Fig. 6.The minimum degree of roughness is σ 2 = 0.03 and the maximum is σ 2 = 1.0.To exclude optically thin clouds, pixels are selected if the CALIOP vertical feature mask product marks total attenuation above ground.No temperature threshold is applied.layer clouds introduce a large bias that brings our estimate out of the range of prescribed parameters.Removal of the multi-layer clouds helps to reduce the number of pixels with very large χ 2 values.

Inference failure in tropics
To gain a better insight into the cause of the long tail in the tropics, a case study is conducted for two cloudy scenes: a typical extratropical scene and a tropical cloud scene with systematically large χ 2 values.Figure 16 displays true color composites from PARASOL with markers indicating the locations of detailed analysis.A green circle is shown where the χ 2 value is less than the 95th percentile of the χ 2 distribution (reasonable deviation from the forward model), and a magenta cross is shown where the χ 2 value is more than the 95th percentile (too far from the forward model).The locations of the magenta crosses in Fig. 16a (typical extratropics) are somewhat systematic; they appear at cloud boundaries or at isolated locations.This may suggest that cloud heterogeneity and cloud 3-D effects cause a small number of inference failures in the extratropics.
In contrast to the typical extratropical scene, magenta crosses are prevalent throughout the tropical scene in Fig. 16b.Since the cloud reflectivity is comparable to the typical extratropical scene, it is not likely that the inference failures are due to contamination by surface reflection.Also failures cannot be fully explained by 3-D effects of clouds as a few green circles appear randomly.Flaws in the assumptions that depend little on the relative location in a cloud, such as cloud particle shape and cloud heterogeneity (e.g.Oreopoulos et al., 2009), or the lack of information content due to the limited scattering angle range are therefore suspected as causes of the inference failure in the tropics.
A close investigation into the correlation of EOF 1 and EOF 2 scores supports the hypothesis that the information content is limited.Figure 17 shows the coefficient of correlation between retrieved EOF 1 and EOF 2 scores.While the distribution is centered at 0 in the extratropics, it peaks at −0.8 in the tropics, indicating limited information content to constrain the parameter space.Also, a validation of retrieved cloud top height using the CALIOP data indicated that the cloud top heights are not properly retrieved in the tropics.The insufficient information content for roughness and cloud top height retrievals is presumably caused by that sampled scattering angles are concentrated near the backscattering direction and zenith angles are small.The directional distribution of the polarized reflectivity is not well captured to constrain the degree of roughness and the spectral contrast of the Rayleigh scattering signal is too weak to infer cloud top height accurately. in Fig. 18 shows the interquartile range of the reconstructed −P 12 which indicates that 50 % of our extratropical observations fall within the shade at a given scattering angle.The blue line is −P 12 for the particle shape used in MODIS Collection 6, and the green line is that for the shape in MODIS Collection 5.Both particle models assume a gamma distribution with effective particle size of 60 µm and effective variance of 0.1.The blue line (Collection 6) is closer to our reconstruction, while the green line (Collection 5) significantly deviates from our reconstruction.This result indicates that the particle habit adopted for MODIS Collection 6 is more consistent with polarimetric observations than the habit mixture used for MODIS Collection 5, for which only one of the habits includes a limited degree of roughness.

Comparison
The reconstructed −P 12 shows stronger side scattering between 80 and 120 • than the MODIS Collection 6 particle model.As the increasing roughness enhances side scattering, weak side scattering of the column aggregate shape may be responsible for the unexpectedly large roughness parameter in the extratropical inferences.By using a shape that has stronger side scattering, it is likely that the degree of roughness that is needed to explain the observations becomes smaller.An example of such a habit mixture is shown by the thick magenta line in Fig. 18.A mixture of two habits (70 % column aggregate particles with roughness parameter of σ 2 = 0.8 and 30 % severely roughened hollow bullet rosette particles with σ 2 = 0.5) included in the scattering property library by Yang et al. (2013), results in a phase function with strong side scattering.

Summary and future directions
In this study, the particle roughness parameter of very cold ice clouds over ocean is inferred by employing a new frame-work that is resilient to the observational error.The distinct feature of the framework is the continuous parameter space that is constructed with an empirical orthogonal function (EOF) analysis.Two EOFs are found to be sufficient to explain the variation of −P 12 with a changing particle roughness parameter, substantially reducing the number of parameters for the forward model.
From unpolarized cloud reflection at a scattering angle of 170 • , the observational error of the PARASOL data is empirically estimated.Supported by the error analysis with parametric bootstrapping, the maximum likelihood method is applied to the inverse problem.The method provides error estimates and correlations for inverted parameters, which are unavailable with the "best-fit" approach used in the previous studies.To correctly incorporate the effect of atmospheric Rayleigh scattering, the cloud-top height is inferred simultaneously.
The application of the present method to cold ice clouds over extratropical oceans results in a roughness parameter of 2.82, implying that the use of the roughened particle model is suitable for cloud property retrievals.By contrasting the distribution of χ 2 values in the tropics and extratropics, we find that the performance of our method needs to be enhanced in the tropics.Possible future technical improvements may be an extension of parameter space to include multiple particle shapes, application to optically thin clouds, and integration with unpolarized radiance observations.The reconstructed −P 12 curve shows better consistency with −P 12 from the particle shape model used in MODIS Collection 6 than −P 12 from MODIS Collection 5.The addition of roughness and a hollow bullet rosette particle shape to the MODIS Collection 6 model further improves the consistency.
Since its launch in 2004, the PARASOL satellite observed global polarimetric reflectivity nearly simultaneously with MODIS for 5 years until leaving the A-train constellation in 2009.A large amount of PARASOL data are available to apply the framework described in this paper.Local variations of the roughness parameter, correlation of the roughness parameter to other meteorological data, and the impact of cloud heterogeneity are to be investigated in our future study.

Data availability
The satellite datasets are available through ICARE Data and Service center, NASA LAADS system, and NASA GSFC GES DAAC.The single scattering property dataset used in this study is available from the author upon request.

Figure 1 .
Figure1.The response of the conventional "best-fit" approach to a synthetic signal with and without random measurement noise.The addition of noise to the synthetic signal results in a distribution of the roughness parameter (hatched bars), from which the true roughness cannot be inferred.This figure is to be compared to Fig.11.

Figure 2 .
Figure 2. Observation density of modified polarized reflectivity (L nmp ) over the Western Pacific during September 2005.L nmp crosses zero at a scattering angle of approximately 170 • .The data in the rectangular box is used to derive the histogram in Fig. 3.

Figure 3 .
Figure 3. Histogram of observed normalized polarized radiance (L np ) from the data in the rectangular box in Fig. 2. The solid line is the simulated error using a parametric bootstrapping method with s = 0.00095.The agreement is sufficient for estimating the noise level.

Figure 4 .
Figure 4. Sum of squared error as a function of standard error (s) of the original sensor noise.The minimum error is achieved when s = 0.00095.

Figure 5 .
Figure 5.The simulated variance of L np as a function of L n .The variance of L np increases as the normalized radiance L n (brightness of a pixel) increases, becoming nearly constant at var(L np ) = 1.35 × 10 −6 once L n reaches L n = 0.2.Insets show that the distribution of L np tends to a normal distribution, justifying the use of a normal distribution as an error distribution of L np for a reflective cloudy pixel.

Figure 6 .
Figure6.The pairs of EOF scores needed to reconstruct the original −P 12 .The EOF 1 score is a monotonic function of particle roughness parameter σ 2 .The EOF 2 score reaches a minimum at particle roughness parameter of σ 2 = 0.1. 

L
Figure 8.The impact of particle roughness parameter change (σ 2 = 0.15 → 0.5) and cloud top pressure change (200 → 500 hPa).The magnitudes of the differences are comparable while the directional patterns are different.In this plot, the solar zenith angle is 54 • and the viewing zenith angle is 30 • .

Figure 10 .
Figure10.The distribution of inferred EOF 1 scores for synthetic data with and without noise.The distribution for the noise-added synthetic data are symmetric about the EOF 1 score corresponding to the true roughness.The median of EOF 1 score is −0.00336, corresponding to roughness parameter of σ 2 = 0.14.

Figure 11 .
Figure11.Frequency distribution of the χ 2 values (variancenormalized residual square sum).The distribution has a peak at about 12, tapering to nearly zero at approximately 40.This is a reasonable distribution because most pixels contain 15 to 24 observations.

Figure 12 .
Figure 12.The distribution of EOF 1 scores obtained from cold ice clouds over extratropical oceans during September 2005.The median of the EOF 1 score is −0.0293, corresponding to a roughness parameter of 2.82.Consistent with previous studies, roughened particles better simulate the measured polarized reflectivity.

Figure 13 .
Figure 13.Distributions of χ 2 values in the tropics and extratropics.The distribution of the χ 2 value in the tropics (b) implies that the forward model is not correctly simulating the reflectivity in the tropics, while the distribution of the χ 2 value in the extratropics (a) indicates successful inversion.

Figure 14 .
Figure 14.Distributions of EOF 1 and EOF 2 scores with different particle shapes.The observation frequency is shaded with color, and the solid line connects the EOF scores for 10 prescribed roughness values (circles).(a) The result of inference with aggregate of columns, (b) hollow column, (c) solid bullet rosette, and (d) aggregate of plates.

Figure 16 .
Figure 16.Comparison of (a) a typical cloud scene in the extratropics and (b) a cloud scene in the tropics where the χ 2 values are much larger than expected.Green circles are inference locations where the χ 2 value is less than the 95th percentile of the χ 2 distribution, whereas magenta crosses are inference locations where the χ 2 value exceeds the 95th percentile.These figures indicate that the causes of a large χ 2 value may be different in the extratropics and tropics.

FrequencyFigure 17 .
Figure17.Histograms of the coefficient of correlation between EOF 1 and EOF 2 scores.Out of 79 192 total inferred pixels, 49 902 pixels are selected by the condition SD(EOF 1 Score) < 0.02, SD(EOF 2 Score) < 0.02.The results in Figs.12-14 are based on the data within the center six bins in these histograms.

Table 1 .
PARASOL pixel and view selection criteria.