Greenland bare-ice albedo from PROMICE automatic weather station measurements and Sentinel-3 satellite observations

The Programme for Monitoring of the Greenland Ice Sheet (PROMICE) provides surface meteorological and glaciological measurements from widespread on-ice automatic weather stations since mid-2007. In this study, we use 105 PROMICE ice-ablation time series to identify the timing of seasonal bare-ice onset preceded by snow cover conditions. From this  collection, we find a bare-ice albedo at ice-ablation onset (here called bare-iceonset albedo) of 0.565 ± 0.109 that has no apparent spatial dependence among 20 sites across Greenland. We then apply this snow-to-ice albedo transition value to measure the variations in daily Greenland bare-ice area in Sentinel-3 optical satellite imagery covering the extremely low and high respective melt years of 2018 and 2019. Daily Greenland bare-ice area peaked at 153 489 km2 in 2019, 1.9 times larger than in 2018 (80 220 km2), equating to 9.0% (in 2019) and 4.7% (in 2018) of the ice sheet area.


Introduction
The recent net loss of Greenland land ice is among the largest contributors to global sea-level rise (Box & Sharp 2017). While warm air advection produces the highest daily ice-ablation observations (Fausto et al. 2016), absorbed sunlight is the largest melt energy source over seasonal time scales (van den Broeke et al. 2008;Box et al. 2012;Fausto et al. 2016). Absorbed sunlight increases during the melt season as surface conditions shift from a highly reflective, dry snow cover, to lower albedo wet snow with larger grains ( Wiscombe & Warren 1980;Brun 1989), and yet lower albedo across the ablation area. Bare ice darkened primarily by ice algae (Stibal et al. 2017;Ryan et al. 2018;Cook et al. 2020;Williamson et al. 2020) plays an important role in peak ice-sheet melt rates. Accurate definition of bare-ice albedo at ice-ablation onset (hereafter bare-ice-onset albedo) has applications in (1) classifying the bare-ice area over large areas of the ice sheet (Ryan et al. 2019;Fausto et al. 2020); (2) constraining Polar regional climate models used to estimate the surface mass balance of the Greenland ice sheet  (3) climate monitoring (e.g. Moon et al. 2020).
Here, we study the surface climate conditions spanning the melt season transition from dry snow to bare ice using PROMICE ground measurements. Our main objective is to determine an albedo value useful in classifying the boundary between seasonal snow cover and bare ice. We proceed to determine the spatial and temporal patterns of bare-ice albedo using spaceborne www.geusbulletin.org observations from the EU Copernicus Sentinel-3 satellite mission applied to a low (2018) and high (2019) melt year to measure the maximum relative differences in bare-ice area.

PROMICE surface measurements
Across the Greenland ice sheet, the sunlight reflectivity of snow and ice, hereafter albedo, and several other surface meteorological and glaciological parameters are measured by more than 20 Automatic Weather Stations (AWSs), operated by The Programme for Monitoring of the Greenland Ice Sheet (PROMICE) since mid-2007 (Ahlstrøm et al. 2008). PROMICE operates AWSs in nine regions around Greenland (Fig. 1) where most locations have a lower and upper AWS, denoted by _L and _U, respectively.
In our analysis of daily average PROMICE AWS data, we used an initial dataset of 225 station years from 26 station locations for air temperature, snow thickness, ice ablation and albedo. Three stations situated in the accumulation area, where the underlying ice did not appear, were excluded from the analysis. The acquisition and/or computation of each variable is described in the following sections.

Air temperature
Air temperature is recorded at PROMICE AWSs using a platinum resistance thermometer in an aspirated shield. The measurement height above the surface varies between 0 and 2.6 m due to snow accumulation, compaction and ablation.

Seasonal snow layer thickness
Snow thickness above the ice surface is obtained from the AWS acoustic recordings of distance from a sonic sensor to the snow or ice surface. The sensor height above bare ice is determined here for each station year from a 20-day average of daily values. This average is computed 10 days after bare-ice onset determined from ice-ablation measurements (see section 2.2) to ensure   snow-free conditions. The sensor height is then subtracted to obtain the snow thickness.

Ice ablation
PROMICE stations measure ice ablation using the pressure of the column of antifreeze over a pressure transducer (Fausto et al. 2012). The transducer is initially drilled 10 to 14 m into the ice and re-installed to avoid complete exhumation. At high ablation (>6 m ice per year) sites like QAS_L and KAN_L, instrument re-installation occurs each year. At most other sites, the re-installation frequency is two to three years.

Surface albedo
Daily average albedo is calculated from 10-minute tiltcorrected (Wang et al. 2015) upward and downward solar irradiance recordings in the 0.3 to 2.5 μm wavelength range. Hourly data are averaged for cases with solar elevation angles above 20°. Daily averages are computed from hourly data between 0 and 1. The daily albedo values are further adjusted after the correction proposed by Aoki et al. (2011) for measurement platform obstruction of the radiometer field of view. This correction increases PROMICE AWS albedo measurements by 0.034 on average.

Determining Greenland bare-ice albedo at ice-ablation onset
In order to study the ablation season albedo as it transitions from snow to bare ice, we use ice-ablation measurements to determine the timing of bare-ice onset. Seasonal snow layer thickness and air temperature further contributed to a better understanding of the evolution of surface conditions. Ice-ablation time series were manually compensated for signal shifts caused by station movement, sensor reinstallation and measurement failure. We then attempted to automatically determine a theoretical candidate date of ice-ablation onset for each station year. To this end, an automatic detection was conducted to identify curve inflexions (Satopaa et al. 2011). However, signal variability as well as heterogeneous transition patterns precluded this approach. Instead, the ice-ablation onset was first identified manually for each station year with supporting data of snow layer thickness and air temperature time series when necessary. Despite careful inspection, roughly 20% of the ice ablation onset dates were flagged as ambiguous because of noisy or complex patterns. Also, distinguishing transient melting days before the ice melt onset from a short sequence of freezing days after the start of ice ablation remained challenging. Measurement uncertainty was not taken into account in this initial and theoretical estimation of ice-ablation onset but helped to better constrain its identification before further refinements.
Following this first step, we excluded time series with interruptions within the period of interest, instrument malfunctions and/or those for which the ice ablation onset could not be identified. Only data within ± 45 days of bare-ice onset were thereafter considered, corresponding to the average time span of the ablation season. The resulting dataset contained 105 station years from 20 stations across the ice sheet (c. 47 and 77% of the initial dataset size, respectively) consisting of more than 9000 daily measurements for each variable.
The raw ice-ablation measurements, consisting of cumulated values after each instrument re-installation, the average instrument recording within 45 days before the initial bare-ice onset were then subtracted from the associated ice-ablation time series for each station year. In this way, the ablation measurements forming the final dataset consisted of values relative to the premelt season and therefore, each station year had a zero average prior to the ice-ablation season. The processing steps are further documented in a Github repository (Wehrlé & Box 2020a).
The timing of bare-ice onset was then refined for each station year by accounting for measurement uncertainty. To this end, a threshold was determined to estimate the first 'significant' ice-ablation value after the manually selected ablation onset; a value for which we have high confidence that bare-ice conditions begin to prevail. A too conservative ice ablation threshold (i.e. above the actual precision of the measurement), would lead to a delayed detection of the bare-ice onset date. Such a delay would cause a dark bias in the determination of the average bareice albedo as the ice may have already been affected by algal darkening (Stibal et al. 2017;Ryan et al. 2018;Cook et al. 2020;Williamson et al. 2020). On the other hand, bare-ice onset would be defined prematurely if the threshold was too restrictive (i.e. below the actual precision of the measurement) leading to a bright bias in bare-ice albedo associated with residual patches of snow. Assessment of the temporal variation in surface conditions spanning the melt season snow-to-ice transition was therefore needed to determine an optimal threshold. For this reason, we computed the average bare-ice albedo from all station years at the date of ice-ablation onset for a range of ice-ablation thresholds. The lower threshold boundary, set to 4 cm, corresponds to the accuracy of the pressure signal in the ice-ablation setup for each measurement. On the other hand, 8 cm which is the sum of the uncertainties of two pressure measurements, represents a maximized theoretical uncertainty of the ice-ablation www.geusbulletin.org measurement. We added 10% of this value to obtain a conservative upper limit threshold of 9 cm. We found a curve inflection in the resulting albedo at a threshold of 6 cm with the slope of the linear regression being 11 times steeper within the 4-6 cm range than within the 6-9 cm range (-0.98 and -0.09, respectively). We contend that this change in data behaviour illustrates a snow-to-ice transition. As the melt of remaining heterogeneous patches of snow rapidly decreases the albedo, the underlying ice appears and further darkens at a slower rate. To be certain of bare-ice prevalence, we therefore defined the day of bare-ice onset as the first day after the manually selected transition with a cumulative ice ablation >6 cm. The 6 cm threshold equals the ice-ablation uncertainty determined in Fausto et al. (2016). The resulting start of significant ice ablation occurs 4 ± 3 days after the manually selected inflection determined in the first step. The time series of each climate variable was then combined to build composites, that is, multi station-year averages synchronised to the determined emergence of bare ice. The average albedo value at the date of bare-ice onset, our bare-ice-onset albedo, was finally extracted.

Results
In this section, we present and analyse surface conditions ± 45 days around the onset of ice ablation with an emphasis on the albedo transition from snow to ice.

Air temperature
Eight days prior to average bare-ice onset, composite average air temperature increases on average by 0.18°C per day before reaching the melting point (Fig. 2a). After four days of relatively constant air temperatures under snow-melt conditions (0.21 ± 0.09°C), composite air temperature increases with an average rate of 0.57°C per day until bare-ice onset. Composite air temperature then stabilises at 2.75 ± 0.43°C, which we suggest is associated with the sensible heat sink effect of the isothermal melting surface.

Seasonal snow layer thickness
Composite snow height decreases at a rate of 7 ± 8 mm per day over the first 31 days, then gradually increases to 16 ± 7 mm per day until bare-ice onset (Fig. 2b). The high variability before bare-ice onset is associated with differences in weather conditions between stations.

Ice ablation
Before melt onset, average near-zero values of ice ablation have a standard deviation of 3.4 cm, associated with signal noise. Average rate of ice ablation is 3.5 ± 0.5 cm per day after bare-ice onset (Fig. 2c).

Surface albedo
The daily average albedo composite is stable (0.794 ± 0.008) until c. 15 days prior to bare-ice onset. Then, average albedo declines by -0.008 ± 0.007 per day until 6 days prior to bare-ice onset (Fig. 2d). The albedo decline rate then increases to -0.029 ± 0.009 per day until bare-ice onset. This steeper decline is partly driven by emergence of darker bare-ice patches and snow metamorphism. At the onset of snow melt, wet snow metamorphism (Brun 1989) causes rapid grain growth, resulting in the reduction of near-infrared snow albedo (Wiscombe & Warren 1980;Brun 1989). The composite bare-ice-onset albedo is 0.565 ± 0.109, which is between the recommended values for superimposed ice and clean ice in Cuffey & Paterson (2010). Values remain stable for 8 days before albedo further declines at an average rate of -0.009 ± 0.004 per day for 10 days. Finally, a slight decrease (-0.002 ± 0.005 per day) brings the composite albedo to its minimum daily average over the entire period (0.454 ± 0.140), 36 days after bare-ice onset. The mean difference of 0.111 between bare-ice onset and minimum albedo may be the result of ice algal growth (see Stibal et al. 2017). The composite average albedo subsequently increases due to a temperature decrease as the melt season comes to its end and seasonal snowfall begins to accumulate.
We further examined a Cloud Cover Probability Index (CPI) based on longwave downward irradiance and air temperature (van As 2011) to assess the influence of cloud cover on albedo. We find a daily standard deviation of 0.013 on composite albedo, albedo depending on the CPI threshold (from 0.1 to 0.9, in 0.01 steps). On the day of bare-ice onset, a standard deviation of 0.007 suggests little influence of cloud cover.
The significant variability associated with the composite bare-ice-onset albedo (0.565 ± 0.109) is a result of a combination of data acquisition factors, method accuracies and regional variations. While the measurement accuracies are estimated, we cannot assess the method accuracy or the regional variations in bareice-onset albedo because of a lack of in situ measurements at field-verified dates of bare-ice onset and with a widespread spatial coverage. Thus, the relative share of responsibility of each of these sources in the total variability remains unknown. Nevertheless, in order to investigate the stability of this bare-ice-onset albedo value, we conducted 10 000 simulations where half the station years (53) were excluded randomly. We found a standard deviation of 0.011 on average bare-ice-onset albedo, which demonstrates www.geusbulletin.org a low sensitivity to sample size. We also found a low variability within the range of ice-ablation thresholds from 4 to 9 cm (daily albedo standard deviation of 0.005 for the period of interest), boundary thresholds being associated with bare-ice-onset albedo values of 0.585 and 0.560, respectively. While a relatively large difference in bare-ice-onset albedo (0.02) is obtained for ice-ablation thresholds between only 4 and 6 cm, a very small difference (0.002) is obtained for a larger difference of 3 cm between ice-ablation thresholds of 6 and 9 cm. These results further support the change in data behaviour observed at an ice-ablation threshold of 6 cm presented in section 2.2, where the evolution in bare-ice albedo switches from a quickly decreasing to a slowly decreasing regime. Finally, low correlations between average bare-ice albedo and station elevation, latitude and longitude (-0.07, 0.25 and 0.04, respectively) suggest no spatial dependence of the bare-ice-onset albedo.

Application to spaceborne observations
The bare-ice-onset albedo determined in the previous section can be used as an upper bound for ice albedo in order to monitor the evolution of the Greenland bareice area throughout the melt season using spaceborne observations. Here, we focus on 2018 and 2019 which are low and high melt years, respectively.

Albedo retrieval from Sentinel-3 observations
The EU Copernicus Sentinel-3 satellite Ocean and Land Color Instrument (OLCI) provides 21 spectral bands from 400 nm in visible wavelengths to 1020 nm in the near-infrared, from October 2016 to present. Here, we computed snow albedo (A Snow ) from OLCI observations using a fast atmospheric correction technique (Kokhanovsky et al. 2018(Kokhanovsky et al. , 2020. Because the extremely heterogeneous bare-ice surface conditions violate assumptions in Kokhanovsky's theory, we determined bare-ice albedo from a simple empirical approach. This approach consists of a fit between 4729 hourly PROMICE albedo measurements and the nearest in time and space OLCI Top of Atmosphere (TOA) reflectances spanning three years (2017)(2018)(2019). In order to compute this fit, OLCI TOA radiances were first converted to reflectances (R) normalised with the Solar Zenith Angle (SZA) by: R = OLCI TOA radiance π / (cos(SZA) S 0 ) Where S 0 is the TOA solar irradiance measured aboard Sentinel-3.

www.geusbulletin.org
We then defined the bare-ice albedo (A Ice ) from a fit to the average of four OLCI TOA reflectances: A Ice = α (R 400 nm + R 560 nm + R 865 nm + R 1020 nm ) / 4+ β (2) where α corresponds to the slope of the linear regression between OLCI TOA and PROMICE albedo measurements, and β is its intercept. Computing an orthogonal distance regression, we found respective α and β coefficients of 1.003 and 0.058 associated with a high correlation of R = 0.899 and a standard error of 0.006. We subsequently built daily mosaics for Greenland albedo using A Snow combined with A Ice for A Snow values below the average bare-ice-onset albedo (0.565) determined in this study. We found an average absolute difference of 0.078 between A Ice and A Snow , ±10% around the bare-ice-onset albedo.
Clouds were detected and thereafter masked in Sentinel-3 imagery using the Simple Cloud Detection Algorithm (SCDA) version 2.0 (Metsämäki et al. 2015;Wehrlé & Box 2021). This algorithm consists of up to six tests on Sea and Land Surface Temperature Radiometer (SLSTR) TOA reflectances (550 and 1600 nm) and brightness temperatures (3.7, 11 and 12 μm).
We further applied a temporal filter based on outlier detection modified after  to remove remaining cloud artifacts, which would otherwise introduce abrupt temporal variations in the albedo time series. This processing consists of a 10-day rolling average applied to each pixel time series, for cases where the central value is within ± 15% of the window median.
To validate the albedo retrievals, we compared the snow and ice OLCI-derived albedo with PROMICE ground measurements. We first used the CPI to exclude albedo ground measurements acquired under cloudy conditions that we associated with a CPI >0.3. Despite the 'collocation problem' of the large difference in footprint size between the ground station (c. 2 m × 2 m footprint) and the 1 km × 1 km OLCI pixel (see Ryan et al. 2017), we found a high correlation (R = 0.885, N = 549), a Root-Mean-Square Error (RMSE) of 0.079 and an insignificant mean bias (0.005) between OLCI-derived albedo and ground measurements. Further, we found a mean absolute difference of 0.069 within ± one standard deviation around the bare-ice-onset albedo (0.565 ± 0.109). We expect that some part of the RMSE is attributable to transient errors in the ground observations, such as water droplets or ice on the radiometer domes, specular reflections from the station or equipment shadowing on the surface or radiometer. The assessed albedo RMSE can also increase due to undetected clouds affecting the spaceborne observations.
Finally, 'gapless' daily 1 km OLCI albedo grids were generated by updating pixel values when an area is considered to be cloud free. The last valid pixel value covering a given area then remains until a new valid value is determined. The different processing and filtering steps are documented in a Github repository (Wehrlé & Box 2020b).

Deriving Greenland bare-ice area
The bare-ice-onset albedo (0.565) was used to compute bare-ice area from 'gapless' daily albedo averages in 2018 (low melt year) and 2019 (high melt year). The 0.565 threshold represents a compromise between the 105 station years of the dataset as it is associated with variability in ice types linked to background impurity concentration and resulting albedo. While the PROMICE network provides unprecedented access to crucial observations of conditions at the ice-sheet surface, the apparent limitations of in situ monitoring at very high resolution still prevent a precise study of these spatial variations using ground measurements. The standard deviation of 0.109 on the bare-ice-onset albedo is therefore associated with the spatial and temporal variability of the parameters measured at the different AWSs. This standard deviation probably overestimates the true variability of bare-ice-onset albedo across the entire ice sheet. The inevitable localised characteristics of field measurements often result in small sampling sizes relative to the area of the region of interest, which leads to a higher influence of outliers compared to larger datasets. On the other hand, we estimate that sampling size has little influence on the average bare-ice albedo itself (see Section 3). Consequently, we decided not to use this standard deviation as a measure of uncertainty for our thresholding to determine bare-ice area. Nevertheless, to assess the sensitivity of the bare-ice area determination, we applied bare-ice albedo values obtained with theoretical ice-ablation thresholds of 4 and 9 cm (0.585 and 0.560, respectively) introduced in Section 2.

Variations in Greenland bare-ice area
The Sentinel-3 satellite-derived albedo maps upscale our analysis and compare the near record high melt year, 2019 (Tedesco & Fettweis 2020), to the low melt year, 2018, across Greenland (Fig. 3). Daily Greenland ice albedo (including peripheral ice caps) was on average 0.030 (-3.7%) lower in 2019 for the melt season defined as 1 May to 15 September, reaching a maximum difference of 0.053 (-6.7%) on 3 August (Fig. 3a). This maximum difference occurred after a 2019 high-melt event, which reduced the 2019 daily albedo to a minimum of 0.738, lower by 0.044 (-5.6%) than the 2018 minimum of 0.782. As estimated from the Moderate Resolution Imaging Spectroradiometer (MODIS) after , albedo averaged over Greenland land-ice was 0.817 and 0.777 from June to August 2018 and 2019, www.geusbulletin.org respectively. This is 0.013 and 0.006 higher than the values determined in this study for the same time span (0.804 and 0.771).
On average from 1 May to 15 September, daily bareice area was 47 329 km² in 2019 (Fig. 3b), 4.5 times larger than in 2018, and reached a maximum difference of 102 943 km² (3.1 times larger than the previous year) on 3 August (Fig. 3c). This deviation in ratios of average and maximum values between 2019 and 2018 is a result of the early 2019 melt onset, while 2018 bare-ice area remained near zero for the same period. A maximum daily bare-ice area of 153 489 km² occurred in 2019, 73 269 km² or 1.9 times larger than the 2018 maximum (80 220 km²). In 2019, the maximum daily bare-ice area corresponded to 9.0% of the Greenland ice sheet, only 4.7% in 2018. Maximum daily bare-ice area occurred 14 days earlier in 2019 than 2018 (3 and 17 August, respectively), while minimum albedo occurred 9 days earlier in 2019 than 2018 (4 and 13 August, respectively).
Using a watershed algorithm (van der Walt et al. 2014) to determine the outer boundaries of the ice sheet, we identified ice caps as ice bodies separated from the ice sheet. While ice caps represent only c. 3.6% of the ice-sheet surface area, their inclusion increases 2019 and 2018 maximum daily bare-ice areas by 17.2% and 15.5%, respectively, increasing calculated bare-ice area ratios to 10.4% and 5.3%.
Despite the extremely low ice-ablation threshold of 4 cm, corresponding to the uncertainty of a single pressure measurement, we found maximum daily bare-ice areas equivalent to ± 11.4% (2018) and ± 7.0% (2019) of the values determined with the selected threshold of 6 cm. By applying the high ice-ablation threshold of 9 cm, we found differences equivalent to ± 2.0% and ± 2.8% (for 2018 and 2019, respectively) further supporting a relatively low ice-ablation threshold sensitivity.

Conclusions
PROMICE AWS time series of air temperature, snow height, ice ablation and albedo provide insights into snow and meteorological processes at the ablation-driven transition from seasonal snow to bare-ice surface conditions. We identified 6 cm as the first significant measured value of ice ablation through a sensitivity analysis, which matches the ice ablation uncertainty determined by Fausto et al. (2016). By applying the 6 cm threshold to

www.geusbulletin.org
identify the date of bare-ice onset for each station year in a semi-automatic procedure, we determined a bareice-onset albedo of 0.565 ± 0.109 for the Greenland ice sheet. This value is between the recommended values for superimposed ice and clean ice reported by Cuffey & Paterson (2010). After bare-ice onset, we found a further albedo decrease of 0.111, which may be the result of ice-algal growth. Average ice ablation was 3.5 ± 0.5 cm per day while average air temperature remained roughly constant, suggesting that stable air temperatures were associated with the sink of heat energy during surface melting. We found no dependence of bare-ice albedo on elevation, latitude or longitude, suggesting that the bare-ice-onset albedo determined here is representative for locations in between PROMICE AWSs. We further combined snow albedo after Kokhanovsky et al. (2018Kokhanovsky et al. ( , 2020 with bare-ice albedo estimated from a fit between Sentinel-3 OLCI TOA reflectances and PROMICE albedo data. In a cross-validation using 4729 daily PROMICE observations, we confirmed a high correlation coefficient (0.885), a RMSE of 0.079 and an insignificant average bias (0.005) with PROMICE ground measurements.
Applying PROMICE-derived bare-ice-onset albedo to the Sentinel-3 imagery we produced quantitative mapping of albedo and bare-ice area variations. The maximum daily bare-ice area was 1.9 times larger in 2019 than in 2018 (153 489 and 80 220 km², respectively), covering 9.0% and 4.7% of the ice sheet. Peripheral ice caps increase bare-ice area estimates by 17.2% (2019) and 15.5% (2018). Thus, the combination of ground and spaceborne observations yields powerful quantitative constraint on snow-cover dynamics across Greenland ice.