Large-scale integrated subglacial drainage around the former Keewatin Ice Divide, Canada reveals interaction between distributed and channelised systems

We identify and map traces of subglacial meltwater drainage around the former Keewatin Ice Divide, Canada from ArcticDEM data. Meltwater tracks, tunnel valleys and esker splays exhibit several key similarities, including width, spacing, their 15 association with eskers and transitions to and from different types, which together suggest they form part of an integrated drainage signature. We collectively term these features ‘meltwater corridors’ and propose a new model for their formation, based on observations from contemporary ice masses, of pressure fluctuations surrounding a central conduit. We suggest that eskers record the imprint of a central conduit and 20 meltwater corridors the interaction with the surrounding distributed drainage system. The widespread aerial coverage of meltwater corridors (5-36 % of the bed) provides constraints on the extent of basal uncoupling induced by basal water pressure fluctuations and variations in spatial distribution and evolution of the subglacial drainage system, which will modulate the ice dynamic response. 25


Introduction 30
Variations in the configuration of subglacial hydrological systems are key to understanding some of the most dynamic ice sheet behaviour at a range of spatial and temporal scales (e.g. Zwally et al., 2002;Das et al., 2008;Joughin et al., 2008;van de Wal et al., 2008;Shepherd et al., 2009;Palmer et al., 2011;Fitzpatrick et al., 2013;35 Doyle et al., 2014). Once water reaches the bed, its impact on ice flow is determined by the hydraulic efficiency of the subglacial hydrological system. Theory developed at alpine glaciers suggests that increasing water pressure results in enhanced ice motion owing to reduced ice-bed contact (Lliboutry, 1968;Bindschadler, 1983) and where sediment is present, enhanced sediment deformation (e.g. Engelhardt et al., 1978;40 Hodge, 1979;Iken and Bindshadler, 1986;Fowler, 1987;Iverson et al., 1999;Bingham et al., 2008). Water pressure at the bed depends on water supply to, storage within and discharge through the subglacial hydrological system (Iken et al., 1983;Kamb et al., 1985;Nienow et al., 1998). Hydraulically efficient drainage can accommodate and evacuate greater meltwater flux with fewer and smaller spikes in basal water pressure, 45 and thus has less impact on ice motion.
In theory, in a steady-state system, water flows from surrounding high pressure distributed regions into lower pressure conduits. Borehole measurements of subglacial water pressure, modelling and ice velocity proxy data (e.g. Hubbard et al., 1995;60 Gordon et al., 1998;Bartholomaus et al., 2008;Werder et al., 2013;Tedstone et al., 2014) suggest, however, that given a sufficiently large and rapid spike in water delivery to a subglacial conduit, the hydraulic gradient can be reversed such that water is forced out and laterally away from the conduit across a zone some metres to several kilometres wide that has been referred to as a Variable Pressure Axis (VPA) (Hubbard, 65 1995). Thus, during significant fluctuations in the delivery of surface water to the bed, conduits interact with their surroundings, supporting the notion of a three-system model of: channelised subglacial drainage connected to the supraglacial hydrological system; distributed subglacial drainage weakly connected to the subglacial channels; and hydraulically isolated subglacial distributed drainage (Hoffman et al., 2016). alongside and associated with the coarse gravelly central ridge (e.g. Cummings et al., 105 2011a;Prowse, 2017). These splays are an order of magnitude wider and more gently sloped than the main ridge (Cummings et al., 2011a). They are proposed to form in both proglacial environments, representing subaqueous outwash fans deposited by sediment laden plumes exiting a subglacial conduit into a proglacial lake (e.g. Powell, 1990;Hoyal et al., 2003;Cummings et al., 2011b), and / or subglacial environments, 110 with sedimentation in subglacial cavities alongside the main esker ridge during periods of conduit over-pressurisation (e.g. Gorrell and Shaw, 1991;Brennand, 1994).

115
Erosional subglacial meltwater channels, or Nye-channels, incised into bedrock or sediment substrate range in size from metres to tens of metres wide (e.g. Sissons, 1961;Glasser and Sambrook Smith, 1999;Piotrowski, 1999) to large tunnel valleys several kilometres in width and tens of kilometres long (e.g. Kehew et al., 2012;van der Vegt et al., 2012;Livingstone and Clark, 2016). Tunnel valleys are observed to 120 occur at various developmental stages from mature and clearly defined to indistinct valleys often associated with hummocky terrain or as a series of aligned depressions (e.g. Kehew et al., 1999;Sjogren et al., 2002). Their formation has been linked to subglacial meltwater erosion at the ice-bed interface (c.f. Ó Cofaigh, 1996;Kehew et al., 2012;van der Vegt et al., 2012) with the assumption that channels transported 125 large volumes of sediment and water. However, their precise mechanism of formation is still debated with the main arguments focussing on i) catastrophic outburst formation with rapid erosion following the release of sub or supraglacially stored water (e.g. Piotrowski, 1994;Cutler et al., 2002;Hooke and Jennings, 2006;Jørgensen and Sandersen, 2006) and ii) gradual steady-state formation with headward erosion of 130 soft-sediments in low water pressure conduits (e.g. Boulton and Hindmarsh, 1987;Mooers, 1989;Praeg, 2003;Boulton et al., 2009).
Here, we use the term meltwater channel to refer to palaeo-evidence of erosional channelised flow left on the ice sheet bed (i.e. the outline of the path the water took) 135 at all scales from N-channels through to tunnel valleys. We use the term conduit to refer to the active channelised flow beneath a contemporary ice mass (i.e. the enclosed (sediment and / or ice walled) pipe carrying water at the ice-bed interface).

140
Detailed mapping in northern Canada and Scandinavia has identified the presence of linear tracks variously termed 'hummock corridors', 'glaciofluvial corridors', 'washed zones' and 'esker corridors', typically a few hundred meters to several kilometres wide and a few kilometres to hundreds of kilometres long (e.g. St-Onge, 1984;Dredge et al., 1985;Rampton, 2000;Utting et al., 2009;Burke et al., 2012;Kerr et al., 2014aKerr et al., , 145 2014bSharpe et al., 2017;Peterson et al., 2017;Peterson and Johnson, 2018;Lewington et al., 2019). These features often contain eskers and hummocks which vary in size, shape and relief (Peterson and Johnson, 2018) as well as 'patches' of glaciofluvial deposits and areas of exposed bedrock. While a subglacial meltwater origin is largely agreed upon, their precise mode of formation is not yet known. These 150 features are collectively termed meltwater tracks hereon in.
Meltwater landforms are typically mapped and interpreted individually (e.g. Clark Brennand, 2000;Storrar et al. 2013;Burke et al., 2015;Livingstone and Clark, 2016;Mäkinen et al., 2017) rather than as a holistic drainage signature (c.f. 155 Storrar and Livingstone, 2017). As such, it is not yet clear whether or how differing expressions of subglacial drainage are interrelated and to what extent variations in drainage and / or background conditions (e.g. bed substrate, geology and local topography) control the preserved geomorphic signature we see today. This study aims to identify and map all discernible evidence of subglacial meltwater drainage 160 across the Keewatin District of northern Canada from the ArcticDEM. We collectively refer to these as meltwater routes (MW Routes ). Producing an integrated map of all visible subglacial meltwater evidence allows us to report on the varying dimensions and geomorphological expressions of these features, to investigate associations between features traditionally treated separately and to explore potential controls on 165 expression and formation.

Study Area
This study focusses on an area approximately 1 million km 2 to the west of Hudson 170 Bay in northern Canada that surrounds the location of the former Keewatin Ice Divide ( Fig. 1) (Lee et al., 1957;McMartin et al., 2004). The area generally exhibits low relief https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. and is underlain by resistant Precambrian bedrock that is either exposed or covered by till ranging from thin and discontinuous to thick and pervasive (e.g. Clark and Walder, 1994).  (Storrar et al., 2013). The Laurentide Ice Sheet extent displayed in the inset is the Last Glacial Maximum (LGM) at 18 14 C ka BP (21.4 cal. ka BP) (Dyke et al., 2003) and the extent of the Precambrian shield is also mapped (Wheeler et al., 1996). (b) Zoomed in location of the study area focussed on the 180 area around the former Keewatin Ice Divide.
Eskers have been widely mapped in northern Canada. Initial mapping was largely undertaken by the Geological Survey of Canada using air photography and field observations (e.g. Aylsworth and Shilts, 1989). This included mapping of 'esker systems' -comprising a series of hummocks or short, flat-topped segments which 195 phase downstream into relatively continuous esker ridges or occasionally beaded eskers -across 1.3 million km 2 of the Keewatin sector of the Laurentide Ice Sheet (LIS) (Aylsworth and Shilts, 1989;Aylsworth et al., 2012). Discontinuous esker ridges are connected to areas of outwash, meltwater channels or belts of bedrock stripped free of drift. More recently, increasing availability of remotely sensed data allowed Storrar 200 et al., (2013) to digitise eskers at an ice sheet scale for the LIS (including the Keewatin sector) using Landsat 7 ETM+ imagery. From this, a secondary dataset was derived by interpolating a straight line between successive aligned esker ridges, creating a continuous pathway, which reflects the location of the major conduits in which the eskers formed (Storrar et al., 2014a). 205

210
ArcticDEM 10 m resolution DEMs (freely available at https://www.pgc.umn.edu/data/arctcidem), generated by applying stereo and autocorrelation techniques to overlapping pairs of high-resolution optical satellite images (Noh and Howat, 2015;Porter et al., 2018), were used in this study to identify and map meltwater landforms. Eskers mapped by Storrar et al. (2013) from 30 m resolution 215 Landsat ETM+ multispectral imagery were used to inform further high-resolution esker mapping from the 10 m resolution DEM. The automatic mapping approach developed in Lewington et al. (2019) was used to create a first pass map of hummock corridors -classified as meltwater tracks here (Appendices Fig. A1) to augment the improved esker map. Together, these were used to create an integrated map of MW Routes by 220 manually mapping centrelines of all visible traces of subglacial meltwater drainage including meltwater tracks, meltwater channels and eskers. Multiple orthogonal hillshades were generated to avoid azimuth bias (Smith and Clark, 2005) and mapping was undertaken at a range of spatial scales to maximise the number of features captured (Chandler et al., 2018). 225 https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.

Classification and morphometry
The MW Routes were used to explore the occurrence and morphology of different types of meltwater landforms. Former ice-margin estimates from Dyke et al. (2003) 230 were used as transects. These transects are spaced approximately 30-40 km apart and in the study area, cover c. 1,000 years of deglaciation between 9.7 and 8.6 ka.
This period encompasses the final stages of deglaciation when the ice sheet was experiencing a strongly negative surface mass balance with associated increasing rates of meltwater production (e.g. Carlson et al., 2008Carlson et al., , 2009. Retreat rates were 235 generally between 100-200 m yr -1 from 13 to 9.5 ka, increasing rapidly between 9.5 and 9 ka to around 400 m yr -1 after which retreat rate decreased briefly before another increase from ~8.5 ka (Dyke et al., 2003).
When a MW Route intersected a transect, an intersection point was added, and 240 the following information recorded: -Landform type (i.e. esker ridge, esker with lateral splay, meltwater track or meltwater channel) -Width of landform (or landforms if an esker ridge was present within a meltwater 245 track, meltwater channel or surrounded by a lateral splay) -Bed substrate and geology (Fulton, 1995;Wheeler et al., 1996).
Spacing between adjacent MW Route centrelines was calculated along each transect with centrelines at the end of each transect and those separated by clear breaks (e.g. 250 due to the coincidence of a lake) discounted. The total length of MW Route centrelines was calculated automatically in a GIS.

255
This study takes a large-scale approach to exploring controls on MW Routes width and expression. While this approach results in a compromise in terms of data resolution available for the surface substrate and geology maps, it also increases statistical confidence in the results due to the larger sample size. Before analysis was https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. undertaken, three test sites were selected from the study area to allow for more 260 detailed mapping and comparisons (Appendices Figs. A2 and A3).
To explore substrate and geological controls on MW Route occurrence, distribution and properties, the overall length of MW Routes overlying each substrate type (Fulton, 1995) and geology (Wheeler et al., 1996) within the three test sites was 265 calculated. The total area of each basal unit within the test sites was also calculated and values converted to percentages. Following this, the percentage area was subtracted from the percentage of MW Routes for each individual substrate and geology type, giving a positive (over-represented) or negative (under-represented) value. Next, MW Routes were split and classified by feature expression (i.e. esker, esker with lateral 270 splay, glaciofluvial material, anabranching sections, meltwater channel and MW Track ).
The above analysis was repeated by feature type to explore whether geomorphological expression is controlled by surface substrate / geology. It is important to note that categorisations along MW Routes were not always independent as the same section was sometimes coded as a meltwater track and an esker with splay 275 / glaciofluvial deposits as often 'positive' features are situated within wider erosional corridors.
It was noted that landform type varies both across adjacent MW Routes and along individual MW Route centrelines. To assess any potential relationship between landform 280 expression and background controls in more detail, individual centrelines were selected and sampled with a higher frequency (5 km intervals). At each sample point the width of the MW Track / meltwater channel, the presence / absence of an esker (and its width if present), surface substrate, bed geology and elevation were recorded.

285
To investigate the relationship between subglacial drainage and bed roughness, the standard deviation of the ArcticDEM (resampled to 50 m) was calculated using a 1 km and 5 km radius circular search window and this was compared to the density of MW Routes. Ice stream locations (Margold et al., 2015) were also qualitatively compared to the distribution of MW Routes . 290 https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.

An integrated drainage signature 295
Mapping all traces of meltwater drainage reveals the ubiquity of former subglacial drainage across the study area (Fig 2). A total of 2,795 MW Routes were mapped over a ~1 million km 2 area with a total length of 51,856 km. The MW Routes exhibit a similar overall pattern to earlier esker maps (e.g. Aylsworth and Shilts, 1989;Storrar et al., 300 2013) radiating out from the former Keewatin Ice Divide. Greater than 90 % of mapped esker ridges in this region are estimated to occur along a MW Route and therefore form part of this wider network of drainage features. In terms of the largescale pattern, there are no obvious trends in MW Route density, width or feature type associated with margin retreat. However, the study area only covers approximately 305 1,000 years, associated with a period of intense meltwater production and rapid retreat.  MW Routes reach a maximum of 340 km in length and 3.3 km in width (Table 1), but are noted to reach up to 760 km when they extend beyond the limits of the study 330 area (Storrar et al., 2014a). Meltwater channels and meltwater tracks are typically an order of magnitude wider (mean width: 900 m) than the eskers which they often contain (mean width: 97 m). MW Routes appear to vary in width around the ice sheet and along individual centrelines but show no clear (temporal) trend from the ice divide towards the margin. Within the study area, adjacent centrelines are spaced on average 335 8 km apart (Table 1). This is at the lower end of the range that has been observed and modelled (e.g. Banerjee and McDonald, 1975;St-Onge, 1984;Shilts et al., 1987;Hebrand and Amark, 1989;Bolduc, 1992;Boulton et al., 2009;Hewitt, 2011). Like variations in width, there appears to be no coherent change in spacing over time within this area (Fig. 3).  This paper extends earlier work which recognises links between eskers and broader traces of subglacial meltwater flow but does not explicitly describe or formally 355 https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. quantify them (e.g. Aylsworth and Shilts, 1989;Storrar et al., 2014a). It is encouraging that despite different datasets and mapping procedures, the overall patterns are similar (Fig. 3). Here we go beyond mapping eskers alone and identify and classify each part of the MW Routes , including those which do not contain an esker. We describe them in detail (section 4.3) and explore the role of background controls in influencing 360 geomorphic expression (section 4.4), providing insight into the drainage system operation (section 5).

365
Landforms along MW Routes exhibit a high degree of geomorphic variability in terms of relief and definition (Fig. 4). Some exhibit strongly negative relief (down to ~30 m below their immediate surroundings) and distinct boundaries, i.e. meltwater channels ( Fig. 4a and 4c), while others exhibit negligible relief and are more difficult to differentiate from their surroundings i.e. meltwater tracks (e.g. Fig. 4d and 4e) or even 370 exhibit positive relief i.e. eskers with lateral splay (e.g. Fig. 4f, 4g and 4h). Channel edges vary from straight to crenulated and from continuous to discontinuous. A variety of forms are found within the meltwater tracks and meltwater channels. These include hummocks of varying size, shape and relief and eskers and associated glaciofluvial material (e.g. Fig. 4e and 4f). In places, till may be entirely eroded revealing patches 375 of bedrock. Eskers display a higher degree of variability along the MW Routes with single, continuous straight-to-sinuous ridges the exception rather than the norm (e.g. Fig. 4a and 4b).
Furthermore, they display a qualitatively similar width range -several hundred meters 410 to ~3 kilometres (Fig. 6) -although, the null hypothesis that the data in each pairing are from the same continuous distribution using the two-sample Kolmogorov-Smirnov test could not be rejected for any pairings (esker, esker with splay, meltwater channel and meltwater track) at the 5 % significance level.

Controls on the width and expression of meltwater landforms
Most subglacial meltwater landforms occur within areas of till (Fig. 7). Meltwater   (Fulton, 1995) and background geology (Wheeler et al., 1996). 'Other' includes marine, 435 lacustrine and glaciofluvial sediments. calculation (1 km and 5 km). Palaeo-ice streams are rare in the Keewatin District region (Stokes and Clark, 2003a,b;Margold et al., 2018), but where they do occur (e.g. Dubawnt Lake and Maguse ice streams), MW Routes are noticeably sparser and more fragmented (Fig. 9). On the bed of the Dubawnt Lake Ice Stream, the MW Routes 445 also exhibit a more dendritic arrangement and extend further towards the ice divide.

455
To explore potential controls governing how meltwater landform expression varies with changing background conditions, detailed measurements of width, feature type and substrate were extracted along individual MW Routes (Fig. 10). Although there is no consistent ratio between esker width and the associated meltwater track / meltwater 460 channel width at sample points across the whole study area, qualitatively there is a positive relationship between the two, with both increasing and decreasing in phase ( Fig. 10). Esker width varies along individual pathways (varying by almost 300 % (Fig.   10b)) and expressions range from a single well-defined ridge to multiple, fragmented and anabranching sections. The MW Route in Fig. 10b suggests that wider sections 465 broadly coincide with higher bed elevations and the presence of more deeply eroded, yet narrower sections with decreasing elevation. In Fig. 10d we observe a large increase in width associated with a rapid increase in elevation, which also coincides with a transition from a strongly negative feature (a meltwater channel) to a positive / https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. depositional feature (esker with lateral splay). This sharp transition is related to the 470 emergence of the MW Route out of the Thelon sedimentary basin.

Discussion
Our new MW Routes map shows that meltwater tracks, meltwater channels and esker splays which flank and connect (in an along-flow direction) esker ridges are a dominant 500 part of the landscape across the former Keewatin sector of the LIS. Mapping complete drainage pathways means we are better able to identify regional meltwater drainage patterns and unravel controls on feature expression.
The large-scale distribution and pattern of meltwater tracks, channels and 505 eskers with lateral splays exhibit several key similarities, including their width, spacing, association with eskers and occurrence within an integrated network characterised by transitions to and from different types along individual MW Routes (Figs. 4-6). Together, this provides strong evidence that these meltwater landforms are varying expressions of the same phenomenon and we therefore collectively group these features with 510 widths in the order of 100s to 1000s of meters and term them meltwater corridors (MW Corridors ) ( Table 2) . This is consistent with previous conceptual work linking meltwater landforms. For example, Sjogren et al., (2002)     https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.
While recognised previously in local case studies (e.g. St-Onge, 1984;Rampton, 535 2000;Utting et al., 2009), we confirm that across this 1 million km 2 area of the former LIS, MW Corridors of varying geomorphic expression are widespread rather than an isolated phenomenon. Within the study area, 87 % of all esker ridges were flanked by a MW Corridor . Esker ridges alone were captured at just 6 % of our sample points while a MW Corridor was recorded at 90 % (7 % of these were lateral splay, 21 % were 540 meltwater channels and 73 % were meltwater tracks) -the remaining 4 % were coded as unclassified. However, we do note that the presence / absence of an esker at the sample point may not be indicative of the entire length of the MW Route as in many cases the esker ridge within a MW Corridor was fragmented. This suggests that the model of Rchannels across the Canadian Shield (e.g. Clark and Walder, 1994) is an 545 oversimplification and the range of landforms observed suggests the presence of various modes of subglacial drainage varying between R-channels entirely cut up into the ice to those incised into the bed with a range in-between. Variations in form and pattern are likely to have been influenced by glaciological conditions, including ice velocity, viscosity, temperature, and thickness, hydraulic potential gradient and water 550 flux, and background conditions such as basal substrate, topography and local scale roughness.

A proposed model of formation: interaction of channelised and distributed drainage elements 555
To interpret palaeo-landforms and reconstruct subglacial meltwater behaviour an understanding of the processes that formed the landforms is needed (i.e. the 'glacial inversion' problem, e.g. Kleman and Borgström, 1996). One approach to understanding glacial processes is through contemporary observations. In this 560 section, we demonstrate how modern observations and modelling of the Variable Pressure Axis (VPA) around a subglacial conduit (e.g. Hubbard et al., 1995) is consistent with the form and distribution of mapped MW Corridors and can explain the range of depositional to erosional signatures observed in the study area.

565
In steady-state conditions, conduits theoretically operate at a lower pressure than the surrounding high-pressure weakly connected system, and therefore will typically draw water in from their surroundings. However, variations in borehole water https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.
pressure measurements observed beneath glaciers in the Alps (e.g. Hubbard et al., 1995;Gordon et al., 1996) and Alaska (e.g. Bartholomaus et al., 2008), ice velocity 570 measurements taken from the Greenland Ice Sheet (e.g. Tedstone et al., 2014) and numerical modelling (e.g. Werder et al., 2013), suggest that large and / or rapid meltwater inputs can cause spikes in conduit water pressure .
This temporarily reverses the hydraulic potential gradient and causes water to flow out of the conduit and into the surrounding weakly-connected drainage system (Fig. 11). 575 This is because the conduit cross-sectional area cannot be expanded rapidly enough by wall melt to accommodate high frequency (e.g. diurnal surface melt) and / or high magnitude (e.g. supraglacial lake drainage) fluctuations in meltwater delivery. As meltwater delivery to the conduit wanes or the conduit cross-sectional area grows at a rate sufficient to accommodate the additional meltwater input, water pressure in the 580 conduit decreases and the conduit will again begin to capture water from its surroundings. Hubbard et al., (1995) note that the pressure perturbation decreases with lateral distance from the conduit like a dissipating wave across the VPA. The distance either side of the conduit influenced by these variations (i.e. the width of the VPA) appears variable over time and space, reaching a maximum of ~140 m in the 585 Alps (Hubbard et al., 1995;Gordon et al., 1996) and modelled to be ~2 km for the western margin of the GrIS (Werder et al., 2013). This tallies well with the spectrum of MW Corridor widths observed in this study, which range from ~100 m to > 3 km.
https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. Here, we propose that MW Corridors are the geomorphic signature of a VPA and 600 represent the interaction between a conduit and the surrounding weakly connected system. Variations in MW Corridor expressions are hypothesised to occur due to variations in discharge and pressure (frequency, magnitude and duration) as well as https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. background glaciological, geological and topographical conditions, which together determine the relative geomorphic activity (i.e. net erosion or deposition). 605 In this model, erosion of corridors with negative relief on the order of metres to 10s of metres depth is likely the result of large magnitude perturbations such as supraglacial or subglacial lake drainages, with the rising limb of the flood event overpressurising the conduit and flooding the entire VPA (Fig. 11a). Flooding of this 610 broader VPA zone is analogous to a narrow sheet flood causing localised flotation or hydraulic lifting (e.g. Brennand, 1994); similar processes have been recorded during jökulhlaups in Iceland (e.g. Russell et al., 2007). Flow across the VPA would enhance erosion by mobilising and entraining unconsolidated sediment in the high velocity flow (e.g. Russell et al., 2006 and references within). As the flood event wanes, the 615 pressure perturbation is reversed causing water and sediment to flow back towards the main conduit where it is rapidly evacuated. Less well defined features that have experienced up to a few metres of erosion (i.e. those with more subdued relief and less well defined boundaries) may be the result of repeated lower magnitude drainage and pressure perturbations (e.g. from diurnal melt or rainfall) forcing water out across 620 the VPA. Instead of completely overwhelming and flooding the system, these events may just fill and expand adjacent cavities (Fig. 11b). Delivery of water to the bed is likely to occur at approximately the same location with lakes and crevasses 'locked' into place by basal topography (e.g. Gudmundsson, 2003;Sergienko, 2013;Ignéczi et al., 2018), resulting in repeated drainage down the same pathways. In this scenario,  Sedimentological investigations suggest that hummocks within meltwater tracks can occur as a result of erosion (e.g. Rampton, 2000) and / or deposition (e.g. Utting et al., 2009) and our proposed model is able to explain each of these. Repeated coupling / uncoupling of the ice-bed interface as water is forced in and out of the VPA 635 (Cowton et al., 2012) and subsequent erosion over seasons or longer may explain the hummocky topography we see today (Fig. 4). Depositional hummocks may be formed https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. by a rapid increase in cross-sectional area associated with a breach in the conduit margin (i.e. it becomes unsealed) and flow across the VPA, facilitating rapid deposition within minor conduits and cavities (e.g. Brennand, 1994). When the conduit returns to 640 being hydraulically isolated from the VPA following a decrease in water flow and pressure, standing bodies of water left within the cavities may also deposit minor fans which build up over time (Brennand, 1994). Alternatively, sediment may be trapped within cavities formed during turbulent sheet flow as water velocities subsequently decrease (Utting et al., 2009). Indeed, it is possible that lateral splays associated with 645 esker ridges also formed due to conduit breaching and subsequent deposition, with outward fining in lateral fans suggesting conduit unsealing and decreasing hydraulic power (Cummings et al., 2011b). Finally, hummocks may be a combination of erosion and deposition. This is not dissimilar to the interpretation of triangular shaped landforms ('murtoos') which are thought to form from subglacial till transported by 650 creep which is then eroded and shaped by subglacial meltwater, and represents high pressure distributed drainage within a 'weakly connected' zone upstream of channelised drainage (Mäkinen et al., 2017;Ojala et al., 2019).

Previous research indicates that esker ridges can be superimposed on hummocks 655
within MW Corridors (e.g. Peterson et al., 2018), but to date, there are no examples of hummocks overlying eskers. Esker ridges do not always sit at the centre of their MW Corridor but instead meander across them, and are recorded alternately as left, central or right at different transect points. This is consistent with eskers being the final depositional imprint of channelised drainage within the large-scale MW Routes network, 660 with formation close to the ice margin (e.g. Hebrand and Amark, 1989;Storrar et al., 2014b;Hewitt and Creyts, 2019;Livingstone et al., in review), while the corridors represent a composite imprint of drainage over a longer period. 665 670 https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.

A brief comment on possible alternative models
While we support the proposed

Exploring potential controls on network patterns and variations in 700 expression of MW Routes
In this section, we explore spatial controls governing the overall pattern of the subglacial hydrologic network, as well as variations in meltwater landform expression (i.e. the balance between erosion / deposition and the resulting geomorphic 705 eskers with lateral splay within wider meltwater channels (Fig. 4H)). This suggests that while spatial controls may be important and exert some control over the relative 'leakiness' of a conduit and therefore the resulting geomorphic signature at any 710 location, there is also a temporal control. This is likely related to short-term variations in the magnitude and rate of subglacial meltwater delivery to the bed and the systems' ability to accommodate it.
We observe a high degree of channelisation across this sector of the ice sheet bed, 715 but this is not uniform, with the densest areas of MW Routes coinciding with the 'roughest' basal topography (Fig. 8). This may be the result of subglacial drainage route fragmentation around bed obstacles, with a greater number of tributaries and broken patterns common in regions of high bed roughness (e.g. Test Site 3). Basal topography also preconditions the large-scale spatial structure of surface drainage 720  and the association between rough areas and dense clusters of MW Routes could be a response to more surface water penetrating to the bed as the result of extensive crevassing. In Greenland crevasses capture a significant amount of surface water -more than moulins or the hydrofracture of surface lakes (Koziol et al., 2017). Surface meltwater inputs are thought to be an important control on the 725 distribution of drainage across the bed (e.g. Gulley et al., 2012;Banwell et al., 2016) and the formation and evolution of subglacial meltwater landforms (e.g. Banerjee and McDonald, 1975;St-Onge, 1984;Hooke and Fastook, 2007;Storrar et al., 2014b;Livingstone et al., 2015;Peterson and Johnson, 2017).

730
Qualitatively, there are markedly fewer MW Routes coinciding with palaeo-ice stream locations -particularly the Dubawnt Lake Ice Stream (Fig. 9). In addition, the network pattern of MW Routes corresponding with the location of the Dubawnt Lake Ice Stream are more dendritic and extend further towards the ice divide. These observations are consistent with Livingstone et al., (2015), who find fewer eskers on palaeo-ice stream 735 beds where modelled subglacial meltwater drainage is greatest. We suggest the scarcity of MW Routes beneath palaeo-ice streams is due to lower ice-surface slopes and hydraulic gradients, which favour distributed rather than channelised drainage (e.g. Kamb, 1987;Bell, 2008). Where channelised drainage does occur beneath palaeo-ice https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. streams, network patterns are typically more dendritic, which may also be the result of 740 shallower hydraulic gradients enabling greater lateral water flow.
Dynamic mass loss via streaming / surging has implications for ice sheet stability (e.g. Bell, 2008;Christianson et al., 2014;Christoffersen et al., 2014). The Keewatin sector of the LIS exhibits a relatively low frequency of ice streams (Margold et al., 745 2015;Stokes et al., 2016). While this may be partially attributed to the resistant bed of the shield (Clayton et al., 1985;Kamb, 1987;Stokes and Clark, 2003a, b), we also suggest that efficient evacuation of meltwater through the dense channelised network that developed in this region during the final stages of deglaciation, as the climate warmed (Storrar et al., 2014b), would have inhibited the development of fast flow and 750 potentially contributed to the shut-down of existing ice streams (Lelandais et al., 2018). This is consistent with modern observations that link decadal-scale ice-flow decelerations with more pervasive and efficient drainage channelisation driven by increased surface meltwater inputs to the bed Tedstone et al., 2014;van de Wal et al., 2015;Davison et al., 2019). We therefore hypothesise that this large-755 scale inverse relationship between drainage channelisation and ice streaming will exist in other palaeo-ice sheet settings. This potential drainage control on ice-sheet activity may also influence the pace of deglaciation; we note slower retreat rates (~230 m yr -1 ) in the northwest of the study area, which coincide with the highest density of MW Routes , compared to much faster retreat rates (~540 m yr -1 ) associated with the 760 sparsest MW Routes . This conclusion is tentative given uncertainty in the regions deglacial chronology (Dyke et al., 2003) and requires further testing.
At a large-scale, there is a general tendency for MW Routes to form on till, which is more easily eroded than bedrock and where geomorphic evidence is likely to be better 765 developed. Eskers are over-represented on harder, more resistant rock (Fig. 7d) where R-channels are more likely to form (Clark and Walder, 1994;Storrar, 2014a), while there is a slight tendency for meltwater channels (i.e. incisional features) to form on the softer, more erodible sedimentary rock (Fig. 7b). Eskers with lateral splay (i.e. depositional features) appear preferentially on till blankets (Fig. 7c) where there is an 770 abundance of sediment that may overwhelm and clog up the conduit (Burke et al., 2015), while isolated esker ridges favour thin till and are under-represented on thick till. Though detailed long-profiles (Fig. 10) hint at local relationships between bed https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.
substrate changes and the resultant landform expression, we caution against the assumption that this is a widespread occurrence rather than an isolated coincidence. 775

Relevance for understanding Greenland's hydrology and associated ice dynamic variations
Western sectors of the contemporary Greenland Ice Sheet are analogous to our 780 study area: both are underlain by resistant Precambrian Shield rocks and both experience(d) rapid retreat and high meltwater production rates. This is also similar to southern Sweden, which lay beneath the palaeo Scandinavian Ice Sheet, where similar geomorphic features to those described here, occur extensively (e.g. Peterson et al., 2017;Peterson and Johnson, 2017). This study therefore has potential 785 implications for our understanding of the impact of subglacial hydrology on overlying ice dynamics and ice flow regulation of past, current and future ice sheets.
The interaction between a subglacial conduit and the surrounding weaklyconnected drainage system (VPA) is believed to be widespread in contemporary 790 glaciological settings (e.g. Hubbard et al., 1995;Gordon et al., 1996;Bartholomaus et al., 2008;Werder et al., 2013;Tedstone et al., 2014) and has been identified as key to understanding ice velocity variations and predicting future ice sheet mass loss (Davison et al., 2019). However, the true extent and influence of VPAs beneath the Greenland Ice Sheet is unknown due to the challenge of observing contemporary 795 subglacial environments. Palaeo-studies, such as this one, offer the potential to reveal new insights into the nature and configuration of the subglacial hydrological system at an ice sheet scale.
Based on our proposed model and the observations within this study we suggest 800 that the VPA is widespread across the Keewatin sector of the LIS, which may be analogous to parts of western Greenland. The drainage footprint is considerably more dense if we consider MW Corridors as well as esker ridges as indicators of areas of the bed which were influenced by subglacial meltwater. The drainage footprint in our study is estimated to cover ~13 % of the bed (using average width and spacing of MW Routes ) 805 but could realistically vary between 5 % (lower quartile width and upper quartile spacing) and 36 % (upper quartile width and lower quartile spacing) if we assume that https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. MW Routes were active at the same time. This represents an area 25 times greater than that of eskers alone which cover ~0.5 % of the bed (using average esker width (100 m) and spacing (18.8 km; Storrar et al., 2014b). The dynamic influence is likely to 810 extend even further beyond the VPA due to lateral stress transfer within the ice (Tedstone et al., 2014). In the longer term, increased channelisation may have additional ice dynamic implications; effectively removing large volumes of meltwater and reducing ice velocity, potentially limiting dynamic mass loss via streaming / surging (e.g. Storrar et al., 2014b). 815 Figure 13. Schematic of the VPA (light blue) beneath a contemporary ice sheet with inputs 820 from the surface and the bed influencing discharge (Q) and therefore the VPA.

825
We identified and mapped all visible traces of subglacial meltwater drainage in the Keewatin sector of the former LIS. We found that wider meltwater features (meltwater tracks, meltwater channels and eskers with lateral splays on the order of 100s to 1000s m) flanking or joining up intervening segments of esker ridges were common. These have previously been termed and described as different features. However, owing to 830 similarities in spacing, morphometry and spatial location (i.e. part of the same https://doi.org/10.5194/tc-2020-10 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.
integrated network), we propose collectively grouping these features under the term MW Corridor (table 2) incorporating the width of the VPA (i.e. the 'weakly connected' drainage system), the drainage footprint in this sector of the ice sheet covered 5-36 % of the bed (on 850 average), which is 25 times greater than previously assumed from esker studies alone, which only account for the central conduit.
Our results suggest that the overall distribution and pattern of drainage is influenced by background topography, with greater relief resulting in denser 855 channelised networks, possibly due to fragmentation of subglacial drainage around basal obstacles and / or enhanced meltwater delivery to the bed through crevasses.
Channelised drainage is relatively rare beneath palaeo ice streams, which may favour distributed drainage configurations due to the lower ice surface slopes and hydraulic gradients. Meltwater drainage may also influence ice dynamics, with the high degree 860 of channelisation observed in the region able to efficiently dewater the bed causing ice-flow deceleration and limiting ice stream activity.
Further research should focus on determining how common this proposed interaction between conduits and the surrounding distributed drainage system is beneath other palaeo and contemporary ice sheets and the controls governing its variability. We hypothesise that where less surface meltwater is delivered to the bed or ice-surface slopes are shallower, for instance, when the LIS was larger and the climate colder, the geomorphic expression will be less extensive and fainter. This is because conduits are less likely to evolve due to lower hydraulic gradients, and their 870 interaction with the surrounding distributed system reduced because of invariant melt supply. Understanding where this interaction and signature occurs will help confirm or refute our proposed model, and develop understanding of how meltwater drainage evolves over long time-scales and influences ice dynamics and mass balance.