Diverse supraglacial drainage patterns on the Devon ice Cap, Arctic Canada

ABSTRACT

The Devon Ice Cap (DIC) is one of the largest ice masses in the Canadian Arctic. Each summer, extensive supraglacial river networks develop on the DIC surface and route large volumes of meltwater from ice caps to the ocean. Mapping their extent and understanding their temporal evolution are important for validating runoff routing and melt volumes predicted by regional climate models (RCMs). We use 10 m Sentinel-2 images captured on 28 July and 10/11 August 2016 to map supraglacial rivers across the entire DIC (12,100 km²). Both dendritic and parallel supraglacial drainage patterns are found, with a total length of 44,941 km and a mean drainage density (D_d) of 3.71 km⁻¹. As the melt season progresses, D_d increases and supraglacial rivers form at progressively higher elevations. There is a positive correlation between RCM-derived surface runoff and satellite-mapped D_d, suggesting that supraglacial drainage density is primarily controlled by surface runoff.

KEYWORDS:

• Supraglacial rivers
• regional climate model (RCM)
• Sentinel-2
• Devon ice Cap
• Arctic

1. Introduction

The Arctic contains numerous large ice masses including the Greenland Ice Sheet, Devon Ice Cap (DIC) and Agassiz Ice Cap, with a total collective area of approximately two million square kilometers (Abdalati et al., 2004; Cook et al., 2019; Gardner et al., 2011). The Arctic is also highly sensitive to global climate change (Screen & Simmonds, 2010), and is currently warming at roughly twice the global average rate, with annual air temperatures for October 2015 to September 2016 averaging +2.0°C warmer than the 1981–2010 average (Richter-Menge et al., 2016). Associated mass loss from Arctic ice masses contributed ∼31% of observed global sea level rise between 1992 and 2017 (Box et al., 2018), with annual ice mass losses averaging 50.2 Gt a⁻¹ in the Canadian Arctic, over twice the pre-1996 average (Noël et al., 2018). It is therefore essential to understand factors that affect the prevalence, spatial pattern, and transport of surface meltwater from Arctic ice masses to oceans.

During 2016, the warmest year in the Arctic since 1900 (Richter-Menge et al., 2016), Arctic ice masses experienced significant surface melt, which formed large and complex supraglacial river networks at their lower elevations in summer (Pitcher & Smith, 2019). Where supraglacial rivers terminate in moulins, large volumes of surface meltwater are transferred to the ice bed, influencing the evolution of subglacial hydrological systems and ice flow dynamics (Chu, 2014; Iken & Bindschadler, 1986; Pitcher & Smith, 2019; Röthlisberger & Iken, 1981; Schoof, 2010; Sundal et al., 2011). Although the spatial structure of these supraglacial river networks is thought to be controlled by the surface topography (Yang et al., 2015a) through the transfer of bed variability to the ice surface at the large scale (Crozier et al., 2018; Ignéczi et al., 2018), their temporal evolution is affected by the rate and magnitude of meltwater production and the permeability of glacier surface throughout the melt season (Irvine-Fynn et al., 2011; Lampkin & Vanderberg, 2014).

In general, large supraglacial river networks have only recently begun attracting scientific attention with most studies focused on the Greenland Ice Sheet (Colgan et al., 2011; Crozier et al., 2018; Gleason et al., 2016; Ignéczi et al., 2018; King et al., 2016; Lampkin & Vanderberg, 2014; Legleiter et al., 2014; Pitcher & Smith, 2019; Poinar et al., 2015; Smith et al., 2015; Smith et al., 2017; Yang & Smith, 2016). Yang et al. (2019a) mapped small areas of the Devon and Barnes Ice Caps (as well as the northwest Greenland Ice Sheet) using 10 m Sentinel-2 visible/near-infrared satellite images and found that the drainage patterns and densities of supraglacial river networks may differ significantly within and between different ice masses.

Previous studies have documented supraglacial river formation on the DIC (Fernandes et al., 2018; Wyatt, 2013; Wyatt & Sharp, 2015; Yang et al., 2019a). Wyatt (2013) manually digitized supraglacial rivers across the DIC from Landsat 7 ETM+ imagery acquired during the 1999 summer, revealing that both surface topography and ice dynamics affect supraglacial drainage patterns. Yang et al. (2019a) mapped a snapshot of supraglacial rivers on the western and southern DIC, demonstrating the strong potential of 10 m Sentinel-2 imagery to investigate supraglacial drainage patterns. However, the spatial distribution, drainage patterns and seasonal evolution of supraglacial rivers across the entire DIC and how their formation is affected by surface runoff remain poorly understood.

Here, we build upon the preliminary study of Yang et al. (2019a) using Sentinel-2 satellite images to map supraglacial river networks and drainage across the entire DIC, their mid- to late-season evolution, and response to surface runoff. First, we extract the supraglacial rivers from satellite images by integrating cross-sectional and longitudinal open-channel morphometry. Then, we calculate the supraglacial drainage density within different elevation bands and drainage basins to investigate the evolution of supraglacial rivers across western, eastern and southern DIC. Finally, we compare the satellite-mapped supraglacial drainage density (D_d) with ice surface runoff modeled from Modèle Atmosphérique Régional (MAR), a regional climate model (RCM), which reveals how surface runoff affects the growth of the supraglacial rivers. We conclude by discussing the implications of the observed supraglacial drainage patterns and the relationship between satellite-mapped D_d and RCM-derived surface runoff estimates.

2. Study area

The DIC is located between 74.5°N and 75.8°N and 80.0°W and 86.0°W in the Canadian Arctic (Figure 1). Its total area is 12,100 km², with an estimated ice volume of 3,980 km³ and maximum ice thickness of 880 m (Dowdeswell et al., 2004). Climate warming has greatly impacted the DIC in recent decades. Increases in summertime temperature and lengthening melt seasons have intensified surface ablation and overall mass loss (Boon et al., 2010; Burgess et al., 2005; Burgess & Sharp, 2004; Burgess & Sharp, 2008; De Jong et al., 2018; Gascon et al., 2013a; Shepherd et al., 2007). Both the total volume and area of the DIC are declining (∼−1.7 km³ a⁻¹; ∼−8.3 km² a⁻¹_, respectively) (Boon et al., 2010; Burgess & Sharp, 2004; Mortimer et al., 2016), and its annual mass loss was six times higher from 2005–2009 than from 1963–2005 (Sharp et al., 2011).

Figure 1. Location of the Devon Ice Cap in the Canadian Arctic.

The main ice cap has a dome-like shape with a maximum elevation of 1921m (Dowdeswell et al., 2004). Eight major drainage basins (471–2,623 km²), obtained from the Randolph Glacier Inventory 6.0 (RGI Consortium, 2017), cover 78% of the total area of DIC (Figure 2) (Table 2). There are strong spatial gradients in mass balance, glacial geometry and ice flow across the ice cap (Burgess & Sharp, 2004; Mair et al., 2005), causing different DIC drainage basins to exhibit distinct melt characteristics and hydrologic conditions (Burgess & Sharp, 2004).

Figure 2. Sentinel-2 images (acquired on 28 July 2016, RGB: bands 4 (Red), 3 (Green), 2 (Blue)) of the entire Devon Ice Cap (excluding the southwest arm). (a–f) show typical areas in different DIC basins (red rectangles in center panel (g)). Black lines mark the boundaries of the major drainage basins (numbered as in Table 2).

Because of the lack of cloud-free image coverage, we only investigated the mid- to late-season evolution of supraglacial rivers on the western, southern and eastern DIC. Only the low-elevation area (below 800 m) on the eastern DIC (covering 5% of eastern DIC and 2% of these three areas) contains crevasses observed by the Sentinel-2 imagery used herein.

3. Data and methods

3.1. Sentinel-2 Multi-Spectral Instrument images

Landsat images have been used for extracting supraglacial rivers on both the Greenland and Antarctic ice sheets (Bell et al., 2017; Kingslake et al., 2017; Lampkin & Vanderberg, 2014; Poinar et al., 2015; Yang & Smith, 2016). However, because supraglacial rivers are normally less than 30 m wide (the spatial resolution of multi-spectral Landsat images) these data can delineate large main stems but miss narrow tributaries. To extract more representative supraglacial river networks, higher spatial resolution images are required (Colgan et al., 2011; Karlstrom & Yang, 2016; Smith et al., 2015; Smith et al., 2017; Yang & Smith, 2013). Ultra high resolution images (meter-scale), such as GeoEye and WorldView-1/2/3/4, can allow accurate delineation of supraglacial rivers, but their high cost and limited spatial coverage reduce their usability in large-scale applications (Yang & Smith, 2013).

As demonstrated by Yang et al. (2019a), the Sentinel-2A/B satellites (hereafter ‘Sentinel-2’) launched by the European Space Agency (ESA) in 2015 and 2017 provide an exciting new opportunity to map supraglacial rivers over large geographic areas with a spatial resolution intermediate between Landsat and WorldView. Sentinel-2 has four spectral bands (blue, green, red and near-infrared) at 10 m spatial resolution. This multi-spectral capability enables calculation of the normalized difference water index (NDWI), which can differentiate surface water bodies from surrounding glacial ice and snow. We use the Sentinel-2 Level-1C product, (orthorectified top-of-atmosphere reflectance) with sub-pixel multispectral registration (Baillarin et al., 2012), provided freely by ESA’s Scientific Data Hub (https://scihub.copernicus.eu/dhus/#/home).

Nine Sentinel-2 images were acquired to extract supraglacial rivers across the DIC from late July through early August, the mid- to late- melt season for Arctic ice masses in this area. These images include six from 28 July 2016 (Figure 2), two from 10 August 2016, and one from 11 August 2016 (Table 1). Elevation contours were generated from the 30 Arc-Second Elevation (GTOPO30) digital elevation model, provided by the USGS EarthExplorer (https://earthexplorer.usgs.gov/). The 2 m ArcticDEM, obtained from Polar Geospatial Center (https://www.pgc.umn.edu/data/arcticdem/), was used to simulate supraglacial drainage patterns and compare with Sentinel-2 mapped supraglacial rivers by using the ArcGIS hydrologic analysis tools.

Table 1. Images used in this study.

Image Type	Acquired Date	Tile Number	Data Source
Sentinel-2A MSI	2016.7.28	T16XEH	ESA
	2016.7.28	T16XEJ
	2016.7.28	T17XMC
	2016.7.28	T16XMD
	2016.7.28	T16XME
	2016.7.28	T16XND
	2016.8.10	T16XEJ
	2016.8.10	T17XMD
	2016.8.11	T17XMD
GTOPO30	-	W100N90	USGS
ArcticDEM	-	27_30	PGC
		27_31
		28_29
		28_30
		28_31

Table 2. Supraglacial hydrological features of the major drainage basins across the DIC on 28 July 2016.

Basin ID	Area (km^2)	Drainage density (km^−1)	RCM-derived surface runoff (mm)
1	1046.46	5.71	9.33
2	765.44	1.63	5.66
3	663.92	2.51	8.26
4	1134.44	2.06	12.03
5	572.17	3.09	13.97
6	471.08	2.46	10.62
7	2622.81	1.94	10.42
8	2170.55	2.19	6.62

3.2. Supraglacial river extraction

The complex physical and spectral characteristics of bare ice make automated extraction of supraglacial rivers from multi-spectral satellite images challenging (Yang & Smith, 2013). Here, we used the method of Yang et al. (2015b) to extract rivers on the basis of their Gaussian-like brightness cross-section and longitudinal open-channel morphometry as follows: First, water features were enhanced by applying NDWI (Mcfeeters, 1996), which differentiates meltwater from the background ice or snow (Figure 3b). Second, low-frequency image background and high-frequency image noise were removed using a band-pass filter ramped between 1/200 m⁻¹ and 1/40 m⁻¹ (Figure 3c). Then, a Gabor filter was applied and the thickness of the Gabor filter was set to 2 to enhance small rivers (river widths smaller than two pixels) (Yang et al., 2015b). The Gabor filter can enhance linear features in different directions and at different scales, which effectively amplifies the river cross-sections (Figure 3d). To improve river connectivity (Figure 3e) by accounting for brightness variations along the river courses and issues with disconnected river segments, a path opening operator was then applied (Heijmans et al., 2005). Path opening is a flexible mathematical morphology operator, which homogenizes the brightness of pixels along linear features in the image by designing a slender, directional path operator with a fixed distance (L_min= 20 pixels) (Figure 3e).

Figure 3. Supraglacial river extraction in Sentinel-2 image acquired on 28 July 2016, following the method of (Yang et al., 2015b). (a) Original Sentinel-2 true color image, (b) NDWI image, (c) NDWI image after band-pass filtering, (d) image after Gabor filtering, (e) image after path opening and (f) final supraglacial river extraction results and manual digitization of the same image.

After Gabor filtering and path opening, the supraglacial rivers become far easier to differentiate from the surrounding ice background. Using a global threshold (T = 5) based on the pixel brightness histogram of the grayscale image, supraglacial rivers are effectively discerned and used to create a binary water mask (Yang et al., 2015b). Finally, ArcScan for ArcGIS was used to generate river vectors from this binary water mask (Figure 3f). To quantify the ability of the automated mapping technique to capture supraglacial rivers, we also manually digitized the supraglacial rivers shown in Figure 3a. The overall spatial pattern was very similar, but the sum length of the automated extraction results is 13.6% shorter than the manual digitization results, due to several gaps in automatically-mapped tributary rivers.

At high elevations, surface meltwater initially percolates into snowpack pore spaces and then forms slush flows when saturation is reached (Irvine-Fynn et al., 2011; Pitcher & Smith, 2019). These slush flows can transform into open meltwater channels once the bare ice surface is exposed (Chu, 2014; Cuffey & Paterson, 2010; Pitcher & Smith, 2019). It is challenging to distinguish slush flows from well-channelized supraglacial rivers in satellite imagery since they are both represented as bright, linear features in NDWI images. Therefore, in this study, we consider these slush flows as initial supraglacial rivers and do not distinguish between slush flows and open meltwater channels.

While we acknowledge that our optical satellite-based study may be missing some key physical processes (e.g. firn stratigraphy) not visible from space, our study does offer an efficient method for mapping temporal and spatial variations in supraglacial river networks from space.

3.3. RCM-derived surface runoff

The Modèle Atmosphérique Régional (MAR) is among the best RCM over the Greenland Ice Sheet (Fettweis et al., 2017; Fettweis et al., 2020). MAR was used to model surface runoff across the DIC by simulating temperature, precipitation, and refreezing variables in 7.5 × 7.5 km grid cells. MAR version 3.8 was used in this study and forced at the lateral boundaries by 6-hourly ERA-Interim reanalysis outputs (Delhasse et al., 2018). MAR integrates a runoff delay function to represent the time lag between surface meltwater production and transportation to nearby supraglacial rivers (Lefebre et al., 2003; Zuo & Oerlemans, 1996). However, the MAR delay function temporally-smooths the MAR runoff signal resulting in less realistic results than the runoff calculated directly from the raw model data (Runoff = Melt + Rain - Refreezing) (Yang et al., 2019b). The uncertainty is about 30% of runoff value (Fettweis et al., 2017). In this study, we therefore calculated ice surface (runoff) directly from MAR on 28 July 2016 and 10/11 August 2016, to yield the total amount of daily runoff generated within each grid cell (mm day⁻¹) based on surface mass balance variables alone with no routing delay imposed. MAR grid cells were intersected with different vector polygon boundaries (elevation bands, drainage basins) and the total runoff computed by summing the intersected runoff values.

Drainage density (D_d) was calculated by dividing the total river length by the area of the corresponding basin/elevation band (Yang et al., 2016). In section 4.3, we compare satellite-mapped D_d with RCM-derived surface runoff above 600 m, because the areas of elevation bands near the ice cap margin (below 600 m) are small and supraglacial D_d below 600 m is influenced by multiple factors (e.g. crevasses and dust zone) rather than only surface runoff (Lea & Brough, 2019).

4. Results

4.1. Mapping of supraglacial rivers across the Devon Ice Cap

Supraglacial rivers were successfully mapped across the entire DIC using nine Sentinel-2 L1C images acquired on 28 July 2016 (Main Map, Figure 4). Diverse supraglacial drainage patterns and densities were found in different DIC basins. In total, 5.5% of the DIC was covered by supraglacial rivers and only 0.1% of the DIC was covered by supraglacial lakes at this time. The total mapped river length was 4.5×10⁴ km and the average D_d was 3.71 km⁻¹. The highest elevation at which supraglacial rivers were mapped was ∼1500 m.

Figure 4. Supraglacial river networks across the entire Devon Ice Cap mapped from 2 m ArcticDEM and six 10 m Sentinel-2 images (acquired on 28 July 2016), for everywhere except the bottom right (southeastern DIC) due to the lack of Sentinel-2 images coverage. (a–f) show examples of supraglacial river networks in different DIC basins (red rectangles in center panel (g)).

The flat western basin of the DIC (Basin 1) was found to have the highest D_d (5.71 km⁻¹), with well-developed supraglacial river networks (Figure 4a). At elevations between 1100 and 1300 m, meltwater is not clearly channelized and represents braided and broadly dendritic supraglacial rivers. As the meltwater flows to lower elevations, the supraglacial rivers become entrenched and comprise of many small rivers parallel to each other (Figure 4a). At elevations below 1000 m, supraglacial drainage becomes organized into fewer large rivers and D_d decreases (Figure 4g, Figure S1).

On the northwest DIC (Basin 2), with the lowest D_d (1.63 km⁻¹), most supraglacial rivers converge into several large supraglacial rivers at the elevation between 1100 and 1400 m, then continue to flow down the ice surface to lower elevation (rather than terminate in supraglacial lakes or moulins) (Figure 4b). Meanwhile, some small sub-parallel rivers, which initiate at ∼1000 m, are distributed from ∼1000 m to ∼600 m near the ice cap margin (Figure 4b, g).

On the northeastern DIC (Basin 3–6), D_d values range from 2.06–3.06 km⁻¹. Supraglacial rivers initiate at ∼1400 m and form dendritic drainage patterns on the main cap, transforming into subparallel arrangements on outlet glaciers at lower elevations, with meltwater routed parallel to the outlet glacier trunk (Figure 4d).

Basins 7 and 8, on the southern DIC, are the two largest drainage basins (2,171 km² and 2,623 km², respectively), and exhibit relatively low D_d (1.94 and 2.19 km⁻¹, respectively). Most supraglacial rivers are dendritic (Figures 4c, e), apart from the western portion of Basin 8, where supraglacial rivers are sub-parallel and have low sinuosity (Figure 4f). Two small supraglacial lakes (area 0.03 and 0.08 km²) were located at ∼1400 m and ∼1250 m in Basin 7 (Figure 4c) and one supraglacial lake (area 0.06 km²) was located at ∼1200 m in Basin 8 (Figure 4e). Supraglacial rivers at the highest elevation were located in Basin 8 at 1500 m, matching previous observations from the 2015 melt season (Fernandes et al., 2018).

Comparing satellite-mapped supraglacial rivers with ArcticDEM-modeled supraglacial rivers, we find that the satellite-mapped supraglacial drainage patterns are similar to ArcticDEM-modeled drainage patterns, which implies that the supraglacial drainage patterns are strongly controlled by surface topography (Figure 4). However, ArcticDEM-modeled drainage networks do not reflect the true distribution of supraglacial rivers. Instead, they tend to overestimate high-elevation rivers (which may be channelized but dry at the time of image capture) and underestimate low-elevation rivers (due to the choice of the accumulation area threshold) (Figure 4).

4.2. Mid- to late-season evolution of supraglacial rivers

Mid- to late-season evolution of supraglacial rivers on the western, eastern, and southern DIC (parts of Basins 1, 7 and 8) was analyzed by comparing the supraglacial drainage densities on 28 July and 10/11 August 2016. Between these dates, D_d of these three areas increased: western DIC increasing from 4.92 km⁻¹–5.99 km⁻¹; southern DIC from 0.98 km⁻¹–3.18 km⁻¹; and eastern DIC from 0.80 km⁻¹–9.59 km⁻¹ (Figure 5). As melting continued, D_d varied markedly with elevation and supraglacial rivers gradually extended toward higher elevations. The snow line increased from ∼1200 m to ∼1500 m during 28 July and 10/11 August 2016 (Figure 6), and the daily runoff was clearly higher on 10/11 August than that on 28 July at 1300–1400 m (Figure 7), which coincided with the largest increase in supraglacial D_d occurred.

Figure 5. Supraglacial rivers evolution. (a–c) show the supraglacial river networks extraction results, and (d–f) are corresponding supraglacial drainage density in different elevation bands.

Figure 6. Snow cover on the western, southern and eastern DIC delineated by thresholding the Sentinel-2 NIR bands (Ryan et al., 2019). The red lines show the boundaries of the snow and bare ice on 28 July and 10/11 August 2016.

Figure 7. Daily RCM-derived surface runoff for 1300–1400 m elevation band of western, southern and eastern DIC during the 2016 melt season.

From 28 July to 10 August 2016, the total D_d of supraglacial rivers on the western DIC (a part of Basin 1) increased slightly from 4.92–5.99 km⁻¹ and the snow covered fraction decreased from 0.39–0.14 (Figures 5a, d, Figure 6). Supraglacial rivers expanded to higher elevations with the peak value of D_d moving from 1200 m to 1400 m (Figure 5d). By the 10 August, snow-covered areas had become slush fields at 1200–1500 m (Pitcher & Smith, 2019), with many poor-channelized supraglacial rivers with low sinuosity forming (Figures 5a, Figure 6), with a maximum D_d of 11.62 km⁻¹ at this elevation band. At 1100–1200 m, the slush field transitioned into relatively impermeable bare ice surface from 28 July to 10 August (Irvine-Fynn et al., 2011), the supraglacial rivers were better channelized and the D_d decreased from 11.32–5.66 km⁻¹. Below 1100 m, the D_d decreased by more than 50% (from 8.23–2.27 km⁻¹). Although the main rivers remained relatively stable across this period, many small supraglacial rivers had either dried up or became narrow, and thus could not be distinguished by Sentinel-2 images (Figures 5a, d).

On the southern and eastern DIC, the evolution of supraglacial rivers differed from that on the western DIC. On the southern DIC (a part of Basin 8), D_d was generally lower, and peak D_d was located at 1300 m. As the snow surface diminished at 1100-1400 m between 28 July and 10 August, a large number of initial supraglacial rivers formed and the D_d increased by a factor of three (from 0.98–3.18 km⁻¹) (Figure 5b, e, Figure 6). This trend was also clearly observed on the eastern DIC (a part of Basin 7), where D_d on 11 August was 12 times higher than that on 28 July (from 1.14–21.29 km⁻¹) (Figure 5c, f) and the snow cover fraction decreased from 0.40–0.05 (Figure 6). The location of the peak D_d extended from 1200 m to 1300 m (Figures 5c, f). Overall, on the southern and eastern DIC, a large number of newly formed supraglacial rivers were observed at both high and low elevations in contrast to the decrease in D_d at low elevations (<1200 m) observed on the western DIC.

4.3. Comparing satellite-mapped D_d with RCM-derived surface runoff

There is a significantly positive linear relationship between the RCM-derived surface runoff (R) and the satellite-mapped D_d, calculated at a range of elevations on 28 July 2016 (Figure 8). Across the entire DIC, the relationship between R and D_d can be expressed by the linear function D_d = 0.20 R + 0.05 (r² = 0.90, p < 0.01). In other words, D_d increases by 0.2 km⁻¹ with each 1 mm increase in R. Both R and D_d gradually increase when moving from the high elevation interior (> 1700m) to lower elevation areas (∼700 m) near the margin of the DIC (Figure 8).

Figure 8. Relationship between daily RCM-derived surface runoff (R) and satellite-mapped drainage density (D_d) on 28 July 2016. All R and D_d were calculated at each 100 m of elevation (from 700 m to 2000m).

Different drainage basins demonstrate variations in the relationship between R and D_d (Figure 9). In Basins 2–8, significantly positive linear relationships (R² > 0.68, p < 0.01) exist between R and D_d (Figures 9b–h). With each 1 mm increase in R, the D_d increases by 0.13–0.37 km⁻¹ in Basin 2–8. However, in Basin 1, where D_d is highest, no direct relationship is apparent between R and D_d (Figure 9a), due to the significant decrease in D_d from ∼1000 m to the ice cap margin.

Figure 9. Relationship between daily RCM-derived surface runoff (R) and satellite-mapped drainage density (D_d) in eight major basins on 28 July 2016. All R and D_d were calculated at each 100 m of elevation (from 700 m to 2000m).

In Basin 3, as R increases, D_d at elevations above 900 m rise faster than at elevations below 900 m (0.52 km⁻¹ mm⁻¹ vs. 0.15 km⁻¹ mm⁻¹) (Figure 9c). Compared to Basin 2 and Basins 4–8, each 1 mm increase in R of Basin 3 causes a greater increase in D_d (0.37 km⁻¹ vs. 0.13–0.24 km⁻¹) (Figure 9c). In Basin 3, the area below 900 m comprises a topographically-constrained outlet glacier approximately 25 km in length, where surface water may be concentrated. Nearly all supraglacial rivers drain into this outlet glacier rather than terminating on land or in supraglacial lakes and moulins. This is in contrast to outlet glaciers in other basins, which only exist below 600 m (Boon et al., 2010) and are therefore not included in the statistics (Figure 2g).

5. Discussion

The overall DIC supraglacial drainage pattern is characterized by high density, dendritic and parallel drainage networks. On both the eastern and southern DIC, dendritic drainage patterns dominate, which is in agreement with patterns observed by previous studies (Dowdeswell et al., 2004; Fernandes et al., 2018; Wyatt, 2013; Wyatt & Sharp, 2015). However, only a few studies have focused on the western DIC (Basin 1) (Wyatt, 2013; Yang et al., 2019a), where we find parallel and dendritic drainage patterns and the highest D_d. Meanwhile, there are significant decreases in D_d at low elevations on western DIC (Basin 1). Ignéczi et al. (2018) suggested that river networks on western Greenland close to the ice margin become fragmented due to crevasses and moulins, and D_d would reduce, resulting in a corresponding reduction in D_d. However, due to the few crevasses that we observed on the western DIC, we attribute the low drainage density mainly to the difficulty in detecting supraglacial rivers in the low albedo zone (Yang et al., 2019a). During the melt season, supraglacial rivers extend toward higher elevations on the DIC as the surface runoff increases at high elevation, similar to the evolution of supraglacial rivers on the southwest Greenland Ice Sheet (Lampkin & Vanderberg, 2014).

We have compared supraglacial rivers delineated in this study with supraglacial rivers delineated in Wyatt (2013). In general, our 10 m Sentinel-2 river map matches well. Notably, however, our Sentinel-2 map delineates numerous small supraglacial rivers that are not captured by Wyatt’s map and yields higher D_d (maximum: 12.5 km⁻¹ vs. 8.7 km⁻¹). This agrees with Yang et al. (2019a) who show that Sentinel-2 images are superior to Landsat images for delineating small supraglacial rivers. In addition, the highest supraglacial D_d is found on the western DIC (Basin 1), but was reported on the eastern DIC by Wyatt (2013). The highest-elevation of satellite-mapped supraglacial rivers is located at 1500 m in Basin 7, which is 80 m higher than that reported by Wyatt (2013) during the 1999 summer. It is probably that with climate change, the occurrence of the further inland melting and the change of the underlying surface (ice/snow/firn) leads to the supraglacial rivers’ expansion (Bezeau et al., 2013; Fisher et al., 2012; Gascon et al., 2013b; Rutishauser et al., 2016; Sharp et al., 2011; Wolken et al., 2009).

We present a positive linear relationship between ice surface runoff and supraglacial D_d, which has not been quantified previously. The slope of the linear function between satellite-mapped D_d and RCM-derived surface runoff in Basin 3 is higher than in other basins, indicating that RCM may underestimate the surface runoff on outlet glaciers. Low-albedo exposed rocks that radiate longwave radiation may increase local melt rates (Kingslake et al., 2017). Further improving parameters of RCM is important for better illustrating the response of supraglacial rivers to surface runoff. Moreover, the positive linear relationship between surface runoff and supraglacial D_d indicate that supraglacial D_d is primarily controlled by surface runoff. We suggest that this finding is important for better quantifying surface runoff. More specifically, runoff is the most difficult variable to qualify in RCMs and is usually calculated as residuals of other climatic variables (Smith et al., 2017); our study provides a potential approach to validate RCM-derived surface runoff using satellite-observed D_d.

6. Conclusion

This paper presents a first complete time-slice map of > 10 m wide supraglacial river networks across the entire Devon Ice Cap using high spatial resolution (10 m) Sentinel-2 L1C satellite images. On 28 July 2016, supraglacial drainage density on the DIC was 3.71 km⁻¹, with both dendritic and parallel drainage patterns evident in different drainage basins of the ice cap, which will improve our understanding of the intensity and efficiency at which meltwater is being transported from the ice cap. The basin with the highest drainage density was located on the western DIC, where peak drainage density was found at elevations between 1200–1400 m. By August 10/11, a large number of new supraglacial rivers had formed at higher elevations, and the region of peak drainage density had moved toward higher elevations. At low elevations (< 1200 m), drainage density decreased on the western DIC but increased in both the southern and eastern DIC. Moreover, we found a positive relationship between satellite-mapped supraglacial drainage density and RCM-derived surface runoff. This indicates that the supraglacial drainage density is primarily controlled by ice surface runoff, which will improve our understanding of supraglacial hydrological processes and mass balance across the Arctic cryosphere. In the future, integrating remotely sensed observations, RCM outputs and high spatial resolution DEMs (e.g. 2 m ArcticDEM) more closely will enable a better understanding of the supraglacial drainage patterns and their response to ice surface runoff.

Software

ESRI ArcGIS 10.6.1 and Adobe Illustrator CC 2018 were used in the production of the main map.

Acknowledgements

Kang Yang acknowledges the support of the National Key R&D Program (2018YFC1406101), the National Natural Science Foundation of China (41871327), and the Fundamental Research Funds for the Central Universities. Laurence Smith acknowledges the NASA Cryospheric Science Program (80NSSC19K0942). Stephen Livingstone and Andrew Sole acknowledge the support of a White Rose University Consortium Collaboration Fund grant. We acknowledge Faye Wyatt for providing a supraglacial river map of the DIC delineated from Landsat 7 ETM+ imagery acquired during 1999 summer. We thank the Associate Editor and four reviewers (Martin Sharp, Caroline Clason, David Burgess and Daniele Cannatella) for their detailed and constructive comments.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the National Key R&D Program [Grant Number 2018YFC1406101], the National Natural Science Foundation of China [Grant Number 41871327], and the Fundamental Research Funds for the Central Universities [Grant Number 0209-14380070].