Photopolarimetric retrievals of snow properties

Abstract. Polarimetric observations of snow surfaces, obtained in the 410–2264 nm range with the Research Scanning Polarimeter onboard the NASA ER-2 high-altitude aircraft, are analyzed and presented. These novel measurements are of interest to the remote sensing community because the overwhelming brightness of snow plagues aerosol and cloud retrievals based on airborne and spaceborne total reflection measurements. The spectral signatures of the polarized reflectance of snow are therefore worthwhile investigating in order to provide guidance for the adaptation of algorithms currently employed for the retrieval of aerosol properties over soil and vegetated surfaces. At the same time, the increased information content of polarimetric measurements allows for a meaningful characterization of the snow medium. In our case, the grains are modeled as hexagonal prisms of variable aspect ratios and microscale roughness, yielding retrievals of the grains' scattering asymmetry parameter, shape and size. The results agree with our previous findings based on a more limited data set, with the majority of retrievals leading to moderately rough crystals of extreme aspect ratios, for each scene corresponding to a single value of the asymmetry parameter.


Introduction
The development of accurate techniques for the retrieval of climatologically relevant parameters in snow-covered regions is of obvious importance for the success of climate model predictions (Comiso, 2006;Zwally and Giovinetto, 2011;Zatko et al., 2013). For instance, the retrieval of grain size and surface temperature leads to information that can be used to map the snow melting state, valuable to understand the mass balance of the ice sheets Lyapustin et al., 2009). Passive optical sensors are employed in this regard to evince the optical and microphysical properties of snow grains from measurements of snow reflectance, and advanced algorithms have been applied to data collected by the MODerate-resolution Imaging Spectroradiometer (MODIS) and sensors with similar capabilities Aoki et al., 2007;Jin et al., 2008;Stamnes et al., 2011). However, the limited information content of such measurements translates into uniqueness problems during the retrieval process, with the result that operational snow products are still limited to snow covered area and albedo (with the exception of the recent addition to MODIS of the dust radiative forcing product as described in Painter et al., 2012). Polarimetric remote sensing (Diner et al., 2007;Cairns et al., 2009;Dubovik et al., 2011;Hasekamp et al., 2011;Chowdhary et al., 2012;Ottaviani et al., 2012a;Knobelspiesse et al., 2012) adds to MODIS-type observations because it can provide information on the ice grain shape Ottaviani et al., 2012b) and an accurate characterization of the optical and microphysical properties of aerosols present in the scene, including possible impurities embedded in the snowpack. The benefits of passive polarimetric observations originate from both their extreme sensitivity to particle microphysics and the specifics of surface polarized reflectance. For conventional land types (soil, vegetation), the latter can be effectively parameterized as a scaled Fresnel reflection function with negligible spectral dependence. As a consequence, the separation of the signal into its surface and atmospheric contributions (i.e., the atmospheric correction procedure) is achieved with improved accuracy.
There certainly is debate on the applicability of the radiative transfer theory to such a densely packed medium as snow. However, the formalism is recognized to model the total reflectance with sufficient accuracy for reasons likely associated with the large number of scattering events taking place within the snowpack. Based on this assumption, the method is expected to perform even better for the polarization component because the polarimetric signatures of a medium are known to originate from its top layer (van Diedenhoven et al., 2013): deeper into the medium (i.e., past the first units of optical depth) multiple scattering randomizes polarization.
The NASA GISS Research Scanning Polarimeter (RSP; Cairns et al., 1999), born as an aerosol research instrument, provides a target's total and polarized reflectance with high accuracy in nine bands between 410 and 2264 nm and at 152 viewing angles per scene. Specifically, the dimensionless reflectances are obtained from the measurements of the first three components of the Stokes vector (I , Q and U ) as follows: where the (annual average) extraterrestrial solar irradiance F 0 is corrected for the Sun-Earth distance r 0 and the solar zenith angle θ s . A recent deployment led to the first published study of the polarimetric signatures of snow reflectance from airborne observations (Ottaviani et al., 2012b, hereafter referred to as O12). In that case, the measurements were acquired over a topologically complex terrain in the Yosemite National Park region and processed with an automated, iterative atmospheric correction scheme purposely created to analyze the spectral dependence of the polarized reflectance. A preliminary retrieval of crystal microphysical and optical properties was also tested. The purpose of this new work is to apply the same atmospheric correction procedure and to extend the surface properties retrieval scheme to a more recent data set of improved quality, in terms of the flatness of the surveyed terrain.
In the next section we describe the data collection, and in Sect. 3 the processing scheme. In Sect. 4 we discuss the results obtained from the analysis. The conclusions are found in Sect. 5.

Data collection
The data analyzed here were collected onboard the highaltitude ER-2 aircraft during the POlarimeter Definition EXperiment (PODEX) mission, a NASA effort to ultimately test and compare the performance of different polarimeters, as part of the Decadal Survey activities and in preparation of the Studies of Emissions and Atmospheric Composition, Clouds and Climate Coupling by Regional Surveys (SEAC 4 RS) campaign. Compared to more traditional aircraft such as the King Air B200 used in O12, the ER-2 offers the advantages of a perspective more similar to that of a satellite. For example, interference caused by cirrus clouds is readily assessed using the RSP 1880 nm channel, which is usually very dark because of the strong absorption by tropospheric water vapor that masks the surface, and is therefore very sensitive to high-altitude clouds. There is in fact interest in exploiting this wavelength as a cirrus screener in a broader context, since cold places tend to have dry atmospheres and the reflectance of snow at 1880 nm is much less than at, say, 1380 nm (the effective strength of absorption by ice at 1880 nm is twice that exhibited at 1380 nm). Also, RSP operations on the ER-2 benefit from the stable attitude intrinsic to stratospheric flights, which enhances data quality since the RSP scans along the flight track (approximately 50 • on either side of nadir) to assemble the multiangle view of a target's reflectance from subsequent scans.
About a dozen snow-covered targets were selected for the 6.8 h flight over Colorado and the Sierra Nevada range in California on 15 March 2012. Well-formed snowpacks with little to no interference from vegetation and flatness of terrain were the main priorities. Three scenes were chosen for this analysis. The Grand Mesa plateau (39 • 08 N, 108 • 12 W, ∼ 3000 m a.s.l.) in western Colorado exhibits remarkable flatness and was covered by a consistent snowpack with a depth of greater than 1 m as reported by nearly simultaneous in situ surveys (Mesa Lakes snow telemetry (SNO-TEL) site, managed by the Natural Resources and Conservation Service of the US Department of Agriculture; M. Skiles, personal communication, 2013). Lake Granby (40 • 10 N, 105 • 53 W, ∼ 2500 m a.s.l.), the third-largest body of water in Colorado and second-largest reservoir, was frozen at the time of the overpass and covered with a thin layer of aged, likely metamorphosed snow (the SNOTEL site at Stillwater Creek recorded a snow depth of ∼ 0.5 m). The scene denoted by Derby Peak (39 • 59 N, 107 • 10 W, ∼ 3500 m a.s.l.) is another high-altitude, rather flat region (in fact belonging to the Flat Tops Wilderness Area). The closest SNOTEL site, at Trapper Lake, 2950 m a.s.l., reported a snow depth of 0.7 m.

Methodology
The method employed for the atmospheric correction procedure is explained in detail in O12, to which the reader can refer if interested in the rigorous mathematical formalism. To provide a brief summary, the signal measured by the instrument is a nonlinear function of the surface reflectance due to the interaction of the radiation reflected by the surface with the atmosphere above. For this reason, a linear regression cannot be applied to infer the model parameters. Nevertheless, a cost function with the departure of the model from the measurements can still be computed and minimized, provided the availability of the derivatives (Jacobians) of the radiative transfer model with respect to the descriptive parameters. The Jacobians, generated without extra-computational cost from the linearization of the doubling-adding radiative transfer code in use at GISS (De Haan et al., 1987), are used in an iterative scheme based on an inversion and optimization procedure of the Gaussian-Newton kind (Rodgers, 2000) to search for a minimum of the cost function. The result is the decoupling of the total signal measured at instrument altitude into its surface and atmospheric contributions, i.e., an automated atmospheric correction. Within this scheme, the surface total reflectance is expressed as the weighted sum of three kernels as in the RossThick-LiSparse form (Wanner et al., 1995;Lucht et al., 2000): where is the scattering angle. The three terms on the righthand side describe different scattering mechanisms: a Lambertian component for isotropic reflection, a volumetric term K vol for scattering from small structures and a geometric term K geo for a macroscopic shadowing from large objects. Although this model was created for reflection from surfaces other than snow, its mathematical form adapts without problems to snowpacks, with the isotropic kernel systematically assuming the largest weight because of the dominance of multiple scattering in the reflection of light from snow. As far as the polarized reflectance is concerned, we instead employ a scaled Fresnel kernel elsewhere used for vegetated surfaces (Nadal and Bréon, 1999;Waquet et al., 2009): with the fractional polarized reflectance F p (θ i ) given by (Born and Wolf, 1999): where in our case n 1 and n 2 are the indices of refraction of pure air and ice, respectively, and the angle of refraction θ t is connected to the angle of incidence θ i by Snell's law. For the retrieval of crystal habit, the polarized reflectance at 864 nm is used because of the very weak ice absorption at this wavelength and the relatively weak Rayleigh scattering contribution. Ice crystals are found to vary within a virtually countless set of shapes, but, limited by the lack of constraints on the correct optical properties, very few studies have departed from the common practice of assuming spherical shapes for the snow grains, an approximation which can negatively impact the quality of the albedo products (Tedesco and Kokhanovsky, 2007;Aoki et al., 2000). Among the available optical models for non-spherical grains are for example those of Baum et al. (2011Baum et al. ( , 2014 and Yang et al. (2013). Incidentally, Räisänen et al. (2015) recently attempted to fit measured phase functions of blowing snow by using phase functions of several complex habits provided by Yang et al. (2013) in addition to phase functions of other habits. However, the selection of habits in their study remained limited, and to fit the data a rather arbitrary mixture of droxtals, aggregates of plates and Koch fractals was needed, leading them to conclude that this mixture "most probably does not represent properly the actual distribution of snow grain shapes in blowing snow (or snow on ground)". In our case, the description of the snowpack is refined by treating it as a compact, optically semi-infinite medium composed of ice crystals modeled as hexagonal prisms with variable aspect ratios and microscale facet roughness.
It should be noted that for the specific task of remote sensing of aerosols over snow, a model capable of accurately reproducing the contribution of the surface to the total signal should be considered successful regardless of its capability to mimic the true shape of the crystals. Our approach is based on the recognition that hexagonal plates and columns can be effectively used as radiative proxies for more complex crystal habits (Fu, 2007;McFarquhar, 2007, 2009;Baran, 2009;Ottaviani et al., 2012a). As demonstrated by van Diedenhoven et al. (2012), matching measured polarized reflectances with a look-up table of values simulated by adjusting the aspect ratios and roughness parameters of simple hexagonal particles yields an estimate of the asymmetry parameter of more complex ice crystals. Our database, initially created to be of use for the retrieval of ice cloud properties, contains 765 different combinations of aspect ratios and roughness parameters and was computationally assembled by running Monte Carlo simulations based on geometric optics (Macke et al., 1996), the performance of which was evaluated by . The retrieved aspect ratios are interpreted as the mean of the aspect ratios of the components of the complex crystals.
The microscale surface roughness is statistically accounted for by tilting, for each interaction with a ray, the normal to the crystal surface by an angle randomly selected in the δ · [0 • , 90 • ] interval, where δ is referred to as the roughness parameter. In our case, δ is limited to 0.7 (i.e., a maximum tilt angle for the facet of 63 • ) because for higher values the probability of unphysical scattering events strongly increases, resulting in progressively larger loss of accuracy of the GO calculations of a given number of rays. Neshyba et al. (2013) have shown that various other definitions and parameterizations of crystal surface roughness lead to similar results. Yang et al. (2008) found that such approaches are efficient yet accurate treatments of microscale roughness. Furthermore, the effects of microscale surface roughness and macroscale crystal distortion were shown to be largely equivalent by Liu et al. (2014). Atmospherically corrected polarized (second row) and total (third row) reflectance for three snowfields (columns) overflown in Colorado, USA. The first row shows real-time imagery from the high-resolution camera onboard the ER-2, captioned by the scene altitude above sea level together with the solar zenith angle and the Sun-sensor relative azimuth at the time of observation. The last row pertains to the search of the optimal ice grain model: the contour plots map the RMSE of the fit to the polarized reflectance at 864 nm when simulating the snowpack with each crystal habit (i.e., each combination of aspect ratio and roughness parameter) in the database. The red circles locate the values of the RMSE falling under a predetermined threshold and are taken to represent the retrieved (optimal) habits.
The retrieval of grain size is instead based on the total reflectance at selected shortwave infrared (SWIR) bands where ice absorption is high. The total reflectance of an optically semi-infinite snow layer at 2264 and 1594 nm is effectively determined by the asymmetry parameter and the single scattering albedo at the respective wavelength, which in turn is mainly determined by the effective diameter of the ice crystals. The grain diameter is not unequivocally defined in the literature, which is understandable since the definition of snow grain itself is often ambiguous. There is a consensus for expressing grain effective diameter as proportional to the ratio of crystal volume to projected area, but the value of the multiplier can vary according to the specific choice of the author Kokhanovsky and Zege, 2004;Aoki et al., 2000). In line with several recent publications (see, for example, Zege et al., 2011;Jin et al., 2008), we define it as 3/2 times the ratio of the crystal volume over its surface area, as outlined in Appendix A, which includes the details of the relationship to geometric size. The effective diameter of snow for a given scene is estimated by matching the measured reflectance with an ice crystal model whose asymmetry parameter is consistent with that retrieved for that scene using the method described above. To this end, various modified gamma size distributions were applied to yield varying effective diameters, and the radiative transfer code was used to simulate the corresponding 1594 and 2264 nm reflectances . For conditions equivalent to those found in snowpacks, Bi et al. (2014) showed that the errors in retrieved effective radius attributable to the use of the conventional geometric optics approach are below 5 %. Another important aspect to consider is that the 1594 and 2264 nm channels experience different penetration depths, with the former weighted more towards the top layer, while the latter probes deeper into the medium; as a consequence, the retrieved sizes contain information on the vertical structure of the snowpack Li et al., 2001;Warren, 1982).
As far as the atmospheric component is concerned, the high elevation of these alpine regions allowed us to assume negligible contributions from aerosols, so that in the radiative transfer code a standard Rayleigh atmosphere was used with a surface pressure extrapolated at the indicated altitude. Cirrus-free conditions were guaranteed from the screening procedure described earlier. The SWIR channels were corrected for the minimal interference caused by absorption from a standard amount of trace gases.

Results and discussion
The core of our findings is condensed in Fig. 1. For each of the three scenes (different columns in the figure) the analysis was performed over a single RSP pixel, corresponding to a spatial resolution of about 225 m. Indeed, averaging over a few adjacent aggregated scans did not show appreciable The "butterfly" pattern of the contour plot for the asymmetry parameter g as a function of aspect ratio and roughness parameter resembles that obtained for the RMSE plots in Fig. 1. The red circles in Fig. 1 (here in green for Grand Mesa, magenta for Lake Granby and blue for Derby Peak) that identify the retrieved optimal crystal habits are found for each scene clustered around a single value of g.
scan-to-scan variability, revealing more or less uniform snow conditions on a scale larger than the instrument instantaneous field of view. Real-time imagery from the high-resolution camera onboard the ER-2, overlaid to Google Earth, is included for context. The second and third rows of panels show the results of the iterative correction scheme applied to the relevant RSP channels. With the exception of O12, the visible and SWIR behavior of the polarized reflectance of snow has not been previously published. The total reflectance exhibits the familiar high spectral albedo in the visible, albeit with some expected scene-toscene variability. Lake Granby in particular shows lower total reflectance than the other two scenes; in situ simultaneous surveys from the personnel at the US Forest Service in Arapaho National Forest reported a snow depth on the order of 5-8 cm. These values are close to the limit of thickness required to visually mask the underlying surface, leaving the possibility that the lake ice underneath the snow cover is darkening the reflectance. Nevertheless, a layer of a few millimeters is sufficient to fulfill the semi-infinite approximation at the longer wavelengths used for the retrieval of grain size: the RSP SWIR channels at 1594 and 2264 nm are very dark because of the strong absorption of ice at these wavelengths.
The polarized reflectance as a function of the scattering angle is very similar among the three scenes. Note that in the case of Derby Peak the flight track was oriented at about 45 • with respect to the principal plane (i.e., the plane defined by the solar azimuth and the normal to the surface), leading to www.the-cryosphere.net/9/1933/2015/ The Cryosphere, 9, 1933-1942, 2015 a smaller range of scattering angles collected by the RSP. It has already been observed in O12 as the land model used for the polarized reflectance converges to zero toward backscatter (scattering angle > 160 • ), and in this region agreement with the data is not to be expected. It is useful to quantify the spectral spread of the channels by examining the differences from the reference 2264 nm channel. The signal in this band is in fact a proxy for the surface contribution, since the Rayleigh (and fine-mode aerosol) contribution becomes more and more negligible as the observation windows shifts to the SWIR. If the spectral differences are small, that same signal can then be subtracted from the shorter RSP wavelengths when aerosol retrievals are attempted. The residuals with respect to the 2264 nm polarized reflectance are found to be within 0.005, which as in O12 is larger than conventional land surfaces but still very small. The polarized reflectance is small too in absolute value. Shifting the attention to the search of a snow model, the plots at the bottom contour the root mean square error (RMSE) of the data fit to each combination of aspect ratio and roughness parameter in the database. Multiple combinations led to a satisfactory fit to each of the three polarized reflectance observations, so we considered as optimal all habits with an RMSE falling below a small threshold value (4.5 × 10 −7 ). These minima, marked with red circles, therefore correspond to those crystal habits that when used to model the snowpack at the bottom of the simulated atmosphere yield the best fits to the data. The analysis shows that the polarized reflectances are consistent with particles that have components that are like thin plates or like thin columns, a fact already observed in O12. Interestingly, the butterfly-like topology of the RMSE was recognized to closely resemble that obtained for the asymmetry parameter g plotted versus the same coordinates (Fig. 2). In each scene, the optimal crystal habits are associated with very similar values of the asymmetry parameter (0.84 for Lake Granby, 0.876 for Derby Peak and 0.90 for Grand Mesa). This result is understood by considering that the asymmetry parameter is the main descriptor for the scattering properties. As demonstrated by  using a large variety of complex ice crystal shapes, all the combinations of aspect ratio and roughness parameter fitting a polarized reflectance measurement are ultimately characterized by similar values of g. Consequently, they will also yield similar particle size retrievals. Figure 3 shows the phase function and the degree of linear polarization, obtained for each scene by averaging the phase matrices associated with the retrieved minima. The lack of halo peaks in the 20-40 • range, typical of the phase functions of pristine crystals, and the general smooth behavior (except from the strong forward peak) are consistent with the moderate-to-high roughness parameters determined in the analysis. All retrieved parameters are reported for comparison in Table 1, together with the results of the size retrievals. Using the simple approach outlined in the previous section, we obtained effective diameters equal to 152, 182 and 144 µm for the Grand Mesa, Lake Granby and Derby Peak scene, respectively, when using the RSP band at 2264 nm. Using instead the 1594 nm channel, diameters typically 30 µm smaller are obtained for each scene, which can be explained by the different penetration depths of the two channels, indicating crystal size increasing with depth as expected Li et al., 2001;Aoki et al., 2000;Warren, 1982). These values, at the lower end of common retrieval ranges, are normally associated with fresh snow conditions as found in studies based on other remote sensors in alpine regions (Negi and Kokhanovsky, 2011;Painter et al., 2009;Dozier and Painter, 2004;Nolin and Dozier, 2000) and polar plateaus (Lyapustin et al., 2009;Hori et al., 2007). However, a direct comparison is challenged by differences in the definition of grain size, assumptions on grain sphericity and different penetration depths achieved by the use of different instrument channels. The larger grain sizes obtained for the Lake Granby scene were attributed to the shallow aged snow over the lake, since metamorphosed snow is recognized to be composed of larger grains (Aoki et al., 2000;Warren, 1982).

Conclusions
We have presented the extension of a novel analysis described by Ottaviani et al. (2012b) to a data set acquired with the high-accuracy RSP airborne instrument overflying highaltitude snowfields. The data were processed within an auto- Table 1. List of ice crystals' parameters retrieved for the three analyzed scenes. The value in parentheses indicates the RSP band used in the retrieval. The effective diameter can be obtained using either the 2264 or the 1594 nm channel (see text), together with a crystal habit characterized by the asymmetry parameter retrieved for the corresponding scene.

Retrieved parameter
Grand Mesa Lake Granby Derby Peak mated retrieval scheme capable of isolating the surface contribution to the total signal. The spectral dependence of the retrieved polarized surface reflectance is slightly larger than for soil or vegetated surfaces, but nonetheless small relative to typical aerosol contributions (cf. Fig. 6 in O12). This fact has important implications for the construction of advanced algorithms which aim at the retrieval of the microphysical properties of aerosol layers located over ice and snowfields. Furthermore, we applied a retrieval procedure which is demonstrated to yield snow microphysical parameters of primary climatological interest. We retrieve snowpacks behav-ing as a collection of grains of extreme geometries (thin plates and/or long columns) and moderate to high microscale roughness. These results reinforce those of our previous study based on a limited data set acquired over very rugged terrain. For each scene, the best fits are obtained for a set of crystal habits characterized by very similar values of the asymmetry parameter and extreme aspect ratios. Size retrievals were also tested, leading to grain effective diameters in the 140-180 µm range at the top of the snowpack.