Introduction

The Greenland Ice Sheet has been losing mass since the 1990s and was responsible for ~20% of global sea-level rise during the 2006–2018 period1. This mass loss has largely been driven by an increase in surface melt and runoff caused by higher air temperatures2,3,4. Radiative feedbacks between melt and ice sheet surface albedo make the ice sheet especially sensitive to increasing air temperatures (also known as the “melt-albedo feedback”)5,6,7. As the climate warms, melt increases, reducing albedo, leading to more shortwave energy absorption and more melting. Constraining the strength of this feedback is therefore critical for accurate modeling of current and future meltwater runoff from the Greenland Ice Sheet to the global ocean.

Several processes (e.g., snowline fluctuations, ice algae, snow grain size) are responsible for ice sheet radiative feedbacks, many of which have been extensively studied6,8,9,10. One process that has received less attention is meltwater ponding. During the Northern Hemisphere summer, vast amounts of meltwater are produced on the Greenland Ice Sheet surface11,12,13. Much of this meltwater is temporarily stored on the surface at a variety of scales, from small cylindrical cryoconite holes with diameters of centimeters to large supraglacial lakes with areas of up to ~10 km2 14,15,16,17. Remote sensing and field measurements demonstrate that supraglacial meltwater features have low albedo and melt rates that are 110–170% higher than nearby non-inundated ice sheet surfaces18,19,20. However, despite evidence that meltwater ponding is a key control on albedo and melt in the ablation zone21,22, these local measurements have not been upscaled to the entire ice sheet. The radiative effect of meltwater ponding, and its importance relative to other processes that reduce ice sheet albedo, therefore remains unconstrained.

The increasing availability of medium-resolution satellite imagery now enables meltwater ponding to be mapped accurately at large spatial scales. Several studies have used Landsat and Sentinel-2 to quantify the extent of supraglacial lakes and rivers in specific sectors of the ice sheet at a spatial resolution of 10–30 m23,24,25,26,27,28. More recently, Landsat and Sentinel-2 were used to map supraglacial lakes29,30 and rivers31 across the entire ice sheet. Finally, one study mapped all surface water contained in lakes, rivers, and water-filled crevasses with 10 m Sentinel-2 imagery32. These maps represent a snapshot of surface water coverage for a cool (2018) and a warm (2019) summer and provide a first opportunity to directly investigate the radiative effect of meltwater ponding.

Here, we provide observational constraints on the radiative effect of meltwater ponding using a combination of satellite and drone remote sensing. We first introduce an empirical method to quantify the reduction of surface albedo (Δα) caused by meltwater ponding. To do this, we intersect surface water maps from Sentinel-232 (from Jul 25–30, 2018, and Jul 29–Aug 5, 2019) with coincident MCD43A3 albedo data derived from the MODerate resolution Imaging Spectroradiometer (MODIS) onboard NASA’s Aqua and Terra satellites. Next, we combine Δα with downward shortwave radiation and near-surface air temperatures from MERRA-2 atmospheric reanalysis to quantify the radiative effect of meltwater ponding and its strength relative to other processes that reduce the ice sheet albedo. Finally, we evaluate our findings using high-resolution (30 cm/pixel) drone imagery from 2015 covering almost 300 km2 of the ice sheet. Our drone imagery reveals many narrow streams and small pools that are not resolved by medium-resolution satellite imagery (i.e., Sentinel-2), implying that satellite-based estimates of the radiative effect caused by meltwater ponding should be considered a lower-bound.

Results

Contribution of meltwater ponding to spatial albedo variability

By intersecting satellite-derived surface water maps with surface albedo products (see Methods), we confirm that meltwater ponding considerably influences spatial albedo variability in late-Jul/early-Aug (Fig. 1a–d). The strongest impact occurs between 1000–2000 m a.s.l., where ponding accounts for up to 51% of albedo variability at short length scales (i.e., <3 km; Fig. 2a). Here “length scale” refers to the side length of the search window used for quantifying albedo variability (see Methods). At longer length scales (25 km), ponding still accounts for 24.2% of albedo variability at 1000 m a.s.l., consistent with a previous study reporting 15% in the western ice sheet at the same elevation (length scale also 25 km)16. At lower elevations (0–600 m a.s.l.), ponding accounts for a smaller proportion of albedo variability (8.1% at 25 km), likely because steep slopes at the margins of the ice sheet promote efficient drainage of meltwater (Fig. 2b)33. Above 2000 m a.s.l., the effect of meltwater ponding on albedo variability diminishes (3.4% at 25 km), likely reflecting less surface melt (Fig. 2b). Our findings therefore suggest that meltwater ponding contributes considerably to spatial albedo variability, especially in the upper ablation and lower percolation zones.

Fig. 1: Meltwater ponding is associated with local minima in surface albedo across the Greenland Ice Sheet.
figure 1

a–c show true color imagery from Landsat 8 with the location of transect in red. The lower panels in a–c show same-day albedo (red line) from the MCD43A3 product and fractional water area (blue fill) from ref. 32. a North Greenland on Jul 30, 2019, b Northeast Greenland on Aug 02, 2019, and c Southwest Greenland on Aug 04, 2019.

Full size image
Fig. 2: Radiative effect of meltwater ponding on the Greenland Ice Sheet during the summers of 2018 and 2019 by elevation.
figure 2

a Proportion of spatial albedo variability caused by meltwater ponding during Jul 29–Aug 5, 2019, colored by length scale. “Length scale” represents the side length of the search window used for quantifying albedo variability (see Methods). b Area of meltwater ponding in a warm (2019) and cool (2018) summer. c Radiative effect of meltwater in a warm (2019) and cool (2018) summer. d Radiative effect of meltwater in a warm (2019) and cool (2018) summer, converted to units of energy absorbed by the ice sheet per day. e Regional radiative effect of meltwater ponding in North and Southwest Greenland in 2019. The dashed lines represent the elevations of our field sites in Southwest Greenland that we surveyed by drone. f Sensitivity of the radiative effect caused by meltwater ponding to near-surface air temperatures.

Full size image

Radiative effect of meltwater ponding

Surface albedo is strongly negatively correlated with the area of meltwater ponding (mean r = −0.63) at short length scales (i.e. ≤ 10 km; Fig. S1). The slopes of these local linear regressions between surface water fraction and albedo average −0.11 ± 0.06 (μ ± σ) at length scales ≤ 10 km, meaning that for every 10% increase in surface water fraction, albedo is reduced by 0.011 ± 0.006 (see Methods). We combine this albedo-lowering effect with satellite-derived surface water maps and downward shortwave radiation from MERRA-2 reanalysis to quantify the radiative effect of meltwater ponding (for a detailed discussion of the uncertainties associated with the use of a static albedo-lowering effect, refer to the Assumptions and Uncertainties section in the Methods). We find that the radiative effect of meltwater ponding exhibits large elevational variations. At 1000–2000 m a.s.l., where large amounts of meltwater accumulate on the surface, the radiative effect averaged 0.30 ± 0.16 W m−2 in the summer of 2019 (Fig. 2b, c). At lower (<600 m a.s.l.) and higher (>2000 m a.s.l.) elevations, the radiative effect of meltwater ponding is smaller, averaging 0.13 ± 0.07 W m−2 and 0.002 ± 0.001 W m−2, respectively (Fig. 2b, c). Summed across the ice sheet, we find that meltwater ponding enhances energy absorption by 5.7 ± 3.1 PJ d−1 in Jul 25–30, 2018 to 15.7 ± 8.6 PJ d−1 in Jul 29–Aug 5, 2019 (Fig. 2d). For context, spatiotemporal albedo variability enhanced ice sheet energy absorption by 2176 PJ d−1 in the summer of 20196. Therefore, our satellite-based analysis indicates that the radiative effect caused by meltwater ponding accounted for 0.7 ± 0.4% of the total ice sheet surface shortwave radiative effect. We note that this value represents mid- to late-summer hydrological conditions and that the contribution of meltwater ponding to total shortwave radiative effect may be higher earlier in the summer.

We find that the radiative effect of meltwater ponding is strongest within the percolation zone34,35. In North Greenland, for example, the radiative effect of meltwater ponding is higher at 1000–1200 m a.s.l. (1.0 ± 0.5 W m−2) than at lower elevations (e.g., 0.1 ± 0.06 W m−2 at 200–400 m a.s.l.; Fig. 2e). Likewise, in Southwest Greenland, the radiative effect of meltwater ponding is highest at 1600–1800 m a.s.l. (1.9 ± 1.0 W m−2; Fig. 2e). This is over six times higher than the radiative effect of meltwater ponding in the dark zone, which is located at 1000–1400 m a.s.l. (0.3 ± 0.2 W m−2; Figs. 1c, 2e). It is plausible that the prevalence of meltwater ponding in the percolation zone may be explained by recent evidence of ice layers or ‘slabs’36,37,38. When infiltrated meltwater refreezes in porous snow or firn, impermeable ice layers form just below the surface that can promote surface runoff into supraglacial streams and lakes31,39. Processes that determine meltwater drainage and infiltration therefore appear to be the primary control on spatial patterns of the radiative effect caused by meltwater ponding.

Interannual variability in the radiative effect of meltwater ponding appears to be largely driven by air temperatures. Satellite-derived surface water maps reveal that more of the ice sheet surface is covered by meltwater during warm summers32 (Fig. 2b, f). Reference 32 found that meltwater covered 4903 km2 in 2018 (a cool summer) and 9988 km2 in 2019 (a warm summer). Given that ice sheet summer (JJA) near-surface (2 m) air temperatures were 1.76 K higher in 2019 than 2018, we estimate that meltwater ponding increased by +1715 km2 K−1, implying that the radiative effect of meltwater ponding is sensitive to air temperatures. We find that the slope of this relationship is +5.7 ± 5.2 PJ d−1 K−1, although we note that the uncertainty is relatively large (Fig. 2f; see Methods). We also note that this relationship may be more complex than identified by our analysis. For example, a melting ice sheet surface will have a larger cooling effect on near-surface air temperatures than a surface that is below freezing. Likewise, the sensitivity of meltwater ponding to air temperatures likely exhibits spatiotemporal variability. For example, the formation of extensive near-surface semi-impermeable ice layers by refrozen meltwater reduces firn meltwater storage capacity, as observed following the extreme 2012 melt season37,39. It is therefore possible that meltwater ponding is more sensitive to air temperatures immediately following extreme melt events, rather than before38. Further mapping of surface water before and after 2012 would be required to investigate such non-linear relationships between air temperatures and meltwater ponding.

Meltwater ponding in high-resolution drone imagery

The accuracy of our findings largely depends on the quality of the satellite-derived water products we used to characterize the area of meltwater ponding. To investigate the accuracy of these products, we surveyed two sites in the ice sheet’s ablation zone at high spatial resolution (0.3 m) using a series of drone surveys in Jul 2015 (Fig. 3a, b). The first site includes both Russell Glacier and Isunguata Sermia (~180 km2) with an elevational range of 150–660 m a.s.l. (Fig. 3c). The second site ( ~ 110 km2) is situated in the dark zone with an elevational range of 1170–1290 m a.s.l. (Fig. 3d). At Russell Glacier and Isunguata Sermia, we find that only 0.001% of the surface was ponded on Jul 11–12, 2015 (Fig. 4a, b; see Methods). This is consistent with ref. 32. who did not identify any supraglacial lakes, ponds, or streams in the same sector using Sentinel-2 imagery in Jul 25–30, 2018 and Jul 29–Aug 5, 2019. Our high-resolution surface water maps therefore confirm that there is negligible meltwater ponding at the margins of the ice sheet.

Fig. 3: Drone surveys across the ablation zone of the Southwest Greenland Ice Sheet provide high-resolution imagery for investigating meltwater ponding.
figure 3

a Locations of drone mapping areas and insets shown in Fig. 4. Background is a Sentinel-2 true color image from Aug 02, 2019. b Location of inset a on the Greenland Ice Sheet. c Survey patterns over Isunguata Sermia and Russell Glacier at the ice sheet margin (150–660 m a.s.l.) on Jul 11–12, 2015. d Survey patterns over the upper ablation zone (1170–1290 m a.s.l.) on Jul 21, 2015.

Full size image
Fig. 4: Examples of meltwater ponding in the Southwest Greenland Ice Sheet ablation zone imaged by drone.
figure 4

Locations of insets can be found in Fig. 3a. a–c There is negligible meltwater ponding (<0.1%) at the margins of the ice sheet due to crevasses and steep slopes. d–e In contrast, the upper ablation zone (1170–1290 m a.s.l.) on Jul 21, 2015 is characterized by vast numbers of small ponds, streams and larger supraglacial lakes and rivers that, in total, cover 3.6% of the imaged area (110 km2, see Fig. 3a). f Tents at our field site, surrounded by small streams and ponds.

Full size image

In contrast, our high-resolution surface water maps reveal substantial meltwater ponding in the ablation zone at elevations higher than our ice-sheet margin field site (Fig. 4d–f). At our dark zone field site, 3.6% of the ice sheet surface was inundated on Jul 21, 2015 – three to four times more than previously estimated using Sentinel-2 imagery in Jul 25–30, 2018 (0.9%) and Jul 29–Aug 5, 2019 (1.2%)32. Some of these differences may be attributed to interannual variability. However, cumulative runoff was similar between Jul 21, 2015 and Jul 25–30, 2018 at the time of mapping, and substantially more runoff occurred before the Jul 29–Aug 5, 2019 mapping (Fig. S3), rendering our conclusion conservative.

A more plausible explanation for the differences in surface water coverage between the drone and satellite products is improved spatial resolution. With a spatial resolution of 0.3 m, we document thousands of meltwater ponds and streams <100 m2 in our drone imagery (Fig. 5). Collectively, these features account for 53.7% of all surface water at our field site. In contrast, the satellite-derived surface water maps (with a spatial resolution of 10 m) resolve none of these <100 m2 ponds32. These maps also appear to underestimate the number of ponds <1000 m2 (Fig. 5). The satellite-derived water maps therefore substantially underestimate the total area of meltwater ponding implying that, in the dark zone, its radiative effect may be stronger than previously thought.

Fig. 5: Histogram of surface meltwater area distributions from drone and satellite imagery at our dark zone field site in Southwest Greenland.
figure 5

Dashed lines represent 100 m2 and 1000 m2. We find that meltwater features (ponds and streams) <1000 m2 account for 63.6% of the total meltwater area while features <100 m2 account for 53.7%. Satellite-derived products underestimate meltwater features <1000 m2 and do not resolve features <100 m2, leading to underestimation of the radiative effect caused by meltwater ponding.

Full size image

Our drone imagery reveals that the radiative effect of meltwater ponding was 0.9 ± 0.5 W m−2 at our field site in the dark zone on Jul 21, 2015. For context, the radiative effect caused by glacier ice albedo change has been previously estimated to be 19.1 W m−2 during the same period6, indicating that 4.7 ± 2.6% of the radiative effect caused by reductions in glacier ice albedo can be attributed to meltwater ponding. This is substantially larger than that predicted by satellite-derived maps of surface water (i.e., 0.2 ± 0.1 W m−2 [1.0 ± 0.5%] in Jul 25–30, 2018 or 0.3 ± 0.2 W m−2 [1.7 ± 0.9%] in Jul 29–Aug 5, 2019). Furthermore, we note that the radiative effect of meltwater ponding is weak in Southwest Greenland at 1250 m a.s.l. relative to the lower percolation zone (1600–1800 m a.s.l.; Fig. 2e). If satellite-derived maps underestimate meltwater ponding by similar magnitudes in the lower percolation zone (i.e., by a factor of four), then the radiative effect caused by meltwater ponding would account for up to ~9% of the total surface shortwave radiative effect in this area of the ice sheet (which is 38.2 W m−2). High-resolution drone imagery in the percolation zone would be required to evaluate this inference. Even without such evidence, our findings demonstrate that it is likely that the radiative effect of meltwater ponding represents a considerable source of energy for ice sheet melt in the upper ablation and lower percolation zones.

Discussion

The prevalence of supraglacial lakes and rivers across the surface of the Greenland Ice Sheet led early studies to hypothesize that meltwater ponding exerts a primary control on albedo patterns22,40 and meltwater runoff21. However, these hypotheses were never directly evaluated outside of a 25 km transect across the western ablation zone16. The strength of the radiative effect caused by meltwater ponding therefore remained unconstrained and excluded from models that forecast the Greenland Ice Sheet’s contribution to sea-level rise. Here, we evaluated the meltwater-albedo hypothesis across the entire ice sheet using satellite-derived maps of surface water coverage and albedo. Our findings demonstrate that (1) meltwater ponding has an important local impact on spatial albedo patterns, (2) the radiative effect of meltwater ponding is strongest in the percolation zone and weaker at the ice sheet margins, (3) in the summer of 2019, meltwater ponding was responsible for at least 0.7 ± 0.4% (15.7 ± 8.6 PJ d−1) of total surface shortwave radiative effect, and (4) the radiative effect of meltwater ponding is sensitive to interannual variation in summer near-surface air temperatures, increasing at a rate of at least +5.7 ± 5.2 PJ d−1 K−1. Given the close relationship between energy absorbed by the ice sheet and meltwater runoff (Fig. S4), and the potential for more meltwater ponding in a warmer climate, we argue that this process should be incorporated into ice sheet and climate models if they are to accurately forecast ice sheet mass loss.

Our satellite-based findings are conservative because high-resolution drone imagery indicates that medium-resolution satellite imagery underestimates the area of meltwater ponding. Although medium-resolution satellite imagery is able to map meltwater features >1000 m2 fairly accurately, we find that ponds and streams smaller than this threshold are substantially underestimated at our field site in the dark zone of the ice sheet (1170–1290 m a.s.l.; Fig. 5). This is problematic because small meltwater features (<1000 m2) account for 63.6% of the total ponded area. The radiative effect of meltwater ponding may therefore be three to four times higher than predicted by satellite-derived maps in the dark zone. High-resolution (<0.5 m) imagery would be required to assess the magnitude of this bias in other areas of the ice sheet. Until then, our satellite-based findings should be considered a lower-bound on the radiative effect of meltwater ponding.

In this study, we focused on meltwater ponding, but note that liquid water also reduces ice sheet surface albedo in other ways. For example, the presence of meltwater can lower snow albedo by increasing effective snow grain size and filling interstitial air-spaces between snow grains41. Snow albedo is therefore strongly dependent on liquid-water content. Suppressed albedo between 20–130 km of the transect in Fig. 1a provides qualitative evidence for this process occurring in a region that has been termed the “slush zone”42. Recently, studies have demonstrated that the slush zone can grow to a much larger area than other hydrological features (e.g., supraglacial lakes) in North and Southwest Greenland during high melt years (e.g., 2019)27,28. While full characterization of the relationship between slush and surface albedo is beyond the scope of our study, we anticipate that its radiative effect could be substantial. For example, at 1600–1800 m a.s.l. in Southwest Greenland, mean albedo from the MCD43A3 product averages 0.62 during Jul 29–Aug 5, 2019, well below that of dry snow (~0.84). Measurements and modeling reveal that light-absorbing particulate (LAP) concentrations are too small to substantially reduce albedo in this sector of the ice sheet43. Likewise, our findings indicate that the area of meltwater ponding is too small to reduce albedo by this magnitude. This implies that most of the albedo reduction in the Southwest Greenland percolation zone during Jul 29–Aug 5, 2019 can primarily be attributed to slush formation and snow grain size growth.

Our findings have implications for modeling the Greenland Ice Sheet’s contribution to sea-level rise. Currently, coupled ice sheet-climate models do not explicitly accumulate meltwater on the surface, so they do not account for shortwave energy absorption caused by meltwater ponding. The exclusion of the radiative effect caused by meltwater ponding may therefore cause these models to underestimate current and future rates of meltwater production. Despite that, incorporating the radiative effect of meltwater ponding into climate models should be feasible since glacio-hydrological models have been available for some time18,44,45. These models use surface elevation data to route and pond meltwater over the Greenland Ice Sheet and would be useful for including the radiative effect of meltwater ponding into projections of ice sheet mass loss. Furthermore, this coupling would allow ice sheet models to retain meltwater in lakes that can refreeze at the beginning of winter46 or evaporate during summer, which would further improve the realism of these models and the accuracy of their projections.

Overall, our study provides a set of observational constraints on the radiative effect of meltwater ponding using a combination of satellite and drone imagery. Our satellite-based analysis describes the spatiotemporal variability of this process at the ice-sheet scale. However, it likely represents a lower-bound on the radiative effect of meltwater ponding since medium-resolution satellite imagery does not resolve the full area of meltwater ponding, much of which is contained in narrow streams and small ponds (Fig. 4d–f). Our drone-based analysis, on the other hand, enables accurate mapping of meltwater features of all sizes. It therefore provides a local upper-bound on the radiative effect of meltwater ponding (+0.9 ± 0.5 W m−2) at 1170–1290 m a.s.l. in Southwest Greenland on Jul 21, 2015. Given that there is more meltwater ponding at elevations higher than our field site (Fig. 2b, e) and that 2015 was a relatively cool summer (Fig. S3), it is likely that the radiative effect of meltwater ponding is even stronger in other areas of the ice sheet (i.e. the lower percolation zone) and in other years (i.e., 2012 or 2019). Coupled climate-ice sheet models should therefore consider including this process for more accurate forecasts of meltwater runoff from the Greenland Ice Sheet.

Methods

Data

We obtained binary surface water coverage maps from ref. 32. which applied a semi-automated water extraction algorithm to 271 Sentinel-2 images to map surface water across the entire Greenland Ice Sheet. These maps represent two distinct time periods: Jul 25–30, 2018, and Jul 29–Aug 5, 2019 when most of the imagery was acquired. We primarily used the 2019 map for our analysis. ref. 32 categorize surface water into three groups: supraglacial lakes, supraglacial rivers, and water-filled crevasses. After comparison with drone imagery from a tidewater glacier, we found that this product substantially overestimates the area of water-filled crevasses, so we decided to exclude water-filled crevasses from our analysis (Fig. S5, S6). Our ice sheet-wide analysis therefore focuses on the radiative effect of meltwater in supraglacial lakes, ponds, rivers, and streams.

We obtained surface albedo data for the same time periods from the MCD43A3 version 6.1 product. This albedo product provides black- and white-sky albedo on a daily basis using a 16-day retrieval window based on MODIS data from Terra and Aqua47. We used black-sky albedo following previous studies that have demonstrated that there is little difference between black- and white-sky albedo for typical solar zenith angles at local noon during summer48,49. We resampled both the surface water and surface albedo data onto the ISMIP6 grid by averaging values within each target grid cell. The ISMIP6 grid has an NSIDC Sea Ice Polar Stereographic North (EPSG:3413) projection and a grid cell resolution of 1 × 1 km50.

We retrieved downward all-sky shortwave radiation data for the periods Jul 25–30, 2018, and Jul 29–Aug 5, 2019, from the MERRA-2 radiation diagnostics product (M2T1NXRAD). This product is provided at an hourly temporal resolution so we averaged these data to daily means. We also statistically downscaled M2T1NXRAD from 0.625 × 0.5° resolution to the ISMIP6 1 × 1 km grid following an elevation-based approach described by previous studies6,51. For each MERRA-2 grid cell, we performed a local linear regression between elevation and downward longwave (or shortwave) radiation using the 24 nearest grid cells (~100 km search radius). We then applied the regression coefficients to higher spatial resolution elevation data (1 × 1 km). If the correlation between elevation and downward longwave (or shortwave) radiation was not significant (p  <  0.1), we did not apply a correction.

We acquired high-resolution aerial imagery using a fixed-wing drone similar to that described by ref. 52. The drone was equipped with a 16.1 MP Sony NEX-5N digital camera and programmed to capture RGB imagery in RAW mode every 50 m of horizontal distance traveled. Each image was georeferenced using the onboard L1 GPS along with attitude (pitch, roll, yaw) data recorded by the accelerometer. The drone collected imagery at two separate field sites. The first site includes both Russell Glacier and Isunguata Sermia (180 km2) with an elevational range of 150–660 m a.s.l53. During six surveys, the drone collected 9732 overlapping images of the glacier surface on Jul 11–12, 2015. The second site (110 km2) is situated in the dark zone with an elevational range of 1170–1290 m a.s.l. Here, the drone collected a total of 3795 images during three surveys between Jul 20–22, 2015. We used Agisoft Metashape Pro v2.2.0 to generate orthomosaics from the overlapping aerial imagery. The final products have a spatial resolution of 0.30 m, offering unprecedented capability for quantifying the area of meltwater ponding on the ice sheet surface.

Local albedo variability and the radiative effect of meltwater ponding

We quantified the impact of meltwater ponding on albedo variability using an empirical approach based on the satellite-derived maps of surface water and albedo. To do this, we compute the standard deviation of albedo within a series of search windows with lengths varying from 1 to 100 km (Fig. S7). In each search window, we compute the standard deviation of albedo for all grid cells (σall) and the standard deviation of albedo for grid cells that do not contain any surface water (σdry). Albedo variability (\({{{{\rm{\sigma }}}}}_{{{{\rm{albedo}}}}}\)) caused by meltwater ponding can then be defined as

$${{{{\rm{\sigma }}}}}_{{{{\rm{albedo}}}}}=1-\frac{{{{{\rm{\sigma }}}}}_{{{{\rm{dry}}}}}}{{{{{\rm{\sigma }}}}}_{{{{\rm{all}}}}}}$$
(1)

This metric is best understood by considering a window in which there is little meltwater ponding. In this window, σdry will be similar to σall. Hence, σalbedo will be low. Alternatively, when there is a lot of meltwater ponding, σall will be larger than σdry, signifying that more of the albedo variability can be attributed to meltwater ponding. We compute this metric for 32,000 randomly selected grid cells that are evenly spread across sixteen elevation bins ranging from 0–3,400 m a.s.l. with an interval of 200 m (i.e., 2000 grid cells per elevation bin). This approach allows us to investigate how σalbedo varies with both length scale and elevation.

We derived the albedo-lowering effect of meltwater ponding using a similar approach. However, instead of computing the standard deviation of albedo in each search window, we linearly regressed albedo (α) against surface water fraction (fwater):

$${{{\rm{\alpha }}}}={{{{\rm{\beta }}}}}_{0}+{{{{\rm{\beta }}}}}_{1}\cdot {{{{\rm{f}}}}}_{{{{\rm{water}}}}}$$
(2)

The slope of this local linear regression (β1) represents the response of albedo to the area of meltwater ponding. We applied this approach to a random sample of grid cells in our coincident maps of surface water and albedo from Jul 25–30, 2018, and Jul 29–Aug 5, 2019 (n = 9079). We find high correlation coefficients between meltwater ponding and surface albedo at length scales < = 10 km (mean r = −0.63), indicating that spatial albedo patterns can primarily be attributed to meltwater ponding at these length scales (Fig. S1 and S2). At greater length scales (i.e., >10 km), other ice sheet properties, such as snow grain size, become more important for spatial albedo patterns, and the impact of meltwater ponding becomes weaker. For example, the correlation coefficients (r) between meltwater ponding and surface albedo decline to −0.56 ± 0.21 at 25 km, −0.45 ± 0.22 at 50 km, and −0.36 ± 0.20 at 100 km. We therefore define the albedo-lowering effect of meltwater ponding as the mean sensitivity of surface albedo to meltwater ponding (β1) at length scales < = 10 km (Fig. S2). The total albedo change caused by meltwater ponding (Δα) for any grid cell can subsequently be computed as

$$\Delta {{{\rm{\alpha }}}}={{{{\rm{\beta }}}}}_{1}\cdot {{{{\rm{f}}}}}_{{{{\rm{water}}}}}$$
(3)

We quantified the radiative effect of meltwater ponding (\({{{\rm{F}}}}\)) for each grid cell as

$${{{\rm{F}}}}=\left(1-\,\Delta {{{\rm{\alpha }}}}\,\right)\cdot {{{\rm{S}}}}{{{{\rm{W}}}}}_{{{{\boldsymbol{\downarrow }}}}}$$
(4)

where SW↓ is downward all-sky shortwave radiation from M2T1NXRAD and \({{{\rm{F}}}}\) is in units of W m−2. Next, we computed the sensitivity of the radiative effect caused by meltwater ponding (\({{{\rm{\lambda }}}}\)) from

$${{{\rm{\lambda }}}}=\frac{\Delta {{{\rm{F}}}}}{\Delta {{{\rm{T}}}}}$$
(5)

where \(\Delta {{{\rm{F}}}}\) is the difference in the radiative effect between the summers of 2018 and 2019, \(\Delta {{{\rm{T}}}}\) is the difference in near-surface air temperature between the summers of 2018 and 2019. \({{{\rm{\lambda }}}}\) is in units of W m−2 K−1. Since we were interested in the radiative effect at the surface, we used 2 m air temperatures from MERRA-2 (M2T1NXSLV). We computed the sensitivity of the radiative effect to air temperature in 200 m elevation bins ranging from 0 to 3400 m a.s.l. and converted \({{{\rm{\lambda }}}}\) into units of energy (in J d−1 K−1) to account for variations in area within each elevation bin.

Meltwater ponding from drone orthomosaics

We mapped meltwater ponding in the orthomosaic produced over the dark zone (1170–1290 m a.s.l.) using a semi-automated classification approach. To do this, we first generated a validation dataset by manually digitizing polygons representing water and non-water. Next, we segmented the orthomosaic using Simple Non-Iterative Clustering (SNIC) to remove noise in the imagery54. We also computed a Normalized Difference Water Index (NDWI) from our orthomosaic using the red and blue bands. The median NDWI value of each segment was then assigned to all pixels within that segment. To find the optimal threshold for distinguishing between water and non-water segments, we iterated through a range of NDWI values (from 0 to 0.3 with a 0.05 interval). The optimal threshold was defined as the threshold with the highest combination of accuracy and precision relative to our validation dataset (Fig. S8). By explicitly accounting for precision (which improves as the number of false positives decreases) during the optimization process, we ensured that our classification does not overestimate meltwater ponding. Once we identified this threshold (0.175), we applied it to every pixel in the segmented orthomosaic to classify water vs. non-water automatically. This approach classified meltwater ponding with an accuracy of 91.5% relative to our validation dataset (Fig. S8).

Our semi-automated approach performed poorly for the Russell Glacier and Isunguata orthomosaic (150–600 m a.s.l.) because self-shadowed pixels caused by crevassing were characterized by high NDWI values. Therefore, for these study sites, we manually digitized ponded areas using the Geo-SAM QGIS plugin, an interactive segmentation tool based on the Segment Anything Model (SAM) foundation AI model55,56. It is possible that the use of different methods to map surface water could introduce some uncertainties. However, based on visual inspection of the imagery and given that both approaches segment water pixels using NDWI values, it is unlikely that these uncertainties are responsible for the large differences in water coverage between the dark zone and margin field sites. We consider the high-resolution (0.3 m/pixel) maps of surface water from the drone imagery to represent “ground-truth” for comparisons with the satellite-derived surface water maps.

Assumptions and uncertainties

In this study, we assume that the albedo-lowering effect of meltwater ponding (β1) is static and can be applied across the ice sheet. This decision is justified by low correlation coefficients between the β1 and latitude (R2 = 0.002) and elevation (R2 = 0.008). Likewise, we find small correlation coefficients between β1 and ponded area (R2 = 0.07), which can be considered a proxy for pond depth. ref. 29, for example, found that pond area is significantly correlated with depth with an R2 = 0.83. We speculate that the lack of correlation between β1 and pond area (depth) is because most of the albedo reduction occurs in the first 20–30 cm of the melt pond18,57. In contrast, albedo exhibits little variation with pond depths >30 cm. ref. 32 find that the mean depth of meltwater features in their dataset 0.66 ± 0.10 m (μ ± σ). There is therefore unlikely to be much variation in β1 with depth for most of the meltwater features in our study, justifying our use of a static value (Fig. S9).

We also assume that meltwater ponding is the only process responsible for the albedo-lowering effects we observe in our local linear regression analysis. However, it is possible that other processes are also responsible for the difference in albedo between inundated and non-inundated ice sheet surfaces. For example, if the concentrations of LAPs are consistently higher in supraglacial lakes than in surrounding areas, then the direct effect of meltwater ponding on albedo will be overestimated. Instead, meltwater ponding would also have an indirect effect on albedo by enhancing the concentration of LAPs. While we are not aware of data to systematically evaluate this assumption, drone imagery indicates that the ice sheet surfaces surrounding supraglacial lakes can be both darker or brighter than the bottoms of supraglacial lakes (Fig. S10; see also Fig. 5c in ref. 16). Therefore, while we cannot dismiss indirect effects, we assume that inundated and non-inundated ice sheet surfaces have similar LAP concentrations and that most of the observed albedo-lowering is attributable to meltwater ponding.