Pre-operational Sentinel-3 Snow and Ice (SICE) Products

This document describes the theoretical basis of the algorithms used to determine properties of snow and ice from the measurements of the Ocean and Land Color Instrument (OLCI) onboard Sentinel-3 satellites within the Pre-operational Sentinel-3 snow and ice products (SICE) project: http://snow.geus.dk/. The code used for the SICE retrieval and its documentation can be found at https://github.com/GEUS-SICE/pySICE. The algorithms were developed after the work from Kokhanovsky et al. (2018, 2019, 2020).


Introduction
This document describes the theoretical basis of the algorithms used to determine properties of snow and ice from the measurements of the Ocean and Land Color Instrument (OLCI) onboard Sentinel-3 satellites within the Pre-operational Sentinel-3 snow and ice products (SICE) project: http://snow.geus.dk/. The code used for the SICE retrieval and its documentation can be found at https://github.com/GEUS-SICE/pySICE. The algorithms were developed after the work from Kokhanovsky et al. (2018Kokhanovsky et al. ( , 2019Kokhanovsky et al. ( , 2020a. Snow is composed of ice crystals in contact with each other and surrounded by air. Snow can include impurities such as dust, soot, algae (e.g., Skiles et al., 2018), here referred to as 'pollution'. Snow can also contain liquid water. The volume concentration of snow grains is usually around 1/3 with 2/3 of the snow volume occupied by air (Proksch et al., 2016). The concentration of pollutants is often low, that is, below 100 ng/g especially in polar regions (Doherty et al., 2010).
The algorithms described here are dedicated to the retrieval of snow optical properties such as snow spectral and broadband albedo and also snow microstructure (snow specific surface area and effective optical grain size). We propose a snow mask based on the Normalized Difference Snow Index (NDSI) and a technique to retrieve the concentration of pollutants in snow, which is possible only for the cases with relatively heavy (above 1ppmv) pollution load (Warren, 2013).

Ocean and Land and Colour Instrument
Ocean and Land and Colour Instrument (OLCI) is a 21 band spectrometer that measures solar radiation reflected by the Earth's atmosphere and surface with a ground spatial resolution of 300 m (see Table 1). The OLCI swath width is 1270 km. OLCI is installed on both Sentinel-3A and Sentinel-3B satellite platforms operated by the European Space Agency (ESA) in service to the European Union Copernicus Programme. 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.

Generated Products
The products of the SICE algorithms are listed in Table 2. Most of retrievals are based on the measurements at 865 and 1020nm, where the influence of atmospheric light scattering and absorption processes on top-of-atmosphere signal as detected on a satellite over polar regions is weak.

Summary of assumptions
The simplified snow layer model represents snow as a: 1. Horizontally homogeneous plane parallel turbid medium; 2. Vertically homogeneous layer; 3. Semi-infinite layer. Therefore, there is no need to account for the reflective properties of underlying surface.
4. Close packed effects are ignored.
5. Geometrical optics can be used to derive local optical snow characteristics.
7. The single light scattering angular pattern is spectrally neutral in the spectral range given in Table 1. 8. Only pixels completely covered by snow are considered, i.e., pixels with ice and/or partially snow pixels are ignored.
9. The effects of slopes and snow roughness are not accounted for.
The output is provided if the OLCI reflectance at 1020nm is larger than 0.1 and the derived diameter of grains is larger than 0.1mm.

Reflectance, spherical and plane albedos
The top-of-atmosphere reflectance is defined as: where, ↑ is the intensity of reflected light, 0 is the solar flux at the top-of-atmosphere. Many satellite instruments simultaneously measure both ↑ and 0 and allow the derivation of the top-of-atmosphere reflectance. In the absence of cloud, the bottom-of atmosphere reflectance or snow reflectance is defined by Eq. (2.1.2) when applied at the bottom of the atmosphere.
The reflectance depends on atmospheric effects due to molecular and aerosol scattering and absorption of solar radiation. For retrieval of surface optical properties, these effects must be removed.
The plane albedo is defined as the integration of bottom-of atmosphere reflectance R across all viewing azimuth and zenith angles: The spherical albedo is found by integration of R over all incident angles 0 :

Retrieval overview
We first convert the top of the atmosphere radiance to reflectance using the SNAP Rad2Refl module. The top of the atmosphere reflectances , are then corrected for ozone and molecular scattering. Retrievals are thereafter approached in two ways, depending on whether they are considered as clean or polluted snow. The test differentiating these two snow conditions is based on the theory presented in Section Clean snow: given the pixel's illumination and viewing geometry, we calculate the theoretical RTOA at band 1 (λ = 400 nm) that this pixel would have if it had a surface reflectance of 1. If the observed RTOA is higher than this virtual RTOA, the pixel is considered as clean snow. Otherwise it is considered as polluted snow.

Clean snow retrieval approach
If the pixel is classified as clean snow, we derive snow spectral albedo in the spectral range 0.3-2.4 micrometres using the two-parameter analytical equation as described by Kokhanovsky et al. (2019) under assumption that atmospheric effects can be ignored at OLCI channels 16 and 21 (λ = 865 nm and λ = 1020 nm). This simple approach has produced good performance when comparing the retrieved snow spectral albedo and broadband albedo to ground observations (Kokhanovsky et al., 2019).

Polluted snow retrieval approach
The atmospheric correction for the polluted snow is treated in two ways depending on the OLCI reflectance at band 21 (λ = 1020 nm).

Case 1: Polluted pixels with RTOA at band 21 over 0.4
If OLCI reflectance at channel 21 is above 0.4, we assume that scattering and absorption of light by surface impurities and atmosphere can be ignored at the wavelengths 865 and 1020nm but not below 865 nm. Therefore, reflectances of bands above 16 are being retrieved based on the OLCI measurements at bands 16 and 21 (λ = 865 nm and λ = 1020 nm) using the analytical equation for the spectral surface albedo presented by Kokhanovsky et al. (2019). For wavelengths below 865 nm, we account for gaseous absorption and light scattering by aerosol in the framework of the theory described in Section Atmospheric correction.
Case 2: Polluted pixels with RTOA at band 21 under 0.4 In the case of polluted and dark (RTOA at band 21 under 0.4) snow and ice pixels, scattering within the atmosphere affect all OLCI channels. We can no longer assume that pollutants and other effects have small influence on OLCI reflectance above 865nm. Here we need an additional assumption on the intrinsic optical properties of these dark surfaces before the retrieval is performed.
In both cases, the albedo inside O2 and water vapor absorption bands is derived using the linear interpolation of results for neighbouring channels.

Correction of the OLCI TOA reflectance for molecular absorption
The top-of-the-atmosphere reflectance is corrected for ozone absorption using the ozone transmittance function 3 : Where the transmittance function is defined as in Rozanov and Rozanov (2010): Where 1 + 1 0 is the air mass factor and ( ) is the ozone vertical optical depth (VOD) at the wavelength , defined as: Here ( , ) is the ozone absorption cross-section at the height z and the wavelength , 3 ( ) is the concentration of the ozone at the height z. Equation (2.3.3) depends on the vertical profile of and 3 but in the absence of such information, we use the simplification by Kokhanovsky et al., (2020b) where the ozone optical depth is expressed as the product the total column ozone and an reference ozone absorption cross-section Eventually, the transmittance can be calculated for each pixel using Equation (2.3.2) and (2.3.4): Where 3 ̅̅̅̅ is the total column ozone concentration provided by the European Centre for Medium-Range Weather Forecast (ECMWF) within the OLCI files.
We note that one should account for the instrument spectral response function because the measurements are usually performed not at a single wavelength but in narrow spectral range . Therefore, the value of , ( ) will differ for different instruments even if measured at the same central wavelength.

Molecular and aerosol scattering of light and effects on the top-of the atmosphere reflectance
The background atmospheric aerosol in Arctic is usually characterized by the low values of aerosol optical thickness and values of single scattering albedo close to one. Therefore, one can neglect light absorption by aerosol and assume that the atmosphere-underlying surface reflectance (due to molecular and aerosol scattering and reflectance from underlying surface) can be presented in the following way: where the surface spherical albedo and snow reflectance are the quantities we want to quantify in this retrieval. But before that, three characteristics of the atmosphere need to be quantified: the atmospheric reflectance , the spherical albedo of atmosphere , and the total atmospheric transmittance from the top-of-atmosphere to the surface and back to the satellite ( 0 , ).
Before the atmospheric reflectance can be derived, several characteristics of the atmosphere should be described: the molecular and aerosol optical depth, which describe how opaque the atmosphere is at a given wavelength, the phase function and its derivatives: the asymmetry parameter and backscattering fraction. In this section we describe how we derive these characteristics and eventually present the atmospheric reflectance calculation in Section:

Molecular and aerosol optical depth
The atmospheric reflectance depends on the atmospheric optical thickness, which can be presented in the following form: The molecular optical depth can be approximated as this assumption does not lead to the substantial errors.

Phase function, asymmetry parameter and backscatter of the atmosphere
The phase function ( ) of a media define the light intensity scattered by the media at a given wavelength and towards the scattering angle . The phase function is normalized so that integrating ( ) for all gives one. The asymmetry parameter of the media is then defined as the intensity-weighted average cosine of the scattering angle (Hansen and Travis, 1974).
It ranges from -1 for completely backscattered light to +1 for entirely forward scattered light and is defined as: In presence of both molecular scattering and aerosol, the phase function can be presented in the following form: is the molecular scattering phase function and ( ) is the aerosol phase function. We shall represent this function as: (2.3.14) Therefore, it follows for the asymmetry parameter: The parameter varies with the location, time, aerosol, type, etc. We shall assume that it can be approximated by the following equation: It should be pointed that the system of equations given above enables the calculation of underlying snow-atmosphere reflectance as a function of the aerosol optical thickness for a known value of the snow spherical albedo.

Calculation of the atmospheric reflectance, transmittance and spherical albedo
The atmospheric reflectance , caused by coupled aerosol-molecular scattering, can be presented within the framework of the Sobolev approximation (Sobolev, 1975) as the sum of the reflectance due to single scattering and the reflectance due to multiple scattering : Both contributions can be expressed as a function of the atmospheric optical depth ( ) and of the atmosphere's phase function ( , ), and its asymmetry parameter ( , ).
The single scattering contribution is expressed as: The multiple light scattering contribution is approximated as  The atmosphere's spherical albedo is found using the approximation proposed by Kokhanovsky et al. (2007): The coefficients of polynomial expansions of all coefficients (a, b, c, , ) with respect to the value of g are given by Kokhanovsky et al. (2005).

Clean snow
The asymptotic radiative transfer theory relates the reflectance of a medium to its spherical albedo and consequently allows for the determination of spherical albedo using reflectance observations for a given observation geometry. See He and Flanner (2020) for a review of current radiative theories. The framework used here was first developed and verified with airborne measurements of albedo and reflectance over a bright cloud field with the spherical albedo in the range 0.8-0.95 (Kokhanovsky et al., 2007). This relationship was later adapted to snow by Kokhanovsky et al. (2018Kokhanovsky et al. ( , 2019Kokhanovsky et al. ( , 2020a to retrieve the snow spherical albedo from the atmosphere-corrected OLCI reflectance as: where 0 is the theoretical reflectance of snow in the absence of absorption (Kokhanovsky et al., 2019, Appendix A), µ 0 is the cosine of the solar zenith angle, µ is the cosine of the viewing zenith angle, and is the angular function (Kokhanovsky et al., 2019) defined as: where µ is either µ 0 or µ.
The spherical albedo for clean snow can be presented in the following form (Kokhanovsky et al, 2018): Where l is the effective absorption length in snow and α(λ) is the bulk absorption coefficient of ice calculated from the wavelength and , the imaginary parts of ice refractive index (Warren and Brandt, 2008) reported in Section Appendix: Data tables used in the retrieval: The plane albedo can be derived eventually from spherical albedo. Namely, it follows (Kokhanovsky et al., 2019): The effective absorption length l does not depend on the wavelength in the OLCI spectral range as demonstrated by Kokhanovsky et al. (2018). The same is true for 0 , the reflectance of non-absorbing snow layer. Therefore, we can derive both parameters from measurements of OLCI reflectance at two wavelengths 1 and 2 as in Kokhanovsky et al., (2018) that satisfy the following two criteria: i) the reflectance at these channels must be sensitive to the parameters of interest, and ii) these channels must be least influenced by atmospheric scattering and absorption processes. Consequently, the OLCI channels centred around 865 and 1020 nm are the best candidates for the retrieval.
The 0 and parameters are subsequently defined as follows. where the parameter depends on the type of snow and shape of grains. We assume that A=0.06 in the retrievals as suggested by Kokhanovsky et al. (2019). The snow specific surface area is derived as = 6 , (2.4.9) where = 0.917 g cm -3 is the bulk ice density.

Polluted snow and ice
As in the case of clean snow, we assumed that scattering and absorption of light by atmosphere and impurities in snowpack can be ignored at the wavelengths 865 and 1020nm.
This makes it possible to derive the parameters 0 , l, d and in the same way as for clean Where is the measured top-of-atmosphere reflectance function, is atmospheric contribution to the measured signal, is the spherical albedo of the atmosphere, is the bottom-of-atmosphere surface reflectance, is atmospheric transmittance from the top-ofatmosphere to the underlying surface and back to the satellite position, 3 is the transmittance of purely gaseous atmosphere. Given that is measured and that , , , 3 , 0 can be calculated following the approach detailed above, the system presented in Eq. (2.4.10) has therefore only one unknown, , which cannot be presented in closed form. We consequently derive iteratively using Simpson's rule.
For the channels 13-15, affected by oxygen absorption, and 19-20, affected by water vapor absorption, the spherical albedos are linearly interpolated between the retrieved spherical albedo at channels 12 and 16 in the first case and channels 18 and 21 in the second.
For the very dark, polluted pixels ( <0.4 at band 21), we assume that the underlying surface is not snow, the application of the equation (2.4.6) relating the snow albedo to the snow grain size is not justified. In this case, it is assumed that the value of reflectance for a non-absorbing surface can be approximated, as discussed by Kokhanovsky et al. (2019), by the following expression: and θ is the scattering angle in degrees. This formulation of 0 is then used when solving Eq. (2.4.10).

Pollutant characteristics
The concentration of pollutants in snow is estimated using the approach described below. It is assumed that the spherical albedo, solved in Eq. (2.4.10), can also be expressed as in Kokhanovsky et al. (2018): where ( ) is spectral absorption coefficient of impurities, is so-called absorption enhancement parameter for ice grains (Kokhanovsky et al., 2019), is the volumetric concentration of ice grains in snowpack and α(λ) is the bulk absorption coefficient of ice defined in Eq. (2.4.4). We shall assume that Babs=1.6 in the retrieval procedure.
In the visible spectrum, the absorption by ice particle can be neglected ( ( ) ≈ 0) and the polluted snow spherical albedo can be presented in the following form: or ( ) = ln 2 ( ( )) (2.4.14) Let us assume that the impurity absorption coefficient ( ) can also be expressed in following form: where 0 = 1 . The Angström absorption coefficient can then be derived from Eq.
(2.4.14) and (2.4.15) evaluated at 1 = 400 nm and , 2 = 412.5 nm where the spherical albedo was previously derived: = ln ( ln 2 ( ( 1 )) ln 2 ( ( 2 )) ) ln ( 2 1 ) (2.4.16) The Angström coefficient can then be used to characterize the type of pollutant present in the snow. Since soot has a typical Angström coefficient around one while dust's Angström coefficient ranges from 3 to 7, we here assume that the snow is polluted by black carbon if m < 1.2 and by dust pollution otherwise.
In the case of soot, ( ) can be approximated as in Kokhanovsky et al., (2018): Here, is the enhancement ( ) = 4 ( ) is the bulk absorption coefficient of soot , ( ) is the imaginary part of soot refractive index, currently set at a constant 0.47.
In the case of dust pollution, we assume: ( 1 ) = 0.01µ −1 to calculate the relative volumetric concentration of pollutants in snow .

General case
The derived spectral albedo is used to integrate the planar and spherical broadband albedo (BBA) ̄, ( 1 , 2 ) over any wavelength interval [ 1 , 2 ]: where ( ) is the incident solar flux at the snow surface, , ( ) is plane (p) or spherical (s) albedo depending plane or spherical BBA, ̄, ( 1 , 2 ), is to be calculated. Currently, only shortwave spherical/plane BBA ( 1 = 300 nm, 2 = 2400 nm) is being retrieved but additional ranges may be added in the future depending on user demand.
Broadband albedo are only weakly sensitive to the variation of ( ). The spectrum of incident solar flux at the snow surface is therefore assumed to be identical in all pixels and is approximated by the following analytical equation: where 0 = 3.238e+1, 1 = -1.6014033e+5, 2 = 7.95953e+3, = 11.71 μm −1 , and = 2.48 μm −1 . The coefficients have been derived using the code SBDART (Ricchiazzi et al., 1998) in the spectral range 300-2400 nm at the following assumptions.

Approximation used for clean snow
In the case of clean snow, the exact integration of BBA (Eq. (2.3.14)) is possible because the spectral reflectance is known for each of OLCI measurement wavelength. However, this integration is time consuming. To speed up the retrieval process, the shortwave spherical albedo can be directly expressed as a function of the retrieved grain diameter:  Table 4.

Approximation used for polluted snow and ice
For polluted snow, the spherical albedo and planar albedo cannot be expressed in a closed form as it is the solution of the system of equation (2.4.10). Nevertheless, , ( and respectively) are known for each of the wavelength corresponding to OLCI channels. To circumvent this issue, we build functions of the wavelength that approximate the retrieved , over three intervals: 1) Over 400-709 nm, we approximate spherical and planar albedo , by a polynomial of the second order fitted to the retrieved , ( = 400 nm), , ( = 560 nm), and , ( = 709 nm).
These assumptions make it possible to derive the value of BBA analytically.

NDSI and NDBI
The normalized difference snow index (NDSI) is calculated as: = (865 ) − (1020 ) (865 ) + (1020 ) . (2.5.1) A pixel is considered snow-covered if NDSI is in the range 0.3-0.4 and R(400nm) is larger than 0.75. The normalized difference bare ice index (NDBI) is calculated as: = (400 ) − (1020 ) (400 ) + (1020 ) . (2.5. 2) The bare ice is classified in two steps. First, dark bare ice is identified where NDBI is less than 0.65 and R (400nm) is less than 0.75. Then for cases the dark bare ice flag is not set, the bare ice flag is equal to one (100% bare icecovered pixel), if NDSI is larger than 0.33. Also is assumed that the dark dirty bare ice flag is equal to one (100% dark dirty bare icecovered pixel), if NDBI is smaller than 0.65 and R (400nm) is smaller than 0.75 and that a land mask is used.