Seasonal evolution of supraglacial lakes and rivers on the southwest Greenland Ice Sheet

Abstract

Supraglacial lakes and rivers dominate the storage and transport of meltwater on the southwest Greenland Ice Sheet (GrIS) surface. Despite functioning as interconnected hydrologic networks, supraglacial lakes and rivers are commonly studied as independent features, resulting in an incomplete understanding of their collective impact on meltwater storage and routing. We use Landsat 8 satellite imagery to assess the seasonal evolution of supraglacial lakes and rivers on the southwest GrIS during the 2015 melt season. Remotely sensed meltwater areas and volumes are compared with surface runoff simulations from three climate models (MERRA-2, MAR 3.6 and RACMO 2.3), and with in situ observations of proglacial discharge in the Watson River. We find: (1) at elevations >1600 m, 21% of supraglacial lakes and 28% of supraglacial rivers drain into moulins, signifying the presence of high-elevation surface-to-bed meltwater connections even during a colder-than-average melt season; (2) while supraglacial lakes dominate instantaneous surface meltwater storage, supraglacial rivers dominate total surface meltwater area and discharge; (3) the combined surface area of supraglacial lakes and rivers is strongly correlated with modeled surface runoff; and (4) of the three models examined here, MERRA-2 runoff yields the highest overall correlation with observed proglacial discharge in the Watson River.

1. Introduction

Extensive fluvial supraglacial hydrologic networks are observed across the southwest Greenland Ice Sheet (GrIS) ablation zone (Smith and others, 2015, 2017; Yang and Smith, 2016; Pitcher and Smith, 2019). These networks, consisting of interlinked supraglacial streams, rivers and lakes, drain large volumes of surface meltwater into the ice sheet via moulins and crevasses (Chu, 2014; Flowers, 2018; Pitcher and Smith, 2019), with corresponding impacts on GrIS subglacial water pressures, and ice flow dynamics (Zwally and others, 2002; Bartholomew and others, 2011; Hoffman and others, 2011; Andrews and others, 2014).

Interconnected supraglacial lakes and rivers are particularly extensive on the southwest GrIS (Fitzpatrick and others, 2014; Smith and others, 2015; Flowers, 2018). In concert with a non-channelized, often porous, near surface weathering crust atop bare ice (Cooper and others, 2018), they control surface meltwater storage and routing on this area of the ice sheet (Smith and others, 2017). Supraglacial lakes are located in surface topographic depressions (Liang and others, 2012) and are estimated to impound 5–13% of instantaneous surface meltwater during the melt season (Leeson and others, 2012; Arnold and others, 2014; Fitzpatrick and others, 2014). While most of these lakes drain laterally into supraglacial rivers over timescales of weeks to months (Tedesco and others, 2013; Smith and others, 2015), some drain vertically via hydrofracture and moulin formation over hours to days (Das and others, 2008), and others retain meltwater overwinter (Koenig and others, 2015; Lampkin and others, 2020). Both rivers and lakes inject meltwater into the ice sheet via moulins and crevasses, thus modulating basal water pressure and evolution of the subglacial hydrologic system (Hoffman and others, 2011; Andrews and others, 2014).

Supraglacial stream/river networks exhibit complex drainage patterns and high drainage densities (Rippin and others, 2015; Smith and others, 2015; Yang and others, 2016). These well-developed fluvial networks route surface meltwater into lakes and moulins (Banwell and others, 2013), and the shape/area of their supraglacial catchments helps control the amplitude and timing of meltwater penetration into the GrIS (Smith and others, 2017).

Despite functioning as interconnected hydrologic networks, supraglacial lakes and streams/rivers are often investigated as independent features. Consequently, the relationship between surface meltwater storage and routing remains poorly understood (Yang and others, 2018). For example, it remains unknown if the areas and volumes of supraglacial lakes and rivers correlate differently with surface melting, and what the relative storage capacities of lakes and river are. Such measures could reduce uncertainties in both surface runoff simulations (Smith and others, 2017; van As and others, 2017; Yang and others, 2018) and ice flow models that incorporate surface meltwater forcing (Andrews and others, 2014; Flowers, 2018).

Current understanding of the coupling between supraglacial lakes and rivers is largely derived from topographic and conceptual models. Banwell and others (2012) developed the Surface Routing and Lake Filling (SRLF) model to simulate surface meltwater routing through the base of a saturated snowpack or over bare ice to fill supraglacial lakes. For a saturated snowpack, supraglacial streams/rivers evolve as the snowline retreats; for a bare-ice surface, SRLF routes supraglacial streams/rivers over a theoretically impervious surface assuming a constant routing network and efficiency (Arnold and others, 1998). However, supraglacial streams/rivers are found to evolve seasonally (Lampkin and VanderBerg, 2014; Yang and others, 2017). This leads to varying spatial partitioning between slow-flowing interfluve (i.e. hillslope) and fast-flowing open-channel zones, and thus changing routing efficiency over the course of the melt season. To investigate this, Yang and others (2018) arbitrarily varied contributing area thresholds (100–5000 m²) in a digital elevation model, but this approach fails to capture the true seasonal evolution of the supraglacial hydrologic networks or the seasonal infilling of lakes.

In this study, we use seven medium-resolution (30 m) Landsat 8 Operational Land Imager (OLI) satellite images to map the seasonal evolution of surface meltwater area and volume throughout the 2015 melt season, a colder than average year. Rather than separating lakes from rivers, we quantify their integrated seasonal evolution. Next, we compare remotely sensed meltwater area and volume with surface runoff simulations from three climate models (MERRA-2, MAR 3.6 and RACMO 2.3). Special attention is paid to high-elevation (>1600 m) areas where well-developed supraglacial hydrologic networks appear to develop despite their location above the equilibrium line altitude (ELA) (~1550 m in this area, van de Wal and others, 2012; Poinar and others, 2015).

2. Study area and data

Two major drainage basins, sourcing the Isortoq and Watson (Qinnguata Kuussua) proglacial rivers, were selected as the study area (×Fig. 1). Each summer, supraglacial lakes and streams/rivers from across the icescapes of these two basins (Fitzpatrick and others, 2014; Smith and others, 2015; Yang and Smith, 2016). Their boundaries were derived using CryoSheds (Pitcher and others, 2016) and their combined basin ice-surface area below 2000 m (a conservative upper limit of melt extent) was calculated as ~11 000 km² (~6000 km² for the Isortoq river basin and ~5000 km² for the Watson river basin, respectively).

Fig. 1. Landsat 8 OLI image of the Isortoq and Watson (Qinnguata Kuussua) river basins, southwest GrIS (acquired 8 July 2015, R5G4B3). Black crosses show the locations of three PROMICE automatic weather stations (AWS) used in this study.

Seven multi-spectral Landsat 8 OLI satellite images covering the study area from 13 June to 25 August 2015 were acquired from the United States Geological Survey (USGS, https://www.usgs.gov/) (×Table 1). While large supraglacial lakes are often mapped using 250 m spatial resolution MODIS images (Liang and others, 2012), Landsat 8 offers superior lake mapping accuracy due to its 30 m spatial resolution and available spectral bands (Yang and others, 2015b; Pope and others, 2016; Williamson and others, 2018a).

Table 1. Seven multispectral Landsat 8 images acquired over the study area (Isortoq and Watson river basins, with a combined ice-surface area of 10 946 km² below 2000 m) during the 2015 melt season

Multi-temporal supraglacial lakes and rivers were mapped and surface meltwater areal extent and storage were calculated

^a Although cloud cover on the 24 July 2015 Landsat 8 image is high, most clouds are located to the north of our study area and have slight impacts on supraglacial lake and river mapping in our study area.

We specifically focus on the 2015 year rather than inter-annual variations between years. Note that surface runoff during 2015 was generally lower than the 2000–2018 average (×Fig. 2). As such, surface melting in 2015 began one month later than is typical for southwest GrIS (early June rather than early May) and runoff during August was similarly below average. Overall, compared to the last two decades, the 2015 melt season can be considered a lower limit of GrIS surface runoff and meltwater extent. Outside of a modest runoff pulse in September, the seven Landsat 8 images analyzed here are broadly representative of the 2015 melt season.

Fig. 2. Modeled (a) daily and (b) cumulative surface runoff of the Isortoq and Watson river basins during the colder-than-average 2015 melt season as simulated by the Modèle Atmosphérique Régionale (MAR) 3.6 regional climate model. For context, modeled values for the extremely warm melt year 2019 and the long-term 2000–2019 average are also presented. Gray uncertainty envelope presents the one std dev. range of surface runoff. Red dots mark the 2015 image acquisition dates for seven Landsat 8 images, while the two blue dots represent the 2019 acquisition dates two higher-resolution images (i.e. a 2 m WorldView-2 image acquired on 25 July 2019 and a 10 m Sentinel-2 image acquired on 4 August 2019, see Figs S2, S3).

Daily ice-sheet surface runoff values for the Isortoq and Watson river basins were extracted from three widely-used climate models (MAR 3.6, RACMO 2.3 and MERRA-2). Detailed descriptions of these climate models are provided by Smith and others (2017). Daily estimates of proglacial discharge in the Watson River were obtained from van As and others (2017), and were compared with daily surface runoff simulations from the climate models. In situ measurements of ice motion were obtained from three Programme for Monitoring of the GrIS (PROMICE) automatic weather stations (AWS) located in the study area. From low to high elevations, they are KAN_L, KAN_M and KAN_U, located at 670, 1270 and 1840 m elevation, respectively (Fig. 1). In addition to their meteorological sensor suites, these stations contain GPS units used to calculate ice surface motion. We used a 5 d smoothing window to eliminate outlier velocity measurements. Additionally, the AWS are equipped with sensors to measure surface accumulation and ablation.

3. Methods

3.1 Surface meltwater area delineation

Supraglacial lakes and rivers were automatically delineated from Landsat 8 images using a Normalized Difference Water Index (NDWI) adapted for ice (NDWI_ice = (BLUE–RED)/(BLUE + RED)) as per Yang and Smith (2013). A high-value, global NDWI_ice threshold (t _ndwi = 0.25) was applied to delineate supraglacial lakes and wide river segments, and narrower river segments were then enhanced from the NDWI_ice images using Gabor filtering and path opening (Yang and others, 2015a). Finally, a second global threshold (t _watermask = 25 out of 255) was applied to classify the enhanced images and to create binary surface meltwater masks.

The two thresholds, t _ndwi and t _watermask were set relatively high, rendering conservative delineation of surface meltwater. We selected t _ndwi = 0.25 because it captures supraglacial lake boundaries reasonably well from all seven Landsat 8 images. This global t _ndwi threshold value also performs well on the Antarctic Ice Sheet (Bell and others, 2017), and is similar to a value of 0.23 previously used to detect southwest GrIS supraglacial lakes (Moussavi and others, 2016). In addition, Moussavi and others (2020) used t _ndwi = 0.19 to identify Antarctic supraglacial lake boundaries and then used additional band calculations to discern lakes from cloud-shadowed snow surfaces. While dynamic thresholding techniques have been used to map supraglacial lakes (Selmes and others, 2011; Williamson and others, 2017), we find that a global threshold t _ndwi = 0.25 is satisfactory given the limited spatial domain and number of images used here.

3.2 Identification of supraglacial lakes, rivers and moulins

Supraglacial lakes and rivers are inseparable in our binary meltwater masks. To broadly partition them, lake polygons and river polylines were created using the ArcScan vectorization tool (ESRI, 2010). The maximum polyline width was set to 5 pixels (150 m) based on visual inspection of small supraglacial lakes. As such, surface meltwater bodies wider than 5 pixels were summarily categorized as supraglacial lakes, and narrower bodies summarily categorized as supraglacial rivers.

This threshold likely miscategorizes some small supraglacial lakes (Williamson and others, 2018a). However, the impact of this on overall surface meltwater area and volume calculations is small. For example, a miscategorized rectangular lake having width × length of 5 × 50 pixels (i.e. 150 m × 1500 m) has a relatively small surface area of 0.225 km². Yang and others (2015b) showed that the summed area of small (<0.25 km²) lakes is <15% of the total lake area on the southwest GrIS. Similarly, Williamson and others (2018a) found that meltwater drained by small (<0.125 km²) lakes represents only 5% of the total meltwater volume drained by lakes on the western GrIS (north of our study area). Finally, numerous small lakes are correctly categorized using a width threshold of 5 pixels (see Fig. S1). We conclude that the impact of small, miscategorized lakes on computed surface meltwater areas and volumes is small.

Small supraglacial streams are ubiquitous on the southwest GrIS yet are largely undetected in 30 m resolution Landsat 8 images (Yang and Smith, 2016; Yang and others, 2019b). Although small streams enhance surface meltwater routing (Decaux and others, 2018; Yang and others, 2018), coarser images of large supraglacial rivers capture the overall drainage pattern of supraglacial hydrologic networks effectively (see Table S1 and Fig. S2). These images are also a primary requirement for identifying terminal moulins (Smith and others, 2015; Yang and Smith, 2016). Therefore, 30 m Landsat 8 imagery is sufficient for capturing the broad-scale seasonal evolution of supraglacial hydrologic networks on the southwest GrIS.

To investigate high-elevation surface-to-bed meltwater connections, we mapped terminal supraglacial river moulins at elevations >1600 m in a 2 August 2015 Landsat-8 image. We selected early August for this analysis as moulins are best developed near the end of melt season (Yang and Smith, 2016), yet surface runoff remains sufficiently high to enable their detection in visible/NIR satellite imagery.

3.3 Surface meltwater volume estimation

Surface meltwater volumes of supraglacial lakes and rivers were estimated from their spectral characteristics. To retrieve water depths, we used the band ratio between the coastal and green bands of Landsat 8 following Pope and others (2016). This approach uses empirical parameters calibrated from in situ spectral and bathymetric data (Pope and others, 2016) as follows: z = a + bX + cX ², X = ln(B1/B3), where a = 0.1488, b = 5.0370 and c = 5.0473, B1/B3 is the ratio of the reflectance of Landsat 8 OLI band 1 (coastal) and 3 (green), and z is the derived lake depth. Meltwater volumes were calculated by multiplying Landsat 8 pixel area (900 m²) by the sum of water depth derived for water pixels (Yang and others, 2019a). Meltwater volume uncertainties were calculated by multiplying RMSE = 0.38 m (Pope and others, 2016) by pixel area and the number of water pixels, as per Williamson and others (2018a).

3.4 Quantitative metrics of supraglacial hydrology

The resultant map products were surface meltwater mask and depth rasters, and lake and river vectors. The areas and depths of supraglacial lakes and rivers exhibit large spatial and temporal variations over the melt season (Chu, 2014; Fitzpatrick and others, 2014; Yang and Smith, 2016; Williamson and others, 2018a). To characterize this variability, we analyzed supraglacial lakes and rivers at 200 m elevation intervals (<1000, 1000–1200, 1200–1400, 1400–1600, 1600–1800 and >1800 m). These six elevation intervals span differing melt intensities and surface conditions (Sundal and others, 2009; Fitzpatrick and others, 2014). For example, the ice surface below 1200 m is heavily crevassed with relatively few lakes and rivers (Colgan and others, 2011), whereas the ice surface between 1200 and 1600 m is smoother, with numerous supraglacial lakes and rivers (Fitzpatrick and others, 2014; Yang and Smith, 2016). Above 1600 m, the ice surface is typically covered by firn/snow, with an end-of-season snowline elevation of 1520 ± 113 m in 2015 (Ryan and others, 2019).

Two quantitative metrics, surface meltwater area proportion and surface meltwater storage, were computed over the 2015 melt season. Surface meltwater area proportion (%) was calculated by dividing the total surface meltwater area by the study area, and surface meltwater storage (mm) was calculated by dividing the total meltwater volume by the study area. For each satellite image, these two metrics were calculated for each elevation interval. Supraglacial river width was not included because, due to mixed pixel effects, 30 m Landsat 8 images overestimate the widths of supraglacial rivers (Fig. S3).

For each elevation interval, a mean supraglacial lake (river) depth was calculated from each Landsat 8 image, yielding times series of instantaneous surface meltwater volumes binned by elevation. Mean depths and std dev. over time were also calculated for each elevation interval to compute seasonally-averaged storage volume and variability, respectively. The coefficient of variation (c _v = std dev./mean) was used to investigate the dispersion of lake and river depth.

We defined a ‘lake disappearance ratio’ as:

(1)$$\left({\displaystyle{{{\rm maximum\ lake\ count}\ndash {\rm lake\ count\ on\ 2\ August}} \over {{\rm maximum\ lake\ count}}}} \right)$$

This ratio quantifies the relative ratio of lakes that drain via moulins or lateral outlet channels by the end of the melt season. We used the supraglacial lake count during late melt season as the denominator because most supraglacial lakes have drained by this time (Selmes and others, 2011). We repeated the calculations for large (0.25 km², the minimum lake area identifiable in MODIS satellite imagery) lake drainage events as they have a large impact on subglacial hydrologic system evolution (Selmes and others, 2011; Andrews and others, 2018).

3.5 Comparison of satellite-observed supraglacial hydrology with modeled surface runoff

We hypothesize that the seasonal evolution of supraglacial lakes and rivers is correlated with the seasonal surface meltwater runoff. To test this hypothesis, we calculated the correlation between remotely sensed quantitative metrics of surface meltwater and modeled surface runoff. The modeled surface runoff was binned into six elevation intervals (see section 3.4) and used to create linear regression models with remotely sensed surface meltwater area, similar to previous studies (Sundal and others, 2009; Fitzpatrick and others, 2014; Yang and Smith, 2016).

Large uncertainties in modeled surface runoff have been reported (Smith and others, 2017; Yang and others, 2019a). Although meltwater retention in snow/firn/weathering crust (Rennermalm and others, 2013; Machguth and others, 2016; Steger and others, 2017; Cooper and others, 2018; Vandecrux and others, 2019) is not considered here, we suggest that remotely sensed surface meltwater reflects ice surface melting and runoff conditions. Therefore, we explored whether satellite-derived supraglacial lake and river maps correlate with modeled surface runoff, and also compared them with in situ measurement of proglacial discharge.

4. Results

4.1 Seasonal evolution of supraglacial lakes and rivers

Supraglacial lakes and rivers progressively migrate from low (<1000 m) to high (>1600 m) elevations in our study area (×Fig. 3), consistent with previous studies (Sundal and others, 2009; Lampkin and VanderBerg, 2014). At low elevations (<1200 m), both supraglacial lake and river area proportions peak on 1 July. Above 1200 m, supraglacial river area proportions peak 7–10 d earlier than supraglacial lakes (Figs 4a, b). This time-lag has not previously been reported, but is not surprising given that supraglacial topographic depressions are small and surface runoff intense below 1200 m (Yang and others, 2015b), resulting in rapid small (mean area 0.10 km²) lake infilling. Above 1200 m, topographic depressions are larger (mean area 0.36 km²) so lake infilling takes longer. In contrast, supraglacial rivers fill quickly because they are highly responsive to surface runoff generation, evacuating runoff from the ice-sheet surface within hours (Smith and others, 2017; Yang and others, 2020).

Fig. 3. Seasonal evolution of supraglacial lakes and rivers on the southwest Greenland Ice Sheet as mapped from seven Landsat 8 OLI satellite images acquired throughout the 2015 melt season. Panels (a)–(g) show blue supraglacial lakes and rivers with the boundaries of the Isortoq and Water river basins overlaid. Panel (h) presents a temporal composite of all seven maps, revealing the temporal evolution of supraglacial lakes and rivers from low to high elevations.

Fig. 4. Seasonal evolution of (a) river area proportion, (b) lake area proportion, (c) surface meltwater area proportion, (d) surface meltwater storage, (e) daily MERRA-2 runoff and (f) cumulative MERRA-2 runoff, (g) daily MAR runoff and (h) cumulative MAR runoff, (i) daily RACMO runoff and (j) cumulative RACMO runoff, (k) ice flow velocity, and (l) surface ablation measured by three PROMICE stations (KAN_L, KAN_M and KAN_U in Fig. 1).

4.2 High-elevation supraglacial lake and river evolution

High-elevation (>1600 m) supraglacial hydrologic networks become well-developed by 17 July (Figs 3, 4a–d). Peak surface meltwater area proportions (lakes and rivers combined) between 1000–1800 m elevation are all ~7% (Fig. 4c). Peak river area proportion is 4.1–5.0% above 1400 m and 4.6–5.9% below 1400 m. Peak lake area proportion at 1400–1600 m (4.5%) is higher than below 1400 m (2.3–3.0%). Peak lake area proportion at 1600–1800 m (2.6%) is ~87% of its low-elevation counterpart at 1200–1400 m (3.0%), more than quadruple than occurred during the 2003 melt season (~20%, Sundal and others, 2009). Moreover, peak surface meltwater storage at 1600–1800 m is ~104% that of at 1200–1400 m (Fig. 4d), more than double that of the Watson river basin for the 2010 melt season (~40%, Fitzpatrick and others, 2014). This confirms that supraglacial lakes and rivers advanced further inland and up-elevation in summer 2015, consistent with model predictions (Leeson and others, 2015) and satellite observed trends from 1972 to 2012 (Howat and others, 2013).

Above 1600 m elevation, six moulins are observed on 2 August 2015 (Fig. 3f). Their combined upstream river length is 302 km, accounting for 28% of the total mapped high-elevation river length (1078 km). This signifies the presence of surface-to-bed meltwater connections at high elevations, consistent with previous reports (Yang and Smith, 2016; Gagliardini and Werder, 2018). A total of 198 km (18%) of high-elevation river length is observed draining into terminal supraglacial lakes, with 579 km (54%) flowing to lower elevations where most rivers terminate in moulins between 1400 and 1600 m elevation (Fig. 3f).

4.3 Surface meltwater area, depth and volume

Total surface meltwater area proportion and storage increase after 13 June, peak on 17 July, then decrease through 25 August. The peak meltwater area proportion and storage are 3.4% and 66.5 mm, respectively (Table 1). Averaged across all elevations and dates, the mean supraglacial lake depth is 1.7 ± 0.7 m and mean supraglacial river depth is 0.9 ± 0.2 m. The coefficient of variation of meltwater depth is 0.41 for supraglacial lakes and 0.22 for supraglacial rivers, suggesting that lake depth is generally more variable than river depth.

Peak meltwater area proportion is primarily controlled by supraglacial rivers, while peak meltwater storage volume is primarily controlled by supraglacial lakes. Both display large variations at different elevation intervals throughout the 2015 melt season (×Figs 4a–d, 5). Other than 1400–1600 m elevation on July 17, the area proportions of rivers are consistently larger than those of lakes (Figs 4a, b). The mean contribution of rivers to peak meltwater proportion is 62%, ranging from 41% at 1400–1600 m to 81% at 1200–1400 m.

Fig. 5. Remotely sensed instantaneous surface meltwater storage in (a) lakes and rivers, (b) lakes only and (c) rivers only, during the 2015 melt season. Error bars are obtained by multiplying RMSE = 0.38 m of water depth estimation using Landsat 8 coastal and green bands (Pope and others, 2016) by the Landsat 8 pixel size (900 m²) and the number of water pixels.

Peak meltwater storage is primarily attributed to supraglacial lakes. Lake storage is larger than river storage in all satellite images except for 1600–1800 m elevation on 17 July and 1200–1400 m elevation on 1 July (Fig. 5). The mean contribution of supraglacial lakes to peak meltwater storage is 68%, ranging from 46% at 1600–1800 m to 81% at 1400–1600 m (Fig. 5b). This is because topographic depressions can store large volumes of surface meltwater (Arnold and others, 2014; Fitzpatrick and others, 2014), whereas supraglacial rivers are shallow and their instantaneous channel storage is small (Legleiter and others, 2014; Gleason and others, 2016). However, for all elevation intervals, supraglacial rivers contain small but significant instantaneous meltwater storage in the early melt season (Figs 3, 5c). We also note that meltwater flux (discharge) is vastly larger in supraglacial rivers, despite their small instantaneous storage volume (Smith and others, 2017).

4.4 Supraglacial lake drainage

We find that small (<0.25 km²), low-elevation supraglacial lakes are more likely to drain, while large (≥0.25 km²), high-elevation supraglacial lakes are more likely to persist throughout the 2015 melt season. Supraglacial lake count evolves seasonally from low to high elevations (×Fig. 6), as do surface meltwater area proportions and volumes (Figs 4a–d). Between 1000 and 1800 m, for example, supraglacial lake counts are higher earlier in the season than their areas (Fig. 6a), indicating large lakes often continue growing while small lakes drain or evolve into narrow channels. The concurrent increase of large lake area (≥0.25 km²) and count further supports this finding (Fig. 6b).

Fig. 6. Seasonal evolution of (a) all; and (b) large (≥0.25 km²) supraglacial lakes during the 2015 melt season. (c) Large (solid line) and all (dashed line) lake disappearance ratios decrease from low to high elevations, while mean lake areas (with one std dev. error bar) increase. Black solid line presents mean lake areas for the same study region during the 2003 melt season as reported in Sundal and others (2009).

From low to high elevations, the supraglacial lake disappearance ratio decreases while mean lake area increases, again indicating that supraglacial lakes at lower elevations are more likely to drain (Fig. 6c). Above 1400 m elevation, supraglacial lake areas are more variable than at lower elevations, suggesting that high-elevation lakes experience greater seasonal expansion. Moreover, 21% of these large lakes drain rapidly via hydrofracture, suggesting the presence of surface-to-bed meltwater connections at high elevations even during a colder-than-average melt season (Yang and Smith, 2016; Gagliardini and Werder, 2018).

4.5 Supraglacial hydrology and modeled surface runoff

Remotely sensed meltwater area and modeled surface runoff positively correlate for all three climate models. The mean correlation coefficient (R ²) between surface runoff and meltwater area proportion across all elevation intervals is 0.90, 0.84 and 0.79 for MERRA-2, MAR and RACMO, respectively. At high elevation intervals (1600–1800 and >1800 m), the corresponding R² values are 0.88, 0.74 and 0.65, respectively (×Fig. 7). Out of the three models, we find that MERRA-2 displays the highest correlation with observed meltwater area. However, due to the intermittent temporal sampling of Landsat 8 images, it is unknown whether time lags exist between peak surface runoff and maximum meltwater area proportion (Figs 4c–j).

Fig. 7. Correlations between MERRA-2 modeled surface runoff and remotely sensed surface meltwater area proportion. Larger circle sizes indicate later image acquisition dates. Except for very high (>1800 m) elevations, positive linear correlations are obtained at all elevations during 17 July to 25 August, indicating that greater (modeled) surface runoff is associated with greater surface meltwater area proportion.

MAR and RACMO runoff progressively declines for elevations >1400 m, whereas MERRA-2 simulates roughly similar runoff between 1000 and 1600 m (Figs 4e–j). MERRA-2 best simulates surface runoff during 8 and 24 July, when supraglacial hydrologic networks are well developed at 1400–1600 m elevation (Figs 4c, e). MERRA-2 cumulative runoff also best matches in situ ablation measurements from the PROMICE stations (Fig. 4l). The gradient of this positive linear relationship is steeper at higher elevations (Fig. 7), suggesting that high elevations are more sensitive to increasing meltwater runoff supply.

Surface meltwater storage accounts for 2.5–14.4% of modeled surface runoff over the 2015 melt season (Table 1), partitioned as 1.7–9.8% in lakes and 0.8–4.6% in rivers. A low percentage on 13 June (2.5%) is mainly attributed to the late start of the 2015 melt season (Fig. 2). While the percentages of lake/river storage (as compared to modeled surface runoff) show large variations at different elevation intervals (Fig. S4), our lake percentages are generally consistent with previously reported values (e.g. 5–13%, Leeson and others, 2012; Arnold and others, 2014; Fitzpatrick and others, 2014).

4.6 Supraglacial hydrology and ice flow

No direct correlation between surface meltwater area and storage and ice surface velocity is found (Fig. 4k). As the 2015 melt season progresses, slight decreases in ice velocity are observed, particularly at lower elevations (Fig. 4k). This is consistent with previous studies demonstrating reduced ice velocities in response to increased surface melt and development of an efficient subglacial hydrologic system (e.g. Bartholomew and others, 2011; Tedstone and others, 2015). Surface meltwater area proportion peaks on 1 July, after which supraglacial lakes continue expanding (Figs 4a–c), while ice flow decreases at 1200–1400 m (KAN_M station). On 24 July, a similar phenomenon occurs at a higher elevation (>1600 m, KAN_U station). This suggests that the well-developed supraglacial hydrologic networks observed at higher elevations drain meltwater into the ice sheet late in the melt season, causing ice flow velocities to accelerate, as occurs earlier in the melt season at lower elevations.

4.7 Supraglacial hydrology and proglacial river discharge

Modeled surface runoff is compared with remotely sensed surface meltwater area proportion and observed proglacial discharge in the Watson (Qinnguata Kuussua) River (×Fig. 8). Two major peaks (~1000 m³/s) in proglacial discharge occur around 8 and 24 July, a time of high surface meltwater area proportion. All three models simulate the 8 July discharge peak, with one-daytime lag between peak surface runoff and peak proglacial discharge. All three models also simulate a second runoff peak around 16–17 July, some 7–8 d before the observed proglacial discharge peak on 24 July. This suggests that the first proglacial discharge peak is driven mainly by low-elevation (<1600 m) meltwater contributions (Figs 4a–d), whereas the second proglacial peak includes high-elevation (>1600 m) meltwater contributions (Figs 4a–d). At higher elevations, a longer lag time is required for surface meltwater to reach the ice edge via subglacial pathways (Chandler and others, 2013; Andrews and others, 2014; van As and others, 2017).

Fig. 8. Daily in situ proglacial discharge in the Watson River at Kangerlussuaq (van As and others, 2017) and corresponding surface runoff as predicted by the MERRA-2, MAR 3.6 and RACMO 2.3 climate models (not including routing delays). Red crosses indicate surface meltwater area proportions derived from Landsat 8 satellite images. Observed proglacial discharge is most closely tracked by remotely sensed surface meltwater area proportions, and MERRA-2 runoff simulations based on simple correlation statistics.

MERRA-2 slightly underestimates observed peak discharge in both cases (by 15.0 and 25.1%, respectively, see ×Fig. 8), and is the only model to simulate a second peak comparable in magnitude to the first peak. MAR and RACMO both simulate smaller second peaks (21.0 and 16.2% lower than the first peak, respectively), and significantly overestimate both peak magnitudes (46.3 and 8.4% by MAR, and 86.2 and 46.4% by RACMO). This finding confirms previous reports of overestimated runoff in MAR and RACMO simulations (Rennermalm and others, 2013; Overeem and others, 2015; Smith and others, 2015, 2017; van As and others, 2017).

Among these three models, MERRA-2 runoff best matches observed proglacial discharge in the Watson River, with a Nash-Sutcliffe model efficiency (Nash and Sutcliffe, 1970) of E = 0.34. In contrast, E for MAR and RACMO are 0 and −1.27, respectively. MERRA-2 matches early-season proglacial discharge before July 19 particularly well (E = 0.88), significantly outperforming the other two models (E = −0.09 for MAR and −2.27 for RACMO). All three models simulate increasing proglacial discharge during July 19–24, just several days too early due to non-representation of routing delays. After July 24, MERRA-2 captures the overall decreasing trend of proglacial discharge but underestimates proglacial discharge by ~50% (Fig. 8).

5. Discussion

5.1 Supraglacial lakes and rivers are well developed at high elevations

This study provides enhanced detail to the temporal and spatial evolution of the southwest GrIS supraglacial hydrologic network and its variation with elevation (Fig. 4). In particular, our finding of well-developed supraglacial lakes and rivers above 1600 m signifies the presence of considerable surface runoff at high elevations. Indeed, we find that peak surface meltwater area proportion at high elevations (7.5%) is slightly larger than at low elevations (7.0–7.4%, Fig. 4d). Peak surface meltwater storage is also slightly greater at high elevations (115 mm) than low elevations (95–110 mm, Fig. 4d), perhaps due to the presence of large supraglacial lakes that persist overwinter (Koenig and others, 2015; Lampkin and others, 2020).

We also document the presence of significant surface-to-bed meltwater connections at high elevations (Figs 5, 6). Because high-elevation supraglacial catchments/depressions are typically large, they can drain significant volumes of meltwater into moulins despite having lower overall ablation rates (Yang and Smith, 2016). The present study firmly asserts the presence of high-elevation (>1600 m) moulins on the southwest GrIS, with 21% of high-elevation supraglacial lakes and 28% of rivers draining into them (Figs 3f, 5c). This affirms the substantial presence of surface-to-bed meltwater connections at high elevations even in a colder-than-average year (Yang and Smith, 2016; Gagliardini and Werder, 2018), contrary to high-elevation meltwater evacuation via very long supraglacial rivers (Poinar and others, 2015).

While the long-term ELA of our study area (~1550 m, van de Wal and others, 2012) is supported by remotely sensed end-of-season snowline elevation (averaging 1520 ± 113 m from 2001 to 2017, Ryan and others, 2019), recent field measurements suggest the ELA is rising and is now situated above 1700 m (Smeets and others, 2018). Our finding of well-developed supraglacial lakes and rivers at elevations >1600 m is therefore not particularly surprising. Future studies should further investigate the impact of these supraglacial hydrologic features on surface mass balance and ice flow at higher elevations in the ice-sheet interior (Poinar and others, 2015; Yang and Smith, 2016; Gagliardini and Werder, 2018).

5.2 Modeled surface runoff is correlated with satellite-observed meltwater area and volume

Positive linear relationships are found between modeled surface runoff and satellite-observed surface meltwater area proportions (Fig. 7). These relationships offer a new tool for understanding GrIS surface runoff where direct observations are difficult to obtain. In addition, we find 2.5–14.4% of surface runoff is stored in supraglacial lakes and rivers, a finding consistent with: (1) a model study for Paakitsoq (north of our study area), which reports values of 5–10% of runoff stored in supraglacial lakes in 2001, 2002 and 2005 (Arnold and others, 2014); (2) a value of 12% modeled for the entire southwest sector of the GrIS in 2003 (Leeson and others, 2012); and (3) remotely sensed values of 7–13% for the Watson river basin in 2010 and 2012 (Fitzpatrick and others, 2014). Supraglacial storage capacity is thus relatively small but stable during both cold (e.g. 2015) and warm (e.g. 2010 and 2012) years. This finding also signifies that supraglacial rivers remove >85% of meltwater runoff from the ice surface, consistent with Smith and others (2015). In this area of the GrIS, supraglacial rivers are therefore the dominant mechanism by which surface meltwater runoff enters the subsurface, not supraglacial lake drainage events.

Because supraglacial lakes and rivers are both fed by surface meltwater runoff (Banwell and others, 2012; Arnold and others, 2014; Chu, 2014; Fitzpatrick and others, 2014; Smith and others, 2017; Pitcher and Smith, 2019; Yang and others, 2019a), it is perhaps not surprising to find a positive linear relationship between surface meltwater area proportion and modeled surface runoff (Fig. 7). However, while numerous positive correlations between numbers/areas/volumes of supraglacial lakes and ice surface temperature/melt/runoff have been reported (McMillan and others, 2007; Sundal and others, 2009; Liang and others, 2012; Howat and others, 2013; Fitzpatrick and others, 2014; Langley and others, 2016; Kingslake and others, 2017; Macdonald and others, 2018; Williamson and others, 2018b; Yang and others, 2019a), supraglacial rivers are rarely considered.

Our results suggest that surface meltwater area proportion increases with modeled surface runoff (Fig. 7). This first empirical demonstration of a direct correlation between these two variables holds promise for calibration and validation of ice-sheet runoff models. It may be possible, for example, to parameterize ice surface drainage density and surface meltwater routing delays and peak discharges to moulins directly from modeled surface runoff. This knowledge would improve estimates of the diurnal timing of surface meltwater delivery to moulins (Arnold and others, 2014; Banwell and others, 2016), which modulates subglacial water pressures and impacts ice flow dynamics (Andrews and others, 2014, 2018; Gagliardini and Werder, 2018; Yang and others, 2020). This exciting potential application of empirical regression models relating observed meltwater area proportion to modeled surface runoff warrants further investigation.

One limitation of this study is the intermittent acquisition of clear-sky Landsat 8 images over our study area. Just a small number of observations (five each) were used to create the regression models presented in Figure 7. Therefore, these regression models are sensitive to the inclusion or removal of data points. Additionally, the evolution of supraglacial lakes and rivers is influenced by ice surface topography (Ádám and others, 2016; Crozier and others, 2018) and crevassing, which are not captured by climate model simulations. Therefore, we do not recommend using climate models to directly compute surface meltwater area proportions. Instead, we suggest that the primary utility of such regressions (e.g. Fig. 7) is to parameterize ice-sheet runoff models (Smith and others, 2017; Yang and others, 2018, 2020).

5.3 Modeled and observed ice-sheet surface runoff

Of the three climate models examined here, MERRA-2 best matches both satellite-observed surface meltwater area and proglacial discharge (Figs 4, 8). This result is consistent with previous reports of reasonable agreement between MERRA-2 and supraglacial river discharge (Smith and others, 2017) and lake infilling rate (Yang and others, 2019a). MERRA-2 replicates proglacial discharge in the Watson River more closely than MAR and RACMO, particularly in the early melt season (Fig. 7). Moreover, MERRA-2 simulates more surface runoff at higher elevations (>1400 m) than lower elevations (<1400 m), which is consistent with our observations of a well-developed supraglacial hydrologic network there (Fig. 4e). The reason for this is unclear and lies beyond the scope of this study, but one possible explanation might be that MERRA-2 uses spatially variable albedo and fractional snow cover (Reichle and others, 2017), possibly better representing GrIS surface energy balance.

Proglacial discharge is commonly compared with modeled surface runoff (Overeem and others, 2015; Smith and others, 2015; van As and others, 2017) but large discrepancies are common due to subglacial delays and other processes (Rennermalm and others, 2013; Overeem and others, 2015; Smith and others, 2015; van As and others, 2017). While our focus here is on meltwater area and storage volume, not discharge, we suggest that remotely sensed meltwater area and storage may correlate with proglacial discharge as much or more than modeled runoff (Fig. 8). Therefore, observations of time-evolving supraglacial hydrologic network should be included in models that simulate proglacial discharge. The contributions of surface meltwater at high elevations (>1600 m), in particular, warrant further study.

5.4 Issues of temporal and spatial scale

This study used a seasonal time series of seven Landsat 8 images to track the evolution of surface meltwater on the southwest GrIS during a colder-than-average melt season. However, seven satellite images are likely intermittent to comprehensively reveal the temporal evolution of a highly dynamic supraglacial hydrologic network (Liang and others, 2012; Chu, 2014; Yang and others, 2017). The 30 m resolution of Landsat 8 also overestimates the widths of large supraglacial rivers due to mixed pixels (Yang and Smith, 2016), while failing to detect numerous small supraglacial streams (Figs S2, S3). We suggest that these mapping uncertainties influence absolute river area proportions, but not their relative pattern (e.g. high-elevation vs. low-elevation comparisons). As the number of higher-resolution satellites increase (e.g. Sentinel-2A/B twin satellites with 10 m spatial resolution; Planet CubeSats with 3–5 m spatial resolution), the prospects for investigating the temporal evolution of supraglacial hydrologic network at finer temporal and spatial scales will improve (Williamson and others, 2018a; Yang and others, 2019b).

6. Conclusions

Using multi-temporal Landsat 8 satellite imagery, we quantified the evolving area and volumetric storage of supraglacial meltwater lakes and rivers on the southwest GrIS during summer 2015, a colder-than-average year. Supraglacial rivers dominate total surface meltwater area and discharge, whereas supraglacial lakes dominate total surface meltwater volume. Owing to their greater storage capacity, lakes require more time to fill and attain maximum area 7–10 d later than supraglacial rivers, confirming that rivers respond rapidly to surface runoff variations. Some supraglacial rivers and lakes drain into moulins at high elevations (>1600 m), signifying the presence of high-elevation surface-to-bed meltwater connections even in a relatively cold year.

Remote sensing of time-varying supraglacial hydrology offers a unique opportunity to validate and improve climate model simulations of ice-sheet runoff. In particular, we present a first quantitative relationship that correlates modeled surface runoff with satellite-observed surface meltwater area proportion, signifying that surface melting intensity drives the formation of supraglacial hydrologic networks. Integrating satellite observations with climate models should therefore improve the parameterization of ice-sheet runoff models, and mechanistic understanding of how surface mass balance drives ice-sheet runoff, moulin inputs and associated subglacial water pressures.