The Determination of the Snow Optical Grain Diameter and Snowmelt Area on the Greenland Ice Sheet Using Spaceborne Optical Observations

Abstract

The optical diameter of the surface snow grains impacts the amount of energy absorbed by the surface and therefore the onset and magnitude of surface melt. Snow grains respond to surface heating through grain metamorphism and growth. During melt, liquid water between the grains markedly increases the optical grain size, as wet snow grain clusters are optically equivalent to large grains. We present daily surface snow grain optical diameters ($d_{opt}$) retrieved from the Greenland ice sheet at 1 km resolution for 2017–2019 using observations from Ocean and Land Colour Instrument (OLCI) onboard Sentinel-3A. The retrieved $d_{opt}$ are evaluated against 3 years of in situ measurements in Northeast Greenland. We show that higher $d_{opt}$ are indicative of surface melt as calculated from meteorological measurements at four PROMICE automatic weather stations. We deduce a threshold value of 0.64 mm in $d_{opt}$ allowing categorization of the days either as melting or nonmelting. We apply this simple melt detection technique in Northeast Greenland and compare the derived melting areas with the conventional passive microwave MEaSUREs melt flag for June 2019. The two flags show generally consistent evolution of the melt extent although we highlight areas where large grain diameters are strong indicators of melt but are missed by the MEaSUREs melt flag. While spatial resolution of the optical grain diameter-based melt flag is higher than passive microwave, it is hampered by clouds. Our retrieval remains suitable to study melt at a local to regional scales and could be in the future combined with passive microwave melt flags for increased coverage.

1. Introduction

The darkening of snow and ice surfaces in the industrial era has caused the increased absorption of solar energy by the Earth surface [1] and accelerated the mass loss from polar ice caps and glaciers and their contribution to sea level rise [2,3,4]. On the Greenland ice sheet, a highly reflective snow cover accumulates on the surface each winter. Surface snow then undergoes melt in the spring and summer, until it either melts away and exposes a dark underlying glacial ice [5,6], or until the melt stops. The surface melt intensity is governed by sensible heat flux and absorption of solar radiation at the surface [7] and the latter is largely controlled by the snow optical properties: the shape and size of its grains, the presence of water and the concentration of light-absorbing impurities [8,9,10,11].

Among the snow optical characteristics, the size of the ice grains within the top layer of snow has a direct impact on the snow albedo and on how much solar energy is absorbed by the snowpack [9]. Larger grains increase the incoming light-path length within ice crystals and therefore increase the probability for the light to be absorbed. Conversely, snow with smaller grains increases the probability of light exiting the snowpack. Beyond its importance for the snow shortwave reflectivity, snow grain size can be important for passive microwave remote sensing, radar and laser altimetry and snow physical and hydrological modelling. Indeed, microwave emissivity of snow can be used to estimate snow surface melt [12] or snow wetness [13] but such retrievals are complicated by the development of large grains at or under the surface [14]. Knowledge of snow grain size is also required to interpret radar altimetry data as the penetration of laser and radar signals into snow depends on grain size [15,16,17]. Eventually, multilayer snow models are usually used in combination with climate models to simulate surface melt and meltwater infiltration and runoff. These models need surface snow grain size at the surface as boundary conditions. Therefore, snow models would benefit from increased knowledge and observation of the surface grain size [18].

Due to the diversity of snow grain shapes found in natural snowpacks, the geometrical snow grain size is ambiguous and difficult to determine in the field [19,20,21]. On the contrary, the effective optical grain diameter ($d_{opt}$) is defined as the diameter of spheres which have the same volume-to-surface ratio as the nonspherical snow particles [22,23]. Those spheres are usually referred to as optically equivalent spheres. $d_{opt}$ is commonly used in remote sensing [24,25] and as a prognostic variable in snow models [26]. $d_{opt}$ is interchangeable with the snow-specific surface area (SSA) through:

$$SSA=\frac{6}{d_{opt} {\rho }_{ice}},$$

(1)

where ${\rho }_{ice}$ is the ice density (917 kg m⁻³). $d_{opt}$ increases and SSA decreases with time through dry and wet metamorphisms [27]. Under the action of wind, large snow grains can be broken into smaller grains resulting in a decrease of $d_{opt}$ and increase in SSA [28]. When surface melt occurs, some meltwater is held between the surface snow grains. Since water and ice have similar refractive index, this sudden appearance of water around the grain reduces air/ice interface scattering and translates into a sharp increase in $d_{opt}$. The potential of using remotely sensed $d_{opt}$ to map surface melt has been mentioned in previous studies (e.g., [24,29]) but no grain-diameter-based surface melt detection method has been presented to date.

Retrievals of snow grain size and albedo from spaceborne multispectral observations usually rely on snow radiative transfer models [30]. These models can directly, or after inversion, calculate the snow grain size and impurity concentration from reflectance measurements at certain wavelengths given certain assumptions about the snow surface. For instance, numerous models assumed spherical snow grains and used Mie scattering theory and radiative transfer equation to retrieve the snow optical characteristics [21,31,32,33,34,35,36]. The spherical grain assumption was motivated by the successful estimation of spectral hemispherical reflectances of snow with nonspherical ice particles when representing them as spherical grains of a similar volume-to-surface-area [37]. While this technique has been widely applied, it has been noted that the spherical assumption was limited in accounting for the directional variation of snow reflectance [38,39,40,41,42,43,44,45] and therefore would lead to errors when used on remotely sensed directional reflectance. Several models using the nonspherical grains assumption have been applied to snow characteristics retrievals, for example on data from the Sea and Land Surface Temperature Radiometer (SLSTR) onboard the Sentinel-3 satellites [43,44] or on data from the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua satellites [29,38]. The asymptotic radiative transfer theory allows the retrieval of the snow albedo and optical grain diameter for snow with nonspherical grains and has been applied to MODIS data [24,30,46,47,48], AATSR and MERIS data from the ENVISAT satellite [47,49] and OLCI data from Sentinel-3 satellites [50,51]. These studies have carefully presented the theoretical background of their snow retrieval algorithms and validated their output against the available in situ measurements. However, little emphasis was brought to the description of estimated $d_{opt}$, its temporal and spatial variations and its capacity to indicate surface melt.

The detection of surface melt on the Greenland ice sheet has been performed using passive microwave remote sensing [12,52,53,54,55]. The evaluation of these surface melt detection algorithms was nevertheless limited by the scarcity of in situ surface melt estimation. Remotely sensed melt maps consequently relied on positive air temperature periods at few automatic weather station locations and regional climate model output for evaluation [12,54]. Additionally, the spacing of passive microwave observations along a scan is coarse (2.5–25 km, [56]) and requires grid enhancement to achieve higher spatial resolution [55]. The use of optical remote sensing to map surface melt, using the response of grain size to the presence of meltwater, can achieve a much higher spatial resolution (10 m–1 km) than from passive microwave observations.

Here, we present the $d_{opt}$ retrieved via the Pre-operational Sentinel-3 Snow and Ice (SICE) toolchain [57] from the OLCI instrument onboard the Sentinel-3A satellite. We present the SICE $d_{opt}$ dataset in Greenland for 2017–2019 and evaluate it against ground observations. We describe the response of $d_{opt}$ to surface warming and melt using Automatic Weather Stations (AWS) observations and build a surface melt flag based on $d_{opt}$ at a 1 km spatial resolution. Eventually, we compare our melt flag with available melt maps derived from passive microwave measurements.

2. Methods

2.1. OLCI Instrument and Data Pre-Processing

The OLCI instrument, onboard the European Space Agency (ESA) Sentinel-3 A and B satellites, is an along-track multispectral imager recording the Earth’s radiance in the visible to near-infrared spectrum at 21 bands ranging from 400 to 1020 nm with a 2% radiometric accuracy [58]. It is composed of 5 cameras arranged in fan with a combined field of view of 68.6°, producing a 1270 km across-track image swath on the ground. To minimize sun glint, OLCI is tilted across track 12.58° away from the sun. The cameras cumulate 3700 detectors across-track, allowing a spatial resolution of ~350 m at nadir. Scenes containing Greenland between years 2017 and 2019 are identified and the corresponding OLCI L1B products [59] are obtained through the Copernicus Open Access Hub (https://scihub.copernicus.eu/, accessed on 3 January 2022).

We convert the satellite top of the atmosphere (TOA) radiance measurement $L\left(\lambda \right)$ into TOA reflectances $R_{TOA}\left(\lambda \right)$ via SNAP software [60] as:

$$R_{TOA}\left(\lambda \right)=\frac{\pi L\left(\lambda \right)z^2}{z_0{}^2{F}_0\left(\lambda \right)cos{\theta }_s},$$

(2)

where the reference solar irradiance $F_0\left(\lambda \right)$ is from [61], adjusted for the wavelength and sensitivity of each of the OLCI detectors and for the Earth–Sun distance [62]. In Equation (2), $F_0\left(\lambda \right)$ which is measured at solar noon, is scaled by the cosine of the solar zenith angle ${\theta }_s$ and adjusted for the Earth–Sun distance z at acquisition day compared to the reference solar irradiance measurement at the reference distance $z_0$ [63]. $F_0\left(\lambda \right){z_0}^2/z^2$ is provided for each pixel in the L1B product.

2.2. Cloud Identification

Cloud masking is performed after the Simple Cloud Detection Algorithm (SCDA) v2.0 that consists of five tests using reflectances at 0.55 µm (R₁) and 1.6 μm (R₅) and brightness temperatures at 11 µm (BT₈), 12 µm (BT₉), 3.7 µm (BT₇) from the Sentinel-3 Sea and Land Surface Temperature Radiometer (SLSTR) after [64], implemented in the SICE toolchain after [65]. A main cloud discrimination test exploits the brightness–temperature (BT) difference between 11 µm and 3.7 µm that yields negative values for clouds given the strong 3.7 µm solar reflection. A radiometric “vicarious calibration” factor of 1.12 was applied to R₅ after [66].

The Normalized-Difference Snow Index is taken as:

$$NDSI=\frac{R1-R5}{R1+R5},$$

(3)

Then, a pixel is considered cloudy if at least one of the following four tests is positive:

$$[1]R_1>0.30and\frac{NDSI}{R_1}<0.8andB{T}_9\le 290\mathrm{K}$$

[2] BT₈ − BT₇ < −13 K and R₁ > 0.15 and NDSI ≥ −0.30 and R16 > 0.10 and BT₉ ≤ 293 K

[3] BT₈ − BT₇ < −30 K

$$[4]{\mathit{BT}}_8-{\mathit{BT}}_7<\mathit{THR}\mathrm{and}\frac{NDSI}{R_1}<S \mathrm{and}-0.02\le \mathit{NDSI}\le 0.75 \mathrm{and} {\mathit{BT}}_9\le 270 \mathrm{K} \mathrm{and} R_1>0.18$$

with

S = 1.1 if R₁ > 0.75; 1.5 otherwise

THR_max = −5.5 K if (R₁ < 0.75 and BT₉ > 265); −8 K otherwise

THR = min (0.5 × BT₉ − 133, THR_max)

The OLCI TOA reflectances are filtered from the identified clouds with an additional buffer radius of 5 km in effort to remove shadows and reflected illumination due to the presence of clouds.

2.3. Mosaic Construction

Swath data are combined into 1 km EPSG:3413 projection daily mosaics over the Greenland ice sheet. When a region is covered by multiple OLCI scenes the same day, the pixel that is cloud-free and presents the minimal solar zenith angle is used in the mosaic. A lookup table identifies the scene ID of each pixel in each daily mosaic. Grids of solar zenith and azimuth angles, viewing zenith and azimuth angles and ozone are extracted from each OLCI scene and assembled according to the same look-up table as the daily reflectance mosaics [57].

2.4. The SICE Retrieval of Snow Albedo and Optical Grain Size

The SICE $d_{opt}$ retrieval uses the asymptotic radiative transfer (ART) theory [50,51,67]. The ART considers a vertically homogeneous, semi-infinite snow layer seen as a horizontally homogeneous plane parallel turbid medium, where geometrical optics can be used to derive local optical snow characteristics. Only pixels completely covered by snow are considered. Impurities (dust, soot, etc.) are assumed to be located outside of ice grains. Given those assumptions, the ART provides an analytical solution to the radiative transfer equation and relates the snow surface reflectance and albedo to the snow-grain diameter $d_{opt}$ and impurity concentration and type [39,45,68,69].

The SICE retrieval uses wavelengths that are not significantly affected by water vapor, ozone or oxygen light-absorption effects. The TOA reflectances measured by OLCI at ${\lambda }_1$ = 865 nm (band 17) and ${\lambda }_2$ = 1020 nm (band 21) wavelengths are well suited for this purpose and can be used as estimation of surface reflectance after the correction for the minor light absorption by ozone:

$$R_{TOA,cor O_3}\left(\lambda \right)=R_{TOA} T_{O_3}\left(\lambda \right),$$

(4)

where the ozone transmittance $T_{O_3}\left(\lambda \right)$ is estimated from the European Centre for Medium-range Weather Forecasts (ECMWF) total column ozone level provided for each OLCI pixel and a reference ozone optical depth spectrum defined by [51].

The pure snow reflectance $R_s$, which for our two wavelengths of interest equals $R_{TOA,cor O_3}$, can be decomposed into the reflectance R₀ of a nonabsorbing snow surface and a factor that accounts for the light absorption in ice, for the shape of the ice crystals, and for the dependence of reflectance on illumination and observation angles [39,45,50,67,68,69,70]:

$$R_s\left(\lambda \right)=R_{0 }exp\left(-\frac{u\left({\mu }_0\right)u\left(\mu \right)}{R_0} y\left(\lambda \right)\right),$$

(5)

where u, the escape function, is evaluated at the cosine of the solar or observation zenith angles (${\mu }_0$ and $\mu ,$ respectively). We use the following approximation [71]

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

(6)

The similarity parameter y can be presented as

$$y=4\sqrt{\frac{k_{abs}/k_{ext}}{3 \left(1-g\right)}},$$

(7)

for weakly absorbing media in the visible and near infrared regions of the electromagnetic spectrum such as snow. In Equation (7), $k_{abs}$ and $k_{ext}$ are the coefficients of absorption and extinction of the snow and g is the asymmetry parameter defined as the average cosine of the scattering angle inside the snowpack. $k_{abs}$ and $k_{ext}$ can be expressed, using the geometrical optics approximation for weakly absorbing grains [69], as: $k_{abs}=B \alpha \left(\lambda \right) C_v$ and $k_{ext}=3 \frac{C_v}{d_{opt}}$ where C_v is the volumetric concentration of ice grains, d_opt is the effective optical diameter, $B$ is the absorption enhancement factor that accounts for the grain shape and $\alpha \left(\lambda \right)$ is the bulk absorption coefficient of ice. $\alpha \left(\lambda \right)$ is calculated at a given wavelength $\left(\lambda \right)$ from the imaginary part of the ice’s refractive index ($\chi$) compiled by [72]:

$$\alpha \left(\lambda \right)=\frac{4\pi \chi }{\lambda }.$$

(8)

In our formulation, the effective optical diameter d_opt is defined as d_opt = 3v/2s [69], where v is the average volume of ice grains and s is their geometrical cross section perpendicular to the incident light direction (equal to A/4 for convex particles in random orientation, where A is the surface area of particles). For monodispersed spherical snow grains, d_opt equals the physical grain diameter while for other shapes d_opt is the diameter of spheres that have similar volume to surface ratio (also called Sauter diameter). Using the definitions of $k_{abs}$ and $k_{ext}$, Equation (7) can be rewritten as:

$$y=\sqrt{\alpha \left(\lambda \right)\frac{16 B}{9 \left(1-g\right)} d_{opt}}=\sqrt{\alpha \left(\lambda \right) l},$$

(9)

where

$$l=\frac{16 B}{9 \left(1-g\right)} d_{opt},$$

(10)

is the effective absorption length as used by [50,51,67]. It is the effective absorption length $l$ that is retrieved from the OLCI reflectance measurements. Indeed, using the ozone-corrected observed surface reflectance $R_s$ and Equation (5) at ${\lambda }_1$ = 865 nm and ${\lambda }_2$ = 1020 nm we can derive $R_0$ and $l$ as [50]:

$$R_0=R_s{\left({\lambda }_1\right)}^{\frac{1}{1-\sqrt{\alpha \left({\lambda }_1\right)/\alpha \left({\lambda }_2\right)}}} R_s{\left({\lambda }_2\right)}^{\frac{1}{1-\sqrt{\alpha \left({\lambda }_2\right)/\alpha \left({\lambda }_1\right)}}},$$

(11)

$$l={[\frac{ln\left(\frac{ R_s\left({\lambda }_2\right)}{R_0}\right)}{\frac{u\left({\mu }_0\right)u\left(\mu \right)}{R_0}}]}^2 \frac{1}{\alpha \left({\lambda }_2\right)}.$$

(12)

In Equations (11) and (12) above, the snow reflectances $R_s$ are derived from satellite measurements, $\alpha$ can be calculated using Equation (8), $u\left({\mu }_0\right)u\left(\mu \right)$ can be calculated using Equation (6). Eventually, the optical grain diameter $d_{opt}$ can be calculated from Equation (12):

$$d_{opt}=l \frac{9 \left(1-g\right)}{16 B}={[\frac{ln\left(\frac{ R_s\left({\lambda }_2\right)}{R_0}\right)}{\frac{u\left({\mu }_0\right)u\left(\mu \right)}{R_0}}]}^2 \frac{1}{\alpha \left({\lambda }_2\right)} \frac{9 \left(1-g\right)}{16 B},$$

(13)

While $l$ in Equation (13) is derived from the OLCI measurements, the fraction $\frac{B}{\left(1-g\right)}$ cannot be retrieved from the observed reflectances. The value of $\frac{B}{\left(1-g\right)}$ depends on the grain shape and ranges from 7.4 for fractal particles to 11.5 for spherical grains [73]. The author of [74] derived the value of $\frac{B}{\left(1-g\right)}$ from an experiment conducted by [75]: simultaneous measurements of shortwave infrared reflectance and specific surface area of snow samples (spanning from fresh dendritic to aged faceted grains) allowed the calculation of both $l$ and $d_{opt}$ and consequently the ratio $\frac{B}{\left(1-g\right)}$, which had an average value of 9.2. We use this value as in [50]. Similar values were also used by [39].

The SICE retrievals also include surface albedo, which variation can then be compared to the variation in retrieved $d_{opt}$. The planar (blue sky) albedo is the integral of R_s for all viewing zenith and azimuthal angles and for a given angle of incident light from the Sun, i.e. for clear sky conditions, when remote sensing is possible. The planar albedo can be simplified to a function of the bulk absorption coefficient $\alpha \left(\lambda \right)$ and absorption length $l$:

$$r_p=exp\left(-\sqrt{\alpha \left(\lambda \right) l} u\left({\mu }_0\right)\right)$$

(14)

If the retrieved grain size is smaller than 0.1 mm, a residual cloud contamination is possible, and we flag the pixel accordingly. Due to limitation of the ART for low sun conditions, pixels that have a solar zenith angle greater than 75° are not considered for retrieval. We also limit our retrieval to the snow covered part of the ice sheet and for the pixels that have been classified as “clean snow” by the SICE retrieval [51]. An extension of this retrieval for polluted snow pixels [50] is also included in the SICE dataset.

2.5. Comparison with Ground Optical Measurements of Snow Grain Diameter

As part of the surface program at EastGRIP in Northeast Greenland [76], the SSA was measured daily using the IceCube instrument [77] in the summers of 2017 (81 days), 2018 (92 days) and 2019 (65 days). The IceCube device measures the reflectance of a 6.0 cm diameter, 2.5 cm depth cylindrical surface snow sample when illuminated with a 1310 nm laser diode underneath an integrating sphere. The reflectance is then converted into SSA and $d_{opt}$ (see Equation (1)). The snow samples were taken every 10 m along a 90 m-long transect, producing 10 daily samples. The IceCube has a reported accuracy of 10–12% for SSA [77]. The IceCube $d_{opt}$ measurements are compared the retrieved $d_{opt}$.

2.6. Surface Melt and Snowfall Detection

2.6.1. Automatic Weather Station Data and Surface Energy Balance Modelling

To study the response of $d_{opt}$ to the meteorological forcing and surface melt, we used data from the PROMICE AWSs [78]. Out of the 20 PROMICE AWSs on the Greenland ice sheet, we could not use: (i) the AWSs located in cloudy regions such as eastern Greenland, because retrievals are infrequent in those regions; (ii) the AWSs that are buried during the winter and that cannot capture the onset of melt; (iii) the AWSs located in regions where either a thin snowpack or intense melt expose the underlying glacial ice shortly after the onset of melt. The highest elevation AWSs (KPC_U, KAN_U, EGP and CEN,

Table 1. Automatic weather station used for the investigation of surface melt.

Station	Latitude (deg. N)	Longitude (deg. E)	Elevation (m a.s.l.)
KAN_U	67.0003	−47.0243	1840
KPC_U	79.8345	−25.1665	870
EGP	75.6247	−35.9748	2700
CEN	77.1826	−61.1127	1886

) were therefore selected and provided data over 3 years (2017–2019) to evaluate the response of $d_{opt}$ to surface melt.

The surface melt rates are calculated using the GEUS surface energy balance (SEB) and firn model [79]. This used, as input, a gapless time series of air temperature, humidity, wind speed, downward and upward shortwave radiation and downward longwave radiation, along with instrument heights and snowfall. To increase the coverage of the AWS during 2017–2019, we gap-filled the AWS data with adjusted data from the nearest cell of the regional climate model RACMO2.3p2 [4]. We did not account for rainfall as no in situ data are currently available. Snowfall was calculated from increments in surface height and adjusted to springtime snow pit measurements as detailed in [79,80].

The SEB model initially developed by [81] was then used and evaluated in [79,80,82]. The model closes the energy budget iteratively by adapting the surface temperature. The budget at the surface is defined as the sum of downward and upward shortwave radiation, downward and upward longwave radiation, latent and sensible heat fluxes, and subsurface conductive heat flux. The first three energy fluxes are given as input while the other energy fluxes are calculated as a function of surface temperature and other input meteorological variables, according to, respectively, the Stefan–Boltzmann Law, Monin–Obukhov similarity theory and Fourier law across the top layer of the firn model. Subsurface shortwave radiation penetration is neglected, as the top layer of the snow model has a minimum water equivalent (w.e.) thickness of 4 cm. If it is not possible to find a subfreezing surface temperature that nullifies the sum of energy fluxes, then surface temperature is set to 0 °C and the sum of all energy fluxes is then used to melt surface snow or ice. The SEB model is coupled to the GEUS firn model, a multilayer snow and firn model that calculates the temperature, density, grain size and water content for each model layer. At each time step, the model column is updated for snow accumulation, temperature diffusion, firn compaction, grain growth and meltwater infiltration. Further SEB and firn model details appear in [79,80,82] and references therein.

2.6.2. Passive Microwave Remote Sensing of Surface Melt and Melt Flag Comparison

To map surface melt over the Greenland ice sheet between 2017 and 2019, we used the MEaSUREs Greenland daily surface melt flag [12,83]. It used, for that period, data from the SSM/I F17 and F19 sensors on board the Defense Meteorological Satellite Program (DMSP) satellites [84] and a snowpack microwave emission model and dynamic threshold adjusted for each year to detect melting/nonmelting surface conditions. It is available at a 25 km resolution on an EASEv2 grid and only fully glaciated pixels are considered.

For each year, between 1st of May and 30th of September, we compared the detected melt from remotely sensed $d_{opt}$ and MEaSUREs melt flag to the reference melt calculated at the PROMICE AWS, using four metrics. The yearly coverage of a melt flag is the fraction of days between 1st of May and 30th of September for which a flag value is available. The accuracy of a melt flag is the fraction of available days that is accurately classified, either as melt or nonmelt. The omission error is the fraction of true melt days, as estimated from AWS data, that are not detected by a given flag. The commission error is the fraction of true nonmelt days, as estimated from AWS data, that are classified as melt by a given flag.

3. Results

3.1. Evaluation of the Grain Size Rretrieval

The SICE retrievals of $d_{opt}$ ranged from 0.10 to 0.47 mm while the measured $d_{opt}$ derived from in situ SSA measurements at EGP ranged from 0.08 to 0.30 mm corresponding to SSA between 22 and 85 m² kg⁻¹ (Figure 1). The retrieved $d_{opt}$ is on average 0.06 mm (30%) larger than the in situ measurements (

Table 2. Comparison statistics between the measured and retrieved snow optical grain diameter at EastGRIP: Mean difference (MD), Root Mean Squared Difference (RMSD), Pearson’s correlation coefficient (r) and number of days used for comparison (N).

Year	MD (mm)	RMSD (mm)	r	N
all	0.06	0.17	0.68	145
2017	0.07	0.20	0.62	59
2018	0.05	0.14	0.59	57
2019	0.05	0.16	0.50	29

, Figure 1). The daily standard deviation of measurements (error bars in Figure 1) relates to the spatial variability of $d_{opt}$ along the 90 m transect and has an average of 0.03 mm. The spatial variability within the 90 m transect is therefore insufficient to explain the difference between the in situ measurements and retrieved $d_{opt}$. This difference could also be due to the deeper penetration of light within the snow at the 865 nm wavelength used for the SICE retrieval (down to ~3 cm, Figure 1 in [85]). The IceCube device uses a laser at a wavelength (1310 nm) corresponding to a sampling depth of less than 2 cm below the surface. In springtime, $d_{opt}$ is minimal at the surface and increases with depth when snow grains have time to undergo metamorphism [29,86,87]. Consequently, a systematic difference of sampling depth in the SICE retrieval and IceCube measurement will lead to a systematic difference in the grain size they estimate. Measurements and retrievals show linear correlation with Pearson’s correlation (r) ranging from 0.50 to 0.62 depending on the year (Table 2).

A short and low intensity (2 mm w.e.) surface melt episode on 22 June 2017 is associated with an increase in both the ground observation and the satellite retrieval optical snow grain diameter (Figure 1). The ground observation and retrieved $d_{opt}$ continued to increase for the 2–3 following days when melt was not calculated but warm temperatures continued to heat the surface. This melt event is small, with less than 2 mm w.e. for the full day. Such surface melt rates represent minor changes at the surface and limited generation of liquid water. Consequently, the $d_{opt}$ increases on 22 June 2017 were only slightly more pronounced than other increases that occurred during periods when melt was not calculated at the weather station. Major snowfall events recorded by the AWS led to a decrease of grain size by a factor of two or more (Figure 1). In some occurrences (e.g., 18 July 2018), the $d_{opt}$ decreased without the AWS recording any precipitation. A possible cause is wind transporting small grains to the site, as also hypothesized by [29]. However, no clear covariation of $d_{opt}$ and wind speed measured at the AWS could be found at EastGRIP.

3.2. The SICE Snow Optical Grain Diameter Dataset

The SICE toolchain provides daily snapshots of optical snow grain diameter at a 1 km resolution. The main challenge for shortwave optical remote sensing is the impossibility of retrieval during cloudy conditions. The average cloud-free coverage was 45% in 2017, 31.8% in 2018 and 40% in 2019 (Figure 2). The ice-sheet-wide snow-covered area-only mean $d_{opt}$ was around 0.3 mm at the beginning of the season and peaked above 0.6 mm around late July or early August, when surface melt is widespread. The smaller standard deviation in May indicates the relative homogeneity of cold, dry snow grains at the surface in the early melt season. As melt starts on the southern and low-lying ice sheet, the standard deviation increases and reaches its maximum synchronously with the average $d_{opt}$, at the peak of the melt season. The warm year of 2019, during which melt was well above average [88], also had the highest Greenland-wide average $d_{opt}$: 0.45 mm compared to 0.36 mm for both 2017 and 2018. That year (2019) also has the maximum daily mean value (1.08 mm on 31 July 2019; Figure 2), as compared to a respective maximum of 0.71 mm and 0.81 mm in 2017 and 2018.

During the melting season, the melt is triggered initially in southern and low-elevation areas before it spreads northward and to more-elevated areas. $d_{opt}$ follows a similar evolution (Figure 3 and Supplementary Figure S1). At the beginning of the season, before the melt onset, the retrieved $d_{opt}$ is relatively low and homogeneous across the ice sheet (e.g., low standard deviation of grain diameters in May 2018, Figure 2) and does not show a visible gradient with elevation (e.g., northern and western regions in Figure 3 and Supplementary Figure S1). Where and when surface warming and melting started, the lowermost areas saw an increase of $d_{opt}$ and of its spatial heterogeneity (Figure 3 and Supplementary Figure S1), while higher areas still had a relatively low and homogeneous $d_{opt}$. In the southern areas and in warm years such as 2019, the grain diameter spatial gradients are visible already in May, when the first retrievals are available (e.g., southwestern regions in Figure 3). As the surface heating and melt progresses to higher elevations, the $d_{opt}$ increase propagates to the ice sheet interior (Figure 3 and Supplementary Figure S1).

3.3. Optical Grain Diameter as an Indicator of Melt Affected Snow

At the four PROMICE AWS (Table 1), the calculated melt rates can be compared to the median $d_{opt}$ retrieved within a 2 km radius of each AWS (Figure 4). The magnitude of melt calculated from the AWSs allows a better understanding of the covariation of $d_{opt}$ and surface melt than at EastGRIP alone (Figure 1). The retrieved $d_{opt}$ reacted closely to the presence and absence of surface melting. At most sites, the $d_{opt}$ had values at or below 0.5 mm at the beginning of the season, when no melt was calculated. At melt onset, $d_{opt}$ jumps above its premelt values at all sites (Figure 4). Some exceptions such as large $d_{opt}$ in the absence of melt or small $d_{opt}$ retrieved during the melt period can be noticed. These punctual decoupling of surface grain diameter and surface melt can be caused by errors in the retrieval, for instance due to cloud contamination or errors in the melt calculation due to inaccurate AWS measurements. A last potential explanation for mismatch cases is that the satellite retrieval may have been taken at an hour of the day when the surface was not melting, (Sentinel 3 overpass time is 10 am local time at nadir), but surface melt occurred later on that day.

The relation between the retrieved snow $d_{opt}$ and the surface melt rate is complex. Pearson’s coefficient of correlation is 0.49, which reveals a co-variation of the two variables although a linear regression is incomplete in explaining the relationship. As an alternative, we look for the $d_{opt}$ value that is most representative of melting conditions and so plot the calculated daily melt as a function of retrieved $d_{opt}$ for our four sites over 2017–2019 (Figure 5). An analysis of variance indicates that a threshold of 0.64 mm in $d_{opt}$ splits the samples into two coherent classes. The samples that have retrieved $d_{opt}$ below 0.64 mm have a median melt of 0 mm w.e. and an average that is not significantly different from zero. The samples that have a retrieved $d_{opt}$ > 0.64 mm have an average daily melt of 7.11 mm w.e., significantly greater than zero and greater than the average of the first class (p-value < 0.001, Figure 5). Additionally, 81.3% of the d_opt retrieval above 0.64 mm corresponds to daily melt greater than 1 mm w.e. and 95.2% of the $d_{opt}$ retrieval below or equal to 0.64 mm corresponds to daily melt lower than 1 mm w.e. This simple threshold approach classifies appropriately the majority of the days, either as melt days, when the retrieved $d_{opt}$ is above that threshold, and non- or low-melt days when the retrieved $d_{opt}$ is below that threshold. This stepwise approach nevertheless misclassifies 18.8% (12 out of 64) of samples as melt days when no or little melt was calculated those days. Similarly, 4.8% (32 out of 660) of days classified as non- or low-melt actually had more than 1 mm of melting calculated from the AWS data.

To assess the robustness of this threshold, we fitted similar stepwise functions to (i) each station-year (except CEN and EGP which did not present sufficient melt), (ii) all samples but removing each station iteratively and (iii) 2017, 2018 and 2019 separately (Supplementary Table S1). The thresholds found in these test cases have an average of 0.61 mm and a standard deviation of 0.10 mm, consistent with the threshold found when pooling all samples. However, the threshold found for these subsets are more sensitive to the errors that occur at each site arising from the difference of footprint between the satellite retrieval and the AWS measurements, erroneous measurements from the AWS that lead to improper melt calculation as well as cloud contamination that leads to improper $d_{opt}$ estimation. We expect that by pooling all these measurements, random noise will be reduced, and a more robust threshold value can be found.

3.4. Application to a Heat Wave in Northeast Greenland and Comparison to the MEaSUREs Melt Flag

We now focus on the 10–18 July 2019 in Northeast Greenland (Figure 6), in the vicinity of the Northeast Greenland Ice Stream (NEGIS). In this area, warm air masses from the northern Atlantic are being pushed onto the Greenland ice sheet through the surface depression created by the NEGIS [90]. As an independent estimate of melting areas, the MEaSUREs melt flag, derived from passive microwave remote sensing, is also considered and compared to our $d_{opt}$-based melt flag.

On the 10th, the majority of the area shows relatively small $d_{opt}$. Only a narrow band below 650 m a.s.l. has a larger $d_{opt}$ (Figure 6). KPC_L is located below that elevation and already showed a too thin snowpack by this point of the season and could not be used for this analysis. KPC_U is located at 870 m a.s.l. and was still snow-covered and untouched by melt as can be seen from the high albedo (~0.8) and the subfreezing air temperature measured before the start of the event (Figure 7). The MEaSUREs melt flag also indicates nonmelting conditions at KPC_U (Figure 6). On the 11th, both melt flags indicate nonmelting conditions at KPC_U. However, the AWS starts to detect minor melt for that day (Figure 7). This amount of melt is either below the detection limit of both melt flags or occurred during a narrow timespan and was missed by both satellites. On the 12th and the following days, the warming continues and the AWS, the $d_{opt}$-based flag and the MEaSUREs flag all indicate melt at KPC_U (Figure 6 and Figure 7). However, different spatial patterns can be seen in the two melt flags. On the 12th, MEaSUREs flag reaches altitudes up to 2000 m a.s.l. The $d_{opt}$-based flag shows a main melting area reaching up to 1600 m a.s.l. and isolated melting patches up to 1800 m a.s.l. On the 13th and 14th, it is in turn the $d_{opt}$-based melt flag that reaches the highest altitudes (up to 2400 m a.s.l.) while the MEaSUREs melting area does not reach higher than 2200 m a.s.l. (Figure 6). On the 18th, the melt stops at higher elevation and continues at below 1600 m a.s.l. According to the MEaSUREs flag and below 1200 m a.s.l. for the $d_{opt}$-based flag. During the entire period, the EGP AWS did not show any sign of melt either as calculated from the AWS data (Figure 1) or according to the two melt flags (Figure 6).

From this specific case, the SICE $d_{opt}$-based flag agrees with MEaSUREs melt flag on the overall spatial variation of the surface melting area. However, substantial differences between the two were highlighted. Unfortunately, no in situ melt data is available in the area where the two melt flags disagree. Additionally, the KPC_U AWS quickly loses its snowpack which limits the number of days for which both $d_{opt}$-based and MEaSUREs melt flags are available. At KAN_U however, melt is frequent, and the snow never melts out, making it an appropriate location to evaluate the two melt flags against the melt estimated from AWS data (Figure 8). As expected, there are more days where retrieval was impossible for the $d_{opt}$-based flag because of cloud cover: retrieval was possible for only 45% of days during May-September of 2017–2019. Considering the AWS-derived melt estimate as ground truth, accurate melt flagging, false positives and false negatives can be found for both melt flags (Figure 8,

Table 3. Comparison statistics of d_opt-based and MEaSUREs melt flags. AWS-derived daily melt greater than 1 mm w.e. is taken as ground reference. The coverage is the fraction of days between 1st May and 30th of September for which a flag value is available. The accuracy is the fraction of days correctly classified by a flag. The omission error is the fraction of true melt days not detected by a flag. The commission error is the fraction of nonmelt days classified as melt by a flag.

	MEaSUREs Melt Flag				${\mathit{d}}_{\mathit{o}\mathit{p}\mathit{t}}$ Melt Flag
Year	Coverage (%)	Accuracy (%)	Omission Error (%)	Commission Error (%)	Coverage (%)	Accuracy (%)	Omission Error (%)	Commission Error (%)
2017	99.2	77.7	65.6	6.7	51.6	90.5	27.8	2.2
2018	94.3	82.6	26.8	12.2	31.1	76.3	81.8	0
2019	98.4	63.3	56.4	0	52.5	60.9	53.2	0
All	97.3	74.4	50.3	7.3	45.1	75.8	51.3	1.1

). The $d_{opt}$-based melt flag has a very low commission error, indicating very few false positives. The higher omission error indicate that some melt events are not flagged, but apart for 2018, omissions of the $d_{opt}$ flag are comparable with the ones of MEaSUREs flag.

4. Discussion

To achieve a reasonable processing time and memory use, we resampled the OLCI images from ~350 m to a 1 km resolution. This spatial resolution is comparable to other snow optical grain diameter or SSA products [24,44]. When using the retrieved d_opt to identify surface melt, the resulting 1 km resolution melt map represents a significant improvement compared to the 25 km resolution of surface melt maps produced from passive microwave data [12,53]. Recent efforts were made to enhance the resolution of the passive microwave observations [55]. Nevertheless, the in situ validation, at one AWS in [55], remained limited and the contamination of coarse passive microwave observations by land near the ice sheet margin remains problematic.

Our approach reveals spatially discontinuous and patchy melt areas (Figure 6). Nevertheless, even at the relatively fine 1 km resolution, melt may be spatially inhomogeneous within a pixel. This spatial heterogeneity complicates the comparison with pointwise AWS observations. Additionally, the nonmelting/melting binary flagging cannot represent the transformation of a melt-free pixel to a pixel where different snow patches can be melting at various intensities and finally to a fully melting pixel. Further, most of the PROMICE stations are located in lower ablation areas where the snowpack rapidly ablates and glacial ice is revealed, reducing the time over which spaceborne retrieval is relevant. CEN and EGP AWS are exceptional by being in the accumulation area and rarely seeing melt. Only two AWS, KPC_U and KAN_U, are located close to the equilibrium line and allow to assess the response of $d_{opt}$ to surface melt. Future in situ estimation of surface melt, from AWS or other techniques, at the equilibrium line altitude could help to further document the relation between $d_{opt}$ and either the presence or the intensity of surface melt.

The threshold that was found to identify melting conditions, $d_{opt}$ = 0.64 mm, corresponds to an SSA of 0.11 m² kg⁻¹. It is consistent with previous values identified for melting surfaces such as 0.6 mm by [29] and ~0.5 mm in [24] (Figure 6 therein). Previous studies also used the optical characteristics of snow, such as albedo, as a proxy for melt and mass loss [91,92]. But such an approach only worked on annual averages. In our study, we have related daily changes in $d_{opt}$ as a response to the presence or absence of meltwater at the surface estimated using the high-temporal-resolution AWS data.

In theory, it is possible that large melt-affected snow grains remain at the surface after the melt stops and consequently that our melt flag would show false positives after the end of the melt event. Although it is a plausible scenario, we do not see this occurring often in our dataset, potentially for two reasons. First, the pervasive action of wind on the ice sheet may quickly erode the surface layer. Small windblown snow grains then cover the melted snow and the $d_{opt}$-based flag is reset to nonmelting. Second, the end of the melt period coincides with cooler temperatures and consequently solid precipitation. The reflective fresh snow deposited at the surface has the potential to shut down the melt. In that situation, once the site becomes cloud-free again, our retrieval will identify fine grained fresh snow and nonmelting conditions. We note again that the presented flag only works for snow surfaces with low light-absorbing particle concentrations. More in situ observation will be needed to assess the performance of our $d_{opt}$ retrieval for melting and non-melting polluted snow and to assess whether the same threshold on $d_{opt}$ can be used to identify melt. Little is known about the optical properties of natural glacier ice surfaces, even less under melting conditions, and further research is needed to constrain radiative transfer models in those areas.

The penetration of solar radiation into surface snow implies that the SICE retrieval samples the first centimeters below the surface. The penetration depth also depends on the grain size [86] and shape [40] and on snow density. As a result, the thickness of the layer that is being sampled by our retrieval evolves throughout the season. Our retrieval assumes a vertically homogeneous snowpack, but since the near-surface snow grain size may vary vertically [21], multilayer snow models may be more suitable [38]. Nevertheless, multilayer radiative transfer models are computationally expensive and not convenient for near-real-time retrievals. Our retrieval also assumes a nonspherical grain shape (characterized by a shape factor $\frac{B}{\left(1-g\right)}=9.2$ as in [50]). Although this assumption has been producing satisfactory results, it is far from the diversity of snow grain shapes in natural snowpacks [20,43,44]. Recent remote sensing efforts have included diverse grain shapes within the radiative transfer model [38,43,44]. Unfortunately, in situ measurements of grain shape at the relevant spatial scale for remote sensing product evaluation are not readily available to this day. In this absence of ground truth, the grain shape could only be used as a tuning parameter within radiative transfer models to minimize the difference between calculated and observed surface reflectances.

A recent report from the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) evaluated the accuracy of the OLCI instruments onboard Sentinel-3A and B when measuring ocean reflectance [93]. They concluded that the OLCI measurements of surface reflectance at wavelength 865 nm were within the 2% accuracy defined in the mission requirements. However, the measured reflectances at 1020 nm were ~9% and ~6% too bright for Sentinel-3A and Sentinel-3B, respectively. These biases have been derived over dark ocean surfaces and cannot be used over bright snowy surfaces. No reference dataset is currently available for a vicarious calibration of OLCI over bright surfaces. The magnitude of this systematic bias should be subject of future research. However, since our melt detection threshold is chosen empirically, the melt maps derived here should remain unchanged by any future change of calibration, as the threshold could be accordingly re-evaluated.

When comparing our retrievals at 1 km resolution to in situ measurements, either for $d_{opt}$ or for surface melt, the spatial representativity of these observations can be questioned. At EastGRIP, the surface type is relatively homogeneous, stemming from a smooth topography: the elevation does not vary by more than 7 m and the slope by more than 0.24° in the 1 km² surrounding the station (derived from the GIMP DEM, [89]) suggesting uniform temperatures and incident shortwave radiation. At this site, the standard deviation of the $d_{opt}$ along the 90 m observation transect, representative of the spatial heterogeneity, was 0.03 mm, which is significantly smaller than the difference between the observed and retrieved $d_{opt}$: 0.06 mm on average. Capturing the spatial heterogeneity in these dry snow regions is therefore less important than addressing the assumptions made in the retrieval procedure (e.g., grain shape or radiometric calibration). In the areas where surface melting is more prevalent, no in situ observation of $d_{opt}$, and of its spatial heterogeneity, is currently available for the Greenland ice sheet. Regarding the detection of surface melt, no observation currently allows bridging between point estimations at AWS sites and microwave melt flags at 12–25 km resolution [12,54]. We consequently compare point observations to gridded retrievals and let to future work bridging these scale gaps such as recently done for surface albedo measurements [94,95].

5. Conclusions

We presented the optical snow grain diameter ($d_{opt}$) retrieved for 3 years of Sentinel-3 OLCI observations over the Greenland ice sheet. After ensuring that the remotely sensed $d_{opt}$ matches with in situ observations at EastGRIP, we analyzed the response of retrieved $d_{opt}$ to surface warming and melt at four PROMICE automatic weather stations where melt can be calculated from meteorological measurements. We established that the retrieved $d_{opt}$ and the in situ estimation of melt covary, which leads to a binary melt flag that can be constructed around the threshold $d_{opt}$ of 0.64 mm with a high likelihood of melting above this threshold. We applied this threshold to our $d_{opt}$ dataset and derived melt maps for daily 1 km mosaics of Greenland. We found that the spatiotemporal evolution of our melt flag compares well with the passive microwave MEaSUREs melt flag during a heat wave in Northeast Greenland and at the KAN_U PROMICE AWS. Although the two flags identify the same general periods and areas as melting, some spatial and temporal mismatches remain. For the dates when $d_{opt}$ and the melt flag could be retrieved, it showed comparable accuracy (fraction of correct flags, 76%), omission (fraction of false negative, 51%) and commission (fraction of false positive, 1%) as the MEaSUREs melt flag (74%, 50% and 7% for the same metrics). Although our retrieval is not possible under cloud cover, which limits the coverage of the dataset, its 1 km spatial resolution is a great improvement compared to the MEaSUREs melt flag. $d_{opt}$ can therefore be used at a regional scale or in combination with passive microwave to describe melt dynamics on the Greenland ice sheet at high spatial resolution.