Glacier Ice Surface Properties in South‐West Greenland Ice Sheet: First Estimates From PRISMA Imaging Spectroscopy Data

Snow and ice melt processes on the Greenland Ice Sheet are a key in Earth's energy balance and are acutely sensitive to climate change. Melting dynamics are directly related to a decrease in surface albedo, amongst others caused by the accumulation of light‐absorbing particles (LAPs). Featuring unique spectral patterns, these accumulations can be mapped and quantified by imaging spectroscopy. We present first results for the retrieval of glacier ice properties from the spaceborne PRISMA imaging spectrometer by applying a recently developed simultaneous inversion of atmospheric and surface state using optimal estimation. The image analyzed in this study was acquired over the South‐West margin of the Greenland Ice Sheet in late August 2020. The area is characterized by patterns of both clean and dark ice associated with a high amount of LAPs deposited on the surface. We present retrieval maps and uncertainties for grain size, liquid water, and algae concentration, as well as estimated reflectance spectra for different surface properties. We then show the feasibility of using imaging spectroscopy to interpret multiband sensor data to achieve high accuracy, frequently repeated observations of changing snow and ice conditions. For example, the impurity index calculated from multiband Sentinel‐3 Ocean and Land Colour Instrument measurements is dependent on dust particles, but we show that algae concentration alone can be predicted from this data with less than 20% uncertainty. Our study evidences that present and upcoming orbital imaging spectroscopy missions such as PRISMA, Environmental Mapping and Analysis Program, Copernicus Hyperspectral Imaging Mission, and the Surface Biology and Geology designated observable, can significantly support research of melting ice sheets.

The most common variable of the cryosphere being monitored from space is the effective snow grain radius. Itexpresses the path lengths of photons that pass through the ice and is also sometimes referred as specific surfacearea. The reflectance in wavelengths of moderate ice absorption depends on these path lengths, so that snowwith larger grains gets darker (Dozier et al., 1981;Warren, 1982). Likewise, the spatial distribution and amount of LAP accumulation can be detected from space. In particular, depositions of algae in snow and glacier ice can be monitored by relying on chlorophyll and carotenoids absorption characteristics (Painter et al., 2001). Algal accumulation can be quantified as concentration in units of cells ml2017; Painter et al., 2001). Finally, the effective grain radius is also an indicator for surface wetness since the crystal size increases due to clustering processes in liquid water enriched snow and ice (Dozier et al., 2009). Alternatively, liquid water content can be expressed as spherical fraction of the snow and ice grains. However, this approach requires a separation of the liquid water and ice absorption lines and can therefore only be pursued by using imaging spectroscopy measurements (Green et al., 2002).
Optical remote sensing of snow and ice surface properties from space was among the earliest geophysical retrieval methods based on satellite missions and is a valuable tool to obtain amount and spatial distribution of different parameters on a global scale with a high temporal resolution (Rango & Itten, 1976). The potential of the near-infrared (NIR) wavelengths to estimate snow grain size was already demonstrated in the early 80's based on measurements from the NOAA Advanced Very High Resolution Radiometer (Dozier et al., 1981). Prominent subsequent missions used to retrieve snow grain size include the Moderate Resolution Imaging Spectroradiometer (MODIS) (Carlsen et al., 2017;Zege et al., 2008Zege et al., , 2011, and the Sentinel-3 Ocean and Land Colour Instrument (S3 OLCI) (Kokhanovsky et al., 2019). The detection of biological LAP on snow and ice surfaces has also been studied in detail and a couple of investigations focused on mapping glacier algal blooms and determining their effects on ice melt on the Greenland Ice Sheet (Cook et al., 2020;Gray et al., 2020;Stibal et al., 2017;Takeuchi et al., 2006;Wang et al., 2018Wang et al., , 2020. These studies applied retrieval algorithms to data from the Satellite Probatoire d' Observation de la Terre, MODIS, S3 OLCI, the Medium Resolution Imaging Spectrometer (MERIS), or Sentinel-2.
In contrast to most of the existing optical satellite missions, imaging spectroscopy can be used to accurately map and quantify snow and ice surface properties using physically based retrievals by modeling characteristic atmospheric and surface absorption features (Painter et al., 2013). So far, this technique has been almost entirely based on airborne spectrometers though, and in particular, on measurements from NASA's Airborne Visible Infrared Imaging Spectrometer (AVIRIS). Approaches to estimate snow grain size from AVIRIS data have been introduced by Nolin and Dozier (1993), and further developed by Nolin and Dozier (2000) and Painter et al. (2013). It has also been demonstrated that concentration of snow algal blooms can be quantified using AVIRIS acquisitions (Painter et al., 2001). The same instrument was used to quantify liquid water (Green et al., 2006). Recently, Bohn et al. (2021) demonstrated a promising potential of spaceborne imaging spectroscopy missions to simultaneously detect and quantify snow and ice grain size, liquid water, and glacier algal accumulation on the Greenland Ice Sheet based on simulated data and AVIRIS measurements. In this context, a new generation of orbital imaging spectroscopy missions is expected to provide much wider coverage on a more regular basis with high resolution of only 30 m. The German Aerospace Center's (DLR) Earth Sensing Imaging Spectrometer (Mueller et al., 2016) and the Italian Hyperspectral Precursor of the Application Mission (PRISMA) (Cogliati et al., 2021) already are 3 of 21 in operation since 2018 and 2019, respectively. Forthcoming missions include NASA's Earth Surface Mineral Dust Source Investigation (EMIT) , the German Environmental Mapping and Analysis Program (EnMAP) (Guanter et al., 2015), the Copernicus Hyperspectral Imaging Mission (CHIME) led by ESA (Rast et al., 2019), and NASA's Surface Biology and Geology (SBG) designated observable (National Academies of Sciences, Engineering, and Medicine, 2018).
In this work, we estimate snow and ice surface properties from existing spaceborne imaging spectroscopy data. We apply a recently developed simultaneous Bayesian inversion of atmospheric and surface state using optimal estimation (OE). The algorithm was introduced by Bohn et al. (2021) and is an extended version of the concept presented in Thompson et al. (2018). It incorporates prior knowledge, measurement noise as well as model uncertainties. We use a data set from the PRISMA instrument in order to map and quantify ice grain size, surface liquid water, and algae mass mixing ratio. The image was acquired over the South-West margin of the Greenland Ice Sheet in late August 2020 capturing the "dark zone" or "k-transect," which is characterized by patterns of clean snow and dark ice featuring high concentration of deposited LAPs (Wientjes et al., 2011). We present retrieval maps and associated posterior uncertainties, as well as estimated reflectance spectra for different surface conditions. We also analyze the optical properties of melt ponds or supraglacial lakes, which are numerous in the selected PRISMA acquisition. In addition to presenting the new spectroscopic retrievals, we finally show how these measurements can be used in concert with multiband data in a comprehensive cryosphere observation system. We demonstrate that simple local regression models applied to multispectral S3 OLCI data can achieve a high degree of alignment with retrieval maps from imaging spectroscopy measurements.

Spectroscopic Snow and Ice Property Retrievals
The algorithm our study is based relies on statistical relationships between surface reflectance spectra and snow and ice properties to estimate the most probable solution state given a particular reflectance. It is based on the principles of OE described by Rodgers (2000) and uses a comprehensive library of reflectance spectra and associated snow and ice surface parameters as a representation of prior knowledge. Bohn et al. (2021) named this approach a "lazy Gaussian" or "lazy prior-driven" inversion since the forward model is a function of the atmospheric state and the surface reflectance, but not of the additional surface parameters. These extra parameters are estimated entirely based on the prior mean and covariance with the surface reflectance. They comprise grain radius, liquid water path length as well as mass mixing ratios of various LAPs. The statistical correlations between reflectance and surface properties are derived from runs of the snow and ice radiative transfer model BioSNICAR-GO.
BioSNICAR-GO simulates surface spectral albedo by combining a bio-optical model with the two-stream multilayer SNow, ICe, and Aerosol Radiation model SNICAR (Cook et al., 2020;Flanner et al., 2007). It facilitates the modeling of ice grains and LAP either as collections of spheres based on Lorenz-Mie theory (Grenfell & Warren, 1999) or as arbitrarily large hexagonal plates and columns using a geometric optics (GO) parameterization from van Diedenhoven et al. (2014). To enable the estimation of surface liquid water, we exploit the coupling of BioSNICAR-GO with a two-layer coated sphere reflectance model for wet snow Green et al., 2002). The model assumes an increased grain radius attributed to a particular liquid water fraction, and is based on a slight shift between the imaginary parts of the spectral refractive index of liquid water and ice (Dozier & Painter, 2004). The coupling with the coated sphere model allows to add liquid water fractions to the Lorenz-Mie simulations, which are then part of the additional surface parameters in the state vector. The derivation of snow and ice surface properties is therefore given implicitly by using BioSNICAR-GO to construct the spectral library of prior surface reflectance and associated snow and ice properties. In case of modeling liquid water in snow, it may probably be more accurate to consider ice-water clusters in the simulation of scattering effects than relying on the coated sphere model from Green et al. (2002). For instance, a weighted mixture might be a better solution compared to a distinct separation of the liquid and solid water phases. However, the setup of the OE-based inversion allows to replace the radiative transfer model used for building the prior spectral library with any other model that is developed for simulating wet snow.
This section presents a brief discussion of the difference in modeling of snow and ice grains, followed by an overview about the forward model and OE in general. We adhere to standard conventions and denote matrices with uppercase boldface letters, and vectors as well as vector-valued functions with a lowercase boldface notation. For in-depth details, the reader is referred to Rodgers (2000), Thompson et al. (2018), and Bohn et al. (2021).

Snow Versus Glacier Ice
Most of the scientific literature on the retrieval of snow and ice surface parameters is focused on snow grain size (see e.g., Kokhanovsky et al., 2019;Nolin & Dozier, 1993, 2000Painter et al., 2013). However, the optical properties of ice crystals are very different compared to snow, which is a mixture of air and ice (Warren, 2019). There is an inner complexity in estimating ice grain dimensions since the transition from snow to glacier ice is a continuum. On ice sheets, snow is compressed by its own weight and with increasing density, present air forms enclosed bubbles. The higher the density and the pressure, the smaller the bubbles get until they finally dissolve to spare pure ice (Warren, 2019).
The most common method to model the shape of snow grains is to assume non-spherical snow particles being arranged as a collection of spheres and to obtain their optical properties from Lorenz-Mie theory (Grenfell & Warren, 1999). This approach is justified by expecting the snow grain radius being much larger than the incident radiation wavelengths. However, this method features clear limitations when applied to surfaces of bare ice since the grains typically appear to be arbitrarily shaped as plates and columns with irregular dimensions (Kokhanovsky & Zege, 2004). To capture this in the modeling, Aoki et al. (2007) proposed to consider length, width, and thickness of the ice crystals instead of the collected-spheres approach. These parameters are likewise the basis of the GO calculations introduced by Kokhanovsky and Zege (2004). Khuller et al. (2021) recently included a specular reflectance term in modeling ice and much earlier, approaches to model lake ice have been presented (Mullen & Warren, 1988). However, the use of modified snow models might also not be appropriate to model glacier ice since they generally consider a transparent medium populated with absorbing and scattering particles. In a glacier, the medium itself is partly absorbing, consisting of both absorbing and transparent scatterers in terms of LAP and air bubbles, respectively. As a consequence, modeling of LAP concentration in glacier ice may be tainted with higher uncertainties than the available models acknowledge.
In this study, we run the "lazy Gaussian" inversion based on both the collected-spheres and the GO method representing the prior distributions. Although the simulated spectra for glacier ice surfaces display the more appropriate prior mean and covariance for our case study, we also applied the Lorenz-Mie based snow spectral library to our PRISMA data set to enable a comparison with the grain radius maps derived from S3 OLCI data. Furthermore, this demonstrates the resulting differences both in spatial distribution and value range of the estimated grain sizes, and therefore, gives an impression of the applicability of the different approaches to model snow and ice grain shape. Figure 1d shows representative surface reflectance spectra of clean snow and dark ice, respectively, with highlighted characteristic absorption features. Abundance of carotenoids and chlorophyll indicates presence of biological impurities on the surface, whereas ice and liquid water absorption bands are used for retrieving grain size as well as liquid water content. The spectra highlight the differences in reflectivity of snow and ice surfaces and thus, confirm the importance of choosing an appropriate prior knowledge for the inversion.

Forward Model
We denote the forward model as a vector-valued function f of the state vector = [ 1, … , n] T yielding the meas- with ϵ representing a random error vector, which in our case includes measurement noise, prior uncertainties in x, and errors due to unknown forward model parameters. Following Thompson et al. (2018), x contains columnar water vapor in g cm −2 and dimensionless aerosol optical thickness (AOT) at 550 nm being an atmospheric part ATM = [ ] T , and the reflectance of each instrument channel as a surface part x SURF . Here, the snow and ice properties are added leading to the extended version SURF = [ 1 , … , m , SURF 1 , … , SURF n ] T . Thompson et al. (2018) use the hemispherical-directional reflectance factor (HDRF) as a representation of the surface reflectance. In contrast, our implementation of the "lazy Gaussian" method optimizes the hemispherically integrated spectral albedo. This approach is limited by the used two-stream snow and ice radiative transfer model BioSNICAR-GO. However, although the HDRF is the more appropriate quantity when modeling measurements of imaging spectrometers (Schaepman-Strub et al., 2006), the use of spectral albedo for applications to the flat parts of the Greenland Ice Sheet can be pursued .
In specific form, f models the wavelength-dependent top-of-atmosphere (TOA) radiance L TOA using a simplified solution of the radiative transfer equation (Chandrasekhar, 1960): where L 0 is the atmospheric path radiance; E dir and E dif are the direct and diffuse solar irradiance arriving at the surface; μ sun is the cosine of the solar zenith angle; T ↑ is the total upward atmospheric transmittance; S is the spherical albedo of the atmosphere; and ρ s is the surface spectral albedo. For simplicity, we assume an infinite, horizontal, and isotropic Lambertian surface as well as clear sky and a plane-parallel atmosphere. Furthermore, surface roughness at subpixel scale tends to make snow reflectance more Lambertian, so that these assumptions ensure validity of using spectral albedo in place of HDRF . The atmospheric flux parameters L 0 , E dir , E dif , T ↑ , and S are functions of x ATM , surface elevation as well as solar and observation geometry. They are derived from radiative transfer simulations using the MODTRAN code (Berk et al., 1989). The prior covariance matrix of x ATM is assumed to be diagonal and unconstrained.
While the first part of the surface state vector, , is expressed by ρ s in Equation 2, the remaining parameters of x SURF , , are not an input to the forward model. They are optimized entirely based on their prior mean and covariance, which are obtained from the prior surface statistics. These statistics are characterized by a multivariate Gaussian distribution of surface reflectance for each instrument channel and the additional surface parameters with a non-diagonal covariance matrix due to expected correlations across (c) the normalized difference snow index calculated from the difference between the visible green and shortwave infrared TOA reflectance; and (d) exemplary surface reflectance spectra estimated from PRISMA TOA radiance data for clean snow and dark ice, respectively. Center wavelengths of characteristic absorption features of carotenoids (Car), chlorophyll (Chl), liquid water (Liq), and ice (Ice) are highlighted with dashed lines. channels. For example, the OE-based inversion assigns a snow surface reflectance to an input TOA radiance spectrum as prior knowledge obtained from the pre-simulated spectral library. This snow reflectance is accompanied by a specific grain size, liquid water content, and algae concentration, which are then optimized along with the reflectance values for each instrument channel by exploiting the prior covariance matrix.

Optimal Estimation
OE acts on two main assumptions: measurement and state vectors as well as the associated errors follow a Gaussian distribution, and the forward model is locally linear. Then, f can be inverted by minimizing the following cost function, which is the negative logarithm of the posterior probability density function: Here, x a is the prior state vector; S a is the prior covariance matrix; and S ϵ is the measurement covariance matrix. The first term of the right-hand side evaluates the difference between prior and solution state by taking into account S a . The second term penalizes the departure of the modeled TOA radiance from the measurement, weighted by S ϵ , which captures both instrument noise, expressed by the noise-equivalent change in radiance, and uncertainties due to unknown forward model parameters. We assume no correlation between the measurement noise of different instrument channels as well as between the unknown parameters, so that S ϵ is diagonal.
x and x a are not weighted equally for all wavelengths and we use loosely constrained priors in the deep shortwave infrared (SWIR) water vapor features around 1,400 and 1,900 nm. This puts less weight on these wavelengths in the optimization. However, this is completely controlled by the diagonal entries of the covariance matrices. The iteration then searches for the solution state ̂ that leads to a local minimum of Equation 3, being the state with the highest probability given the measurement and the prior state. In this work, we find ̂ using a Gauss-Newton iteration scheme that typically converges in less than 30 iterations. In fact, the optimization may approach a local minimum that in some cases is not equal to the global one. Investigations of correct convergence are underrepresented in existing literature so far (Cressie, 2018), but we assume our approach to be stable as the initial guess tends to be within ±10% of the full probabilistic solution Thompson et al., 2018).
Besides the converged solution state, the OE retrieval scheme reports the posterior predictive uncertainty for each ̂ :̂= where K is the Jacobian of the forward model with respect to ̂ . To facilitate an interpretation of the posterior uncertainties, ̂ can be normalized leading to an error correlation matrix (Govaerts et al., 2010).

Sentinel-3 OLCI Snow Property Retrievals
Measurements from S3 OLCI can be used to derive several snow properties including spectral and broadband albedo, snow specific surface area, snow extent, and snow grain size (Kokhanovsky et al., 2019). Additionally, several multiband indices have been developed for identifying impurities on snow and ice surfaces from instruments such as MERIS, MODIS, or S3 OLCI, including different chlorophyll indices and the impurity index Wang et al., 2018Wang et al., , 2020. In this section, we briefly introduce the S3 OLCI grain size retrieval algorithm as well as the impurity index, as results from both are used for comparison with retrieval maps from PRISMA data.
The snow grain radius is estimated from S3 OLCI data using the following relation (Kokhanovsky et al., 2019): where l is the effective ice absorption length, and A is derived from a scaling constant depending both on snow type and grain shape. Kokhanovsky et al. (2019) suggest A = 0.06 based on findings from various studies, which analyze the scaling constant (see Di Mauro et al., 2015;Kokhanovsky, 2006;Libois et al., 2014). The absorption 10.1029/2021JG006718 7 of 21 length l is derived from the exponential approximation of the asymptotic radiative transfer theory and is calculated by: where R 1 and R 2 are the OLCI TOA reflectance at 865 and 1,020 nm, and α 2 is the ice absorption coefficient at 1,020 nm. f is an angular function that depends on the theoretical reflectance of a non-absorbing snow layer R 0 and on the escape function u(μ): Here, u(μ) is approximated by a linear function applied both to the solar zenith angle μ 0 and to the viewing zenith angle μ 1 : The important assumptions of this approach are that R 1 and R 2 have to be sensitive to the snow grain radius and the illumination angle, and least influenced by atmospheric absorption and scattering (Kokhanovsky et al., 2019).
For more details about the algorithm, the reader is referred to Kokhanovsky et al. (2019Kokhanovsky et al. ( , 2020. The impurity index was introduced by Dumont et al. (2014) and exploits the much higher sensitivity of the visible (VIS) wavelengths to impurity content compared with the NIR spectral region. It is calculated by the ratio of the natural logarithms of green and NIR surface reflectance at 560 and 865 nm, respectively: Dumont et al. (2014) showed that i imp is almost non-sensitive to the ice grain size, whereas it can be affected by atmospheric aerosols in case of biased atmospheric correction results. An accurate surface reflectance retrieval is therefore needed prior to calculating i imp . Furthermore, Di Mauro et al. (2017) demonstrated that i imp is also sensitive to mineral dust and black carbon concentration on ice surfaces. Typical values of the impurity index are 0.2-0.5 for bare ice, 0.7-0.9 for low to moderate chlorophyll content, and more than 0.9 for high chlorophyll concentration (Wang et al., 2020). Its values can reach up to 1.2 for high loads of impurities and cryoconite on bare ice (Di Mauro et al., 2017).

Study Area
Our study area is located at the South-West margin of the Greenland Ice Sheet at 66°−68°N and 48°−50°W. It belongs to the Kangerlussuaq transect (k-transect) and is characterized by patterns of clean snow and dark ice. Especially in the summertime, that is, July and August, the k-transect features a low surface albedo forming a zone of dark ice (Alexander et al., 2014;Ryan et al., 2018). This process is highly correlated with meltwater production and runoff as well as with associated occurrences of algal blooms on the ice surface Cook et al., 2020;Wang et al., 2018). As shown by previous studies, the predominant species of biological impurities during the melt season in the dark zone are Mesotaenium berggrenii and Ancylonema nordenskioldii (Williamson et al., 2018;Yallop et al., 2012). In fact, these eukaryotic species are known to dominate the supraglacial environment both in Greenland and elsewhere . Additionally, the large amount of meltwater production leads to the development of several widespread melt ponds (Diamond et al., 2021).

PRISMA Data
PRISMA is an Italian satellite mission led by the Italian Space Agency (ASI) (Cogliati et al., 2021). The instrument was launched in March 2019 and provides on-demand data for most of the Earth. It features 239 spectral bands covering the wavelength region from 400 to 2,500 nm with a spectral sampling interval (SSI) less than 12 nm. The ground sampling distance (GSD) is 30 m, while the swath is 30 km.
For our study, we selected an acquisition from 30 August 2020, covering a part of the k-transect. Figures 1a  and 1b show a true-color representation of the scene and its location on the Greenland Ice Sheet. The image contains representative examples of both clean snow and dark ice at the end of the melting season. Several melt ponds are also displayed.
After converting PRISMA L1 TOA radiance data to reflectance, we calculated the normalized difference snow index (NDSI) (Dozier, 1989), which is visualized in Figure 1c. We mostly obtain an NDSI beyond 0.8 with 0.74 being the minimum value of the entire image, which clearly indicates that the surface is covered with snow and ice (Dozier & Painter, 2004). We can also observe some smooth structures toward the East showing lower values of NDSI, which might be some thin clouds not easily detectable in the true-color image. Stillinger et al. (2019) have shown that the NDSI of dark clouds can be high enough to cause misclassification. In order to improve the radiometric and spectral quality of the selected PRISMA data, we applied a suite of preprocessing tools, including a spectral smile correction and a radiometric radiance correction (Chlus et al., 2022).
To obtain the individual noise-equivalent change in radiance for each PRISMA spectrum needed by the OE-based inversion, we use an estimation of the signal-to-noise ratio based on a discrete cosine transform and scale the results assuming a photon shot noise square root dependence with the radiance (Gorroño & Guanter, 2021).

Sentinel-3 OLCI Data
OLCI is a moderate resolution imaging spectrometer installed on the Sentinel-3 satellite, which was launched in 2016. The instrument provides 21 spectral bands spanning 400-1,020 nm with an SSI between 2.5 and 40 nm. With 1,270 km and 300 m, it features much larger swath and GSD, respectively, than the PRISMA imaging spectrometer. OLCI was specifically designed for retrieving chlorophyll content, primarily over ocean surfaces, which is highly facilitated by its large footprint (Malenovský et al., 2012).
For the comparison with our PRISMA data set, we selected an OLCI acquisition from the same date, that is, 30 August 2020, and almost the same time of overpass, that is, approximately 15:00 GMT-2. The scene covers large parts of the western shore of the Greenland Ice Sheet and part of the Canadian arctic. It includes our study area in the k-transect of southwest Greenland and shows a slightly larger cloud fraction, which is mainly located over water surfaces though. We used the OLCI L1B product providing radiometrically calibrated TOA radiances and converted the data both to TOA and surface reflectance using the S3 OLCI Snow and Ice Properties Processor (SICE). Details on SICE can be found in Kokhanovsky et al. (2020). Subsequently, we produced a snow grain size map and calculated the impurity index for each pixel using OLCI bands 6 at 560 nm and 17 at 865 nm.

Snow and Ice Parameter Maps
The panels on the left hand side of Figure 2 quantify the spatial distribution of ice grain radius, ice liquid water path length, and glacier algae mass mixing ratio from the PRISMA data using the glacier ice spectral library as prior knowledge.
Comparing the maps with the true-color image in Figure 1b, it is obvious that the darker the surface, the larger are the estimated ice grains and the algae concentration since high amounts of both quantities lead to decreasing reflectance in the VIS . Likewise, the liquid water path length detected on the ice surface is significantly larger for the dark zone in the western part. The algae map is calculated from the sum of retrieved values of the species M. berggrenii and A. nordenskioldii, and conforms to values measured in the field (Cook et al., 2020). Figure 3 illustrates these findings by showing spatial transects of ice grain radius, liquid water path length, and algae mass mixing ratio at 67.14°N (see dashed red lines in the panels on the left hand side of Figure 2).
We selected this particular latitude as this transect not only covers the dark zone and clean ice and snow, but also the large dark melt pond located in the north-eastern part of the image. Between 48.5° and 48.8°W, the transect can generally be characterized as transition area from the dark zone near the coastline toward the clean ice at higher elevated parts of the Ice Sheet. This transition area is interrupted by some small-scale accumulations of algae around melt ponds, which typically cause algal disposition in the surrounding area. In contrast, we observe constantly large ice grain radii and liquid water path lengths as well as high algae concentration within the dark zone. The discrete spike in all transects between 48.8° and 48.9°W originates from a small and shallow melt pond, whose brighter reflectance properties are most likely influenced by underlying bare ice featuring a smaller grain radius and very low algae concentration. On the other hand, the dark melt pond is characterized by large estimated ice grains and high concentration of algae (see Section 4.5 for a detailed analysis). Finally, the region of clean ice and snow shows small grain radii, less liquid water on the ice surface, and almost no biological impurities. Overall, the reported value ranges for the various parameters coincide with findings from previous studies Cook et al., 2020). Especially the comparison with samples of algal field measurements collected and provided by Cook et al. (2020) proves a similar value range of mass mixing ratios remotely retrieved from PRISMA data (see box plot in Figure 3). In fact, the concentrations observed in the field are slightly lower, but this is probably due to an earlier sampling date within the melting season, that is, mid of July instead of late August (acquisition time of the PRISMA data). However, the results from the PRISMA data represent average values of 30 × 30 m pixels instead of point measurements. Thus, the mean of the algal field measurements is the more appropriate quantity to compare with. The panels on the right hand side of Figure 2 present the estimated maps for snow grain radius, snow liquid water path length, and snow algae mass mixing ratio using the Lorenz-Mie based snow spectral library as prior knowledge. In contrast to the retrieved ice grain size, a correlation with surface brightness is not observable for the snow grain radius. In fact, this retrieval ideally works for sphereshaped snow grains, so that the reported values for the dark ice surface have to be treated carefully. Toward the most eastern part, the map features smaller grain radii potentially related to the increasing surface elevation, which rises from 1,000 to 1,500 m above sea level in our PRISMA image and leads to lower air temperatures when moving landwards. Under these conditions, generally dry snow with small grain size is found on the surface (Warren, 2019). Several studies well describe the spatial distribution of snow grain size including its decline on the uplifted parts of the Greenland Ice Sheet (see e.g., Kokhanovsky et al. (2019) or Bohn et al. (2021)). The estimated snow liquid water path lengths confirm the retrieved snow grain size map since the highest values can be observed for pixels with large grain radii of up to 800-1,000 μm. The grains in liquid water enriched wet snow tend to form clusters, which behave as larger grains with the respective optical properties (Dozier et al., 2009). Finally, the snow algae map in Figure 2f points out the importance of selecting appropriate priors for the inversion. Applying the snow spectral library to retrievals on glacier ice surfaces obviously leads to less realistic results of algae concentration. The estimated mass mixing ratios do not correlate with surface brightness and show artificially high values when compared to the field observations of Cook et al. (2020) (see Figure 3). This may be due to the different approaches of modeling the shape of both snow and ice grains and the algal cells. The size of the latter might be an even bigger source of uncertainty than the shape by acting as an equivalent radius. However, relying on prior knowledge based on GO calculations significantly enhances the retrieval results since it simulates existing conditions on glacier ice surfaces more appropriately (Cook et al., 2020).

Posterior Error Correlation
We present posterior error correlation matrices for selected atmosphere and surface parameters to show how retrieval uncertainties of particular state vector elements affect each other. We calculated the mean coefficients from the posterior predictive uncertainties for all ̂ of the PRISMA image. Depending on the used surface prior spectral library, Figures 2g and 2h divides into glacier ice and snow surface parameters.
Although we do not analyze the retrieval of the atmospheric state x ATM in this study, we take a look at potential effects of water vapor and AOT on the additional surface parameters. Whereas water vapor uncertainties are clearly uncorrelated with all other parameters over glacier ice, a negative correlation between errors in the ice grain retrieval and the AOT estimation can be observed. Scattering and absorption by atmospheric aerosols show similar effects on the reflectance shape and magnitude in the VIS as increasing ice grain radii. Thus, corresponding retrieval uncertainties are introduced, which was one of the key findings in Bohn et al. (2021). Furthermore, posterior errors for the glacier algae species are strongly negatively correlated since their absorption features are similar (Cook et al., 2020). However, we report the sum of both in the retrieval maps, so that potential inaccuracies compensate for each other. Figure 2h illustrates the positive correlation between uncertainties in the snow grain size retrieval and errors in the liquid water estimation. This is most likely due to the similarities between liquid and ice absorption shapes (Green et al., 2006). Even errors in the solution state for atmospheric water vapor can be little affected by posterior uncertainties in the surface liquid water estimation. However, the respective correlation coefficient is only −0.13 and likewise, the remaining values of the matrix are more or less close to 0. Overall, Figures 2g and 2h confirm the independence of most state vector parameters and therefore, our ability to estimate them with the "lazy Gaussian" inversion.

Top-of-Atmosphere Radiance Fits
We present a comparison between PRISMA L1 data and the respective TOA radiance fits, modeled by Equation 2. As an example, the upper panels of Figure 4 show three selected spectra of different Ice Sheet surface types as highlighted in Figure 3.
The left panel represents a clean snow surface in the eastern part of the image featuring small snow grains, a smaller liquid water path length, and no algae accumulation. In contrast, the spectrum in the middle panel originates from a dark ice pixel in the ablation zone having large ice crystals, a large liquid water path length, and a high glacier algae mass mixing ratio. Finally, the right panel emphasizes the radiative and reflective properties of the dark melt pond located on the spatial transect drawn in the panels on the left hand side of Figure 2. While showing almost no residuals in the SWIR, all spectral fits illustrate discrepancies of up to 2 μWnm −1 cm −2 sr −1 in the VIS/NIR wavelength region. Generally, the modeled TOA radiance rather underestimates the measured PRISMA L1 data, except for the NIR part of the melt pond spectrum. However, we observe slightly different spectral regions of largest error occurrence. The radiance fit for the dark ice surface almost exclusively deviates from the PRISMA measurement between 400 and 750 nm, where the TOA radiance signal is strongly affected by the scattering of atmospheric aerosols. An explanation is directly presented in Figure 2g. Here, we notice a correlation coefficient of −0.69 between errors in the ice grain radius retrieval and the AOT estimation. Therefore, the AOT value reported in the solution state for the dark ice spectrum might be overestimated due to an underestimated ice grain radius. This reduces the modeled radiance in the VIS. Additionally, the AOT estimation is biased by a missing first guess retrieval prior to the inversion.
The fit for the clean snow spectrum shows less influences by the AOT retrieval in the VIS. Here, we observe the largest model discrepancies in the NIR wavelength region. As the inversion reports a much smaller snow grain radius, but remarkably higher relative liquid water fraction, the residuals might be explained by error correlation in-between the three phases of water, that is, atmospheric water vapor, surface liquid water, and snow grain radius. Figure 2h confirms this assumption since we note correlation coefficients of 0.60 between snow grain and liquid water retrieval uncertainties as well as at least −0.13 for errors in water vapor estimation and the reported liquid water fraction.
Finally, the radiance fit for the melt pond spectrum slightly deviates from the PRISMA L1 data in the blue VIS region, but shows larger differences in the NIR wavelength range. The former is most likely caused by uncertainties in the AOT estimation, while the latter might be explained by insufficient surface prior knowledge. The applied spectral libraries of snow and ice reflectance do not include simulations for melt pond surfaces and consequentially, the prior state vector does not cover these characteristics in the inversion. This is also reflected in the estimated ice crystal size for this spectrum. The inversion reports a disproportionately large radius of 23,212 μm, although we rather find open water than ice-covered surface in this pixel. Here, the solution state of the ice crystal size is clearly guided by the relatively low radiance, which is commonly observed for water surfaces.
Overall, the discrepancies in modeled TOA radiance may also originate from too strong constraints on the surface reflectance priors. The optimization then attends less to the measurement part of the cost function and consequentially, models y with a higher associated uncertainty. Increasing the surface reflectance diagonal of the prior covariance matrix may improve the performance of our forward model. Also, uncertainties introduced by the radiometric calibration of the instrument itself might be another source of errors influencing the TOA radiance fits.
Finally, we presume that at least the amount of algae accumulation on the ice surface has less effects on the fitted TOA radiance. Bohn et al. (2021) have shown that the information content of the radiance measurement is almost unaffected by biological impurities. However, small errors might still remain in the TOA radiance fits.

Estimated Surface Reflectance
Since the "lazy Gaussian" inversion is embedded in an atmospheric correction algorithm and the spectral albedo for each instrument channel are elements of the state vector, the evaluation of the retrieved surface reflectance is an essential part of our analysis. Although we lack appropriate field measurements for validation, a qualitative comparison with the official PRISMA L2C product is informative. PRISMA provides two surface reflectance products (L2) that are either in sensor geometry (L2C) or in orthorectified map geometry (L2D). It is an independent estimate of the surface reflectance using a traditional sequential atmospheric correction algorithm based on the MODTRAN radiative transfer code and does include both a water vapor and AOT estimation, but no additional surface properties (ASI, 2020). Since our resulting reflectance map is in similar sensor geometry as the PRISMA L1 product, we use the L2C data for comparison instead of the final orthorectified L2D product.
The lower panel of Figure 4 shows results for the same pixels as analyzed in Section 4.3. For clarity, we excluded reflectance values from instrument channels located within the deep SWIR water vapor absorption features around 1,350 and 1,850 nm, where the solar radiation is almost entirely absorbed by the atmosphere. Even marginally biased simulations of atmospheric water vapor transmission could lead to artificially high reflectance values at these wavelengths. Again, we evaluate spectra of clean snow, dark ice, and a melt pond surface. Overall, we see a good agreement with PRISMA L2C spectra. The results from the "lazy Gaussian" inversion feature less spikes and a smoother reflectance gradient especially in the VIS. This emphasizes the capabilities of OE, which enables a less noisy reflectance estimation by incorporating the prior distribution in the surface model (Thompson et al., 2018). However, all spectra show deviations from the PRISMA data in the same spectral ranges as illustrated by the upper panel of Figure 4. This confirms the assumptions of the previous Section 4.3. On the other hand, further studies are needed to assess the quality of PRISMA L2C spectra and if they can serve as validation targets (Cogliati et al., 2021). Instead, an accurate evaluation of the retrieval results from the "lazy Gaussian" inversion would require field measurements of surface reflectance. Figure 5 shows selected melt pond reflectance spectra representing different water types.

Melt Ponds
Additionally, the estimated glacier algae accumulation for the respective pixels is given in the plot. When comparing with snow or ice surfaces, the reflectance spectrum of melt ponds is characterized by a missing peak at 1,100 nm. The reflectance beyond 900 nm is typically low due to strong liquid water absorption in these wavelengths, with any signal due only to Fresnel reflection (Malinka et al., 2018). Spectra (a) and (b) in Figure 5 only show a marginal peak in the NIR indicating a melt pond without ice cover. Shape and magnitude of both spectra conform with field spectrometer measurements of dark and light-blue ponds presented in Malinka et al. (2018). However, while the inversion reports no present algae accumulation for spectrum (a), the estimated mass mixing ratio of 71 μg/g ice is comparatively high for spectrum (b). Here, we most likely observe the influence of cryoconite on the bottom of the pond, which has been interspersed with melt water. In contrast, spectra (c) and (d) exhibit absorption features in the VIS spectral region caused by abundance of biological impurities on the surface. This assumption is confirmed by retrieved glacier algae mass mixing ratios of 38 and 154 μg/g ice , respectively. Even a distinction between different species of algae is enabled by the retrieval result since both spectra hold different characteristic absorption features. The plots on the right hand side of Figure 5 present a closer look at carotenoid and chlorophyll absorption between 400 and 700 nm present in spectra (c) and (d) in the middle panel of Figure 5. We observe a mixture of phycoerythrin and chlorophyll absorption around 620 nm in spectrum (c) (Bryant, 1982), pointing to green algae or blue colored cyanobacteria, which are commonly found on the Greenland Ice Sheet (Di Gray et al., 2020;Wientjes et al., 2011;Yallop et al., 2012). In contrast, spectrum (d) can be distinguished by a broad carotenoid feature around 500 nm indicating the presence of red or purple algae (Hoham & Remias, 2020). They are found in large quantities on the Greenland Ice Sheet (Cook et al., 2020), which is underlined by the relatively high retrieved concentration of 154 μg/g ice . Present reflectance peaks at 1,100 nm in spectra (c) and (d) suggest though that the respective ponds seem to be either partly covered with ice or to consist of a mixture of water and ice grains (Malinka et al., 2018). This is further endorsed since both spectra (c) and (d) resemble the shape of spectrum (e), which is retrieved from a frozen pond featuring almost clean ice without significant algae accumulation.
Overall, the results demonstrate that the "lazy Gaussian" inversion is able to report meaningful results from PRISMA data for glacial melt ponds by quantifying different brightness of water surfaces, distinguishing turbid and clear water as well as identifying potential ice cover. Furthermore, we show that even weak chlorophyll absorption can be resolved by PRISMA data. To our knowledge, this is the first time that this small absorption is observed from a spaceborne imaging spectrometer, which opens a valuable perspective for the life detection on snow and ice using imaging spectroscopy data.

Comparison With Sentinel-3 OLCI
We present results from the S3 OLCI snow properties retrieval and show a comparison with the PRISMA retrieval maps. In particular, we demonstrate the potential of snow and ice surface parameters derived from imaging spectrometers to develop regression models for multispectral data.  Figure 6 shows S3 OLCI snow grain radius and impurity index calculated according to Equations 5 and 9, respectively.
It is important to note that the OLCI grain size algorithm assumes a spherical grain shape, so that the retrieval rather reports radii of snow grains than dimensions of arbitrarily shaped ice crystals (Kokhanovsky et al., 2019). We masked out non-snow covered pixels to save processing time and complemented the plot with a true-color image of the S3 acquisition. When looking at the eastern part of the scene, we observe a distinct spatial pattern of both parameters having the largest values toward the edge of the ice sheet in a stripe parallel to the coastline. Moving landwards, snow grain radius and impurity index significantly decline. Both their value range and spatial distribution coincide with reported values in, for example, Kokhanovsky et al. (2019) or Wang et al. (2020), and are in line with the seasonal conditions to be found at the end of the melting season in late August (Alexander et al., 2014). As a next step, we generated spatial subsets from the S3 OLCI retrieval maps to match the geographic extent of the PRISMA acquisition. Figure 7 shows a visual comparison of retrieved snow grain radius from both instruments as well as the S3 impurity index and estimated PRISMA glacier algae concentration.
First of all, the maps derived from PRISMA data reveal finer spatial structures and patterns on the surface due to the much smaller GSD. Nevertheless, both distribution and value range of snow grain radius are very similar. We observe a broader stripe of larger radii of up to 1,000 μm extending from North to South in the eastern part of the image and a distinct decrease toward the most north-eastern corner with values of around 200 μm. The impurity index likewise follows the spatial distribution of retrieved glacier algae accumulation. However, the PRISMA glacier algae map yields a clearer distinction of high algae accumulation spots, which is especially demonstrated by the patterns in the middle of the image with mass mixing ratios of up to 160 μg/g ice , and the large melt pond toward the North showing algae concentration both on the water surface and at the shoreline. It is important to note that the impurity index is not only sensitive to biological impurities but also to inorganic LAP such as mineral dust, black carbon, and cryoconite (Di Mauro et al., 2017;Dumont et al., 2014;Wang et al., 2020). Consequentially, deposits of these particles on the ice surface might influence the value of i imp , and thus, explain a part of the variability in the comparison. Moreover, this underlines the importance of estimating the broadband albedo, which characterizes the energy balance of the surface regardless of the type of LAP.
We assess the before-mentioned spatial correlation of S3 and PRISMA snow grain radius as well as impurity index and algae concentration by showing scatter plots in Figures 7c and 7f. To enable a per pixel comparison, we resampled the PRISMA surface parameter maps to 300 m GSD by taking the mean values of 10 × 10 pixel aggregates. Estimated snow grain radii show a remarkable consistency. While we achieve an R 2 of 0.61 and an RMSE of the square roots of 1.79 μm, the values retrieved from multispectral S3 data spread over a larger range reaching 1,000 μm. In contrast, the estimated grain radii for the most north-eastern part of the image are much smaller when applying our proposed approach to the PRISMA data. Here, we observe values even lower than 200 μm. The impurity index seems to be less correlated with glacier algae mass mixing ratio, although featuring an R 2 of 0.76. It is obvious that most of the correlation is influenced by two clusters in the scatter plot, one at i imp around 0.6-0.7 and mass mixing ratios of 100-140 μg/g ice , and another at concentrations below 40 μg/g ice with corresponding i imp of 0.2-0.5. When only considering high glacier algae mass mixing ratio, the impurity index does not significantly increase and remains almost constant at values of around 0.7. This is an indicator that i imp is in fact able to detect algae accumulation on the ice surface, but is less appropriate for describing fine-scale variations of higher amounts of concentration (Wang et al., 2020). Finally, the scattering of points in both plots may also be due to a geometric mismatch, so that a correction for geolocation of the PRISMA image may improve the regression. However, our results demonstrate sufficient potential of the correlation between impurity index and glacier algae mass mixing ratio derived from PRISMA spectra to build predictors for S3 OLCI data. Figure 8 presents predicted glacier algae concentrations for the S3 OLCI acquisition using two different regression methods.
First, we applied the linear regression derived from Figure 7f, y = 227.04x − 41.43, to each pixel of the S3 OLCI image. Then, we fit a Gaussian process regressor (GPR) with a constant kernel to the data from the subset and predicted the glacier algae mass mixing ratios for the complete data set. We selected these two regression approaches as examples for both a simple and a more complex method in order to show the manifold choice of well performing algorithms in the field of supervised learning. Figures 8c and 8d illustrate the performance of both regressors when applied to the S3 OLCI subset covering the same extent as the PRISMA image. We observe almost identical R 2 values of 0.76 and 0.75, respectively, with a larger RMSE of 36.12 μg/g ice for the GPR. Furthermore, the Gaussian kernel densities suggest that a larger fraction of the values predicted by the linear regression is located on the 1:1-line. The respective prediction maps in Figures 8a and 8b indicate that both methods are able to locate the dark zone of high glacier algae accumulation at the edge of the ice shield. However, the linear regression leads to smoother transitions toward lower concentrations, whereas the GPR can better reproduce high amounts of algae on the surface. Nevertheless, for GPR prediction quality, learning the kernel is critically important, and the results could be improved by a detailed investigation and a careful selection of the covariance function and the optimizer of the kernel parameters (Rasmussen & Williams, 2006). Overall, our results provide a promising basis for future exploitation of spectroscopic retrievals to be used as predictors for multispectral data. Different instrument revisit times and the possibility to use imaging spectroscopy data for re-calibration purposes of multiband sensors are other potential synergies. However, a detailed analysis of uncertainty quantification would require concurrent field measurements for a validation of estimated quantities of ice surface parameters.

Scaling to a Global Cryosphere Product
With the setup presented in this study, the "lazy Gaussian" inversion can appropriately be applied to snow and ice surfaces without significant topographic characteristics under sufficient illumination conditions, that is, solar zenith angles not significantly exceeding 50°-60° . This holds true for many parts of the Greenland Ice Sheet during summertime. However, forthcoming orbital imaging spectroscopy missions will deliver high-resolution data on a global scale, which requests for independently applicable retrieval algorithms (Cawse-Nicholson et al., 2021). Especially the SBG designated observable and ESA's CHIME mission are expected to record large data volumes covering a wide range of different snow and ice surface conditions spanning over almost all latitudes.
The results from PRISMA data demonstrate that our approach for mapping snow and ice surface properties hasthe potential for providing a robust cryosphere product based on orbital imaging spectroscopy. However, themethod still faces some challenges that need to be confronted prior to a global application. So far, the inversion uses simulations of spectral albedo by a two-stream snow and ice radiative transfer model as prior knowledge, which does not account for directional effects in the reflectance. Likewise, geometric characteristics ofthe surface such as slope, aspect, sky view factor, or shadow fraction are not considered in the forward model.In order to achieve accurate retrieval results over mountainous areas with complex terrain as well as varyingillumination and observation geometries, simulations of directional reflectance based on multi-stream radiativetransfer models such as DISORT have to be considered as prior knowledge (Lamare et al., 2020). Furthermore, Equation 2 needs to be extended by some additional terms accounting for surface topography. However, by applying an OE-based simultaneous atmosphere and surface inversion scheme our approach provides the basis for a straightforward implementation of these requirements. This will enable a global mapping of snow and ice surface properties corrected for latitudinal and topographic biases including a rigorous quantification of uncertainties.

Conclusion
We present first results from the recently introduced "lazy Gaussian" inversion to infer glacier ice surface properties from a PRISMA imaging spectroscopy data set acquired over the Greenland Ice Sheet. It is the first time that PRISMA data are used for studying the cryosphere and it serves as a precursor to the global availability of spaceborne imaging spectroscopy data, which will allow to detect and quantify snow and ice variables with unprecedented accuracy. The algorithm maps grain radius, liquid water path length, and algae mass mixing ratio, and reports associated posterior predictive uncertainties. Additionally, we show a comparison with multispectral data from the S3 OLCI instrument to detect potential synergies and to reveal how these data can be complimented by satellite spectroscopy observations.
Our results demonstrate that spectroscopic observations from space will play a crucial rule in upcoming research of the Greenland Ice Sheet. We show that these data can be used to detect and quantify patterns of LAP accumulated on the surface in areas such as the dark zone or k-transect. Maps of algae accumulation, surface liquid water, and melt pond evolution provided on a regular basis can support the ongoing investigations of ice sheet melt processes and the resulting sea level rise.
Furthermore, we evidence that glacier algae maps derived from the PRISMA imaging spectrometer can be used to predict the same surface parameter from simple multiband indices such as the impurity index. This opens new possibilities of producing multi-year time series of glacier algae mapping on the Greenland Ice Sheet based on multispectral data sets acquired by instruments such as Landsat or Sentinel-2 and 3. High-frequency observations may not be possible even from the next generation of imaging spectrometers due to their global charter and the high fraction of cloud cover over the Arctic. In contrast, multiband sensors like Sentinel-3 have far greater temporal coverage, but lack imaging spectrometer's sensitivity to subtler snow and ice parameters. Under such circumstances, a hybrid approach can capture the best of both, with sparse imaging spectroscopy data being used to build local models for a more complete interpretation of the multiband data. At the same time, this can fill the gap of missing spectroscopic observations from space during the past four decades. A multitude of upcoming missions such as EnMAP, EMIT, CHIME, and SBG will lead to an unprecedented availability of high-resolution data on a global scale and daily, and thus, improve our understanding of snow and ice surface processes and facilitate the monitoring of glacier ice changes over time.