Articles | Volume 17, issue 11
Research article
09 Nov 2023
Research article |  | 09 Nov 2023

Seasonal evolution of the supraglacial drainage network at Humboldt Glacier, northern Greenland, between 2016 and 2020

Lauren D. Rawlins, David M. Rippin, Andrew J. Sole, Stephen J. Livingstone, and Kang Yang

Supraglacial rivers and lakes are important for the routing and storage of surface meltwater during the summer melt season across the Greenland Ice Sheet (GrIS) but remain poorly mapped and quantified across the northern part of the ice sheet, which is rapidly losing mass. Here we produce, for the first time, a high-resolution record of the supraglacial drainage network (including both rivers and lakes) and its seasonal behaviour at Humboldt Glacier, a wide-outlet glacier draining a large melt-prone hydrologic catchment (13 488 km2), spanning the period 2016 to 2020 using 10 m spatial resolution Sentinel-2 imagery. Our results reveal a perennially extensive yet interannually variable supraglacial network extending from an elevation of 200 m a.s.l. to a maximum of  1440 m a.s.l. recorded in 2020, with limited development of the network observed in the low-melt years of 2017 and 2018. The supraglacial drainage network is shown to cover an area ranging between 966 km2 (2018) and 1566 km2 (2019) at its maximum seasonal extent, with spatial coverage of up to 2685 km2 recorded during the early phases of the melt season when a slush zone is most prominent. Up-glacier expansion and the development of an efficient supraglacial drainage network as surface runoff increases and the snowline retreats is clearly visible. Preconditioning of the ice surface following a high-melt year is also observed, with an extreme and long-lasting 2019 melt season and over-winter persistence of liquid lakes, followed by low snow accumulation the following spring, culminating in earlier widespread exposure of the supraglacial drainage network in 2020 compared to other years. This preconditioning is predicted to become more common with persistent warmer years into the future. Overall, this study provides evidence of a persistent, yet dynamic, supraglacial drainage network at this prominent northern GrIS outlet glacier and advances our understanding of such hydrologic processes, particularly under ongoing climatic warming and enhanced runoff.

1 Introduction

The Greenland Ice Sheet (GrIS) has experienced significant mass loss throughout the 21st century and currently represents the largest single cryospheric component of global sea level rise, contributing an estimated 10.6 ± 0.9 mm since 1992 (The IMBIE Team, 2020). Over the last 2 decades, GrIS mass loss has become increasingly dominated by surface mass balance (SMB) processes, accounting for 60 % of ice loss annually since 1991, with the remainder attributed to dynamical mass losses from marine-terminating glaciers along the ice sheet periphery (van den Broeke et al., 2016). Such SMB losses are being increasingly revealed by the magnitude and spatial extent of seasonal surface melting and runoff (Trusel et al., 2018), attributed to climate-driven atmospheric warming (Hanna et al., 2012, 2021), summertime atmospheric circulatory behaviour (i.e. Greenland Blocking Index; Hanna et al., 2012, 2021; McLeod et al., 2016; van den Broeke, 2017), and the ongoing expansion (Noël et al., 2019) and darkening (Tedesco et al., 2016; Ryan et al., 2018, 2019; Riihelä et al., 2019) of the bare ice zone. Between 2011 and 2020, runoff was 21 % higher than any of the preceding 3 decades (Slater et al., 2021).

Surface runoff is transported by an expansive and complex supraglacial drainage system that is activated during the summer season (Pitcher and Smith, 2019). This drainage system, made up of an ephemeral network of interconnected supraglacial rivers and lakes, transports and stores large volumes of surface meltwater on the ablating ice surface (Rippin and Rawlins, 2021). Such runoff can become intercepted by crevasses and moulins that provide connections to the ice sheet bed where the timing and delivery of such water has been shown to affect ice velocity (Zwally et al., 2002; Bartholomew et al., 2010, 2012; Hoffman et al., 2011; Sole et al., 2011; Andrews et al., 2014; Nienow et al., 2017). In particular, meltwater delivery into an inefficient subglacial configuration, such as linked cavities (Kamb, 1987), which typically occurs during the early period of the melt season, can temporarily overwhelm the subglacial hydrologic system, increasing water pressure and enhancing subsequent sliding (Andrews et al., 2014; Davison et al., 2019). In some regions where moulins and crevasses are absent, supraglacial rivers can extend undisturbed for tens of kilometres across the bare ice surface, flowing directly into the proglacial zone (Yang et al., 2019a; Li et al., 2022). Ultimately, much of this meltwater will end up in the ocean, contributing directly to global sea level rise (Pitcher and Smith, 2019).

Whilst many remote sensing studies have examined components of the supraglacial drainage network in-depth across the largest melt producing western and southwestern sectors of the GrIS (Smith et al., 2015; Gleason et al., 2016, 2021; Yang et al., 2021), it is only recently that other significant ice sheet sectors have begun to be mapped (Gledhill and Williamson, 2018; Macdonald et al., 2018; Yang et al., 2019a; Schröder et al., 2020; Turton et al., 2021; Lu et al., 2021; Boghosian et al., 2023). Focus has only recently shifted to the rapidly changing northern regions of the GrIS, with evidence of inland expansion of supraglacial lakes observed in northeastern Greenland (Turton et al., 2021) that aligns with climate model projections (Leeson et al., 2015) and the existence of a widespread supraglacial network (Lu et al., 2021). This study utilises Sentinel-2 satellite imagery to map the supraglacial drainage network, including both rivers and lakes, on a major northern outlet glacier of the GrIS – Humboldt Glacier (7923.86 N, 6420.60 W), here denoted as HG – to examine its seasonal behaviour at high spatial (10 m) and temporal resolution over 5 consecutive melt years (2016–2020).

2 Study location

The drainage basins of outlet glaciers in the northern sector of the GrIS comprise  14 % of the total ice sheet area, with 82 % of the northern sector predominately drained by 12 marine-terminating glaciers which together hold a sea level equivalent of 93 cm (Mouginot et al., 2019). Since 1990, this sector has experienced some of the most pronounced changes in surface melt and runoff, attributed to the rapid expansion of the ablation (46 %) and bare ice (33 %) zone at rates twice as fast as in southern Greenland, with this trend expected to continue with ongoing climatic warming (Noël et al., 2019). Of the northern outlet glaciers that drain the GrIS, HG, also known as Sermersuaq Glacier, is the widest marine-terminating outlet glacier ( 91 km wide) in Greenland and is responsible for draining  5 % of the ice sheet alone northwestward into the Kane Basin (Hill et al., 2017; Rignot and Kanagaratnam, 2006; Rignot et al., 2021; Fig. 1). Since the late 1990s, HG has experienced rapid rates of retreat ( 162 m yr−1) attributed to increases in mean summer air temperatures and sea ice decline (Carr et al., 2015). Holding an ice volume equivalent of 19 cm of sea level rise, HG is the fourth-largest Greenland glacial contributor to sea level rise (Rignot et al., 2021), having lost 161 Gt since 1972 and 311 km2 of its area between 2000 and 2010 (Box and Decker, 2011).

Figure 1(a) Study location of HG, northern Greenland, and its highlighted drainage basin (shaded red). The dashed box shows the inset for other figure boxes (b–d). (b) The four Sentinel-2 tiles used for extraction of the supraglacial drainage network across the study region. (c) True colour (R: band 4; G: band 3; B: band 2) Sentinel-2 image of HG acquired from 25 July 2020 courtesy of the Copernicus Open Access Hub (, last access: 15 December 2022). HG denoted drainage basin and 200 m contour lines derived from ArcticDEM (10 m) are shown. The 1500 m a.s.l. contour denotes the maximum melt extent. The shaded section up to 200 m a.s.l. also shows the heavily crevassed zone that exists within the northern sector of the terminus. (d) Bed topography of HG and the surrounding area from BedMachine version 4 (Morlighem et al., 2021), ice velocity contours via NASAs MEaSUREs ITS_LIVE project (Gardner et al., 2019), and the division of Humboldt's north and south sectors (dashed line) as per Carr et al. (2015).

Until recently, few glaciological studies had focussed on HG (Joughin et al., 1996, 1999; Carr et al., 2015; Livingstone et al., 2017; Hill et al., 2017, 2018; Mouginot et al., 2019; Gray, 2021; Rignot et al., 2001, 2021, Hillebrand et al., 2022). Studies that have examined HG identified a distinctive ice velocity divide between the northern and southern sectors (Rignot et al., 2001, 2021; Carr et al., 2015; Fig. 1d); the northern sector has up to 4 times faster ice flow than the south. In terms of surface hydrology, several studies have noted the presence of supraglacial lakes (SGLs; Joughin et al., 1996; Selmes et al., 2011; Carr et al., 2015); however, none to-date have examined the overall drainage system, including both rivers and lakes, in detail.

3 Data and methodology

3.1 Data sources

Earth-observing satellites enable the study of supraglacial drainage features with broad spatial and temporal coverages (Rennermalm et al., 2013; Yang and Smith, 2012; Chu, 2014). Over the last 4 decades, the Landsat Program has provided a wealth of remotely sensed data for the mapping and quantification of a number of supraglacial features such as SGLs (Lampkin and VanderBerg, 2011; Banwell et al., 2014; Pope et al., 2016; Williamson et al., 2017; Gledhill and Williamson, 2018; Williamson et al., 2018; Yang et al., 2019b; Dell et al., 2022; Otto et al., 2022), as well as for exploring the generalised configuration (i.e. main river stems) of the supraglacial drainage system (Lampkin and VanderBerg, 2014; Yang et al., 2021). However, its spatial resolution in the visible spectrum (30 m) precludes the reliable delineation of numerous smaller supraglacial rivers (Yang et al., 2019a). This has resulted in these complex networks being unmapped and underrepresented (Chu, 2014). The application of the Multispectral Instrument (MSI) on Sentinel-2 satellites (Sentinel-2A and -2B), which launched in 2015 and 2017, respectively, offers a higher-resolution (10 m) perspective of such systems (Yang et al., 2019a; Lu et al., 2020, 2021). Sentinel-2 imagery enables the detection and delineation of both wide, main-stem river channels, which have high stream orders and are perennially reoccupied (Pitcher and Smith, 2019), as well as narrower (one pixel, or 10 m) tributary-style channels that are of lower order and shallower in depth (Smith et al., 2015; Fig. 2). Sentinel-2 imagery has also been shown to better resolve supraglacial networks in general for mapping purposes at a glacier-wide scale, particularly in terms of river continuity (Yang et al., 2019a), which contributes to its preferred use in this study.

Figure 2Example of the supraglacial drainage features found across the study region of HG, represented in panel (a), from a Sentinel-2 image taken on 23 June 2020 (RGB). (b) A supraglacial lake with a central ice lid feeding an outlet supraglacial river, with evidence of river advection transverse to ice flow. (c) Braided supraglacial rivers flowing between supraglacial lakes, known as “connector” lakes. (d) A large canyonised supraglacial river with a 90 bend terminating abruptly in a moulin. (e) Deep, canyonised supraglacial rivers flowing off the ice edge. (f) Narrow supraglacial rivers with rill features seen in the higher-elevation regions of HG, flowing and coalescing into braided tributary rivers. (g) The slush zone, with dense areas of saturated snow (light blue) culminating in rill-like channels. The Sentinel-2 image is courtesy of the Copernicus Open Access Hub (

For the years 2016 to 2020, a total of 176 Sentinel-2 Level-1C (orthorectified top-of-atmosphere reflectance) images with sub-pixel multispectral registration (Baillarin et al., 2012) were acquired over HG (Fig. 1b; Table S1 in the Supplement) obtained from ESA's Scientific Data Hub (, last access: 15 December 2022). These images covered the entirety of the study area on 44 d between the months of May and September across the study period, equating to 1–2 images per month, allowing us to gain a full melt season perspective of supraglacial drainage evolution for the HG drainage catchment. For scenes with cloud cover below a 20 % threshold, cloud coverage was typically restricted to the Kane Basin waterway (Fig. 1c) or ice interior locations beyond the melt extent, and thus they did not pose any significant problems for river and lake mapping.

The HG drainage catchment was generated using the ArcticDEM 10 m mosaic product obtained from the Polar Geospatial Centre (, last access: 15 December 2022) and delineated following the method of Karlstrom and Yang (2016). Elevation contours at 100 m intervals were defined across the HG basin up to 1500 m a.s.l. (maximum limit of the melt-prone zone), equating to a size of 13 488 km2. Daily surface meltwater production and runoff for the study area were extracted from the Modèle Atmosphérique Régional (MAR) regional climate model (RCM) v3.11 (available at, last access: 15 December 2022; Fettweis et al., 2017, 2020, 2021). MAR is among the best RCM to cover the GrIS as it explicitly models important polar processes (i.e. SMB) forced with ERA5 reanalysis data and has been extensively evaluated against in situ automatic weather station and satellite data (for a detailed description of MAR v3.11, see Amory et al., 2021; Fettweis et al., 2020, 2021). MAR has now been widely used in other GrIS-based supraglacial hydrologic studies (Smith et al., 2017; Yang et al., 2019b; Lu et al., 2020) for quantifying the relationship between modelled runoff and satellite-derived meltwater metrics for RCM accuracy validation. The version MAR 3.11 used in this study was run at high spatial (6 km) and temporal (daily) resolution to generate estimates of daily meltwater production and runoff (R) in millimetres of water equivalent per day (mm w.e. d−1) within each grid cell to assess the spatial and temporal distribution of meltwater against mapped supraglacial rivers and lakes, similar to other studies (Yang et al., 2019b; Lu et al., 2020, 2021; Yang et al., 2021). MAR grid cells were also sampled at each 100 m interval to assess elevational gradients in both mapped drainage and runoff. A MAR uncertainty value of ±15 % was also calculated (Fettweis et al., 2020).

3.2 Supraglacial river and lake extraction

To effectively delineate supraglacial rivers from remotely sensed imagery, an automatic linear enhancement method developed by Yang et al. (2015) was used, which characterises supraglacial rivers according to their Gaussian-like brightness cross sections and longitudinal open-channel morphology. Firstly, a normalised difference water index was performed following McFeeters (1996) to differentiate active surface meltwater from the background ice and snow (Lu et al., 2020; Li et al., 2022). This equation (Eq. 1) utilises Band 3 (“Green”) and Band 8 (“NIR”) from Sentinel-2 imagery as follows:

(1) NDWI = ( Green - NIR ) ( Green + NIR ) .

An ice-derived spectral index, NDWIice, has been widely applied to supraglacial mapping studies in recent years (Yang and Smith, 2012; Moussavi et al., 2016; Williamson et al., 2018; Yang et al., 2021), as it produces fewer false classifications of blue ice and slush areas. We use NDWI (McFeeters, 1996) here because of its successful implementation in other studies (Stokes et al., 2019; Lu et al., 2020, 2021; Corr et al., 2022) and its ability to map all active surface melt of interest, including slush zones. Additionally, preliminary testing of a small sample area found using NDWI (McFeeters, 1996) was able to map 16.3 % more regularly connected supraglacial river channels compared to those mapped by the alternate NDWIice (Supplement Fig. S1). Next, an ice mask was applied, created from manually digitising the HG terminus from the latest, end-of-season image (late August–early September) in each of the study years to extract ice-only regions and remove surrounding land, rocky outcrops and the ocean of the Kane Basin. A separate crevasse mask was delineated from manually identifying the heavily crevassed zone known to extend up to 25 km from the northern sector of the terminus, with particular prevalence across the 7 km floating section (Carr et al., 2015). This mask was applied, and this section of the terminus removed to avoid the erroneous delineation of crevasses and crevasse shadows in this area during image processing and to reduce the effects of drainage overestimation in these ice marginal regions in further calculations (Ignéczi et al., 2018). A global NDWI threshold of 0.4 (t0.4) was then applied to the masked NDWI image as it effectively captures SGL boundaries and wide, main-stem river segments (Lu et al., 2021). To aid in delineating narrower river segments and obtain a complete and continual supraglacial network, the automatic river detection algorithm for linear enhancement (Yang et al., 2015) was then applied. This involved the removal of the low-frequency image background and high-frequency image noise using a band-pass filter ramped between 1/200 and 1/40 m−1 (Yang et al., 2019a), before the application of a Gabor filter to enhance the cross sections of smaller rivers (< 2 pixel width). This was followed by a parsimonious path opening (PPO) operator, which is a flexible mathematical morphological operator, to stabilise linear brightness across river lengths and preserve connectivity, with a minimum path length of 20 pixels.

After Gabor–PPO filtering, the supraglacial river network becomes easier to differentiate and delineate from the surrounding icy background. A global pixel brightness threshold of 5 (out of 255) denoted t5 was used to extract supraglacial rivers of varying widths from Gabor–PPO opened filtered images (Lu et al., 2020), i.e. from tributary-style rivers to main-stem, large river channels. Hydrologically connected slush zones, which is where the pore space of snow becomes entirely water saturated when temperatures permit melting (above freezing, 0 C) and forms expansive fields of ponded surface water, were retained as they play an important role in the initial mobilisation of melt as slush flows or within rill-type channels and the inland expansion of the melt-prone zone as summer progresses (Holmes, 1955; Marston, 1983; Cuffey and Paterson, 2010; Chu, 2014; Rippin and Rawlins, 2021). Slush zones can be spectrally distinguished in true-colour satellite images as dense, light blue patches on the surface as snow becomes water saturated (Fig. 2g), and they can be partially distinguished in NDWI images as bright, dense features similar to that of individual linear river channels (Lu et al., 2021). For mapping conducted in this study, slush zones are mapped inclusively within the drainage network and not treated independently due to their overall hydrologic importance and spectral similarity to other hydrological components. Whilst dynamic thresholding techniques have been used in other studies, particularly in the identification and mapping of supraglacial lakes across independent dates and/or years (Selmes et al., 2011; Williamson et al., 2017), applying these two separate thresholds (t0.4 and t5) to all images across this study period is reasonable for exploring the seasonal behaviour over the available dates used. Finally, these masks were vectorised into separate river channel polylines and lake polygons for analysis, with a size threshold of 0.1 km2 applied to SGLs. To assess the ability of the automated mapping technique to capture supraglacial rivers at HG, a small area was manually digitised (Fig. 3). The overall spatial pattern extracted was similar, and the difference between manual vs. automated river area was 5.4 %, with automated rivers tending to have several gaps in the network compared to manual results. A similar finding was also produced by Lu et al. (2020), who found a difference of 13.6 % between manual vs. automated supraglacial rivers also using Sentinel-2 imagery.

Figure 3(a) RGB image of a sample area of HG with both main-stem and tributary supraglacial river channels present. (b) A performance accuracy assessment comparing the automatic river detection algorithm (blue) used in this study with manually digitised networks (yellow). Overall, automated rivers networks were 5.4 % shorter than those manually derived. The Sentinel-2 image in (a) is courtesy of the Copernicus Open Access Hub (

3.3 Supraglacial river and lake quantification

To characterise the extracted supraglacial drainage system, metrics were calculated for both rivers and lakes. These metrics are summarised in Table S2 and include meltwater area (km2), meltwater area fraction (MF, %), drainage density (Dd), the number of supraglacial lakes (Ln) and supraglacial lake area (La). MF is defined as the percentage total meltwater area across the drainage catchment below a conservative upper melt limit of 1500 m a.s.l. for each date mapped, which is then also further divided into separate feature ratios including river area fraction (RF) and lake area fraction (LF). MF was also calculated across 100 m elevation contours from above the heavily crevassed zone at 200 m a.s.l., to the maximum melt extent at 1500 m a.s.l. To explore the relationship between remotely mapped drainage and modelled RCM MAR runoff (R), a Spearman's rank correlation (rs) was performed and linear regression analysis undertaken with subsequent R2, rs and P values reported.

4 Results

4.1 Spatial characteristics of the supraglacial drainage network

Supraglacial rivers and lakes were mapped from a total of 44 dates across the lower, melt-prone 13 488 km2 HG drainage basin from the melt seasons of 2016 to 2020. The mapped supraglacial drainage network across HG is shown extend up to 1440 m a.s.l., with well-developed, main-stem river channels occurring up to 1000 m a.s.l., which we characterise as the persistent zone, and an ephemeral network of tributary-style rivers and slush zones extending beyond 1000 m a.s.l. in a transient zone up to maximum extent (Fig. 4a). Active supraglacial rivers and lakes form progressively upglacier from low elevations (200 m a.s.l.) to a maximum of 1440 m a.s.l. as the melt seasons progress across the study period, with interannual variability observed.

Figure 4(a) Map showing the recurrence frequency of the supraglacial drainage network across the study period (2016 to 2020) at HG. Dark blue shades denote a higher frequency of occurrence, which is highly prominent in the persistent zone below 1000 m a.s.l., with rivers and lakes typically reforming in the same location each year. Red arrows denote the location and direction of two distinct parallel drainage structures. The black arrow represents ice flow direction. The background Sentinel-2 image is provided courtesy of the Copernicus Open Access Hub ( (b) Close-up panel showing highly persistent rivers in the persistent zone with a prominent 90 angle in their channel form. The red arrow denotes some channel reconfiguration due to ice advection. (c) Close-up panel showing the dendritic supraglacial drainage structure typical in the northern sector of HG. (d) Close-up panel of the transient zone (> 1000 m a.s.l.), where yearly river persistence is lower and characterised by lengthy tributary river channels. (e) Close-up panel showing the more parallel form of supraglacial drainage structure in the southern sector of HG, with red arrows denoting some channel advection with ice motion.

Within the lower-elevation regions (< 400 m a.s.l.) of HG, collectively the supraglacial drainage network is largely fragmented, with many short (< 750 m long), supraglacial rivers observed alongside small SGLs with an average size of 0.23 km2, which is the smallest observed across all elevation bands. At greater elevations, beyond 400 m a.s.l., we observe large, main-stem supraglacial rivers, some with incised-canyon features (Fig. 2b), interconnected with increasingly larger SGLs. At elevations > 1000 m a.s.l., average SGL size is 0.41 km2, with a maximum SGL size of 2.08 km2. SGLs and rivers parallel to ice flow tend to be highly persistent year-on-year across the study period. In Fig. 4b, we also see some evidence of a potential main-river reconfigurations, with the northwestward advection of a river channel that runs transverse to ice flow. In the upper parts of the catchment, the supraglacial drainage network becomes increasingly dense, especially between 800–1000 m a.s.l., where we see a 120 % increase in average meltwater area (94.5 km2) compared to 200–400 m a.s.l. (42.9 km2), a 79 % increase compared to 400–600 m a.s.l. (54.2 km2) and a 10 % increase compared to 600–800 m a.s.l. Not only are persistent main-stem rivers still present up to 1000 m a.s.l., but an extensively connected tributary river system is also observed within this persistent zone.

The drainage network up to 1000 m a.s.l., in particular within the central and northern reaches of HG, is dendritic in nature (Fig. 4c). This type of drainage, however, is not uniform across HG, with the supraglacial network to the south exhibiting a more sub-parallel drainage style (Fig. 4e), with this configuration extending beyond the persistent zone and into the more transient zone with slush zone development and inland evolution a key part of the temporal aspect of the network at higher elevations (Fig. 4a). Even though density remains high until 1300 m a.s.l., vastly transient tributary rivers and slush zones dominate, feeding meltwater from headwater regions downstream. In terms of SGLs, the maximum recorded elevation was 1346 m a.s.l. both in 2019 and 2020. Fewer SGLs are observed in the transient zone, accounting for 16 % of total SGL area across the study period. However, average SGL area increases with elevation, with lakes in this zone being 54 % larger (0.4 km2) on average than those found below 1000 m a.s.l. (0.26 km2). The single largest SGL size of 2.08 km2 was recorded in 2016 at 1150 m a.s.l.

A key feature that is particularly prominent in the supraglacial drainage network is the presence of two parallel lines that track across glacier from the southwest to the northeast (Fig. 4a). Many supraglacial rivers and lakes are aligned along these two features, which appear as depressions in the ice surface, with the abrupt termination of many rivers indicating meltwater capture via moulins, indicating a strong structural element influencing drainage configuration.

4.2 Temporal evolution of the supraglacial drainage network

Typically, the supraglacial drainage network becomes active in early June with the onset of melt production and runoff in the region, with only a small number of large-stem supraglacial rivers becoming active and subsequently recorded (MF < 3.2 %) during this time within the mapped elevation bands of 200–950 m a.s.l. By late June a widespread (500 to 1150 m a.s.l.) slush zone develops and advances upglacier as the melt season progresses runoff increases, with MF ranging between 11.4 % and 19.9 %. As bare ice is exposed below this slush zone, the drainage system becomes increasingly channelised (Fig. 5). The formation of the slush zone at the end of June typically coincides with maximum melt storage in SGLs (both numbers and size of lakes). The network with the largest expanse of slush zone and number of SGLs in this study was observed on 30 June 2019 (Fig. 5), with 2685 km2 (19.9 %) of the HG ice surface comprised of a hydrologically connected, unchannelised system and 111 SGLs recorded (total area 27.4 km2). As the season progresses, the slush zone shifts upglacier while reducing in size, with average June MF decreasing by up to 33 % before stagnating at a maximum inland extent, ranging between 1050 m a.s.l. (2018) and 1440 m a.s.l. (2020) across the study period. At this elevation, the slush zone operates thereafter as headwaters, feeding the complex, transient, tributary systems below and further supplying the larger, well-defined supraglacial rivers at lower elevations (< 1000 m a.s.l.).

Figure 5Satellite-derived mapping of the temporal evolution of the supraglacial drainage network, including rivers and lakes, across the HG drainage basin from 2 melt years during the study period, with 2018 showing the typical behaviour of a low-melt year and 2019 the behaviour during a high-melt year. The background image is a Sentinel-2 image courtesy of the Copernicus Open Access Hub ( For all mapped study dates, please refer to Fig. S3.

Towards the end of the melt season, despite melt and runoff cessation, the supraglacial drainage network remains (Fig. 5).

The interannual variability in seasonal behaviour of the supraglacial drainage network between 2016 and 2020 (Fig. 5) corresponds to the length and intensity of the melt season. Drainage within the melt seasons of 2016, 2019 and 2020 follow a similar pattern characterised by a rapid increase and peak in MF in late June, yielding values of 11.4 % (R= 7.1 mm d−1), 19.9 % (R= 19.4 mm d−1) and 12.1 % (R= 19.6 mm d−1), respectively (Fig. 6), concurrent with early season melt production and runoff. These high MF values are largely associated with widespread slush zone initiation, with a subsequent peak in MF increasing the drainage network area by 267 % in 2016 (28 June) and 322 % in 2019 (30 June), with its spatial extent observed below 1000 m a.s.l. at this time (Fig. S3). The beginning of the melt season in 2020 is an exception, with an anomalously high MF (11.6 %) recorded on 15 June reaching 1100 m a.s.l., as well as a high number (57) and cumulative area of SGLs (11.2 km2; Fig. 6). This high is followed by a 66 % reduction in MF by 23 June, after which the network behaviour is similar to that of other seasons, with a subsequent MF increase of 211 % to its peak at the end of June (28 June; Fig. 6). Despite high rates of melt production and runoff throughout July across these three melt seasons, MF plateaus and reduces, reaching a steady state of between 7.6 % and 13.9 %. The number of SGLs also reduces on average between 20 % and 27 % throughout July, with cumulative SGL area reducing between 9 % and 38 %.

Figure 6The satellite-derived water metric meltwater area fraction (MF, %) for each mapped date across the study period (2016–2020) alongside melt production and runoff values (mm w.e d−1) derived from MAR v3.11 for the HG catchment. The ±15 % uncertainty envelopes are provided for MAR-derived values (Fettweis et al., 2020).


By the end of July, the supraglacial drainage network consistently extends to 75–80 km inland (1440 m a.s.l.) at its maximum areal extent. As melt and runoff reaches declines into August, the drainage network reduces between 53 % (2016) and 9 % (2020). In 2019, persistent high rates of melt production and runoff result in persistence of the drainage networks maximum inland extent and even a late-season increase of 54 % in MF from the 5 August to 13 August. The year 2019 is registered as an exceptional year at HG, with an average meltwater area up to 75 % greater in June compared to 2016, between 25 % and 92 % greater in July compared to 2016 and 2020, and up to 3 times greater (300 %) in August compared to both 2016 and 2020. Additionally, 129 % and 86 % more lakes are recorded compared to 2016 and 2020, respectively, with SGLs persisting much later into the melt season (mid-August) than the other years.

The supraglacial drainage network in 2017 and 2018 behaves quite differently to the other melt years (Fig. S4). Until late July, MF remains low (< 7 %; Fig. 6) with runoff predominately occurring within the lower 700 m of the HG basin (Fig. S3). The number and area of SGLs is also low during this period, with an average of 15 lakes observed in June and July across 2017 and 2018 with a cumulative average area of 21.5 and 31.1 km2, respectively; 280 % lower than average lake counts and 218 % lower than cumulative average area for the same monthly periods across 2016, 2019 and 2020. Unlike the other study years, peak MF does not occur until much later in the melt season, reaching 10.2 % on 22 August 2017 and 7.2 % on 20 August 2018 with a 690 km2 (197 %) and 521 km2 (135 %) increase in meltwater area to its prior mapped date, respectively; concurrent with peak runoff and the widespread occurrence of a late slush zone. Similarly, SGL numbers and area also peak on these dates, with 54 (15.8 km2) and 74 (20 km2) lakes recorded, comprising 32 % and 39 % of total lakes recorded per year, respectively. Examination of the maximum extent of the supraglacial drainage network is also revealed to be limited, with a maximum elevation of  1150 m in 2017 and  1050 m in 2018, equivalent to between 48 and 51 km inland.

When comparing what could be considered high-melt vs. low-melt year seasonal patterns, the supraglacial drainage network is twice as large during high-melt years (2016, 2019 and 2020) and is seen to form at elevations  300 m higher than in low-melt years (2017, 2018; Figs. S2 and S3). In terms of the average number of SGLs, there is almost double the number of SGLs observed (92 %) in the high-melt years compared to low-melt years, with average SGL area also being 111 % higher, showing large year-on-year variability in this system.

The contribution of rivers and SGLs to the supraglacial drainage network follows the same general seasonal trend as MF. RF and LF (%) peak at the end of June in the high-melt years of 2016, 2019 and 2020 and in mid-August for the low-melt years of 2017 and 2018 (Fig. 7), i.e. when hydrologically connected slush zones are most prominent and SGL numbers are shown to peak. Overall, the network is largely composed of supraglacial rivers, with RF accounting for an average of 6.2 % ± 5.6 % (mean ± SD) of the HG catchment across the study dates and is higher during high-melt years (7.8 % ± 6.0 %) compared to low-melt years (3.5 % ± 3.4 %). The highest RF recorded occurred on the 7 July 2019 and measured 24.5 %; the highest RF measured in a supraglacial hydrologic study. For SGLs, average LF was much lower, accounting for an average of 0.07 % ± 0.06 % of the melt-prone area of the catchment, with slightly higher LF in high-melt years (0.08 % ± 0.07 %) compared to low-melt years (0.05 % ± 0.04 %). In terms of the overall contribution of supraglacial rivers and SGLs to peak meltwater area fraction, 98 % of the network is dominated by supraglacial rivers across all study years, with SGLs playing a less dominant role in HG's drainage network. In comparison, rivers contribute  62 % to the drainage network in southwestern Greenland, with SGLs contributing  38 % (Yang et al., 2021).

Figure 7River area fraction (RF, %) and lake area fraction (LF, %) across the melt-prone zone of HG for each mapped date across the study period. On average, 98 % of the supraglacial drainage network across the study period at HG is comprised of RF.


4.3 Supraglacial hydrology and MAR runoff

A strong positive linear relationship was identified between satellite-derived MF and regional climate model MAR surface R for the HG catchment across the study period 2016 to 2020 (Fig. 8; R2= 0.77, Rs= 0.91, p=< 0.001) up until peak MF values and rapid surface runoff decline at the end of the melt season. Both MF and R increased concurrently each year as the melt season progressed, with peak runoff often coinciding with the existence of expansive slush fields across the upper part of the catchment. Runoff remained high until maximum extent occurred, particularly for high-melt years: 2016 (2 July, R= 7.5 mm d−1), 2019 (13 August, R= 8.2 mm d−1) and 2020 (30 July, R= 11.4 mm d−1). For low-melt years, runoff remained relatively high until early August (6 August 2017, R= 4.0 mm d−1; 2 August 2018, R= 7.1 mm d−1), with maximum extent occurring within 2 weeks (the next mapped date). This relationship between MF and R shows the reliability of simulated variations in seasonal surface meltwater runoff in capturing the behaviour of the supraglacial drainage network via satellite-derived water metrics, particularly during high-melt years and until runoff declines each melt season.

Figure 8Linear relationship between the satellite-derived water metric meltwater area fraction (MF, %) and RCM-derived runoff from MAR v3.11 for each year studied up until runoff declines.


5 Discussion

5.1 Spatial characteristics of the supraglacial drainage network

Our satellite-derived mapping reveals a complex supraglacial drainage network at HG that reaches a maximum elevation of 1440 m a.s.l., recorded in 2020. This maximum elevation is 400 m lower than supraglacial drainage observations in southwestern Greenland ( 1800 m a.s.l.; Yang et al., 2021), consistent with observations by Ryan et al. (2019), who found the maximum snowline elevation to be 400 m lower in northern Greenland than southwestern Greenland. Up to elevations of 1000 m a.s.l., the supraglacial drainage network is highly persistent across the study period, with a variable system beyond this elevation observed within the transient zone; defined in this study as a high-elevation region where drainage is transient in time but not necessarily space. The recurrence frequency of meltwater pixels in the supraglacial drainage network (Fig. 4) demonstrates the stability of SGLs in particular, which re-occupy the same location year-on-year despite ice advection (Smith et al., 2015; Pitcher and Smith, 2019). SGLs increase in size with elevation, with higher-elevation lakes persisting for longer (i.e. incapable of draining) and undergoing greater season expansion (Gledhill and Williamson, 2018; Yang et al., 2021). Larger lakes at higher elevations are the result of longer wavelength bed undulations being favourably transferred to the ice surface through thicker ice, creating large, shallow surface depressions (Ng et al., 2018). Unlike SGLs found < 400 m a.s.l., which are smaller and have likely reached their maximum available melt area, higher-elevation lakes are not yet topographically constrained (Krawczynski et al., 2009). This is similar to observations in southwestern Greenland where SGL size was larger and more variable > 1400 m a.s.l., with 21 % of these lakes draining via hydrofracture (Yang et al., 2021).

Whilst SGL location is known to be largely controlled by bed topography (Lampkin and VanderBerg, 2011; Ignéczi et al., 2018), this study also notes that many well-established rivers that are longitudinal to ice flow, including many with canyonised features, also reoccupy locations. Supraglacial rivers that are transverse to ice flow or have a transverse element to their channel, however, may be less stable in some areas by up  300 m (Fig. 4a) over the study period, probably due to the impact of ice advection. In the transient zone, the recurrence frequency is reduced, with tributary rivers and slush zones dominating at higher elevations. Here, their persistence is highly reliant on there being enough melt at higher elevations to initiate and sustain channel formation, which in 2017 and 2018 was limited, and therefore drainage occurrence was much reduced in this zone. Also, tributary rivers are typically lower order, with narrower channels and shallower depths (Smith et al., 2015; Pitcher and Smith, 2019), meaning their form has the potential to migrate, close and reform quickly if melt is available. The transient zone is therefore not only influenced by melt availability over time but also the potential for migration in space. In agreement with previous studies (Joughin et al., 2013; Poinar et al., 2015), we also show that rivers tend to be longer at higher elevations (> 25 km long), consistent with observations in SW GrIS (> 40 km long), likely due to the basal transfer of only long-wavelength basal undulations to the surface due to thicker ice and the reduced presence of surface crevassing (Gudmundsson, 2003; Lampkin and VanderBerg, 2011; Crozier et al., 2018).

Drainage patterns are also shown to vary across HG, with a dendritic style of drainage observed in the northern sector and a parallel style of drainage observed in the southern sector. These differing patterns highlight not only variations across different hydrologic catchments of the GrIS but also intra-catchment variations, which may stem from local variations in surface topography via the transmission of basal topography (Raymond and Gudmundsson, 2005; Ng et al., 2018) controlled by bed roughness and structure, wavelength transfer, and differing ice flow regimes (Gudmundsson et al., 1998; Gudmundsson, 2003; Lampkin and VanderBerg, 2011; Crozier et al., 2018; Ignéczi et al., 2018). This has been shown to play an important role in mapped rivers and lake hydromorphology at both the 79 N Glacier (Lu et al., 2021) and across the Devon Ice Cap (Wyatt and Sharp, 2015).

In the northern sector of HG, we observed short supraglacial rivers and small SGLs at lower elevations (200–400 m a.s.l.) and a prominent dendritic-style drainage pattern of interconnected rivers and lakes up to 1000 m a.s.l., with some larger rivers abruptly terminating. Such characteristics suggest the interception of runoff by crevasses and moulins, with such capture directed to the englacial and subglacial system with the potential for pronounced impacts on localised flow rates (Catania et al., 2008; Schoof, 2010; Mejia et al., 2022). This drainage style is typical of that observed within the western and southwestern sectors of GrIS, whereby structurally controlled drainage of supraglacial rivers flowing between SGLs promote shorter river channels, while high rates of crevassing at lower elevations (< 1000 m a.s.l.) means virtually all surface meltwater is captured and diverted before reaching the ice edge. This more compact style of drainage is likely to create a “flashier” response in hydrographs, with a greater runoff peak and shorter rising limb (Smith et al., 2017). Additionally, as shown by Carr et al. (2015) and Rignot et al. (2021), this sector of HG sits within a 475 m deep basal trough that extends  45 km wide and > 70 km inland and is characterised by fast rates of flow (200–600 m yr−1). Faster basal sliding has the ability to promote the more efficient transfer of basal topography to the surface and can subsequently precondition the large-scale spatial structure of the surface drainage system, including drainage fragmentation due to higher rates of crevassing and moulin formation, which can further supply meltwater to this faster-flowing sector (Crozier et al., 2018; Ignéczi et al., 2018).

In comparison, the sub-parallel drainage structure of the supraglacial network in the southern sector differs greatly to that of the northern sector of HG and observed drainage in western and southwestern regions of the GrIS. Drainage largely consists of continuously flowing rivers that drain surface meltwater from the slush zone at  1500 m a.s.l. to much lower elevations (200 m a.s.l.), with some rivers directly terminating off the ice sheet edge. This suggests limited opportunities for meltwater to penetrate to the ice sheet bed, with meltwater having longer transport times to travel to the catchment outlet (i.e. proglacial zone/Kane Basin) with the hydrograph expected to have a more subdued and longer rising limb (Karamouz et al., 2013; Yang et al., 2019a). Within this southern sector of HG, ice velocity is significantly slower (< 100 m yr−1; Rignot et al., 2021) than its northern counterpart, with relatively thick ice contributing to the absence of crevasses and moulins (Oswald and Gogineni, 2011; Yang et al., 2019a; Andrews et al., 2022), as well as controlling the hydromorphology of the drainage network found here. Similar drainage hydromorphology was also mapped at the neighbouring glacier at Inglefield Land (Yang et al., 2019a; Li et al., 2022), with supraglacial rivers flowing uninterrupted into the proglacial zone, with discharge into terrestrial rivers directly reflecting the timing and intensity of surface meltwater runoff from the outlet glacier catchment without modification from englacial and/or subglacial processes.

The supraglacial drainage configuration is also further influenced by significant structural elements that were identified by Livingstone et al. (2017) via the Moderate-Resolution Imaging Spectroradiometer (MODIS) mosaic of Greenland (Haran et al., 2013). Two linear structures expressed as depressions on the ice surface that run in a southwest to northeast direction across HG are clearly visible on the mapped glacier surface in this study (Figs. 4a, 9), with many supraglacial rivers and SGLs aligned-to or terminating at them. Some longitudinal rivers are also shown to suddenly change direction when encountering these structures, with subsequent channels diverting at a 90 angle, transverse to ice flow (Fig. 4a, b). It is at the intersection of such structures we observe some channel advection with spacing of 300 m (Fig. 4a), broadly representing the ice displacement over the study period. Other basal structures are also reflected within the supraglacial drainage system, including many V- and X-shaped patterns clearly controlled by depressions in the bed. There is a strong glacier-wide emphasis here on the influential control of these structures on drainage, which are shown here in Fig. 7 based on Livingstone et al. (2017) from MODIS imagery (MOG2015, Haran et al., 2018) and also within bed topography data from BedMachine (version 4; Morlighem et al., 2017, 2021). This identification provides independent confirmation of the existence of these depressions in the bed and their subsequent expression on the surface, as well as how they significantly control the multi-year surface drainage configuration within the vicinity of such structures.

Figure 9(a) Combined map of the supraglacial drainage network at its maximum extent across 2016 to 2020 showing the two parallel structures orientated southwest to northeast identified by dashed red lines. Other V- and X-shaped structures in the drainage network are also highlighted by dashed purple lines. (b) Moderate-Resolution Imaging Spectroradiometer (MODIS) Mosaic of Greenland (MOG2015, Haran et al., 2018) showing evidence of the these structures on the ice surface. (c) Bed topography via BedMachine (Morlighem et al., 2021) showing the structures as depressions within the bed. (d) The supraglacial drainage network as presented in (a) overlain on top of bed topography presented in (c) showing the overall influence bed topography has on the surface drainage structure at HG.

The supraglacial drainage network at HG is also shown to be largely dominated by supraglacial rivers (98 %), consistent with findings by Lu et al. (2021) in northeastern Greenland where supraglacial rivers also dominate the surface drainage network (83.8 %). Such river contributions in northern Greenland are significantly higher than southwestern GrIS, whereby supraglacial rivers account for an average of 62 % of the network (Yang et al., 2021). The difference in contribution of rivers and SGLs to the drainage network between northern and southwestern Greenland is likely the result of (i) reduced lake distributions and storage > 800 m a.s.l. in northern Greenland due to surface topographic constraints (i.e. reduced number of available depressions due to ice thickness; Ignéczi et al., 2016; Lu et al., 2021) and (ii) steeper northern surface slopes promoting meltwater runoff over storage, with gentler elevation gradients in southwestern Greenland promoting a wider, more susceptible ablation zone to melt and meltwater ponding (Mikkelsen et al., 2016; van As et al., 2017). At HG, the maximum supraglacial RF recorded for a single date across the melt-prone zone reached 24.5 %, the largest RF value recorded to date from a supraglacial drainage mapping study (e.g. 5.3 % in northeastern Greenland, 5.5 % at the Devon Ice Cap and  14 % for southwestern Greenland). A reason for this large value is likely due to the inclusion of expansive hydrologically connected slush zones within hydrologic mapping, which are prevalent during the early-to-middle melt season when snow becomes water saturated. This slush zone progressively moves upglacier and is removed to reveal the bare ice below. Additionally, a lower river threshold (t5) was used in this study compared to other studies (e.g. t20, Lu et al., 2021), allowing narrower supraglacial rivers and slush zones to be captured, increasing the overall supraglacial drainage and river components. Whilst supraglacial rivers are likely to dominate the drainage network here, similar to northeastern Greenland (Lu et al., 2021), regardless of the threshold used due to reduced lake presence, care should be taken when directly comparing drainage network statistics between mapping studies from differing regions of the GrIS.

5.2 Temporal evolution of the supraglacial drainage network

Across the study period (2016–2020), we observe the seasonal development and inland evolution of the supraglacial drainage network as the melt season progresses and runoff increases up to the maximum melt extent, which typically occurs at the end of July. During the early melt season and at higher elevations as melt progresses, we observe the growth and inland migration of a large, poorly channelised slush zone (Greuell and Knap, 2000). Within this zone, slush flows can occur as surface melt percolates and saturates the snowpack, promoting slush mobilisation into topographical lows and initiating the reopening of perennially occupied channels (Cuffey and Paterson, 2010; Irvine-Fynn et al., 2011). On higher-draining slopes, such mobilisation can form slush-filled rills, which coalesce into networks of arborescent tributary channels, efficiently transporting melt to larger, primary river channels and subsequently the transportation of melt downglacier (Marston, 1983; Cuffey and Paterson, 2010; Chu, 2014; Rippin and Rawlins, 2021). The seasonal development of the drainage network is shown to transform from a system initially dominated by water percolation to one dominated by channelised, efficient flow (Fig. 5), further confirming behaviour identified across multiple supraglacial drainage mapping studies across the GrIS (Lu et al., 2021; Yang et al., 2021; Li et al., 2022).

The rate and extent of the spatial and temporal evolution of the supraglacial drainage network is highly variable between years. Several years within the last decade have been characterised by high air temperatures and extreme melt events, including 2 years represented within this study: 2016 and 2019. Both years, in particular 2019, experienced a strong negative North Atlantic Oscillation and simultaneously a positive East Atlantic Index and Greenland Blocking Index phase, which are associated with persistent, anticyclonic conditions over Greenland driving enhanced surface mass loss (Lim et al., 2016; Cullather et al., 2020; Zhang et al., 2022). Mass loss during summer 2019 in particular was promoted by enhanced solar radiation, reduced cloud cover and the northwestward advection of warm, moist air from the western margins as a result of such atmospheric variability (Hanna et al., 2021; Cullather et al., 2020; Tedesco and Fettweis, 2020; Elmes et al., 2021). Combined with low snow accumulation in the 2018/19 winter (Sasgen et al., 2020), extensive melting occurred along much of the Greenland coast, with surface melt experienced in the north being the highest on the record since 1948 (Tedesco and Fettweis, 2020). It was during this exceptional and long melt year that we observed the highest MF values (19.9 % or 2685 km2 recorded on the 30 June 2019) and second highest areal extent of supraglacial drainage network (1375 m a.s.l.; Fig. 5). Ablation continued throughout September (Sasgen et al., 2020; Tedesco et al., 2019; Tedesco and Fettweis, 2020); however, this was beyond our mapped time frame. The year 2019 also recorded the largest number of SGLs (527) and cumulative lake areas (151.8 km2). The hydrologic expansion of the drainage network was also rapid, in-line with a record early melt season event that, combined with low snow accumulation (Tedesco and Fettweis, 2020), promoted rapid snowpack warming, disintegration and exposure of the bare ice zone, resulting in an enhanced melt–albedo feedback mechanism (Tedesco and Fettweis, 2020). Similar findings by Turton et al. (2021) and Hochreuther et al. (2021) at 79 N Glacier recorded the largest SGL numbers and extents of their study periods in 2019, indicating the widespread impact of this extreme melt event, particularly across the GrIS northern sector.

Similar seasonal and multi-annual behavioural patterns during 2017 and 2018 were again also observed by Turton et al. (2021) and Lu et al. (2021) on the northeastern glacier of 79 N and Otto et al. (2022) on the northern Ryder Glacier. These findings all record a slow rate of SGL increase and late area peak (early-August) in 2017 and 2018, with Lu et al. (2021) confirming the late area peak in combined supraglacial drainage mapping of both lakes and rivers during August 2017. SGL mapping by Turton et al. (2021) also identified the delayed development and lower SGL presence during the 2018 melt season, with SGLs largely limited to < 900 m a.s.l. Such findings are consistent with observations in this study, with the onset and inland evolution of the supraglacial network, including both rivers and lakes, delayed by  1 month compared to other the other study years (2016, 2019, 2020) and the limited areal development of the network (< 1150 m a.s.l.). Both the melt seasons of 2017 and 2018 recorded below-average melt (1981–2010 reference period) and melt extents (32.9 % and 44 %, respectively) across the GrIS (Tedesco et al., 2017, 2018; Sasgen et al., 2020). There was heavy springtime snowfall and late surviving snow in bare-ice areas (Tedesco et al., 2017, 2018), consistent with a strongly positive average summer (JJA) North Atlantic Oscillation and a negative Greenland Blocking Index, hypothesised to inhibit surface melt and promote increased summertime snowfall (Ruan et al., 2019; Sasgen et al., 2020), leading to these anomalously cold summers.

The anomalously early spike in satellite-derived meltwater area recorded in early June 2020 raises questions as to how extreme melt years, such as 2019, may precondition the ice surface (Cullather et al., 2020) and affect surface conditions and subsequent surface hydrologic behaviour the following year (Culberg et al., 2021). Some SGLs on the GrIS have been found to persist throughout the winter months due to insulation from a layer of ice and/or snow (Koenig et al., 2015; Law et al., 2020; Schröder et al., 2020). This lake persistence includes the winter of 2019/2020, when late summer surface melt and high autumnal temperatures (August–November) are believed to have increased subsurface firn temperatures, delaying and even decreasing the ability for subsurface meltwater freezing in northern Greenland, contributing to higher totals of liquid-buried SGLs (Dunmire et al., 2021). This, alongside a drier than average (1981–2020) winter and spring (Tedesco and Fettweis, 2020; Moon et al., 2020), has the potential to cause the rapid disintegration of limited snow present and the subsequent swift exposure of the bare ice zone the following summer. This swift exposure would also include that of perennial rivers and lakes much earlier in the melt season than expected, leading to/ the increased MF value observed in this study despite low MAR-derived melt production and runoff values.

5.3 Satellite-derived MF and runoff simulations

We found a positive linear relationship between satellite-derived MF and MAR simulated R before runoff declined (Fig. 8), showing how the MF–R relationship can be used to reliably simulate seasonal surface meltwater variation and provide further understanding into how runoff is routed and stored, at least up to peak melt events. This finding supports other studies that have used satellite-derived meltwater metrics and RCM-modelled runoff which have focussed on the southwestern GrIS (Yang et al., 2021), Northern GrIS (Lu et al., 2021; Li et al., 2022) and the Devon Ice Cap (Lu et al., 2020). In terms of post-peak melt events, surface meltwater can still cover a substantially large area even after surface runoff has reduced and ceased, as seen in this study. This is known as the “delay” effect (Lu et al., 2021), whereby meltwater may continue to be routed or stored on the ice surface via the slow routing of melt out of snowpack or firn at higher elevations, the stagnation and subsequent preservation of transported meltwater in large supraglacial rivers, or the storage of melt in SGLs. Therefore, while this MF–R relationship is promising in providing comparative assessments between satellite observations and RCM-modelled runoff at HG, calculated runoff volume via satellite-based (Yang et al., 2021) or field-based measurements (Smith et al., 2017) are required to provide further validation of such relationships, particularly over space and time across a full melt season.

5.4 Future implications

Substantial changes have taken place at HG over the last 2 decades driven by atmospheric and oceanic change (Carr et al., 2015; Rignot et al., 2021). It is therefore important to consider HG and the overall northern region sensitivity to such warming under present climate scenarios. Northern Greenland is expected to undergo the greatest warming of the 21st century across the GrIS (Hill et al., 2017), and given its already low rates of winter accumulation compared to other ice sheet sectors (Goelzer et al., 2013), this region is likely to become ever more sensitive to climatic change in the future. Mapping performed within this study illustrates the multi-annual persistence of the supraglacial drainage network within this high-latitude region (Fig. 4a) and the rapid and extensive response of this system to high-melt years (Fig. 5). This response, particularly in the extreme melt year of 2019, can precondition the ice surface for the following melt season, resulting in earlier but widespread hydrologic activity and longer-lasting melt season. This preconditioning and subsequent behaviour is likely to become increasingly normalised as melt events and atmospheric variability, such as persistent blocking events, increase in frequency and intensity (Rahmstorf and Coumou, 2011; McLeod and Mote, 2016).

Inland migration of the supraglacial drainage network is also projected with continued warming (Leeson et al., 2015), with recent work already showing ablation area expansion and amplification of melt and runoff post-1990 across northern Greenland (Noël et al., 2019). Many more surface depressions for future SGL locations are present above the current northern equlibrium line altitude (Ignéczi et al., 2016), with the potential to accumulate high volumes of meltwater and feed lengthening overflow supraglacial rivers that extend tens of kilometres downstream to non-local, low-elevation moulins (Poinar et al., 2015). For ponded water, if ice becomes thin enough (Poinar et al., 2015) or localised ice columns become vulnerable to fracture from refrozen ice complexities within firn (ice blobs; Culberg et al., 2022), new hydrofracture events will bring such meltwater to isolated areas of the ice sheet bed. This could have knock-on impacts on ice flow, with the likely delivery of water and heat to a persistent inefficient subglacial system, where thicker, flatter ice may prohibit the development of an efficient subglacial system (Dow et al., 2014), enhancing ice flow (Christoffersen et al., 2018).

Persistent low-permeability ice slabs that block vertical percolation have continued to thicken over time in the lower-accumulation zone. Ice slabs are expected to enhance runoff from Greenland's interior, particularly in consecutive warm summers (MacFerrin et al., 2019). Enhanced runoff and inland expansion of the supraglacial drainage network will impact meltwater feedback processes, not only in driving overall SMB decline (Noël et al., 2021) but further impacting dynamical behaviour, including hydrofracture potential. At HG, a particular concern is the vulnerability of its northern terminus to increased hydrofracture events from greater melt runoff (Carr et al., 2015). Such events at HG have the potential to promote future rapid runaway retreat of HG, especially if the northern sector of the terminus retreats beyond its pinning point into the deepened bed (below sea level) in which is sits (Carr et al., 2015; Hillebrand et al., 2022).

As the northern region mainly consists of fast-flowing, marine-terminating outlet glaciers like HG that drain a large proportion of the GrIS, further understanding of the mechanisms that drive their dynamical behaviour, in particular those related to enhanced runoff, are required for predicting their future contribution to GrIS mass loss and subsequent sea level rise.

6 Conclusion

In this study, we mapped and quantified for the first time the spatial and temporal evolution of the supraglacial drainage network, including both rivers and lakes, using 10 m Sentinel-2 images from the melt seasons of 2016 to 2020 at Humboldt Glacier, northern Greenland. We identify an extensive supraglacial drainage network at Humboldt Glacier that is particularly prevalent up to 1000 m a.s.l. with a further variable transient zone extending up to  1440 m a.s.l. The seasonal evolutionary behaviour of this network migrates upglacier in response to increasing runoff as air temperatures rise throughout the melt season, with the network transforming from an inefficient system dominated by water percolation and slush flows to one dominated by channelised, efficient flow. Interannual variability in the extent and behaviour of the system is associated with high- and low-melt years across the study period, with the low-melt years of 2017 and 2018 having both limited and delayed spatial development. The extreme melt year of 2019 showed the extensive development and persistence of the supraglacial drainage network into September, which followed by low snow accumulation during the subsequent winter and spring, preconditioned the ice sheet for earlier hydrologic activity in 2020. This behaviour may become more representative with extreme melt events and longer-lasting melt seasons into the future. This work ultimately contributes to advancing our understanding of supraglacial hydrologic processes across the Greenland Ice Sheet by expanding detailed drainage mapping to other understudied regions of the ice sheet, in particular to Greenland's rapidly changing northern region, aiding in projections of future mass loss as enhanced runoff continues with climatic warming.

Data availability

Supraglacial river and lake shapefiles can be requested by contacting the lead author. Regional climate model MAR v3.11 datasets (Amory et al., 2021; Fettweis et al., 2021) were provided by Xavier Fettweis and are freely available via (, Amory et al., 2021, Fettweis et al., 2021). Sentinel-2 imagery is available from the Copernicus Open Access Hub (, Copernicus, 2022). The automatic river detection algorithm (Yang et al., 2015, 2017) to generate supraglacial network datasets is freely available on the GitHub site (, njuRS, 2022).


The supplement related to this article is available online at:

Author contributions

LDR designed the study, conducted data collection and analysis, and prepared the manuscript. KY provided the source code for the automatic detection algorithm for automatic river mapping in MATLAB. AJS aided with MAR data extraction for the Humboldt Glacier catchment. DMR, AJS, SJL, and KY provided comments on draft versions of the manuscript produced by LDR.

Competing interests

At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


Lauren D. Rawlins acknowledges financial support from a Natural Environmental Research Council (NERC) Doctoral Training Partnership (grant no. NE/L002450/1). The authors would like to acknowledge the Copernicus Open Access Hub for providing free access to the satellite imagery used in this study and also thank Xavier Fettweis for the availability of MAR data. We also thank the Polar Geospatial Center for the availability of the ArcticDEM. The authors would like to thank the two anonymous reviewers for providing constructive comments that led to improvements to the manuscript and to the editor (Caroline Clason) for handling the manuscript.

Financial support

Lauren D. Rawlins was supported by a PhD studentship via the Natural Environment Research Council (NERC) Doctoral Training Partnership (grant no. NE/L002450/1).

Review statement

This paper was edited by Caroline Clason and reviewed by two anonymous referees.


Amory, C., Kittel, C., Le Toumelin, L., Agosta, C., Delhasse, A., Favier, V., and Fettweis, X.: Performance of MAR (v3.11) in simulating the drifting-snow climate and surface mass balance of Adélie Land, East Antarctica, Geosci. Model Dev., 14, 3487–3510,, 2021 (data available at: 

Andrews, L. C., Catania, G. A., Hoffman, M. J., Gulley, J. D., Lüthi, M. P., Ryser, C., Hawley, R. L., and Neumann, T. A.: Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet, Nature, 514, 80–83,, 2014. 

Andrews, L. C., Poinar, K., and Trunz, C.: Controls on Greenland moulin geometry and evolution from the Moulin Shape model, The Cryosphere, 16, 2421–2448,, 2022. 

Baillarin, S., Meygret, A., Dechoz, C., Petrucci, B., Lacherade, S.,Trémas, T., Isola, C., Martimort, P., and Spoto, F.: Sentinel-2 level 1 products and image processing performances, IEEE international geoscience and remote sensing symposium, Munich, Germany, 22–27 July 2012, 7003–7006,, 2012. 

Banwell, A. F., Cabellero, M., Arnold, N., Glasser, N., Cathles, L. M., and MacAyeal, D.: Supraglacial lakes on the Larsen B Ice Shelf, Antarctica, and Paakitsoq Region, Greenland: a comparative study, Ann. Glaciol., 55, 66,, 2014. 

Bartholomew, I., Nienow, P., Mair, D., Hubbard, A., King, M. A., and Sole, A.: Seasonal evolution of subglacial drainage and acceleration in a Greenland outlet glacier, Nat. Geosci., 3, 408–411,, 2010. 

Bartholomew, I., Nienow, P., Sole, A., Mair, D., Cowton, T., and King, M. A.: Short-term variability in Greenland Ice Sheet motion forced by time-varying meltwater drainage: Implications for the relationship between subglacial drainage system behavior and ice velocity, J. Geophys. Res.-Earth, 117, F03002,,, 2012. 

Boghosian, A. L., Pitcher, L. H., Smith, L. C., Kosh, E., Alexander, P. M., Tedesco, M., and Bell, R. E.: Development of ice-shelf estuaries promotes fractures and calving, Nat. Geosci., 14, 899–905,, 2021. 

Box, J. E. and Decker, D. T.: Greenland marine-terminating glacier area changes: 2000–2010, Ann. Glaciol., 52, 91–98,, 2011. 

Carr, J. R., Vieli, A., Stokes, C., Jamieson, S., Palmer, S., Christoffersen, P., Dowdeswell, J., Nick, F., Blankenship, D., and Young, D.: Basal topographic controls on rapid retreat of Humboldt Glacier, northern Greenland, J. Glaciol., 61, 137–150,, 2015. 

Catania, G. A., Neumann, T. A., and Price, S. F.: Characterizing englacial drainage in the ablation zone of the Greenland ice sheet, J. Glaciol., 54, 567–578,, 2008. 

Christoffersen, P., Bougamont, M., Hubbard, A., Doyle, S. H., Grigsby, S., and Pettersson, R.: Cascading lake drainage on the Greenland Ice Sheet triggered by tensile shock and fracture, Nat. Commun., 9, 1064,, 2018. 

Chu, V. W.: Greenland ice sheet hydrology: A review, Prog. Phys. Geogr., 38, 19–54,, 2014. 

Copernicus: Copernicus Open Access Hub,, last access: 15 December 2022. 

Corr, D., Leeson, A., McMillan, M., Zhang, C., and Barnes, T.: An inventory of supraglacial lakes and channels across the West Antarctic Ice Sheet, Earth Syst. Sci. Data, 14, 209–228,, 2022. 

Crozier, J., Karlstrom, L., and Yang, K.: Basal control of supraglacial meltwater catchments on the Greenland Ice Sheet, The Cryosphere, 12, 3383–3407,, 2018. 

Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, Academic Press, ISBN 978-0-12-369461-4, 2010. 

Culberg, R., Schroeder, D. M., and Chu, W.: Extreme melt season ice layers reduce firn permeability across Greenland, Nat. Commun., 12, 2336,, 2021. 

Culberg, R., Chu, W., and Schroeder, D. M.: Shallow Fracture Buffers High Elevation Runoff in Northwest Greenland, Geophys. Res. Lett., 49, e2022GL101151,, 2022. 

Cullather, R. I., Andrews, L. C., Croteau, M. J., Digirolamo, N. E., Hall, D. K., Lim, Y. K., Loomis, B. D., Shuman, C. A., and Nowicki, S. M.: Anomalous circulation in July 2019 resulting in mass loss on the Greenland Ice Sheet, Geophys. Res. Lett., 47, e2020GL087263,, 2020. 

Davison, B. J., Sole, A. J., Livingstone, S. J., Cowton, T. R., and Nienow, P. W.: The influence of hydrology on the dynamics of land-terminating sectors of the Greenland ice sheet, Front. Earth Sci., 7, 10,, 2019. 

Dell, R., Banwell, A. F., Willis, I., Arnold, N., Halberstadt, A. R. W., Chudley, T. R., and Pritchard, H.: Supervised classification of slush and ponded water on Antarctic ice shelves using Landsat 8 imagery, J. Glaciol., 68, 401–414,, 2022. 

Dow, C., Kulessa, B., Rutt, I., Doyle, S. H., and Hubbard, A.: Upper bounds on subglacial channel development for interior regions of the Greenland ice sheet, J. Glaciol., 60, 1044–1052,, 2014. 

Dunmire, D., Banwell, A. F., Wever, N., Lenaerts, J. T. M., and Datta, R. T.: Contrasting regional variability of buried meltwater extent over 2 years across the Greenland Ice Sheet, The Cryosphere, 15, 2983–3005,, 2021. 

Elmes, A., Levy, C., Erb, A., Hall, D. K., Scambos, T. A., DiGirolamo, N., and Schaaf, C.: Consequences of the 2019 greenland ice sheet melt episode on albedo, Remote Sensing, 13, 227,, 2021. 

Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033,, 2017. 

Fettweis, X., Hofer, S., Krebs-Kanzow, U., Amory, C., Aoki, T., Berends, C. J., Born, A., Box, J. E., Delhasse, A., Fujita, K., Gierz, P., Goelzer, H., Hanna, E., Hashimoto, A., Huybrechts, P., Kapsch, M.-L., King, M. D., Kittel, C., Lang, C., Langen, P. L., Lenaerts, J. T. M., Liston, G. E., Lohmann, G., Mernild, S. H., Mikolajewicz, U., Modali, K., Mottram, R. H., Niwano, M., Noël, B., Ryan, J. C., Smith, A., Streffing, J., Tedesco, M., van de Berg, W. J., van den Broeke, M., van de Wal, R. S. W., van Kampenhout, L., Wilton, D., Wouters, B., Ziemen, F., and Zolles, T.: GrSMBMIP: intercomparison of the modelled 1980–2012 surface mass balance over the Greenland Ice Sheet, The Cryosphere, 14, 3935–3958,, 2020. 

Fettweis, X., Hofer, S., Séférian, R., Amory, C., Delhasse, A., Doutreloup, S., Kittel, C., Lang, C., Van Bever, J., Veillon, F., and Irvine, P.: Brief communication: Reduction in the future Greenland ice sheet surface melt with the help of solar geoengineering , The Cryosphere, 15, 3013–3019,, 2021 (data available at: 

Gardner, A. S., Fahnestock, M. A., and Scambos, T. A.: ITS_LIVE Regional Glacier and Ice Sheet Surface Velocities, National Snow and Ice Data Center [data set],, 2019. 

Gleason, C. J., Smith, L. C., Chu, V. W., Legleiter, C. J., Pitcher, L. H., Overstreet, B. T., Rennermalm, A. K., Forster, R. R., and Yang, K.: Characterizing supraglacial meltwater channel hydraulics on the Greenland Ice Sheet from in situ observations, Earth Surf. Proc. Land., 41, 2111–2122,, 2016. 

Gleason, C. J., Yang, K., Feng, D., Smith, L. C., Liu, K., Pitcher, L. H., Chu, V. W., Cooper, M. G., Overstreet, B. T., and Rennermalm, A. K.: Hourly surface meltwater routing for a Greenlandic supraglacial catchment across hillslopes and through a dense topological channel network, The Cryosphere, 15, 2315–2331,, 2021. 

Gledhill, L. A. and Williamson, A. G.: Inland advance of supraglacial lakes in north-west Greenland under recent climatic warming, Ann. Glaciol., 59, 66–82,, 2018. 

Goelzer, H., Huybrechts, P., Fürst, J. J., Nick, F. M., Andersen, M. L., Edwards, T. L., Fettweis, X., Payne, A. J., and Shannon, S.: Sensitivity of Greenland ice sheet projections to model formulations, J. Glaciol., 59, 733–749,, 2013. 

Gray, L.: Brief communication: Glacier run-off estimation using altimetry-derived basin volume change: case study at Humboldt Glacier, northwest Greenland, The Cryosphere, 15, 1005–1014,, 2021. 

Greuell, W. and Knap, W. H.: Remote sensing of the albedo and detection of the slush line on the Greenland ice sheet, J. Geophys. Res.-Atmos., 105, 15567–15576,, 2000. 

Gudmundsson, G. H.: Transmission of basal variability to a glacier surface, J. Geophys. Res.-Sol. Ea., 108, 2253,, 2003. 

Gudmundsson, G. H., Raymond, C. F., and Bindschadler, R.: The origin and longevity of flow stripes on Antarctic ice streams, Ann. Glaciol., 27, 145–152,, 1998. 

Hanna, E., Mernild, S. H., Cappelen, J., and Steffen, K.: Recent warming in Greenland in a long-term instrumental (1881–2012) climatic context: I. Evaluation of surface air temperature records, Environ. Res. Lett., 7, 045404,, 2012. 

Hanna, E., Cappelen, J., Fettweis, X., Mernild, S. H., Mote, T. L., Mottram, R., Steffen, K., Ballinger, T. J., and Hall, R. J.: Greenland surface air temperature changes from 1981 to 2019 and implications for ice-sheet melt and mass-balance change, Int. J. Climatol., 41, E1336–E1352,, 2021. 

Haran, T., Bohlander, J., Scambos, T., Painter, T., and Fahnestock, M.: MEaSUREs MODIS Mosaic of Greenland 2005 (MOG2005) Image Map, Version 1, Boulder, Colorado, NSIDC: National Snow and Ice Data Center [data set] ,, 2013. 

Haran, T., Bohlander, J., Scambos, T., Painter, T., and Fahnestock, M.: MEaSUREs MODIS Mosaic of Greenland (MOG) 2005, 2010, and 2015 Image Maps, Version 2, Boulder, Colorado USA, NASA National Snow and Ice Data Center Distributed Active Archive Center [data set],, 2018. 

Hill, E. A., Carr, J. R., and Stokes, C. R.: A review of recent changes in major marine-terminating outlet glaciers in Northern Greenland, Front. Earth Sci., 4, 111,, 2017. 

Hill, E. A., Carr, J. R., Stokes, C. R., and Gudmundsson, G. H.: Dynamic changes in outlet glaciers in northern Greenland from 1948 to 2015, The Cryosphere, 12, 3243–3263,, 2018. 

Hillebrand, T. R., Hoffman, M. J., Perego, M., Price, S. F., and Howat, I. M.: The contribution of Humboldt Glacier, northern Greenland, to sea-level rise through 2100 constrained by recent observations of speedup and retreat, The Cryosphere, 16, 4679–4700,, 2022. 

Hochreuther, P., Neckel, N., Reimann, N., Humbert, A., and Braun, M.: Fully automated detection of supraglacial lake area for Northeast Greenland using sentinel-2 time-series, Remote Sensing, 13, 205,, 2021. 

Hoffman, M., Catania, G., Neumann, T., Andrews, L., and Rumrill, J.: Links between acceleration, melting, and supraglacial lake drainage of the western Greenland Ice Sheet, J. Geophys. Res.-Earth, 116, F04035,, 2011. 

Holmes, G. W.: Morphology and hydrology of the Mint Julep area, southwest Greenland, in: Project Mint Julep Investigation of Smooth Ice Areas of the Greenland Ice Cap, 1953, Part II Special Scientific Reports, Arctic, Desert, Tropic Information Center, Research Studies Institute, Air University, 1955 

Ignéczi, Á., Sole, A. J., Livingstone, S. J., Leeson, A. A., Fettweis, X., Selmes, N., Gourmelen, N., and Briggs, K.: Northeast sector of the Greenland Ice Sheet to undergo the greatest inland expansion of supraglacial lakes during the 21st century, Geophys. Res. Lett., 43, 9729–9738,, 2016. 

Ignéczi, Á., Sole, A. J., Livingstone, S. J., Ng, F. S., and Yang, K.: Greenland Ice Sheet surface topography and drainage structure controlled by the transfer of basal variability, Front. Earth Sci., 6, 101,, 2018. 

Irvine-Fynn, T. D., Hodson, A. J., Moorman, B. J., Vatne, G., and Hubbard, A. L.: Polythermal glacier hydrology: A review, Rev. Geophys., 49, RG4002,, 2011. 

Joughin, I., Kwok, R., and Fahnestock, M.: Estimation of ice-sheet motion using satellite radar interferometry: method and error analysis with application to Humboldt Glacier, Greenland, J. Glaciol., 42, 564–575,, 1996. 

Joughin, I., Fahnestock, M., Kwok, R., Gogineni, P., and Allen, C.: Ice flow in the Humboldt, Petermann, and Ryder Glaciers, North Greenland, J. Glaciol., 45, 231–341,, 1999. 

Joughin, I., Das, S. B., Flowers, G. E., Behn, M. D., Alley, R. B., King, M. A., Smith, B. E., Bamber, J. L., van den Broeke, M. R., and van Angelen, J. H.: Influence of ice-sheet geometry and supraglacial lakes on seasonal ice-flow variability, The Cryosphere, 7, 1185–1192,, 2013. 

Kamb, B.: Glacier surge mechanism based on linked cavity configuration of the basal water conduit system, J. Geophys. Res.-Sol. Ea., 92, 9083–9100,, 1987. 

Karamouz, M., Nazif, S., and Falahi, M.:. Hydrology and hydroclimatology: principles and applications, CRC Press, ISBN 978-1-4665-1220-7,, 2013. 

Karlstrom, L. and Yang, K.: Fluvial supraglacial landscape evolution on the Greenland Ice Sheet, Geophys. Res. Lett., 43, 2683–2692,, 2016. 

Koenig, L. S., Lampkin, D. J., Montgomery, L. N., Hamilton, S. L., Turrin, J. B., Joseph, C. A., Moutsafa, S. E., Panzer, B., Casey, K. A., Paden, J. D., Leuschen, C., and Gogineni, P.: Wintertime storage of water in buried supraglacial lakes across the Greenland Ice Sheet, The Cryosphere, 9, 1333–1342,, 2015. 

Krawczynski, M. J., Behn, M. D., Das, S. B., and Joughin, I.: Constraints on the lake volume required for hydro-fracture through ice sheets, Geophys. Res. Lett., 36, L10501,, 2009. 

Lampkin, D. and VanderBerg, J.: A preliminary investigation of the influence of basal and surface topography on supraglacial lake distribution near Jakobshavn Isbrae, western Greenland, Hydrol. Process., 25, 3347–3355,, 2011. 

Lampkin, D. and VanderBerg, J.: Supraglacial melt channel networks in the Jakobshavn Isbræ region during the 2007 melt season, Hydrol. Process., 28, 6038–6053,, 2014. 

Law, R., Arnold, N., Benedek, C., Tedesco, M., Banwell, A., and Willis, I.: Over-winter persistence of supraglacial lakes on the Greenland Ice Sheet: results and insights from a new model, J. Glaciol., 66, 362–372,, 2020. 

Leeson, A., Shepherd, A., Briggs, K., Howat, I., Fettweis, X., Morlighem, M., and Rignot, E.: Supraglacial lakes on the Greenland ice sheet advance inland under warming climate, Nat. Clim. Change, 5, 51–55,, 2015. 

Li, Y., Yang, K., Gao, S., Smith, L. C., Fettweis, X., and Li, M.: Surface meltwater runoff routing through a coupled supraglacial-proglacial drainage system, Inglefield Land, northwest Greenland, Int. J. Appl. Earth Obs., 106, 102647,, 2022. 

Lim, Y.-K., Schubert, S. D., Nowicki, S. M., Lee, J. N., Molod, A. M., Cullather, R. I., Zhao, B., and Velicogna, I.: Atmospheric summer teleconnections and Greenland Ice Sheet surface mass variations: Insights from MERRA-2, Environ. Res. Lett., 11, 024002,, 2016. 

Livingstone, S. J., Chu, W., Ely, J. C., and Kingslake, J.: Paleofluvial and subglacial channel networks beneath Humboldt Glacier, Greenland, Geology, 45, 551–554,, 2017. 

Lu, Y., Yang, K., Lu, X., Smith, L. C., Sole, A. J., Livingstone, S. J., Fettweis, X., and Li, M.: Diverse supraglacial drainage patterns on the Devon ice Cap, Arctic Canada, J. Maps, 16, 834–846,, 2020. 

Lu, Y., Yang, K., Lu, X., Li, Y., Gao, S., Mao, W., and Li, M.: Response of supraglacial rivers and lakes to ice flow and surface melt on the northeast Greenland ice sheet during the 2017 melt season, J. Hydrol., 602, 126750,, 2021. 

Macdonald, G., Banwell, A., and MacAyeal, D.: Seasonal evolution of supraglacial lakes on a floating ice tongue, Petermann Glacier, Greenland, Ann. Glaciol., 59, 56–65, https://, 2018. 

MacFerrin, M., Machguth, H., As, D. v., Charalampidis, C., Stevens, C., Heilig, A., Vandecrux, B., Langen, P., Mottram, R., and Fettweis, X.: Rapid expansion of Greenland's low-permeability ice slabs, Nature, 573, 403–407,, 2019. 

Marston, R. A.: Supraglacial stream dynamics on the Juneau Icefield, Ann. Assoc. Am. Geogr., 73, 597-608,, 1983. 

McFeeters, S. K.: The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features, Int. J. Remote Sens., 17, 1425–1432,, 1996. 

McLeod, J. T. and Mote, T. L.: Linking interannual variability in extreme Greenland blocking episodes to the recent increase in summer melting across the Greenland ice sheet, Int. J. Climatol., 36, 1484–1499,, 2016. 

Mejia, J., Gulley, J., Trunz, C., Covington, M. D., Bartholomaus, T., Breithaupt, C., Xie, S., and Dixon, T. H.: Moulin density controls the timing of peak pressurization within the Greenland Ice Sheet's subglacial drainage system, Geophys. Res. Lett., 49, e2022GL100058,, 2022. 

Mikkelsen, A. B., Hubbard, A., MacFerrin, M., Box, J. E., Doyle, S. H., Fitzpatrick, A., Hasholt, B., Bailey, H. L., Lindbäck, K., and Pettersson, R.: Extraordinary runoff from the Greenland ice sheet in 2012 amplified by hypsometry and depleted firn retention, The Cryosphere, 10, 1147–1159,, 2016. 

Moon, T. A., Tedesco, M., Box, J., Cappelen, J., Fausto, R., Fettweis, X., Korsgaard, N., Loomis, B., Mankoff, K., and Mote, T.: Arctic Report Card 2020: Greenland Ice Sheet, National Snow and Ice Data Center (U.S.), University of Colorado Boulder, United States, National Oceanic and Atmospheric Administration. Office of Oceanic and Atmospheric Research,, 2020 

Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., and Dorschel, B.: BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multibeam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, 11051–11061,, 2017. 

Morlighem, M., et al.: IceBridge BedMachine Greenland, Version 4. [Bed Topography], Boulder, Colorado USA, NASA National Snow and Ice Data Center Distributed Active Archive Center [data set],, 2021 

Mouginot, J., Rignot, E., Bjørk, A. A., Van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.: Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, P. Natl. Acad. Sci. USA, 116, 9239–9244,, 2019. 

Moussavi, M. S., Abdalati, W., Pope, A., Scambos, T., Tedesco, M., MacFerrin, M., and Grigsby, S.: Derivation and validation of supraglacial lake volumes on the Greenland Ice Sheet from high-resolution satellite imagery, Remote Sens. Environ., 183, 294–303,, 2016. 

Ng, F. S., Ignéczi, Á., Sole, A. J., and Livingstone, S. J.: Response of surface topography to basal variability along glacial flowlines, J. Geophys. Res.-Earth, 123, 2319–2340,, 2018. 

Nienow, P., Sole, A., Slater, D. A., and Cowton, T.: Recent advances in our understanding of the role of meltwater in the Greenland Ice Sheet system, Current Climate Change Reports, 3, 330–344,, 2017. 

Noël, B., van de Berg, W. J., Lhermitte, S., and van den Broeke, M. R.: Rapid ablation zone expansion amplifies north Greenland mass loss, Science Advances, 5, eaaw0123,, 2019. 

Noël, B., van Kampenhout, L., Lenaerts, J. T. M., van de Berg, W. J., and van den Broeke, M. R.: A 21st century warming threshold for sustained Greenland ice sheet mass loss, Geophys. Res. Lett., 48, e2020GL090471,, 2021. 

njuRS: River_detection, GitHub [data set],, last access: 15 December 2022. 

Oswald, G. K. and Gogineni, S.: Mapping basal melt under the northern Greenland Ice Sheet, IEEE T. Geosci. Remote, 50, 585–592,, 2011. 

Otto, J., Holmes, F. A., and Kirchner, N.: Supraglacial lake expansion, intensified lake drainage frequency, and first observation of coupled lake drainage, during 1985–2020 at Ryder Glacier, Northern Greenland, Front. Earth Sci., 10, 978137,, 2022. 

Pitcher, L. H. and Smith, L. C.: Supraglacial streams and rivers, Annu. Rev. Earth Planet. Sci., 47, 421–452,, 2019. 

Poinar, K., Joughin, I., Das, S. B., Behn, M. D., Lenaerts, J. T., and Van Den Broeke, M. R.: Limits to future expansion of surface-melt-enhanced ice flow into the interior of western Greenland, Geophys. Res. Lett., 42, 1800–1807,, 2015. 

Pope, A., Scambos, T. A., Moussavi, M., Tedesco, M., Willis, M., Shean, D., and Grigsby, S.: Estimating supraglacial lake depth in West Greenland using Landsat 8 and comparison with other multispectral methods, The Cryosphere, 10, 15–27,, 2016. 

Rahmstorf, S. and Coumou, D.: Increase of extreme events in a warming world, P. Natl. Acad. Sci. USA, 108, 17905–17909,, 2011. 

Raymond, M. J. and Gudmundsson, G. H.: On the relationship between surface and basal properties on glaciers, ice sheets, and ice streams, J. Geophys. Res.-Sol. Ea., 110, B08411,, 2005. 

Rennermalm, A. K., Smith, L. C., Chu, V. W., Box, J. E., Forster, R. R., Van den Broeke, M. R., Van As, D., and Moustafa, S. E.: Evidence of meltwater retention within the Greenland ice sheet, The Cryosphere, 7, 1433–1445,, 2013. 

Rignot, E. and Kanagaratnam, P.: Changes in the velocity structure of the Greenland Ice Sheet, Science, 311, 986–990,, 2006. 

Rignot, E., Gogineni, S., Joughin, I., and Krabill, W.: Contribution to the glaciology of northern Greenland from satellite radar interferometry, J. Geophys. Res., 106, 34007–34019,, 2001. 

Rignot, E., An, L., Chauche, N., Morlighem, M., Jeong, S., Wood, M., Mouginot, J., Willis, J. K., Klaucke, I., and Weinrebe, W.: Retreat of Humboldt Gletscher, North Greenland, driven by undercutting from a warmer ocean, Geophys. Res. Lett., 48, e2020GL091342,, 2021. 

Riihelä, A., King, M. D., and Anttila, K.: The surface albedo of the Greenland Ice Sheet between 1982 and 2015 from the CLARA-A2 dataset and its relationship to the ice sheet's surface mass balance, The Cryosphere, 13, 2597–2614,, 2019. 

Rippin, D. and Rawlins, L.: Supraglacial River Networks, In International Encyclopedia of Geography, edited by: Richardson, D., Castree, N., Goodchild, M. F., Kobayashi, A., Liu, W., and Marston, R. A.,, 2021. 

Ruan, R., Chen, X., Zhao, J., Perrie, W., Mottram, R., Zhang, M., Diao, Y., Du, L., and Wu, L.: Decelerated Greenland Ice Sheet melt driven by positive summer North Atlantic oscillation, J. Geophys. Res.-Atmos., 124, 7633–7646,, 2019. 

Ryan, J., Smith, L., Van As, D., Cooley, S., Cooper, M., Pitcher, L., and Hubbard, A.: Greenland Ice Sheet surface melt amplified by snowline migration and bare ice exposure, Science Advances, 5, eaav3738,, 2019. 

Ryan, J. C., Hubbard, A., Stibal, M., Irvine-Fynn, T. D., Cook, J., Smith, L. C., Cameron, K., and Box, J.: Dark zone of the Greenland Ice Sheet controlled by distributed biologically-active impurities, Nat. Commun., 9, 1065,, 2018. 

Sasgen, I., Wouters, B., Gardner, A. S., King, M. D., Tedesco, M., Landerer, F. W., Dahle, C., Save, H., and Fettweis, X.: Return to rapid ice loss in Greenland and record loss in 2019 detected by the GRACE-FO satellites, Commun. Earth Environ., 1, 1–8,, 2020. 

Schoof, C.: Ice-sheet acceleration driven by melt supply variability, Nature, 468, 803–806,, 2010. 

Schröder, L., Neckel, N., Zindler, R., and Humbert, A.: Perennial supraglacial lakes in Northeast Greenland observed by polarimetric SAR, Remote Sensing, 12, 2798,, 2020. 

Selmes, N., Murray, T., and James, T.: Fast draining lakes on the Greenland Ice Sheet, Geophys. Res. Lett., 38, L15501,, 2011. 

Slater, T., Shepherd, A., McMillan, M., Leeson, A., Gilbert, L., Muir, A., Munneke, P. K., Noël, B., Fettweis, X., and van den Broeke, M.: Increased variability in Greenland Ice Sheet runoff from satellite observations, Nat. Commun., 12, 6069,, 2021. 

Smith, L. C., Chu, V. W., Yang, K., Gleason, C. J., Pitcher, L. H., Rennermalm, A. K., Legleiter, C. J., Behar, A. E., Overstreet, B. T., and Moustafa, S. E.: Efficient meltwater drainage through supraglacial streams and rivers on the southwest Greenland ice sheet, P. Natl. Acad. Sci. USA, 112, 1001–1006,, 2015. 

Smith, L. C., Yang, K., Pitcher, L. H., Overstreet, B. T., Chu, V. W., Rennermalm, Å. K., Ryan, J. C., Cooper, M. G., Gleason, C. J., and Tedesco, M.: Direct measurements of meltwater runoff on the Greenland ice sheet surface, P. Natl. Acad. Sci. USA, 114, E10622–E10631,, 2017. 

Sole, A. J., Mair, D. W. F., Nienow, P. W., Bartholomew, I., King, M., Burke, M. J., and Joughin, I.: Seasonal speedup of a Greenland marine-terminating outlet glacier forced by surface melt–induced changes in subglacial hydrology, J. Geophys. Res.-Earth, 116, F03014,, 2011. 

Stokes, C. R., Sanderson, J. E., Miles, B. W. J. Jamieson, S. S., and Leeson, A. A.: Widespread distribution of supraglacial lakes around the margin of the East Antarctic Ice Sheet, Sci. Rep., 9, 13823,, 2019. 

Tedesco, M. and Fettweis, X.: Unprecedented atmospheric conditions (1948–2019) drive the 2019 exceptional melting season over the Greenland ice sheet, The Cryosphere, 14, 1209–1223,, 2020. 

Tedesco, M., Doherty, S., Fettweis, X., Alexander, P., Jeyaratnam, J., and Stroeve, J.: The darkening of the Greenland ice sheet: trends, drivers, and projections (1981–2100), The Cryosphere, 10, 477–496,, 2016. 

Tedesco, M., Box, J. E., Capplelen, J., Fausto, R. S., Fettweis, X., Hansen, K., Mote, T., Sasgen, I., Smeets, C. J. P. P., van As, D., van de Wal, R. S. W., and Velicogna, I.: NOAA Arctic Report Card 2018: Greenland Ice Sheet in Arctic Report Card 2017, (last access: 20 November 2022), 2017. 

Tedesco, M., Box, J. E., Cappelen, J., Fausto, R. S., Fettweis, X., Anderson, J. K., Mote, T., Smeets, C. J. P. P., van As, D., and van de Wal, R. S. W.: NOAA Arctic Report Card 2018: Greenland Ice Sheet in Arctic Report Card 2018, (last access: 20 November 2022), 2018. 

Tedesco, M., Moon, T., Anderson, J. K., Box, J. E., Cappelen, J., Fausto, R. S., Fettweis, X., Loomis, B., Mankoff, K. D., Mote, T., Smeets, C. J. P. P., van As, D., and van de Wal, R. S. W.: Greenland Ice Sheet in Arctic Report Card 2019, (last access: 21 November 2022), 2019. 

The IMBIE Team: Mass balance of the Greenland Ice Sheet from 1992 to 2018, Nature, 579, 233–239,, 2020. 

Trusel, L. D., Das, S. B., Osman, M. B., Evans, M. J., Smith, B. E., Fettweis, X., McConnell, J. R., Noël, B. P., and van den Broeke, M. R.: Nonlinear rise in Greenland runoff in response to post-industrial Arctic warming, Nature, 564, 104–108,, 2018. 

Turton, J. V., Hochreuther, P., Reimann, N., and Blau, M. T.: The distribution and evolution of supraglacial lakes on 79° N Glacier (north-eastern Greenland) and interannual climatic controls, The Cryosphere, 15, 3877–3896,, 2021. 

van As, D., Bech Mikkelsen, A., Holtegaard Nielsen, M., Box, J. E., Claesson Liljedahl, L., Lindbäck, K., Pitcher, L., and Hasholt, B.: Hypsometric amplification and routing moderation of Greenland ice sheet meltwater release, The Cryosphere, 11, 1371–1386,, 2017. 

van den Broeke, M., Box, J., Fettweis, X., Hanna, E., Noël, B., Tedesco, M., van As, D., van de Berg, W. J., and van Kampenhout, L.: Greenland ice sheet surface mass loss: recent developments in observation and modeling, Current Climate Change Reports, 3, 345–356,, 2017. 

van den Broeke, M. R., Enderlin, E. M., Howat, I. M., Kuipers Munneke, P., Noël, B. P. Y., van de Berg, W. J., van Meijgaard, E., and Wouters, B.: On the recent contribution of the Greenland ice sheet to sea level change, The Cryosphere, 10, 1933–1946,, 2016. 

Williamson, A. G., Arnold, N. S., Banwell, A. F., and Willis, I. C.: A Fully Automated Supraglacial lake area and volume Tracking (“FAST”) algorithm: Development and application using MODIS imagery of West Greenland, Remote Sens. Environ., 196, 113–133,, 2017. 

Williamson, A. G., Banwell, A. F., Willis, I. C., and Arnold, N. S.: Dual-satellite (Sentinel-2 and Landsat 8) remote sensing of supraglacial lakes in Greenland, The Cryosphere, 12, 3045–3065,, 2018. 

Wyatt, F. R. and Sharp, M. J.: Linking surface hydrology to flow regimes and patterns of velocity variability on Devon Ice Cap, Nunavut, J. Glaciol., 61, 387–399,, 2015. 

Yang, K. and Smith, L. C.: Supraglacial streams on the Greenland Ice Sheet delineated from combined spectral–shape information in high-resolution satellite imagery, IEEE Geosci. Remote Se., 10, 801–805,, 2012. 

Yang, K., Smith, L. C., Chu, V. W., Gleason, C. J., and Li, M.: A caution on the use of surface digital elevation models to simulate supraglacial hydrology of the Greenland ice sheet, IEEE J. Sel. Top. Appl., 8, 5212–5224,, 2015. 

Yang, K., Smith, L. C., Sole, A., Livingstone, S. J., Cheng, X., Chen, Z., and Li, M.: Supraglacial rivers on the northwest Greenland Ice Sheet, Devon Ice Cap, and Barnes Ice Cap mapped using Sentinel-2 imagery, Int. J. Appl. Earth Obs., 78, 1–13,, 2019a. 

Yang, K., Smith, L. C., Fettweis, X., Gleason, C. J., Lu, Y., and Li, M.: Surface meltwater runoff on the Greenland ice sheet estimated from remotely sensed supraglacial lake infilling rate, Remote Sens. Environ., 234, 111459,, 2019b.  

Yang, K., Smith, L. C., Cooper, M. G., Pitcher, L. H., Van As, D., Lu, Y., Lu, X., and Li, M.: Seasonal evolution of supraglacial lakes and rivers on the southwest Greenland Ice Sheet, J. Glaciol., 67, 592–602,, 2021. 

Zhang, Q., Huai, B., van Den Broeke, M. R., Cappelen, J., Ding, M., Wang, Y., and Sun, W.: Temporal and Spatial Variability in Contemporary Greenland Warming (1958–2020), J. Climate, 35, 2755–2767,, 2022. 

Zwally, H. J., Abdalati, W., Herring, T., Larson, K., Saba, J., and Steffen, K.: Surface melt-induced acceleration of Greenland ice-sheet flow, Science, 297, 218–222,, 2002. 

Short summary
We map and quantify surface rivers and lakes at Humboldt Glacier to examine seasonal evolution and provide new insights of network configuration and behaviour. A widespread supraglacial drainage network exists, expanding up the glacier as seasonal runoff increases. Large interannual variability affects the areal extent of this network, controlled by high- vs. low-melt years, with late summer network persistence likely preconditioning the surface for earlier drainage activity the following year.