# Retrieval of Snow Properties from the Sentinel-3 Ocean and Land Colour Instrument

Alexander Kokhanovsky, Maxim Lamare, Olaf Danne, Carsten Brockmann, Marie Dumont, Ghislain Picard, Laurent Arnaud, Vincent Favier, Bruno Jourdain, Emmanuel Le Meur, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Alexander Kokhanovsky, Maxim Lamare, Olaf Danne, Carsten Brockmann, Marie Dumont, et al. Retrieval of Snow Properties from the Sentinel-3 Ocean and Land Colour Instrument. Remote Sensing, 2019, 11 (19), pp.2280. 10.3390/rs11192280. meteo-03657917

HAL Id: meteo-03657917

https://meteofrance.hal.science/meteo-03657917

Submitted on 3 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Article

# Retrieval of Snow Properties from the Sentinel-3 Ocean and Land Colour Instrument

Alexander Kokhanovsky ¹*, Maxim Lamare ²,³, Olaf Danne ⁴, Carsten Brockmann ⁴,†,
Marie Dumont ³, Ghislain Picard ²,†, Laurent Arnaud ²,†, Vincent Favier ²,†, Bruno Jourdain ²,
Emmanuel Le Meur ², Biagio Di Mauro ⁵, Teruo Aoki ⁶,⁷,†, Masashi Niwano ⁷,†,
Vladimir Rozanov ⁸, Sergey Korkin ⁹, Sepp Kipfstuhl ¹⁰, Johannes Freitag ¹⁰, Maria Hoerhold ¹⁰,
Alexandra Zuhr ¹⁰,¹¹, Diana Vladimirova ¹², Anne-Katrine Faber ¹³,†,
Hans Christian Steen-Larsen ¹³,†, Sonja Wahl ¹³,†, Jonas K. Andersen ¹⁴, Baptiste Vandecrux ¹⁴,†,
Dirk van As ¹⁴,†, Kenneth D. Mankoff ¹⁴,†, Michael Kern ¹⁵, Eleonora Zege ¹⁶ and Jason E. Box ¹⁴

¹ VITROCISET Belgium SPRL, Bratustrasse 7, D-64293 Darmstadt, Germany

² Institut des Géosciences de l'Environnement (IGE), Université Grenoble Alpes, CNRS, UMR 5001, 38041 Grenoble, France; maxim.lamare@meteo.fr (M.L.); ghislain.picard@univ-grenoble-alpes.fr (G.P.); laurent.arnaud@univ-grenoble-alpes.fr (L.A.); vincent.favier@univ-grenoble-alpes.fr (V.F.); bruno.jourdain@univ-grenoble-alpes.fr (B.J.); emmanuel.lemeur@univ-grenoble-alpes.fr (E.L.M.)

³ Université Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Etudes de la Neige, 38400 Grenoble, France; Marie.Dumont@meteo.fr

⁴ Brockmann Consult, Max Planck Strasse 2, 21502 Geesthacht, Germany; olaf.danne@brockmann-consult.de (O.D.); Carsten.brockmann@brockmann-consult.de (C.B.)

⁵ Department of Earth and Environmental Sciences, University of Milano-Bicocca, Piazza della Scienza, 1 20126 Milan, Italy; biagio.dimauro@unimib.it

⁶ Arctic Environment Research Center, National Institute of Polar Research, Tachikawa 190-8518, Japan; aoki.teruo@nipr.ac.jp

⁷ Meteorological Research Institute, Japan Meteorological Agency, Tsukuba 305-0052, Japan; mniwano@mri-jma.go.jp

⁸ Institute of Environmental Physics, University of Bremen, D-28359 Bremen, Germany; rozanov@iup.physik.uni-bremen.de

⁹ USRA GESTAR, 7178 Columbia Gateway Drive, Columbia, MD 21046, USA; skorkin@usra.edu

¹⁰ Alfred-Wegener-Institute, 27570 Bremerhaven, Germany; Sepp.Kipfstuhl@awi.de (S.K.); johannes.freitag@awi.de (J.F.); Maria.Hoerhold@awi.de (M.H.); alexandra.zuhr@awi.de (A.Z.)

¹¹ Institute of Earth and Environmental Sciences, University of Potsdam, 14476 Potsdam, Germany

¹² Center for Ice and Climate, Copenhagen University, 1165 Copenhagen, Denmark; diana.vladimirova@nbi.ku.dk

¹³ Geophysical Institute and Bjerknes Centre for Climate Research, University of Bergen, 5007 Bergen, Norway; akfaber@uib.no (A.-K.F.); hans.christian.steen-larsen@uib.no (H.C.S.-L.); Sonja.Wahl@uib.no (S.W.)

¹⁴ Geological Survey of Denmark and Greenland (GEUS), 1350 Copenhagen, Denmark; jkvand@space.dtu.dk (J.K.A.); bav@geus.dk (B.V.); dva@geus.dk (D.v.A.); kdm@geus.dk (K.D.M.); jeb@geus.dk (J.E.B.)

¹⁵ ESTEC/ESA, Science, Applications and Climate Department (EOP-SME), 2201 Noordwijk, The Netherlands; Michael.Kern@esa.int

¹⁶ Institute of Physics, National Academy of Sciences of Belarus, 220072 Minsk, Belarus; ellazege@gmail.com

<figure><img src="image_4.png" /></figure>

check for
updates

Received: 19 July 2019; Accepted: 16 September 2019; Published: 29 September 2019

**Abstract:** The Sentinel Application Platform (SNAP) architecture facilitates Earth Observation data processing. In this work, we present results from a new Snow Processor for SNAP. We also describe physical principles behind the developed snow property retrieval technique based on the analysis of Ocean and Land Colour Instrument (OLCI) onboard Sentinel-3A/B measurements over clean and polluted snow fields. Using OLCI spectral reflectance measurements in the range 400–1020 nm, we derived important snow properties such as spectral and broadband albedo, snow specific surface

<footnote>* Correspondence: a.kokhanovsky@vitrocisetbelgium.com</footnote>
area, snow extent and grain size on a spatial grid of 300 m. The algorithm also incorporated cloud
screening and atmospheric correction procedures over snow surfaces. We present validation results
using ground measurements from Antarctica, the Greenland ice sheet and the French Alps. We find
the spectral albedo retrieved with accuracy of better than 3% on average, making our retrievals
sufficient for a variety of applications. Broadband albedo is retrieved with the average accuracy
of about 5% over snow. Therefore, the uncertainties of satellite retrievals are close to experimental
errors of ground measurements. The retrieved surface grain size shows good agreement with ground
observations. Snow specific surface area observations are also consistent with our OLCI retrievals.
We present snow albedo and grain size mapping over the inland ice sheet of Greenland for areas
including dry snow, melted/melting snow and impurity rich bare ice. The algorithm can be applied
to OLCI Sentinel-3 measurements providing an opportunity for creation of long-term snow property
records essential for climate monitoring and data assimilation studies—especially in the Arctic region,
where we face rapid environmental changes including reduction of snow/ice extent and, therefore,
planetary albedo.

**Keywords:** snow characteristics; optical remote sensing; snow grain size; specific surface area; albedo; Sentinel 3; OLCI

**1. Introduction**

Snow is an essential component of the cryosphere, and plays an important role in the hydrological cycle and also in climate change studies [1]. In particular, it has been found that the global snow cover and snow albedo are decreasing with time [1]. If this trend persists, this may lead to important climatic consequences in the future. In particular, just 1% reduction in snow albedo may lead to a diurnal radiative forcing of 4 W m⁻² locally, similar to that of carbon dioxide doubling [2]. This leads to a substantial climate warming as discussed in [1].

Therefore, it is very useful to develop ground-based, airborne and satellite snow remote sensing
techniques to monitor the evolution of snow surfaces and their properties: snow extent, snow depth,
snow water equivalent, snow specific surface area, snowmelt areas, snow grain size, pollution load
and albedo [3,4]. The World Meteorological Organization has defined snow albedo and snow extent
as essential climate variables (ECVs) that must be monitored on a daily basis globally using satellite
observations. The retrieval of snow optical properties is currently done using various space-borne
instrumentation including optical sensors such as MODerate-resolution Imaging Spectroradiometer
(MODIS) [5–12], Global Imager (GLI) [13,14], Medium Resolution Imaging Spectrometer (MERIS) [15],
Visible Infrared Imaging Radiometer Suite (VIIRS) [16], MultiSpectral Imager (MSI)/Sentinel-2
(S-2) [17], etc.

In this study, we propose a new technique to determine snow optical and microphysical properties
using Ocean and Land Colour Instrument (OLCI) on board ESA Sentinel-3 mission [18]. In particular,
we provide snow products, such as snow specific surface area, snow grain size, snow extent, snow
spectral and broadband albedo.

Currently, two Sentinel-3 satellites are in operation, S3A and S3B, launched, respectively, 16 February 2016 and 25 April 2018, each with design lifespans of seven years. ESA has planned to launch additional Sentinel-3 satellites (S3C, S3D) with identical OLCI instruments on board during next five years (2021 and 2023). Therefore, there is now the possibility to create long-term ECV records of snow properties using identical and inter-calibrated spaceborne instrumentation for at least 15 years with the current constellation of S3A and S3B satellites. The swath of OLCI is 1270 km. The instrument has 21 channels in the spectral range 400 nm to 1020 nm with spatial resolution of 300 m on ground and is suitable for monitoring snow extent/albedo and also other snow properties using backscattered sunlight.
## 2. Materials and Methods

### 2.1. Theory

Snow is composed of transformed ice crystals that precipitate from the atmosphere. The crystals accumulate on the ground and metamorphose in place. The transformed snow may ultimately melt, slide and sublimate. The volumetric concentration of ice in snow is about 30%, although this number varies depending on the snow type and age. The effective grain size (diameter) *d* of crystals in snow is much larger as compared to the wavelength *λ* of incident radiation in the visible and near infrared regions of electromagnetic spectrum. Therefore, the geometrical optics approach is suitable for the description of light scattering, extinction and absorption in a snowpack. The processes of light reflection and refraction occur at the surfaces of crystals. Additionally, light absorption by bulk ice occurs after penetration of photons to the interior of ice crystals. Through a semi-infinite snow layer approximation, we assume that snowpack is so optically thick that the influence of the underlying surface can be neglected. This is a valid assumption in the visible part of the electromagnetic spectrum for snow layers having a thickness larger than about 10–30 cm, depending on snow type, pollution level and underlying surface type. The penetration depth is much smaller for the near-infrared radiation (below 2–5 cm depending on the wavelength) mostly used for the retrievals discussed in this paper. This means that we retrieved snow surface parameters (e.g., grain size at the top of snow layer). We also ignored vertical and horizontal inhomogeneities of snow variables and assumed that the snow is a plane-parallel medium. The retrieval of properties of snowpack covering the tilted surfaces is considered in the framework of the simplified theory.

In our theory, the reflection function of snow depends on just two inherent snow properties: single scattering albedo (SSA) $ω_0$ (or, alternatively, the probability of photon absorption (PPA) ($β = 1 - ω_0$)) and the probability of photon scattering (phase function $p(θ)$) in a given direction (specified by the scattering angle $θ$) by a unit volume of snowpack [19,20]. The phase function $p(θ)$ depends on the scattering angle being larger in the forward scattering region, where the diffraction of light on ice crystals also takes place. In the visible spectral domain and in the absence of pollutants in snow, the single scattering albedo is spectrally neutral and close to one due to extremely weak absorption of radiation by bulk ice [21,22]. This explains the brightness and white color of freshly accumulated dry snow. The theoretical value of fresh clean snow albedo (planar and spherical) is close to one in the visible wavelengths. Its dependence on the size and shape of ice crystals can be neglected in the visible. For instance, it has been found from spectral radiance measurements [23] that the snow albedo in Antarctica is in the range 0.96–0.98 in the visible range of the electromagnetic spectrum almost independently of the snow grain size, shape and illumination conditions. The small difference of albedo from one (0.02–0.04) can be attributed to the snow contamination and measurement errors. Indeed, in the absence of absorption processes, all photons penetrating through the snowpack surface must return back toward the atmosphere and outer space, resulting in albedo equal to 1.0. This is true independently on the type of local scattering law as long as the medium can be treated as semi-infinite. The fresh snow bidirectional reflection distribution function (BRDF) does not depend on the grain size in the visible wavelengths because the light scattering takes place in the geometrical optics domain and the absorption processes by bulk ice in the visible are weak. Although it depends on the shape of particles. It is known that the single scattering pattern (phase function) reaches its asymptotic shape as $d/λ→∞$. This asymptotic pattern is different for absorbing and non-absorbing particles. Furthermore, the pattern is influenced by the shape of grains. Therefore, for instance, the assumption of spherical particles cannot be used in snow optics because this assumption leads to BRDFs never observed for snow [24]. Experimental measurements demonstrate that the snow BRDF can be quite accurately modeled by assuming that the snow phase function is flat in the backward hemisphere. The characterization of semi-infinite non-absorbing snow BRDF using theoretical calculations is a complicated matter because the snow is composed of large and irregularly shaped close-packed ice
grains. Internally and externally mixed pollutants (such as mineral dust and black carbon) may be present in a snowpack, strongly affecting its optical properties [25,26].

Let us introduce the snow reflectance $R = \pi * BRDF = \pi I^\uparrow / \mu_0 F_0$, where $I^\uparrow$ is the reflected light intensity, $\mu_0$ is the cosine of the solar zenith angle and $F_0$ is the solar flux incident on the area perpendicular to the solar beam. $R$ is equal to 1.0 for the ideal Lambertian surface with albedo 1.0 at any observation/illumination geometry. Although visible snow is close to the Lambertian surface with albedo 1.0, the deviations from the Lambertian surface model are significant and cannot be ignored [27–29].

An important result of the asymptotic radiative transfer theory (see Appendix A) is the fact that the snow reflectance at no absorption $R_0$ (e.g., in the visible in absence of pollutants) can be related to the same snow (with the same phase function) reflectance with single scattering albedo not equal to 1.0. In particular, it follows in the framework of the exponential approximation [30–32] that

$$R = R_0 \exp(-f y). \qquad (1)$$

where

$$f = \frac{u(\mu_0)u(\mu)}{R_0}. \qquad (2)$$

$$y = 4 \sqrt{\frac{\beta}{3(1-g)}}. \qquad (3)$$

where $g$ is the average cosine of scattering angle in snow, which is defined as

$$g = \frac{1}{2} \int_{0}^{\pi} p(\theta) \cos(\theta) \sin(\theta) d\theta. \qquad (4)$$

The angular function $f$ satisfies to the reciprocity principle and depends on the viewing/observation geometry via the reflectance of non-absorbing snow layer $R_0$ and via the escape function, which can be approximated by the linear function of the cosine of the viewing zenith angle $\mu$ [33]

$$u(\mu) = \frac{3}{7}(1 + 2\mu). \qquad (5)$$

The angular function $f$ also depends on the type of snow/shape of ice crystals (via reflectance/phase function), but not on the snow grain size.

The snow plane albedo $r_p$ can be derived from Equation (1). In particular, it follows by definition that

$$r_p(\mu_0) = 2 \int_0^1 \bar{R}(\mu_0, \mu) \mu \, d\mu. \qquad (6)$$

where

$$\bar{R}(\mu_0, \mu) = \frac{1}{2\pi} \int_0^{2\pi} R(\mu_0, \mu, \phi) d\phi. \qquad (7)$$

where $\phi$ is the relative azimuthal angle. The double integral (6) can be evaluated numerically, but we will use the accurate analytical approximation described below. As $y \to 0$, it follows from Equation (1) that

$$R = R_0 - u(\mu_0)u(\mu)y. \qquad (8)$$

Because $R_0$ is referred to as a reflecatace for a non-absorbing snow, the double integral (6) with $R = R_0$ is equal to the plane albedo of non-absorbing snow. Therefore, this integral is equal to 1.0 by definition. Then, it follows from Equation (8) (see also Equations (5)–(7)) that

$$r_p(\mu_0) = 1 - y u(\mu_0). \qquad (9)$$
We will use the exponential approximation similar to Equation (1)

$$r_p(\mu_0) = \exp(-y u(\mu_0)). \qquad (10)$$

Equation (10) coincides with Equation (9) as $y \to 0$. However, it can be used in a wider range of the parameter $y$ as compared to Equation (9).

The spherical albedo is defined as

$$r_s = 2 \int_0^1 r_p(\mu_0) \mu_0 d\mu_0. \qquad (11)$$

It follows from Equations (5) and (9) that

$$r_s = 1 - y. \qquad (12)$$

Therefore, one derives in the framework of the exponential approximation(see Equations (9) and (10))

$$r_s = \exp(-y). \qquad (13)$$

Equation (12) can be used to make an interpretation of the similarity parameter $y$. Namely, it follows for the absorptance $A$ of the snow layer under the diffuse illumination conditions and small values of the probability of photon absorption

$$A = 1 - r_s \equiv y. \qquad (14)$$

This means that the parameter $y$ coincides with the value of the snow layer absorptance at small values of the probability of photon absorption in snow. The value of $A$ increases as $g$ (see Equations (3) and (14)) increases as one may expect (stronger forward scattering and, therefore, light penetration to a layer, where light can be effectively absorbed).

It follows from the above equations that the reflectance and plane albedo of snow can be expressed via spherical albedo in the approximation under study

$$R = R_0 r_s^f. \qquad (15)$$

and

$$r_p(\mu_0) = r_s^{u(\mu_0)}. \qquad (16)$$

It follows from Equations (5) and (16) that the plane and spherical albedo coincide at the cosine of the solar zenith angle equal to 2/3. Equation (15) allows the estimation of spherical (and, therefore, plane albedo (see Equation (16)) from the reflectance measurements at fixed observation geometry avoiding integration procedure as shown in Equations (6), (7), (11). It should be stressed that most of optical instruments do not sample all observation directions; therefore, the theory presented here gives an important hint on how to avoid poor sampling of angular measurements in the procedures for the determination of snow albedo.

## 2.2. The Accuracy of Exponential Approximation

In this section, we study the accuracy of the presented approximations at various values of the single scattering albedo using the exact radiative transfer calculations. The Henyey-Greenstein phase function

$$p(\eta) = \frac{1 - g^2}{(1 + g^2 - 2g\eta)^{3/2}}. \qquad (17)$$
where η is the cosine of the scattering angle, with the average cosine of scattering angle g = 0.75 similar to that for snow and crystalline clouds, is assumed in all calculations. As a matter of fact, it is not the specific phase function, but rather the average cosine of scattering angle that determines the relationship between R and R₀ to be used in this work (see Equation (1)). We assumed an optical thickness of 5000 in our numerical solution of the integro-differential radiative transfer equation for determination of snow reflectance and in the evaluation of the approximations presented in this work. Such an assumption is needed to have the case of a semi-infinite medium considered in the previous section. In this paper, we made use of Intensity and POLarization calculation radiative transfer code (IPOL [34]). IPOL is a Fortran 90/95, BLAS/LAPACK based, discrete-ordinates matrix-operator vector radiative transfer code that computes intensity and polarization (including ellipticity) of the solar radiation reflected from various turbid media. The code also computes planar and spherical albedos by numerical integration of the azimuthally averaged intensity. Recent intercomparison of radiative transfer codes [35] revealed high accuracy of IPOL, including the case of cloud layers, which are in many respects analogous to snowpack.

The results of inter-comparison of exact and approximate results for the escape function are illustrated in Figure 1. It follows that the approximation can be used for solar zenith angles under 80 degrees. The accuracy for the reflectance is illustrated in Figure 2 at the nadir observation and for various values of single scattering albedo. It follows that the reflectance decreases for the cases of high solar zenith angles. The reflectance increases for the overhead Sun position. It follows that the approximation is highly accurate even for low (0.2) reflectance values (and PPA equal to 0.05), when the dependence of reflectance on the solar zenith angle is low. The errors are also small (see Figure 3) for the plane albedo in the range 0.3–1.0 typical for clean and polluted snow. In Figure 3, the accuracy of two approximations for the plane albedo is shown (via the similarity parameter *y* as shown in Equation (10) and also using Equation (16) with spherical albedo derived from Equation (15)). Both approximations produce acceptable results. The accuracy of approximation for the spherical albedo is illustrated in Figure 4. It follows that the approximation provides accurate results even at comparatively large values of similarity parameter (*y* = 2.5). The summary of relative errors for various parameters is illustrated in Figure 5 for the nadir observation case. The spread along the vertical axis is due to solar zenith angle. We conclude that errors are smaller than 5% for the majority of cases. Therefore, the derived approximations are highly accurate and can be used for the solution of the inverse problem: determination of snow properties using optical measurements. Similar conclusions have been reached using a ray tracing model at the grain scale for a wavelength of 1300 nm [36]. A slight discrepancy appeared only at 1550 nm, where the ice absorption is extremely strong [37].

The indirect confirmation of the applicability of the approximate theory also follows from Figure 6,
where we present the results of plane and spherical albedo calculations for the various values of single
scattering albedo (as in Figure 2). The spherical albedo and planar albedo curves indeed intersect at
μ₀ = 2/3, as suggested by the theory presented above. In this case, the solar zenith angle (SZA) is equal
to 48 degrees.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** The intercomparison of escape functions calculated using exact theory and Equation (5).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** The dependence of the snow reflectance calculated using Henyey-Greenstein phase function at $g = 0.75$ (semi-infinite snow case) and various values of single scattering albedo (SSA) as a function of the solar zenith angle at the nadir observation (exact results: crosses, approximation: boxes).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** The plane albedo as the function of the cosine of the solar zenith angle at various values of single scattering albedo for the semi-infinite snow layer with the asymmetry parameter $g = 0.75$. The exact results are presented together with two approximations (boxes-Equation (13), circles-calculation of plane albedo using Equations (15) and (16)).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** The spherical albedo derived using Equation (13) and exact radiative transfer calculations.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** The relative errors of approximations.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** The dependence of plane albedo on the solar zenith angle according to exact radiative transfer calculations for various values of single scattering albedo (the same as in Figure 5). The spherical albedo is shown by dashed lines. The spherical albedo coincides with the plane albedo at $\mu_0 = 2/3$.</figcaption></figure>
3. Results

3.1. The Approximate Solution of the Inverse Problem: Clean Snow

Let us consider first the case of fresh, dry and unpolluted snow. In the 400 nm to 1020 nm spectral range spanned by the OLCI instrument, one can ignore the ice crystal size dependence of the value of the average cosine of scattering angle g.

The parameter β is the ratio of absorption to extinction coefficients

$$
\beta = \frac{k_{abs}}{k_{ext}}. \qquad (18)
$$

The absorption and extinction coefficients are defined as

$$
k_{abs} = N\langle C_{abs} \rangle, \quad k_{ext} = N\langle C_{ext} \rangle. \tag{19}
$$

where N = c/<V> is the number of grains in unit volume of snow, <V> is the average volume of grains and c is the volumetric concentration of grains. Because grains are much larger than the wavelength of the incident light, one may assume [38]

$$
\langle C_{ext} \rangle = 2\langle \Sigma \rangle. \tag{20}
$$

where <Σ> is the average cross section of grains (perpendicular to incident light). It is equal to <S>/4 for spherical particles and also convex particles in random orientation [38]. Here, <S> is the average surface area of grains. Assuming that grains are weakly absorbing, which is the case in the visible and near-infrared regions of the electromagnetic spectrum, the following approximation holds [19,39,40]

$$
\langle C_{abs} \rangle = B\alpha\langle V \rangle. \tag{21}
$$

where α = 4πχ/λ is the bulk ice absorption coefficient, χ is the imaginary part of ice refractive index and the absorption enhancement factor B depends on the shape of grains but not on the size of grains. The slight dependence of B on the spectral refractive index in the visible and near-infrared can be neglected.

It follows

$$
\beta = q\alpha. \tag{22}
$$

where $q = B<V>/2<Σ>$. Additionally, the following important relationship for the similarity parameter $y$ can be derived using equations given above

$$
y = \sqrt{\alpha l}. \tag{23}
$$

where

$$
l = \frac{q}{\Lambda} \quad \Lambda = \frac{3(1-g)}{16}. \qquad (24)
$$

This equation is highly accurate in the visible and near-infrared wavelengths covered by OLCI.
However, it is less accurate at larger wavelengths, where the spectral variation of the asymmetry
parameter shall be taken into account. The same is true for Equation (21). The extension of equations
to the longer wavelengths in the case of spherical particles has been given in [41] (see also [42]). Taking
into account Equation (1), one derives for OLCI channels

$$
R = R_0 \exp(-f \sqrt{\alpha l}). \qquad (25)
$$

Therefore, one concludes that the spectral snow reflectance in the visible and near infrared is
determined by the spectral bulk ice absorption coefficient and just two additional parameters: the
reflectance of a non-absorbing snow and the absorption length $l$. The value of $f$ in Equation (25) is expressed via $R_0$ and the function given by Equation (5). These two parameters can be determined from the measurements at two wavelengths (say, at the OLCI wavelengths 865 nm ($\lambda_1$) and 1020 nm ($\lambda_2$), see Equation (25)) using [19]

$$R_0 = R_1^\epsilon R_2^{(1-\epsilon)} \qquad (26)$$

$$l = \frac{1}{\alpha_2 f^2} \ln\left(\frac{R_2}{R_1}\right). \qquad (27)$$

where $\epsilon = (1-b)^{-1}$, $b = \sqrt{\alpha_1/\alpha_2}$ and indices signify the wavelengths used. The error budget for the retrieved parameters $R_0$ and $l$ is presented in Appendix B.

The derived value of $l$ can be used to find the spectral PPA with Equation (22) and the relation (see Equation (24))

$$q = \Lambda l. \qquad (28)$$

The determination of the parameters given by Equations (26) and (27) makes it possible to find the spectral reflection function, spectral planar and spherical albedos using equations given above (see Equations (13), (16), (23), (25)) in the broad spectral range 300–2400 nm. This is the main idea of this technique. The main point is that the retrieved albedo values are not influenced by the assumption on the shape of particles. The shape effects are accounted for by the value of $l$ found from the measurements.

The snow broadband albedo (BBA) $r_b$ can be determined using the integration of the spectral albedo. In particular, it follows

$$r_b = \frac{\int_{\lambda_s}^{\lambda_f} r_{p(s)}(\lambda)F(\lambda)}{\int_{\lambda_s}^{\lambda_f} F(\lambda)d\lambda}. \qquad (29)$$

where $F(\lambda)$ is the incident solar flux at the snow surface, $r_{p(s)}$ is plane or spherical albedo depending if plane or spherical BBA is of interest. The indices signify the wavelengths used. We have derived the incident solar flux at the snow surface using the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) code [43] in the 300 nm to 2400 nm range with the solar zenith angle grid of 1 degree and the following assumptions:

* the water vapor column: 2.085 g/m;

* the ozone column: 0.35 atm-cm;

* the tropospheric ozone: 0.0346 atm-cm;

* the aerosol model: rural;

* the vertical optical depth of boundary layer at 550 nm: 0.1;

* the altitude: 825 m;

* and the surface set to snow, calculated in SBDART using the Mie/δ-Edington model [44] with a snow grain diameter of 0.25 mm.

Generally, the results are only weakly sensitive to the variation of the function $F(\lambda)$. Therefore, we have used the same solar flux at the snow surface for the retrievals at different geographical locations. We define the following parameters with Equation (29):

* BBA visible (300–700 nm),

* near IR BBA near IR (700–2400 nm),

* BBA shortwave (300–2400 nm).

BBA visible is the most accurate because the approximations have higher accuracy at small values of the probability of photon absorption. However, the errors of the determination of BBA near IR and BBA shortwave are small as well because the main contribution to integral (29) comes from the visible and near-infrared regions of the electromagnetic spectrum, where derived equations are very accurate.
The values of albedo are retrieved without any assumptions with respect to the particle size or shape. By assuming the snow crystal shape, one may also derive the snow grain size using theory presented above. Namely, let us introduce the effective optical grain diameter

$$d = \frac{3\langle V \rangle}{4\langle \Sigma \rangle}. \qquad (30)$$

Then, the effective diameter can be related to the retrieved value of the absorption length (see Equation (24)) using

$$d = \frac{9l}{16\gamma}. \qquad (31)$$

where $\gamma = B/(1 - g)$ is the so-called scaling constant. It has been found in [39] from the analysis of experimental data of [45] that the scaling constant is in the range 7.4–11.5 and advised to use the scaling constant value of 9.2. We adopt this value in our study. Then, it follows: $d = \sigma l$, where $\sigma \approx 0.06$. The values of the pair $(B, g)$ depend on the type of snow/shape of ice crystals [25,46]. In particular, if one assumes that $B = 1.6$ (as suggested in [40]), one can derive that $1 - g = 0.17$ at $\gamma = 9.2$, which is smaller as compared to the value of $1 - g = 0.25$ used in modern ice cloud remote sensing research [47,48]. This may be due to roundness of grains in snow. Clearly, the retrieved value of $d$ will be biased if the shape of snow grains in the studied snowpack differs from that assumed in the retrieval. Importantly, the dependence on the snow grain shape is influenced by the scaling constant and not by separate values of $B$ and $1 - g$. Clearly, more attention must be given to experimental studies of the variability of scaling constant for various types of snow.

We can also retrieve the snow specific surface area. It is related to the diameter $d$ by the following simple equation [36,37,49]

$$s = \frac{6}{\rho d}. \qquad (32)$$

where $\rho$ is the bulk ice density and it has been assumed that the surface area of ice grains is equal to $4\langle\Sigma\rangle$.

## 3.2. The Determination of Polluted Snow Properties

### 3.2.1. Albedo

In the case of polluted snow, we used the above equations, except that the spectral albedo from 400 nm to 1020 nm was derived directly from OLCI measurements. The determination of PPA also needed to be modified. In particular, the following expression for the spherical albedo was used (see Equation (35))

$$r_{s(\lambda)} = \left[ \frac{R_{BOA}(\lambda)}{R_0} \right]^{1/f}. \qquad (33)$$

where $R_{BOA}(\lambda)$ is the bottom of atmosphere snow reflectance. It is derived from the top-of-atmosphere OLCI reflectance accounting for molecular scattering, water vapor and ozone absorption. Aerosol effects are neglected. Equation (33) is applied at the wavelengths 400, 560 and 1020 nm. The value $R_0$ is derived using equations presented in Appendix C. The planar albedo is found from Equations (16) and (33).

For the BBA calculation, we needed values of spectral albedo in the 300 nm to 2400 nm range. We approximated the derived spherical albedo spectra (see Equation (33)) in the range 300 nm to 1020 nm for polluted snow by the quadratic polynomial using the measurements at 400 nm, 560 nm and 1020 nm least influenced by gaseous absorption (except ozone at 560 nm). The derived polynomial was then used for the determination of PPA for the case of polluted snow (see Equations (3) and (13))

$$\beta(\lambda) = \Lambda \ln^2(r_{s(\lambda)}). \qquad (34)$$
The albedo values at the wavelengths longer than 1020 nm needed for the evaluation of BBA were found using Equations (10) and (13) as for clean snow. This is possible because the influence of pollutants on snow reflectance spectra is less pronounced or negligible in the near-infrared, where ice grains themselves are strong absorbers of radiation.

In the case of clean snow, only the wavelengths 865 nm and 1020 nm were used in our retrieval. Therefore, the errors due to uncertainties of atmospheric correction [50] were minimized. In the case of polluted snow, albedo was derived using the measured OLCI top-of-atmosphere reflectance also at channels 400 and 560 nm, which were influenced by the optical properties of atmosphere to a much greater extent (e.g., molecular and atmospheric aerosol scattering/absorption), see also [20]. It is expected that the derived parameters for the polluted snow have larger errors as compared to those for unpolluted snow due to the problems of atmospheric correction over snow-covered areas.

### 3.2.2. Snow Grain Size and Absorption Coefficient of Pollutants

The snow grain size and absorption coefficient of pollutants in the case of polluted snow is found as explained below. It is assumed that the reflection function of polluted snow is described by Equation (15)

$$R = R_0 r_s^f. \qquad (35)$$

with (see Equation (13))

$$r_s = \exp(-y). \qquad (36)$$

where

$$y = 4 \sqrt{\frac{\beta}{3(1-g)}}. \qquad (37)$$

$$\beta = \frac{k_{abs}}{k_{ext}}. \qquad (38)$$

We shall assume that the scattering and extinction of light by pollutants is much smaller than that of ice grains and absorption coefficient of pollutants can be approximated as

$$k_{abs, pol}(\lambda) = k_{abs, pol}(\lambda = \lambda_0) \left(\frac{\lambda}{\lambda_0}\right)^{-m}. \qquad (39)$$

where $\lambda_0 = 1$ µm. Then, it follows for the polluted snow reflection function

$$R = R_0 \exp \left( -4f \sqrt{\frac{(B\alpha c + \kappa \tilde{\lambda}^{-m})d}{3(1-g)c}} \right). \qquad (40)$$

where $\kappa \equiv k_{abs, pol}(\lambda_0)$, $\tilde{\lambda} = \frac{\lambda}{\lambda_0}$. We have used the geometrical optics result for the snow extinction coefficient $k_{ext} = \frac{3c}{d}$ (see Equations (19), (20) and (30)) and the following relationship for the absorption coefficient of polluted snow in the framework of the impurity external mixture approximation

$$k_{abs} = k_{abs, ice} + k_{abs, pol}. \qquad (41)$$

where (see Equations (19), (21) and (30))

$$k_{abs, ice} = B\alpha c. \qquad (42)$$

Using the absorption length $l$ (see Equation (24)), one can derive from Equation (40)

$$R = R_0 \exp \left( -f \sqrt{(\alpha + \tilde{k} \tilde{\lambda}^{-m})l} \right). \qquad (43)$$
where $\tilde{k} = \frac{k}{B_c}$.

Equation (43) can be used to solve the inverse snow remote sensing problem using spectral reflectance measurements. There are four unknown parameters to be found using Equation (43): $l$, $m$, $R_0$, $\tilde{k}$. They can be found numerically (e.g., using optimal estimation approach) or analytically under the assumption that pollutants do not influence reflectance in the near infrared (say, at 865 and 1020 nm for OLCI). We used the second approach. Then, one can derive [20]

$$R_0 = R_3^\epsilon R_4^{(1-\epsilon)}, \quad l = \frac{1}{\alpha_4 f^2} \ln^2\left(\frac{R_4}{R_0}\right), \quad m = \frac{\ln\left(\frac{p_1}{p_2}\right)}{\ln\left(\frac{\lambda_2}{\lambda_1}\right)}, \quad \tilde{\kappa} = \frac{p_1 \tilde{\lambda}_1^m}{f^2 l}. \quad (44)$$

where

$$\epsilon = \frac{1}{1-b'}, \quad b = \sqrt{\frac{\alpha_3}{\alpha_4}}, \quad p_k = \ln^2\left(\frac{R_k}{R_0}\right). \quad (45)$$

Symbols 1–4 signify the wavelengths used (400, 560, 865 and 1020 nm, respectively, for OLCI). These equations are used to derive the snow grain diameter (see Equation (31)) and also the absorption coefficients of pollutants (see Equation (39)) normalized to the concentration of snow grains (snow density). The error budget for the determination of parameters given by Equations (44) is discussed in Appendix D.

### 3.3. Snow Properties Processor for SNAP

The Snow Properties Processor (SPP) has been developed and implemented under the framework of the European Space Agency (ESA) Sentinel-3 for Science (S34Sci) Land Study 1: Snow that delivered a SNAP plug-in. The processor version 2.3 provided the implementation of the algorithms as described above for the various snow properties. As input, the processor requires an OLCI L1b product that is atmospherically corrected and cloud screened in pre-processing. The output was a set of snow properties, defined by the user via processing parameters. Atmospheric correction was performed ignoring effects of aerosol scattering and absorption. Molecular scattering and absorption was fully accounted for as provided by the SNAP by default.

When the Snow Properties Processor was called from its menu entry in the SNAP desktop application, the processor Graphic User Interface was displayed. It contains two tabs ‘I/O Parameters’ and ‘Processing Parameters’ (Figure 7). Here, I/O stands for Input/Output.

The most important processing parameters, which can be set by the user, are as follows:

* Select OLCI wavelengths for spectral snow quantities: for the selected wavelengths, the computed spectral snow quantities (i.e., spectral snow albedo) were written as the corresponding band into the target product.

* Name of binary mask band in cloud classification product (if present): see details in the software user manual at https://readthedocs.org/projects/s3tbx-snow/.

* Consider snow mask based on the OLCI Normalized-Difference Snow Index (NDSI): if selected, an NDSI value will be computed from atmospherically corrected reflectances at wavelengths 865 nm and 1020 nm (NDSI = (R(865 nm) – R(1020 nm))/(R(865 nm) + R(1020 nm))). If this value exceeded 0.03 (and reflectance at 410 nm > 0.5), the given pixel was regarded as a snow pixel. The normalized bare ice index (NDBI) was also provided in the output. It was defined as NDBI = (R(410 nm) – R(1020 nm))/(R(410 nm) + R(1020 nm)). It was assumed that the values of NDBI below 1/3 corresponded to snow and the ice was assumed for values of NDBI above 1/3. The dark bare areas correspond to values of NDBI above 2/3. OLCI snow mask is given in output (snow mask was equal to 1 for 100% snow pixels. Otherwise, the snow mask was zero).

* Consider snow pollution: if selected, a retrieval for polluted snow was applied for the considered snow pixel.
* Compute PPA: if selected, the spectral probability of photon absorption (PPA) was written to the target product for each selected OLCI wavelength.

* OLCI reference wavelength: reference wavelength used in the snow property algorithms.

* OLCI gain for band n: OLCI system vicarious calibration (SVC) gain for the relevant bands (n = 1, 6, 17, 21 (400, 560, 865, 1020 nm)) used in the retrieval algorithms outlined above.

More details on this can be obtained from the SNAP processor help documentation and from the software user manual (see above). The result product from the snow properties processor contains various snow properties, depending on the processing parameters specified by the user. The bands which can be generated are presented in Table 1.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** The Sentinel Application Platform (SNAP) Snow Properties Processor: Input/Output (I/O) and processing parameters tabs.</figcaption></figure>
<table><caption>**Table 1.** Bands in the final snow properties product.</caption><thead><tr><th>Name in Product</th><th>Units, When Not Dimensionless</th><th>Description</th></tr></thead><tbody><tr><td>albedo_bb_sphericalдан</td><td></td><td>Spherical albedo in the broadband visible range from 300 nm to 700 nm</td></tr><tr><td>albedo_bb_spherical_nir</td><td></td><td>Spherical albedo in the broadband near infrared range from 700 nm to 2400 nm</td></tr><tr><td>albedo_bb_spherical_sw</td><td></td><td>Spherical albedo in the broadband shortwave range from 300 nm to 2400 nm</td></tr><tr><td>albedo_bb_planarдан</td><td></td><td>Planar albedo in the broadband visible range from 300 nm to 700 nm</td></tr><tr><td>albedo_bb_planar_nir</td><td></td><td>Planar albedo in the broadband near infrared range from 700 nm to 2400 nm</td></tr><tr><td>albedo_bb_planar_sw</td><td></td><td>Planar albedo in the broadband shortwave range from 300 nm to 2400 nm</td></tr><tr><td>albedo_spectral_spherical</td><td></td><td>Spectral spherical albedo in 21 OLCI bands</td></tr><tr><td>albedo_spectral_planar</td><td></td><td>Spectral planar albedo in 21 OLCI bands</td></tr><tr><td>rBRR</td><td></td><td>Estimated bottom of atmosphere OLCI reflectance in 21 OLCI bands</td></tr><tr><td>ppa_spectral</td><td></td><td>PPA in 21 OLCI bands</td></tr><tr><td>grain_diameter</td><td>mm</td><td>Snow grain diameter</td></tr><tr><td>snow specific area</td><td>m<sup>2</sup> kg<sup>-1</sup></td><td>Snow specific surface area</td></tr><tr><td>ndbi</td><td></td><td>Bare ice indicator</td></tr><tr><td>pollution_mask</td><td></td><td>Pollution mask</td></tr><tr><td>k̃</td><td>mm<sup>-1</sup></td><td>Normalized snow impurity absorption coefficient at 1 µm</td></tr><tr><td>l</td><td>mm</td><td>Effective absorption length</td></tr><tr><td>m</td><td></td><td>Absorption Angstrom parameter</td></tr><tr><td>r_0</td><td></td><td>Reflectance of non-absorbing snow layer</td></tr><tr><td>f_rel_err</td><td></td><td>Relative error of normalized snow impurity absorption coefficient at 1 µm</td></tr><tr><td>l_rel_err</td><td></td><td>Relative error of parameter l</td></tr><tr><td>m_rel_err</td><td></td><td>Relative error of parameter m</td></tr><tr><td>r_0_rel_err</td><td></td><td>Relative error of parameter r_0</td></tr><tr><td>ndsi</td><td></td><td>Normalised differential snow index</td></tr><tr><td>ndsi_mask</td><td></td><td>NDSI mask for snow identification</td></tr><tr><td>quality_flags</td><td></td><td>L1b quality flags</td></tr><tr><td>pixel_classif_flags</td><td></td><td>Pixel classification flags</td></tr></tbody></table>

### 3.4. Cloud Identification

Cloud screening of Antarctic and Col du Lautaret (French Alps) OLCI scenes used for the SPP validation studies is achieved by rejecting mostly cloudy scenes identified by visual inspection. The cloud cover in the vicinity of the field location (approximately 10 × 10 pixels) was visually assessed using RGB true-color composites generated from OLCI bands 4, 6 and 8, and classified in three categories: cloud-free, partly-cloudy and fully-cloudy. Owing to the availability of the scenes (multiple daily acquisitions), only the cloud-free images were retained for this study.

For Greenland, we have found that the majority of clouds pixels can be identified in OLCI imagery using a simple thresholding of snow grain diameter for values less than 0.1 mm or OLCI band 21
(1020 nm) bottom of atmosphere reflectance exceeding 0.72. Residual cloud effects include cloud shadows and thin clouds.

Cloud conditions are constrained for Greenland validation studies with PROgramme for Monitoring Greenland ICE Sheet (PROMICE) ground observations using a cloud probability index (CPI) estimated from PROMICE automatic weather station longwave downward irradiance (L↓) after [51]. CPI is approximated from the bi-modality of L↓, i.e., clear or cloudy, and its partial impact on near surface temperature T. The dependence of L↓ on T is evaluated at each PROMICE site [52] and showed agreement with the blackbody theory, i.e., L↓ = εσTQ. For overcast conditions, ε was found to be 0.993 and Q was 4, approaching near blackbody radiative dependency. The parameter σ is the Stefan-Boltzmann constant (5.67 × 10⁻⁸ Wm⁻²K⁻⁴). For clear conditions the best fit to the fifth percentile resulted in variable values of ε and Q. Linear interpolation between the clear and cloudy estimates produces CPI values from zero to one for pairs of hourly L↓ and T measurements. Values below and exceeding the range of theoretical L↓ values were assumed to have been recorded respectively in clear and overcast conditions, i.e., CPI = 0.0 and CPI = 1.0.

Polar aerosols have typically low optical thickness (below 0.1 at the at 550 nm), have little influence on the top-of-atmosphere reflectance at the near-infrared wavelengths beyond ~865 nm [53], and thus, were currently neglected.

## 3.5. Validation of Snow OLCI Products

### 3.5.1. Snow Spectral Albedo

The spectral planar albedo products derived from OLCI scenes were compared to ground measurements of surface albedo. In a first step, the accuracy of the OLCI albedo product and its ability to capture the temporal evolution of the surface, and thus seasonal trends, was assessed using a dataset of sub-hourly observations of spectral surface albedo obtained at a fixed location at Dome C, Antarctica (75°06′00″S, 123°19′58″E, 3233 m elevation). Instrumentation details are provided below. In a second step, the evaluation of the OLCI spectral planar albedo retrievals at multiple locations spreading across the Antarctic plateau with varying snow surface properties was performed using portable spectral albedo measurements obtained along a 1371 km transect south of Dumont d'Urville (66°39′47″S, 140°00′10″E, 202 m). In a third step, the OLCI spectral planar albedo retrievals of snow containing impurities (atmospheric mineral aerosol deposits) were compared to portable spectral albedo measurements performed at Col du Lautaret, in the French Alps (45°02′07″N, 6°24′20″E, 2000 m), after a Saharan dust deposition event.

#### Data Processing and Instrumentation

An OLCI database was created for the different validation sites, applying the following steps. For each of the validation datasets, the OLCI L1c full-resolution products overlapping the location of the field measurements for the dates corresponding to the ground-based acquisitions were gathered. The Snow Properties Processor was run for all the remaining OLCI scenes, and the values for the pixels overlapping the field sites were extracted from the product. For all OLCI scenes, the SPP was run for the 21 OLCI bands, using band 21 (1020 nm) as the reference wavelength. In this study, the OLCI SVC gains (described in Section 3.3, see Figure 7) were not applied. For the scenes in Antarctica, the SPP was run without applying the snow pollution algorithm, whereas for the Col du Lautaret scenes, the snow pollution algorithm was used.

Spectral albedo measurements of the snow surface were acquired at Dome C using the Autosolex5 instrument [22,54,55], an automated spectral-radiometer that has been acquiring surface albedo measurement every 12 min between October and March since 2012. The instrument, shown in Figure 8, is described in detail in [22], and therefore only a succinct description is provided here.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Left: Autosolexs, the automated spectral albedometer at Dome C, Antarctica; right: measuring the surface albedo with the Solalb instrument during the ASUMA traverse in Antarctica.</figcaption></figure>

Autosolexs is composed of two measurement heads located approximately 1.5 m above the surface, each containing two upward and downward facing cosine receptors. The downward facing collectors were estimated to receive 75% of the reflected signal from the surface from an area of approximately 3.5 m in radius [22]. The receptors are connected to an Ocean Optics Maya2000 Pro spectrometer with fiber optics passing through an optical switch. Autosolexs measures the downward and upward (reflected) solar radiation. Using the calibration and processing steps described in [22], the spectral albedo was calculated as the ratio of the upwelling to downwelling spectral irradiances, providing usable data across the 400–1050 nm range. Owing to the negligible amount of diffuse illumination, the bi-hemispherical reflectance measured by the instrument (here called albedo) was directly compared to the directional–hemispherical reflectance (planar albedo) retrieved from OLCI images. For the validation purposes of this study, the datasets of the 2016–2017 and 2017–2018 summer seasons were used. The comparison between the ground observations and the satellite retrievals of spectral albedo was performed in several steps. Over the two summers, 1151 OLCI scenes were available at the Dome C location. For each scene, the closest Autosolexs albedo measurement in time was selected. To avoid a bias due to the albedo changing with the SZA and possibly changing snow physical properties, Sentinel-3 scenes more than 15 min apart from an Autosolexs measurement were automatically discarded. A second filter was applied to remove the measurements performed with a SZA larger than 75°, as it was shown in [22] that large SZA values introduce significant variations in the Autosolexs measurements. The remaining measured albedo spectra were smoothed with a 5 pixe-wide averaging sliding window to remove excessive noise due to instrument oversampling. For each measurement the spectra from the two heads of the instrument were averaged, allowing to diminish the effect of small-scale local surface roughness features [56], that may be located beneath the instrument, introducing a bias in the measurements owing to the small footprint of the instrument (approximately 40 m²). After the processing steps described above, 252 OLCI scenes were retained for comparison with the Autosolexs instrument over the two seasons. Although the full measured spectra are shown in Figure 9, in order to compare the ground values with each individual Sentinel 3 OLCI band the measured albedo was spectrally integrated over each bandwidth using the average spectral response provided by ESA.
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Intercomparison of the retrieved spectral planar albedo from Sentinel-3 OLCI using the Snow Properties Processor (SPP) (black squares) and Autosolexs albedo (blue lines) measured at Dome C, Antarctica for four different dates across the 2016–2017 and 2017–2018 summer seasons. The shaded area around the measurements represents the optimum uncertainty target of ±2%. The absorption bands of oxygen and water vapor are marked here as grey squares.</figcaption></figure>

For the spatial validation of the satellite albedo retrievals, the OLCI planar spectral albedo was compared to measurements obtained along a 1371 km transect on the Antarctic plateau south of the French Dumont d'Urville research station, over a period of 34 days between December 2016 and January 2017. The measurements were obtained within the framework of the ASUMA program (Improving the Accuracy of the SURface Mass balance of Antarctica). Albedo measurements were taken at 23 locations along the transect (Table 2) using the Solalb instrument (shown in Figure 8, a handheld version of Autosolexs). The instrument is composed of a 3 m boom at the end of which a light-collector is connected. It is characterized by a cosine response. The light-collector is connected to the same spectrometer model used in Autosolexs (Maya200pro, Ocean Optics) via an optic fiber. The instrument measures incoming radiation over the wavelength range 200–1300 nm with a usable range of 400–1050 nm. To measure albedo, the light-collector was first placed horizontally, in an upward-looking position to collect the incident irradiance, then in a downward-facing position to collect the upwelling reflected solar radiation from the snow surface. Following the processing chain developed for Autosolexs [22], the albedo was calculated as the ratio of reflected to incident solar radiation. For each of the 23 locations, multiple albedo measurements were performed at random positions within a radius of approximately 100 m, and the coordinates were recorded with a handheld GPS. The series of individual measurements at each location were averaged and compared to the corresponding OLCI retrieval for the overlapping pixel. The variance of the measurements was considered to represent the sub-pixel spatial variability of the pixel. For each set of field measurements, the closest cloud-free OLCI overpass in time was selected. If no OLCI cloud-free scenes were available on the day of the field measurements, the dataset was discarded. Because of the difference in time between the satellite overpass and the ground measurements, and thus the difference in SZA (Table 2), the measured albedo was converted to planar spectral albedo and corrected for the difference in SZA based on the asymptotic analytical radiative transfer theory presented in [19]. A full description of the
processing is available in Appendix E. After these processing steps, 11 OLCI scenes were retained for comparison with the Solalb measurements.

The OLCI albedo retrievals for polluted snow were validated using portable spectral albedo measurements obtained with the Solalb instrument previously described, at the Col du Lautaret field site with neighboring slopes of ~4°. Measurements were made on the 17, 19 and 24 April 2014 after Saharan mineral dust deposition events 14–16 and 22–23 April 2018. To minimize the effects of topography, the ground measurements acquired on 17 April were collected over a relatively flat (average slope of 3.2°) area of the plateau. Five spectral albedo measurements were obtained at random within a radius of approximately 10 m, with an average SZA of 55.3°. The averaged measurements were compared to the overlapping OLCI pixel (SZA = 41.2°). On 19 April, six surface albedo measurements were obtained every 40 m along a 240 m transect, to account for the spatial variability of the snow surface (slopes, snow physical properties, deposited dust). The average slope along the transect was 8.6° and the average SZA at the time of the field acquisitions was 40.9°. The entire transect was contained within an OLCI pixel, allowing a comparison of all the measured points with the satellite retrieval (SZA = 39.2°) obtained on the 20 April owing to the lack of Sentinel 3 data on 19 April. The comparison was considered valid, owing to the stable weather conditions during the two days. On 24 April, measurements were obtained at 24 locations along a 420 m transect on the plateau covering three OLCI pixels. However, two of the pixels contained only three and four measurements, respectively, located along the edge of the pixels, and thus, they were discarded, since the measurements were not considered representative of the pixels. Therefore, 17 measurements were retained for the comparison with the overlapping OLCI pixel. The slopes measured along the transect on the 24 April (average slope: 10.6°) were more pronounced than the slopes measured on the other dates. The SZA at the time of the field measurements was similar to the SZA at the time of the OLCI overpass with values of 36.5° and 36.9°, respectively.

Once processed, the albedo spectra were corrected for the difference in configuration between the measurement setup and the satellite retrievals. The ground spectral albedo acquisition was performed over sloping terrain, with a horizontal sensor, whereas the theory used in the retrievals assumes a horizontal surface. It is known that the surface albedo measurements are affected by the departure of the surface from horizontal [54], and therefore, the ground measurements need to be corrected for a direct comparison with the satellite spectral albedo. A full description of the algorithm appears in Appendix E. Furthermore, the difference in SZA between the ground acquisitions and the satellite overpass was accounted for, as described for the processing of the ASUMA dataset. Finally, for each date the albedo measurements were averaged, and compared to the closest spectral planar albedo retrieval in time for the overlapping OLCI pixel.

## Clean Snow

The spectral performance of the snow albedo retrieval algorithm is illustrated in Figure 9, where the OLCI snow spectral planar albedo retrievals at all 21 bands are compared to the concurrent. Autosolexys spectral albedo measurements for four dates selected across the 2016/17 and 2017/18 seasons. We have applied a scaling factor to the Autosolexys data to compensate measurement artifacts [22,54]. The scaling factor is obtained by fitting a two-parameter model to the observed spectrum in the 700–1050 nm range [22]. This leads to the ground measured plane albedo variability, which may be induced by surface roughness.

An excellent agreement across the wavelength range was observed, particularly for the two retrievals at smaller SZAs (<70°), where the satellite albedo values were within 2% of the ground measurements for all OLCI bands. For the retrievals with a SZA close to 70°, despite an overall good agreement, slightly more variability between observed and retrieved albedo was detected at larger wavelengths. This can be related to surface roughness effects, or problems with the atmospheric correction in case of large values of SZA.
The four dates illustrated in Figure 9 were selected to show the high radiometric accuracy of the retrievals, and are mostly representative of the two summer seasons at Dome C. It follows from Table 3 that the mean albedo bias over the two years was less than 1% for all bands below 1020 nm and it was smaller than 2% at 1020 nm. Considering the individual comparisons between the OLCI retrievals and the Autosolexs measurements at two wavelengths (865 and 1020 nm) in Figure 10, the values were well distributed along the identity line at 865 nm, and the satellite retrievals were mostly within 2% of the ground measurements for all dates, highlighting the ability of the retrieval algorithm to capture temporal variations. At 1020 nm, although the average bias is low, the retrievals tended to underestimate albedo in comparison to the ground measurements. Dependency of the satellite retrievals to the SZA is evident, with 90% of the points underestimating the albedo by more than 2% having SZA over 65°.

The validation of the OLCI spectral albedo retrievals for varied snow surface conditions across the Antarctic plateau is illustrated in Figure 11, where the average albedo measured with the Solalb instrument at each site and its standard deviation, representing that the subpixel variability is shown for four wavelength bands. Although the inter-site variability is low across the 1371 km transect, the wide error bars provide an indication of the strong sub-pixel spatial heterogeneity at each of the individual sites. Since the retrieved satellite albedo and the field measurements are not directly comparable owing to the differences (2.6 to 15.6°) of the SZA between the satellite overpasses and the measurements (Table 2), the measured albedo corrected for the SZA at the time of the overpass is also illustrated in Figure 11, with the assumption that the snow surface conditions did not change during the day. The data suggest an excellent radiometric accuracy between the spectral albedo and the retrieved satellite spectral albedo across the sites. Differences between the albedo estimated from the ground measurements and the satellite retrievals occur at 1020 nm, which we suspect is the result of the heterogeneity of the physical properties of the snow surface and could be explained by the under-representation of the pixel by the measured points (not shown here).

The comparison of OLCI spectral planar albedo with spectral surface albedo measurements at the Dome C site over 2 summer seasons shows that the OLCI snow products can be used to accurately monitor the albedo of clean snow surfaces and their temporal evolution over flat surfaces (see Table 3). However, the data have also shown that biases are introduced in the albedo retrieval algorithm for large SZAs. A part of uncertainty in the validation results remains due to differences introduced by the point-to-pixel comparison which does not account for the spatial heterogeneity (albeit low in the case of Dome C) and other scale effects.

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Scatterplot of the relationship between Sentinel-3 OLCI spectral planar albedo retrievals and the concurrent spectral albedo measured with the Autosolexs instrument at 2 OLCI bands (865, and 1020 nm). The thin dashed lines running on each side of the identity line represent $y = x \pm 0.01$, and the thin dotted lines in the right panel represent $y = x \pm 0.05$.</figcaption></figure>
<table><caption>**Table 2.** Location, Sentinel-3 acquisition dates and times, and SZA values for the sites along the ASUMA traverse used for the validation of the spectral albedo products.</caption><thead><tr><th rowspan="2">Site Name</th><th rowspan="2">Sentinel-3 OLCI Acquisition Date/Time (UTC)</th><th colspan="2">SZA/°</th><th rowspan="2">Average Latitude/Longitude of the Field Measurements</th></tr><tr><th>S3 Overpass</th><th>Measurement</th></tr></thead><tbody><tr><td>st1</td><td>2016/12/01, 22h21</td><td>59.14</td><td>64.1</td><td>-66.881150, 139.280006</td></tr><tr><td>st2</td><td>2016/12/02, 23h36</td><td>52.81</td><td>47.2</td><td>-67.167791, 138.968926</td></tr><tr><td>st6</td><td>2016/12/06, 21h51</td><td>63.25</td><td>68.5</td><td>-68.872819, 135.303749</td></tr><tr><td>st10</td><td>2016/12/10, 23h28</td><td>55.23</td><td>67.6</td><td>-69.636122, 135.279596</td></tr><tr><td>st16</td><td>2016/12/16, 22h32</td><td>58.60</td><td>70.8</td><td>-69.953844, 138.551435</td></tr><tr><td>st17</td><td>2016/12/17, 22h06</td><td>60.74</td><td>70.8</td><td>-69.953799, 138.553538</td></tr><tr><td>st21</td><td>2016/12/21, 22h02</td><td>60.01</td><td>75.6</td><td>-69.786595, 141.972969</td></tr><tr><td>st22</td><td>2016/12/22, 23h17</td><td>54.20</td><td>64.9</td><td>-69.786595, 141.972969</td></tr><tr><td>st25</td><td>2016/12/25, 23h40</td><td>53.39</td><td>64.7</td><td>-69.376932, 139.016187</td></tr><tr><td>st26</td><td>2016/12/26, 23h13</td><td>55.53</td><td>67.4</td><td>-69.060033, 138.225079</td></tr><tr><td>st27</td><td>2016/12/27, 22h47</td><td>57.90</td><td>55.3</td><td>-68.747273, 137.443734</td></tr></tbody></table>

<table><caption>**Table 3.** Comparison between the Sentinel-3 OLCI spectral planar albedo retrievals and the Autosolexks spectral albedo measurements over the 2016–2017 and 2017–2018 summers at Dome C, Antarctica. The standard deviation is calculated from the measurements or OLCI retrievals performed over the 2 years, representing the variability of the snow conditions. The mean albedo bias column shows that the error of retrievals is below 1% for all channels below 1020 nm and approximately 2% at 1020 nm.</caption><thead><tr><th rowspan="2">Sentinel-3 OLCI Band</th><th colspan="2">Sentinel-3 OLCI Spectral Planar Albedo Retrievals</th><th colspan="2">Autosolexks Spectral Albedo Measurements</th><th rowspan="2">Mean Albedo Bias (Absolute):</th><th rowspan="2">Mean Albedo Bias (%):</th></tr><tr><th>Average</th><th>Std</th><th>Average</th><th>Std</th></tr></thead><tbody><tr><td>01 (400 nm)</td><td>0.9988</td><td>0.0002</td><td>0.9899</td><td>0.0210</td><td>0.0089</td><td>0.9446</td></tr><tr><td>02 (421.5 nm)</td><td>0.9987</td><td>0.0002</td><td>0.9941</td><td>0.018</td><td>0.0046</td><td>0.4991</td></tr><tr><td>03 (442.5 nm)</td><td>0.9981</td><td>0.0003</td><td>0.9962</td><td>0.0136</td><td>0.0019</td><td>0.2099</td></tr><tr><td>04 (490 nm)</td><td>0.9955</td><td>0.0006</td><td>0.9971</td><td>0.0086</td><td>0.0016</td><td>0.1572</td></tr><tr><td>05 (510 nm)</td><td>0.9938</td><td>0.0008</td><td>0.9965</td><td>0.0061</td><td>0.0027</td><td>0.2651</td></tr><tr><td>06 (560 nm)</td><td>0.9889</td><td>0.0015</td><td>0.9929</td><td>0.0044</td><td>0.0039</td><td>0.3951</td></tr><tr><td>07 (620 nm)</td><td>0.9819</td><td>0.0024</td><td>0.9893</td><td>0.0040</td><td>0.0075</td><td>0.7550</td></tr><tr><td>08 (665 nm)</td><td>0.9749</td><td>0.0034</td><td>0.9812</td><td>0.0052</td><td>0.0062</td><td>0.6345</td></tr><tr><td>09 (673.75 nm)</td><td>0.9742</td><td>0.0034</td><td>0.9793</td><td>0.0052</td><td>0.0050</td><td>0.5133</td></tr><tr><td>10 (681.25 nm)</td><td>0.9731</td><td>0.0037</td><td>0.9776</td><td>0.0052</td><td>0.0045</td><td>0.4581</td></tr><tr><td>11 (708.75 nm)</td><td>0.9663</td><td>0.0046</td><td>0.9706</td><td>0.0057</td><td>0.0043</td><td>0.4466</td></tr><tr><td>12 (753.75 nm)</td><td>0.9575</td><td>0.0057</td><td>0.9574</td><td>0.0063</td><td>0.0001</td><td>0.0112</td></tr><tr><td>13 (761.25 nm)</td><td>0.9536</td><td>0.0062</td><td>0.9584</td><td>0.0072</td><td>0.0047</td><td>0.4905</td></tr><tr><td>14 (764.375 nm)</td><td>0.9524</td><td>0.0064</td><td>0.9544</td><td>0.0072</td><td>0.0020</td><td>0.2072</td></tr><tr><td>15 (767.5 nm)</td><td>0.9508</td><td>0.0066</td><td>0.9510</td><td>0.0073</td><td>0.0003</td><td>0.0270</td></tr><tr><td>16 (778.75 nm)</td><td>0.9452</td><td>0.0073</td><td>0.9455</td><td>0.0078</td><td>0.0002</td><td>0.0248</td></tr><tr><td>17 (865 nm)</td><td>0.9213</td><td>0.0104</td><td>0.9226</td><td>0.0099</td><td>0.0013</td><td>0.1444</td></tr><tr><td>18 (885 nm)</td><td>0.9055</td><td>0.0124</td><td>0.9056</td><td>0.0118</td><td>0.0001</td><td>0.0050</td></tr><tr><td>19 (900 nm)</td><td>0.8992</td><td>0.0132</td><td>0.8975</td><td>0.0127</td><td>0.0017</td><td>0.1906</td></tr><tr><td>20 (940 nm)</td><td>0.8875</td><td>0.0146</td><td>0.8908</td><td>0.0147</td><td>0.0032</td><td>0.3576</td></tr><tr><td>21 (1020 nm)</td><td>0.7939</td><td>0.0254</td><td>0.8075</td><td>0.0258</td><td>0.0136</td><td>1.6661</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Comparison of the retrieved Sentinel-3 OLCI spectral planar albedo and the measured albedo for the sites along the ASUMA transect in Antarctica. The error bars of the spatialised measurements at each site are an indicator of the sub-pixel variability. The orange markers represent the planar albedo calculated using the field measurements for the SZA at the time of the satellite acquisition.</figcaption></figure>
## Polluted Snow

The OLCI retrievals of spectral planar albedo using the snow pollution branch of the retrieval algorithm for three dates at the Col du Lautaret validation site (see Figure 12) are illustrated in Figure 13. The OLCI retrievals were compared to ground measurements from a single location on 17 April 2018 and along transects on the 19 and 24 April 2018. The calculated planar albedo accounting for the slope under the sensor and the SZA at the overpass time, and the spatial variability of the measurements are also illustrated in the figure. For the three dates, the figure highlights the ability of the algorithm to retrieve the spectral signature of the snow containing dust, i.e., with enhanced absorption in the 400 nm to 650 nm range [25]. The measurements obtained at a single location on 17 April 2018 and along a transect across a homogenous slope on 19 April 2018 exhibited little sub-pixel variability, indicating similar snow physical properties and dust concentrations across the site (see Figure 12). On the other hand, the transect performed on 24 April 2018 across a larger range of slopes, suggests a greater variability of the snow surface, probably owing to the differences in snow properties between the north and south facing slopes. When taking into account the slope and difference in SZA between the measurement and the satellite retrieval, an excellent agreement was observed on 17 April 2018. On 19 April 2018, a good agreement was also observed, with an overestimation of the OLCI retrieval between 400 and 510 nm and between 940 and 1020 nm. Nevertheless, the satellite retrieval is able to estimate the overall spectral shape of the polluted snow surface. On 24 April, the OLCI retrieval overestimated the albedo at all wavelengths, with differences between ground and OLCI albedo of 8% at 400 nm and 16% at 1020 nm. The disagreement between the OLCI retrieval and the measurements on the last date may be explained by the larger slopes across the pixel (10.6°), the higher sub-pixel variability in the measurements, and the location of transect on the plateau, closer to surrounding topographic features than the sites measured on the two previous days. Indeed, the measurements made on 17 and 19 April 2018 were obtained in flatter terrain (average slopes of 3.2° and 8.6°, respectively) and further away from large surrounding topographic features (cliffs, peaks), which have been shown to significantly impact satellite measurements due to the interactions between the solar radiation and the topographic features [54]. Therefore, in the configuration of 24 April 2018 the effect of the slope can not be accounted for using the algorithm described in Appendix E.

<figure><img src="image_16.png" /><figcaption>**Figure 12.** The Solalb light collector shown over a polluted snowpack at Col du Lautaret, French Alps.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Comparison of the retrieved Sentinel-3 OLCI spectral planar albedo for polluted snow (black markers) and the measured albedo (black lines) at Col du Lautaret for three dates in April 2018, after a dust deposition event. The planar albedo accounting for the slope under the sensor calculated for the SZA at the time of the satellite overpass based on the measurements is shown in blue. The blue shaded area represents the sub-pixel variability.</figcaption></figure>

### 3.5.2. Broadband Albedo

#### Data Processing and Instrumentation

In Greenland, snow and ice albedo was monitored by the Programme for Monitoring of the Greenland Ice Sheet automatic weather stations (AWSs) located at 24 sites on the Greenland ice sheet (Figure 14). The PROMICE AWS record shortwave irradiance upward and downward (300 nm to 2400 nm) each 10 min using a Kipp CM3 pyranometer as part of a Kipp CNR-1 or CNR-4 radiometer bundle (see Figure 13). The CM3 thermopile is covered by a single glass dome and complies with ISO 9060 s-class specifications, i.e., ±10% accuracy for daily totals. Under conditions with no frost or water on domes and a level station, as is common here, the CM3 performs better than ±5%. PROMICE AWS tilt measurements (see Figure 14) were used to exclude cases having tilt greater than ±1.5 degrees. Furthermore, to avoid larger measurement uncertainty, albedo values include only those instantaneous
measurements for which sunlight incidence angle for both the ice sheet surface and the upper dome of
the sensor was larger than 20 degrees.

<figure><img src="image_18.png" /></figure>

1: solar and infrared radiation, 2: tilt sensor, 3: satellite antenna, 4: wind speed & direction, 5: snow/ice surface height, 6: air temperature & humidity, 7: ice ablation 'hose', 8: solar panel, 9: data logger, barometer and GPS, 10: battery, 11: ice temperature profile (8 levels)

<caption>**Figure 14.** PROMICE_map, (left) location of Programme for Monitoring of the Greenland Ice Sheet (PROMICE) automatic weather stations (AWSs). (Right) photo of a PROMICE AWS with instrumentation indicated with numbering.</caption>

In an effort to neutralize PROMICE station shortwave radiation obstruction effects by the AWS
measurement platform, we apply the correction used in [56] (equation 16). An assumption was that the
OLCI BBA<sub>shortwave</sub> for dry clean snow had no bias, as shown for clean snow in the Antarctica spectral
albedo validation. In the correction, we assumed a constant station (mast, enclosure, solar panel)
average reflectivity of 0.4, and a fraction of station obstruction coming from above and below the
pyranometer. We used a variable [56] measurement frame obstruction "f" parameter to produce zero
average bias with the OLCI retrieval for dry snow, i.e., when PROMICE hourly average albedo is
above 0.875. Experimenting with correction factors leads to the need for variable f parameter on a
site by site basis, corresponding with differences in the relative spacing and orientation of the AWS
solar panel, data logger enclosure and whether snow persists on the surface (covering the battery box).
The correction increases the average PROMICE albedo by 0.027-0.081. We did not compare the results
for all PROMICE stations because for the Sentinel-3 years, the excluded sites did not have useful data.
These are the best-behaved results.

The nearest OLCI 300 m pixel values were extracted for Greenland PROMICE network locations
during five-month periods spanning 1 May to 31 September in years 2017 and 2018. For some sites and
days, multiple images are available. Data from the closest hour are compared. The NASA Moderate
Resolution Imaging Spectroradiometer (MODIS) MOD10A1 product [7] collection 6 [12] provides daily
broadband albedo (BBA<sub>shortwave</sub>) retrievals. The BBA<sub>shortwave</sub> retrieved using OLCI measurements is
compared with those derived from both MODIS and PROMICE measurements. The results are given
in the next section.
Broadband Albedo Validation

Comparison of shortwave (300 nm to 2400 nm) BBA<sub>OLCI</sub> and BBA<sub>MOD10A1</sub> with BBA<sub>PROMICE</sub>
demonstrate that both the OLCI-derived albedo and the independent PROMICE ground observations
record a common climate signal, capturing a wide range of albedo values (0.25 to 0.85), despite their
footprint size differing by a factor of ~104 (Figure 15). A realistic BBA<sub>OLCI</sub> is evident over a wide
range of albedo. Pearson's correlation coefficient for BBA<sub>OLCI</sub> versus BBA<sub>PROMICE</sub> average 0.870 when
comparing the cases with PROMICE radiometer tilt under 1.5 degrees and cloud probability index
under 30% (Table 4). For these cases, the comparison had an average BBA root mean squared difference
(RMSD) of 0.050 and an insignificant average difference, i.e., a bias of 0.011. Regression slopes average
less than unity, due to a wider range of ground observations, which may well be the case given the
large difference in spatial scale between the 300 m OLCI pixel and the ~6 m² footprint size of the
BBA<sub>PROMICE</sub> data. BBA<sub>PROMICE</sub> captures a more patchy surface while the OLCI (and MODIS) retrievals
average out patch and sastrugi-scale surface variations. The regression slope departure from unity is
statistically insignificant with some PROMICE stations having a large (~25%) slope departure from
unity. BBA<sub>OLCI</sub> correlations with BBA PROMICE are usually higher than BBA<sub>MOD10A1</sub> correlations.
Additionally, the RMS differences are lower for OLCI retrievals of BBA<sub>shortwave</sub> (see Figure 15).

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Comparison of OLCI and MODIS derived broadband albedo for a selection of Greenland ice sheet PROMICE ground station locations.</figcaption></figure>
<table><caption>**Table 4.** Broadband albedo comparison statistics between OLCI and Greenland ice sheet PROMICE stations.</caption><thead><tr><th>Site Name</th><th>Latitude, Degrees N</th><th>Longitude, Degrees W</th><th>Elevation, m</th><th>Regression Slope</th><th>Regression Constant</th><th>Correlation</th><th>Average Difference, OLCI-PROMICE</th><th>RMSD</th><th>N obs</th></tr></thead><tbody><tr><td>KPC_L</td><td>79.908</td><td>24.080</td><td>366</td><td>0.734</td><td>0.217</td><td>0.936</td><td>0.053</td><td>0.082</td><td>221</td></tr><tr><td>SCO_U</td><td>72.394</td><td>27.259</td><td>988</td><td>0.890</td><td>0.094</td><td>0.921</td><td>0.024</td><td>0.046</td><td>22</td></tr><tr><td>QAS_L</td><td>61.031</td><td>46.849</td><td>270</td><td>0.900</td><td>0.046</td><td>0.979</td><td>-0.008</td><td>0.055</td><td>23</td></tr><tr><td>KAN_M</td><td>67.067</td><td>48.838</td><td>1268</td><td>0.861</td><td>0.081</td><td>0.868</td><td>-0.017</td><td>0.072</td><td>66</td></tr><tr><td>KAN_U</td><td>67.000</td><td>47.027</td><td>1842</td><td>0.571</td><td>0.347</td><td>0.707</td><td>-0.003</td><td>0.022</td><td>100</td></tr><tr><td>UPE_U</td><td>72.887</td><td>53.585</td><td>929</td><td>1.106</td><td>-0.033</td><td>0.918</td><td>0.030</td><td>0.066</td><td>62</td></tr><tr><td>EGP</td><td>75.625</td><td>35.973</td><td>2660</td><td>0.842</td><td>0.134</td><td>0.760</td><td>0.001</td><td>0.007</td><td>74</td></tr><tr><td>mean</td><td>70.8</td><td>-40.5</td><td>1189</td><td>0.843</td><td>0.127</td><td>0.870</td><td>0.011</td><td>0.050</td><td>81</td></tr><tr><td>std.</td><td>6.3</td><td>11.5</td><td>839</td><td>0.164</td><td>0.124</td><td>0.100</td><td>0.025</td><td>0.027</td><td>67</td></tr></tbody></table>

### 3.5.3. Specific Surface Area

#### Measurements in Antarctica

The Specific Surface Area (SSA) retrievals based on OLCI images were evaluated using the Autosolexs time-series at Dome C over the 2016–2017 and 2017–2018 summer seasons. The temporal comparison allows assessing the ability for the OLCI retrieval to capture daily to seasonal variations of the snow surface physical properties. The albedo dataset provided by the Autosolexs instrument was used to derive the SSA of the snow surface at Dome C. The SSA was retrieved by fitting a theoretical model to the observed albedo, based on the analytical asymptotic radiative transfer theory [19,22].

The seasonal and daily variations in SSA retrieved from the ground data were well reproduced by the satellite retrievals, with an underestimation of the ground SSA obtained from OLCI (Figure 16). The systematic bias present in both years, and much more pronounced for the 2017–2018 summer, was consistent with the under-estimation of the albedo in the near-infrared illustrated in Figure 9. The range of SSA values observed with both instruments, nevertheless, was in agreement with the range observed by [55] between 2012 and 2014 at Dome C. Although it was more pronounced during summer 2016–2017, Figure 16 illustrates the decrease in SSA over both observed seasons due to the increased metamorphism caused by the summertime enhanced solar radiation and temperatures. Furthermore, the satellite SSA retrievals are able to capture the annual cycle of SSA, and its daily variations, which are due to snowfall or clear sky ice crystal precipitation, i.e., diamond dust events [55]. However, a solar zenith angle dependency was observed in the OLCI retrievals (as was illustrated in Figure 9 for spectral albedo). Indeed, unrealistic and significant variations in the SSA were observed for retrievals obtained on the same day at different times (and thus SZA). Although such sub-daily variations were observed in the Autosolexs SSA dataset, the variations are much smaller for the ground-based retrievals, which led us to think that the OLCI retrievals are much more sensitive to the SZA than the Autosolexs measurements.

The comparison of the OLCI SSA retrievals with the SSA retrieved from the Autosolexs instrument indicates that the algorithm is suitable to track the temporal evolution of variables describing the physical properties of the snowpack in flat polar regions, which is critical to parameterize and validate snow models. Indeed, exploiting the high revisit frequency over the Poles, the algorithm was able to resolve the short term (weekly, and sometimes daily) variations of the snow conditions, as well as long term trends.

It should be pointed out that the OLCI snow specific surface area was determined from the derived value of the snow grain diameter. Therefore, measurements at the wavelength 1020 nm were used. The ground measurements of SSA were based on spectral albedo measurements. Multiple wavelengths in the framework of optimal estimation approach were used (assuming $g = 0.845$). Because the retrievals of plane albedo from OLCI and ground were very close, the difference in observed SSA can be attributed to different algorithms of the SSA determination using reflectance measurements. It was of importance to validate OLCI SSA retrievals based on non-optical methods of SSA determination. This will be the subject of our future work. It should be pointed out that the optimal estimation approach based
on ground spectral planar albedo measurements has been tested in the range 10–40 m²kg⁻¹ using an independent technique [57]. An excellent agreement with the measurements based on independent methane adsorption technique was found.

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Time series of SSA calculated using the Autosolexs spectral albedo measurements, and the Sentinel-3 OLCI spectral planar albedo retrievals over Dome C for the 2016/17 and 2017/18 summer seasons. The gaps in the time series are due to periods of time with cloud cover.</figcaption></figure>

## Measurements in Greenland

Daily snow specific surface area (SSA) was measured in two consecutive field seasons at the EastGRIP Greenland site (EGP in Figure 14) using the IceCube instrument, in accordance with [37]. IceCube illuminates a 5.0 cm diameter × 2.5 cm depth cylindrical surface snow sample underneath an integrating sphere with a 1310 nm laser diode. Snow samples were taken every 10 m along a 90 m long transect, producing 10 daily samples. The observed SSA values range from below 20 m²kg⁻¹ to above 70 m²kg⁻¹ after snowfall. The standard deviation is typically in the range of 10 m²kg⁻¹, but the distribution was heavily skewed to high values (skewness = 1.4) occurring briefly due to snowfall. The average observed SSA along the transect is compared with the nearest 300 m OLCI SSA retrieval. The comparison covers the middle of the sunlit season; early May to mid-August in 2017 and 2018. The comparison is illustrated in Figure 17 (see also Table 5). Agreement was better below ~35 m²kg⁻¹. When SSA is extremely high, it is a condition of higher scattering.

The difference in footprint size between the ground measurement and the OLCI retrieval is a source of uncertainty. The ground observers note that some samples are taken from areas of snow deposition (dunes and lee slopes of sastrugi) and others from areas of snow erosion (sastrugi). This character of spatial variability is found immediately after a fresh snow event due to wind redistribution of the snow surface at the meter-scale. Spatial variability then decreases with time. The OLCI footprint is on a scale that is above the meter-scale variations captured by the 90 m transect of 1 m samples.

In the other EastGRIP comparison, three observations were made within 2 min of satellite overpass times and during the tandem Sentinel-3A and Sentinel-3B phase of the Sentinel-3B commissioning. This work was part of a JAXA Global Change Observation Mission-Climate (GCOM-C)/Second generation Global Imager (SGLI) field campaign by T. Aoki, S. Matoba, M. Niwano, and R. Shimada. In this campaign, SSA was also observed using the IceCube instrument and the data were converted to
optically equivalent snow grain diameter using Equation (32). The OLCI-derived snow grain diameter retrievals are extremely close to these coincident observations (Table 4). The absolute difference (12.3 µm for S3A and 2.5 µm for S3B) is much smaller than the ~20 µm uncertainty in the field data. This finding must be confirmed for a larger dataset, which is currently under investigation.

<figure><img src="image_21.png" /><figcaption>**Figure 17.** Comparison of OLCI retrievals of snow Specific Surface Area (SSA) from the average of 10 ground observations each day along a 90 m transect at the EGP location on Greenland.</figcaption></figure>

<table><caption>**Table 5.** Comparison of OLCI surface (0–1 cm depth) snow grain size retrieval in comparison with in-situ observations.</caption><thead><tr><th>Satellite</th><th>Date</th><th>Overpass Time, UTC</th><th>Observation Time, UTC</th><th>OLCI-Retrieved Diameter of Grains, mm</th><th>In-Situ Derived Diameter of Grains, mm</th><th>Difference in Diameters (OLCI-In-Situ), µm</th><th>% Bias</th></tr></thead><tbody><tr><td rowspan="3">S3A</td><td>8 July 2018</td><td>1414</td><td>1416</td><td>0.138</td><td>0.145</td><td>-6.8</td><td>0.95%</td></tr><tr><td>9 July 2018</td><td>1529</td><td>1528</td><td>0.209</td><td>0.206</td><td>3.5</td><td>1.02%</td></tr><tr><td>13 July 2018</td><td>1525</td><td>1523</td><td>0.164</td><td>0.124</td><td>40.1</td><td>1.32%</td></tr><tr><td rowspan="3">S3B</td><td>8 July 2018</td><td>1414</td><td>1416</td><td>0.135</td><td>0.145</td><td>-10.1</td><td>0.93%</td></tr><tr><td>9 July 2018</td><td>1528</td><td>1528</td><td>0.176</td><td>0.206</td><td>-30.1</td><td>0.85%</td></tr><tr><td>13 July 2018</td><td>1525</td><td>1523</td><td>0.173</td><td>0.124</td><td>47.6</td><td>1.38%</td></tr><tr><td></td><td></td><td></td><td>average</td><td>0.161</td><td>0.158</td><td>2.5</td><td>1.02%</td></tr></tbody></table>

## 4. Discussion

The results given above confirm that the developed technique is robust and can be used to derive various snow properties from OLCI observations. Example of OLCI retrieval maps (see Figures 18 and 19) for the land ice areas of Greenland during July illustrate a wide range of values. We chose July for the examples because snow metamorphic and ice ablation effects are strong at this time of year from high solar insolation, high temperatures and snow/ice ablation effects. In image mosaics, retrieval value difference along image boundaries (a.k.a. “fringes”) are the result of some combination of real physical properties changes and the effect of differing illumination and viewing geometry.

Due to low temperature and therefore low snow grain metamorphism in the high elevation (above 2500 m and up to 3250 m) dry snow area, small (~0.15 mm) grain diameters (and high SSA above
10 m²kg⁻¹, not shown) are typical (Figure 18a). As elevation decreases, retrieved grain diameters increase above ~1.8 mm, where snow metamorphic effects are strongest. Inversely proportional to grain diameter, the retrieved SSA reaches minimum values of ~4 m²kg⁻¹, not shown. A highlighted area of higher grain diameter values along the northeastern ice sheet is where the effect of localized precipitation else heating were observed after 25 July 2017, which persisted through 31 July 2017.

Broadband albedo on this date (Figure 18a) is uniformly high (above 0.78) in the dry snow area where snow metamorphism is low. The impact on broadband albedo of the snow grain size increase in the circled high elevation area of the northeastern ice sheet is subtle. However, in longer wavelengths, the grain size impact is substantial. This covariation of snow albedo and grain size indicates the importance of retrieving, not only the surface albedo, but also the characteristics that are driving changes of snow reflectivity. The lowest elevation where snow persists, commonly referred to as snowline, has retrieved albedo of 0.6, consistent with observations. Below snowline, albedo drops sharply. The bare ice has a variable concentration of impurities, the most absorptive of which are microbiological in origin (e.g., [58]). Bare ice albedo retrievals have minimum values below 0.35 in agreement with field measurements (e.g., [59]).

Daily retrievals or monthly averages are calculated from pixels not identified as cloudy (Figure 19). Monthly sample counts indicate that as few as 1/3 of possible temporal coverage occurs in areas, where surface based ice fog occurrence is high from an environment dominated by negative net radiation. Absolute variance in the averages is at a maximum around the ice sheet periphery where surface snow conditions are more temporally variable due to transient snowfall, snow drift and melting events thereby affecting the surface properties. These maps of surface snow characteristics can be used to validate the albedo calculated by climate models (e.g., [60]) or derived from other remote sensing platforms or can even be assimilated by regional climate models to prescribe more reliable snow characteristics to their snow module (e.g., [61]).

<figure><img src="image_22.png" /><figcaption>**Figure 18.** Example Sentinel-3 daily mosaics of (a) snow grain diameter for areas classified as snow using broadband albedo above 0.55 and (b) OLCI-derived broadband albedo. The chosen date on 28 July 2017 was a date with little cloud cover over the entire ice sheet.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 19.** Example monthly statistics from OLCI-derived snow physical properties; (a) grain diameter; (b) specific surface area; (c) planar albedo at 1020 nm; (d) broadband albedo (300–2500 nm).</figcaption></figure>

**5. Conclusions**

This work presents first retrievals of snow optical and microphysical properties using the OLCI instrument onboard the Sentinel-3A and B satellites. The retrieval of spectral snow albedo for the case of clean snow is performed in two steps: (1) the absorption length and non-absorbing snow
reflectance are derived using OLCI measurements at 865 and 1020 nm, least affected by atmospheric effects; (2) the fact that spherical albedo of clean snow (in the range 400–1020 nm covered by OLCI) is uniquely defined by the absorption length is used. This enables also the determination of plane albedo accounting for the solar zenith angle and also broadband albedo extrapolating the retrieved spherical albedo to other wavelengths in the range 300–2400 nm. The retrieval for the polluted snow case is more involved and requires a priori knowledge of non-absorbing snow reflectance, which is derived from exact radiative transfer calculations for the case of clouds composed of ice crystals.

The work also emphasizes validation. The underlying theory, data processor plugin for the open source SNAP and validation of various snow products document an advanced snow data delivery tool planned for data exploitation for snow covered areas globally. The algorithm is not based on the look-up table approach and can be used with small modifications for other optical instruments orbiting our planet (MODIS, S-GLI, etc.). The algorithm may be biased at specific observation conditions (forward/backward scattering, high contamination of snow by pollutants, vertically inhomogeneous snow and sloping terrain) as discussed by [54,62–64].

Nevertheless, the validation steps presented here demonstrate that the algorithm developed to derive spectral albedo of clean dry snow over flat expanses from OLCI images is suitable to estimate surface spectral albedo with 2–3% accuracy, for solar zenith angles under 70 degrees, which is sufficient for many climate studies. For an alpine site with high aerosol mineral dust loading, ground validation suggests reliability in retrievals of the spectral albedo. However, with rugged terrain, field measurements cannot be directly compared to the retrievals without accounting for the terrain effects.

Shortwave broadband (300 nm to 2400 nm) albedo from OLCI and MODIS MOD10A1 product compared with ground observations from PROMICE.org automatic weather stations demonstrate that both the OLCI-derived albedo and the independent ground observations record a common climate signal, accurately capturing a very wide range (0.25 to 0.85) of albedo, despite their footprint size differing by a factor of ~104. Thus, we can conclude a realistic broadband albedo from our OLCI retrievals for dry snow, wet snow and impurity rich bare glacier ice.

For broadband shortwave albedo, OLCI-retrieval correlations are typically higher than those from MODIS (the MOD10A1 product) versus the PROMICE AWS albedo measurements and the root mean squared (RMS) differences are lower. Thus, it appears that OLCI-retrieved broadband albedo is possibly of higher accuracy than the MODIS-products, paving the way for a multi-decade reliable climate data record of an essential climate variable. Where surface slope and terrain shading occur, additional compensation of radiometric effects come into focus to achieve accurate results for areas off ice sheets, where surface slopes often exceed a few degrees.

**Author Contributions:** Conceptualization, A.K., J.E.B., C.B., G.P. and M.L.; methodology, M.L., M.D., J.E.B., A.K., G.P.; software, M.L., O.D. and A.K.; validation, L.A., V.F., B.J., E.L., B.D.M., T.A., M.N., S.K. (Sepp Kipfstuhl), J.F., M.H., A.Z., D.V., A.-K.F., H.C.S.-L., S.W., J.K.A., B.V., D.v.A. and K.M.; formal analysis, A.K., V.R., E.Z., S.K. (Sergey Korkin), J.E.B., M.L., G.P., M.D., M.K. and C.B.; investigation, A.K., M.L., J.E.B. and O.D.; resources, O.D., C.B. and J.E.B.; data curation, M.L. and J.E.B.; writing—Original draft preparation, A.K., M.L., J.E.B., H.C.S.-L. and C.B.; writing—Review and editing, A.K., J.E.B., M.L., H.C.S.-L. and B.V.; visualization, K.M., J.B. and B.V., supervision, M.K., J.B., C.B., G.P. and M.D.; project administration, M.K., J.E.B., C.B., G.P. and M.D.; funding acquisition, J.B., C.B. and A.K.

**Funding:** This research has been funded by ESRIN contract 4000118926/16/I-NB and the ESRIN contract 4000125043 - ESA/AO/1-9101/17/I-NB EO SCIENCE FOR SOCIETY.

**Acknowledgments:** Agency (ESA) studies; the Scientific Exploitation of Operational Missions (SEOM) Sentinel-3 Snow (Sentinel-3 for Science, Land Study 1: Snow), ESRIN contract 4000118926/16/I-NB and the ESRIN contract 4000125043 - ESA/AO/1-9101/17/I-NB EO SCIENCE FOR SOCIETY. Additional support came from The Program for the Monitoring of the Greenland Ice Sheet (PROMICE), part of the Danish Energy Agency through the DANCEA program. The work of Sergey Korkin is supported by grant "NASA Science for Terra, Aqua and SNPP" (PI: Alexei Lyapustin, 17-TASNPP17-0116; solicitation NNH17ZDA001NTASNPP). EGRIP is directed and organized by the Center of Ice and Climate at the Niels Bohr Institute. It is supported by funding agencies and institutions in Denmark (A. P. Møller Foundation, UCPH), USA (US National Science Foundation, Office of Polar Programs), Germany (AWI), Japan (NIPR and ArCS), Norway (UiB and BFS), Switzerland (SNF), France (IPEV, IGE) and China
(CAS and BNU). This work was conducted as a part of the GCOM-C/SGLI Project supported by the Japan Aerospace Exploration Agency (JAXA). The authors acknowledge the support from Agence Nationale de la Recherche Scientifique for the scientific traverse in Antarctica and the associated research: project ANR-14-CE01-0001 (ASUMA). The authors are grateful to S. Delwart and S. Warren for important discussions related to peculiarities of snow remote sensing from space and snow optical properties.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. The Derivation of Main Equations

Let us consider the case of absorbing snow. The reflectance can be presented as

$$R(\beta) = \sum_{n=1}^{\infty} a_n (1 - \beta)^n. \qquad (A1)$$

where $\beta = 1 - \omega_0$ is the probability of photon absorption by unit volume of snow. In the case of nonabsorbing snow it follows

$$R(0) = \sum_{n=1}^{\infty} a_n. \qquad (A2)$$

Let us define the ratio

$$\mathfrak{R} = \frac{R(\beta)}{R(0)}. \qquad (A3)$$

Then, we have

$$\mathfrak{R} = \frac{\sum_{n=1}^{\infty} a_n (1 - \beta)^n}{\sum_{n=1}^{\infty} a_n}. \qquad (A4)$$

Expanding the nominator of Equation (A4), we derive at small $\beta$

$$\mathfrak{R} \approx 1 - \beta \langle n \rangle + \frac{\beta^2}{2} \langle n^2 \rangle - \frac{\beta^3}{6} \langle n^3 \rangle + \dots \approx \langle \exp(-\beta n) \rangle. \qquad (A5)$$

where

$$\langle n^k \rangle \equiv \sum_{n=1}^{\infty} f_n n^k, \quad \langle \exp(-\beta n) \rangle \equiv \sum_{n=1}^{\infty} f_n \exp(-\beta n). \qquad (A6)$$

$$f_n = \frac{a_n}{\sum_{n=1}^{\infty} a_n}. \qquad (A7)$$

and we assumed that $n(n-1) \approx n^2, n(n-1)(n-2) \approx n^3 \dots$ in our derivations. This is possible because $\beta$ is close to zero and the number of scattering events is high. For the same reason, we have

$$\langle \exp(-\beta n) \rangle \approx \int_{0}^{\infty} f(n) \exp(-\beta n) dn. \qquad (A8)$$

This integral can be evaluated assuming the function $f(n)$. In particular, it follows from the random walk theory that the probability of a particle (photon) appearing at a given place, time and direction after large number of iterations $n$ can be presented as [65]

$$f(n) = \sqrt{\frac{\eta}{\pi}} n^{-3/2} \exp\left\{-\frac{\eta}{n}\right\}. \qquad (A9)$$

where the parameter $\eta$ depends on the process studied. The substitution of Equation (A9) into Equation (A8) gives

$$\mathfrak{R} = \exp(-2\sqrt{\eta\beta}). \qquad (A10)$$

Therefore, we can write:
$$R(\beta) = R_0 \exp(-\sqrt{s\beta}). \qquad (A11)$$

where $R_0 = R(0)$, $s = 4\eta$. This equation has been proposed for the first time in [30] (see also [31]). It shows how the spectral snow reflectance depends on the probability of photon absorption $\beta$. The parameter $s$ depends on the scattering and not on absorption processes and, therefore, one may assume that it does not depend on the wavelength for snow composed of large snow grains in contact. It should be pointed out that we have derived Equation (A9) not using the radiative transfer equation (RTE). In fact, the applicability of RTE for closely-packed media such as snow is in question. It should be pointed out that Equation (A11) is very general and can be applied to many geophysical media such as snow, leaves, etc. The main shortcoming of this equation as applied to natural snow is the fact that it can fail in areas close to forward scattering (glint) and backward scattering, where snow has enhanced brightness, e.g., due to oriented snow crystal surfaces or thin ice films on the top of snowpack. Such directions must be avoided in the snow property retrieval process. Additionally, this equation can be applied only for the horizontal snow surfaces. Otherwise, a modification is needed.

The value of $s$ can be related to the asymmetry parameter $g$ of ice grains using asymptotic results of RTE valid at small values of the probability of photon absorption. Then, it follows from Equation (A9)

$$R(\beta) = R_0(1 - \sqrt{s\beta}). \qquad (A12)$$

moreover, from analytical radiative transfer theory [33]

$$R(\beta) = R_0 - yu(\mu_0)u(\mu). \qquad (A13)$$

where

$$y = 4 \sqrt{\frac{\beta}{3(1-g)}}. \qquad (A14)$$

and

$$u_0(\mu_0) = \frac{3}{4}(\mu_0 + \varphi(\mu_0)). \qquad (A15)$$

$$\varphi(\mu_0) = 2 \int_0^1 \langle R_0(\mu_0, \mu) \rangle \mu^2 d\mu. \qquad (A16)$$

$$\langle R_0(\mu_0, \mu) \rangle = \frac{1}{2\pi} \int_0^1 R_0(\mu, \mu_0, \psi) d\psi. \qquad (A17)$$

One can show [33] that the following approximation holds

$$\varphi(\mu_0) = \frac{4}{7} \left[ 1 + \frac{\mu_0}{4} \right]. \qquad (A18)$$

and, therefore

$$u(\mu_0) = \frac{3}{7} [1 + 2\mu_0]. \qquad (A19)$$

Comparing Equations (A12) and (A13), one derives

$$s = \frac{16f^2}{3(1-g)}. \qquad (A20)$$

where

$$f = \frac{u(\mu)u(\mu_0)}{R_0}. \qquad (A21)$$
## Appendix B. The Retrieval Errors and Uncertainties

The absolute error of the retrieved parameter $x_j$ can be presented in the following way

$$ \Delta x_j = \sqrt{\sum_{n=1}^{N} \left(\frac{\partial x_j}{\partial y_n}\right)^2 (\Delta y_n)^2}. \qquad (A22) $$

where $N$ is the number of measurements and $y_n$ are the results of the measurements (e.g., snow reflectivity at $N$ channels as measured from space) performed with absolute errors $\Delta y_n$. Equation (A22) can be also presented in the following way for the relative errors

$$ \frac{\Delta x_j}{x_j} = \sqrt{\sum_{n=1}^{N} \Xi_{jn}^2 \left[ \frac{\Delta y_n}{y_n} \right]^2}. \qquad (A23) $$

where

$$ \Xi_{jn} \equiv \frac{y_n}{x_j} \left( \frac{\partial x_j}{\partial y_n} \right). \qquad (A24) $$

These coefficients can be found analytically for the pair ($R_0, l$) (see Equations (26) and (27)). Namely, it follows

$$ \frac{\Delta R_0}{R_0} = \sqrt{\varepsilon^2 \left(\frac{\Delta R_1}{R_1}\right)^2 + (1-\varepsilon)^2 \left(\frac{\Delta R_2}{R_2}\right)^2}. \qquad (A25) $$

$$ \frac{\Delta l}{l} = \sqrt{\varepsilon^2 \left(\frac{\Delta R_1}{R_1}\right)^2 + (2-\varepsilon)^2 \left(\frac{\Delta R_2}{R_2}\right)^2}. \qquad (A26) $$

where

$$ \varepsilon = (1-b)^{-1}, \quad b = \sqrt{\frac{\alpha_1}{\alpha_2}}, \quad \alpha_j = \frac{4\pi\chi_j}{\lambda_j}, \quad \sigma = 2(\varepsilon-z), \quad z = \ln^{-1}\left(\frac{R_2}{R_1}\right). \qquad (A27) $$

where $\chi_j$ is the bulk pure ice refractive index measured the j-th wavelength and $R_j$ is the measured reflectance at the j-th wavelength. The accuracy of the determination of the pair is determined by the accuracy of reflectance measurements at two wavelengths in near-infrared. In particular, the wavelengths 865 and 1020 nm can be used. The error enhancement coefficients depend on the ratio of the bulk ice absorption coefficients at the two selected wavelengths. We conclude that to have accurate results one should select the pair of the wavelengths, which are not too close one to another.

In the case of the retrievals based on the wavelength pair 865 and 1020 nm, one derives that $\varepsilon = 1.55$. Let us assume that the relative errors in both channels are the same and equal to $\Delta$. Then, it follows

$$ \frac{\Delta R_0}{R_0} = k_1 \Delta, \quad \frac{\Delta l}{l} = k_2 \Delta. \qquad (A28) $$

where

$$ k_1 = \sqrt{1 + 2\varepsilon(\varepsilon - 1)}, \quad k_2 = 2\sqrt{1 + 2(\varepsilon - z)(\varepsilon - 1 - z)}. \qquad (A29) $$

The error enhancement factors given by Equations (A29) are presented in Figures A1 and A2. We find that the first enhancement factor is equal to 1.6 at $\varepsilon = 1.55$. This means that the error of the determination of the parameter reflectance of nonabsorbing snow is better than 4% in case the reflectances are determined with the accuracy 2%. The dependence of the second factor on the ratio of reflectances $x = R_2/R_1$ is given in Figure A2. It follows that the accuracy increases for smaller values of $x$ (e.g., for large ice particles). For a typical value of $x = 0.6$, it follows: $k_2 = 8$. This means that the error of the retrieval of the parameter absorption length is about 16% in this particular case.
<figure><img src="image_24.png" /><figcaption>**Figure A1.** The dependence of the error enhancement coefficient on the parameter $ε$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure A2.** The dependence of the enhancement factor on the value of the ratio of reflectance at two channels.</figcaption></figure>
## Appendix C. The Reflection Function of a Semi-Infinite Non-Absorbing Snow

It is assumed that the reflection function of a semi-infinite snow layer at single scattering albedo equal to one can be approximated by the following expression, proposed in [66]:

$$R_0 = \frac{A + B(\mu_0 + \mu_0) + C\mu_0\mu + p}{4(\mu_0 + \mu)} \quad (A30)$$

where $A = 1.247$, $B = 1.186$, $C = 5.157$,

$$p = 11.1 \exp(-0.087\theta) + 1.1 \exp(-0.014\theta). \quad (A31)$$

and $\theta$ is the scattering angle in degrees. It holds:

$$\cos\theta = -\mu_0\mu + ss_0\cos\varphi. \quad (A32)$$

Here, $s$ and $s_0$ are the sines of the viewing zenith angle (VZA) and SZA, respectively, $\varphi$ is the relative azimuthal angle and $(\mu_0, \mu)$ are the cosines of the SZA and VZA, respectively. The OLCI relative azimuth angle $\upsilon$ must be transformed as: $\varphi = \text{abs}(180 - \upsilon)$ to be used in the equations given above.

## Appendix D. The Error Budget in the Case of Polluted Snow

The absolute error of the retrieved parameter $x_j$ can be presented in the following way

$$\Delta x_j = \sqrt{\sum_{n=1}^{N} \left( \frac{\partial x_j}{\partial y_n} \right)^2 (\Delta y_n)^2}. \quad (A33)$$

where $N$ is the number of measurements and $y_n$ are the results of the measurements (e.g., snow reflectivity at $N$ channels as measured from space) performed with absolute errors $\Delta y_n$. The symbol $j$ signifies the parameter under study (see Table A1 for actual snow parameters retrieved from snow observations). Equation (A33) can be also presented in the following way for the relative errors

$$\frac{\Delta x_j}{x_j} = \sqrt{\sum_{n=1}^{N} \Xi_{jn}^2 \left[ \frac{\Delta y_n}{y_n} \right]^2}. \quad (A34)$$

Let us find the matrix elements

$$\Xi_{jn} \equiv \frac{y_n}{x_j} \left( \frac{\partial x_j}{\partial y_n} \right). \quad (A35)$$

for the snow parameters mentioned in Table A1.

<table><caption>**Table A1.** The snow parameters retrieved from satellite observations ($y_n = R_n$).</caption><thead><tr><th>N</th><th>Parameter</th><th>Equation</th><th>Comments</th><th>Physical Meaning</th></tr></thead><tbody><tr><td>1</td><td>x<sub>1</sub> ≡ R<sub>0</sub></td><td>R<sub>3</sub><sup>ε</sup>R<sub>4</sub><sup>1-ε</sup></td><td>ε = (1 - b)<sup>-1</sup>,<br/>b = √(α<sub>3</sub>/α<sub>4</sub>)</td><td>reflectance of nonabsorbing snow layer</td></tr><tr><td>2</td><td>x<sub>2</sub> ≡ l</td><td>(R<sub>0</sub><sup>2</sup>/α<sub>4</sub>u<sub>0</sub><sup>2</sup>(μ<sub>0</sub>)u<sub>0</sub><sup>2</sup>(μ)) ln<sup>2</sup>[R<sub>4</sub>/R<sub>0</sub>]</td><td></td><td>effective absorption length</td></tr><tr><td>3</td><td>x<sub>3</sub> ≡ m</td><td>(ln[p<sub>1</sub>/p<sub>2</sub>]) / (ln[λ<sub>2</sub>/λ<sub>1</sub>])</td><td>p<sub>k</sub> = ln<sup>2</sup>(R<sub>k</sub>/R<sub>0</sub>)</td><td>absorption Angström parameter</td></tr><tr><td>4</td><td>x<sub>4</sub> ≡ κ</td><td>(R<sub>0</sub><sup>2</sup>λ<sub>1</sub><sup>m</sup>/u<sub>0</sub><sup>2</sup>(μ<sub>0</sub>)u<sub>0</sub><sup>2</sup>(μ)) ln<sup>2</sup>[R<sub>1</sub>/R<sub>0</sub>]</td><td></td><td>the snow impurity absorption coefficient at the wavelength λ<sub>0</sub> (normalized to the snow grain concentration and also to the value of the absorption enhancement factor B)</td></tr></tbody></table>
We can also present the relative errors in the following way:

$$ 
\begin{aligned}
\frac{\Delta R_0}{R_0} &= \sqrt{\varepsilon_1^2(\frac{\Delta R_3}{R_3})^2 + \varepsilon_2^2(\frac{\Delta R_4}{R_4})^2}, \quad \frac{\Delta l}{l} = \sqrt{v_1^2(\frac{\Delta R_3}{R_3})^2 + v_2^2(\frac{\Delta R_4}{R_4})^2} \\
\frac{\Delta m}{m} &= \sqrt{\sum_{n=1}^4 w_n^2(\frac{\Delta R_n}{R_n})^2}, \quad \frac{\Delta \tilde{\kappa}}{\tilde{\kappa}} = \sqrt{\sum_{n=1}^4 h_n^2(\frac{\Delta R_n}{R_n})^2}.
\end{aligned}
\quad (A36) $$

The corresponding derivatives can be found analytically. The answer for the coefficients is

$$ 
\begin{aligned}
\varepsilon_1 &= (1 - b)^{-1}, \\
\varepsilon_2 &= 1 - \varepsilon_1, \\
v_1 &= -2\varepsilon_1 z_4, \\
v_2 &= -v_1, \\
w_1 &= \frac{2z_1}{m\ln(\lambda_2/\lambda_1)}, \\
w_2 &= -\frac{2z_2}{m\ln(\lambda_2/\lambda_1)}, \\
w_3 &= -(w_1 + w_2)\varepsilon_1, \\
w_4 &= -(w_1 + w_2)\varepsilon_2, \\
h_1 &= 2(1 + s)z_1, \\
h_2 &= -2sz_2, \\
h_3 &= -\varepsilon_1(h_1 + h_2) - v_1, \\
h_4 &= -\varepsilon_2(h_1 + h_2) - v_2,
\end{aligned}
\quad (A37) $$

where

$$ s = \frac{\ln \tilde{\lambda}_1}{\ln(\lambda_2/\lambda_1)}, z_i = \ln^{-1}(R_i/R_0), p_i = \ln^2(R_i/R_0), b = \sqrt{\frac{\alpha_3}{\alpha_4}}. \quad (A38) $$

These coefficients determine the matrix coefficients introduced in Equation (A35). The expressions given above enable to make the following conclusions:

* The parameters ($R_0, l$) can be found from the measurements at two wavelengths in the near-infrared. Therefore, the accuracy of their determination is determined by the accuracy of reflectance measurements at two wavelengths in near-infrared. In particular, the wavelengths 865 and 1020 nm can be used. The error enhancement coefficients depend on the ratio of the bulk ice absorption coefficients at the corresponding wavelengths. Both coefficients increase, if the wavelengths selected are to close one to another (see Appendix C).

* The accuracy of the pair ($\tilde{\kappa}, m$) determination depends on the accuracy of measurements at all four channels (two-in the visible, two- in the near infrared) used in the retrieval process.

The derived equations can be applied to plan the experimental measurements and assess the relative errors of the retrieved parameters for a given set of the wavelengths.

## Appendix E. Converting Measured Spectral Albedo Over Sloping Terrain to Spectral Planar Albedo in a Flat Terrain Configuration

The hand-held field measurements, performed with a horizontal sensor (which is not parallel to the surface for sloping surfaces), are not in the same configuration as the satellite spectral planar albedo retrievals that consider a flat surface (plane-parallel radiative transfer theory). In order to compare the ground measurements with the satellite retrievals, the following algorithm converting the field measurements into plane-parallel spherical or planar albedo was applied. As shown by [54], the measured albedo with a horizontal sensor over a surface (which may be sloping) can be approximated as follows for small slopes

$$ r_{\text{meas}}(\lambda) = D(\lambda, \theta_s)r_s(\lambda) + c(\theta_s, \tilde{\theta}_s)(1 - D(\lambda, \theta_s))r_p(\lambda, \tilde{\theta}_s). \quad (A39) $$
where $r_s(\lambda)$ is the spherical albedo (albedo under diffuse illumination, i.e., white sky albedo) and $r_p(\lambda, \tilde{\theta}_s)$ is the plane albedo (albedo under direct illumination, i.e., black sky albedo), $D(\lambda, \theta_s)$ is the diffuse-to-total irradiance ratio for a wavelength $\lambda$, and

$$ \zeta = \frac{\cos \tilde{\theta}_s}{\cos \theta_s}. \qquad (A40) $$

where the effective solar zenith angle is $\tilde{\theta}_s$ [54], which represents the effective direction of the sun with respect to the tilted surface, and $\theta_s$ is the solar zenith angle. In this study, the diffuse-to-total irradiance ratio, $D(\lambda, \theta_s)$, was calculated using the SBDART detailed radiative transfer code [43]. For the calculations performed on the ASUMA and Col du Lautaret datasets, the solar geometries at the time of the measurements, and the altitude at which the measurements were performed were used as inputs to the model. Assuming a clear sky, a standard Sub-Arctic summer atmosphere model with no boundary layer aerosols, and an optical depth of 0.05 at 550 nm for the background stratospheric aerosols, was selected for the ASUMA calculations. For the Col du Lautaret measurements, the diffuse-to-total irradiance calculations were performed using a standard mid-latitude winter atmosphere, with a rural boundary layer aerosol type, and an optical depth of 0.15 at 550 nm. The atmosphere and aerosol models are fully described in the SBDART documentation.

The spectral albedo measured in the field is converted to spectral spherical albedo and solved iteratively

$$ r_s^{i+1}(\lambda) = \frac{r_{\text{meas}}(\lambda) + \zeta(\theta_s, \tilde{\theta}_s)(1 - D(\lambda, \theta_s))(r_s^i(\lambda) - (r_s^i(\lambda))^n)}{\zeta(\theta_s, \tilde{\theta}_s)(1 - D(\lambda, \theta_s)) + D(\lambda, \theta_s)}. \qquad (A41) $$

where the $r_s^i(\lambda)$ is initialised with $r_{s,\text{meas}}(\lambda)$ in the first iteration. Following Equation (16) in this paper, the spectral planar albedo can be expressed as

$$ r_p(\lambda, \tilde{\theta}_s) = r_s^u(\lambda). \qquad (A42) $$

where

$$ u = \frac{3}{7}(1 + 2 \cos \tilde{\theta}_s). \qquad (A43) $$

To compare spectral planar albedo measurements and retrievals acquired at different solar zenith angles over flat terrain (or corrected for the slope following the previous equations), the albedo can be estimated for a given solar zenith angle based on the theory described in Section 2.1 of this paper (Equation 16)

$$ r_p(\lambda, \theta_2) = r_p(\lambda, \theta_1)^{\frac{(1+2\cos(\theta_2))}{(1+2\cos(\theta_1))}}. \qquad (A44) $$

with $\theta_1$ is the solar zenith angle at the time of the input measurement or retrieval of spectral planar albedo and $\theta_2$ is the solar zenith angle for which the estimation is made.

## References

1. Hansen, J.; Nazarenko, L. Soot climate forcing via snow and ice albedos. Proc. Natl. Acad. Sci. USA **2004**, *101*, 423–428. [CrossRef] [PubMed]

2. Warren, S.G. Optical properties of ice and snow. Philos. Trans. R. Soc. A **2019**, *377*, 20180161. [CrossRef] [PubMed]

3. Frei, A.; Tedesco, M.; Lee, S.; Foster, J.; Hall, D.K.; Kelly, R.; Robinson, D.J. A review of global satellite—derived snow products. Adv. Space Res. **2012**, *50*, 1007–1029. [CrossRef]

4. Tedesco, M. (Ed.) Remote Sensing of the Cryosphere; Wiley: New York, NY, USA, 2015.

5. Sirguey, P.; Mathieu, R.; Arnaud, Y. Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: Methodology and accuracy assessment. Remote Sens. Environ. **2009**, *113*, 160–181. [CrossRef]
6. Klein, A.; Stroeve, J. Development and a validation of a snow albedo algorithm for the MODIS instrument. *Ann. Glaciol.* **2001**, *34*, 45–52. [CrossRef]

7. Hall, D.K.; Riggs, G.A.; Salomonson, V.V. Development of methods for mapping global snow cover using moderate resolution imaging spectroradiometer data. *Remote Sens. Environ.* **1995**, *54*, 127–140. [CrossRef]

8. Hall, D.K.; Riggs, G.A.; Solomonson, V.V.; De Girolamo, N.E.; Bayr, K.J. MODIS snow cover products. *Remote Sens. Environ.* **2002**, *83*, 181–194. [CrossRef]

9. Lyasputin, A.; Tedesco, M.; Wang, Y.; Aoki, T.; Hori, M.; Kokhanovsky, A.A. Retrieval of snow grain size over Greenland from MODIS. *Remote Sens. Environ.* **2009**, *113*, 1976–1987. [CrossRef]

10. Painter, T.H.; Rittger, K.; McKenzie, S.; Slaughter, P.; Davis, R.E.; Dozier, J. Retrieval of subpixel snow covered area, grain size, and albedo from MODIS. *Remote Sens. Environ.* **2009**, *113*, 868–879. [CrossRef]

11. Zege, E.P.; Katsev, I.L.; Prikhach, A.; Heygster, G.; Wiebe, H. Algorithm for retrieval of the effective snow grain size and pollution amount from satellite measurements. *Remote Sens. Environ.* **2011**, *115*, 2674–2685. [CrossRef]

12. Riggs, G.A.; Hall, D.K. MODIS Snow Products Collection 6, User Guide, 2015. Available online: https://nsidc.org/sites/nsidc.org/files/files/MODIS-snow-user-guide-C6.pdf (accessed on 23 September 2019).

13. Stamnes, K.; Li, W.; Eide, H.; Aoki, T.; Hori, M.; Storvold, R. ADEOS-II/GLI Snow/ice Products—Part I: Scientific Basis. *Remote Sens. Environ.* **2007**, *111*, 258–273. [CrossRef]

14. Li, W.; Stamnes, K.; Chen, B.; Xiong, X. Snow grain size retrieved from near-infrared radiances at multiple wavelengths. *Geophys. Res. Lett.* **2011**, *28*, 1699–1702. [CrossRef]

15. Kokhanovsky, A.A.; Rozanov, V.V.; Aoki, T.; Odermatt, D.; Brockmann, C.; Krüger, O.; Bouvet, M.; Drusch, M.; Hori, M. Sizing snow grains using backscattered solar light. *Int. J. Remote Sens.* **2011**, *32*, 6975–7008. [CrossRef]

16. Appel, I. Retrieval and Validation of VIIRS Snow Cover Information for Terrestrial Water Cycle Applications. In *Remote Sensing of Terrestrial Water Cycle*; Laksmi, V., Alsdorf, D., Anderson, M., Biancamaria, S., Cosh, M., Entin, J., Huffman, G., Kustas, W., van Oevelen, P., Painter, T., et al., Eds.; AGU: Washington, DC, USA, 2014; pp. 175–197. [CrossRef]

17. Gascoin, S.; Grizzonet, M.; Bouchet, M.; Salgues, G.; Hagolle, O. Theia snow collection: High resolution operational snow cover maps from Sentinel-2 and Lansat-8 data. *Earth Syst. Sci. Data* **2019**, *11*, 493–514. [CrossRef]

18. Donlon, C.; Berruti, B.; Buongiorno, A.; Ferreira, M.-H.; Féménias, P.; Frerick, J.; Goryl, P.; Klein, U.; Laur, H.; Mavrocordatos, C.; et al. The Global Monitoring for Environment and Security (GMES) Sentinel-3 mission. *Remote Sens. Environ.* **2012**, *120*, 37–57. [CrossRef]

19. Kokhanovsky, A.A.; Zege, E.P. Scattering optics of snow. *Appl. Opt.* **2004**, *43*, 1589–1602. [CrossRef]

20. Kokhanovsky, A.; Lamare, M.; Di Mauro, B.; Picard, G.; Arnaud, L.; Dumont, M.; Tuzet, F.; Brockmann, C.; Box, J.E. On the reflectance spectroscopy of snow. *Cryosphere* **2018**, *12*, 2371–2382. [CrossRef]

21. Warren, S.G.; Brandt, R.E. Optical constants of ice from the ultraviolet to the microwave: A revised compilation. *J. Geophys. Res.* **2008**, *113*, 3681–3691. [CrossRef]

22. Picard, G.; Libois, Q.; Arnaud, L. Refinement of the ice absorption spectrum in the visible using radiance profile measurements in Antarctic snow. *Cryosphere* **2016**, *10*, 2655–2672. [CrossRef]

23. Grenfell, T.C.; Warren, S.G.; Mullen, P.C. Reflection of solar radiation by the Antarctic snow surface at ultraviolet, visible, and near-infrared wavelengths. *J. Geophys. Res.* **1994**, *99*, 18669–18684. [CrossRef]

24. Dumont, M.; Brissaud, O.; Picard, G.; Schmitt, B.; Gallet, J.-C.; Arnaud, Y. High-accuracy measurements of snow Bidirectional Reflectance Distribution Function at visible and NIR wavelengths—Comparison with modelling results. *Atmos. Chem. Phys.* **2010**, *10*, 2507–2520. [CrossRef]

25. Di Mauro, B.; Fava, F.; Ferrero, L.; Garzonio, R.; Baccolo, G.; Delmonte, B.; Colombo, R. Mineral dust impact on snow radiative properties in the European Alps combining ground, UAV, and satellite observations. *J. Geophys. Res. Atmos.* **2015**, *120*, 6080–6097. [CrossRef]

26. He, C.; Liou, K.-N.; Takano, Y.; Yang, P.; Qi, L.; Chen, F. Impact of grain shape and multiple black carbon internal mixing on snow albedo: Parameterization and radiative effect analysis. *J. Geophys. Res.* **2018**, *123*, 1253–1268. [CrossRef]

27. Kokhanovsky, A.A.; Aoki, T.; Hachikubo, A.; Hori, M.; Zege, E.P. Reflective properties of natural snow: Approximate asymptotic theory versus in situ measurements. *IEEE Trans. Geosci. Rem. Sens.* **2005**, *43*, 1529–1535. [CrossRef]
28. Ding, A.; Jiao, Z.; Dong, Y.; Qu, Y.; Zhang, X.; Xiong, C.; He, D.; Yin, S.; Cui, L.; Chang, Y. An assessment of the performance of two snow kernels in characterizing snow scattering properties. *Int. J. Remote Sens.* **2019**, *40*, 6315–6335. [CrossRef]

29. Jiao, Z.; Ding, A.; Kokhanovsky, A.; Schaaf, C.; Bréon, F.; Dong, Y.; Wang, Z.; Liu, Y.; Zhang, X.; Yin, S.; et al. Development of a snow kernel to better model the anisotropic reflectance of pure snow in a kernel-driven BRDF model framework. *Remote Sens. Environ.* **2019**, *221*, 198–209. [CrossRef]

30. Rosenberg, G.V. Optical characteristics of thick weakly absorbing scattering layers. *Dokl. Akad. Nauk* **1962**, *6*, 775–777.

31. Zege, E.P.; Katsev, I.L.; Ivanov, A.P. *Image Transfer Through Scattering Media*; Springer: Berlin/Heidelberg, Germany, 1991.

32. Kokhanovsky, A. Statistical properties of photon gas in random media. *Phys. Rev.* **2002**, *E66*, 037601. [CrossRef]

33. Sobolev, V.V. *Light scattering in planetary atmospheres: International series of monographs in natural philosophy*; M.: Nauka; Elsevier: Pergamon, Turkey, 1975.

34. Korkin, S.; Lyapustin, A. Matrix exponential in C/C++ version of vector radiative transfer code IPOL. *J. Quant. Spectrosc. Radiat. Transf.* **2019**, *227*, 106–110. [CrossRef]

35. Emde, C.; Barlakas, V.; Cornet, C.; Evans, F.; Korkin, S.; Ota, Y.; Labonnote, L.C.; Lyapustin, A.; Macke, A.; Mayer, B.; et al. IPRT polarized radiative transfer model intercomparison project—Phase A. *J. Quant. Spectrosc. Radiat. Transf.* **2015**, *164*, 8–36. [CrossRef]

36. Picard, G.; Arnaud, L.; Dominé, F.; Fily, M. Determining snow specific surface area from near-infrared reflectance measurements: Numerical study of the influence of grain shape. *Cold Reg. Sci. Technol.* **2009**, *56*, 10–17. [CrossRef]

37. Gallet, J.-C.; Domine, F.; Zender, C.; Picard, G. Measurement of the specific surface area of snow using infrared reflectance in an integrating sphere at 1310 and 1550 nm. *Crysosphere* **2009**, *3*, 167–182. [CrossRef]

38. Van de Hulst, H.C. *Light Scattering by Small Particles*; Dover: New York, NY, USA, 1981.

39. Kokhanovsky, A.A. Scaling constant and its determination from simultaneous measurements of light reflection and methane adsorption by snow samples. *Opt. Lett.* **2006**, *31*, 3282–3284. [CrossRef] [PubMed]

40. Libois, Q.; Picard, G.; Dumont, M.; Arnaud, L.; Sergent, C.; Pougatch, E.; Sudul, M.; Vial, D. Experimental determination of the absorption enhancement parameter of snow. *J. Glaciol.* **2014**, *7*, 714–724. [CrossRef]

41. Kokhanovsky, A.A.; Zege, E.P. Local optical parameters of spherical polydispersions: Simple approximations. *Appl. Opt.* **1995**, *34*, 5513–5519. [CrossRef] [PubMed]

42. Dombrovsky, L.A.; Kokhanovsky, A.A. The influence of pollution on solar heating and melting of a snowpack. *J. Quant. Spectrosc. Radiat. Transf.* **2019**, *233*, 42–51. [CrossRef]

43. Ricchiazzi, P.; Yang, S.; Gautier, C.; Sowle, D. SBDART, A research and teaching tool for plane-parellel radiative transfer in the Earth’s atmosphere. *Bull. Am. Meteorol. Soc.* **1998**, *79*, 2101–2114. [CrossRef]

44. Wiscombe, W.J.; Warren, S.G. A Model for the Spectral Albedo of Snow. I: Pure Snow. *J. Atmos. Sci.* **1980**, *37*, 2712–2733. [CrossRef]

45. Domine, F.; Salvatori, R.; Legagneux, L.; Salzano, R.; Fily, M.; Cassacchia, R. Correlation between the specific surface area and the short wave infrared (SWIR) reflectance of snow. *Cold Reg. Sci. Technol.* **2006**, *46*, 60–68. [CrossRef]

46. Räisänen, P.; Kokhanovsky, A.; Guyot, G.; Jourdan, O.; Nousiainen, T. Parameterization of single-scattering properties of snow. *Cryosphere* **2015**, *9*, 1277–1301. [CrossRef]

47. Hori, M.; Aoki, T.; Stamnes, K.; Chen, B.; Li, C. Preliminary validation of the GLI cryosphere algorithms with MODIS daytime data. *Polar Meteorol. Glaciol.* **2001**, *15*, 1–20.

48. Garrett, T.J. Observational quantification of the optical properties of cirrus cloud. In *Light Scattering Reviews*; Kokhanovsky, A., Ed.; Praxis-Springer: Berlin/Heidelberg, Germany, 2008; Volume 3, pp. 1–26.

49. Linow, S.; Horhold, M.W.; Freitag, J. Grain-size evolution of polar firn: A new empirical grain growth parameterization based on X-ray microcomputer tomography measurements. *J. Glaciol.* **2012**, *58*, 1245–1252. [CrossRef]

50. Dozier, J.; Frew, J. Atmospheric corrections to satellite radiometric data over rugged terrain. *Remote Sens. Environ.* **1981**, *11*, 191–205. [CrossRef]

51. Van As, D. Warming, glacier melt and surface energy budget from weather station observations in the Melville Bay region of northwest Greenland. *J. Glaciol.* **2011**, *57*, 208–220. [CrossRef]
52. Van As, D.; van den Broeke, M.R.; Reijmer, C.; van de Wal, R. The summer surface energy balance of the high Antarctic plateau. *Bound. Layer Meteorol.* **2005**, *115*, 289–317.

53. Tomasi, C.; Kokhanovsky, A.; Lupi, A.; Ritter, C.; Smirnov, A.; O'Neill, N.T.; Stone, R.S.; Holben, B.N.; Nyeki, S.; Wehrli, C.; et al. Aerosol remote sensing in polar regions. *Earth Sci. Rev.* **2015**, *140*, 108–157. [CrossRef]

54. Dumont, M.; Arnaud, L.; Picard, G.; Libois, Q.; Lejeune, Y.; Nabat, P.; Voisin, D.; Morin, S. In situ continuous visible and near-infrared spectroscopy of an alpine snowpack. *Cryosphere* **2017**, *11*, 1091–1110. [CrossRef]

55. Libois, Q.; Picard, G.; Arnaud, L.; Dumont, M.; Lafaysse, M.; Morin, S.; Lefebvre, E. Summertime evolution of snow specific surface area close to the surface on the Antarctic Plateau. *Cryosphere* **2015**, *9*, 2383–2398. [CrossRef]

56. Arnaud, L.; Picard, G.; Champollion, N.; Domine, F.; Gallet, J.C.; Lefebvre, E.; MFily; Barnola, J.M. Measurement of vertical profiles of snow specific surface area with a 1cm resolution using infrared reflectance: Instrument description and validation. *J. Glaciol.* **2011**, *57*, 17–29. [CrossRef]

57. Aoki, T.; Kuchiki, K.; Niwano, M.; Kodama, Y.; Hosaka, M.; Tanaka, T. Physically based snow albedo model for calculating broadband albedos and the solar heating profile in snowpack for general circulation models. *J. Geophys. Res. Atmos.* **2011**, *116*. [CrossRef]

58. Stibal, M.; Box, J.E.; Cameron, K.A.; Langen, P.L.; Yallop, M.L.; Mottram, R.H.; Khan, A.L.; Molotch, N.P.; Chrismas, N.A.M.; Quaglia, F.C.; et al. Algae drive enhanced darkening of bare ice on the Greenland ice sheet. *Geophys. Res. Lett.* **2017**, *44*, 11463–11471. [CrossRef]

59. Wientjes, I.G.M.; Van de Wal, R.S.W.; Reichart, G.J.; Sluijs, A.; Oerlemans, J. Dust from the dark region in the western ablation zone of the Greenland ice sheet. *Cryosphere* **2011**, *5*, 589–601. [CrossRef]

60. Van Angelen, J.H.; Lenaerts, J.T.M.; Lhermitte, S.; Fettweis, X.; Kuipers Munneke, P.; van den Broeke, M.R.; van Meijgaard, E.; Smeets, C.J.P.P. Sensitivity of Greenland Ice Sheet surface mass balance to surface albedo parameterization: A study with a regional climate model. *Cryosphere* **2012**, *6*, 1175–1186. [CrossRef]

61. Langen, P.L.; Fausto, R.S.; Vandegrux, B.; Mottram, R.H.; Box, J.E. Liquid water flow and retention on the Greenland ice sheet in the regional climate model HIRHAM5: Local and large-scale impacts. *Front. Earth Sci.* **2017**, *4*, 110. [CrossRef]

62. Shao, D.; Xu, W.; Li, H.; Wang, J.; Hao, X. Reconstruction of remotely sensed snow albedo for quality improvements based on a combination of forward and retrieval models. *IEEE Trans. Geosci. Remote Sens.* **2018**, *56*, 6969–6985. [CrossRef]

63. Lv, Y.; Wu, D.; Sub, Z. Effect of black carbon concentration on the reflection property of snow: A comparison with model results. *IEEE Trans. Geosci. Remote Sens.* **2018**, *56*, 6823–6839. [CrossRef]

64. Saito, M.; Yang, P.; Loeb, N.G.; Kato, S. A novel parameterization of snow albedo based on a two-layer snow model with a mixture of habits. *J. Atmos. Sci.* **2019**, *76*, 1419–1436. [CrossRef]

65. Chandrasekhar, S. Stochastic problems in physics and astronomy. *Rev. Mod. Phys.* **1943**, *1*, 1–89. [CrossRef]

66. Kokhanovsky, A.A. Reflection of light from particulate media with irregularly shaped particles. *J. Quant. Spectrosc. Radiat. Transf.* **2005**, *96*, 1–10. [CrossRef]

<figure><img src="image_26.png" /></figure>