# Evidence for subglacial flooding in labyrinthine channels on Devon Island, Nunavut, Canada

Simona F Ruso, Gordon R Osinski, Anna Grau Galofre, Antero Kukko

► To cite this version:

Simona F Ruso, Gordon R Osinski, Anna Grau Galofre, Antero Kukko. Evidence for subglacial flooding in labyrinthine channels on Devon Island, Nunavut, Canada. 2026. hal-05455670

HAL Id: hal-05455670

https://hal.science/hal-05455670v1

Preprint submitted on 13 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Evidence for subglacial flooding in labyrinthine channels on Devon Island, Nunavut, Canada

Simona F. Ruso<fnref n="1" />,<fnref n="2" />, Gordon R. Osinski<fnref n="1" />,<fnref n="2" />, Anna Grau Galofre<fnref n="3" />, Antero Kukko<fnref n="4" />

*Correspondence to:* Simona F. Ruso (sruso@uwo.ca)

**Abstract.** Subglacial drainage systems route glacial meltwater to the ice margin either via efficient, channelized systems or inefficient, distributed systems. The interplay between channelized and distributed drainage systems varies spatially and temporally, governed by meltwater supply and abundance, bed roughness and topography, ice sliding velocity, and ice driving stress. Subglacial channel formation and evolution are therefore affected by variability in meltwater supply to subglacial conduits, and these changes may be recorded in the geomorphology of these channels. The formation of subglacial bedrock channels is attributed to higher energy and/or higher magnitude discharge events, such as the episodic release of meltwater in the form of either subglacial or proglacial floods, in comparison to the energy or discharge required to excavate channels in soft sediment. Common features of landscapes modified by meltwater floods include anastomosing channels and multiple erosive surfaces, wherein the pre-existing drainage system is inundated, resulting in the incision of new channels that reconnect downstream. Devon Island in the Canadian Arctic Archipelago was covered by the thin (<1000 m), cold-to-polythermal based Innuitian Ice Sheet over the course of at least three glacial expansions during the last glacial cycle. Despite this, there is a conspicuous lack of typical glacial landforms, and instead, the inland plateau region of the island is incised by ubiquitous subglacial and lateral meltwater channels. Some sets of bedrock subglacial channels on Devon Island bear a striking resemblance to the morphology of The Labyrinth in Antarctica, which formed by the episodic drainage of a subglacial lake. The characteristics, topology, and morphology of these channels, referred to as 'Labyrinthine channels' hereafter, together with two subglacial channel networks make the focus of this study. We argue that, within both labyrinthine and other subglacial channel networks on Devon Island, the presence of distinct erosional surfaces, anastomosing channels, and profile slope breaks imply formation by short-lived locally intense episodes of erosion. The presence of well-defined erosional surfaces suggests floods progressively incised into lower elevations where meltwater was captured by pre-existing or incipient channels. Moreover, steep contacts between erosional surfaces, termed here as "slope breaks," are similar to fluvial knickpoints and hanging valleys found in other notable landscapes caused by flooding, such as the Channeled Scablands, possibly indicating

<footnote n="1">¹Department of Earth Sciences, University of Western Ontario, London, N6A 3K7, Canada</footnote>

<footnote n="2">²Institute for Earth and Space Exploration, London, N6A 3K7, Canada</footnote>

<footnote n="3">³Laboratoire de Planetologie et Geosciences, Nantes, 44322, France</footnote>

<footnote n="4">⁴Finnish Geospatial Research Institute, National Land Survey of Finland, 02150 Espoo, Finland</footnote>
<figure><img src="image_3.png" /></figure>

channel headward erosion in response to pulses of intense erosion. Overall, we suggest that the presence of discrete erosional surfaces implies multiple flooding events, and that changing flow conditions during these events are evidenced by slope breaks. Multiple erosional surfaces, scabland-type landscapes, anastomosing bedrock channels, and hanging valleys with steep slope breaks are not consistent with ice marginal melt, demanding large discharge conditions and pulses of activity, and pointing at subglacial rather than marginal or proglacial environments of formation. This work aids in enhancing the current understanding of the role and dynamics of meltwater drainage systems operating under the cold-to-polythermal based Innuitian Ice Sheet, perhaps shedding light into its retreat dynamics, and bolstering the interpretation of glacial dynamics on Devon Island.

# 1 Introduction

Subglacial meltwater channels are integral components of glacial hydrological systems and can provide insight into ice sheet dynamics (e.g., Dewald et al., 2022; Hogan et al., 2022; Jansen et al., 2014; Kleman & Borgström, 1996; Simkins et al., 2023). Subglacial channels form drainage systems wherein meltwater, whether from a supra-, en-, or subglacial source, is efficiently routed to the ice margin in pressurized conduits. Low water pressures in conduits concentrate meltwater along these pathways and effectively increase effective pressure (ice overburden minus water pressure) at the ice base, decelerating glacial sliding (Flowers, 2015; Lelandais et al., 2018; Schoof, 2010). A channelized drainage system differs from less efficient “distributed” systems, which may be composed of linked cavities, thin water sheets, in-sediment macroporous pathways, or soft-sediment canals, all acting to increase water pressures at the ice base and promote sliding (Flowers, 2015; Kamb, 1970).

The interplay between these two drainage systems is governed by meltwater supply rate and abundance, bed roughness and topography, ice sliding velocity, and ice driving stress, and, therefore, varies at both temporal and spatial scales (Cuffey and Paterson, 2010; Flowers, 2015; Grau Galofre et al., 2022; Lewington et al., 2020; Schoof, 2010). Variability in meltwater supply to subglacial conduits affects channel inception and evolution, and these changes may be recorded in the geomorphology of these channels. While tracer experiments have been successful in improving the current understanding of varying subglacial discharge and its effect on ice velocity over diurnal and seasonal timescales (e.g., Andrews et al., 2015; Chandler et al., 2013; Nienow et al., 1996), the formation, evolution, and topology of subglacial conduits themselves is poorly constrained. Direct observations of subglacial drainage under ice sheets are difficult to conduct due to accessibility constraints. Some researchers argue for time-transgressive formation of tunnel valleys, km-scale subglacial conduits, by incision of diurnal or seasonal meltwater (Kehew et al., 2012; Ó Cofaigh, 1996; Van der Vegt et al., 2012), formation via outburst flooding events (Jørgensen and Sandersen, 2006; Lewis et al., 2006), or mega-floods (Shaw, 2002).

Of particular interest because of their long-time preservation potential are subglacial meltwater drainage pathways carved into bedrock, termed N-channels at the mm- to m-scale (Nye, 1976), or tunnel valleys at the km-scale (Kehew et al., 2012). Considering the relative hardness of bedrock compared to soft sediment, the formation of subglacial bedrock channels is attributed to higher energy and/or higher magnitude discharge events such as the episodic release of meltwater in the form of
<figure><img src="image_4.png" /></figure>

either subglacial (e.g., Gombiner & Lesemann, 2024; Kirkham et al., 2019; Lewis et al., 2006; Sugden et al., 1991) or proglacial megafloods (Baker, 2009; O'Connor et al., 2020). A common feature of landscapes modified by meltwater floods are anastomosing channels, wherein the pre-existing drainage system is inundated, resulting in the incision of new channels that reconnect downstream. In rivers, anastomosis forms as older channel segments are bypassed or by flooding which results in overbank flow and incision of multiple new channels at once (Makaske, 2001; Schumm et al., 1996). While the former exists in long-lived channels across an entire floodplain, the latter is typically associated with discrete events or stages of river evolution (Makaske, 2001). In subglacial environments, water flow through conduits is driven by effective pressure gradients (Cuffey and Paterson, 2010; Flowers, 2015). As such, meltwater pressure in the conduit must exceed that of the overburden to create space for overbank flow and incision of anastomosing channels (e.g., Gombiner & Lesemann, 2024). It is for this reason that anastomosing channels in glacial environments invoke large outburst flood hypotheses wherein glaciers can provide ample meltwater to drive high-energy, large-magnitude floods (Alley et al., 2019). However, recent modelling predicts efficient subglacial bedrock erosion can occur over multiple moderate-magnitude events such as seasonal meltwater discharge or the episodic drainage of supra-, en-, or subglacial meltwater reservoirs (Beaud et al., 2018).

This study focuses on subglacial bedrock channels on central Devon Island, Nunavut, Canada (Fig. 1A), a site occupied by Innuitian cold- to polythermal-based glaciation with several advances during the last ~100 14C ka BP (Dalton et al., 2022), and with a last recorded episode retreating around ~8 14C ka BP from the central plateaus around our study region (England et al., 2006). Specifically, we focus on investigating the origin of anastomosing subglacial meltwater channels incised into predominantly dolomite bedrock on the inland plateaus of Devon Island (Fig. 1B). Some sets of channels bear a striking resemblance to the so-called ‘Labyrinth’ in Antarctica, a set of channels incised subglacially in carbonate bedrock during the high intensity, episodic drainage of a subglacial lake (Lewis et al., 2006; Sugden et al., 1991), and thus we refer to these Devon Island channel networks as “labyrinthine channels.” We show that the dominant anastomosing channel patterns, trapezoidal cross-sections, and undulous longitudinal profiles of The Labyrinth (Sugden et al., 1991) are comparable – albeit on a smaller scale – to our observations of the labyrinthine channels on Devon Island. We also expand on previous investigations of bedrock subglacial channels (Grau Galofre et al., 2018, Ruso et al., 2024) and document other instances of subglacial channel anastomosis on Devon Island (Fig. 1B). We propose that the formation of subglacial channels on Devon Island occurred via multiple episodes of large meltwater discharge events, perhaps characterizing the retreat dynamics of a largely cold-based Innuitian ice sheet and compare our study sites to other examples of landscapes modified by glacial outburst floods to refine our understanding of the nature of glacial and meltwater dynamics of the former Innuitian ice sheet on Devon Island.
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 1. Satellite imagery of study area (Google Earth; MAXAR and Earthstar Geographics). A) Devon Island with surrounding islands and interisland channels. B) Context image showing the study locations in north-central Devon Island. Note the Haughton impact structure and major fjords, valleys, and outlets to Thomas Lee Inlet (top right corner). C) Deglacial isochrons illustrating retreat of the Innuitian Ice Sheet on Devon Island from 11 ka 14C BP to present day (redrawn from Dyke, 1999).</figcaption></figure>

## 2 Study area

Devon Island is a 55,247 km² uninhabited island located between ~ 77 – 74°N and ~ 97 – 79°W in the Canadian Arctic Archipelago that hosts a remarkably well-preserved record of glacial and subglacial landforms, including numerous fjords along the coast and a combination of ice marginal and subglacial drainage pathways on the inland plateaus (Dyke, 1999; Grau Galofre et al., 2018; Lee et al., 1999; Ruso et al., 2024). During the course of at least three glacial expansions, Devon Island was repeatedly covered by the Innuitian Ice Sheet (IIS), although the extent and timing of these glaciations remain largely unconstrained (Dalton et al., 2022). Indeed, the glaciation of North America leading up to the Last Glacial Maximum (LGM) was dynamic; both the Laurentide and the Innuitian Ice Sheets likely experienced multiple cycles of advance and retreat from
<figure><img src="image_7.png" /></figure>

115 – 25 ka (Dalton et al., 2022). However, the Innuitian Ice Sheet (IIS), which covered the Queen Elizabeth Islands (QEI) of northern Nunavut, is poorly understood in terms of advance and retreat cycles. Modelling suggests that the IIS experienced at least two episodes of advance: significant growth of the IIS occurred from 70 – 50 ka, followed by decay around 45 ka, then reaching maximum extent at ~ 25 – 19 ka (Dalton et al., 2022; England et al., 2006). Minimum IIS extent estimates glaciated highland areas on Ellesmere, Axel Heiberg, and Devon Island, comparable to modern-day ice cover, while the maximum extent of the IIS is proposed to have covered all the QEI and interisland channels, likely reaching the continental shelf (Dalton et al., 2022). Deglaciation of the QEI began at 11 14C ka BP and culminated to approximately its present-day configuration by around 7.5 14C ka BP (England et al., 2006). On Devon Island, it was proposed that ice retreated inland towards the Devon Ice Cap, the final remnants of the IIS (Fig. 1C) (Dalton et al., 2022; Dyke, 1999). While the IIS has been recognised in glaciochronological (e.g., Dalton et al., 2022; England et al., 2006), geomorphological (e.g., Dyke et al., 2002; England et al., 2006; Grau Galofre et al., 2018; Ruso et al., 2024), and glacial isostatic modeling (e.g., Dyke, 1999; Simon et al., 2015), the extent, dynamics, meltwater content, and derived geological record of IIS glaciation on Devon Island and other islands in the QEI has yet to be fully explored.

105

110

115

The record of glaciation on Devon Island is sparse; glacial erosion features are rare and have only been documented along coastal areas (Dyke, 1999), while subglacial and lateral meltwater channels have been documented and explored in the inland plateau (Fig. 1B) (Dyke, 1999; Grau Galofre et al., 2018; Ruso et al., 2024). Ice marginal and subglacial drainage pathways on Devon Island have both been identified in the island plateaus (Dyke, our work). Marginal and subglacial drainage pathways on Devon Island are differentiated from each other based primarily on relative channel size and the degree of channel incision in bedrock; identified subglacial channels are incised in the sub-horizontal Ordovician Allen Bay dolomite formation, are heavily anastomosing, present different erosive surfaces, and have stripped the channel headwaters of soil and fine sediment (scabland), as discussed in following sections. Lateral meltwater channels are observed perched on valley walls and typically following elevation contours or draining into a main channel flowing oblique to topography. They incise into cryofractured dolomite to a much lesser degree, and do not present bare bedrock sections or different erosional levels, and they have a relatively smaller channel size; marginal channels are significantly shorter (~300 – 730 m), shallower (~1 – 4 m), and narrower (~20 – 40 m) (Ruso et al., 2024) than subglacial channels discussed in this study. Moreover, consistent spacing is a diagnostic feature of subglacial meltwater channels (Hewitt, 2011; Schoof, 2010) and is observed in subglacial channels throughout Devon Island in contrast to the spacing of marginal meltwater channels which is not as consistent. We discuss this aspect in more detail in Section 4.

120

125

130
<figure><img src="image_8.png" /></figure>

# 3 Methods

This work presents results following field campaigns to Devon Island in the summers of 2022 and 2024. Prior to fieldwork, sites of interest were identified and characterized using satellite imagery from Maxar (30 cm/pixel, available in open format from Google Earth) and stereo-derived topography data from the ArcticDEM (2m/pixel), as well as existing field image and topography data from previous campaigns (Grau Galofre et al., 2018, Ruso et al., 2024, supplementary information).

## 3.1 Field data acquisition

Fieldwork was conducted on foot, all-terrain vehicle, and by helicopter in July 2024 and included the study areas SG2, SG3, and the labyrinthine channels (Fig. 1B). Data collection for SG1 took place in 2016, 2017, and 2022 and is summarized in section 4.1 and previous publications (Grau Galofre et al., 2018 and Ruso et al., 2024). Channels in the western section of SG2 were traversed from the headwaters to the downstream reaches of a tributary channel and a portion of the main stem was explored in a 2 km long circuit (Fig 2D). SG3 was briefly traversed, however, significant snow cover and poor weather conditions inhibited detailed exploration of the area. In SG3, one tributary channel was traversed starting from ~ 300 m from the headwaters for a distance of ~ 700 m downstream (Fig. 2F). A 650 m portion of the labyrinth terrain as traversed on foot and by ATV over two days (Fig. 8B).

Photos were taken within channels and interfluves and qualitative observations such as lithology, grain size, scale, and morphology of channels were noted. Upstream or downstream changes in morphology and sedimentology were also documented.

We acquired image data inside and around the channels, and aerially from an uncrewed aerial vehicle (UAV). Geotagged images include observations of channel cross-sections, tributary divergent sections and junctions, nature of the bed (covered by frost-shattered dolomite or bedrock exposed, scabland), paleo-bedload buried under the post-glacial infill veneer, evolution of channel depth and width, comparison with lateral meltwater channels in the vicinity, etc. Lithology and sediment type, size, and characteristics were also noted, together with downchannel changes in any of these parameters.

A DJI Matrice 300 RTK Pilot Uncrewed Aerial Vehicle (UAV) equipped with a DJI Zenmuse L1 LiDAR camera (maximum 3 cm resolution, 10 cm horizontal accuracy, 5 cm vertical accuracy) was utilized in the field to obtain LiDAR scans and imagery of the labyrinth terrain and SG3. LiDAR-derived Digital Elevation Models (DEMs) of SG1 and SG2 were obtained in a previous fieldwork campaign (see Ruso et al., 2024). LiDAR scans were set to an altitude of 100 m. Scans of localized areas within the labyrinth terrain and SG2 were completed using a Akhka-R4 backpack kinematic LiDAR equipped with a Riegl MiniVUX-3UAV and VUX-1HA laser scanners. Data presented here was obtained from the Riegl MiniVUX-3UAV with a ranging accuracy of 15 mm with 10 mm precision and a maximum range of about 290 – 170 m. The surveys were carried out by walking along the channel floors and topographic highs to ensure maximum coverage of the area.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

## 3.2 Data processing

LiDAR data from the UAV was processed using DJI Terra to develop high-resolution DEMs. Data was first processed into point cloud files (.las) in DJI Terra before processing in ArcGIS Pro, where point cloud data was converted into a raster form. The workflow consists of the “Make LAS Dataset Layer” tool, and the “LAS to Raster” tool. LiDAR data from the kinematic LiDAR scans was output as compressed LAS files (.laz) which were then converted to LAS and then raster datasets (workflow: “Convert LAS” and “LAS to Raster” tools in ArcGIS Pro). Final DEM resolution for the labyrinth and anastomosing section of SG2 is 0.25 m. Hillshaded DEMs were produced from the raster and overlain by the original with colourized elevation to assist in visualization.

165

Channels were manually delineated in ArcGIS Pro and longitudinal profiles were extracted from channel polylines and plotted using the “Profile Graph” tool. Contours (1 m) were processed and overlain on the DEMs to identify clear slope breaks and areas of relatively consistent elevation; hereafter referred to as “terraces”. Erosional Surfaces (ES) were identified using 1 m and 5 m contour data overlain on the DEMs; contacts for ES were drawn at clear slope breaks. Where no clear slope breaks were identified, ES were delineated visually as areas of similar elevation ranges and similar morphological characteristics. ES were delineated manually on the DEMs as polygon shapefiles and notable slope breaks features were marked in the DEM as a point shapefile.

170

175
<figure><img src="image_11.png" /></figure>

# 4 Results

## 4.1 Subglacial Channel Network 1

Subglacial channel network 1 (SG1) is located roughly 16.5 km S/SW of Thomas Lee Inlet (75°20′16.23″ N, 89°23′50.76″ W) and ~ 2.5 km to the southeast of the Haughton impact structure (Figs. 1A-B, 2A-B). This network is incised into the dolomite bedrock of the Lower Allen Bay Formation within the Arctic Platform sedimentary succession (Osinski et al., 2005; Thorsteinsson and Mayr, 1987). SG1 covers an area of ~ 2.3 km² (Fig. 2A, B). SG1 is composed of a network of tributary channels oriented S-N that drain into a main stem channel running approximately SW-NE (Figs. 2B, C). The main stem of SG1 drains into large valleys that follow the crater rim northwards (Fig. 2A), eventually connecting to the same fjord as the labyrinth terrain (see section 4.4). The elevation of SG1 is 315 – 235 m ASL and regional slope is gently NNW.
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /><figcaption>Figure 2. A) Satellite image (Google Earth/Maxar/Landsat) of region surrounding SG1 (B-C), SG2 (D-E), and SG3 (F) on Devon Island. Transect from 2024 annotated as black dashed line in B). Channels are delineated in white. Magenta circles in D) and F) annotate locations of photos in Figures 4-5 and Figure 7, respectively. Resolution of DEMs in C) and E) are 0.41m/pixel and 0.28 m/pixel, respectively. Lines on DEMs are gaps in LiDAR data.</figcaption></figure>

We investigated the headwaters of SG1 along a transect that spanned the entire width of the network (Fig. 2B). Details are provided in Ruso et al. (2024). In brief, SG1 is comprised of a series of curved tributary channels that originate from the eastern end of the network and relatively shallow anastomosed pathways in the western end (Figs. 2B, C). Tributary channels grade smoothly into the bedrock headwaters, and they originate as meter wide and flat-bottomed scabland channels (bedrock is exposed on their floors), grading into anastomosing channels that narrow and deepen downstream before merging into the main stem channel (Fig. 2C), which connects to a larger valley and ultimately drains into a fjord (Fig. 2A). No post-glacial infill resulting from fluvial transport is observed in the tributary channels of SG1; channels contain angular gravels, cobbles, and boulders of the Allen Bay Formation dolomite, formed as a result of in-situ cryofracturing of the channel floors and walls, resulting in a talus covering parts of the channel floor.

Following our 2024 fieldwork, three erosional surfaces (ES) were identified in SG1 using drone imagery from 2022 and based on field observations from other sites visited in 2024 (SG2, SG3, labyrinthine channels); erosional surfaces were not described in previous work by Ruso et al. (2024). ES1 is the uppermost ES and ranges in elevation from 315 m – 260 m in our site (Fig. 3A) corresponding to minimum tributary channel floor elevation and is comprised of channels up to 1350 m in length that incise into local topographic highs. Channel dimensions are up to 120 m wide and 16 m deep. Some channels anastomose at the headwaters and downstream, and interfluves are elongated streamlined features, carved into bedrock (Figs. 2B, C). Both the elevation gradient and channel orientation are towards the main stem channel of the network, approximately S-N (Figs. 2C, 3B). Towards the northern edge of ES1 there are a set of short, shallow channels that are ~3 m deep, <10 m wide, maximum 200 m long, and exhibit minor anastomosis (Fig. 3A). The contact between ES1 and ES2 is gradational with no major slope breaks observed. Occupying the western end of SG1, anastomosing tributary channels oriented S-N comprise the second erosive surface (ES2) (Fig. 3A), which has an elevation range of 281 m – 256 m (minimum elevation corresponds to minimum channel floor elevation). Channels in ES2 are up to 112 m wide and 11 m deep. Channel anastomosis is common and occurs from the headwaters to the downstream reaches before meeting the main stem channel (Fig. 3A). Interfluves are always streamlined, yet they appear shorter and wider than those in ES1. At 210 m elevation, there is a slope break where ES2 meets ES3 (Figs. 3A, B). Erosive surface 3 (ES3) comprises the lowest surface with the main stem channel of the network (Fig. 3A), which has maximum dimensions of 172 m wide and 36 m deep.
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure 3. A) Erosional surfaces (ES) of SG1 overlain on LiDAR-based DEM (0.41 m/pixel resolution). B) ES overlain by 5 m contour intervals to illustrate relatively flat “terrace” areas, slope breaks, and channel incision into topographic highs.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

## 4.2 Subglacial channel network 2

Subglacial channel network 2 is located at 75°17′21.01″ N, 89°11′16.80″ W ~17 km south of Thomas Lee Inlet and ~10.5 km southeast of the Haughton impact structure (Figs. 1A-B, 2A), also incised into the dolomite bedrock of the Lower Allen Bay Formation (Grau Galofre et al., 2018, Osinski et al., 2005; Thorsteinsson and Mayr, 1987). The SG2 study site covers an area of ~0.85 km² (Figs. 2A, D) and has an elevation range of 352 – 282 m ASL. Regional slope is gently NNW. The channel network covers an area of around 4 km², with a gentle slope of 3% towards the NW. The topology of this network is linear, highly directional and very crudely dendritic, containing a set of 15 tributaries oriented SSE-NNW, with heavy bedrock anastomosis in the headwaters that becomes less prominent downflow. The western section of SG2 is characterized by a series of pathways that originate as individual evenly spaced tributary channels (cf., Grau Galofre et al., 2018) that grade into anastomosing channels and finally into the main stem channel (Figs. 2D, E). SG2 tributary channels drain directly into a main stem channel oriented approximately S-N (Figs. 2D, E). A set of lateral meltwater channels are found perched on the western side of the main stem, offering a comparative between the morphology and characteristics of the two meltwater drainage styles. Lateral meltwater channels are loosely incised, are perched on the valley sides, have a much smaller scale than observed subglacial channels, and do not display signs of clear bedrock erosion and direct exposure on their course, particularly on their headwaters (see annotation in Fig. 2). The main stem of the subglacial network drains into a major canyon flowing northwards, and ultimately into Thomas Lee Inlet, a large fjord to the north (Fig. 2A). In its course, the major northern canyon splits west to feed the Haughton River Valley, running through the crater, and to the east towards a cluster of small lakes that are interconnected by small (m-scale) undulating channels, draining towards the fjord and Thomas Lee Inlet (Fig. 2A). We investigated the network of channels located at the western headwaters of SG2 over the course of the 2024 (Fig. 2D); the traverse covered a length of 1.2 km from the headwaters to the downstream reaches of a tributary channel. A portion of the main stem was characterized by means of the kinematic LiDAR scan. The main stem connects to large valleys to the north and south of the network (Fig. 2A). The headwaters of SG2 tributary channels (Fig. 4A, B) are morphologically distinct from the downstream reaches and main stem (Figs. 4 C, D, 4).

We characterized the granulometry and nature of the sediment cover and bedrock exposures found inside and around the channels. Near the channel headwaters (Figure 4), we observed channel beds consisting of barren (scabland) slabs of subhorizontal dolomite bedrock (Fig. 4A,B,D), indicating very high energy flows or flood conditions right at the onset. Bedrock outcrops occur along the channel walls, revealing an originally boxy shaped cross-section (Fig. 4A). Cryofracturing and consequent accumulation of talus at the base of channel walls was minimal at channel headwaters but became more prominent downstream to the point of obscuring the original shape of the channel wall shape (Fig. 4A, D). Hence, channel cross-section was notably box-like with near flat bedrock slabs making the channel floors, gently dipping (~5°) and steep walls at the headwaters (Fig. 4), evolving into trapezoidal shapes downflow as talus material accumulated (Fig. 5).
<figure><img src="image_22.png" /></figure>

At the headwaters, channels gently graded with the top of the plateaus (Fig. 4B), without displaying chutes or steep junctions that would later characterize tributary intersections. This smooth origin speaks of water flow continuity with the top of the plateaus, discarding water sourced from top-down melting of the ice sheet.

<figure><img src="image_23.png" /><figcaption>Figure 4. Geometry and nature of channel headwaters, black arrow indicates water flow direction. (A) Cross-section of a tributary looking upstream towards the origin, graded towards the plateaus. The original boxy shape of the cross-section is visible when considering the bedrock outcrops, notably to the right of the channel (top white arrow). Cryofractured scabland bedrock slabs are visible on the bed (arrows). (B) Channel headwaters grading into the plateaus, highlighting the in-situ dolomite slabs of clean bedrock (lower white arrow). (C) Irregular cross-sectional shape and longitudinal profile of a tributary, always bedrock carved. Blue arrow shows the position of a pothole. (D) Cross-section of the same tributary as (B), here looking downstream before a chuted junction (orange dashed line marks junction site). Note the bedrock slabs on the bed and the larger cross-section when compared to (B), no input of water from tributaries is noted between B and D.</figcaption></figure>

Approximately 500 m downstream from the origin, the channel deepened and widened into the main stem (Fig. 5A; the channel floor (Figs. 4, 5) became increasingly covered by cryofractured bedrock, which itself can be observed in outcrops along channel walls and at topographic highs (Fig. 5). Despite difficult digging conditions into cryofractured cobbles, we found evidence for glaciofluvially reworked, rounded and sub-rounded cobbles of ~10 cm in diameter under a ~15 cm veneer of cryofractured

260

265

270
<figure><img src="image_24.png" /></figure>

rocks, indicating glaciofluvial transport prior to the current channel cover. We observed very little evidence for sediment deposition in relation to subglacial channel networks in the study area. The only evidence for sediment deposition were sets of small (~1 – 10m scale) and localized patches of unconsolidated silt to silty sand, which are visible in Figure 5A as yellowish spots on the darker substrate that comprises the channel floor. However, we could not determine if they formed in situ or are the result of post-glacial deposition and periglacial reworking (i.e., loess, mass wasting, etc.), and it is, therefore, assumed that minimal to no channel deposition occurred before infill by post-glacial erosion. Current fluvial modification of the study area corresponds to small snowmelt channels with transport of fines and accumulation of small levees, but discharges are insufficient to transport cryofractured dolomite fragments (~10 cm), consistent with the low amounts of precipitation on Devon Island.

<figure><img src="image_25.png" /><figcaption>Figure 5. Field photos of SG2 and associated erosional surfaces. See locations of each photo in Figure 2D. A) Channel cross-section ~1 km down valley. Note size difference and level of degradation (from post-glacial cryofracturing) compared to headwaters. B)</figcaption></figure>

275

280

285
<figure><img src="image_26.png" /></figure>

Downstream channel reaches near main stem channel looking upstream from channel floor towards hanging valleys (black arrows). C) Hanging valley (black arrow) along wide and deep main stem channel looking downstream towards a section of significant anastomosis. D) Close up of a hanging valley.

290

SG2 in planform is highly intertwined (anastomosing), and channel junctions are usually defined by hanging valleys merging through steep chutes, which define slope breaks between different erosive surfaces. Consistent with these previous works, we observed slope breaks, akin to river knickpoints (e.g., Crosby & Whipple, 2006; Howard et al., 1994), in SG2 as hanging valleys where a channel at a higher elevation than the channel floor converged with the channel at a lower elevation (Figs. 5B-D). Slope break heights increased downflow from ~1 m near the headwaters to ~ 10 m at the downstream end of our field site (Fig. 5D), as channels became deeper. Graded channel junctions formed wide, flat areas that narrowed into a single channel downstream (Fig. 5B).

295

LiDAR DEMs (25 cm) of an ~ 0.67 km² area of SG2 (Fig. 6A) reveal well-defined slope breaks in undulating longitudinal profiles along select channels. These slope breaks typically define contacts between erosional surfaces and are characterized by a sharp drop in elevation along a channel. In the field, slope breaks occurred as hanging valleys or steeply graded channel junctions (Fig. 5B-D). In SG2, Profile 1 (Fig. 6B) was identified along a prominent incised channel adjacent to the main stem and oriented in the same ~ SW-NE direction, along the boundary between two erosional surfaces (ES2 and ES3, discussed below). This channel is ~ 310 m long with a concave longitudinal profile. Profile 1 contains slope breaks on both the downstream and upstream ends of the channel, which are steeper on the upstream (SW) side (Fig. 6B). Elevation changes for the south and north slope breaks are 10 m and 17 m, respectively. A longitudinal profile spanning the entire length of the traversed section of SG2 (Profile 2, Fig. 6C) exhibits localized concavity along the first ~ 775 m followed by a more linear, slightly convex profile downstream. Profile 2 displays at least two minor slope breaks up to 10 m in height (Fig. 6C) followed by relatively flat or low slope gradient areas.

300

305

310
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>
<figure><img src="image_33.png" /></figure>

**Figure 6.** A) Kinematic LiDAR-derived DEM of anastomosing channels in SG2. Relative location of DEM on satellite image in left panel (Google Earth MAXAR Imagery). B-C) Longitudinal profiles along two transects in heavily anastomosed section of SG2, corresponding to ES2. D) Erosional surfaces (ES) of SG2 overlain on LiDAR-based DEM (0.28 m/pixel resolution). E) ES overlain by 5 m contour intervals to illustrate relatively flat “terrace” areas, slope breaks, and channel incision into topographic highs. Red box in E) illustrates area coverage of the kinematic LiDAR scan that produced A).

Figure 6D-F shows the three erosive surfaces (ES) identified in SG2, separated by slope breaks in the LiDAR and ArcticDEM data. Erosive surface 1 (hereafter named ES1) is defined as the uppermost surface with an elevation range of 352 m – 312 m containing tributary bedrock channels incised into the top of the plateaus, oriented approximately SE-NW (Fig. 6D). ES1 channel dimensions are up to 82 m wide and 18 m deep, and anastomosis is common after a distance of ~300 m from channel headwaters. Most anastomosis occurs through obliquely connecting tributaries, delineating streamlined interfluves between channels, although one channel ~53 m long is incised perpendicular to the surrounding channels in ES1 (Fig. 6D, black box), perhaps following structural controls. Erosive surface 3 is defined as the bottom elevation surface, containing the main stem at the end of the network. The contact between ES1 and ES3 (Fig. 6D near black arrow), is defined through two slope breaks approximately 220 m apart (along the main stem) at 300 m and 308 m elevation (Fig. 6D, black arrow is pointing to 308 m slope break). A similar series of slope breaks at 300 m and 305 m define the contact between the middle erosive surface, ES2, and ES3 (Fig. 6D; corresponding to downstream slope break in Profile 1 (Fig. 6B); however, the contact between ES1 and ES2 is gradational. ES2 ranges in elevation from 324 m – 305 m. Channels in ES2 are up to 46 m wide and 11 m deep, relatively short, and anastomosis is common (Fig. 6D). Channels are also oriented roughly SW-NE, almost perpendicular to those in ES1 (Fig. 6D). ES3 contains the bottom of the sequence comprised by the main stem channel (Fig. 6D), which is a maximum of 148 m wide and 27 m deep. Anastomosis of this main stem occurs downstream from our study site, with interfluves that appear to lose the streamlined shape and become increasingly circular (Fig. 6D).

## 4.3 Subglacial channel network 3

Subglacial channel network 3 (SG3) is located approximately 16 km south of Thomas Lee Inlet and ~12 km to the southeast of the Haughton impact structure (Figs. 1A-B, 2A), incised into the Lower Allen Bay Formation dolomite (Osinski et al., 2005; Thorsteinsson and Mayr, 1987). SG3 covers an area of ~3.4 km² (Figs. 2A, F) and has an elevation range of 384 – 316 m ASL. The regional slope is gently NNW. SG3 is comprised of slightly curved and highly directional tributary channels oriented SE-NW, draining into a main stem oriented roughly S-N (Fig. 2F). The main stem of SG3 drains into a valley to the north before reaching a fjord that drains into Thomas Lee Inlet (Fig. 2A). The channels that make SG3 originate as evenly spaced tributary channels, with a slight funnel-shaped origin, (Grau Galofre et al., 2018) grading gently from the top of the plateaus. that display anastomosis at the headwaters and the first ~300 m downstream (Fig. 2F); channels grade into individual, evenly spaced channels before draining into the main stem channel. Two sets of lateral meltwater channels are perched on the northern and eastern sides of the main stem channels (Fig. 2F). These channels, similar to those observed perched on the main stem of
<figure><img src="image_34.png" /></figure>

SG2, are significantly smaller in length, width, and depth compared to the network's tributary channels and display minimal evidence of bedrock erosion (Fig. 2F). The two sets of lateral meltwater channels are oriented obliquely to each other and to the tributary channels. Field observations of SG3 are presented in Grau Galofre et al. (2018). Channel anastomosis develops predominantly in the western tributaries around 300 m downslope from the headwaters (Fig. 7A), defining streamlined bedrock interfluves between anastomosing tributaries (Fig. 7B) and exposed bedrock outcrops (Fig. 7C). Similar to SG2, bedrock exposures observed in the field are heavily degraded from post-glacial cryofracturing; consequently, channel floors and the bases of channel walls are covered by angular ~5 – 10 cm clasts of cryofractured dolomite fragments. Towards the main stem, channels appear to widen downstream of junctions (Figs. 7D-F), however, cross-section evolution is dominated by increasing depth rather than width downflow (Grau Galofre et al., 2018). Tributaries join the main stem through stepped or chuted junctions, defining a similar pattern of slope breaks in the longitudinal profile as in SG2.
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /><figcaption>Figure 7. Field photos of SG3. A) Wide and flat junction where anastomosing channels merge and diverge. B) Channel near junction in panel A, looking downstream. C) Exposed bedrock spire between two channels that diverged at panel A) and merge downstream beyond the bedrock spire. D) Flat-bottomed channel upstream of a channel junction, looking upstream. E) Wide, flat-bottomed</figcaption></figure>
<figure><img src="image_42.png" /></figure>

channel downstream of channel junction, looking downstream. F) Same channel as in panel E), looking upstream towards bedrock spire in panel C). Geologist in the distance circled in red for scale.

4.4 The Devon Island Labyrinth

The labyrinthine channels that make the so-called 'Devon Island Labyrinth' are located ~ 20 km southwest of the Thomas Lee Inlet and northwest of the 31 Ma-old 23 km-diameter Haughton impact structure (Osinski et al., 2005), within the faulted rim region of the crater (Figs. 1A-B, 2A). The area of study is ~ 4.7 km² with an elevation range of 274 – 160 m ASL and a gentle regional slope towards the north. The study area comprises a series of heavily anastomosing channels incised into the dolomite bedrock of the Upper Member Allen Bay Formation (interbedded massive resistant dark brown dolomite) within the Arctic Platform sedimentary succession (Fig. 8) (Osinski et al., 2005; Thorsteinsson and Mayr, 1987). The Devon Island Labyrinth headwaters are found in a local, relatively flat-lying depression of ~1.6 km wide by 2 km long (Fig. 8A "recessive unit"), comprised of silt- and sand-size dolomitic sediment from a recessive unit of the Lower Member Allen Bay formation, and punctuated by meter-scale erratic blocks (Fig. 9H). Anastomosing channels appear on the western side of the recessive unit (Fig. 8A). To the north, this unit connects to a cluster of lakes and eventually to a major fjord that drains into the western side of Thomas Lee Inlet (Fig. 8A). To the east, the system of heavily interconnected bedrock channels we refer to as the Devon Island Labyrinth drains to the SW-NE, before flowing into major valleys and fjords to the west and ultimately reaching Thomas Lee Inlet (Figs. 8A-C). Specifically, we did not visit the western system of anastomosing channels in the field or its relationship with the depression and were consequently unable to collect data for this specific location, and therefore field observations focus only on the eastern system. Channel disjunctions and anastomosing patterns start in this system less than 10 m downflow from headwaters, producing crudely streamlined and blocky interfluves that range in size from 1-10 m at the headwaters, evolving to larger inter-channel islands (100-500 m long) downflow.

365

370

375

380
<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /><figcaption>Figure 8. A) Satellite image (Google Earth/Maxar/Landsat) showing context of region surrounding the labyrinth terrain (white box) on Devon Island. B) Satellite image (Earthstar Geographics) of the labyrinth terrain. Channels are delineated in white. Magenta circles show photo locations for Figure 9. C) Digital Elevation Model of the labyrinth terrain (0.25 m/pixel resolution). For ease of comparison, the blue star in B) shows the same spot in C).</figcaption></figure>

As with other subglacial channels visited, channels in the Devon Island Labyrinth similarly appear to occupy three distinct levels or erosional surfaces (ES, discussed below) defined by distinct elevation ranges (Figs. 9A, B) that are defined by either stepped slopes or by steep channel junctions (Fig. 8D), which resemble hanging valleys observed in SG2. Channels

385

390
<figure><img src="image_45.png" /></figure>

in the upper erosional surface, corresponding to the top of the plateaus, are wide and grade gently into the depression described above (Figs. 9A-C) and channel interfluves are irregularly shaped topographic highs (Figs. 8C, 9B, 10A,E). In the middle erosive surface, channels are wider and deeper than those in the first (Figs. 9E, F). Channel floors are flat, and channels converge into smooth (graded) junctions that are particularly wide, flat areas with a distinct trapezoidal shape (Fig. 9F, G). Channels in the lowest erosional surface are the widest and deepest of those observed (Fig. 9G). Channel floors are covered by a veneer of angular, cryofractured fragments of dolomite bedrock. In-situ bedrock was observed near the tops of some interfluves, similarly to other subglacial channels (Fig. 9D). Some periglacial modification of channels floors is observed as sorted circles containing primarily cobble-sized clasts and silt, with some rare boulders.

395
<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /><figcaption>Figure 9. A) Overview of the Devon Island labyrinth terrain. B) Overview of the different erosive surfaces from the bottom of the main stem, the bottom of the tributary channels, and the top of the mesas. Black arrows denote examples of slope breaks. C) Wide and shallow flat-bottomed channel within the upper erosional surface (ES1; see discussion below). ATV tracks are approximately 1 metre wide for scale. Looking towards slope break to ES2. D) Example of shallow slope break or knickpoint in the Labyrinth (black arrow). Note bedrock outcrops at the top of the channel wall. E-F) Channel in middle erosional surface (ES2); note the difference in width and depth with channel in A. Geologists and ATV for scale. G) Lowest level erosional surface (ES3). Note progressive increase</figcaption></figure>

400

405
<figure><img src="image_55.png" /></figure>

in channel width and depth towards the lower erosional surfaces. H) Large boulder interpreted as a potential glacial erratic in the flat-lying area adjacent to the labyrinth terrain. Note absence of other landforms or boulders in the vicinity.

High-resolution DEMs (25 cm) produced from backpack kinematic LiDAR scans of ~0.36 km² of the Labyrinth (Fig. 10A) reveal well-defined slope breaks in longitudinal profiles. We once again take these as the defining distinction between erosional surfaces. In the Labyrinth, the most prominent slope break occurs at 227 m elevation (Profile 3, Fig. 10B) defining the contact between erosional surfaces ES2 and ES3. The channel connected to Profile 3 is relatively flat for ~200 m eastward before a near vertical drop of ~10 m in elevation into a wider and relatively flat section of the trunk channel (Figs. 10A, B). Another slope break in the Labyrinth (Profile 4, Fig. 10C) occurs at 230 m and drops ~20 m in elevation. The channel in Profile 4 displays an upward sloping concave longitudinal profile, is notably oriented almost perpendicular to that of Profile 3 (Fig. 10A) and is located along the same contact between erosional surfaces. The slope gradient of Profile 4 is steep but more gradual than Profile 3, and there is a m-scale depression below the slope break (Fig. 10C). Profile 5 (Fig. 10D) is an example of an incised interfluve between an anastomosed section of the trunk channel (Fig. 10A). The incision is ~40 m long and is oriented ~N-S. The southern elevation gradient is significantly steeper and deeper (5 m elevation change) than the northern side (3.5 m elevation change) (Fig. 10D) and connects to a wide and flat upstream section of the trunk channel (Fig. 10A). The southern side of Profile 5 is also situated at the contact between two erosional surfaces (ES2 and ES3).
<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /></figure>

<figure><img src="image_58.png" /><figcaption>Figure 10. A) Kinematic LiDAR-derived DEM of labyrinth terrain. Relative location of DEM on satellite image in top panel (Google Earth/MAXAR Imagery). B-D) Longitudinal profiles along three transects in the labyrinth. Note the prominent slope breaks observable in the profile and corresponding locations on the DEM. E) Erosional surfaces (ES) overlain on LiDAR-generated DEMs of the labyrinth terrain. F) DEM of study areas overlain by 5 m contour intervals to illustrate relatively flat “terrace” areas, slope breaks, and channel incision into topographic highs. Red box in F) illustrates area coverage of the kinematic LiDAR scan that produced A).</figcaption></figure>

Three erosive surfaces (ES) were identified in the Labyrinth based on the location of slope breaks (sharp knickpoints) in the longitudinal profiles as well as by observing flat surfaces in the topography, as shown in Figures 10E-F. The contact between ES1 and ES2 is gradational, whereas a very well-defined slope break is observed at the contact between ES2 and ES3, at 227 m elevation (Fig. 10E, black arrow). ES1 ranges in elevation from 274 m – 237 m and is primarily flat with a 400 m by 270 m area of anastomosed channels, which are oriented approximately E-W (Fig. 10E). These channels are ~5 – 10 m deep and ~25 – 30 m wide and appear to be more degraded or weathered when compared to other channels in the Labyrinth. ES1

425

430

435
<figure><img src="image_59.png" /></figure>

also appears as irregularly shaped plateaus within ES2 (Fig. 10E). The next surface, ES2, has an elevation range of 238 m – 227 m and spans the entire length of the Labyrinth (Fig. 10E). Significant anastomosis of channels is observed in the western half of ES2, where channels are approximately 80 – 120 m wide and ~ 7 – 15 m deep. Channel interfluves appear as crudely streamlined mesas lacking the elongation and definition of SG1 and SG2. Towards the eastern edge of ES2, lateral meltwater channels are identified as sets of short, parallel and straight channels merging on both sides of the main stem (Fig. 10E), significantly smaller in size and without the bedrock anastomosing observed in the Devon Island Labyrinth. The bottom surface ES3 ranges in elevation from 227 m – 160 m and is comprised of the main stem (Fig. 10E) that is ~ 30 – 37 m deep and ~ 70 – 120 m wide (widening downflow). The eastern end of ES3 contains short, shallow channels ~ 7 m deep and ~ 10 – 25 m wide that form broad, flat interfluves (Fig. 10E). These channels, while near those at the eastern end of ES2, do not resemble the identified lateral meltwater channels and appear morphologically distinct.

# 5 Discussion

## 5.1 Erosional Surfaces: Evidence for subglacial floods

Landscapes carved by glaciogenic meltwater flood events, subglacial and proglacial, include soil barren, scabland-type landscapes, mesa and butte topography, streamlined islands and interfluves, complex topologies including anastomosing and multiply connecting tributaries, hanging valleys and chuted junctions, and incision in bedrock. High-energy depositional features include large sandurss, streamlined tails following islands, megaripples, and boulder pavements (Baker, 2009; Denton and Sugden, 2005; Gombiner and Lesemann, 2024). In our study sites we observe several of these characteristics. Notably, anastomosing is noted near (~10-100 m) the channel headwaters, where channels additionally display scabland type bedrock incision exposing slabs of bedrock at the channel beds and formation of crudely streamlined interfluves. Anastomosis here refers to interconnected drainage pathways incised in bedrock, in contrast to fluvial anastomosis which typically occurs downstream of long-lived, stable rivers (Makaske, 2001). In this context, anastomosis has been proposed to record flood events that temporarily overwhelmed the drainage system, leading to the development of a complex landscape of interconnected channels (e.g., Gombiner & Lesemann, 2024; Lewis et al., 2006; O'Connor et al., 2020; Sugden et al., 1991). This morphology has been used to invoke the glacial flood hypothesis in well-known landscapes such as the Antarctic Dry Valleys (Fig. 11A) (e.g., Denton & Sugden, 2005; Lewis et al., 2006; Sugden et al., 1991) and proglacial megafloods in the Channeled Scablands of Washington, USA (Baker, 2009; Bretz, 1923). We thus interpret the bedrock-carved anastomosis at the headwaters of channels in this study (Figs. 2, 8) as evidence of high-energy flows and overwhelmed subglacial drainage systems.
Another characteristic of these channels, which stands also as a key difference to lateral meltwater channel systems, is the existence of multiple erosive surfaces, as described in section 3. Lewis et al. (2006) describe discrete erosional surfaces within The Labyrinth (Fig. 10A) that resemble those we documented here on Devon Island (compare Figs. 3, 6D-E, and 9E-F with
<figure><img src="image_60.png" /></figure>

Fig. 10B). These authors discuss the progressive incision of the landscape over multiple subglacial floods that occurred as a result of the repeated failure of subglacial meltwater dams along the margin of the East Antarctic Ice Sheet, associated with fluctuating climate during the Miocene.

Multiple, discreet erosive levels are also found in the Channeled Scablands, a landscape carved over the course of multiple proglacial Quaternary megafloods located in southeast Washington State (Baker, 2009; O'Connor et al., 2020). The Channeled Scablands are comprised of an anastomosing network of large, steep-walled troughs ("coulees") incised into basalt bedrock, forming streamlined interfluves, and are also characterized by hanging valleys of over 100 m in height formed at the junctions of tributary channels into larger valleys, due to progressive widening and deepening following catastrophic floods (Baker, 2009). These anastomosing channels have been interpreted to have formed by subglacial flooding as a result of rapid drainage of either supraglacial or subglacial meltwater reservoirs (Gombiner and Lesemann, 2024). Importantly, the Channeled Scablands channels contain adverse slopes, cross drainage divides, and connect smaller channels to larger trunk channels, similar to the morphologies observed in this study (e.g., Figs. 2, 8, 6B-C, 10B-E). The floods that formed the Channeled Scablands were likely large in magnitude ($10^4$ to > $10^6$ m³/s; Kirkham et al., 2019) and potentially occurred multiple times over multiple glaciations (Gombiner and Lesemann, 2024). Subglacial scabland erosion has also been modelled under the Cordilleran Ice Sheet at the Fraser Plateau (B.C., Canada) following the drainage of a supraglacial or subglacial meltwater reservoir wherein the existing subglacial drainage network could not accommodate the sudden influx of water (Burke et al., 2012). In this case, meltwater flowed along a broad zone of inefficient subglacial drainage, eroding the scabland before being funnelled into existing ice tunnels (Burke et al., 2012).
<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /><figcaption>Figure 11. Satellite image (Google Earth/Maxar) of anastomosing channels of A) The Labyrinth in Antarctica and B) labyrinth terrain on Devon Island.</figcaption></figure>

490

By analogy, the presence of distinct erosional surfaces, anastomosing channels, and slope breaks found both in the longitudinal profile and as tributary hanging valleys suggests that multiple erosive events were involved in the formation of the
<figure><img src="image_63.png" /></figure>

anastomosing portions of the Devon Island subglacial channel networks and labyrinthine channels (Figs. 3, 6, 8, 10). We suggest that multiple events of this nature resulted in the layering of an erosional landscape; that is, each flooding event developed or modified a distinct erosional surface that eventually formed a composite landscape of heavily anastomosed tributary channels and large, well-developed main channels (Fig. 12A). We suggest that the uppermost erosional surface, ES1, corresponds to the original top of the plateau in all three study areas. This surface in SG1 (Fig. 3A) and SG2 (Fig. 6D) contains the network's regularly spaced channel headwaters. Within ES2, flat-bottomed, bedrock incised and scabland-like anastomosing channels (Figs. 3, 6D-E, 10E-F), the occurrence of slope breaks at ES2 boundaries (Figs. 3A, 6B, 9B-D), and the large size of bedload cobbles found (Fig. 10H) all suggest that recurring high-energy, erosive subglacial meltwater flows were at play. Multiple high energy events would have led to well-developed anastomosis in older erosional surfaces by reusing meltwater pathways until overspill flow incises new channels in a younger ES (Fig. 12A). This suggests that the strongly anastomosing character of ES2 (Figs. 3, 6D-E, 9E-F) likely represents erosion during particularly high-energy flood events or that the boundary between ES1 and ES2 was not breached often, and floodwaters were confined to ES2. This may also represent a period of more common flooding events, perhaps related to a rapid retreat or collapse of the Innuitian ice sheet from the plateaus.
<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /><figcaption>Figure 12. A) Generalized cartoon illustrating the formation of erosional surfaces beginning with i) an incipient series of regularly spaced subglacial meltwater channels funneling subglacial meltwater (perhaps subglacial water sheets) from high-intensity water supply events, from which overbank flow inundates the surrounding bedrock surface, incising channels to form ES1. Further flooding events in ii) further develop ES1 (streamlined interfluves are shown color-coded with respect to erosive surface number) and eventually carve a new erosional surface (iii) once the incipient channel is sufficiently wide and deep to prevent overbank flow onto highest elevation surface (ES1). B) Generalized cartoon illustrating headward growth of tributary channels and associated slope breaks due to either a retreating ice margin or changing flow conditions at the point of meltwater input (e.g., episodic drainage of supraglacial lakes).</figcaption></figure>

We interpret ES3, which contains the main stems of SG1 and SG2 (Figs. 2B-F) and the trunk channel of the labyrinth (Figs. 8B, C), as a surface containing developed subglacial meltwater channels. These erosional surfaces contain the widest and deepest channels in each study area and likely captured meltwater from the surrounding smaller channels (e.g., Fig. 8D). While channels in The Labyrinth and Channeled Scablands are larger than those documented in this study (> 100 m wide), the

510

515

520
<figure><img src="image_66.png" /></figure>

anastomosing character and apparent funnelling of meltwater into larger valleys is directly comparable to the configuration of the labyrinth terrain on Devon Island (Fig. 10), as well as localized areas of anastomosis in channel networks (Figs. 2, 3, 6, 7). This leads us to conclude that anastomosis at our study sites on Devon Island (Figs. 2, 8) is likely associated with subglacial flood events, perhaps over multiple episodes, that consistently incised new channels as the drainage system was overwhelmed and enlarged main or trunk channels by funnelling floods from areas of anastomosis to wider, deeper channels (Fig. 12A). While the extent of upward erosion into the ice is unknown, it is possible that channel anastomosis is more localized in SG1 and SG2 (Figs. 3A, 6E) compared to the Devon Island labyrinth due to a higher number of pre-existing pathways for efficient meltwater drainage (i.e., tributary channels: Figs. 2B, D, F). This may also provide an explanation for the clearer “layering” of erosional surfaces in the labyrinth terrain (Figs. 10E) where flooding was focussed consistently along the trunk channel (Fig. 12A) rather than spread throughout a channel network. However, it is also possible that the Labyrinth reflects a better-developed, more mature flood landscape, either because it has sustained more flooding cycles or because it has undergone multiple glacial cycles. Indeed, in terms of evolution, the striking regular spacing, funnel-shaped heads, smooth grading into the top of the plateaus, and consistent elevation across all tributaries in a same network (SG1, SG2, SG3, and the Labyrinth) imply that tributaries formed synchronously and competed for water supply (e.g., Hewitt et al., 2011), potentially delivered as subglacial sheet flow events (see below). Development of erosive surfaces, as shown in Fig. 12, requires flow water input into all tributaries to fluctuate between flood, overspill conditions leading to anastomosis and quieter conditions concentrating flow into the main stem. At a longer evolutionary scale, the main stem and tributary channels may have either formed during a single glacial cycle or during the course of multiple glaciations and continued to widen and deepen over time following Innuitian advances and retreats (Dalton et al., 2022).

Today, Devon Island's interior plateau is incised with many channel networks like SG1 and SG2 (Figs. 1B, 2A), suggesting that this process was not localized but appears to have occurred across the entire plateau. The Labyrinth presents a distinct morphology from SG1-SG3 that is worth noting. The notable absence of tributary channels draining into the trunk channel in the labyrinth (Fig. 8B) may be a result of a lack of seasonal meltwater input in this area; rather, the labyrinth may have experienced more frequent or larger magnitude flooding followed by periods of relative inactivity of the subglacial drainage system and drainage closure (i.e., minimal seasonal meltwater input). Additionally, the enlarged dimensions of main stem and trunk channels may represent a convergence to an equilibrium geometry that has been interpreted to represent subglacial channel formation over multiple glacial cycles (e.g., Kirkham et al., 2020).

### 5.1.2 Potential meltwater sources

Previous work at SG1 and SG2 has suggested that meltwater was sourced from supraglacial reservoirs and routed to the ice-bed interface via moulins or crevasses (Grau Galofre et al., 2018, Ruso et al., 2024). Episodic drainage of supraglacial lakes to the bed is commonly observed in other glaciated environments such as Greenland (e.g., Bartholomew et al., 2012; Das
<figure><img src="image_67.png" /></figure>

et al., 2008; Zwally et al., 2002) and Antarctica (e.g., Langley et al., 2016; Stokes et al., 2019; Tuckett et al., 2019). Several characteristics in SG1, SG2, and SG3 hint to their origin. First, the regular spacing between all tributaries in a network implies that all tributaries formed synchronously, sourcing from the same water reservoir and competing with each other for meltwater supply. Regular spacing has been observed in eskers across Canada (e.g., Storrar et al., 2015), and the regular spacing in synchronous subglacial drainage routes is theoretically predicted (e.g., Schoof, 2010, Hewitt et al., 2011). Next, the funnel morphology of channel heads at the headwaters (see Fig 5, panels A and D), and the smooth continuation of the elevation profile from the top of the plateaus into the channels (see grading in figure 4) seems to suggest that subglacial channelized drainage self-organized sourcing from a subglacial sheet flow, with funnel heads potentially corresponding to the transitional morphology between distributed and channelized subglacial drainage, eroded backwards as channel formation progressed. Headwaters anastomosing incision into bedrock, and cobble-sized bedloads, both suggest large discharge, flood events: it is unclear whether these correspond to punctuated but cyclical events (i.e., seasonal) or give an insight on how the otherwise cold-based Innuitian ice sheet may have collapsed without transitioning into a dynamic, warm-based ice mass for which there is no morphological evidence (e.g., Dyke, 1999, Grau Galofre et al., 2018, Ruso et al., 2024), but water sources feeding a sheet flow were likely supraglacial, perhaps through moulins or crevasses.

The labyrinth on Devon Island could have also formed via subglacial flooding sourced from supraglacial meltwater; however, the presence of a depression at the western headwaters of the labyrinth (Fig. 8A) presents an alternative hypothesis. The recessive unit in the Allen Bay Formation (Lower Member) is more readily eroded than the surrounding massive thickly-bedded dolomite bedrock (Upper Member), offering a naturally low-lying depression that could have accumulated subglacial meltwater. Evidence for this lake comes in the topography of the basin divide, which is breached in five different locations at different elevations and connected to cross-divide channels flowing to the WSW and NE, and E (Fig. 8A). This is inconsistent with a proglacial lake setting, where once a breach is emplaced at a given elevation, further breaching in different directions is unlikely as the original breach defines the low point and outlet. In this way, there is potential for formation of the labyrinth terrain by impoundment and release of meltwater from a subglacial reservoir, analogous to the Antarctic Labyrinth (Lewis et al., 2006; Sugden et al., 1991).

## 5.2 Slope Breaks: Evidence for changing subglacial flow properties

Slope breaks in the studied subglacial channel networks (Figs. 3A, 6B-C) and labyrinth terrain (Figs. 10B-D) are interpreted to indicate pulses of erosion as the result of a change in flow dynamics (meltwater overpressure or supply rate) or glacial dynamics (ice margin position). These slope breaks are associated with relatively flat stretches of channel floor referred to as “terraces” (Figs. 3A, 6D, 10E) due to the resemblance to fluvial strath terraces, erosional terraces carved into bedrock (Antoine et al., 2007; Schanz et al., 2018; Wegmann and Pazzaglia, 2009). We suggest the slope breaks described would have worked as fluvial knickpoints do, albeit probably reflecting abrupt changes in water erosive capability rather than related to long-term
<figure><img src="image_68.png" /></figure>

tectonics or climate signals (e.g., Boulton, 2020; Crosby & Whipple, 2006; Marrucci et al., 2018; Whittaker & Boulton, 2012)
In rivers, knickpoint formation can occur due to base level change or tectonics (e.g., Wegmann and Pazzaglia, 2009; Whittaker
and Boulton, 2012), due to variations in discharge (climate pulses, e.g., Leonard and Whipple, 2021), changes in substrate
composition, or due to stream capture (Crosby and Whipple, 2006).

In subglacial environments, temporal changes in driving stresses could very well include, in addition to those cited above, varying effective pressure gradients following changes in meltwater input rate, changing drainage conditions coupled to ice sliding dynamics (Schoof, 2010), changing organization of the drainage system (connectivity, tortuosity), and evolving ice margin position, to cite but a few possible processes. For example, the development of large-scale hanging valleys (> 100 m above valley floors) in the Channeled Scablands is attributed to the enlargement of main valleys by flooding events that rapidly dropped valley floor elevation (Baker, 2009; Lehnigk and Larsen, 2022). Repeated high-discharge events significantly increased the shear stress of the flow which resulted in progressive deepening and headward migration of these hanging valleys (Lehnigk and Larsen, 2022). Channels at all three locations in this study feed into wide and deep channels downstream (main stems of SG1 and SG2, trunk channel of labyrinth terrain; Figs. 2C, 8C, E, F) and, while smaller in magnitude, hanging valleys are observed in SG1 and SG2 at the junctions of tributary channels to the main stem (Figs. 3A, 5B-D, 6E). Capture of meltwater into increasingly larger, lower elevation pathways (e.g., Fig. 12B) therefore provides a potential trigger for local base level changes and consequent fluctuations in flow conditions.

To form slope breaks in a glacial environment, perturbation to the system would either occur as a result of changing ice margin position, equivalent to a changing base level in a river, or at the point of meltwater input (e.g., moulin or crevasse) (Fig. 12B), as opposed to propagating from the downstream reaches as with fluvial systems (Howard et al., 1994). Burke et al. (2012) describe a scenario that could have led to the development of a meltwater corridor under the Cordilleran Ice Sheet during the LGM (their Figure 9). This scenario invokes rapid drainage of a supraglacial lake via moulins or crevasses, resulting in the formation of scabbed terrain and headward erosion of subglacial “knickpoints” at the headwaters of pre-existing (pre-flood) meltwater channels (Fig. 7A), while downstream focussing into efficient drainage pathways increased the velocity of the flow (Burke et al., 2012). The shallow box-like channel headwaters of SG1 and SG2 (Fig. 4) provide evidence for headward growth of tributary channels and funnelling of meltwater into progressively wider and deeper channels (e.g., main stem and trunk channels). Moreover, anastomosis and slope breaks in SG1, SG2, and the labyrinth (Figs. 3, 6, 10) offer further evidence for multiple events where the drainage system was temporarily overwhelmed.

A second possible explanation for the formation of slope breaks can be envisaged using the modelling by Braun et al. (1999), who predicted that within a series of glacial-interglacial cycles glacial erosion increases to a steady-state while fluvial erosion increases with time. If we assume negligible glacial erosion in the interior plateau region of Devon Island while occupied by a dominantly cold-based ice sheet (Fig. 1C) – or that it had reached a terminal/maximum glacial erosion rate largely confined to valleys while ice was frozen to divides and interfluves – then fluvial erosion becomes increasingly more prominent in
<figure><img src="image_69.png" /></figure>

shaping the landscape. Fluvial erosion decreases during glaciations but increases during interglacial periods (Jansen et al., 2011), incising existing glacial features such as U-shaped valleys and moraines (Braun et al., 1999). This pulse in fluvial erosion, potentially enhanced by the mobilization of unconsolidated glacial sediment, results in “glacial knickpoint retreat” (Braun et al., 1999). This pulse is likely short-lived as post-glacial “knickpoint” retreat slows as sediment is evacuated and the erosive power of the flow is consequently reduced (Jansen et al., 2011). Moreover, during periods of significant climatic transition such as glacial-interglacial transitions, 2 – 3 m of incision can occur on a timescale of less than 1 ka (Antoine, 1997; Vandenberghe, 1995). Since fluvial knickpoints translate changes in downstream flow to the upstream reaches as a propagating feature, the lack of knickpoints in a river suggests that it is in equilibrium or steady-state flow (i.e., no significant perturbations to the system). Therefore, preservation of slope breaks in our study implies short-lived events or changes in flow conditions that did not have sufficient time or energy to propagate entirely through the system.

In this way, we suggest that the slope breaks identified in the subglacial channel networks (Figs. 3A, 6B, D) and labyrinth terrain (Figs. 10B-E) on Devon Island represent “subglacial knickpoints,” the erosive response to periods or pulses of intense subglaciofluvial erosion responding to changing flow conditions, which can include anything from changing ice margins, time varying meltwater input rates, or subglacial drainage rerouting. It is worth noting that these features could have survived multiple glaciations and could have formed at the beginning of a glacial period or over multiple glacial-interglacial cycles where topography is inherited from the previous glacial or interglacial period and then subsequently modified. Indeed, Devon Island has experienced multiple cycles of glaciation and inter-glacial periods (Dalton et al., 2022) and, as a result, has been modified by both glacial, fluvial, and periglacial processes, and it is possible that the landscapes here described are palimpsests of this long history.

## 6 Conclusions

In this study we present detailed morphological characterization of bedrock subglacial drainage channels on Devon Island, including the first description of subglacial labyrinthine channels on Devon Island, Nunavut, Canada, drawn from field and remote sensing data are presented. Morphological analysis of the labyrinthine channels and similarly highly anastomosing sections of two bedrock subglacial channel networks suggests they formed over multiple episodes of erosive, high-energy flows. In particular, distinct erosional surfaces identified by the presence of slope breaks in longitudinal profiles, as well as highly anastomosing bedrock channels, both suggest formation by short-lived, locally intense episodes of subglacial flood erosion. The presence of well-defined erosional surfaces, particularly in the labyrinth terrain, suggests floods progressively incised into lower elevations where meltwater was captured by pre-existing or incipient main stem or trunk channels.

Comparison to The Labyrinth (subglacial flood landscape, Wright glacier valley, Antarctica) and the Channeled Scablands (proglacial flood landscape, Washington, USA), two well-known glacial outburst flood systems, reveals a resemblance in the anastomosing character of channels. This suggests that anastomosis of channels in SG1 and SG2 and significant anastomosis
<figure><img src="image_70.png" /></figure>

throughout the labyrinth terrain, is also likely associated with subglacial outburst floods. The timing and magnitude of these floods is currently unknown; however, based on channel dimensions they are associated with relatively lower-magnitude flooding compared to The Antarctic Labyrinth and the Channeled Scablands. Slope breaks identified in all three study areas are comparable to hanging valleys observed in the Channeled Scabland that formed as a result of glacial outburst flooding and interpreted as a result of subglacial flooding and consequent headward growth of channels. In comparison to fluvial knickpoints, we suggest that the presence of slope breaks in longitudinal profiles indicate short-lived erosive events as the pulse of erosion or headward channel growth did not fully propagate through the entire system.

## Data Availability

Data for this project is available upon request from the lead author, Simona F. Ruso (sruso@uwo.ca).

## Author contributions

SFR, GRO, and AGG devised and discussed project idea and timeline. SFR and GRO collected drone data and field observations in 2022. SFR, GRO, AGG, and AK collected drone data and field data in 2024. AK compiled and developed DEMs using LiDAR data from kinematic LiDAR scans. SFR compiled and analysed drone data, field notes and photos, and developed DEMs, obtained measurements, and carried out morphometric analyses. Literature review and manuscript writing was carried out by SFR. AK provided detailed description of kinematic LiDAR methodology. GRO and AGG provided insight and edits throughout manuscript writing.

## Competing interests

The authors declare that they have no conflict of interest.

## Acknowledgements

Funding was provided by a Canadian Space Agency (CSA) Flights and Fieldwork for the Advancement of Science and Technology (FAST) grant, and a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant Northern Supplement to GRO. SFR thanks the Northern Scientific Training Program (NSTP), The Explorers Club, and the Royal Canadian Geographical Society for funding. A.G.G was supported by the European Commission in the program H2020 under the MSCA grant agreement n° 101027900, the Institut national des sciences de l'Univers – PNP project EGV E&M, and the French polar institute IPEV, under project SUBGLACE-2. Logistical support from the Polar Continental Shelf Program
<figure><img src="image_71.png" /></figure>

(NRCan) is gratefully acknowledged. SFR thanks Dr. Sam Anderson (SFU) and Axel Noblet (UWO) for their assistance with field work.

## References

Alley, R. B., Cuffey, K. M., and Zoet, L. K.: Glacial erosion: Status and outlook, Ann Glaciol, 60, 1–13, https://doi.org/10.1017/aog.2019.38, 2019.

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

Antoine, P.: Évolution Tardiglaciaire et début Holocène des vallées de la France septentrionale: nouveaux résultats, Comptes Rendus de l'Académie des Sciences-Series IIA-Earth and Planetary Science, 325, 35–42, 1997.

Antoine, P., Limondin Lozouet, N., Chaussé, C., Lautridou, J. P., Pastre, J. F., Auguste, P., Bahain, J. J., Falguères, C., and Galehb, B.: Pleistocene fluvial terraces from northern France (Seine, Yonne, Somme): synthesis, and new results from interglacial deposits, Quat Sci Rev, 26, 2701–2723, https://doi.org/10.1016/j.quascirev.2006.01.036, 2007.

Baker, V. R.: The channeled scabland: A retrospective, https://doi.org/10.1146/annurev.earth.061008.134726, May 2009.

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

Beaud, F., Venditti, J. G., Flowers, G. E., and Koppes, M.: Excavation of subglacial bedrock channels by seasonal meltwater flow, Earth Surf Process Landf, 43, 1960–1972, https://doi.org/10.1002/esp.4367, 2018.

Boulton, S. J.: Geomorphic Response to Differential Uplift: River Long Profiles and Knickpoints From Guadalcanal and Makira (Solomon Islands), Front Earth Sci (Lausanne), 8, https://doi.org/10.3389/feart.2020.00010, 2020.

Braun, J., Zwartz, D., and Tomkin, J. H.: A new surface-processes model combining glacial and fluvial erosion, Ann Glaciol, 28, 1999.

Bretz, J. H.: Glacial drainage on the Columbia Plateau, Bulletin of the Geological Society of America, 34, 573–608, 1923.

Burke, M. J., Brennand, T. A., and Perkins, A. J.: Evolution of the subglacial hydrologic system beneath the rapidly decaying Cordilleran Ice Sheet caused by ice-dammed lake drainage: Implications for meltwater-induced ice acceleration, Quat Sci Rev, 50, 125–140, https://doi.org/10.1016/j.quascirev.2012.07.005, 2012.

Chandler, D. M., Wadham, J. L., Lis, G. P., Cowton, T., Sole, A., Bartholomew, I., Telling, J., Nienow, P., Bagshaw, E. B., Mair, D., Vinen, S., and Hubbard, A.: Evolution of the subglacial drainage system beneath the Greenland Ice Sheet revealed by tracers, Nat Geosci, 6, 195–198, https://doi.org/10.1038/ngeo1737, 2013.
<figure><img src="image_72.png" /></figure>

Crosby, B. T. and Whipple, K. X.: Knickpoint initiation and distribution within fluvial networks: 236 waterfalls in the Waipaoa River, North Island, New Zealand, Geomorphology, 82, 16–38, https://doi.org/10.1016/j geomorph.2005.08.023, 2006.

705

Cuffey, K. M. and Paterson, W. S. B.: The Physics of Glaciers, Fourth., edited by: Butterworth-Heinemann, Elsevier, 2010.

Dalton, A. S., Stokes, C. R., and Batchelor, C. L.: Evolution of the Laurentide and Innuitian ice sheets prior to the Last Glacial Maximum, Earth Sci Rev, 224, 103875, 2022.

710

Das, S. B., Joughin, I., Behn, M. D., Howat, I. M., King, M. A., Lizarralde, D., and Bhatia, M. P.: Fracture propagation to the base of the Greenland Ice Sheet during supraglacial lake drainage, Science (1979), 320, 778–781, https://doi.org/10.1126/science.115336, 2008.

Denton, G. H. and Sugden, D. E.: Meltwater features that suggest Miocene ice-sheet overriding of the Transantarctic Mountains in Victoria Land, Antarctica, Geografiska Annaler, Series A: Physical Geography, 87, 67–85, https://doi.org/10.1111/j.0435-3676.2005.00245.x, 2005.

715

Dewald, N., Livingstone, S. J., and Clark, C. D.: Subglacial meltwater routes of the Fennoscandian Ice Sheet, J Maps, 18, 382–396, https://doi.org/10.1080/17445647.2022.2071648, 2022.

Dyke, A. S.: Last Glacial Maximum and deglaciation of devon island, arctic canada: support for an innuitian ice sheet, Quat Sci Rev, 18, 393–420, https://doi.org/10.1016/S0277-3791(98)00005-5, 1999.

720

England, J., Atkinson, N., Bednarski, J., Dyke, A. S., Hodgson, D. A., and Ó Cofaigh, C.: The Innuitian Ice Sheet: configuration, dynamics and chronology, Quat Sci Rev, 25, 689–703, https://doi.org/10.1016/j.quascirev.2005.08.007, 2006.

Flowers, G. E.: Modelling water flow under glaciers and ice sheets, https://doi.org/10.1098/rspa.2014.0907, 8 April 2015.

Gombiner, J. and Lesemann, J.-E.: Okanogan lobe tunnel channels and subglacial floods into Moses Coulee, Channeled Scabland, northwestern United States, Geology, https://doi.org/10.1130/g52005.1, 2024.

725

Grau Galofre, A., Mark Jellinek, A., Osinski, G. R., Zanetti, M., and Kukko, A.: Subglacial drainage patterns of Devon Island, Canada: Detailed comparison of rivers and subglacial meltwater channels, Cryosphere, 12, 1461–1478, https://doi.org/10.5194/tc-12-1461-2018, 2018.

Hewitt, I. J.: Modelling distributed and channelized subglacial drainage: The spacing of channels, Journal of Glaciology, 57, 302–314, https://doi.org/10.3189/002214311796405951, 2011.

Hogan, K. A., Arnold, N. S., Larter, R. D., Kirkham, J. D., Noormets, R., Ó Cofaigh, C., Golledge, N. R., and Dowdeswell, J. A.: Subglacial Water Flow Over an Antarctic Palaeo-Ice Stream Bed, J Geophys Res Earth Surf, 127, e2021JF006442, https://doi.org/10.1029/2021JF006442, 2022.

730

Howard, A. D., Dietrich, W. E., and Seidl, M. A.: Modeling fluvial erosion on regional to continental scales, J Geophys Res, 99, 971–984, 1994.

Jansen, J. D., Fabel, D., Bishop, P., Xu, S., Schnabel, C., and Codilean, A. T.: Does decreasing paraglacial sediment supply slow knickpoint retreat?, Geology, 39, 543–546, https://doi.org/10.1130/G32018.1, 2011.

735
<figure><img src="image_73.png" /></figure>

Jansen, J. D., Codilean, A. T., Stroeven, A. P., Fabel, D., Hättestrand, C., Kleman, J., Harbor, J. M., Heyman, J., Kubik, P. W., and Xu, S.: Inner gorges cut by subglacial meltwater during Fennoscandian ice sheet decay, Nat Commun, 5, 3815, https://doi.org/10.1038/ncomms4815, 2014.

Jørgensen, F. and Sandersen, P. B. E.: Buried and open tunnel valleys in Denmark-erosion beneath multiple ice sheets, Quat Sci Rev, 25, 1339–1363, https://doi.org/10.1016/j.quascirev.2005.11.006, 2006.

Kamb, B.: Sliding motion of glaciers: theory and observation, Reviews of Geophysics, 8, 673–728, 1970.

Kehew, A. E., Piotrowski, J. A., and Jørgensen, F.: Tunnel valleys: Concepts and controversies - A review, Earth Sci Rev, 113, 33–58, https://doi.org/10.1016/j.earscirev.2012.02.002, 2012.

Kirkham, J. D., Hogan, K. A., Larter, R. D., Arnold, N. S., Nitsche, F. O., Golledge, N. R., and Dowdeswell, J. A.: Past water flow beneath pine Island and Thwaites glaciers, west Antarctica, Cryosphere, 13, 1959–1981, https://doi.org/10.5194/tc-13-1959-2019, 2019.

Kirkham, J. D., Hogan, K. A., Larter, R. D., Arnold, N. S., Nitsche, F. O., Kuhn, G., Gohl, K., Anderson, J. B., and Dowdeswell, J. A.: Morphometry of bedrock meltwater channels on Antarctic inner continental shelves: Implications for channel development and subglacial hydrology, Geomorphology, 370, 107369, https://doi.org/10.1016/j geomorph.2020.107369, 2020.

Kleman, J. and Borgström, I.: Reconstruction of palaeo-ice sheets: The use of geomorphological data, Earth Surf Process Landf, 21, 893–909, https://doi.org/10.1002/(SICI)1096-9837(199610)21:10<893::AID-ESP620>3.0.CO;2-U, 1996.

Langley, E. S., Leeson, A. A., Stokes, C. R., and Jamieson, S. S. R.: Seasonal evolution of supraglacial lakes on an East Antarctic outlet glacier, Geophys Res Lett, 43, 8563–8571, https://doi.org/10.1002/2016GL069511, 2016.

Lee, P., Rice, J. W. Jr., Bunch, T. E., Grieve, R. A. F., Mckay, C. P., Schutt, J. W., and Zent, A. P.: Possible analogs for small valleys on Mars at the Haughton Impact Crater site, Devon Island, Canadian High Arctic, in: Lunar and Planetary Science Conference, https://doi.org/20000092052, 1999.

Lehnigk, K. E. and Larsen, I. J.: Pleistocene Megaflood Discharge in Grand Coulee, Channeled Scabland, USA, J Geophys Res Earth Surf, 127, https://doi.org/10.1029/2021JF006135, 2022.

Lelandais, T., Ravier, É., Pochat, S., Bourgeois, O., Clark, C., Mourgues, R., and Strzerzynski, P.: Modelled subglacial floods and tunnel valleys control the life cycle of transitory ice streams, Cryosphere, 12, 2759–2772, https://doi.org/10.5194/tc-12-2759-2018, 2018.

Lewington, E. L. M., Livingstone, S. J., Clark, C. D., Sole, A. J., and Storrar, R. D.: A model for interaction between conduits and surrounding hydraulically connected distributed drainage based on geomorphological evidence from Keewatin, Canada, Cryosphere, 14, 2949–2976, https://doi.org/10.5194/tc-14-2949-2020, 2020.
<figure><img src="image_74.png" /></figure>

Lewis, A. R., Marchant, D. R., Kowalewski, D. E., Baldwin, S. L., and Webb, L. E.: The age and origin of the Labyrinth, western Dry Valleys, Antarctica: Evidence for extensive middle Miocene subglacial floods and freshwater discharge to the Southern Ocean, Geology, 34, 513–516, https://doi.org/10.1130/G22145.1, 2006.

Makaske, B.: Anastomosing rivers: a review of their classification, origin and sedimentary products, Earth Sci Rev, 53, 149–196, https://doi.org/10.1016/S0012-8252(00)00038-6, 2001.

Marrucci, M., Zeilinger, G., Ribolini, A., and Schwanghart, W.: Origin of knickpoints in an alpine context subject to different perturbing factors, stura valley, maritime alps (North-Western Italy), Geosciences (Switzerland), 8, https://doi.org/10.3390/geosciences8120443, 2018.

Nienow, P. W., Sharp, M., and Willis, I. C.: Velocity-Discharge Relationships Derived from Dye Tracer Experiments in Glacial Meltwaters: Implications for Subglacial Flow Conditions, Hydrol Process, 10, 1411–1426, 1996.

Nye, J. F.: Water flow in glaciers: jokulhlaups, tunnels and veins, Journal of Glaciology, 17, 181–207, https://doi.org/10.3189/S002214300001354X, 1976.

Ó Cofaigh, C.: Tunnel valley genesis, Prog Phys Geogr, 20, 1–19, 1996.

O'Connor, J. E., Baker, V. R., Waitt, R. B., Smith, L. N., Cannon, C. M., George, D. L., and Denlinger, R. P.: The Missoula and Bonneville floods—A review of ice-age megafloods in the Columbia River basin, https://doi.org/10.1016/j.earscirev.2020.103181, 1 September 2020.

Osinski, G. R., Lee, P., Spray, J. G., Parnell, J., Lim, D. S. S., Bunch, T. E., Cockell, C. S., and Glass, B.: Geological overview and cratering model for the Haughton impact structure, Devon Island, Canadian High Arctic, Meteorit Planet Sci, 40, 1759–1776, https://doi.org/10.1111/j.1945-5100.2005.tb00145.x, 2005.

Ruso, S. F., Grau Galofre, A., and Osinski, G. R.: Geomorphological and morphometric characterization of subglacial channels on Devon Island, Nunavut, Canada, Geomorphology, 462, https://doi.org/10.1016/j.geomorph.2024.109345, 2024.

Schanz, S. A., Montgomery, D. R., Collins, B. D., and Duvall, A. R.: Multiple paths to straths: A review and reassessment of terrace genesis, Geomorphology, 312, 12–23, https://doi.org/10.1016/j.geomorph.2018.03.028, 2018.

Schoof, C.: Ice-sheet acceleration driven by melt supply variability, Nature, 468, 803–806, https://doi.org/10.1038/nature09618, 2010.

Schumm, S. A., Erskine, W. D., and Tilleard, J. W.: Morphology, hydrology, and evolution of the anastomosing Ovens and King Rivers, Victoria, Australia, Geol Soc Am Bull, 108, 1212–1224, https://doi.org/10.1130/0016-7606(1996)108<1212:MHAEOT>2.3.CO;2, 1996.

Shaw, J.: The meltwater hypothesis for subglacial bedforms, Quaternary International, 90, 5–22, 2002.

Simkins, L. M., Greenwood, S. L., Winsborrow, M. C. M., Bjarnadóttir, L. R., and Lepp, A. P.: Advances in understanding subglacial meltwater drainage from past ice sheets, Ann Glaciol, 1–5, https://doi.org/10.1017/aog.2023.16, 2023.
<figure><img src="image_75.png" /></figure>

Stokes, C. R., Sanderson, J. E., Miles, B. W. J., Jamieson, S. S. R., and Leeson, A. A.: Widespread distribution of supraglacial lakes around the margin of the East Antarctic Ice Sheet, Sci Rep, 9, 13823, https://doi.org/10.1038/s41598-019-50343-5, 2019.

Sugden, D. E., Denton, G. H., and Marchant, D. R.: Subglacial meltwater channel systems and ice sheet overriding, Asgard Range, Antarctica, Geografiska Annaler, Series A, 73 A, 109–121, https://doi.org/10.1080/04353676.1991.11880335, 1991.

Thorsteinsson, R. and Mayr, U.: The sedimentary rocks of Devon Island, Canadian arctic archipelago, Supply and Services, 1987.

Tuckett, P. A., Ely, J. C., Sole, A. J., Livingstone, S. J., Davison, B. J., Melchior van Wessem, J., and Howard, J.: Rapid accelerations of Antarctic Peninsula outlet glaciers driven by surface melt, Nat Commun, 10, 4311, https://doi.org/10.1038/s41467-019-12039-2, 2019.

Vandenberghe, J.: Timescales, climate and river development, Quat Sci Rev, 14, 631–638, 1995.

Van der Vegt, P., Janszen, A., and Moscariello, A.: Tunnel valleys: current knowledge and future perspectives, Geological Society, London, Special Publications, 368, 75–97, 2012.

Wegmann, K. W. and Pazzaglia, F. J.: Late Quaternary fluvial terraces of the Romagna and Marche Apennines, Italy: Climatic, lithologic, and tectonic controls on terrace genesis in an active orogen, Quat Sci Rev, 28, 137–165, https://doi.org/10.1016/j.quascirev.2008.10.006, 2009.

Whittaker, A. C. and Boulton, S. J.: Tectonic and climatic controls on knickpoint retreat rates and landscape response times, J Geophys Res Earth Surf, 117, https://doi.org/10.1029/2011JF002157, 2012.

Zwally, H. J., Abdalati, W., Herring, T., Larson, K., Saba, J., and Steffen, K.: Surface melt-induced acceleration of Greenland ice-sheet flow, Science (1979), 297, 218–222, https://doi.org/10.1126/science.1072708, 2002.