Comparison of Multiple Methods for Supraglacial Melt-Lake Volume Estimation in Western Greenland During the 2021 Summer Melt Season

Abstract

Supraglacial melt-lakes form and evolve along the western edge of the Greenland Ice Sheet and have proven to play a significant role in ice sheet surface hydrology and mass balance. Prior methods to quantify melt-lake volume have relied upon Landsat-8 optical imagery, available at 30 m spatial resolution but with temporal resolution limited by satellite overpass times and cloud cover. We propose two novel methods to quantify the volume of meltwater stored in these lakes, including a high-resolution surface DEM (ArcticDEM) and an ablation model using daily averaged automated weather station data. We compare our methods to the depth-reflectance method for five supraglacial melt-lakes during the 2021 summer melt season. We find agreement between the depth-reflectance and DEM lake infilling methods, within +/−15% for most cases, but our ablation model underproduces by 0.5–2 orders of magnitude the volumetric melt needed to match our other methods, and with a significant lag in meltwater onset for routing into the lake basin. Further information regarding energy balance parameters, including insolation and liquid precipitation amounts, is needed for adequate ablation modelling. Despite the differences in melt-lake volume estimates, our approach in combining remote sensing and meteorological methods provides a framework for analysis of seasonal melt-lake evolution at significantly higher spatial and temporal scales, to understand the drivers of meltwater production and its influence on the spatial distribution and extent of meltwater volume stored on the ice sheet surface.

1. Introduction

During the summer months (May–August), near-surface air temperature often exceeds 0 °C, melting accumulated snow and firn to produce surface meltwater in the Greenland Ice Sheet (GrIS) ablation zone. The meltwater is routed through a complex surface network of saturated crevasses [1] and supraglacial streams [2,3,4,5] and routed into topographic depressions forming supraglacial melt-lakes (SGLs) [2]. The primary drivers of the production of surface meltwater are incoming shortwave radiation and surface air temperature [6,7,8]. Supraglacial melt-lakes that form along the margins of the Greenland Ice Sheet have been observed to form in the same coordinate location across multiple years [2,9,10,11,12], and their locations are strongly modulated by basal topography [13,14]. Surface meltwater stored in melt-lakes has been observed to quickly (<2 h) transfer a significant amount (~44 × 10⁶ m³) of meltwater deep into the ice sheet [15] in discrete drainage events. The rapid injection of vast amounts of meltwater into the ice sheet, known as hydrofracture events, can decouple the ice sheet from the bedrock below, reducing basal friction, and significantly impact glacier velocity [16,17,18,19,20], with Maier et al. (2023) identifying drainage events taking place during the winter months [20]. The frequent injections of meltwater into the ice sheet’s base contributes to the observed ice sheet mass loss, such as the 327 ± 40 Gt mass loss from Sermeq Kujalleq alone between 1972 and 2018, making it vital to understand melt-lake evolution. Furthermore, Krawczynski et al. (2009) [13] identified lakes with a minimum area (0.049 km²) and volume (0.098 × 10⁶ m³) to reach a threshold for hydrofracture. Across the ablation zone of western Greenland, numerous melt-lakes that reappear annually are found to meet the theoretical threshold for hydrofracture in multiple melt seasons [21]. Because meltwater production on the GrIS is heavily modulated by surface air temperature and insolation receipt, the climatological warming trend (0.23 °C per decade) along the western GrIS over recent decades [22], along with increases in the presence of high pressure in the region (known as Greenland Blocking) during the summer months [23], supports the observed increase in the count and areal coverage of melt-lakes across the GrIS [8,24]. Future climate scenarios suggest significant warming across the GrIS [25], which would augment the production of surface meltwater, yielding a higher number of supraglacial melt-lakes that have the propensity to hydrofracture, that is, drain catastrophically, in turn having a profound impact on ice sheet dynamics and the ice sheet’s contribution to sea level rise. Fang et al. (2023) identified a linear relationship between surface temperature and ice sheet mass balance, such that a 1 °C increase in surface air temperature would result in a nearly 100 Gt/year decrease in ice sheet mass balance, nearly all of it (94%) stemming from glacier surface processes, such as drainage of supraglacial melt-lakes [26].

Remote sensing methodologies have been developed to determine melt-lake area and volume amounts in the ablation zone, relying on sensors, such as the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [27,28,29], the Moderate Resolution Imaging Spectroradiometer (MODIS) [9,10,11,12,30,31,32,33], and Landsat sensors [5,9,10,27,34]. Estimation of supraglacial melt-lake area is relatively straightforward as it requires observation of the melt-lake surface, which can be acquired from cloud-free satellite imagery. Melt-lake volume, however, is more challenging as it requires information on the lake’s bathymetry to generate depth values and thus produce melt-lake volume. The remoteness and harsh climate of the ablation region of Greenland make collecting in situ observations on melt-lake depth extremely challenging, which has resulted in few available datasets on melt-lake depth; these datasets acquire depth information on a small number of melt-lakes and at one instance in time [11,29,35]. In order to expand beyond the few instances of direct depth measurement, a significant effort has gone into coregistering depth values to at-satellite reflectivity, thus developing depth-reflectance algorithms [11,29,36,37,38,39]. Most recently, Pope et al. (2016) [40] applied a depth-reflectance schematic for several bands aboard multiple Landsat satellites and found that two operational land imager (OLI) bands on the Landsat-8 satellite—red and panchromatic—are optimal for SGL depth retrieval. Melling et al. (2024) [41] applied a combination of Sentinel-2 satellite imagery, high-resolution DEM, and IceSat-2 laser altimetry to estimate volume and found their DEM- and laser-altimetry-derived depths to be in close agreement, but that further work to more strongly correlate depth-reflectance and elevation-based (e.g., via the DEM) approaches was needed.

Additionally, methods that model the production of surface meltwater and the routing of it to fill melt-lakes often rely on a surface energy balance model [42,43]. While these models are robust and rely on a physical basis for the estimation of surface melt, they require a large number of input variables, many of which are not available at the necessary spatial or temporal resolution for this study. Here, we apply a conceptually simple Positive Degree Day (PDD) model outlined by Braithwaite (1985) [44] and Hock (2003) [45], which has been improved upon by Tsai and Ruan (2018) [46], to rout surface melt across the watersheds of each of our five lakes for the study period. PDD models rely on fewer input variables—namely near-surface air temperature and, in the case of Tsai and Ruan (2018) [46], percolation depth of meltwater, which serves as a proxy for the depth to which energy is distributed in the upper part of the ice. Deeper percolation will generally produce less surface melting, as incident energy is diluted through a larger ice volume. Utilizing an ablation model that relies on fewer atmospheric variables increases our ability to model the production of surface meltwater across space and time, especially as there are frequent data gaps at in situ weather stations due to the harshness of the environment.

For this study, we compare three methods, two of them new, that quantify melt-lake volume and rely on imagery from the Landsat-8 satellite (30 m spatial resolution), elevation data from the ArcticDEM dataset, and atmospheric data from the Greenland Climate Network (GC-NET) to monitor the evolution of five supraglacial melt-lakes during the 2021 summer melt season (Figure 1). The area is remote and thus has little in situ data on melt-lake depth and volume. Melt-lakes in this region both develop and drain during a short melt season, making it difficult to fully understand the controls on meltwater production and routing, and the formation and filling of melt-lake basins. However, such an understanding is vital for quantifying surface meltwater in the ablation zone and modeling the movement of meltwater through the karst-like environment of the ablation zone [47]. Our work enhances the spatial resolution of the Pope et al. (2016) methodology by using a high-resolution (2 m) DEM to determine lake bathymetry, resolving bathymetric features by an order of magnitude. This work also improves upon the coarse temporal resolution of remote sensing methods by applying a meteorological ablation model to calculate total volume of meltwater within the watershed of each lake. The daily resolution of the ablation model improves upon the coarse resolution of Landsat imagery, which is limited by overpass times and the presence of cloud cover. In addition to increasing spatial and temporal resolution to model melt-lake evolution from our individual methods, this work compares results from these remote-sensing-based and meteorologically forced approaches, to yield a more holistic understanding on the drivers and patterns of seasonal melt-lake evolution. Improvements upon the spatial and temporal resolutions of surface melt-lake volume monitoring, such as we propose in this study, provide a better assessment of the meltwater inputs into the subglacial hydrologic system. This in turn can help provide a better accounting of glacier mass balance across the ablation season, with better information on when, where, and how much melt is leaving the glacier.

2. Study Area and Data

Sermeq Kujalleq is a fast-flowing (~12.5 km yr⁻¹, [48]) glacier in west-central Greenland and has a high number [31] and concentration [49] of supraglacial melt-lakes. At Sermeq Kujalleq’s calving front, this 15 km-wide glacier is responsible a for a significant portion of the ice sheet’s mass loss (~7%) [50,51], making it highly sensitive to its regional surface hydrology (e.g., hydrofracture events). The presence of numerous lakes that meet the hydrofracturing threshold and reappear over multiple years make this region optimal for the monitoring of melt-lake evolution at intra-seasonal and inter-annual time scales.

In the Sermeq Kujalleq Ablation Region (SKAR), melt-lakes form across the region and vary in size, with smaller and deeper lakes at lower elevation due to the relative steepness of the ice sheet surface and larger (area), shallower lakes at higher elevation, due to the GrIS having a gentler slope [49]. Our study area (Figure 1, red box) is a small region (70 km × 20 km) within the SKAR that includes numerous melt-lakes. Within this area, we studied five supraglacial melt-lakes over the 2021 melt season (28 May–23 August) in a narrow elevation band (~150 m) ranging from 1036 to 1185 m a.s.l. (Figure 1), which are located 45–115 km south of the glacier.

Specific lakes were chosen for availability of cloud-free Landsat 8 imagery over the majority of the ablation season, availability of high-resolution Arctic DEM, and similar elevation to that of available AWS data. High-elevation lakes were also intentionally chosen to avoid any hydrologic impact of significant upstream lake formation and drainage.

In this work, we use a combination of Landsat-8 satellite optical imagery, a satellite-derived digital elevation model (DEM), part of the ArcticDEM repository provided by the Polar Geospatial Center, and automated weather station (AWS) data from the Swiss Camp weather station, which is part of the Greenland Climate Network (GC-Net), provided by the Geologic Survey of Denmark and Greenland (GEUS), for the 2021 summer melt season.

2.1. Remote Sensing

2.1.1. Landsat-8 Imagery

Five Landsat-8 Collection 2 Level-2 scenes (

Table 1. Five Landsat-8 images acquired for the 2021 melt season in west-central Greenland, including their path/row identifier.

Image Date	Path/Row
28 May	008/012
22 June	007/012
8 July	007/012
31 July	008/012
23 August	009/012

) spanning the 2021 melt season were acquired from the Google Earth Engine (GEE) platform [52]. These data have a spatial resolution of 30 m and a nominal temporal resolution of 16 days. However, optical imagery of the surface is often obscured by cloud cover, decreasing the effective temporal resolution of the imagery. For this work, we utilized the top-of-atmosphere (TOA) scenes for the red (0.64–0.67 µm) and panchromatic (0.50–0.68 µm) bands to generate the depth-reflectance depth and volume estimates (see Section 3.1). Surface reflectance scenes for the red (0.64–0.67 µm) and blue (0.45–0.51 µm) bands were acquired to compute the Normalized Differenced Water Index for ice (see Section 3.2). Our first image date of the 2021 season—28 May—showed a surface that is snow-covered across our region, with no visible melt-lakes. We therefore determine this date to be the beginning of our melt season such that any lake that formed after 28 May is assumed to have no melt-lake volume. This determination was supported by the AWS data (see Section 2.2), as near-surface temperatures were sufficiently below the melting point during this period.

2.1.2. ArcticDEM

To capture the topographic complexity of the surface lake feature, we use a high-resolution (2 m) digital elevation model (DEM) mosaic provided by the Polar Geospatial Center [53]. The ArcticDEM repository is an open-access dataset where Surface Extraction from TIN-based Searchspace Minimization (SETSM) software (version 4.3.0) is applied to stereoscopic image pairs of high-resolution imagery from the Worldview (-1, -2, and -3 missions) and GeoEye satellites to generate high-resolution DEMs [54].

Because the ArcticDEM dataset is compiled from optical satellite imagery, there are environmental challenges (e.g., cloud cover) in generating DEMs with high spatial and temporal coverage. Additionally, because the melt-lakes form, grow, and either drain or refreeze, careful consideration is made to identify DEMs that cover melt-lakes that are found in several Landsat-8 images during a single melt season, to detect lake growth and evolution. Due to the low density of available ArcticDEM imagery in the region, along with the challenges of finding multiple, clear-sky Landsat-8 images, we selected the ArcticDEM derived from 21 July 2021 Worldview-1 imagery (0.5 m spatial resolution) with an area of ~2100 km².

2.2. Automated Weather Station

Level-1 Automated Weather Station (AWS) data were acquired for the Swiss Camp AWS from 1 May to 31 August 2021 from the Geological Survey of Denmark and Greenland (GEUS). The Swiss Camp station is part of the Greenland Climate Network (GC-Net) and is located at (69°33′53″ N, 49°19′51″ W), with an elevation of 1176 m a.s.l. [55]. We utilize daily mean near-surface (2 m) air temperature to run the ablation model.

3. Methods

We apply three methods, wo remote-sensing-based and one meteorologically based ablation model, to estimate melt-lake volume for our five melt-lakes during the 2021 melt season. The depth-reflectance method, first outlined by Pope et al. (2016), is used as a control comparison to the other two methods, which are new to this study.

3.1. Depth-Reflectance Lake Volume

Supraglacial melt-lake depth estimation using physically based and empirically derived methods is well documented [29,39,40,56]. These approaches account for the attenuation of light as it enters and exits water and carry several implicit assumptions about melt-lake characteristics, including homogenous lake-bottom albedo, flat surface water slope, and minimal dissolved or suspended matter within the water column. Equation (1) uses the formulation originated by Philpot (1989) [56] to solve for our optically derived melt-lake depths, using TOA reflectance data from the Landsat-8 imagery:

$$z={\displaystyle \frac{\left[\ln \left(A_d-R_{\infty }\right)-\ln \left(R_{lake}-R_{\infty }\right)\right]}{g}},$$

(1)

where A_d represents lake bottom albedo, $R_{\infty }$ is the reflectance of optically deep water, R_lake is the reflectance of lake pixels, and g is an attenuation coefficient that is dynamic based on sensor and wavelength, accounting for upward and downward losses of light energy as it passes through the melt-lake water column [29,39,40]. Lake bottom albedo (A_d) is represented by the reflectance values of pixels on the inside edge of the melt-lake where we assume the water column to be at its shallowest. Values for $R_{\infty }$ are determined on each date by acquiring an average reflectance value of optically deep water—typically from the adjacent fjord in our Landsat-8 scenes. The input for Equation (1) is our red and panchromatic band imagery, or R_lake. We borrow the g parameter from Pope et al. (2016), wherein g = 0.7507 and g = 0.3817 for Landsat OLI bands 4 (red) and 8 (panchromatic), respectively [40]. To provide the most accurate melt-lake depth estimates, we follow the Pope et al. (2016) approach by averaging the melt-lake depths calculated, using Equation (1), for the red and panchromatic bands. The geospatial methods outlined above were conducted with the terra package in R and ArcGIS Pro.

3.2. ArcticDEM Lake Basin

Perimeters for our five melt-lakes were delineated using the Normalized Difference Water Index for ice (NDWI_ICE) [57].

$${NDWI}_{ICE}={\displaystyle \frac{blue-red}{blue+red}}$$

(2)

Yang and Smith (2013) [57] sought a range of NDWI_ICE thresholds used to delineate supraglacial streams and lakes in southwest Greenland and found a range of 0.011, 0.008, and 0.021. We applied Equation (2) to our Landsat-8 imagery and identified a range of NDWI_ICE values that closely resemble melt-lake perimeters from Landsat-8 imagery from visual inspection. After applying Equation (2) to determine whether pixels were inundated, that is, part of a melt-lake or not, we needed to further determine whether an individual inundated pixel was part of a melt-lake, that is, not a small pond (i.e., less than 10 Landsat-8 pixels) or a supraglacial stream. To do this, we followed the method of melt-lake sieving outlined by Pope (2016) to remove outlier pixels [58]. Using the patches function of the terra package in R [59], to determine a melt-lake we set a minimum pixel width of 1 and minimum pixel count of 10 and checked images for connected water pixels using 4D connectivity. In other words, sets of connected water pixels < 9000 m² in area are omitted. For consistency across our study period, we used the 0.25 threshold for NDWI_ICE, such that those that do not correspond to a filled melt-lake are removed from the raster.

Using the pixels classified as melt-lake (from the NDWI_ICE method outlined above), the NDWI_ICE output was overlaid onto our DEM, and an extraction of melt-lake elevations was made using the ArcGIS Pro Version 2.8 software package [ESRI], yielding a total of 20 individual melt-lake elevation rasters (one for each lake for each Landsat-8 image date). In ArcGIS Pro, perimeter pixels of our extracted melt-lake elevation rasters were isolated and averaged to determine the melt-lake surface elevation. Using these melt-lake elevations, per-pixel depths were calculated by subtracting raster elevation from the melt-lake surface height and summed to yield melt-lake volume for each melt-lake per image date.

3.3. Ablation Model Lake Basin Filling

To estimate the volume of water produced by atmospheric forcing, we combined remote sensing methods to delineate melt-lake watersheds with automated weather station data run in a meteorologically forced ablation model to generate watershed-wide melt. Here, we assume that all melt generated within the watershed from the model goes into lake basin filling, which, while perhaps an impractical assumption in reality due to evaporation and retention of water in supraglacial streams and interstitial spaces in the firn, nevertheless provides a theoretical maximum for melt-lake volume.

3.3.1. Generating the Melt-Lake Watershed

We delineate individual watersheds for our five melt-lakes by applying a set of hydrology tools (Figure 2) provided by ArcGIS Pro to the ArcticDEM (Figure 2). As there are inherent anomalous elevation values in a stereoscopically derived DEM, to remove individual pixels with anomalous sinks and peaks from the DEM raster, we apply the Fill tool. Identified anomalous pixels are replaced with interpolated elevation values from their neighboring pixels. With our ‘filled’ raster, we apply the Flow Direction tool to determine the direction of water flow for each pixel based on the steepest elevational decline by nearest neighbor. Each pixel is thus assigned a cardinal directional value (e.g., north, southeast) corresponding to the direction in which that pixel would flow based on its steepest neighboring pixel.

We use the Flow Accumulation tool to represent the modeled accumulation of routed surface water from one pixel to its steepest descending neighboring pixel. The accumulation—or cumulative sum of pixels that flow into a particular pixel—increases as water flows downstream. This allows us to visualize the modeled location of surface stream flow, and thus stream order, into and out of melt-lakes. We superimpose our flow accumulation raster onto the Landsat scene of max melt-lake area (8 July) and, using the Landsat-observed melt-lakes, visually identify the location of stream outflow from each melt-lake within the study area (including melt-lakes not studied). Each melt-lake in the area is assigned a pourpoint, which represents the outflow point for each lake within the ArcticDEM area coverage. Providing pourpoints for each melt-lake allows us to use the Watershed function to determine the watershed of each melt-lake identified in the ArcticDEM area. Using this watershed, we assume that during the summer melt season, water produced through surface melt will be routed toward the pourpoint, filling up the melt-lake depression as delineated in the DEM.

3.3.2. Temperature-Driven Ablation Model

The Positive Degree-Day Model (PDD) has been used to estimate surface ablation, using some empirical data [44] and relationships between near-surface temperature and melting of ice and snow across regions [45]. PDD models typically rely on one input parameter (temperature), which is most closely associated with surface melting. While many other parameters are known to affect the occurrence and quantity of surface melt, e.g., insolation [8] and precipitation type and intensity [60], these data are often impossible to find at a spatial and temporal scale appropriate for polar research. The lack of appropriate model forcing renders full energy balance models impractical. Tsai and Ruan (2018) [46] improved upon the simplicity of the PDD model by providing a physical explanation for surface ablation accounting for earlier over-estimates of ice melt. Antecedent cold conditions will delay the onset of melt for several days; therefore, the Tsai and Ruan (2018) model will generally not produce surface melt until several days of above-freezing temperatures have been observed. The ablation rate will vary based on the percolation depth, with faster melting observed coincident with shallower percolation depths. Thus, the Tsai and Ruan (2018) model strikes a balance between a cumbersome full energy balance model and a rudimentary PDD model by accounting for at least one other input parameter aside from temperature. In this work, we utilized the ablation model outlined by Tsai and Ruan (2018) and ran the MATLAB code they provide in their supplemental materials.

Daily mean near-surface temperature (2 m) was collected from the Swiss Camp weather station as input into the ablation model; the station is located 45–120 km north of our five melt-lakes but is comparable in elevation (+/−125 m), and thus should experience similar weather conditions. To account for uncertainty in local air temperatures in each watershed as compared to a more distant AWS, as well as variability in the fracture structure across the ice sheet, several model runs were conducted with temperatures as-is from each weather station as well as temperatures 1 °C, 2 °C, and 3 °C both above and below what was recorded at the AWS. Furthermore, we ran the model with varying percolation depth (2 m, 5 m, 10 m, and 50 m). Varying both of the input parameters is our attempt to account for any possible regional climatic differences between the location of the AWS and the melt-lakes used in our study. The results presented here are derived from our as-is temperature data and a 2 m percolation depth. We report the sum total of the meltwater that is generated on each of the same four dates for which Landsat imagery was available.

4. Results

4.1. Comparison of Remote-Sensing-Based Methods

Volume estimates from our two remote sensing methods (depth-reflectance and DEM) are calculated as the product of our derived depths and the pixel’s surface area. We find general agreement in melt-lake bathymetry between these two methods in that we observe the steepest lake-bed slopes along the perimeter of our melt-lakes and the higher-depth pixels in similar locations within each melt-lake (Figure 3). Figure 3 shows melt-lake depths plotted using the depth-reflectance (left) and DEM (right) methods for Lake C for our peak melt-lake volume date, 31 July. Similar to what we see in Figure 3, a majority of our melt-lake depths over the entire melt period show mean depth (d_µ) derived from the higher-resolution (2 m) DEM method having higher maximum melt-lake depths (d_max) than the Landsat-derived (30 m) depth-reflectance method (Figure 3). For additional comparison, we downscaled the high-resolution ArcticDEM to match the Landsat imagery for Lake C, that is, resampling the DEM from 2 m to 30 m resolution. We found that the downscaled volume was ~3600 m³ less than the volume derived from the original 2 m DEM, or a 0.08% difference in lake volume. Differences in maximum melt-lake depths and overall melt-lake volumes between the DEM and depth-reflectance methods are likely due to the order-of-magnitude difference in spatial resolution between the two datasets (Landsat-8 vs. ArcticDEM) such that the DEM dataset is better able to resolve small-scale undulations at the lake bottom, generating the observed higher mean and maximum melt-lake depths. Additional difference in maximum melt-lake depths is likely attributable to the depth-reflectance method’s dependence on the attenuation of light, thus affecting the A_d, R_∞, and g terms in Equation (1) (see Section 5). Similar findings by Pope et al. (2016) and Melling et al. (2024) noted a plateauing effect in their depth-reflectance approaches for deeper melt-lake depths, which is likely due to the rapid absorption of light in the red band, suggesting that depth-reflectance methods are better suited to shallower melt-lakes.

The smallest difference in melt-lake depth estimates between these methods occurred on 31 July, when peak melt-lake volumes were achieved by four of our melt-lakes (Figure 4). Across all lakes on this date, the depth-reflectance d_µ (DR d_µ) was 0.09 m greater (on average) than the DEM d_µ, whereas on the other image dates, the DEM d_µ across all lakes averaged 0.41 m greater. The date on which there was the least agreement was 8 July, when the DEM d_µ was 0.67 m greater than the DR d_µ. Lakes A and C, which were the deepest and largest in extent, also saw the highest differences in lake depth between the two methods. Comparing lake depths between the two methods is not as straightforward as it might appear. Because the methods offer estimates at different spatial resolutions (30 m vs. 2 m), one would expect substantially higher maximum depths with a finer spatial resolution, and only generally comparable average depths between the two methods.

For each lake, the pattern across the ablation season is the same: low depths in June, maximum depth in July, and drainage of the lake by the end of August. This is also consistent with visible satellite imagery (Figure 5), showing small or no lake extent in June, maximum extent in July, and drainage of the lakes by the end of August. In particular, the DEM method seemed more consistent with lake extent seasonal patterns as seen by visible satellite imagery. For Lakes A, D, and E, the DEM method better captures the change in lake extent from 8 July to 31 July than the depth-reflectance method. For Lakes B and C, the two methods capture the changes in July equally well.

4.2. Performance of the Ablation Model Method

For each of the five melt-lakes, the volumes calculated through the ablation model correlated poorly to those of the other methods. Furthermore, it is clear from Landsat-8 imagery that lakes are present and of significant size despite the ablation model showing little meltwater having been generated (Figure 6). Because we are reporting the sum total of all the meltwater generated from 1 May through the time of each image date, the ablation model should serve as a theoretical maximum for the volume of water in each of our melt-lakes. While we know that there is evaporation and percolation of meltwater into the ice, this is not represented by the Tsai and Ruan (2018) model [46]. However, for four of the five lakes in this study, the ablation model shows nearly an order of magnitude less melt than either of the remote-sensing methods for 8 July. The ablation model generally shows larger lake volumes for the late-August time period, but again, because we are not accounting for sudden lake drainage via hydrofracturing, we are only able to obtain the theoretical maximum lake volume by that date. Still, the maximum lake volumes generated by late August fall far below the volumes that we calculate at maximum lake size via these other methods.

The one lake for which the ablation method was somewhat consistent with the remotely sensed methods was Lake D. By 31 July, the ablation model had generated enough melt to fill the lake at a volume nearly equal to the depth-reflectance method. However, there was still not enough melt generated by 8 July to fill the lake to the spatial extent observed in the Landsat imagery.

We report here model runs that assume a minimal percolation depth of 2 m and temperatures that are equal to what was observed at the Swiss Camp AWS. Because the melt-lake watersheds are roughly 100 km south of Swiss Camp, and some of the watersheds are somewhat lower in elevation, we repeated the model runs with temperatures +1 °C, +2 °C, and +3 °C from the AWS data. Each successive run, while producing somewhat more melt, still did not produce a pattern that was consistent with a large lake forming by early July. Additionally, we forced the ablation models with temperatures colder than the AWS station, and with percolations depths up to 50 m, but these model runs were even less likely to generate enough melt to form the lakes observed in the satellite imagery.

During the study period, we found that four of our five melt-lakes reached the Krawczynski et al. (2009) [13] threshold for hydrofracturing, suggesting that these melt-lakes had the potential for hydrofracture during the 2021 melt season.

5. Discussion

We present two novel methodologies and compare them to a depth-reflectance method to calculate changes in melt-lake volume for five study lakes in the SKAR through the 2021 melt season. Our results show a range of melt-lake volumes between these three methods. The evolutional pattern of melt-lake volumes was detected by both the depth-reflectance and DEM–lake-filling methods, in that melt-lakes started off smaller (22 June) and grew and peaked in our July image dates, and then, by our late-August date, melt-lakes either slowly or catastrophically drained. We found differences in melt-lake maximum depths between our two remote sensing methods, such that the DEM-derived depths produced higher maximums depths for most of our lakes in all image dates. This is likely due to the depth-reflectance method’s reliance on upwelling light in the red and panchromatic bands. Others [40,41] have found similar results when using the depth-reflectance method derived from the red band. As there is rapid attenuation of red light within a water column, using red band reflectance can on one hand resolve the shallow melt-lake depths well, but after reaching a limit, becoming a hinderance in accurately capturing deeper melt-lake depths. Another source for the higher maximum depths by the ArcticDEM dataset can arise from it being at a higher spatial resolution—relative to the Landsat imagery—which would better resolve lake bottom undulations that are not captured in the coarser Landsat dataset. For the ablation model method, we found that it did not realistically or adequately show the same evolution of melt-lake extent and volume.

The formation of surface melt-lakes has important implications for when and where pulses of meltwater leave the Greenland Ice Sheet, and a better understanding of surface ablation is important for understanding the regional and global-scale implications of global climate change, specifically with respect to sea level rise and the strength of Atlantic Meridional Overturning Circulation. While further work will be necessary to upscale the contribution of melt from individual lakes to the scale of the entire Greenland Ice Sheet, this work represents an important step forward in constraining the spatial and temporal contribution of surface melt.

Considerations for the differences in our volume estimates between these methods, as well as known limitations of each method, are discussed below.

5.1. Considerations in Our Remote Sensing Approaches

5.1.1. Cloud Cover

Often-relied-upon satellites, such as MODIS, Sentinel-2, ASTER, and Landsat, orbit high above the earth’s atmosphere and at fixed temporal resolution. Capturing imagery above the altitude of clouds can obscure the ablation zone surface, effectively reducing the temporal resolution of the remote-sensing-based methods presented in this study. Additionally, partial cloud coverage within an image may be radiometrically affected by thin cirrus clouds that are capable of escaping cloud masking techniques. The presence of such clouds affects the detected at-sensor energy of the melt-lake. Cloud presence thus modifies the top of atmosphere pixel reflectance by increasing the reflectance value and affecting the estimated inundation depth. This leads to a bias in inundation depth by generating shallower melt-lakes and decreasing volume totals.

5.1.2. Floating Ice

A common phenomenon that occurs in the ablation zone is the presence of floating ice inside of a melt-lake or marginal freezing, that is, refreezing of surface water along the margins of the melt-lake. The spectral reflectance of these objects mirrors that of nearby surface ice outside of the melt-lake, so that optically derived methods are unable to estimate melt-lake depth below their corresponding pixels. As the method for ArcticDEM generation is reliant on optical imagery, this too affects the derived surface elevation corresponding to pixels associated with the floating ice or marginal freezing. The inability of these methods to account for inundation depth below these features will yield lower melt-lake volume.

5.1.3. NDWI Thresholding

Lake perimeter delineation is defined using NDWI-defined thresholds from the Landsat-8 imagery. We found little area change, that is, the number of pixels classified as melt-lake, using several NDWI_ICE thresholds (0.18 ≤ NDWI_ICE ≤ 0.25) that fall within the parameters described by Pope (2016) [57]. However, due to the order-of-magnitude difference in spatial resolution between Landsat-8 imagery (30 m) and the ArcticDEM (2 m), small changes in the Landsat-derived lake boundaries had a substantial effect on melt-lake surface height as extracted by the ArcticDEM, strongly modulating overall lake volume.

5.2. Considerations in Our Ablation Model Approach

5.2.1. Temperature Anomalies from AWS to Melt-Lake Watersheds

Perhaps the most obvious reason that a meteorologically forced ablation model would poorly correspond to observed melt is that the temperatures used to force the model could be significantly different from those experienced in the field. With the Swiss Camp AWS being roughly 45–120 km from each of the watersheds, this was an initial concern. However, temperature anomalies alone cannot account for the significant discrepancy between modeled and observed ablation rates, particularly early in the season. We attempted to account for this variability by increasing the temperatures in the model inputs to levels beyond what one might expect for temperature variability (+3 °C) and were still unable to model melt rates consistent with what was observed in satellite imagery. We also used AWS data from the KAN-M station to force the ablation model. The amount of melt generated was similar to that from the forcing from the Swiss Camp AWS. We then looked at other physically based constraints that might be more plausible reasons for the lack of modeled melt. Furthermore, when Tsai and Ruan (2018) calibrated their model, they also contended with a limited surface temperature record—using NCEP/NCAR Reanalysis [61], which is available only at a 2.5° × 2.5° spatial resolution.

5.2.2. Limitations of PDD Ablation Model over Full Energy Balance Ablation Models

Insolation (incoming solar radiation) and temperature have been shown to be the primary drivers of melt-lake formation in Greenland [8]. While the PDD model does account for temperature, insolation is more difficult to model. The microscale extent of cloud cover can have a significant effect on insolation, and because AWS data, let alone insolation, were unavailable at the scale of individual lake watersheds, this constrains observed surface melting. While insolation can be obtained from reanalysis output such as the NCEP/NCAR Reanalysis [61] or North American Regional Reanalysis [62], this information is still only available at a fairly coarse spatial resolution as compared to watershed and lake size.

Furthermore, the modified PDD model does not account for the significant influence of rainfall on surface melting events. Rain-on-snow events have been observed to increase the temperature of the firn by as much as 23 °C in one day [60], which would rapidly overwhelm the antecedent cold ice conditions and generate melt far more quickly than would be generated by the modified PDD model. While Tsai and Ruan (2018) did apply their model to a glacier in southwestern Greenland (Qamanarssup Sermia), their calibration year was chosen specifically for a more complete early-season ablation record. The chosen year was 1986, which was significantly cooler than the year in our study, 2021.

5.2.3. Non-Meteorological Constraints on Lake Volume

It is also certainly possible that non-meteorological processes are accelerating the filling of the melt-lakes beyond what is generated in the modified PDD ablation model. While we can constrain the watersheds at the surface via high-resolution DEM, we do not have an accurate picture of whether these watersheds are connected under the ice surface. Connectivity of watersheds could conceivably conduct excess meltwater into the lakes; however, it could just as likely reduce the amount of available water through transport out of the watershed. We also considered that the modified PDD ablation model did not account for the time that it takes for water to travel across the watershed and into the melt-lakes. But again, this would likely lead to a decrease in meltwater observed in each lake, as water would be lost to evaporation and downward percolation over the course of travel. We also have considered that perhaps the incoming solar radiation is concentrated in a much shallower layer of ice than was estimated by our modified PDD model. This would be possible in cases where there is little accumulated snow, perhaps in areas with less snowfall during the cold season or ice exposed by high winds and snow blown out of the watershed.

6. Conclusions

We applied several methods to estimate surface meltwater storage in five melt-lakes in the ablation zone south of Sermeq Kujalleq over the 2021 melt season. These methods utilized a combination of satellite-derived and meteorological data to represent the rapid growth of melt-lake volume across our four image dates. While the DEM infilling method is in general agreement with the depth-reflectance method, our ablation model method significantly underestimates lake volume. The DEM infilling method provides information about lake volume at a higher spatial resolution (one order of magnitude higher) than is currently possible using the depth-reflectance methodology. We found that during this period, four of our five melt-lakes achieved sufficient melt-lake volume to exceed the Krawczynski et al. (2009) threshold for hydrofracture. This suggests that these melt-lakes had the potential to hydrofracture, that is, catastrophically drain, which can affect nearby ice glacial velocity and regional ice flow dynamics. Further work will be necessary to rectify the spatial scale mismatch between higher-resolution remote sensing methodology and coarser-resolution regional and global climate models.

Challenges in the selection of melt-lakes for study included the lack of in situ data on melt-lake depth to more precisely correlate our remote-sensing-based estimates, along with the need for overlapping data from Landsat-8 overpasses and the ArcticDEM dataset. Despite these challenges, we observed similar trends in melt-lake volume from our two remote sensing methods (depth-reflectance and DEM) such that lakes began small at the melt season onset, grew during July, and then decreased in size by late August. More work is necessary to determine the validity of the DEM infilling method with respect to deeper, smaller lakes that are closer to the ice sheet margin. Additionally, these lakes have a more complicated watershed which may fill and/or drain at a temporal scale shorter than the effective temporal resolution of either the depth-reflectance or DEM infilling methods. As such, it is possible that modeling of lake volume may provide the appropriate temporal resolution to study these features.

However, when comparing the volumes of each of the remote-sensing-based methods to those estimated by the modified PDD ablation model method, we found that there was a significant delay in the production of meltwater. The melt that was produced was not sufficient to represent melt-lake areas as observed in the Landsat-8 imagery until after the lakes were observed to have drained at the end of the ablation season. This suggests that a PDD-based ablation model, without benefit of additional energy balance parameters such as insolation or rainfall, is likely insufficient on its own to model supraglacial lake volume in this region. It is possible that the influence of rainfall was greater in 2021 than it might have been in other years, owing to the observed rainfall event at the summit station in mid-August. This excess rainfall may have contributed to an unusual amount of water in each of the lakes. However, the disparity was so great between the large size of observed melt-lakes and the paucity of modelled meltwater that a more thorough energy balance model would seem to be necessary to adequately capture melt dynamics at this scale.

Here, we described sources of agreement and disagreement between our three approaches in the hope of facilitating better estimation of supraglacial lake volume, as the dynamic production, routing, and storage of surface meltwater on the ice sheet have significant impact on ice sheet mass balance. Each of our three methods has an advantage for melt-lake volume estimation. Depth-reflectance approaches rely on satellite imagery with repeat coverage and can provide longitudinal studies on melt-lake evolution at decadal time scales. Our DEM method is not sensitive to attenuation of light, particularly at red wavelengths, and can more accurately model deeper melt-lakes, which contain larger volumes. The ablation model is not sensitive to obstructive cloud cover that would hinder the acquisition of satellite imagery and can provide melt-lake development at much higher temporal resolution. However, this particular approach will require a more intensive energy balance model to account for rapid increases in lake volume during rain-on-ice events. Despite the disagreements in melt-lake volume estimates between our three methods, this work provides a framework for integrating remote sensing and meteorological approaches to melt-lake evolution to yield a more comprehensive understanding—at higher spatial and temporal resolutions—of the drivers of meltwater production and its influence on the spatial distribution and extent of supraglacial melt-lakes in the ablation region.