1. Introduction
The Arctic is undergoing rapid warming (Rantanen and others, Reference Rantanen2022), which has been particularly pronounced in the Canadian Arctic Archipelago (Sharp and others, Reference Sharp, Kargel, Leonard, Bishop, Kääb and Raup2014). From 2000 to 2023, the Canadian Arctic Archipelago experienced the highest rate of glacier mass loss after Alaska, averaging 53.6 ± 4.7 Gt a−1 (Zemp and others, Reference Zemp2025). Surface melt is the primary driver of this mass loss, accounting for 90% of the total mass losses in the northern Canadian Arctic Archipelago since 2005, leading to the formation and expansion of supraglacial hydrological networks (Millan and others, Reference Millan, Mouginot and Rignot2017; Lu and others, Reference Lu2020). Supraglacial meltwater channels play a crucial role in a glacier’s hydrological system by influencing surface mass balance (SMB) and routing meltwater that can eventually access englacial and subglacial pathways, thereby serving as an important connection between the surface and internal drainage systems (Pitcher and Smith, Reference Pitcher and Smith2019). When supraglacial channels terminate in sinks such as moulins and large crevasses, large volumes of water can transfer to the glacier bed, impacting basal water pressure and the development of subglacial hydrological systems (Lu and others, Reference Lu2020), in turn impacting ice motion (Schoof, Reference Schoof2010).
Advances in remote sensing (Yang and others, Reference Yang2019; Bash and others, Reference Bash2023), hydrological modeling (King and others, Reference King, Hassan, Yang and Flowers2016) and field studies (St Germain and Moorman, Reference St Germain and Moorman2019) have significantly enhanced our understanding of supraglacial hydrological networks. Despite these advances, supraglacial systems remain among the least studied hydrological components on Earth (Yang and others, Reference Yang2019). Recent studies have begun to bridge this gap through the delineation of supraglacial hydrological features, which has been facilitated by the availability of high-resolution (e.g. WorldView-2 [WV2], 0.5 m; Yang and others, Reference Yang, Li, Liu, Cheng, Huang and Chen2015a) to moderate-resolution (e.g. Sentinel-2, 10 m; Yang and others, Reference Yang, Karlstrom, Smith and Li2017; Lu and others, Reference Lu2020) satellite imagery and digital elevation models (DEMs) such as the ArcticDEM (Rippin and others, Reference Rippin, Pomfret and King2015; King and others, Reference King, Hassan, Yang and Flowers2016; Yang and others, Reference Yang2018; Bash and others, Reference Bash2023). Recent studies in the Canadian Arctic Archipelago have contributed to the development or refinement of methodological frameworks for supraglacial hydrology mapping (Yang and others, Reference Yang2019; Bash and others, Reference Bash2023), assessed the evolution of the supraglacial drainage system over a single melt season (Lu and others, Reference Lu2020), or have explored the relationship between surface melt and glacier dynamics for single points in time (Wyatt and Sharp, Reference Wyatt and Sharp2015; Yang and others, Reference Yang2019). Despite these advances, research on supraglacial channels in the Canadian Arctic Archipelago remains limited, with most large-scale studies having taken place on the Greenland Ice Sheet (Yang and Smith, Reference Yang and Smith2013; Yang and others, Reference Yang, Li, Liu, Cheng, Huang and Chen2015a; King and others, Reference King, Hassan, Yang and Flowers2016; Yang and others, Reference Yang, Karlstrom, Smith and Li2017; Yang and others, Reference Yang2019; Lu and others, Reference Lu2020).
In the Canadian Arctic Archipelago, St Germain and Moorman (Reference St Germain and Moorman2019) made initial strides in addressing the gap in knowledge of supraglacial hydrological networks by studying the evolution of specific streams and rivers on Fountain Glacier, Bylot Island, although their research focused on a relatively short period (2010–17), capturing only a snapshot of their evolution. Although Kochtitzky and others (Reference Kochtitzky2023) mapped the development of supraglacial channels over the last ∼60 years on Thores Glacier, northern Ellesmere Island, their study focused solely on major rivers in the ablation area and so did not capture fluctuations across the entire supraglacial system.
Recent findings from Greenland show that supraglacial drainage systems are highly dynamic, responding strongly to interannual variations in surface melt and the seasonal retreat of the snowline to higher elevations, which drives up-glacier expansion of drainage networks (Rawlins and others, Reference Rawlins, Rippin, Sole, Livingstone and Yang2023; Glen and others, Reference Glen2025). Network morphology also evolves with time: changes in discharge influence channel incision and the transition from ephemeral surface streams to more deeply incised rivers and canyons. When combined with local glacier geometry, this will govern whether channels become perennial features or remain seasonally transient (St Germain and Moorman, Reference St Germain and Moorman2019).
In this study, we examine how long-term climatic changes, particularly those influencing SMB and the position of the Equilibrium Line Altitude (ELA), have impacted the development of supraglacial hydrological features in the northern Canadian Arctic Archipelago. We aim to determine whether multi-decadal trends reveal up-glacier expansion of drainage networks and an increasing tendency for supraglacial channels to become perennial, rather than remaining primarily governed by annual variations in melt intensity. To address this, we employ a combination of optical satellite imagery, DEMs and historical aerial photographs to manually map supraglacial channel networks and assess their multi-decadal evolution across selected glaciers on Ellesmere Island in the northern Canadian Arctic Archipelago. We quantify their temporal changes over approximately 60 years by analyzing shifts in supraglacial drainage patterns and the composition of different channel types within the drainage system. We also evaluate spatial differences among glaciers to better understand the magnitude of change in supraglacial systems located in varied environmental settings. Together, these analyses provide a multi-decadal quantification of the supraglacial hydrological networks’ evolution in this region, contributing novel insights into their response to climate warming.
2. Study area
Ellesmere Island (Umingmak Nuna), Nunavut (Fig. 1a), is the largest island in the northern Canadian Arctic Archipelago. Approximately 40% (∼77 516 km2) of its area is covered by glaciers, of which ∼53% are land-terminating and ∼47% are marine-terminating (Millan and others, Reference Millan, Mouginot and Rignot2017; White and Copland, Reference White and Copland2018). The island experiences a pronounced air temperature gradient from warmer mean annual temperatures in the south to colder temperatures in the north. There is also a precipitation gradient extending from southeast to northwest, with precipitation highest along the coastlines facing Baffin Bay and decreasing with elevation and distance inland (Koerner, Reference Koerner1979).
Figure 1. Location of the main icefields and ice caps, as well as the five studied glaciers, on Ellesmere Island. Projection: Canada Lambert Conformal Conic (main map), WGS UTM 16N (Sydkap Glacier) and WGS UTM 18N (other glaciers). Data: Panel (a)—Statistics Canada, 2016 Census—Provinces/territories—Cartographic Boundary File (statcan.gc.ca). Natural Resources Canada—Lakes, Rivers and Glaciers in Canada—CanVec Series—Hydrographic Features (Open Government Portal). Panels (b–f) —Glacier outlines are manually corrected Randolph Glacier Inventory v7 outlines. Satellite imagery: PlanetScope for panels (b, d, f) and Sentinel-2 for panels (c, e).
The current study focuses on five glaciers (Fig. 1b–f), which were selected based on i) geographical distribution across the ∼830 km latitudinal extent of Ellesmere Island, ii) data availability and iii) the presence of a perennial supraglacial drainage system (except for Sydkap Glacier, which was included for comparison). From north to south, these are:
1. Unnamed 1 Glacier (82.998°N, 72.855°W; Fig. 1b) is the smallest of the studied glaciers, with a basin area of approximately 18 km2 and elevations ranging from ∼30 to 980 m above sea level (a.s.l.). This glacier lies southeast of Ward Hunt Island near Canada’s northernmost point and has never been studied.
2. Henrietta-Nesmith Glacier (81.857°N, 73.059°W; Fig. 1c), which drains the central eastern region of the Northern Ellesmere Icefield, covering ∼1043 km2 and with elevations ranging from ∼100 to 2500 m a.s.l. Flowing from near the highest point on Ellesmere Island (Hattersley-Smith, Reference Hattersley-Smith1961), it contributes substantially to Lake Hazen’s hydrological input (Cavaco and others, Reference Cavaco2019).
3. John Evans Glacier (79.650°N, 74.187°W; Fig. 1d) is a predominantly cold polythermal glacier, which covers an area of ∼79 km2 and has a length of about 15 km, with elevations ranging from ∼100 to 1500 m a.s.l. Previous studies have documented a well-developed supraglacial drainage network below the ELA, which forms during a melt season lasting about 60 days (Bingham and others, Reference Bingham, Nienow, Sharp and Copland2006). The surface hydrological network is strongly linked to the subglacial drainage system (Bingham and others, Reference Bingham, Nienow, Sharp and Boon2005), with marked seasonal variability in glacier dynamics (Copland and others, Reference Copland, Sharp and Nienow2003b; Bingham and others, Reference Bingham, Nienow, Sharp and Copland2006).
4. Unnamed 2 Glacier (78.793°N, 76.538°W; Fig. 1f) is a large land-terminating glacier that drains the northern part of the Prince of Wales Icefield, with a basin size of ∼242 km2 and an elevation range from ∼30 to 1800 m a.s.l. Research conducted on nearby glaciers has reported rapid rates of mass loss since the mid-20th century (Bergsma and others, Reference Bergsma, Svoboda and Freedman1984), with rates nearly doubling between 1959–2001 and 2001–19, and the ELA rising by 77 m between 1974–82 and 2014–19 (Curley and others, Reference Curley, Kochtitzky, Edwards and Copland2021).
5. Sydkap Glacier (76.686°N, 85.461°W; Fig. 1e) is a large marine-terminating glacier with a basin area of ∼475 km2, draining the southern side of Sydkap Ice Cap. It lost 36.41 km2 of ice due to calving between 1959 and 2019 (Smith and others, Reference Smith2020), resulting in a terminus retreat exceeding 9 km, one of the greatest in the Canadian Arctic Archipelago (Copland and others, Reference Copland, Sharp and Dowdeswell2003a).
3. Data and methods
3.1. Time-series analysis
Various recent optical satellite images (Sentinel-2, PlanetScope) were compared to ASTER and SPOT images, and historical aerial photographs, to document temporal and spatial changes in supraglacial hydrology across Ellesmere Island since the late 1950s (Tables A1 and A2, both in Appendix A). Two different kinds of analyses were undertaken, based on image quality and availability:
a) Qualitative analysis: a series of images with approximately decadal separation was compiled for the five studied glaciers to qualitatively assess supraglacial drainage evolution, although image availability, particularly before 2000, did not always allow for this. This included Sentinel-2, ASTER and SPOT imagery, with resolutions of ≥10 m, which could only resolve channels ∼10 m wide or larger. While unsuitable for precise quantification, these images provided valuable insights into the evolution of major channels.
b) Quantitative analysis: detailed mapping was conducted for 1959 and 2020 for the four northernmost glaciers. In 2020, supraglacial drainage networks were mapped in detail by manually delineating channels on 3 m-resolution PlanetScope images following specific criteria outlined in Section 3.6 (Table A2). Mapping for 1959 used ∼0.6 m-resolution historical aerial photographs (Table A1). To maintain consistency, only channels ≥3 m wide (five pixels in the air photos) were included, aligning with PlanetScope resolution. While images were examined at various scales to interpret features, all delineation was performed at scale of 1:2500 to ensure consistency across glaciers and years, as well as to enhance the differentiation of tones and textures (Bash and others, Reference Bash2023). We note two important limitations in this mapping. First, full mapping was limited on Henrietta-Nesmith Glacier, due to poor aerial photograph image quality and snow cover in some images, particularly those taken up-glacier. As a result, mapping focused on larger channels in the lower ablation area of this glacier. Second, quantitative analysis was not carried out on Sydkap Glacier because it does not exhibit a well-defined perennial drainage system.
For all analyses, cloud-free images were preferentially acquired in July and August, when snow cover was minimal, ensuring supraglacial drainage networks were fully developed and mapped under consistent surface conditions (Tables A1 and A2).
3.1.1. Equilibrium line altitude
The 2020 ELA was estimated from the average elevation of the snowline in cloud-free, late summer (July–August) Sentinel-2 and PlanetScope satellite imagery spanning 2016–20, using the methodology of Kochtitzky and others (Reference Kochtitzky2023; Table A2). Although this method does not account for the presence of superimposed ice that may exist below the ELA, Casey and Kelly (Reference Casey and Kelly2010) found that the late-summer snowline offers the most reliable ELA estimate, showing strong agreement between remote sensing and in situ measurements. For each year, the snowline was manually delineated for all five glaciers by tracing the boundary between snow-covered and snow-free glacier surfaces in the late summer imagery. Except for Unnamed 1, which only had images available in 2017 and 2020, the ELA on all other glaciers was traced on four to five images. Elevation values along each polyline were extracted from 2 m-resolution ArcticDEM strips from 2020, and an average ELA was calculated from these values (Table A3 in Appendix A). A contour line at the average elevation was subsequently drawn across each glacier surface to represent the ∼2020 ELA. The standard deviation of the annual average snowline positions from 2016 to 2020 was used to represent the uncertainty in the ELA estimates, capturing interannual variability in late-summer snowline elevation.
The 1959 ELAs over the studied glaciers were estimated using regional reconstructions of the lowest ELAs in 1960, first conducted by Miller and others (Reference Miller, Bradley and Andrews1975) and later recreated by Wolken and others (Reference Wolken, England and Dyke2008). More recently, White and Copland (Reference White and Copland2018) digitized the recreated map of Wolken and others (Reference Wolken, England and Dyke2008) for northern Ellesmere Island and overlaid it on a Landsat 8 mosaic from 2015/16. This facilitated our estimate of the ELA for the smallest glacier, Unnamed 1, located near the northern coast of Ellesmere Island.
3.2. Historical aerial photographs
The earliest available imagery of the glaciers was 1959 panchromatic 1:60 000 scale vertical aerial photographs archived at the National Air Photo Library in Ottawa, Canada (Table A1). These images were acquired by the Royal Canadian Air Force and taken from an altitude of 9000 m a.s.l. using a Wild RC5A aerial film camera and Aviogon lens (15 AG). Images over sections of glaciers where channels were present were scanned at a resolution of 2033 dpi, resulting in a ground resolution of ∼0.6 m. Lower-resolution photos scanned at 300 dpi with a ground resolution of ∼4 m were used to provide coverage of the entire glacier surface for the quantitative analyses described later in the methods.
3.3. Satellite imagery
Satellite scenes used in both the qualitative and quantitative analyses were acquired from various sources (Table A2). SPOT 1-5 imagery covering the period from 1986 to 2015 was downloaded from the SPOT World Heritage data site (https://regards.cnes.fr/user/swh/modules/60). ASTER images, available since 1999 and acquired from the United States Geological Survey EarthExplorer database (https://earthexplorer.usgs.gov/), were used to fill gaps in the time series when SPOT images were not available. Sentinel-2 and PlanetScope imagery covered the most recent period since 2015, with Sentinel-2 images acquired from the Copernicus Open Access Hub (https://www.copernicus.eu/en) and the PlanetScope imagery downloaded through the Education and Research Program, courtesy of Planet Images (https://www.planet.com/). In addition, two WorldView-2 (WV2) and one WorldView-3 (WV3) image from 2022 covering the ablation area of Unnamed 2 and the entirety of Unnamed 1 glacier, respectively, were acquired through the European Space Agency’s Third-Party Missions program and used as validation for the mapping of sinks.
3.4. Digital elevation models
To assist with mapping in shadowed and low contrast regions in satellite imagery, a shaded relief raster (hillshade model) was used as a support. Specifically, strips from the 2 m ArcticDEM version 4.1 from the Polar Geospatial Center (https://www.pgc.umn.edu/data/arcticdem/) were used for this purpose (Table A3). The selected strips were closely aligned with the dates of the PlanetScope images for each glacier, ensuring consistency between datasets and reducing the risk of misclassifying channels with less well-defined banks (Bash and others, Reference Bash2023).
The hillshade model was generated in the Quantum Geographic Information System (QGIS) v 3.24 utilizing the Hillshade tool within the GDAL library. To enhance the visibility of depressions, a solar altitude angle of 45° was chosen, and the lighting direction was aligned to be perpendicular to the predominant flow of the glacier, given that channels typically align parallel to the glacier flow (Bash and others, Reference Bash2023). The hillshade model was produced without including shading from adjacent terrain, highlighting channel depressions with contrasting light and dark areas on either side (Bash and others, Reference Bash2023).
3.5. Georeferencing
Manual georeferencing was performed in QGIS for each historical aerial photo, utilizing the 2020 PlanetScope scenes of each glacier as the master imagery. Apart from Unnamed 1 Glacier, which was fully covered by a single photo, all other glaciers required two to six photos to cover their entire surface. In such cases, the aerial photo that captured the glacier’s terminus was initially georeferenced to its corresponding 2020 PlanetScope image. Subsequently, overlapping sections of the other aerial photos from that glacier were aligned to the previously georeferenced photo, while non-overlapping sections were directly aligned to the Planet imagery, ensuring improved alignment at the boundary of overlapping photos. Processing statistics of the final georeferenced products are provided in Table 1; the mean error reflects the average deviation between actual ground control point positions and their positions predicted by the transformation.
Table 1. Georeferencing statistics for historical aerial photographs
3.6. Channel mapping criteria
Although various automated and semi-automated methods have been developed to map supraglacial features, such as using spectral indices to detect water (Yang and Smith, Reference Yang and Smith2013; Lu and others, Reference Lu2020; Bash and others, Reference Bash2023) or DEM-based hydrological modeling (Rippin and others, Reference Rippin, Pomfret and King2015; King and others, Reference King, Hassan, Yang and Flowers2016; Yang and others, Reference Yang2019) to infer flow paths, their application to long-term studies remains challenging due to the scarcity of high-quality data and inconsistencies in the spatial and spectral resolution between contemporary and historical imagery. Consequently, we manually identified and delineated supraglacial channels in this study. While manual mapping introduces potential uncertainty due to user bias, subjectivity was minimized by applying a consistent set of morphological and textural criteria (described below) and by having all delineations performed by a single analyst.
Following the terminology of Bash and others (Reference Bash2023), channels were identified based on their curvilinear shape, smooth texture and distinct blue/white tones in multispectral PlanetScope imagery. A curvilinear shape, with smooth, sinuous outlines and few sharp angles, distinguishes natural water channels from features like crevasses. Channels also exhibited a homogeneous reflectance pattern, appearing smooth in texture in the multispectral imagery. While maintaining the same morphology, Lampkin and VanderBerg (Reference Lampkin and VanderBerg2014) also note that supraglacial channels show high contrast against the surrounding ice in panchromatic imagery (historical air photos), appearing darker. In this study, high mapping confidence was ensured by only including continuous channels, excluding features with repetitive breaks.
Once mapped, channels were classified according to terminology developed by St Germain and Moorman (Reference St Germain and Moorman2019), categorizing supraglacial channels into three main classes: surface, incised and canyons. Surface streams were identified by a valley depth-to-width ratio of less than one, meaning they lacked a distinct depression in the 2 m ArcticDEM but remained visible in PlanetScope imagery (Fig. 2 and Fig. A1 in Appendix A). Since no hillshade model was available for 1959, surface streams in historical aerial photographs were identified using image interpretation techniques (Bash and others, Reference Bash2023) and relied on the concept that channel width and depth are positively correlated (St Germain and Moorman, Reference St Germain and Moorman2019). The threshold beyond which depressions were no longer discernible in the 2020 hillshade model (∼6 m or two pixels in PlanetScope imagery) was used as a reference. Channels of this width in the 1959 aerial photographs, lacking clearly defined banks, were classified as surface streams. Well-defined channel banks suggest an incised feature where the valley depth-to-width ratio would likely be more than one, indicating that the channel is not of surface nature.
Figure 2. Decision tree for classifying supraglacial channels in imagery. Sections of the channel network are examined, and channels are classified as surface streams, incised rivers or canyons.
In contrast to surface streams, canyons are readily discernible in both satellite imagery and aerial photographs due to their substantial size and perennial presence. These deep, well-defined channels typically contain flowing water that occupies only a small portion of their valley volume (Fig. 2 and Fig. A1). Their prominent banks, often casting shadows, further distinguish them in satellite imagery.
All remaining channels, those deeper than surface streams but smaller than canyons, were classified as incised rivers. These channels exhibited strong contrast on either side due to well-defined banks and appeared as continuous features in the hillshade model. The complete mapping workflow is summarized in Fig. 2, and examples of each channel class, as they appear in the imagery and hillshade models, are provided in Fig. A1 (Appendix A).
3.6.1. Uncertainty quantification
To quantify uncertainty in channel length measurements, we assumed a positional error of ±1 pixel at both the start and endpoints of each digitized channel. When a channel flowed into another channel, its endpoint was excluded from the calculation to avoid double-counting. The cumulative uncertainty for each glacier was then estimated by multiplying the total number of considered start and endpoints by the pixel size. This approach was also applied to estimate the uncertainty in the total length of the individual channel classes.
3.7. Hydrography quantification
Several components of hydrography were derived from the final channel polylines, including total channel length, proportions of total length (PTL), drainage density (Dd) and sinuosity.
Dd standardizes channel length by glacier area, making it useful for comparing supraglacial networks across glaciers with varying surface areas and assessing channel abundance at different elevations. Here, Dd was calculated as the total channel length (km) divided by the glacier surface area (km2) for each glacier and for each stream class. To obtain glacier surface area for 1959 and 2020, we used the Randolph Glacier Inventory v7 outlines, which were manually corrected at glacier margins and drainage divides using the same ArcticDEM strips employed to create the hillshade models (Table A3). These strips were then clipped to match the extent of each glacier outline, and elevation polygons were generated at 100 m increments for both years to calculate PTL and Dd by elevation. PTL highlights channel type distributions independently of length changes, allowing for year-to-year comparisons.
Sinuosity, the ratio of channel length to straight-line distance, was also calculated for canyons on each glacier. Incised and surface channels were excluded due to their narrow width and/or transient nature, making consistent mapping difficult. The straight-line distance was drawn following the valley centerline (midpoint between channel walls), capturing directional changes along different reaches. The Mann–Whitney U test with a statistical significance level of p < 0.05 and boxplots were then used to assess whether sinuosity distributions statistically differed between years.
3.8. Moulins and sink points
Moulins typically require substantial and consistent water input, making them hydrologically significant when connected to active perennial supraglacial channels (Bash and others, Reference Bash2023). Similarly, crevasses can disrupt supraglacial drainage by capturing channel flow, storing water or transferring it to the other glacial hydrological subsystems (St Germain and Moorman, Reference St Germain and Moorman2019). Due to the difficulty of differentiating these features in satellite and aerial imagery, both hydrologically significant moulins and crevasses are referred to as sink points in this study. These were manually identified where channels terminated abruptly in the 1959 aerial photographs and 2020 PlanetScope imagery. The higher resolution of the historical aerial photographs often allowed direct identification of moulins as rounded, tone-distinct features, beyond which channels did not extend, which was not always possible in the PlanetScope imagery. To address this, manual validation was conducted on two glaciers using the high-resolution (0.5 m) panchromatic WV2 and WV3 images from summer 2022, covering the ablation area of Unnamed 2 Glacier and all of Unnamed 1 Glacier (Table A2).
3.9. Surface mass balance
To better understand the drivers behind the changes observed in various supraglacial systems, we utilized the daily modeled SMB dataset for the northern Canadian Arctic Archipelago from Noël and others (Reference Noël, Van De Berg, Lhermitte, Wouters, Schaffer and Van Den Broeke2018), updated to cover the period 1959 to 2020. SMB estimates are derived from the 11 km resolution regional climate model RACMO2.3, by summing individual daily components (precipitation, melt, runoff, refreezing, sublimation and drifting snow erosion). These outputs are then statistically downscaled to a finer 1 km resolution to better capture smaller ice masses. For additional information on the downscaling procedure, see Noël and others (Reference Noël, Van De Berg, Lhermitte, Wouters, Schaffer and Van Den Broeke2018).
In the current study, the 2020 ELA was used to divide the glaciers into their ablation and accumulation areas, and the specific SMB was calculated separately for each.
4. Results
4.1. Contemporary drainage patterns
The supraglacial drainage system of glaciers north of Sydkap generally follows a similar pattern, where channels are not clearly defined in areas of high topographic variability, such as within accumulation zones. Instead, they form a broad dendritic drainage network primarily composed of smaller streams and incised rivers with some terminating in sinks on both Unnamed 2 and John Evans (Fig. 3). Smaller streams and rivers that do progress down-glacier transition into sub-parallel configurations before converging into larger, entrenched incised rivers, which eventually develop into canyons flowing parallel to the glacier’s direction of flow, a pattern particularly evident on Henrietta-Nesmith Glacier (Fig. 3b). Smaller dendritic and sub-parallel tributary systems, primarily composed of streams and incised rivers, flow alongside these main channels, such as on Unnamed 1 Glacier where they cover most of the surface. This pattern is especially pronounced given that Unnamed 1 no longer possesses a distinct accumulation zone (Fig. 3a).
Figure 3. Canyon, incised and surface channels, along with sink areas, on (a) Unnamed 1 Glacier, (b) Henrietta-Nesmith Glacier, (c) John Evans Glacier and (d) Unnamed 2 Glacier, delineated on historical aerial photographs from 1959 and PlanetScope images from 2020. The solid black line represents the estimated ELAs for 1959 and 2020 (Note: no ELA is shown for Unnamed 1 in 2020 as it is now above the highest elevation of the glacier). The pie charts indicate the percentage contribution of each channel class to the total channel length in both years.
In the ablation zones of both Unnamed 2 and John Evans glaciers, where most of the supraglacial drainage system is concentrated (84.7% and 62.5%, respectively), meltwater transport is disrupted by prominent sink fields. These features, located at approximately 740 m elevation on Unnamed 2 and 600 m on John Evans, intercept flow from major supraglacial rivers, preventing continuous downstream transfer of meltwater from the upper to the lower ablation areas. In these lower sections, channels eventually reappear, fed by other streams and rivers originating from the southern portions of the glaciers’ accumulation areas (Fig. 3c and d). While a limited number of channels, particularly rivers, terminate in sinks near the mapped glaciers’ termini (only one sink on Unnamed 1), most flow off the glaciers’ surfaces (Fig. 3).
4.2. Qualitative multi-decadal changes
The qualitative time series analysis reveals substantial changes in the supraglacial hydrological networks of the studied glaciers over the past ∼60 years, with trends indicating a recent shift towards more established perennial drainage systems. These changes are observed in various aspects such as the number (e.g. Fig. 4 and Fig. B1 in Appendix B), length (e.g. Fig. 4), width (e.g. Fig. 4, Fig. B2 in Appendix B and Fig. 5), incision (e.g. Fig. 5) and sinuosity of the channels (e.g. Figs. B1 and B2), with some changes accelerating in the last decade. Additionally, some perennial channels have extended further up-glacier over the study period (e.g. Fig. 4 and Figs. B1 and B2).
Figure 4. Channel evolution on John Evans Glacier between 1959 and 2020. (a) Overview of the glacier, with the ablation area shown in panels (b–f) outlined by the black box. The red line represents the estimated 2020 ELA at 789 ± 52 m a.s.l. (b) ∼0.6 m resolution historical aerial photo, (c) 20 m resolution SPOT-1 image, (d) 15 m resolution ASTER image, (e) 5 m resolution SPOT-5 image, (f) 3 m resolution PlanetScope image. Colored dots show down-glacier progression of a channel over time and orange/pink cross-section lines in 1959 and 2020 mark the area where channel width was compared over time. See Tables A1 and A2 for image dates.
Figure 5. Channel evolution on Unnamed 1 Glacier between 1959 and 2020. (a) ∼4 m resolution historical aerial photo, (b) 20 m resolution SPOT-1 image, (c) 15 m resolution ASTER image, (d) 3 m resolution PlanetScope image. The red arrow is at the same location in each image and marks a canyon river that has seen particularly important changes in width and incision, with an increase in incision rate over the last decade. The orange cross-section line marks the area where channel width was compared over time. See Tables A1 and A2 for image dates.
John Evans Glacier exemplifies many of these changes. For example, a new perennial river formed in the northern part of the ablation area of the glacier between 1959 and 1987 (Fig. 4b and c). Channels have also lengthened over time, such as one channel that extended from the upper ablation area to lower elevations, connecting the upper supraglacial system with the system located in the terminus region (dots in Fig. 4c–f). This channel extended by ∼0.3 km between 1959 and 1987, barely extended between 1987 and 2002, extended by ∼0.5 km between 2002 and 2011, and eventually reached the downstream system sometime between 2011 and 2020. Channel banks and valley walls have also widened over time, especially along the now larger canyon rivers. For example, between 1959 and 2020, the channel cross-sections represented by the orange and pink lines widened by several tens of meters (∼99 and 34 m, respectively). In addition, larger rivers have increased in importance in the southern portion of the accumulation area, especially since the early 2010s (Fig. 4e and f).
On the northern Prince of Wales Icefield, Unnamed 2 Glacier has experienced significant changes in its supraglacial drainage network over time, particularly in channel sinuosity near the terminus. In this area, larger meanders have developed and subsequently cut through their necks, indicating dynamic channel evolution. As a result, channel banks have widened considerably (e.g. by 90 m since 1959 at the orange cross-section on the main channel near the terminus; Fig. B2b in Appendix B). Such changes are apparent in the main canyon at the terminus in the early 1990s but appear to have accelerated in the 21st century, especially with the expansion of additional canyons (Fig. B2). On Henrietta-Nesmith Glacier, channel dynamism is also observed within older, larger perennial canyons. These large rivers have also extended up-glacier, and numerous new perennial channels have emerged near and within the terminus region, especially since the start of the 21st century (see arrows in Fig. B1 in Appendix B). In far northern regions, such as Unnamed 1 Glacier, incision and channel bank widening have become increasingly pronounced across much of the glacier, particularly over the past decade, with canyon valley width approximately increasing by 77 m along the orange-marked cross-section since 1959 (Fig. 5).
In contrast to the other glaciers, the supraglacial hydrology of Sydkap Glacier has exhibited a different temporal evolution. A large crevasse field dominates the terminus region and extends ∼7.5 km up-glacier, inhibiting the formation of a perennial hydrological system (Fig. 6b). Further up-glacier, channels have remained poorly defined throughout the study period, indicating that they may change their course frequently (Fig. 6c–f). Only near the equilibrium line and within the accumulation area are larger channels found, but these disappear relatively quickly as they progress down-glacier, terminating in sinks, crevasses or smaller distributed systems (Fig. 6d–f).
Figure 6. (a) PlanetScope image of Sydkap Glacier, showing the absence of a perennial supraglacial drainage system; boxes indicate the location of images in (b–f), and the red line represents the estimated 2020 ELA at 831 ± 32 m a.s.l. (b) Close-up of the terminus, highlighting the upper limit of the crevasse field (green line), which inhibits the uninterrupted transport of meltwater across the glacier surface ∼7.5 km from the terminus. (c) Close-up showing an area further up-glacier, where supraglacial channels remain limited even in the absence of crevasses. (d–f) Channel evolution in the accumulation area of Sydkap Glacier from 1959 (∼4 m resolution historical aerial photo) to 2012 (10 m resolution SPOT-5 image) to 2020 (10 m resolution Sentinel-2 image). See Tables A1 and A2 for image dates.
4.3. Quantitative drainage system evolution between 1959 and 2020
Since 1959, the general structure of the supraglacial drainage systems has remained similar (i.e. dendritic networks at higher elevations feeding more streamlined systems at lower elevations), but each system has expanded, particularly into accumulation areas (Fig. 3). Over the past ∼60 years, total channel length and consequently total Dd, have increased across all glaciers. By 2020, Unnamed 2 and John Evans Glacier had total channel lengths of 474 ± 1.3 km and 283 ± 1.4 km, respectively, reflecting increases of 35.8% ± 0.4% and 39.4% ± 0.7% since 1959. In contrast, Unnamed 1 Glacier, the northernmost, saw a modest 7.1% ± 0.8% increase to 106 ± 0.8 km. Consequently, total Dd rose from 1.38 ± 0.001, 2.45 ± 0.002 and 5.21 ± 0.01 in 1959, to 1.96 ± 0.01, 3.57 ± 0.02 and 5.83 ± 0.04 km km−2 in 2020 on these three glaciers, respectively, with Dd remaining lowest on Unnamed 2 and highest on Unnamed 1 over the years (Fig. 7). Henrietta-Nesmith Glacier also experienced drainage system expansion, with total canyon length increasing by 204.8% ± 1% to 128 km by 2020.
Figure 7. Contribution of each channel class to the total drainage density (Dd; km km−2) for each glacier for the years 1959 and 2020.
The relative contribution of channel classes to the total length has also shifted across all glaciers towards increased channelization and incision (Fig. 7). In 1959, surface streams accounted for approximately 78.3%, 60.2% and 56.8% of the total channel length on Unnamed 2, John Evans and Unnamed 1 glaciers, respectively. By 2020, incised rivers had become the dominant channel type, comprising approximately 53.6%, 54.9% and 60.0% of the total length. Canyons, which only made up approximately 1.6%, 15.8% and 17.0% of the total length of channels on these glaciers in 1959, made up approximately 3.9%, 19.9% and 17.9% by 2020 (pie charts in Fig. 3).
While the overall trend of changing channel types is evident across all three glaciers, the absolute rate and relative shift toward a more perennial drainage system vary between them. On Unnamed 2 Glacier, canyon length tripled from 6 ± 0.03 km in 1959 to 18 ± 0.4 km in 2020, expanding at an average rate of approximately 0.20 ± 0.01 km a−1. However, canyons still contribute less to the total channel length compared to Unnamed 1 Glacier, despite the latter experiencing a slower growth rate of 0.03 ± 0.003 km a−1 (pie charts in Fig. 3). A similar pattern is observed for incised rivers on these two glaciers, contributing to the persistence of perennial features, defined here as the combined length of canyons and incised rivers, which remains highest on Unnamed 1 Glacier compared to all other glaciers. Although John Evans and Henrietta-Nesmith Glaciers exhibited higher absolute rates of canyon expansion (0.39 ± 0.003 and 1.41 ± 0.01 km a−1, respectively) than Unnamed 2, their relative increases were still lower in comparison.
4.3.1. Spatiotemporal variability of supraglacial channels
Between 1959 and 2020, total Dd increased across all elevations on both John Evans Glacier and Unnamed 2 Glacier, with the most pronounced rise occurring within the first 300 m bands from the terminus (Fig. 8c, d, g and h). In contrast, on Unnamed 1 Glacier, the most significant changes occurred above the 1959 ELA, beyond 400 m, with peak total Dd shifting to a higher elevation (Fig. 8a and e).
Figure 8. Drainage density (Dd; km km−2) by channel class as a function of elevation for the four mapped glaciers in 1959 and 2020. Total Dd combining all classes (black solid line) is shown on the secondary axis. Graphs are organized by latitude based on the location of each glacier: (a and e) Unnamed 1, (b and f) Henrietta-Nesmith, (c and g) John Evans and (d and h) Unnamed 2. Elevation is capped at 1400 m, which represents the highest elevation band at which channels were mapped. The black and red dotted lines show the position of the 1959 and 2020 ELAs, respectively, with the red shaded areas representing the uncertainty range for the 2020 ELAs. No ELA is shown for Unnamed 1 in 2020 as it is now above the highest elevation of the glacier at ∼980 m a.s.l.
These changes in total Dd are primarily driven by the development and expansion of perennial (incised) channels, as surface stream Dd has decreased across all glacier surfaces (Fig. 8). On Unnamed 2 Glacier, the PTL of canyons, which was previously concentrated within the 100 m band, is now more evenly distributed between the 100 and 200 m bands (Fig. 9d and h), with canyon Dd increasing across all bands up to 500 m, where they were previously absent (Fig. 8d and h). Incised rivers are also more prominent at higher elevations, with a particularly pronounced increase in their PTL observed between the 600 and 800 m bands by 2020 (Fig. 9d and h). Their Dd has increased substantially at and beyond the 800 m band, with rivers extending up-glacier past their 1959 maximum extent within the 1100 m band, to now reach the 1400 m band (Fig. 8d and h).
Figure 9. Proportion of total length (PTL; %) by channel class as a function of elevation for the four mapped glaciers in 1959 and 2020. Graphs are organized by latitude based on the location of each glacier: (a and e) Unnamed 1, (b and f) Henrietta-Nesmith, (c and g) John Evans and (d and h) Unnamed 2. Elevation is capped at 1400 m, which represents the highest elevation band at which channels were mapped. The black and red dotted lines show the position of the 1959 and 2020 ELAs, respectively, with the red shaded areas representing the uncertainty range for the 2020 ELAs. No ELA is shown for Unnamed 1 in 2020 as it is now above the highest elevation of the glacier at ∼980 m a.s.l.
A similar trend is evident on John Evans Glacier, where canyon distribution has expanded, and Dd has increased between the 100 and 700 m bands and into the 800 m band where they were previously absent (Fig. 8c and g). Incised rivers, which had a PTL formerly concentrated in the 600 m band, now peak within the 800 m band, with their maximum elevation extending from the 900 to 1100 m band by 2020 (Fig. 9c and g). The Dd of incised rivers has nearly doubled or more between the 400 and 900 m bands, with the most pronounced increases occurring in the upper ablation zone near 700 m (Fig. 8c and g).
Henrietta-Nesmith Glacier exhibits a comparable pattern, with canyon length extending from a maximum elevation within the 800 m band in 1959 to the 1100 m band in 2020 (Fig. 9b and f ). Concurrently, peak Dd has shifted from the 200 to 400 m band, reflecting the overall increase across all elevation bands (Fig. 8b and f).
As a result, canyons have expanded across a broader elevation range on these three glaciers compared to 1959, now dominating the lower ablation areas. Meanwhile, incised rivers have become more prevalent further up-glacier, with peak Dd values shifting to higher elevations.
In contrast to these glaciers, Unnamed 1 Glacier exhibits more limited changes in canyon distribution (PTL) (Fig. 9a and e), despite a rise in Dd at almost all elevations and a threefold rise in the 500 m band (Fig. 8a and e). Similarly, the Dd of incised rivers has more than doubled across much of the glacier, with the most notable increases occurring beyond the 1959 ELA, which has now retreated beyond the glacier’s highest elevation (Fig. 8a and e). This channel class now dominates most of the glacier’s surface, even in elevation bands where canyons are most abundant.
4.3.2. Changes in sinuosity
In 1959, sinuosity values were generally higher for most glaciers, except for John Evans, compared to 2020. This is particularly evident on Unnamed 1 Glacier, where the median, interquartile range and a single outlier from 1959 show substantially higher sinuosity values than those recorded in 2020 (Fig. 10). Despite the overall decrease, some glaciers exhibit considerable variability, with multiple outliers in both years indicating sporadic instances of higher sinuosity. This variability is especially notable for glaciers other than Unnamed 1, where the number of outliers increased, and maximum values rose for Unnamed 2 and John Evans (Fig. 10). Interquartile ranges for 1959 are larger for northern glaciers (Henrietta-Nesmith, Unnamed 1), while southern glaciers (Unnamed 2, John Evans) show increased variability in 2020, indicating a trend towards greater variability in the south and reduced variability in the north (Fig. 10). Despite these observations, p values do not yield statistically significant results when comparing the distributions of sinuosity values in 1959 and 2020.
Figure 10. Box plot showing the median, interquartile range and outlier sinuosity values for Unnamed 2, John Evans, Henrietta-Nesmith and Unnamed 1 glaciers for 1959 and 2020.
4.4. Changes in SMB
Modeled monthly specific SMB for the ablation areas of the five studied glaciers reveals significant trends in mass losses since the mid-1960s (p < 0.05) (Fig. 11a–e and Table C1 in Appendix C). These losses are primarily associated with the melt season, which typically occurs from June to August. Since the mid-1960s, snowfall (i.e. SMB gains) in the ablation areas of the five studied glaciers have shown generally low variability (Fig. 11a–e), with significant decreases observed only on Sydkap and Henrietta-Nesmith (Table C1).
Figure 11. Annual surface mass balance (SMB) components for the ablation (a–e) and accumulation (f–i) areas of the five studied glaciers from 1959 to 2020: Gains (blue), losses (pink) and net annual (purple), with trendlines from 1964 depicted in dark blue (gains), red (losses) and violet (net annual). Cumulative net specific SMB (SSMB; black solid line) is shown on the secondary axis. SMB values were derived from the regional climate model RACMO2.3 downscaled to 1 km resolution, and specific glacier values were computed by summing all raster values within the glacier area and dividing by the glacier’s surface area.
All glaciers show a significant negative trend in net annual specific SMB, with mass losses increasing between 48.14% and 127% between the periods 1964–99 and 2000–20 in the ablation area of all studied glaciers (Fig. 11a–e and Table 2). The magnitude and rate of these losses vary by latitude, with the highest values observed on Sydkap Glacier (furthest south) and the lowest on Unnamed 1 Glacier (furthest north) (Fig. 11a–e and Table C1). Cumulative summer SMB losses between 1964 and 2020 reached −51.45 m w.e. on Sydkap compared to −23.72 m w.e. on Unnamed 1 (Fig. 11a and e). Linear regression slopes also highlight this contrast, indicating faster loss rates on Sydkap (−13 mm w.e. a−1) than on Unnamed 1 (−6.24 mm w.e. a−1) (Table C1).
Table 2. Average losses (mm w.e. a−1) used as a proxy for surface melt across the ablation (Ab) and accumulation (Ac) areas of the five studied glaciers over the 1959–2020 period, derived from the regional climate model RACMO2.3
In the accumulation area, losses also exhibit a significant negative trend from 1964 to 2020 (p < 0.05), with mass gains also showing significant negative trends on some glaciers (Fig. 11f–i and Table C1). However, losses are increasing at a much faster rate than gains are decreasing (Fig. 11f–i and Table C1). Between 1964–99 and 2000–20, accumulation area losses more than doubled on Sydkap and John Evans, more than tripled on Unnamed 2, and increased over elevenfold on Henrietta-Nesmith (Table 2). These relative changes in losses are much greater in the accumulation areas compared to their respective ablation zones, highlighting the increasing vulnerability of higher-elevation regions (Table 2).
Following this pattern of increasing losses, the ELA has risen on all mapped glaciers since 1959, though the magnitude of change varies considerably. On northern Ellesmere Island, the ELA now lies above the highest point of Unnamed 1 Glacier, an increase of more than 580 m a.s.l., effectively eliminating its accumulation zone (Fig. 12a). On Henrietta-Nesmith Glacier, the ELA rose by approximately 53 ± 77 m between 1959 and 2020, yet the impact on its Accumulation Area Ratio (AAR) was minimal, which remains high at ∼0.81 (Fig. 12b) (Dyurgerov and others, Reference Dyurgerov, Meier and Bahr2009). Further south, John Evans and Unnamed 2 glaciers experienced ELA increases of 89 ± 52 and 145 ± 50 m, respectively, resulting in notable AAR reductions of ∼12.5% and ∼20.2% (Fig. 12c and d).
Figure 12. Glacier hypsometry for (a) Unnamed 1, (b) Henrietta-Nesmith, (c) John Evans and (d) Unnamed 2 glaciers. The ELA for each glacier in 1959 and 2020 is shown as the blue and red dashed lines, respectively, with the red shaded areas representing the uncertainty range for the 2020 ELAs. Only the 1959 ELA is included for Unnamed 1 as the highest point on this glacier is now below the regional ELA.
5. Discussion
5.1. Supraglacial drainage systems
Studies on supraglacial drainage patterns in the northern Canadian Arctic Archipelago remain limited. However, Lu and others (Reference Lu2020) reported patterns similar to those observed in this study on Devon Ice Cap, where high-elevation supraglacial channels formed broad dendritic networks that became increasingly channelized and subparallel at lower elevations, eventually developing into regularly spaced, parallel mainstems within outlet glacier trunks. Consistent with our findings and studies on polythermal glaciers in the northern and southern Canadian Arctic Archipelago, Greenland and Svalbard (Rippin and others, Reference Rippin, Pomfret and King2015; St Germain and Moorman, Reference St Germain and Moorman2016, Reference Cavaco2019; Yang and others, Reference Yang2019; Lu and others, Reference Lu2020), sinuosity was highest within the mainstems and increased towards the glacier terminus. Charlton (Reference Charlton2007) notes that natural rivers adjust their slope, typically decreasing it, by increasing sinuosity. On Fountain Glacier, Bylot Island, St Germain and Moorman (Reference St Germain and Moorman2019) found that supraglacial channels remained relatively straight until reaching the sloped terminus, where lateral meanders developed at rates of 0.2–1.9 m a−1, depending on channel size.
Although our study did not specifically examine meander growth, sinuosity values for canyon rivers show no significant pattern of change over the study period. This likely reflects the dynamic nature of these channels: continuous bank erosion from melting leads to meander cutoffs and channel straightening, after which new meanders begin to form as the channels readjust their slope (Rippin and others, Reference Rippin, Pomfret and King2015; St Germain and Moorman, Reference St Germain and Moorman2019). This process is clearly visible in the time series for Unnamed 2 through Henrietta-Nesmith glaciers (Figs. 4, 5, B1, B2). This is supported by St Germain and Moorman (Reference St Germain and Moorman2019), who show that canyon sinuosity is strongly correlated with slope and watershed area, which influence discharge and cause meander development rates to vary between canyons.
5.1.1. System evolution
While drainage patterns were broadly consistent across glaciers and over time, total channel length and Dd differed significantly, with the highest Dd on the northernmost glacier (Unnamed 1) and the lowest on the southernmost (Unnamed 2) in both 1959 and 2020 (Fig. 7). Glaciers with extensive cold ice are less dynamic and have smoother, less crevassed surfaces, limiting surface meltwater access to englacial and subglacial pathways (Rippin and others, Reference Rippin, Pomfret and King2015). As a result, supraglacial channels on such glaciers are often widespread, draining meltwater over large areas (Irvine-Fynn and others, Reference Irvine-Fynn, Hodson, Moorman, Vatne and Hubbard2011; Ryser and others, Reference Ryser, Lüthi, Blindow, Suckro, Funk and Bauder2013; Rippin and others, Reference Rippin, Pomfret and King2015). Bash and others (Reference Bash2023) found that channels were 98% more common on Thores Glacier on Northern Ellesmere Island, where ice is likely frozen to the bed (Kochtitzky and others, Reference Kochtitzky2023), than on Nàłùdäy (Lowell Glacier) in the Yukon, suggesting broader channel development on colder ice masses, consistent with this study. In addition, no moulins were found on Thores Glacier (Kochtitzky and others, Reference Kochtitzky2023), consistent with our observations that prominent sink fields were present on both Unnamed 2 and John Evans glaciers to the south, whereas only limited sinks were mapped on Henrietta-Nesmith and Unnamed 1 to the north (Fig. 3).
This pattern echoes findings by Yang and others (Reference Yang2019), who noted that moulins are rare in northwest Greenland but vital to supraglacial drainage in the southwest, attributing this difference to colder, thicker and less permeable ice in the north, along with slower flow and lower strain rates. Similar conditions may explain their absence on northern glaciers in our study area, as Kochtitzky and others (Reference Kochtitzky2023) reported nearby Thores Glacier to be slow-moving and cold-based, with modeled basal temperatures ranging from −7°C to −12°C. Previous research on John Evans Glacier has shown that sink fields are linked to seasonal increases in surface motion, highlighting its sensitivity to variations in surface meltwater input (Bingham and others, Reference Bingham, Nienow and Sharp2003). The emergence of new sinks in the western accumulation area of Unnamed 2 Glacier suggests the development of a subglacial drainage system, which could lead to enhanced basal flow (Fig. 3d). At a similar latitude on Axel Heiberg Island, Thomson and Copland (Thomson and Copland, Reference Thomson and Copland2017, Reference Thomson and Copland2018) observed that the upglacier development of moulins since the 1960s on White Glacier enabled basal sliding to occur at higher elevations than previously recorded. These findings suggest that southern glaciers exhibit stronger hydrological coupling, allowing more meltwater to reach the bed, whereas weaker coupling in northern glaciers retains water at the surface, promoting more widespread supraglacial drainage systems.
While changes in Dd on Unnamed 2 and John Evans glaciers were primarily driven by increases in total channel length, rather than reductions in glacier area, approximately 50% of the increase in Dd on Unnamed 1 can be attributed to surface area loss. This aligns with the relatively small increase in channel length on Unnamed 1 and the fact that the glacier now lies entirely within the ablation zone, indicating a negative mass balance across its entire surface (Fig. 3a).
One consistent pattern observed across all glaciers is the shift toward a more perennial drainage system, marked by an increasing contribution of incised rivers and canyons over time (pie charts in Fig. 3). This trend highlights an important transformation: while canyons and incised rivers have long played a significant role at lower elevations, they have become increasingly widespread across the glacier surfaces. This includes canyon development on Henrietta-Nesmith, which has seen an increase in both PTL at higher elevations and Dd over the study period (Figs. 8b and f and 9b and f).
On Unnamed 1 Glacier, the retreat of the ELA beyond the glacier’s highest point likely contributed to a concurrent rise in Dd and PTL of incised rivers at higher altitudes (Figs. 8a and e and 9a and e). In contrast, John Evans and Unnamed 2 glaciers showed the greatest changes in Dd at lower elevations, where increases were mainly driven by the rapid expansion of canyon rivers into the upper ablation zone (Fig. 8c, d, g and h). While Dd did not increase as significantly at higher elevations on John Evans and Unnamed 2, incised rivers still showed substantial development there (Fig. 8c, d, g and h).
Our observations suggest that the surface drainage system of Unnamed 1 Glacier adjusted to prevailing climatic conditions earlier than those of the other glaciers, likely due to its smaller watershed area and already high total Dd in 1959. This is further supported by the limited increase in total channel length over the study period, as well as the lowest absolute and relative rates of change in perennial rivers observed among all glaciers. These findings indicate that a well-established supraglacial drainage network was already in place during the earlier period, unlike Unnamed 2 Glacier, which, despite having the lowest total Dd, exhibits the greatest relative change in perennial rivers, suggesting rapid hydrological reorganization. As with Unnamed 1, the watershed area likely influenced the absolute canyon growth rates on both John Evans and Henrietta-Nesmith glaciers. Henrietta-Nesmith, the largest glacier and John Evans, with its wide ablation zone, have geometries that can potentially support the formation of large individual canyon watersheds, which may enhance meltwater discharge and contribute to canyon development (Fig. 3b and c). These findings suggest that while northern glaciers may exhibit greater stability and higher drainage densities due to broader environmental conditions, glacier-specific characteristics, such as geometry, may play a key role and should be considered in future studies.
Overall, the most prominent changes in incised rivers occurred in the transition zone between ablation and accumulation areas, typified by glaciers Unnamed 2 through Unnamed 1, reflecting the upward shift of the ELA since 1959. This transition has resulted in higher melt rates and increased meltwater availability, promoting enhanced erosion and the formation of new incised rivers throughout this zone, as well as within the initial elevation bands of the accumulation areas, particularly on John Evans Glacier (Fig. 8c and g).
To date, only St Germain and Moorman (Reference St Germain and Moorman2019) have examined the mechanisms driving the development of different supraglacial channel types in the Canadian Arctic Archipelago. Their findings, which align with earlier work by Knighton (Reference Knighton1981) in Norway and Marston (Reference Marston1983) in Alaska, indicate that higher meltwater discharge and steeper glacier gradients are strongly associated with faster incision rates. This suggests that the combination of favorable glacier gradients and high discharge over the study period likely facilitated substantial erosion of channel beds and banks, ultimately enabling the formation and expansion of incised rivers and canyons across our studied glaciers.
5.2. Hydrological implications of SMB changes
We observe small or insignificant trends in gains in the ablation zones of the studied glaciers since the mid-1960s, suggesting that the overall decline in net annual SMB is primarily driven by heightened summer melt. These intensified melt conditions are reflected in the evolution of the supraglacial hydrological systems, which have developed differently depending on each glacier’s melt regime.
In the far northern regions, lower surface melt rates and generally slower rates of change in losses since the mid-1960s have supported the long-term persistence of a well-developed supraglacial drainage system (Rippin and others, Reference Rippin, Pomfret and King2015) on Unnamed 1. In contrast, the more southerly glaciers, such as on Unnamed 2 and John Evans, have experienced stronger melt intensities that have accelerated supraglacial system development over the past ∼60 years (Fig. 11 and Table 2). On Sydkap Glacier, the combination of high surface ablation and extensive crevassing is likely responsible for the lack of a perennial system, as ablation exceeds incision, limiting channel development beyond the crevasse field.
Relative changes in losses have become particularly pronounced in the accumulation areas of the studied glaciers, where mass losses are increasing at a much faster rate than mass gains are declining (Fig. 11f–i and Table C1). This suggests that supraglacial system development at higher elevations is also primarily melt-driven. The accelerated relative rate of change in the accumulation areas, compared to their respective ablation zones, is consistent with the observed elevation-dependent warming across the glaciated regions of the northern Canadian Arctic Archipelago since the start of the 21st century. Specifically, the linear rate of change of clear-sky mean summer land surface temperature (°C a−1) has been 40% greater at elevations above 1000 m a.s.l. and 76.9% greater at elevations above 1400 m a.s.l. (Mortimer and others, Reference Mortimer, Sharp and Wouters2016).
The accelerating ice loss in the accumulation areas of the studied glaciers underscores their vulnerability and highlights the potential for expansion of supraglacial channel networks at higher elevations. Since 1959, such expansion has been observed through the development of incised rivers in regions that were previously part of the accumulation zone but have since retreated (Fig. 3). Even relatively minor changes in the ELA can shift large portions of glaciers from the accumulation zone to the ablation zone, facilitating the observed growth of the supraglacial channel networks. This demonstrates the high vulnerability of glaciers with large areas at elevations near the ELA and allowed for rapid expansion of the supraglacial drainage system at higher elevations. If the ELA continues to rise at the observed rates for these glaciers, Unnamed 2 is likely to be particularly vulnerable, due to the large amount of ice concentrated between 1000 and 1200 m a.s.l., potentially accelerating the development of perennial supraglacial features over a larger surface area (Fig. 12d). Unnamed 1 Glacier, with its surface already below the regional ELA, is expected to continue retreating and will ultimately disappear, following trends observed in smaller glaciers and ice caps across Ellesmere Island (Serreze and others, Reference Serreze, Raup, Braun, Hardy and Bradley2017; White and Copland, Reference White and Copland2018; Medrzycka and others, Reference Medrzycka, Copland and Noël2023).
5.3. Limitations and future outlook
While this study provides the first multi-decadal perspective on supraglacial drainage evolution in the Canadian Arctic Archipelago, several limitations constrain the interpretation of our results. The most significant is the absence of a DEM for 1959. As a result, channel mapping for that year relied exclusively on image interpretation. Although the high resolution of the historical aerial photographs (∼0.6 m) allowed for the distinction between surface and incised channels, thanks to the visibility of clear channel banks, mapping in 1959 was inherently more subjective than in 2020, where DEM-supported hillshade models enhanced our mapping confidence.
A second major limitation is the scarcity of high-resolution imagery between 1959 and 2020, which restricted our quantitative analysis to just those two years. This temporal gap limits our ability to determine what proportion of supraglacial hydrological development reflects sustained decadal trends versus transient surface hydrological conditions driven by year-to-year variations in melt intensity, particularly for the more ephemeral surface streams. Our decadal qualitative analysis helped partially mitigate this uncertainty by documenting general trends in channel expansion and persistence, but higher temporal resolution would enable a more nuanced understanding of when and how key hydrological shifts occurred.
Future research should aim to link hydrological changes to climatic drivers such as temperature gradients, precipitation regimes, glacier elevation and proximity to the coast. Incorporating regional climate datasets and modeling efforts will help refine predictions of how these systems will continue to evolve in a warming Arctic. Moreover, the significant changes in supraglacial systems presented in this study should be examined within the broader context of glacier dynamics to better understand how glaciers will respond to ongoing and future climatic changes.
6. Summary and conclusion
In this study, we provided the first detailed, multi-decadal analysis of supraglacial drainage development on five glaciers across Ellesmere Island, offering insights into long-term hydrological changes in the northern Canadian Arctic Archipelago. Between 1959 and 2020, Dd increased on glaciers north of Sydkap, primarily due to the expansion of incised and canyon rivers at higher elevations as ELAs retreated up-glacier. These changes reflect a broader reorganization of supraglacial drainage networks, with perennial, incised channels becoming dominant. Nonetheless, variations in the rate of change of supraglacial hydrological systems exist with northern glaciers such as Unnamed 1 displaying a more stable and less dynamic system compared to glaciers further south such as John Evans and Unnamed 2.
Modeled increases in surface melt support these observations, with glaciers at lower latitudes experiencing faster rates of supraglacial system evolution. Despite this, high surface melt rates on Sydkap Glacier combined with heavy crevassing have inhibited the formation of a persistent surface drainage system there. Although Unnamed 1 has maintained a relatively stable drainage system due to lower melt rates and the likely persistence of cold surface ice, the recent rise of its ELA above the glacier’s highest elevation signals ongoing retreat and long-term vulnerability. As ELAs continue to rise, glaciers with most of their ice at mid elevations, such as Unnamed 2, are particularly susceptible to rapid hydrological expansion at higher elevations and continued mass loss.
Acknowledgements
This research was generously supported by funding from the Natural Sciences and Engineering Research Council (NSERC) of Canada, ArcticNet Network of Centres of Excellence Canada, Northern Scientific Training Program (NSTP), Polar Continental Shelf Program (PCSP), Amundsen Science, Weston Family Foundation, Fonds de Recherche du Québec (FRQ), University of Manitoba and University of Ottawa. We also extend our thanks to the community members of Grise Fiord for their support throughout this research, and to the Nunavut Research Institute for permission to conduct fieldwork. Brice Noël is a Research Associate of the Fonds de la Recherche Scientifique de Belgique—F.R.S.-FNRS.
Appendix A: Datasets and example visualization of supraglacial channel mapping
Table A1. Historical aerial photographs used for the qualitative (multidecadal time series) and quantitative (channel delineation) assessments of the supraglacial drainage systems on all five studied glaciers in 1959
Table A2. Satellite imagery used for the qualitative and quantitative time series analyses, including supraglacial channel delineation, sink identification and delineation of the 2020 ELAs
Table A3. ArcticDEM strips used to generate hillshade models and determine the 2020 ELAs for each examined glacier
Figure A1. Surface, incised and canyon channels as they appear in the 1959 historical aerial photographs and 2020 PlanetScope imagery, along with a comparison of the same channels in the hillshade model (inset).
Appendix B: Qualitative time series analysis of glaciers Henrietta-Nesmith and Unnamed 2
Figure B1. Channel evolution on Henrietta–Nesmith Glacier between 1959 and 2020. (a) Overview of the glacier, with the lower ablation area shown in parts (b–f) outlined by the black box. The red line represents the estimated 2020 ELA at 1253 ± 77 m a.s.l. (b) ∼0.6 m resolution historical aerial photo, (c) 20 m resolution SPOT-1 image, (d) and (e) 15 m resolution ASTER image, (f) 3 m resolution PlanetScope image. See Tables A1 and A2 for image dates.
Figure B2. Channel evolution on Unnamed 2 Glacier between 1959 and 2020. (a) Overview of the glacier, with the terminus area shown in parts (b–f) outlined by the black box. The red line represents the estimated 2020 ELA at 1045 ± 50 m a.s.l. The orange cross-section line marks the area where channel width was compared over time. (b) ∼0.6 m resolution historical aerial photo, (c) 20 m resolution SPOT-2 image, (d) 15 m resolution ASTER image, (e) 10 m resolution SPOT-5 image, (f) 3 m resolution PlanetScope image. See Tables A1 and A2 for image dates.
Appendix C: Trends and significance values for surface mass balance components
Table C1. Rate of change (slope, θ, mm a-1) and significance values (p) for gains, losses and net annual specific SMB over the period 1964–2020 across the (a) ablation and (b) accumulation areas of the five studied glaciers, derived from the regional climate model RACMO2.3
Note: Unnamed 1 was excluded from the accumulation area analysis as its surface area is now completely within the ablation zone.
Open access
