The scientific interest in quantifying the surface energy and mass balance (SEMB) of the ice sheets is obvious when considering the potential impacts on global climate of variations in the brightness of such regions (Hansen and Nazarenko, 2004; Fettweis et al., 2008; van den Broeke et al., 2011; Rae et al., 2012; Van Angelen et al., 2012; Tedesco et al., 2013; Alexander et al., 2014; Colgan et al., 2014). On the other hand, upper estimates of the projected sea level rise are still very uncertain because of poor knowledge of ice sheet dynamics and the SEMB in the polar areas (Van Angelen et al., 2012), which demand a better characterization of their albedo.

As the albedo of snow-covered regions depends on the microphysical and optical properties of the constituent ice crystals (Wiscombe and Warren, 1980; Aoki et al., 2000; Bartelt and Lehning, 2002; Flanner and Zender, 2006; Bougamont et al., 2007; Dang et al., 2016; He et al., 2018), and on the content of light-absorbing impurities (LAIs) (Warren and Wiscombe, 1980; Hansen and Nazarenko, 2004; Dumont et al., 2014), a better knowledge of such properties and their evolution is a high-priority objective for the remote sensing of the cryosphere (Tedesco et al., 2013; Dumont et al., 2014). In developing techniques for the retrieval of these climatologically relevant parameters (Krabill et al., 2004; Comiso, 2006; Zwally and Giovinetto, 2011; Zatko et al., 2013), particular attention has been given to the Greenland ice sheet (GrIS) (Hori et al., 2007; Tedesco et al., 2011, 2013; Box et al., 2012), for which quantitative explanations of the ascertained mass losses are still lacking (Velicogna, 2009). Aerosol (especially black carbon) transport (Dumont et al., 2014; Tedesco et al., 2016) and deposition (Doherty et al., 2014), consolidation of LAIs (including black and organic carbon, and dust) through melt and scavenging, exposure of bare ice, and biogenic contributions cause large variations in the albedo of the GrIS and deeply affect its SEMB, especially in the topographically complex coastal region known as the ablation zone (Alexander et al., 2014; Schmale et al., 2021). A reduction in surface albedo can increase the likelihood of wide swaths of melt, which can fundamentally change the structure of the ice sheets, thus leading to further changes in the albedo and meltwater retention. It is, therefore, essential to understand these factors and how well these changes are captured in the climate models used to simulate and project the GrIS SEMB. Advanced Radiative Transfer (RT) algorithms have been developed to extend the current knowledge to operational products of climatological relevance (Nolin and Dozier, 2000; Stamnes et al., 2007; Jin et al., 2008; Zege et al., 2008; Lyapustin et al., 2009; Painter et al., 2009; Kokhanovsky et al., 2011; Stamnes et al., 2011; Stroeve et al., 2013). For example, the retrieval of grain size in conjunction with standard thermal infrared techniques to infer surface temperature such as split-window algorithms (Key and Haefliger, 1992), gives information on the snow age and metamorphic state and, therefore, an indication of the onset of melting. Many of the retrieval schemes have been applied to the reflectance data from the MODerate resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua satellites, but the information content is often insufficient to unambiguously resolve the vertical distribution of grain size, nonsphericity, and LAI content (Warren, 2013) in the snow medium, critical to accurately estimate the albedo and its recent trends (Xie et al., 2006; Tedesco et al., 2013; Dumont et al., 2014; Dang et al., 2016; Tedesco et al., 2016).

Poorly constrained by the correct optical properties, many RT studies model snow as a collection of spherical grains. It is now established that such an assumption biases the derived albedo low by a few percent (Xie et al., 2006; Tedesco and Kokhanovsky, 2007; Kokhanovsky et al., 2011; Libois et al., 2013; Tedesco et al., 2013; Dumont et al., 2014; Dang et al., 2016; Räisänen et al., 2017). Such biases can be largely amplified by the feedback processes, measurably altering the modeled forcing in polar regions. Furthermore, since the crystal habit determines the asymmetry parameter, the assumption of spherical grains also impacts grain size retrievals, with additional inaccuracies in determining shortwave absorption. In fact, errors in the albedo due to wrong assumptions on grain shape can be offset by adjusting the grain size (Dang et al., 2016; Räisänen et al., 2017).

Decoupling the competing contributions of grain size and LAIs in determining absorption within the snowpack is an unresolved issue of undoubted urgency (Dumont et al., 2014; Tedesco et al., 2016), also affecting the altimetric measurements (Kwok, R., priv. comm.). Differentiating LAIs in the snow from atmospheric aerosols, of paramount climatological importance (Hansen and Nazarenko, 2004; Dumont et al., 2014), has proven to be very difficult (Warren, 2013).

A most striking feature of ice is the eightfold-order-of-magnitude increase in absorption as the observational wavelength shifts from the visible part of the spectrum toward the Near-InfraRed (NIR) (Warren, 2019). For this reason, the total reflectance in the NIR is virtually unaffected by the presence of extraneous LAIs (Wiscombe and Warren, 1980), while in the visible even a minimal amount of LAIs is sufficient to cause a drastic reduction in the albedo. This phenomenon has received much attention since the darkening of the snowpack increases the absorption of radiation in the snow, potentially generating a runaway feedback mechanism (Nolin and Frei, 2003; Hansen and Nazarenko, 2004).

In terms of grain shape, the first natural approach to snow modeling is that of assuming a collection of spheres. Due to the lack of information about crystal shape, some authors (Stamnes et al., 2007) have employed the Henyey–Greenstein phase function (Henyey and Greenstein, 1941), after running Mie computations on spheres of the desired effective radii to determine the asymmetry parameter. The H–G phase function is an ad hoc way to mimic particle roughness by eliminating the resonant structures (rainbow and glory) associated with scattering from perfect spheres, not observed in snow. A growing body of research points out the limitations of using spherical shapes in the modeling of optical properties of snow grains (Xie et al., 2006; Tedesco and Kokhanovsky, 2007; Kokhanovsky et al., 2011; Libois et al., 2013; Tedesco et al., 2013; Dang et al., 2016; Räisänen et al., 2017; He et al., 2018). A database of ice crystal optical properties was therefore constructed at the NASA Goddard Institute for Space Studies (GISS), based on ray tracing and Geometric Optics (GO) (Macke et al., 1996; van Diedenhoven et al., 2012). Although snow and ice crystals can assume complex aggregated shapes, we used the fact that hexagonal columns and plates have been shown to be adequate radiative proxies for more elaborated structures (Fu, 2007; van Diedenhoven et al., 2012, 2014a). The results for hexagonal prisms were tabulated as a function of the ratio of their dimensions (the Aspect Ratio, AR, with AR>1 for columns and AR<1 for plates) and the microscale roughness of the crystal facets (Macke et al., 1996; Yang et al., 2008). Analogous to the Cox and Munk variance used for the water-wave slopes (Cox and Munk, 1956), the latter numerically represents the standard deviation of the distribution of angles used to randomly perturb the orientation of a crystal facet when the beam encounters it in the GO simulations. For a given combination of aspect ratio and roughness, the asymmetry parameter does not depend on the crystal size for the weakly absorbing wavelengths (van Diedenhoven et al., 2014a, 2012) such as 864 nm. This database yielded successful polarimetric retrievals for parameters descriptive of the crystals forming ice clouds (van Diedenhoven et al., 2014b).

The implementation of these aspects into a retrieval scheme addresses the long-standing challenges which limit the accurate determination of the relationship between the snow grain size and shape with the surface temperature, LAI content, location, and season. Retrievals exploiting the GO database of nonspherical ice crystals have been tested on airborne measurements over snow surfaces collected with the NASA GISS Research Scanning Polarimeter (RSP), and confirmed the superior capabilities of polarimetry to retrieve the snow grain size, shape, and microscale roughness (Ottaviani et al., 2015, 2012, hereafter referred to as O15 and O12). At the heart of such retrievals, and also used here, is a forward RT model based on the general doubling-adding formalism described by De Haan et al. (1987), which is capable to simulate the stokes vector of the light field at any desired set of illumination and viewing angles and altitude above the surface. Both the light field and the partial derivatives with respect to the individual parameters can be determined from a single run of the code without the additional computational cost, should the linearization of the model be desired in order to feed an inversion scheme (Rodgers, 2000). The code includes a rigorous treatment of gaseous absorption and multiple scattering processes by molecules, aerosols, and clouds. Bimodal, log-normal aerosol size distributions (Hansen and Travis, 1974) are typically employed for aerosols. The code has been validated against the benchmark results (see Appendix in De Haan et al. (1987)), and is routinely used at GISS to model the RSP measurements over land and ocean for the retrieval of aerosol, cloud, and surface properties (e.g., Waquet et al. (2009); Knobelspiesse et al. (2011, 2012); Ottaviani et al. (2019)).

In a preliminary study, we applied a root-mean-square-error search through the GO database to the surface polarized reflectance signal isolated after an accurate atmospheric correction from the measurements acquired during a survey of alpine snowfields in the Sierra Nevada range (O12). The collection of these measurements, the first of its kind, was personally designed during the Carbonaceous Aerosol and Radiative Effects Study (CARES) mission. The results pointed to snow composed of, or behaving as, a collection of ice crystals with roughness values of ∼0.3 and above, a range of quite extreme aspect ratios, and correspondingly high asymmetry parameters (g ∼ 0.85). Subsequently, a dedicated science flight over a dozen targets in Colorado and California was designed during the POlarimeter Definition EXperiment (PODEX) campaign, which served as a training activity for the inter-agency Studies of Emissions and Atmospheric Composition, Clouds, and Climate Coupling by Regional Surveys (SEAC4RS) mission, with the RSP deployed on the ER-2 aircraft. With 95% of the atmospheric mass residing below the nominal cruising altitude, this stratospheric platform provides a perspective close to that of a satellite. The analysis (O15) replicated the results from the previous limited RSP dataset, confirming the solidity of the method. Building on the improved estimates of the asymmetry parameter, the results were also augmented with grain size retrievals from the RSP shortwave infrared measurements in a MODIS-like fashion (Lyapustin et al., 2009).

The O12 and O15 studies did not include the basic plots showing the sensitivity of the polarization signal to the crystal habits spanned by the GO dataset. For this reason, Figure 1 presents the simulations of the radiances emerging from a 1 m-thick, single-layer snowpack (density = 0.1 g/cm³, optically semi-infinite), with an exponential distribution of hexagonal prism ice grains of different ARs, and with effective radii of 50 μm typical of fresh snow conditions. The roughness parameter was set to 0.3 as retrieved from the real data (O15, O12). The upper and lower rows correspond to the total reflectance and the polarized reflectance, respectively, defined as R _p = Q 2 + U 2 , where Q and U are the second and third components of the Stokes vector. The columns correspond to different wavelengths in the visible and NIR (469, 864 and 1,589 nm), which give similar results (with the obvious exception of the total reflectance at the absorbing wavelength of 1,589 nm) because the scattering from ice crystals belongs to the geometrical optics regime and the medium is optically semi-infinite. These simulations are performed for an RSP-like-instrument scanning along the principal plane of reflection, so as to collect a wide range of scattering angles (top axis).

Since the intention of Figure 1 is to highlight the effects of crystal shape on the surface contribution (i.e., on a signal that has hypothetically been atmospherically corrected), no atmosphere is included. The sensitivity to crystal shape emerging from these simulations is an important result, given that the total reflectance at absorbing wavelengths is used to estimate the grain size (in the top layer of a centimeter or so depending on the specific penetration depth), and demonstrates that the retrievals can be improved by including polarization to achieve a better determination of the grain shape. Albeit low, the magnitude of the polarized reflectance does not cause detection problems given the polarimetric accuracy and the signal-to-noise ratio of heritage (POLDER (Fougnie et al., 2007; van Diedenhoven et al., 2012)) and upcoming polarimeters, such as those slated for launch onboard the NASA Plankton, Aerosol, Ocean, and Ecosystem (PACE) and the European MetOP-SG (Lang et al., 2019) satellite missions. More importantly, the polarized signal carries specific angular signatures. A substantial variation in the behavior occurs for crystal plates across the AR value of 0.5, mirrored by a jump in the asymmetry parameter as large as 0.05. This situation is replicated to a lesser extent by columns as the AR varies across the reciprocal value (i.e., AR = 2).

To simulate the sensitivity of satellite observations to impurities in or above the snow, a crystal habit with AR = 1 (blue curves in Figure 1) was used to model the bottom layer of a model atmosphere. Contamination was added in the form of plausible, increasing concentrations of 1) impurities embedded in the snowpack (0, 0.01, 0.1, and 10 ppmw of soot as in Fig. 7 of Warren and Wiscombe (1980)), with a refractive index n = 1.8–0.6i and density of 2 g/cm³) and 2) tropospheric aerosol similar to the “Boreal Forest” smoke class by Dubovik et al. (2002), with a refractive index n = 1.5–0.01i and the optical depths at 555 nm equal to 0, 0.1, 0.4, 0.8, and 1.6. The impurities in the snowpack are considered as externally mixed to the snow grains, and their optical properties are calculated by Mie calculations internal to the RT code, as done for the aerosols. The results of these simulations are shown in Figures 2, 3, where the columns pertain to the same wavelengths as in Figure 1. The TOA polarized reflectance at 469 nm is dominated by Rayleigh scattering. At all wavelengths, comparing the two figures, it is evident that the model predicts very different TOA reflectances depending on whether the contamination is embedded in the snowpack or suspended above it (tropospheric aerosols), because the overwhelming amount of multiple scattering in the snowpack suppresses any polarization effect due to the small fraction of LAIs. A complete sensitivity study on an extended set of atmospheric aerosols is outside the scope of this article, and does not invalidate the basic finding that the polarization signatures of aerosols overwhelm those of the impurities, which are de facto null.

We also modeled (not shown) the impurities embedded in the snowpack with the optical properties typical of dust (n = 1.5–0.01i, r _eff = 1 μm, v _eff = 0.3, density = 2.6 g/cm³ and concentrations of 1, 10, 100, 1,000 ppmw), obtaining the same insensitivity of the polarized reflectance as obtained for soot. These findings confirm that it is impossible to differentiate between the two kinds of contaminants as long as their refractive indices, although different, are taken as spectrally invariant. To this regard, note that the recent efforts have aimed at estimating the absorption properties of these contaminants further into the UV region where dust absorption increases dramatically (Müller et al., 2008; Moosmüller et al., 2009; Wagner et al., 2012), since some of the current space-borne sensors such as the Thermal And Near-infrared Sensor for carbon Observation (TANSO), the Cloud and Aerosol Imager (CAI) on GOSAT, and the JAXA Second-generation GLobal Imager (SGLI) on GCOM-C1 have moderate (1 km) resolution observations at 380 nm.

The simulations presented in this study confirm the importance of accounting for the grain shape in remote sensing of snow properties, and illustrate the benefits of polarimetry in correctly partitioning LAIs between the snowpack and the atmosphere. For instruments equipped with polarization and multiangular capabilities, they suggest that a comprehensive database of ice crystals’ shapes can be exploited within a simultaneous inversion scheme to largely improve the description of both the surface and atmospheric parameters, and avoid unphysical correlations stemming from a lack of information content. After a tailored atmospheric correction procedure applied to dedicated airborne measurements collected by the RSP sensor and previously published (O15 and O12), we are in the process of constructing a general inversion scheme, and we prospect in the near future to validate the model simulations presented here starting from the augmentation of MODIS-like observations with the data from the POLDER sensor, insofar the only space-borne polarimeter with the sufficient capabilities ever to furnish a reliable data record.

The snow grain shape retrievals hinge on the exploitation of multi-angle, polarimetric measurements at NIR wavelengths. Together with traditional, MODIS-like techniques for the retrieval of grain size and impurity content based on VIS and SWIR channels, the results are therefore of special interest for the development of novel algorithms to be applied to the next generation of polarimeters planned for upcoming space missions, such as PACE and those identified in the Decadal Survey Report, or the European 3MI (Lang et al., 2019) sensor onboard MetOP-SG. The results from polarimeters with limited angular capabilities, such as JAXA’s SGLI sensor, are to be evaluated with rigorous information content studies.