Temporal Variability of Surface Reflectance Supersedes Spatial Resolution in Defining Greenland’s Bare-Ice Albedo

Ice surface albedo is a primary modulator of melt and runoff, yet our understanding of how reflectance varies over time across the Greenland Ice Sheet remains poor. This is due to a disconnect between point or transect scale albedo sampling and the coarser spatial, spectral and/or temporal resolutions of available satellite products. Here, we present time-series of bare-ice surface reflectance data that span a range of length scales, from the 500 m for Moderate Resolution Imaging Spectrometer’s MOD10A1 product, to 10 m for Sentinel-2 imagery, 0.1 m spot measurements from ground-based field spectrometry, and 2.5 cm from uncrewed aerial drone imagery. Our results reveal broad similarities in seasonal patterns in bare-ice reflectance, but further analysis identifies short-term dynamics in reflectance distribution that are unique to each dataset. Using these distributions, we demonstrate that areal mean reflectance is the primary control on local ablation rates, and that the spatial distribution of specific ice types and impurities is secondary. Given the rapid changes in mean reflectance observed in the datasets presented, we propose that albedo parameterizations can be improved by (i) quantitative assessment of the representativeness of time-averaged reflectance data products, and, (ii) using temporally-resolved functions to describe the variability in impurity distribution at daily time-scales. We conclude that the regional melt model performance may not be optimally improved by increased spatial resolution and the incorporation of sub-pixel heterogeneity, but instead, should focus on the temporal dynamics of bare-ice albedo.


Introduction
Melt and runoff from the Greenland Ice Sheet has been greater this century than during the previous 350 years [1]. In recent decades, the ice sheet has experienced intensified and prolonged melt seasons [2], more anticyclonic air circulation and reduced cloud cover [3,4], a greater proportion of precipitation falling as rain [5,6], and expansion of the ice sheet's ablation area [7][8][9]. Each have contributed to record runoff, ice sheet mass loss, and global sea-level rise [10]. The bare-ice ablation area contributes up to 80% of total runoff [11], with the majority of melt driven by shortwave radiation fluxes [12,13]. As a result, Greenland's surface melt and mass balance is highly sensitive to seasonal albedo variations [14,15]. Recent satellite-based studies have reported long-term darkening of the ablation zone [16][17][18], underpinned by spatially variable inter-and intra-annual evolution of bare-ice reflectance [14,19,20]. However, the detail of short-term spatiotemporal reflectance dynamics in Greenland's bare-ice regions remains poorly resolved and crudely represented in regional climate models. This, therefore, represents an important research goal for accurate projection of the ice sheet's surface mass balance and contribution to sea-level rise [21,22].
Greenland's seasonally bare-ice zone comprises a wide variety of ice types and foliation, interspersed by crevasses, fractures, residual snow patches, meltwater features, and scattered distributed impurities (including dust, aerosols, emergent particles and debris, mineral-microbe aggregates, algae, and other microbiota) [23][24][25][26][27][28][29][30]. Drivers of change in the optical properties of Greenland's bare-ice have previously been investigated at seasonal timescales and illustrate the importance of reductions in superimposed ice [19,20,23] and the growth of ice algae [25,27,31]. These processes darken the bare-ice surface over the summer, while autumn freeze-up by the creation of superimposed ice [32,33] and freezing of cryoconite holes [34], along with new snowfall, raises the surface reflectance. Reflectance variability is often attributed to ongoing melting across a complex and spatially heterogeneous bare-ice zone where several inter-dependent melt processes may operate simultaneously [23,25,35,36]. These include: melt-driven creation of a porous weathering crust, which may then be affected by subsurface freezing or removed by rainfall or high turbulent energy exchange [26,37]; transient surface meltwater ponding in the weathering crust or crevasses [30,[37][38][39]; evolving surface meltwater channel networks [40]; and, the distribution and geometry of impurities and cryoconite holes in response to relative fluxes of radiative and turbulent energy [41][42][43]. Collectively, these concurrent, spatiotemporally diverse processes, driven by micro-to synoptic-scale meteorological conditions, result in a quasi-stochastic variability in the reflectance of the ice sheet's bare-ice surface.
Recent studies have attempted to investigate the relative impact of these processes on bare-ice reflectance variability at the satellite sub-pixel scale [25,26,30,31,35,[43][44][45][46]. These studies show that the presence of impurities causes a positive skew in the distribution of surface reflectance over 250-500 m Moderate Resolution Imaging Spectrometer (MODIS) pixel scales [26,47], and begin to demonstrate a complex and dynamic interplay between the bio-physical characteristics that underpin the spatiotemporal variability of reflectance. The presence of these variabilities at very fine spatiotemporal scales are further complicated by both the evolution of surface topography [48,49] and its associated anisotropic characteristics [50]. The subsequent retrieval of bare-ice albedo is therefore particularly challenging due to the requirement of appropriate bidirectional reflectance distribution functions (BRDFs) for the many and varied surface types in the bare-ice zone [51,52].
In the face of these challenges, regional surface mass balance models have applied a range of reflectance values across the entire Greenland Ice Sheet: from invariant reflectance [53][54][55], to spatially variable but temporally static values derived from satellite observations [56], to composite satellite products with coarse spatial resolution (~1 km) and temporal (8-16 day) variability [57]. A recent assessment using the MODIS MOD10A1 product in the Canadian Arctic Archipelago [58] demonstrated that, while accounting for spatial variation, employing a temporally interpolated average albedo yielded improved model performance relative to MODIS-derived minimum [59] or seasonally-averaged [60] bare-ice albedo. However, the 500 m resolution MOD10A1 product fails to capture spatial surface albedo variability [26], and there is now an increasing focus on developing higher spatial resolution ice sheet surface mass balance models [22,56,61]. There is, therefore, an outstanding requirement to explore high resolution/fine-scale heterogeneity in bare-ice reflectance to resolve the impact of spatiotemporal variations in the optical properties of Greenland's bare-ice over daily to seasonal timescales.
To address this research gap, we investigate short term (diurnal-to-synoptic) and seasonal variability in bare-ice surface reflectance using multi-scale observations from a site in western Greenland, obtained during the 2016 ablation season. Our data, captured over an area of~0.13 km 2 (or 13 ha), and at spatial resolutions from 0.01 to 10 m, are used to examine variations in the distribution and mean surface reflectance over a bare-ice region. We use a simplified energy balance model to illustrate the importance of adequately describing the spatial and temporal variability in local areal mean reflectance or albedo for modelling bare-ice ablation.

Study Site
The study site is located on the Kangerlussuaq (K-) Transect in western Greenland, 38 km from the ice margin and proximal to the Institute for Marine and Atmospheric Research Utrecht (IMAU) S6 Automatic Weather Station (AWS) (67 • 04.78 N, 49 • 24.08 W, 1020 m a.s.l., Figure 1a). This location is within the Greenland's ablation area where 90% of the net June-August (JJA) energy fluxes are defined by the surface albedo and an annual ice loss of~1.96 m w.e. is typical [12]. Measurements were collected during the ablation season between 16 July and 17 August 2016 (day of year (DOY) 198-230) during which time the site was characterized by bare-ice [20]. The site lies at the western edge of the 'dark zone', where dust, cryoconite, and algae reduce the surface albedo [25,27,29,30,36,44]. A supraglacial region of interest (ROI) was demarcated over~0.13 km 2 and ablation measurements were taken following standard protocols [62] at daily to 3-day intervals at 12 regularly distributed stakes.
Local meteorological data, including shortwave incoming (SW in ) and reflected (SW out ) radiation and air temperature (T a ), were recorded at the IMAU S6 AWS, located at a similar elevation approximately~1 km west of the ROI ( Figure 1a) and exhibiting equivalent, representative bare-ice surface characteristics (see Smeets et al., 2018; Figure 1a). Time periods influenced by cloud were identified from the incident radiation record. Ice surface reflectance at S6 was derived following Smeets et al. [63] using a 24-h moving 'accumulated albedo' ratio of |SW out |: SW in .

Satellite Reflectance Products
To examine how bare-ice albedo evolved over the study period, we acquired the 500 m MOD10A1 (Collection 6) albedo product [64] for the MODIS pixel containing the centerpoint of ROI for the summer melt season months (30 May-7 September: DOY 150-250) 2000 and 2016. The MOD10A1 product provided the black-sky broadband albedo estimate for the visual and near infrared (VNIR) wavelengths. These observations were taken under clear-sky conditions and therefore the black-sky estimate, assuming no diffuse illumination, was considered the most appropriate. These MOD10A1 data were filtered for quality issues using the Snow Cover Basic QA band bitmask and, following the 11-day running median and standard deviation methods proposed by Box et al. [14], we constrained albedo values to between 0.21 [23] and 0.84 [14]. The accuracy of the albedo product over bare-ice is thought to be <0.07 [46,65]. For comparison, the daily MOD09GQ (at 250 m pixel size) and the 8-day composite MOD09A1 reflectance products were also acquired for the pixel containing the ROI during the 2016 study period. To enable comparison with medium resolution data, Sentinel-2 multi-spectral imagery (at 10 m horizontal resolution) was acquired for the same MODIS pixel area. A total of fifteen Sentinel-2 cloud-free scenes were acquired for the JJA ablation period in 2016 and were processed to surface reflectance (including topographic normalization) using the 6S atmospheric model [66], implemented using the open-source Python library ARCSI [67]. The 6S model was parameterized with an aerosol optical depth (AOD) of 0.05 using the 'maritime' aerosol profile, the atmospheric water vapor and ozone were defined using the 6S standard 'subarctic summer' atmospheric profile. The mean VNIR reflectance was calculated using the 10 m resolution Sentinel-2 VNIR bands for the analysis and comparison with our other measurements. We opt to explore 'reflectance' as a proxy for albedo, owing to the multi-scale nature of our study and the absence of the appropriate BRDFs for heterogenous bare-ice across the scale ranges.

Field Spectroscopy Observations
Field-spectroscopy was conducted using a StellarNet Blue-Wave (SNBW) spectrometer, measuring reflectance over visible wavelengths (400-700 nm) within which >50% of the total incident solar radiation is received, and where most of the variability in ice surface reflectance takes place [68]. The spectrometer was fitted with a 10 • mat black collimating cowl giving an effective measurement footprint diameter of 0.1 m (or measurement area of 0.008 m 2 ) from the acquisition height of~0.7 m. Following a consistent operator approach, hemispheric-conical reflectance factor (HCRF) measurements [42,69] were collected in natural illumination under constant cloud conditions within two hours of local solar noon at a spectral resolution of 1 nm. The sensor was positioned at nadir and oriented into the principal solar plane. For each survey, HCRF data were recorded at 30 points distributed over the ROI (Figure 1b) in a grid layout at~50 m spacing to provide a spatially representative sample of surface reflectance [46].
At each of the 30 points, individual sample surfaces were identified in a circle of radius 1 m around the central location, sampling at~45 compass degrees. Each sample measurement was a mean of 10 replicates made with the sensor remaining in position and with a consistent integration time (<1 s). The reflectance of a horizontal white Halon reference panel (>97% reflectance) was made immediately before and after each sample measurement set to estimate the incident irradiance, with all samples at a point measured within a~2-min period. This sample measurement sequence was repeated ten times across the ROI between DOY 198 to 230. As the survey points were guided by hand-held GPS locations with an estimated horizontal precision of 3 m, and the ice flows to the west at 70 m a −1 , measurements were not made repeatedly at identical locations throughout the survey period but remained within the ROI. Post-processing of the SNBW data involved the discarding of corrupt sample spectra (totaling <7% of all sample measurements), the transformation of individual spectra to irradiance using the instrumental calibration, and calculation of reflectance defined as the ratio of each ice surface sample irradiance to the mean measurement set of the white reference panel irradiance.

Uncrewed Aerial Vehicle Imagery and Reflectance Measurements
Upscaling from discrete spectroscopic measurements, yielding data at a finer spatial resolution, was achieved using a DJI Phantom Pro 3 drone (hereafter, DPP3) following the uncrewed aerial vehicle (UAV) methodologies described by Ryan et al. [30,35,70]. Between DOY 203 and 230, the DPP3 was flown periodically over the study site at 70 m altitude above the ice surface (Figure 1b), using a pre-defined flight-plan. Optical red-green-blue (RGB) images (12 Mpix) were taken with the onboard FC300X camera with a 4 mm focal length and ISO-100. For each survey flight, ice surface characteristics across the ROI were reconstructed using Agisoft PhotoScan Pro (v1.2). Photogrammetric point clouds contained between 180 and 250 million points for each scene. Following mesh construction, orthomosaics and point clouds were generated at ground resolutions of between 2.5 and 2.7 cm pixel −1 . These datasets were re-projected (using nearest-neighbor resampling) to 2.5 cm pixel resolution using fixed markers on the ice surface demarcating the ROI, to enable the comparison of imagery acquired on different survey dates. Geolocation of the images included feature mapping and co-registration to ensure the same surface area extent was included in the final orthomosaics, and locational artefacts arising from ice flow were avoided. In the absence of ground control point data for each survey, horizontal uncertainties of~1 m in the co-located orthomosaics were estimated from feature matching and accuracy of DPP3 s onboard GPS.
To produce comparable datasets over JJA, here, the photogrammetrically derived RGB orthomosaics of the ROI were transformed to surface reflectance through an empirical line calibration [71,72] using in-situ reflectance measurements of matt white, grey and black calibration panels collected using the StellarNet spectrometer. This post-processing normalization approach assumed illumination across each image set (collected over a < 18-min window) remained constant, and was employed to account for irradiance variations between the DPP3 flights, and internal image color adjustment during the photogrammetric orthomosaic processing itself. The calibration panels (each~0.1 m 2 or 16 pix per panel) were placed within the ROI, and nadir reflectance at each panel was measured coinciding with the time of each DPP3 survey. Spectral resampling was undertaken to ensure correspondence between the spectrometer readings and the RGB data recorded by the DPP3: spectra were resampled to RGB band discretization by fitting Gaussian functions to the onboard camera sensor responses and integrating these at 1 nm intervals to enable the total reflectance over each RGB band to be estimated from the spectral radiometer readings [73]. Empirical linear relationships for each waveband were derived between the spectrometer reflectance measurements and the associated orthomosaic Digital Numbers (DNs) from the three calibration panels by squaring the original DN values (yielding r 2 > 0.992). This resulted in surface reflectance for each of the RGB wavebands, and while linear methods exist for the narrow-to-broadband albedo conversion [74], these have not been verified for the DPP3 low altitude imaging and on-board camera. Therefore, here, as a proxy for surface reflectance across the visible spectrum, we used the mean of the RGB band values for each pixel [75].

Numerical Modelling of Ice Ablation
To explore the impact of dynamism in surface reflectance on ice melt (ablation), we employed a distributed melt model [76] based on a simplified surface energy balance [77]. Ablation is defined by the total energy available for melt (Q m ), given by: where α is local albedo or reflectance, SW in is incident shortwave radiation, T a is air temperature at 2 m, and c 1 and c 2 are constants. The combined components c 1 T a + c 2 estimate the turbulent heat fluxes and longwave radiation. The model was driven by the hourly meteorological input data: SW in and T a collected at the IMAU S6 weather station, and Q m converted to ice ablation, in m w.e. The in-situ ice ablation measurements across the ROI, at daily time-steps, were used to internally calibrate c 1 and c 2 . Spatial sensitivity was explored using a 600 × 600 cell model domain that included an arbitrary topography taken from an excerpt of the photogrammetrically derived resampled point cloud for 23 July (DOY 205). During the reference model run (which was based on the ice reflectance distribution obtained from the DPP3 data on DOY 205) we calibrated our model and optimized the bulk exchange coefficient to account for the sensible heat flux. This coefficient was kept constant in all consecutive model runs where synthetic reflectance distributions were implemented over the model domain to describe albedo (α) for the ROI. Then, to assess the relevance of temporal sensitivity to areal mean reflectance, we re-optimized the model for a single point for differing albedo (α) parameterizations over the summer ablation season (JJA). Eleven α parameterizations from previous studies were used (Table 1) to simulate temporal variations in the areal mean surface albedo. Of the five temporally invariant parameterizations, we used: a fixed bare-ice mean [56]; MODIS minimum α [59]; and seasonal means derived from filtered MOD10A1 data [60] and following a linear interpolation to fill missing daily values [58]. Rather than average the lowest 5% of values extracted from a 16-day average composite product (MOD43 or MCD43A3) over a multi-year period [60], we opted for the lowest 10% of values in the MOD10A1 record (here, n = 5) for 2016 to define a constant α min5 .
Continuing to use the MOD10A1 data series to provide a temporally coarse variability in ice surface characteristics, we firstly simulated an 8-day composite albedo series (α MOD8dc ) by extracting the values for the days that provided the MOD09A1 composite data points (cf. [57]). Examination of the MOD10A1 records for the ROI from 2000 to 2015 highlighted a perennial, broadly sinusoidal pattern of albedo across the summer season, as reported by Box et al. [14]. In recognition of the cosine approach used in early positive degree-day models of ice sheet surface mass balance [78], we fitted a cosine function centered on mid-July to the MOD10A1 albedo records for 150 days from May to September (MJJAS: DOY121-271); the cosine function was limited by α < 0.84 [14,79] and defined by day-of-year (time, t). The ordinary least squares coefficient of determination (r 2 ) was used to iteratively solve values for the amplitude of variation around the mean (α amp ) and constants d 1 and d 2 , where d 1 − d 2 represents the day of the onset of albedo decline as the seasonal snow begins to metamorphose and melt. For each of the years from 2000 to 2015, this cosine approach revealed a mean r 2 of 0.68 (0.39 < r 2 < 0.89) between the observed and parameterized albedo, with the spring onset of the albedo decline ranging from DOY 83 to 132. To represent more physically based bare-ice albedo parameterizations [39,80], we redefined our cosine model, exchanging the independent variable, time (t), for cumulative Positive Degree Days (ΣPDD), where PDD was defined as the daily average temperature if >0 • C. As a comparative alternative, we also examined a parabolic fit for α against ΣPDD, where constant b 2 is the ΣPDD at which α reaches its minimum (b 3 ) ( Table 1). Table 1. Description of constant and time-variant albedo parameterizations explored for prediction of cumulative melt over the ROI, including the parameter values and Nash-Sutcliffe, Willmott, mean absolute and root mean square error goodness-of-fit metrics (respectively, NSE, d, MAE and RMSE) to observed ablation over the JJA season. The associated goodness-of-fit metrics for the average hourly melt rate between on-site ablation observations are shown in parentheses. NSE and d are dimensionless, MAE and RMSE are in m w.e. Goodness-of-fit metrics are displayed with a red relative heatmap shading, ranging from red for weaker to white for stronger model performance.  For each temporally varied albedo parameterization, we calibrated the ablation model coefficients c 1 and c 2 using the observed average ablation measurements (in m w.e. hr −1 ) defined over the time-steps between ablation surveys. These coefficients were then kept constant over time, and SW in and T a records for 2016 were used to drive the ice ablation model simulations.

Results
We focus here on the nature and impact of spatiotemporal variability in reflectance characteristics, describing reflectance distribution statistics within the ROI over time. Although our analysis precludes the direct comparison of absolute values due to differences in sensor configuration, after co-registration of the UAV orthomosaics with a co-located systematic ground-sampling strategy within the ROI, our approach ensures comparability and consistency in the spatiotemporal patterns detected at each scale.

Bare-Ice Reflectance Variability
The 2016 melt season was darker than average [20], with MOD10A1 albedo persistently below the mean albedo recorded for the locality between 2000 and 2015 ( Figure 2a). The long-term accumulated albedo proxy recorded at S6 suggested broadly stable surface reflectance typified the locality during the summer period: reflectance averaged~0.47 (± 0.06) with intermittent, temporary decreases or increases (Figure 2c). These temporal variations qualitatively showed an association with the synoptic changes evidenced by periods of cloud cover which decreased SW in (e.g., DOY 199-202, 215, 222-224 and 229; Figure 2b). Bare ice persisted until DOY249, when snowfall increased reflectance at S6 to >0.7. Notwithstanding the uncertainty in representativeness of point measurements [46], the discrepancy between the HCRF reflectance approximations and true albedo [42,52,81], and the locational difference between our ROI and S6, the MOD10A1 albedo record for the 2016 bare-ice season (DOY180-250) correlates with the equivalent daily mean S6 reflectance time-series (r = 0.48, p < 0.001).
Significant correlations between the 10 pair-wise combinations of the MOD10A1, S6 reflectance, MOD09GQ and Sentinel-2 time-series (0.33 < r < 0.93, p < 0.01) highlighted the statistical similarity in the seasonal (JJA) pattern of local albedo and reflectance at the locality (Figure 2c-e). This overarching pattern was also observed in the ground-based SNBW and DPP3-derived reflectance data (Figure 2e), with means of 0.385 and 0.468, respectively. The cause of the lower SNBW average reflectance value is likely due to the sample locations for the SNBW measurements missing the brightest areas of the ice surface (see Figure 1b). Owing to the absence of concurrent measurements and a dearth of data pairs, retrieving correlations for these fine-scale observations with the satellite-derived time-series is not possible. Cloud cover can influence observed glacier surface albedo by up to 0.08 [23]; this likely results from the combined influence of factors, including the spectral attenuation of irradiance [82] on the heterogenous surface, and the changing optical properties of emergent, ablating ice [42]. During cloudy periods, typically the ratio of shortwave irradiance to turbulent energy fluxes decreases and this degrades and/or removes the more reflective, less consolidated and porous decimeter-scale surface ice layer that is created by the subsurface melting that is particularly associated with clear-sky conditions [26,83]. However, the DPP3 reflectance data collected on DOY224 and SNBW measurements on DOY 200 and 223 were consistent with similar records from clear-sky days. This suggests that, owing to the brevity and nature of the cloudy phases during our observation period, the effect of such influences was slight. With both reflectance datasets being generated from nadir downward-facing ratios and calibrated against the same white reference, these data are sufficiently representative. As the spatial resolution of observations increases, the range of reflectance values also rises (Figure 2f). The coarser resolution of the Sentinel-2 product thus masks the finer resolution spatial variability through pixel scale spatial averaging. At the finest resolution of the DPP3 product, discrete and localized points of high or low reflectance are recorded. This variability is demonstrated in standard deviations of reflectance (σ; Figure 2g), particularly in the latter half of the period. Despite a high range, DPP3 shows a lowered σ, indicative of reduced variability around the mean, with a high degree of spatial autocorrelation [46]. The SNBW surveys show a smaller number of observations than the DPP3 spatial maps, and sample points at distances that are larger than the spatial autocorrelation, and so reveal more variability around the mean.
Changes in distribution between each measurement set are most clearly illustrated by the skewness (s) and kurtosis (k) (Figure 2h,i). The Sentinel-2, SNBW and DPP3 datasets all exhibit positive (right) skew, as previously reported for bare-ice using MODIS and Sentinel-2 data [26,45,47]. While the day-to-day variability in the distribution function for Sentinel-2 is masked as a result of the temporal resolution and the spatial averaging at the pixel scale, the higher spatial resolutions show that non-Gaussian distributions are persistent and that these distributions vary at the daily timescale. Such variability is suggestive of short-term changes and dynamics in the ice surface properties and characteristics.
Temporal variability is investigated further by examining spatial reflectance frequency distributions themselves (Figure 3a-c, refer also to Figure 1b). The spread of values depends on the dataset resolution: for example, Sentinel-2 shows narrower reflectance distributions compared to SNBW and DPP3, a result of spatial averaging within the pixel footprint, while the higher spatial resolution datasets better detail the range of reflectance seen over the heterogeneous bare-ice. With features, such as cryoconite holes, occupying~8% of the surface area in the locality [84], the coarser resolution of Sentinel-2 eliminates observation of these discrete, 'dark' (<0.2) reflectance points and reduces the presence of tails within the distribution, with values exhibiting a closer fit around the mean. In contrast, the SNBW and DPP3 show a wider range, with spatial resolution able to capture the larger spread of reflectance values and, consequently, a lower mean. Over time, as shown by the distribution descriptors (Figure 2), the variations in the histogram form were not systematic. The change in distribution between each observation, for the three contrasting scale reflectance datasets, revealed changes of up to 20% (Figure 3d-f). This is indicative of a rapidly evolving and topographically complex bare-ice environment.
Our multi-scale observations support the previously identified need to examine subpixel variability in bare-ice surface reflectance in west Greenland [26,85], and to explore the relative precision of reflectance and/or albedo derived from platforms, particularly satellite platforms, at a range of scales. The discrepancies between in-situ HCRF, Sentinel-2 and MODIS products, as shown here, are unresolved yet well-known and often ascribed to differences in spatial aggregation and calibration of conversion methods [26,45,46,85]. While current methods often simplify the conversion of reflectance to albedo, at fine resolutions uncertainties in BRDFs, especially for lower reflectance surfaces common to heterogenous bare-ice environments typical of the ice sheet's ablation zone, result in errors in albedo evaluation of up to 11% [52].
Our analysis shows that despite differences in sensor configuration, there are variations in the spatiotemporal patterns (Figure 3) of surface reflectance that can be resolved at each scale of measurement. This is critical, as it undermines the use of near-surface measurements, such as field spectroscopy and fixed-albedometers, as providing perfect benchmark satellite-validation tools [35]. Moreover, it also precludes the use of coarseresolution satellite imagery to fully elucidate the processes operating across the ice surface that cause albedo to vary.

Spatial Variability
The ice melt (ablation) model [76] was used to examine the relevance of the spatial reflectance variability observed in this study. To begin, model runs employed a Gaussian reflectance distribution for variable mean albedo (Figure 4a), randomly assigned over the model domain (Figure 4b). Although linear conversions are feasible [74,86], in the absence of validated BRDF models at the necessary spatial resolution for impurity contaminated glacier ice, we use reflectance as a proxy for albedo (α). Comparison between model runs utilizing the arbitrary topography and a flat, unvaried surface yielded changes in the total melt of <0.5%; this result indicates that small-scale topography (slope and aspect variations) makes minimal difference to seasonal melt production. The mean reflectance varied while retaining a constant standard deviation (σ), skewness (s) and kurtosis (k), and revealed a linear response of surface melt to changes in albedo: a 0.01 change in bare-ice albedo increases or decreases melt by 2% (Figure 4e). As expected, retaining a constant mean (0.468, the average derived from the eight DPP3 surveys), and sequentially adjusting σ, s and k independently, resulted in negligible changes (<1%) in the melt production across the model runs, likely to be associated with the randomized population.  In (a,c), distributions with the lowest and highest mean α are shown with solid and dashed vertical lines, respectively; in (b,d), albedo is in greyscale, from black (α = 0) to white (α = 1). (e) Illustrates the linear increase in melt following reduction in mean surface reflectance for the normal distribution (circles) and for the beta distributions (triangles) focused on the window 0.35 < α < 0.6. The filled circles and triangles in (e) represent the reflectance distributions and associated means shown in (a,c).
As we have shown, variability in the reflectance distribution across the ROI exhibits simultaneous changes to the mean, σ, s and k (Figure 2). To simulate these changes, a beta distribution was employed to describe α over the model domain. The standard beta probability density function, f (x), defined over the interval 0 to 1, is parameterized by two positive shape parameters (A and B > 0) and is well-suited to modelling the dynamic distributions of observations of a percentage variable, X (e.g., albedo measurements over space) [87]: where x is a realization of X and β is the beta function. Here, again using the ice reflectance distribution obtained from the arbitrarily selected DPP3 survey on DOY205 as the 'control', we fitted a probability density function to the observed data, which defined parameters A and B as 2.108 and 2.583, respectively, (or A:B = 0.82; Figure 4c,d). As the beta distribution can also be defined by the sum of the shape parameters, we generated synthetic surface reflectance distributions for the model domain defined by A + B = 4.7; this provided positively and negatively skewed distributions with co-varied mean, σ, s and k (Figure 4c). The output from this second set of model runs highlighted a statistically similar inverse relationship between the mean reflectance and total melt (Figure 4e): 95% confidence intervals for the linear regressions showed significant similarity in the slope of the albedomelt relationships. The mean daily variation of 0.025 exhibited in the MOD10A1 data series for 2016 implies that, for ablation controlled by albedo, daily changes in meltwater production were ±4.75%. These simple simulations show that subtle changes in distribution of reflectance revealed by our nested examination of sub-MODIS-pixel variability in 2016 are not primary controls on melt production, at least at the scale of our ROI. It remains unclear as to whether coarser-resolution products (derived from MODIS or Sentinel-2) adequately and robustly represent pixel area mean reflectance through the assumed linear mixing models. Rather, our data suggest that for a given area, short-term (daily) changes in the areal mean reflectance may be an important control on melt production.

Temporal Variability
To assess the importance of temporal variation in areal mean reflectance of the bare-ice zone of western Greenland, we explored the efficacy of simulations of cumulative ablation against observed ablation (stake measurements) over the ROI by reporting the mean absolute error (MAE), the root mean squared error (RMSE) and the Nash-Sutcliffe Efficiency criterion (NSE: [88]) and Willmott Index of Agreement (d: [89]) which are predictive skill metrics commonly used in hydraulic modelling (Table 1). These metrics were also calculated for the average hourly melt rate each model predicted compared to that for the observed ablation measurements. Over the JJA period, characterized by clear skies and daily peak irradiance that was typically >500 Wm −2 , the daily average air temperature remained above 0 • C throughout the summer (except for DOY182, 228, 229) until the autumn transition following DOY240 (Figure 5a). The variety of albedo parameterizations (Figure 5b) all performed similarly well in simulating cumulative ablation with goodness-of-fit metrics NSE and d > 0.9, and typical mean absolute errors of <0.10 m. The projections of cumulative bare-ice ablation of between 2.01 and 2.45 m w.e. (Figure 5c) fit well with the local long-term annual mean surface mass balance of −1.96 m w.e. at the IMAU S6 weather station [12], particularly in light of 2016 being a relatively moderate, unexceptional melt year [8].
The models that performed least well in predicting cumulative ablation across the goodness-of-fit metrics were the invariant fixed (α fix ), minimum (α min ) and mean observed MOD10A1 (α MODmean ) albedo parameterizations (Table 1); in contrast, all the models employing a dynamic albedo input yielded higher NSE and d indices, and lower MAE and RMSE (Table 1). Intermediate model performance was found for the linearly interpolated daily MOD10A1 series (α MODint ), the lowest 10% MOD10A1 values (α min5 ), and the composite product approximation (α MOD8dc ). At the average hourly melt time-scale, despite modelled ablation using α fix resulting in high NSE and low MAE scores, across the model skill metrics, the dynamic albedo models consistently outperformed other formulations, and particularly the α min5 and α MOD8dc parameterizations. All the model runs exhibited an underestimation of melt between DOY199 and 217 (Figure 5d), when reflectance values were low (Figure 2). With seasonal minimum albedo typically reached~36 days after bare-ice onset [90], we speculate that this underestimation of melt is a function of ice algae: discrete algal blooms not only darken the ice surface over length-scales below that of the resolution of satellite observation, but also exhibit non-linear growth following a brief lag phase after the onset of bare-ice conditions [25][26][27]44]. These biological early melt season processes may not be fully captured by either satellite products or air temperature-derived estimations of albedo. Following DOY218, all of the model variants fall within the range of observed ablation measurements, which also illustrate the cumulative uncertainties associated with the ground-based measurements themselves. Examination of our modelled ablation during the observation period highlights how the temporally dynamic α simulations performed almost identically (Figure 5d): all yielding RMSEs < 0.013 m w.e. and mirroring the α MODint output. These latter models all tracked the ROI's median ablation well, and show that the use of poorly defined, invariant albedo values can result in over-or under-estimation of cumulative summer ablation by ±11%.

Discussion
Our multi-scale analysis reveals a spatial and temporal heterogeneity of bare-ice reflectance in the western section of the Greenland Ice Sheet. The daily MODIS MOD10A1 albedo record for the pixel containing our study site during the 2016 summer JJA melt season exhibited a decreasing trend followed by an increase, over a range of 0.26, superimposed with day-to-day variation up to 0.09. Even excluding the end of summer period from DOY240, the seasonal bare-ice albedo variation of 0.16 exceeds the 0.11 reported for bare-ice across the ice sheet [90]; late season increases in surface reflectance is not solely a result of late season snowfall accumulation [90], but rather can also arise independently from the formation of ice lids and/or freeze-up of cryoconite holes [34] and initial autumn superimposed ice formation [13,33]. The daily variability, however, is comparable to or below the uncertainties reported from assessments of the MOD10A1 product of between 0.03 to 0.1 [14,45,46,91], with the greater magnitudes of uncertainty likely in the JJA period owing to bare-ice heterogeneity [35]. The MOD10A1 product failed to resolve albedo increase [26], despite Sentinel-2 and fine-scale multi-spectral measurements showing substantial increases in clean ice over the pixel scale. Consequently, confidence in the short-term variations in the MOD10A1 record may be questioned. Nonetheless, similar seasonal and marked daily-to-synoptic scale variations in the daily and composite MOD09 reflectance products and the weekly Sentinel-2 observations for the bare-ice period during the melt season are observed.
Mean summer reflectance in the 0.25 km 2 ROI study site from Sentinel-2 showed the same seasonal cycle as that observed in the MOD10A1 time series, albeit over a smaller reflectance range (~0.13). However, Sentinel-2 s coarse temporal resolution precludes identification of shorter-term variability and surface reflectance value tend towards the higher end of in-situ measurements (e.g., our DPP3 and SNBW data, see Figure 3), and darker and brighter surfaces and over-and under-estimated, respectively [92]. This is due to areal averaging of sub-pixel (<10 m) scale heterogeneities. Spatial aggregation of mean surface characteristics at both the MODIS (250 or 500 m) and Sentinel-2 (10 m) pixel scale is enhanced due to some impurities (e.g., cryoconite holes that are less than 1 m in diameter; Figure 1b) being obscured from non-zenith illumination and non-nadir observation angles [35,45,46] and spatial changes occurring at sub-pixel length scales [25,26].
Despite showing similar, relatively small (<0.13) daily to multi-day changes in mean reflectance, as observed by Sentinel-2, our finer scale data displayed increased standard deviation (σ > 0.1) and range (>0.95) in the reflectance distribution. Such resolution-dependent distribution has been reported elsewhere [52]. Cloud cover is known to adjust retrieved reflectance derived from ground-and UAV-based measurement platforms owing to changes in the downwelling and upwelling spectral energy distributions [85,93]. However, only three days of our fine-resolution observations were influenced by cloud cover, and deviations in reflectance distribution were not particularly marked; changes in spatial mean reflectance were characteristically less than the maximum of 0.08 identified for cloud-effects on HCRF broadband reflectance for Greenland ice types [23]. These findings suggest that processes causing variations in the apparent optical properties of bare-ice at daily timescales under clear sky conditions may have a similar magnitude impact on surface albedo as changes in cloud cover. Therefore, we argue that the reflectance records presented here demonstrate a seasonal and daily variability in bare-ice reflectance distribution and, by inference, a spatiotemporal heterogeneity in the albedo of the bare-ice surface, as implied by the MOD10A1 data record.
The processes underlying the dynamics in the optical characteristics of the bare-ice environment of our ROI are expressed in not only the mean surface reflectance but its changing distributions. Such seasonal and synoptic transitions in reflectance and albedo distributions have previously been reported in western Greenland [26,47]. At the start of the bare-ice phase, higher reflectance (>0.6) surfaces, such as residual seasonal snow and superimposed ice, result in either bimodal or negatively skewed unimodal reflectance distributions [45,47]. However, our fine-resolution data does not capture these early season conditions, and despite the elevated mean reflectance (DOY183-193) shown by Sentinel-2, its spatial resolution and sampling bias [92] likely masks the true expression of the reflectance distributions. Once early summer season bright surfaces are eliminated, with increased exposure of darker bare-ice, and associated impurities, meltwater and ice algal growth, the mean reflectance decreases and reflectance distributions tend to be those with a positive skew with variable kurtosis, as evident in our multi-scale datasets. Subsequently, as ice ablation is reduced later in the season, mean reflectance increases and skewness declines with a clear trend toward negative skew ( Figure 2h); we acknowledge, however, that due to the synoptic event evident for DOY236-243 and the duration of our fine-resolution observations, our datasets do not fully capture the freeze-up of the bare-ice surface.
At the scale of a MODIS pixel, skew in reflectance distribution may result in a 2% increase in hourly surface melt predictions [26]. Based on our model of ablation utilizing a beta distribution to simulate a skewed surface albedo distribution, we demonstrate that the nature of the albedo distribution skew is not itself critical for simulating melt at daily to seasonal timescales if the areal mean surface albedo is accurately characterized. Recognition of the temporal evolution of the areal mean albedo of the bare-ice surface is, instead, crucial. Approaches utilizing temporally (and spatially) constant representations of bare-ice albedo will be subject to greater uncertainties: arbitrary bare-ice albedo or an albedo defined by the lowest observed value(s) for a location may, respectively, under-or over-predict seasonal ablation by~10%. Certainly, for Greenland's dark zone and other marginal areas, assumptions that bare-ice albedo lies between 0.5 and 0.55 [54], which is closer to the bare-ice onset of 0.565 [90], will underestimate ice melt. Similarly, use of satellite-derived albedo minima will overestimate melt through lacking recognition of the seasonal changes in the optical properties of the ice surface.
Our ablation model sensitivity analysis study shows how the use of an albedo parameterization that portrays short-term (daily) temporal variability, rather than fixed values, results in improved ablation estimates. Here, ablation model performance at seasonal and daily (or average hourly) time-scales was improved whether accounting for fine scale temporal variability by using a mean albedo calculated from linearly interpolated daily time-series of the MOD10A1 record [58] or by employing more sophisticated periodically evolving [57], time-dependent, or physically based parameterizations. Consequently, our results support the use of melt-dependent parameterizations [39], as employed within the Modèle Atmosphérique Régionale (MAR) [18,94]; although in its current form, with α limited to >0.4, MAR likely underestimates ablation, particularly in the darker regions of the Greenland Ice Sheet. We therefore argue that exploration, development and testing of cyclical or physically based bare-ice albedo dynamics at moderate to high temporal resolution are required for further forecasting and sensitivity studies.
Irrespective of the data source uncertainties in observing true bare-ice albedo, this study emphasizes the need to account for the short-term spatial variation in the optical properties of the bare-ice surface in modelling ice ablation. Specifically, there is a need to constrain the representativeness and accuracy of the spatial mean albedo at numerical model grid resolutions. While dynamic non-Gaussian reflectance distributions, primarily arising from ice surface processes, are themselves not an important consideration for ablation, at least at daily timescales, definition of representative spatial mean reflectance observations, and variability therein, is essential for accurate predictions of bare-ice melt at daily to seasonal timescales. To achieve this, it is critical that the inherent surface variability induced by the biologically active, ablating ice interface is recognized. With some regions in Greenland recently experiencing increased clear-sky conditions [3,4], high-pressure blocking [95,96], rising air temperatures [2] and increasing bare-ice areas [8,9], adequately characterizing bare-ice albedo, and its spatiotemporal evolution, is critical to improving projections of the ice sheet's surface mass balance.

Conclusions
Despite its critical role in amplifying Greenland's ice ablation, bare-ice albedo remains poorly characterized and parameterized in numerical models. Here, we demonstrate substantial spatial and temporal variability in the reflectance of the bare-ice environment using a multi-scale analysis that uniquely draws together and compares finer-resolution ground-based spectrometry, UAV-and Sentinel-2 derived reflectance observations and the coarse resolution MODIS reflectance and albedo products. Our analysis highlights daily-to-seasonal synoptic variability in bare-ice reflectance, evident in both the spatially averaged reflectance for our 1300 m 2 study area along with the actual reflectance distributions measured in-situ with a field-spectrometer. At the highest spatial resolution, varying but consistently positive skewed reflectance distributions were observed throughout the 2016 peak melt season. Modelling bare-ice melt with these reflectance distributions described using a beta function yields no significant difference compared to a Gaussian distribution, but highlights the importance of accurately describing the areal mean albedo. Further assessment exploring the temporal variability in bare-ice albedo emphasizes the critical importance of parameterizations at daily-to-synoptic resolution, with inappropriate schemes yielding large errors in seasonal ablation of up to 11%. Given that summer high-pressure, clear-sky blocking conditions have been increasing in the last 30 years, with concomitant air temperature warming and expansion of ice sheet ablation areas, we argue that future research efforts must provide robust albedo evaluations for heterogenous bare-ice across the range of available observational platform scales, and further refine evolving bare-ice albedo parameterizations using temporally and physically based descriptors.

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