An Improved Retrieval of Snow and Ice Properties Using Spaceborne OLCI/S-3 Spectral Reflectance Measurements: Updated Atmospheric Correction and Snow Impurity Load Estimation

Abstract

We present an update of the Snow and Ice (SICE) property retrieval algorithm based on the spectral measurements of Ocean and Land Color Instrument (OLCI) onboard Sentinel-3 satellites combined with the asymptotic radiative transfer theory valid for weakly absorbing turbid media. The main improvements include the introduction of a new atmospheric correction, retrieval of snow impurity load and properties, retrievals for partially snow-covered ground and also accounting for various thresholds to be used to assess the retrieval quality. The technique can be applied to various optical sensors (satellite and ground-based) operated in the visible and near infrared regions of electromagnetic spectra.

1. Introduction

This paper presents a fast and accurate approach for the determination of snow properties from the measurements of the Ocean and Land Color Instrument (OLCI) onboard Sentinel-3 satellites (see also http://snow.geus.dk/, accessed on 10 November 2022). The algorithm targets retrieval of snow and ice optical properties, such as snow spectral and broadband albedo, important for climatic studies [1] and snow microstructure (snow specific surface area and effective optical grain size (see [2]). We propose a technique to retrieve the concentration of snow impurities (dust, soot, see [3]), which is possible only for the cases of relatively heavy (above 1 ppmv) impurity loading [4] as the concentration of snow impurities is often low, below 100 ng/g, especially in polar regions [5]. The detailed description of the previous versions of the retrieval algorithm is given in [6,7,8] and an application for Greenland ice sheet broadband albedo and bare ice area appears in [9]. The main updates presented in this work are:

-

The retrievals of snow fraction and sub-pixel snow albedo;

-

The improvement of the atmospheric correction technique [10];

-

Characterization of snow impurities in terms of their type (dust, soot) and load;

-

Retrieval of total ozone column as discussed in [11];

-

Introduction of retrieval quality metrics based on the root-square-mean-difference (RSMD) between measured and retrieved OLCI spectra;

-

Updated value of the calibration coefficient for the relation between the effective absorption length and the effective snow grain diameter;

-

Accounting for the recent results related to OLCI gains [12].

In this work, we describe an updated version of SICE algorithms, validate the retrieval results with ground in situ observations and demonstrate the capabilities of the technique using retrievals over Greenland and Antarctica. The algorithm is made available as either a Python or Fortran package at https://github.com/GEUS-SICE/pySICE, accessed on 10 November 2022, and archived in [13].

2. Retrieval of Snow Properties Using OLCI Observations

2.1. Ocean and Land and Colour Instrument

The Ocean and Land and Colour Instrument is a 21 band imaging spectroradiometer that measures solar radiation reflected by the Earth’s atmosphere and surface with a ground spatial resolution up to 300 m (see

Table 1. Band characteristics of the Sentinel-3 Ocean and Land Colour Instrument (OLCI).

Band	λ Centre (nm)	Width (nm)	Band	λ Centre (nm)	Width (nm)	Band	λ Centre (nm)	Width (nm)
1	400	15	8	665	10	15	767.5	2.5
2	412.5	10	9	673.75	7.5	16	778.75	15
3	442.5	10	10	681.25	7.5	17	865	20
4	490	10	11	708.75	10	18	885	10
5	510	10	12	753.75	7.5	19	900	10
6	560	10	13	761.25	2.5	20	940	20
7	620	10	14	764.375	3.75	21	1020	40

https://sentinel.esa.int/web/sentinel/user-guides/sentinel-3-olci/resolutions/radiometric (accessed on 1 December 2022). Channels affected by oxygen and water vapor absorption are given in red color.

). The OLCI swath width is 1270 km. OLCI is installed on both Sentinel-3A and Sentinel-3B satellites operated by the European Space Agency (ESA) in service to the European Union Copernicus Program. The Sentinel-3 A and B orbit at 802 km altitude, 98.6 orbital inclination and a 10:00 UTC sun-synchronous equatorial crossing time. A Sentinel-3 C platform is planned for launch in 2024. The mission is useful for snow monitoring because spectral signatures of clean and polluted snow in the spectral range of OLCI and Sea and Land Surface Temperature Radiometer (SLSTR), another instrument installed on Sentinel-3, vary depending on the level of snow pollution and snow properties, such as effective optical snow grain size.

2.2. Definitions

2.2.1. Geometry of the System

The angles describing the solar illumination and observation direction (Figure 1) from which we obtain the cosine of the solar zenith angle ${\mu }_0$, the cosine of the viewing zenith angle $\mu$ and the single scattering angle $\theta$ for the reflected light:

$$\begin{gathered} {\mu }_0=\cos {\theta }_0, \mu =\cos {\theta }_v \\ cos \theta =-cos {\theta }_0 cos {\theta }_v+sin {\theta }_0 sin {\theta }_v cos {\varphi }_r \end{gathered}$$

(1)

It follows from Figure 1 that relative azimuthal angle (RAA) is equal to zero in the glint region (${\theta }_0={\theta }_v$). RAA is equal to π in the exact backscattering direction and at ${\theta }_0={\theta }_v$. Such a definition is usually used in the radiative transfer studies [14]. RAA derived from OLCI data differs from that given above by π. This difference must be accounted for while applying the theory described below to OLCI data.

Figure 1. Definition of the solar zenith angle ${\theta }_0$, solar azimuthal angle ${\varphi }_0$, viewing zenith angle ${\theta }_v$, viewing zenith angle ${\varphi }_v$ and relative azimuth angle ${\varphi }_r$. Illustration after Hudson et al. (2006) [15].

Figure 1. Definition of the solar zenith angle ${\theta }_0$, solar azimuthal angle ${\varphi }_0$, viewing zenith angle ${\theta }_v$, viewing zenith angle ${\varphi }_v$ and relative azimuth angle ${\varphi }_r$. Illustration after Hudson et al. (2006) [15].

2.2.2. Reflectance, Spherical and Plane Albedo

The top-of-atmosphere (TOA) reflectance:

$$R_{TOA}=\frac{\pi I^{\uparrow }}{{\mu }_0{F}_0},$$

(2)

where $I^{\uparrow }$ is the intensity of reflected light, $F_0$ is the spectral solar flux at the top-of-atmosphere incident on the unit area perpendicular to the solar light beam. Many satellite instruments simultaneously measure both $I^{\uparrow }$ and $F_0$ and allow for the derivation of the top-of-atmosphere reflectance. The reflectance depends on atmospheric effects due to molecular, cloud and aerosol scattering and absorption of solar radiation. For accurate retrieval of surface optical properties, these effects must be removed.

The plane albedo $r_p$ is defined as the integral of reflectance R across all viewing azimuth angles ${\varphi }_v$ (assuming that the solar azimuth angle is equal to zero) and zenith ${\theta }_v$ angles:

$$r_p\left({\mu }_0\right)=2{{\displaystyle \int }}_o^1\frac{1}{2\pi }{{\displaystyle \int }}_0^{\pi }R\left({\mu }_0,{\mu }_v,{\varphi }_v\right)d{\varphi }_v{\mu }_{v}d{\mu }_v,$$

(3)

The spherical albedo $r_s$ is found by integration of $r_p\left({\mu }_0\right)$ over all incident angles ${\theta }_0$:

$$r_s=2{{\displaystyle \int }}_o^1{r}_p\left({\mu }_0\right){\mu }_{0}d{\mu }_0.$$

(4)

The plane and spherical albedo can be calculated for the complete atmosphere—underlying surface system. However, we are interested in this work in the albedo of underlying snow surface [16]. In this case the value of $R$ in Equations (3) and (4) corresponds to the bottom-of-atmosphere (snow) reflectance. In some cases, yet another terminology is used. The spherical albedo is also called white sky albedo (WSA) and planar albedo is called black sky albedo (BSA) [17]. The snow albedo measured on the ground becomes independent of the solar zenith angle under cloudy conditions and corresponds to WSA. The snow planar albedo corresponds to the ground reflected solar light measurements for an idealized case, when the contribution of the atmosphere can be ignored, e.g., at 1020 nm for OLCI band 21 for the case of cloudless low-aerosol polar atmosphere. Generally speaking, neither WSA nor BSA, but rather their linear combination is observed on the ground for the case of a cloudless atmosphere. The weighting coefficient between WSA and BSA is often approximated as the diffuse light fraction [18].

2.2.3. Assumptions

The snow surface often has complex morphology from, e.g., sastrugi that resembles frozen water waves or from variations of snow specific surface area and pollution load both in horizontal and vertical directions. Nonetheless, we derive snow properties using the following simplifying assumptions:

Horizontally homogeneous, plane-parallel snow layer;

Vertically homogeneous snow layer;

The slope effects and snow roughness are neglected;

The snow is a semi-infinite layer. Therefore, there is no need to account for the reflective properties of an underlying surface, such as ice or soil;

Close-packed media effects are ignored;

Asymptotic radiative transfer theory valid for weakly absorbing turbid media can be applied to the interpretation of measurements in the OLCI spectral range;

Geometrical optics under weak light absorption approximation can be used to derive local optical snow characteristics in the OLCI spectral range;

Impurities (dust, soot, algae, etc.) are located external to ice grains;

The single light scattering by snow grains is spectrally neutral in the OLCI spectral range.

2.2.4. Pre-Processing

The retrieval process consists of the following pre-processing steps:

Spectral OLCI top of the atmosphere radiance to reflectance conversion using the SNAP Rad2Refl module;

Cloud screening;

Ground scene classification.

Cloud screening is performed using Simple Cloud Detection Algorithm (SCDA) algorithm [19,20]. A land mask is used. Therefore, ground pixels containing oceanic surfaces are not processed. If R (400 nm) is smaller than 0.2, the retrievals are not performed because ground scenes are dark, i.e., not covered by snow or bare ice. The normalized difference bare ice index (NDBI) is calculated as:

$$NDBI=\frac{R\left(400 \mathrm{nm}\right)-R\left(1020 \mathrm{nm}\right)}{R\left(400 \mathrm{nm}\right)+R\left(1020 \mathrm{nm}\right)}$$

(5)

Bare ice is classified in two steps where NDBI < 0.65 and R (400 nm) < 0.75 with a bare ice index set to 2. Then, for the cases without the bare ice flag = 2, the bare ice flag is assumed to be equal to 1 (100% bare ice covered pixel), occurring when the normalized difference snow index (NDSI) > 0.33. The NDSI is:

$$NDSI=\frac{R\left(865 \mathrm{nm}\right)-R\left(1020 \mathrm{nm}\right)}{R\left(865 \mathrm{nm}\right)+R\left(1020 \mathrm{nm}\right)}$$

(6)

In summary, the bare glacier ice index is, respectively, 1 (2) for the case of clean (polluted) bare ice. Otherwise, the bare ice flag = 0.

The OLCI spectral index (K):

$$K=\frac{R\left(1020 \mathrm{nm}\right)}{R\left(400 \mathrm{nm}\right)}$$

(7)

indicates the range of variability of the snow OLCI reflectance, being closer to 1 for snow surfaces as compared to bare ice. A snow index (I):

$$I=\frac{1-K}{1+K},$$

(8)

where a pixel is considered as a snow-covered one (snow index is equal to 1) if NDSI < 0.1 and R (400 nm) > 0.75. Otherwise, snow index = 0.

The indices (I, K, N) are used for the scene classification. They are provided in the output of the algorithm and calculated after processing of OLCI pixels with respect to the determination of their microstructure, albedo, and impurity characterization. The indices are used for the assessment of the quality of the retrievals.

2.2.5. Post Processing and Retrieval Quality Checks

The results for ground pixels with root-mean-square deviation (RMSD) between initial and modelled (using retrieved parameters) OLCI reflectances (at 16 OLCI channels almost free of gaseous absorption, see Table 1) are larger than 5% and are not provided in output in the postprocessing stage. They are assumed to be covered by other types of surfaces or partially covered by clouds. The same is true for the pixels, where retrieved total ozone column (TOC) differs by more than 25% from the TOC provided in OLCI files. The diameters of ice crystals in clouds and snow occupy different sized bins. If the effective grain diameter (EGD) is smaller than the threshold value (THV) of 0.14 mm, then it is assumed that the ground scene is covered by clouds or diamond dust. All values of THVs mentioned in this section can be changed depending on the user requirements.

2.3. Snow Property Retrieval

2.3.1. The Retrieval of Snow Grain Size and Albedo for Completely Snow-Covered Ground Scenes in Absence of Impurities

Here, we describe the snow property retrievals procedure that accounts for atmospheric scattering effects and can be used both for clean and polluted snow. A priori knowledge on the snow pollution level is not required. It is assumed that the ground scene is completely snow-covered and the measured top-of-atmosphere reflectance at the wavelength λ can be presented as:

$$R_{meas}=R{T}_g,$$

(9)

$$R=R_a+\frac{T_a{r}_s}{1-r_a{r}_s},$$

(10)

where R is the reflectance of atmosphere without account for the atmospheric gaseous absorption, $R_a$ is the atmospheric reflectance for the case of black underlying surface, $r_a$ is the atmospheric spherical albedo for the black underlying surface, $r_s$ is the snow spherical albedo, $T_g$ is the gaseous transmittance, $T_a$ is total atmospheric transmittance for black underlying surface. The functions $R_a, T_a, r_a$ can be calculated using the radiative transfer theory for a given distribution of aerosol particles in atmosphere with account for multiple light scattering effects as described in [10], where the approach for the parameterization of $T_g$ is also given.

Equation (9) is an approximation. While the gaseous absorption effects and multiple light scattering in vertically inhomogeneous atmosphere should be considered in the framework of the solution of the radiative transfer equation [21], we concentrate our snow property retrievals on the OLCI channels that are only weakly affected by the gaseous absorption.

Equation (10) follows from the first principles in the case of Lambertian underlying surfaces [22]. As snow is not an exactly Lambertian reflector, we substitute $r_s$ in the nominator by the value of the snow reflection function $R_s$ which partially accounts for the snow reflectance anisotropy and provides a correct limit for the measured reflectance for the idealized case of absence of atmosphere between satellite and ground ($R$ = $R_s$ by definition in this case). Therefore, we shall assume that

$$R=R_a+\frac{T_a{R}_s}{1-r_a{r}_s}$$

(11)

for an optically thin polar atmosphere with underlying snow surface. The snow reflectance function can be approximated at OLCI channels in the following way [6,23,24,25]:

$$R_s(\mu ,{\mu }_0,{\varphi }_r) = R_0(\mu ,{\mu }_0,{\varphi }_r)r_s^{\xi }$$

(12)

where

$$r_s=exp\left\{-\sqrt{\left(\alpha +\gamma {\tilde{\lambda }}^{-m}\right)L} \right\}$$

(13)

is the snow spherical albedo. The meaning of all parameters and functions in Equations (12) and (13) is explained in

Table 2. Description of parameters in Equations (12) and (13).

Parameter	Meaning	Equation
$\mu$	Cosine of viewing zenith angle
${\mu }_0$	Cosine of solar zenith angle
${\varphi }_r$	Relative azimuthal angle
$R_0$( $\mu ,{\mu }_0,{\varphi }_r$)	The reflectance of snow at the idealized assumption that there is no light absorption in snow
$\xi$	Anisotropy function	$u\left({\mu }_0\right)u\left(\mu \right)/R_0$($\mu ,{\mu }_0,\psi$)
$u\left({\mu }_0\right)$	The escape function	$u\left({\mu }_0\right)=\frac{3}{5}{\mu }_0+\frac{1+\sqrt{{\mu }_0}}{3}$
$\tilde{\lambda }$	The normalized wavelength (${\lambda }_0=1000 \mathrm{nm}$)	$\lambda /{\lambda }_0$
$\gamma$	The impurity load parameter
m	The impurity absorption Angström parameter
α	Bulk ice absorption coefficient	$4\pi \chi /\lambda$
$\chi$	The imaginary part of ice refractive index
L	The effective absorption length

.

It follows from Equations (11)–(13) that the underlying snow–atmosphere reflectance outside gaseous absorption bands is determined by the four underlying snow and impurity parameters ($L,R_0,m,\gamma$), illumination/observation directions and atmospheric characteristics ($R_a,T_a,r_a)$. The snow parameters listed above do not depend on the wavelength [23]. On the other hand, the atmospheric characteristics strongly depend on the wavelength, primarily due to the spectral dependence of atmospheric molecular and aerosol optical thickness, which decreases with wavelength. Both aerosol and molecular optical thickness is smaller than 0.01 at the OLCI wavelengths above 865 nm because of typically clean polar atmospheres. For example, it has been found [26] that the molecular optical thickness is smaller than 0.01 in polar regions at wavelengths larger than 850 nm. The aerosol optical thickness is smaller than 0.01 at 870 nm at Dome C in Antarctica [27].

Therefore, in a reasonable approximation one can assume that OLCI reflectance at channels 865 nm and 1020 nm are insensitive to atmospheric light scattering effects. These channels are outside major gaseous absorption bands. Importantly, light absorption by snow grains is much larger as compared to light absorption by small amounts of pollutants in polar snowpack. Therefore, we can assume that the second term in Equation (13) can be ignored in the near IR (above 800 nm wavelengths). With neglecting aerosol effects for near IR, the parameters ($L,R_0$) can be obtained analytically from Equations (12) and (13). We can thus write for the measured reflectance at the pair of near IR OLCI channels only weakly affected by the atmospheric scattering and absorption effects in clean polar atmospheres as:

$$R_{meas,1}=R_0 exp(-\xi \sqrt{{\alpha }_{1}L}), R_{meas,2}=R_0 exp(-\xi \sqrt{{\alpha }_{2}L}),$$

(14)

where indices signify the OLCI bands. In this work, we have used the OLCI channels at ${\lambda }_1=865 \mathrm{nm} \mathrm{and} {\lambda }_2=1020 \mathrm{nm}$ that are almost free of gaseous absorption. The spectral calibration can be applied for all OLCI channels ahead of retrievals (see Appendix A).

It follows from Equations (14) that neglect weak atmospheric scattering effects:

$$R_0=R_{meas,1}^{ϵ}{R}_{meas,2}^{1-ϵ}, L=W{\xi }^{-2}l{n}^2(\frac{R_{meas,2}}{R_0})$$

(15)

where $ϵ=\frac{1}{1-b}, b=\sqrt{\frac{{\alpha }_1}{{\alpha }_2}}$, W = 1/${\alpha }_2$. Taking into account the data for ice refractive index [28] one derives that $ϵ=1.55$ and $W=36.08 \mathrm{mm}$.

Thus, OLCI measurements at two near infrared channels enable determining the parameters $R_0$ and $L$ and in turn one can also estimate spectral spherical albedo of clean snow at any wavelength using Equation (13) at γ = 0. Further, one can determine the spectral bottom of-atmosphere (BOA) clean snow surface reflectance (see Equation (12)) and the plane albedo $r_p$ can be derived as [23,25]:

$$r_p=r_s^{u\left({\mu }_0\right)}$$

(16)

Conveniently, the approach for the determination of the snow albedo described above and introduced in [6] does not require simultaneous or daily/weakly (at different observation and illumination conditions) measurements of angular TOA reflectance patterns. The approximation greatly simplifies the retrieval algorithm, enabling computationally efficient estimation of snow properties.

The effective absorption length (EAL) L can be used to determine the effective grain diameter [2,23,29,30] and the snow specific area (SSA) [7,31,32]. In particular, it follows for the effective grain diameter [7] that $d=\frac{9L}{16S},$ where S is the scaling constant, which is equal to the ratio of the absorption enhancement parameter B and co-asymmetry parameter 1-g. The dependence of the pair (B, g) and the constant S on the wavelength and size of ice grains in the visible and NIR [23] can be neglected. However, both parameters depend on the shape of ice grains and, therefore, on snow type, i.e., dry or wet, old or fresh, including combinations. The value of S = 9.2 for the natural snow covers has been proposed [29]. Considering the corresponding uncertainties, we assume the value of S = 9 in this work. This leads to a simple relationship between the values of L and d: L = 16d (see also

Table 3. The relationships between the values of EGD, SSA and EAL.

Quantity	Symbol	Equation	Parameters
EGD	d	pL	p = 0.0625
SSA	σ	q/L	q = 0.1047 m3/kg

). The snow specific surface area σ can be related to the value of the snow grain diameter d [23]: $\sigma =\frac{6}{\rho d},$ where ρ = 0.917 $\mathrm{g}/{\mathrm{cm}}^3$ is the ice density. We can also write: $\sigma =\frac{96}{\rho L}$ (see Table 3).

The broadband planar albedo is defined as the integral of the reflected spectral radiance over the 300–2400 nm wavelength range, divided by the integral of the solar flux over the same range. For a given wavelength, the reflected radiance is equal to the incoming solar flux multiplied by the planar albedo, which itself depends on the absorption length $L$ and on the value of the escape function for the viewing angle $u\left({\mu }_0\right)$. Integrating this product for all pixels is computationally expensive, so the broadband planar albedo of clean snow was parameterized by the following equation as proposed in [33]:

$$r_b=a+b exp(-u\left({\mu }_0\right)\sqrt{vL})$$

(17)

where $a=0.5271, b=0.3612, \mathrm{and} v=0.0235 {\mathrm{mm}}^{-1}$. Therefore, the technique specified above makes it possible to derive all essential clean snow parameters for the case of 100% snow covered ground OLCI scenes.

2.3.2. The Retrieval of Snow Impurity Properties

Let us concentrate on the determination of the pair (m, γ) in the case of snow containing impurities. In the case of polluted snow, the parameters (m, γ) can be derived fitting OLCI measurements in the visible (outside gaseous absorption bands) and calculations according Equations (11)–(13). Importantly, the visible OLCI bands are influenced not only by snow properties but also by atmospheric scattering and absorption effects. We can assume that the molecular and aerosol optical thicknesses can be presented after [26,34,35,36] as:

$${\tau }_{mol}\left(\lambda \right)={\tau }_{mol}\left({\lambda }_0\right){(\lambda /{\lambda }_0)}^{-n}$$

(18)

$${\tau }_{aer}\left(\lambda \right)={\tau }_{aer}\left({\lambda }_0\right){(\lambda /{\lambda }_0)}^{-k}$$

(19)

where ${\tau }_{mol}\left({\lambda }_0=1000 \mathrm{nm}\right)$ = ζexp(−h/$h_0$), ζ = 0.008735, n = 4.08, h is the underlying surface height, $h_0$ is the assumed scale height.

Since the aerosol optical thickness tends to be small in the polar regions, it is difficult to retrieve its value accurately. As an approximation, we use a 550 nm aerosol optical thickness and Angström parameter typical for polar regions [8]. Alternatively, one may apply an aerosol climatology [37].

The only unknown parameter in Equation (11) is the snow spherical albedo $r_s$, which can be derived solving the transcendent equation:

$$a{r}_s^{\xi }+b{r}_s-c=0$$

(20)

where $a=T_a, b=r_a\left(R-R_a\right), c=R-R_a$. Equation (20) follows from Equation (11). This equation can be solved analytically at ξ = 1:

$$r_s=\frac{c}{a+b}$$

(21)

and also at ξ = 2:

$$r_s=\frac{b}{2a}\left[\sqrt{1+4ac{b}^{-2}}-1\right]$$

(22)

We solve Equation (20) numerically for the case of arbitrary value of ξ at the wavelengths ${\lambda }_3$ = 400 nm and ${\lambda }_4$ = 490 nm most not influenced by gaseous assumption, when parameters of aerosol and molecular scattering effects can be estimated in advance [10]. Then, we use Equation (13) to derive the parameters (m, γ). It follows from Equation (13), under assumption that light absorption by snow in the visible is negligible:

$$r_{s,3}=exp(-\sqrt{\gamma {\left({\lambda }_3/{\lambda }_0\right)}^{-m}L}), r_{s,4}=exp(-\sqrt{\gamma {\left({\lambda }_4/{\lambda }_0\right)}^{-m}L})$$

(23)

If the derived value of $r_s$ at 400 nm > 0.99, it is assumed that snow is pollution-free and parameters (γ, m) are not retrieved. If $r_s$ < 0.99, then the parameters (γ, m) are retrieved from Equation (23) under assumption that L is known. L is retrieved from measurements at 865 and 1020 nm as explained above. It follows:

$$m=\frac{2lnz}{lnp},$$

(24)

$$\gamma =L^{-1}{\left[{\lambda }_3/{\lambda }_0\right]}^{m}l{n}^2{r}_{s,3}$$

(25)

where $z=\frac{ln{r}_{s,3}}{ln{r}_{s,4}}, p=\frac{{\lambda }_4}{{\lambda }_3}$. The parameters ($R_0,L,$m, γ) can be derived from OLCI measurements at four wavelengths almost not affected by gaseous absorption (400, 490, 865, 1020 nm), enabling derivation of:

Spectral spherical albedo;

Spectral planar albedo;

Spectral bottom of atmosphere reflectance at any observation and illumination geometry;

Effective grain diameter;

Snow specific surface area.

In addition, the polluted snow broadband albedo can be determined using integration of spectral snow albedo derived directly from Equation (20) as described in [11]. This makes it possible to avoid an assumption on the type of pollutants in the retrieval procedure.

Importantly, the value of $m$ can be used to identify the type of pollution. The value of $m$ for black carbon (BC) is close to 1 [5,38,39]. We shall assume that the pollutant is BC, if 0.9 ≤ m ≤ 1.2. Otherwise, the pollution by dust is assumed [40,41,42,43,44]. Given these considerations, the concentration of pollutants can be derived using the following approximate formula [45]:

$$c=\frac{B\gamma }{k\left({\lambda }_0\right)},$$

(26)

where the normalized impurity absorption coefficient at the wavelength 1 μm γ is derived as explained above, c =$\frac{c_p}{c_i}$ is the relative volumetric concentration of snow impurities, $c_p$ is the volumetric concentration of snow pollutants, $c_i$ is the volumetric concentration of ice grains, $B=1.8$ is the absorption enhancement factor [24,30] and $k\left({\lambda }_0\right)$ is the volumetric absorption coefficient of impurities at the wavelength ${\lambda }_0=1000 \mathrm{nm}$ defined as the ratio of average absorption cross section of particles to their average volume [23]. The quality of retrievals of the concentration of pollutants c depends on the availability information on the value of the volumetric absorption coefficient of impurities at the selected wavelength ${\lambda }_0$ = 1000 nm. The derived concentration of pollutants is equal to the ratio of volumetric concentrations of pollutants and ice grains in snow. To determine the relative mass concentration of pollutants $c_m$, one must multiply the value of c by the ratio ζ of densities of pollutants $\left({\rho }_p\right)$ and ice $({\rho }_i)$:

$$\zeta =\frac{{\rho }_p}{{\rho }_i},$$

(27)

where ${\rho }_i=0.917 \mathrm{g}/{\mathrm{cm}}^3$ and we assume that ${\rho }_p=1.9 \mathrm{g}/{\mathrm{cm}}^3$ for black carbon and ${\rho }_p=2.65 \mathrm{g}/{\mathrm{cm}}^3$ for dust, which gives: ζ = 2.1 for BC and ζ = 2.9 for dust. Namely, it follows:

$$c_m=\frac{B\zeta \gamma }{k\left({\lambda }_0\right)}$$

(28)

The parameter ζ can vary depending on the type of BC and dust. In the case of BC, one can assume [23]:

$$k\left({\lambda }_0\right)=\frac{4\pi {\chi }_{BC}\left({\lambda }_0\right)}{{\lambda }_0}D$$

(29)

where the value of imaginary part of BC refractive index ${\chi }_{BC}$ varies depending on the origin and age of BC. We shall assume that ${\chi }_{BC}=0.47$ in the visible. The parameter D depends on the real part of BC particles and their shape. We assume that this parameter is equal to 1.3 for BC [23]. We can now introduce volumetric mass absorption coefficient:

$$k_m\left({\lambda }_0\right)=\frac{k\left({\lambda }_0\right)}{{\rho }_s},$$

(30)

where ${\rho }_s$ is the soot density. The derivation of the value of $k\left({\lambda }_0\right)$ for the case of dust impurities is much more involved. In particular, the value of $k\left({\lambda }_0\right)$ depends on the shape and size distribution of dust particles in snow. It depends on the position of dust particles (inside of ice grains, outside, or both). In addition, refractive index and colour of atmospheric dust depends on the origin of dust. This makes it more difficult to retrieve the dust concentration as compared to the BC concentration retrieval. The value of $k_0={\rho }_d{k}_{abs}^m$$({\lambda }_0)$ for dust aerosol at ${\lambda }_0=1000 \mathrm{nm}$ has been parameterized with respect to the absorption Angström parameter m [45]:

$$k_0={{\displaystyle \sum }}_{i=0}^2{b}_i{m}^i,$$

(31)

where $k_0$ is in 1/mm, $b_0=10.916, b_1=-2.0831 , b_2=0.5441$ and the model of spherical dust particles is assumed. The size of dust particles can be estimated from the value of m as well [45]:

$$d_{ef}={\displaystyle \sum }_{i=0}^2{h}_i{q}^i$$

(32)

where $h_0=39.7373, h_1=-11.8195, h_2=0.8235,$ and $d_{ef}$ is expressed in micrometers.

2.3.3. The Retrievals for Partially Snow-Covered Ground Scenes

The case of partially snow-covered ground scenes can be reduced to the case described above using the linear mixing approximation (LMA). It follows in the framework of the LMA that the bottom of atmosphere reflectance $\Re$ can be presented as a linear combination of reflectances from various portions of the ground scene:

$$\Re ={\displaystyle \sum _{j=1}^N{a}_j}{R}_j,$$

(33)

where $R_j$ is the reflectance of the j-th portion of the ground scene and $a_j$ is the weight of this portion to the total signal and N is the number of various types of surfaces in the ground scene. The value of $a_j$ can be approximated by the ratio of the surface area of the j-th portion to the surface area of the whole ground scene. Let us assume that N = 2. Then, it follows:

$$\Re =f{R}_s+(1-f)R_b$$

(34)

where f is the snow fraction, $R_s$ is the snow reflectance, $R_b$ is the reflectance of the snow-free background surface. We shall consider a special case of partially water or sediment—covered glacier with snow patches. In this special case the last term in Equation (34) can be neglected and one derives from Equation (34):

$$\Re =f{R}_s.$$

(35)

The value of f can be estimated using OLCI measurements at 400 nm: f = $\Re /R_0$ where $R_0$ is the unabsorbing snow reflectance [46]. The derived value of the snow fraction f at 400 nm makes it possible to determine the snow reflectance of the fraction of ground scene completely covered by snow at any wavelength:

$$R_s=f^{-1}\Re .$$

(36)

Then, the clean snow parameters can be derived from the measurements at 865 and 1020 nm (EAL, $R_0$) as discussed above. The derived value of EAL can be used to estimate the albedo of partially snow covered surface. We follow yet another approach in this work and derive the snow spectral albedo (for a portion of pixel covered by snow) using Equation (20) after assumption that the snow reflectance is given by Equation (36).

We assume that the surface type index M is equal to 3 (partially snow-covered ground scene) in the case f is smaller than 0.99. The value of M is equal to 1 for clean 100% snow-covered ground scenes (f is larger than 0.99 and spherical albedo at 400 nm is above 0.98) and it is equal to 2 for polluted (f is larger than 0.99 and spherical albedo at 400 nm is smaller than 0.98) 100% snow-covered ground scenes. The retrievals for partially snow-covered pixels are performed in case the reflectance R_meas (400 nm) is smaller than 0.75. This THV can be adjusted depending on the requirement of users. We do not make an attempt for the retrieval of snow pollution in the case of partially dirty snow-covered ground scenes.

2.3.4. Quality of Retrievals

The quality of retrievals can be assessed making the intercomparison of the retrieved and measured OLCI spectra. The retrieved OLCI spectrum can be calculated as discussed above. Namely, it follows:

$$R_{retr}=(R_a+\frac{f{T}_a{R}_s}{1-r_a{r}_s})T_g.$$

(37)

The quality of retrievals is assessed from the analysis of the values of the RMSD of the OLCI spectral reflectance determination:

$$\delta \left(N\right)=\frac{1}{N}\sqrt{{{\displaystyle \sum }}_{j=1}^N{\left[R_{meas}\left({\lambda }_j\right)-R_{retr}\left({\lambda }_j\right)\right]}^2},$$

(38)

where N is the either number of all OLCI channels (N = 21) or the number of selected OLCI channels almost free of gaseous absorption (N = 16, all OLCI channels except three channels at oxygen A-band and two channels at water vapor band, see Table 1). The relative spectral RMSD $sRMS{D}_{rel}=\frac{\delta \left(N\right)}{{\overline{R}}_N}$ is reported in the output of the retrieval routine. Here, ${\overline{R}}_N$ is the average value of the reflectance. The output is provided if $sRMS{D}_{rel}\le 5\% \mathrm{at} N=16$. In addition, we assume that the retrieved value of EGD $\ge 0.14 \mathrm{mm}$, which enhances the cloud mask routine. We also derive the total ozone column (TOC) using OLCI observations as described in [11]. If the retrieved TOC differs from that provided in OLCI files by more than 12%, the output of snow parameters is not provided (the contamination of pixel by clouds is supposed). The threshold values for the corresponding parameters can be changed, if needed.

3. The Validation of the Algorithm

3.1. Spectral Albedo

We have compared the spectral albedo retrieved using the algorithm described above with ground albedo measurements in the vicinity of DOME C in Antarctica (69.3762S, 139.0162E) on 25 December 2016 using the Autosolex instrument [32,44,47]. The spectral ice refractive index as measured by Picard et al. [47] in the visible and that presented by Warren and Brandt [28] at other wavelengths have been used (see also [23]). A good correspondence of measured and retrieved albedo has been found (see Figure 2). The satellite measurements have been performed at the solar zenith angle 61.5 degrees. However, ground measurements have been performed at the SZA equal to 64.8 degrees. Therefore, the satellite-retrieved snow plane albedo has been corrected for this small difference in SZA using the following equation [23]:

$$r_{p2}\left({\mu }_2\right)={r_{p1}}_{}^{u\left({\mu }_2\right)/u\left({\mu }_1\right)}.$$

(39)

Here, ${\mu }_1$ is the cosine of SZA at the time of satellite measurements and ${\mu }_2$ is the same parameter except at the time of ground measurements. The retrieved EGD for this case is 0.36 mm and BBA is equal to 0.8. Both parameters are characteristic for clean snow surfaces in the vicinity of Dome C station in Antarctica.

The second example of the intercomparison of satellite and ground-based retrievals of the snow plane albedo using the SICE-3 algorithm is shown in Figure 3 for the case of dust-loaded snow at Col du Lautaret site in the French Alps (45°2′4″N, 6°24′18″E). The measurements were performed on 17 April 2018 at SZA equal to 55.27 degrees and were part of the EBONI campaign running from 2016 to 2019 [48,49]. Snow pollution originating from the Sahara has been detected at the site at the moment of measurements [45]. The ground measurements are corrected for the solar zenith angle (41.25 degrees) corresponding to the satellite measurements and also for the slope effects (3 degrees slope). The circles represent satellite measurements of snow plane and spherical albedo with OLCI S-3A gains not applied. The plane albedo retrieved with assumption on OLCI gains [12] is given by stars (see Appendix A). One can see that the retrieved plane albedo (with and without gains) is close to the ground measured albedo with results including OLCI gains closer to the ground measurements (the largest difference of retrieved and measured plane albedo occurs at 1020 nm and it is equal to 6% in the case of gains applied).

The difference can occur due to imperfect atmospheric correction (the aerosol optical thickness has been assumed to be equal to 0.07 at 550 nm, which is a realistic assumption for the elevation of the site (2000 m a.s.l.) and Angström exponent has assumed to be equal 1.3) and also due to the scene inhomogeneity both with respect to the snow grain size and pollution load. The albedo of polluted snow decreases with decrease of the wavelength in the spectral range 400–600 nm as seen in Figure 3. This decrease can be attributed to various pollutants present in snowpack. The parameters of pollutants and snow retrieved by us for the case shown in Figure 3 are listed in the first two lines of

Table 4. The retrieved snow and pollutants parameters at two dates and two locations in European Alps. The results derived for the French Alps have been obtained using spaceborne observations. The results for the Italian Alps [45] have been derived from ground-based spectral albedo observations. The retrievals with the gains applied are given in brackets.

Date	Location	m	γ, 1/mm	L, mm	${\mathit{k}}_{\mathbf{0}}, \mathbf{1}/\mathbf{mm}$	d, mm	${\mathit{c}}_{\mathit{m}}, \mathbf{ppm}$	${\mathit{d}}_{\mathit{e}\mathit{f}}, \mathsf{\mu }\mathbf{m}$
17.04.2018	Col du Lautaret (French Alps)	3.04 (2.16)	$1.53\times {10}^{-4}$ ($3.74\times {10}^{-4}$)	17.5 (23.9)	9.61 (8.96)	1.1 (1.5)	82.6 (217.0)	11.5 (18.1)
17.05.2018	Torgnon (Italian Alps)	2.5	$1.52\times {10}^{-4}$	25.6	9.11	1.6	77.4	15.2

(with and without OLCI gains applied). The last line contains the parameters derived using ground measurements (one month later) in the Italian Alps. One can see that the retrieved parameters are similar. This points to the same aerosol source location (Sahara as suggested by [42,44]). The spectral dust mass absorption coefficient derived by us from spaceborne observations of polluted snow is given in the first two lines of

Table 5. The retrieved spectral mass absorption coefficient (MAC) of dust as retrieved with or without gains by SICE as well as measured in [40,50].

Wavelength, nm	Dust MAC, 10−3 m2 g−1
	SICE without Gains	SICE with Gains	Denjean et al. [50] Caponi et al. [40]
428	47.7	21.1	37
532	24.6	13.2	20
660	12.8	8.3	11 *
850	5.9	4.8	5 *
1000	3.6	3.4	3 *

* derived using the data for the MAC at 428 nm and the absorption Angström exponent for the Saharan dust equal to 2.9 provided in Table 5 of the paper Caponi et al. [40].

. One can see that our data are consistent with those reported in [40]. In particular, we derived that the absorption Angström exponent is 2.2–3.0 depending on the application of OLCI gains, which is close to the value of 2.9 suggested in [40] for the dust from the Sahara.

3.2. Broadband Albedo

In this section we describe the validation of the algorithm for the broadband albedo retrieval using OLCI comparing the retrieval results with local BBA measurements in the framework of the Program for Monitoring of the Greenland Ice Sheet (PROMICE) [51]. The locations of PROMICE stations are given in Figure 4. Horizontally levelled up- and down-facing Kipp & Zonen CNR1 or CNR4 pyranometers record incoming and upcoming solar radiation (in W m⁻²). These data are used to calculate the shortwave broadband albedo. The measurement height is at the sensor boom level of 2.7  m over the ice surface. Short-wave radiation is measured by the pyranometers within plastic meniscus domes, allowing minimal water droplet adhesion. The manufacturer reports that sensor uncertainty is 10 %. In practice, this sensor uncertainty has been found to be 5 % for daily totals in Antarctica [52,53]. Short-wave radiation measurements are corrected for sensor tilt inpost-processing [2]. The intercomparison of the broadband albedo measured by ground-based instrumentation and derived from OLCI measurements is shown in Figure 5. The measurements have been collocated both in space and time domains. We also show the results derived from MODIS (MOD10A1 boadband albedo product, collection 6) and also from the previous version of the algorithm. Only cloudless cases (PROMICE cloud index below 0.33) have been used in the comparisons. One can see that OLCI provides more stable values of BBA as compared to MODIS. The differences of PROMICE and new algorithm BBA values are generally smaller as compared to the MODIS BBA values and also those derived from the previous version of the algorithm (see

Table 6. Broadband albedo (BBA) statistics at three PROMICE AWS (EGP, SCO-U, KAN-U) and evaluation of the BBA retrievals from the current algorithm (SICE_new), previous version (SICE_prev) and MODIS for 2017, 2018, and 2019. N gives the number of collocated measurements by PROMICE, OLCI/S-3 and MODIS/Terra/Aqua (free of clouds).

Year	N	Yearly Average BBA				Yearly Average Relative Difference Relative to PROMICE (%)			Yearly Average Absolute Difference Relative to PROMICE (-)
		PROMICE	SICE_New	SICE_Prev	MODIS	SICE_New	SICE_Prev	MODIS	SICE_New	SICE_Prev	MODIS
EGP
2017	119	0.78	0.79	0.83	0.84	1.2	5.9	7.0	0.01	0.04	0.05
2018	35	0.79	0.79	0.83	0.85	0.1	5.4	7.5	0.00	0.04	0.06
2019	48	0.79	0.79	0.83	0.85	−0.6	5.1	7.9	0.00	0.04	0.06
SCO-U
2017	88	0.51	0.53	0.58	0.54	2.8	13.7	5.3	0.02	0.07	0.03
2018	55	0.61	0.62	0.69	0.66	0.1	13.0	8.1	0.00	0.08	0.05
2019	69	0.60	0.60	0.64	0.62	−0.6	6.7	1.3	−0.01	0.04	0.02
KAN-U
2017	63	0.78	0.77	0.81	0.80	−1.2	3.9	2.9	0.00	0.03	0.02
2018	25	0.79	0.78	0.82	0.83	0.0	4.7	5.8	−0.01	0.04	0.04
2019	42	0.81	0.81	0.78	0.76	−4.7	−1.6	−4.2	−0.06	−0.03	−0.05

and also Figure 5). The average bias of OLCI BBA retrievals and PROMICE ground observations is smaller than 0.05 for almost all cases studied. This difference can occur not only due to experimental errors of ground-based and satellite instruments but also due to the assumptions used in the process of the retrieval and scene inhomogeneity (300 m OLCI measurements versus point PROMICE pyranometer measurements). The presence of the instrument itself may also lead to some decrease of ground-measured BBA. Overall, the absolute accuracy of satellite retrievals is good. The differences between ground-based and OLCI retrievals are inside of the interval of accuracy required for climate models (0.02–0.05) [54,55] on a global scale. A smaller uncertainty of $\pm 0.02$ is required for regional climate simulations [56]. Such an uncertainty is difficult to achieve considering typical experimental errors of measurements and processing software. Nevertheless, an accuracy better than $\pm 0.02$ is achieved for almost all cases specified in Table 5 for OLCI BBA retrievals using new version of the code (without account for OLCI gains).

The sites EGP and KAN_U are located far from the ocean (see Figure 4) and characterized by quite a large value of the average surface albedo (0.8) as derived from different types of measurements (see Figure 5A,B,E,F). The elevation of the KAN_U site is 1840 m a.s.l. The EGP PROMICE station is located at an even higher level (2660 m a.s.l.). We show the absolute BBA differences (satellite–ground) for the three PROMICE sites in Figure 5. We conclude that temporal changes of surface albedo are less pronounced for the EGP site. This is due to its higher location, larger distance to the ocean as compared to the other sites and also due considerable difference in latitude (approximately 9 degrees difference with EGP being more to the North, see Figure 4). Interestingly, PROMICE data provide almost the same BBA (0.8) for the EGP and KAN-U sites. Peng et al. (2018) have found that PROMICE data is biased low as compared to VIIRS satellite albedo retrievals. We have arrived at the same conclusion as far as OLCI retrievals for retrievals over the EGP site are concerned. KAN-U and SCO-U BBA biases could be both positive and negative.

The intercomparison of satellite and ground BBA measurements for the SCO_U site (elevation 670 m a.s.l.) with smaller yearly average ground albedo (0.5–0.6) is shown in Figure 5C,D. During wintertime the site is covered by snow and BBA is close to 0.8. The melt season starts early June at this site, which is well captured by both satellite and ground measurements. The lowest values of BBA (0.3–0.4) are reached in July with further increase in September due to the drop of temperatures and subsequent snowfall events. As seen from Figure 5, the largest disagreement between satellite and ground measurements occur during melting season, when the scene is very heterogeneous and may contain ice/snow, water ponds, crevassing, supraglacial water, cryoconite, algae and bare ground in various proportions. One can see that the new version of the SICE algorithm provides more accurate results for a darker ground as compared to the previous versions with the average biases below 0.02 for years 2017–2019 (0.08 for the previous version of the algorithm). Additionally, biases are smaller as compared to the MODIS BBA retrievals (0.05). Such an accuracy is sufficient for running global climate models and agrees with requirements for other satellite missions with respect to the snow BBA determination (e.g., Joint Polar Satellite System requirement [57]).

The time series of OLCI/S-3 BBA retrievals for all PROMICE stations are given in Figure 6. It follows that the algorithm is robust and produces results for the BBA similar to the ground measurements with correct identification of seasonal trends in ground albedo and mixed pixels. Some differences can be explained by the different spatial measurement scales (especially relevant in melt season) and also by the residual clouds not cleared out by our cloud screening procedure.

4. Examples of Retrievals

4.1. Antarctica

We have applied the algorithm to the OLCI/S-3B data for the selected area in Antarctica (81–82S, 120–130E, 31 January 2020, 22:26UTC). The results of the application of the algorithm and the statistical distributions of the derived parameters are shown in Figure 7. The average values of the retrieved parameters for valid retrievals are given in Table 6. As expected, for the interior of the Antarctic ice sheet, no snow pollution has been found using the algorithm applied. The retrievals are selected using the following retrieval quality criteria: the maximum per-pixel relative spectral RMSD between observed and retrieved TOA reflectance (sRMSD_rel) was set at 5% (at N = 16); the maximum per-pixel relative difference (Δ_TOC) between the retrieved total ozone column and the TOC value from ECMWF provided with the OLCI files was set to 12%, and only the cases with the retrieved EGD > 0.14 mm were considered. Unfortunately, we cannot validate the results given in this section using ground observations. However, it is known that the effective diameter of grains is almost uniformly in the range 0.2–0.4 mm in the interior of ice sheet in Antarctica [58,59]. Our results (see Table 6 and Figure 7) confirm the conclusions of the Gay et al. [58] study. It has been found that the coefficient of variance of EGD is approximately 13% for the huge area studied (more than 100 km along latitude and longitude). The average snow specific area is approximately 22 m² kg⁻¹. It corresponds to the case of aged snow (without distinct ice crystals) according to the classification given in [31]. Such type of snow is common in the interior of the Antarctic ice sheet. Pirrazzini [60] reports the values of BBA equal to 0.8 at Concordia (approximately 75S, 123E) on 26–29 January 1997 and Kuipers Munneke et al. [61,62] report BBA equal to 0.81 at 75S and 0E (Amundsen station at the East Antarctic Plateau), which is similar to the value 0.81 reported in

Table 7. Average, standard deviation and coefficient of variance of the retrieved parameters with the updated SICE algorithm (SICE_New) and with the previous version of the algorithm (SICE_Prev) presented in [8] for the area presented in Figure 7.

Parameter	Average		Standard Deviation		Coefficient of Variance (%)
	SICE_New	SICE_Prev	SICE_New	SICE_Prev	SICE_New	SICE_Prev
EGD (mm)	0.31	0.36	0.04	0.06	13.10	16.70
SSA (m2 kg−1)	21.75	18.78	2.82	4.10	13.00	21.80
BBA (-)	0.81	0.84	0.00	0.02	0.60	2.40

. Therefore, one can see that our findings are in accordance with previous research based on ground observations.

4.2. Greenland

We also retrieved snow properties using old and new versions of the SICE software over Greenland for 27 July (day 210), 2019. The pixels have been selected in a random way in the region 27–67W, 62–80N. The comparisons between retrievals have been performed only for the grain diameters in the range 0.2–0.6 mm. Lower values of EGD could be due to the contamination of pixels by clouds. Larger values could be due to 3-D effects/glitter. The differences in the retrieved values of EGD, SSA, and BBA are presented in Figure 8 as the function of the longitude. It follows from Figure 8 and also from data presented in

Table 8. The average values of the retrieved parameters according to new and old versions of the algorithm.

Parameter	New Version	Old Version	Relative Difference, %
EGD, mm	0.33	0.40	−16.1
SSA, ${\mathrm{m}}^2/\mathrm{kg}$	20.95	17.55	17.5
BBA	0.79	0.81	−2.9

that the derived value of EGD and BBA are reduced (by 16 and 3%, respectively). The value of SSA is larger by 17.5% as compared to the previous version. These differences are due to update in the coefficient of the transformation from EAL to EGD and also due to update of the equation for the BBA. The retrieval results of main parameters (BBA, SSA, EGD) in the vicinity of the EGP site on 2 August 2019 are shown in Figure 9. It follows that all parameters are within the expected range of variability. The coefficient of variance (CV) of retrieved and measured TOA reflectance (for channels not affected by gaseous absorption) is also given in Figure 9. It follows that the CV is below 5%, which shows that our retrievals make it possible to predict not only BOA but also TOA reflectance. The difference between retrieved total ozone column and that provided in OLCI files is given in Figure 9. The areas with large differences of both TOC datasets could be due to the presence of clouds or violations of assumptions used in the algorithm. Therefore, these differences are used in post-processing.

5. Conclusions

This work aimed to present the updated version of the Snow and ICE (SICE-2022) retrieval algorithm based on OLCI Copernicus S-3 observations. The algorithm is based on the asymptotic radiative transfer theory, which provides accurate results as compared to the numerical solutions of the integro-differential radiative transfer equation [7]. The technique can be applied to various optical instruments operating in the visible and near IR regions of the electromagnetic spectrum. The example of retrievals using the Multispectral Imager on board Sentinel-2 is given by Kokhanovsky et al. [63]. The retrievals of spectral and broadband albedo highly relevant to climate change studies [1]. They have been validated using ground-based observations both for polluted and clean snow cases. Generally, the accuracy of observations is higher for clean snow and clean atmosphere. The errors of retrievals increase for polluted snow and also for the case of overlying atmosphere with a large load of aerosol particles. The snow grain size retrievals have not been validated. However, the retrieved SICE grain size is close to that derived by the previous version of the code. The validation results for the grain size for the previous version of the code have been presented in a separate study by Vandecrux et al. [2]. The multi-year SICE retrieval results for the Greenland Ice Sheet derived from OLCI/S-3A/3B observations can be found at http://snow.geus.dk/ (accessed on 1 December 2022). The retrieval algorithm also includes the atmospheric correction scheme over snow, which is of particular importance due to the increased aerosol concentrations in the High Arctic attributable to changing atmospheric transport patterns [64]. The algorithm is available in both Python and Fortran implementation at https://github.com/GEUS-SICE/pySICE (accessed on 1 December 2022) [13]. The list of generated snow products is given in Appendix B.