A GIS‐based multi‐proxy analysis of the evolution of subglacial dynamics of the Quebec–Labrador ice dome, northeastern Quebec, Canada

The glacial records of the inner‐core regions of the Laurentide Ice Sheet (LIS) document complex yet coherent patterns reflecting ice‐sheet change (e.g. ice‐divide migration), providing unique insights into past glacial conditions. This study develops a conceptual model of subglacial dynamics evolution within a major ice‐dispersal centre of the LIS in northeastern Quebec, Canada using a GIS‐based analysis of the surficial geologic record. Multiple proxies of subglacial conditions (subglacial streamlined landforms, lake density and lake area over thin drift/bedrock) were analysed through grid‐overlay techniques and then classified based on different proxy variables ranging from highly mobile warm‐based to immobile cold‐based conditions. An additional proxy (till blanket) was used to identify areas of thick till deposition, but with few proxies (few lake or landform metrics). Based on local ice‐flow reconstructions, the most ‘relict’ glacial terrain zone (GTZ1) has warm‐based conditions over 66% of its area and is remarkably well preserved, suggesting laterally extensive warm‐based conditions during the oldest identified ice‐flow phase. This relict glacial terrain is partially overprinted by two subsequent ice‐flow phases in spatially restricted zones in the northeast (73% warm‐based), east‐central (41% warm‐based), and northwest (33% warm‐based) of the study area. A zone of more sluggish conditions (only 3% warm‐based) was identified in the highlands at the centre of the study area, characterized by thin till cover, few landforms, yet with large patches of relatively abundant small lakes, indicative of areal scouring. No clear evidence of sustained cold‐based conditions (i.e. high chemical index of alteration values or high 10Be abundances) was found in the study area. These results suggest that warm‐based conditions (active erosion and/or deposition) were uniformly widespread during the earliest ice‐flow phase, later becoming more spatially restricted with broader sluggish ice conditions. These spatially restricted regions of warm‐based subglacial regimes were likely controlled by surrounding and down‐flow ice streaming. © 2020 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd


Introduction
Inner regions of the Laurentide Ice Sheet (LIS) show crosscutting subglacial landforms indicative of multiple ice-flow phases that several studies have interpreted to reflect ice-divide migration caused by internal ice dynamic shifts and/or other processes that take place far from the ice margins or closer to the ice margins during the final stages of deglaciation (Boulton and Clark, 1990;Clark et al., 2000;McMartin and Henderson, 2004;Trommelen et al., 2012;Hodder et al., 2016;Gauthier et al., 2019). Portions of the glacial landscape in these regions thus contain subglacial landforms correlated to old ice-flow phases. In some places, these old features form extensive coherent landscapes that were largely preserved under unmodified 'frozen-bed' and 'ephemeral' patches (Kleman and Glasser, 2007). Their occurrence indicates a shift in thermomechanical conditions from warm-based (flowing ice) to cold-based (frozen at the bed/not flowing) at some point during the last or a preceding glaciation. Despite their fragmentary nature, these patches of relict glacial terrains offer a window into previous ice-sheet configurations and subglacial conditions (Kleman and Glasser, 2007). In addition, the distribution and degree of overprinting of these relict glacial beds by younger ice-flow phases provide key information about how glacial dynamics evolved over time. These inner regions also contain evidence of more sustained cold-based ice, where subglacial landscape features are lacking (Dyke, 1983;Edmund, 1991;Kleman and Stroeven, 1997;De Angelis and Kleman, 2005;Davis et al., 2006). Therefore, inner core regions of ice sheets often form subglacial mosaics of different Glacial Terrain Zones (GTZs) (Trommelen et al., 2012), and investigating these records can improve our understanding of the long-term subglacial conditions and overall ice-sheet dynamics of the LIS (e.g. Gauthier et al., 2019).
One of the largest ice dispersal centres of the LIS occupied northern Quebec and is commonly referred to as the New Quebec dome (Hillarie-Marcel et al., 1981;Vincent, 1989;Roy et al., 2015), or the Quebec-Labrador dome (QLD) Jansson et al., 2003;Occhietti et al., 2011;Stokes et al., 2012). The QLD has been estimated to be 2-3km thick at the Last Glacial Maximum (LGM) (e.g. Ganopolski et al., 2010;Peltier et al., 2015), with ice-flow dispersal patterns spreading across and west of Hudson Bay during various phases (Prest et al., 2000;Trommelen and Ross, 2014), leaving evidence of the QLD's influence on the overall behaviour of the LIS. The QLD is characterized by a complex ice-flow record (Veillette et al., 1999). At the landscape scale (Figure 1), crosscutting flow sets indicate a partially preserved, long temporal record that has been overprinted locally by late glacial features associated with dynamic adjustments in regional ice-flow configuration (Clark et al., 2000;Jansson et al., 2003;Margold et al., 2015Margold et al., , 2018. The landform record of the QLD is characterized by a horseshoe-shaped unconformity (Clark et al., 2000) separating landforms converging towards Ungava Bay from landforms radially oriented away from Ungava Bay ( Figure 1).
The southeastern branch of this horseshoe unconformity was recently the focus of ice-flow reconstruction and deglacial studies (Dubé-Loubert and Roy, 2017;Dubé-Loubert et al., 2018;Rice et al., 2019). The subglacial record of our 7000km 2 study area is thought to reflect ice-divide migration, as well as changes in paleo-ice stream catchment, which together suggest a complex evolution characterized by transient polythermal subglacial conditions.
The goal of this study is to improve our understanding of the net effect of the changing ice-sheet dynamics on the glacial landscape by analysing the subglacial record of an inner core region of the QLD that experienced complex ice-flow history. Our investigation combines various proxies typically used to analyse patterns and intensity of bedrock erosion and/or glacial landform development (e.g. drumlins) related to actively flowing ice. By placing this spatial analysis within the context of the local ice-flow history and related fragmented GTZs, we get insights into the relative strength of subglacial activity related to the older ice-flow landscapes and younger overprinting phases. The results are further assessed using independent proxies of glacial erosion, specifically the chemical index of alteration (CIA) and abundance of 10 Be in till and bedrock, to determine whether there is evidence of weathering and/or inherited exposure in areas where the other proxies indicate limited net glacial imprint on the landscape or where 'relict' glacial landscapes of unknown age are particularly well preserved. These proxies are then used to interpret changes in subglacial ice-flow dynamics and related thermomechanical conditions throughout the established ice-flow history, including an ice-divide migration within the study area. This work builds on our understanding of how the interior regions of ice sheets evolve from a subglacial landscape perspective throughout one or more glaciation(s).

Study area
The study area is located in northcentral Quebec, Canada (Figure 1) and encompasses the border between two large lithotectonic blocks: the Core Zone in the east and the New Quebec Orogen in the west (Figure 2). The Core Zone is a 300km-wide suite of Archean rocks reworked during the collision of the Superior and Nain cratons, which forms the eastern edge of the Canadian Shield ( Figure 2A) that has been identified as a region of mineral potential with recent bedrock mapping projects (e.g. Corrigan et al., 2018). Different bedrock lithologies can have significant controls on landscape evolution (Phillips et al., 2010), glacial erosion rates (Dühnforth et al., 2010), and glacial sediment properties (McClenaghan and Paulen, 2018), therefore considering the lithological variation across the study area is essential. Bedrock within the study FIGURE 1. A digital elevation model of northern Quebec and Labrador with select regional landforms annotated from Fulton (1995) and the horseshoe unconformity of Clark et al. (2000). The approximate locations of the ice divide identified by Dyke and Prest (1987) are also annotated. Several notable ice streams are also identified (Margold et al., 2015;Rice et al., 2020); thick lines represent large ice divides and thinner lines represent smaller divides (see Dyke and Prest, 1987). The current study area is outlined in red. (Inset) The approximate extent of the LIS at 18ka (Dyke, 2004) area is subdivided into three general domains: the Rachel-Laporte of the larger New Quebec Orogen in the west, the George River in the centre, and the Mistinibi-Raude in the east, both part of the larger Core Zone block ( Figure 2A). The Paleoproterozoic Rachel-Laporte domain is comprised of metasedimentary and felsic plutonic rocks (Sanborn- Barrie, 2016;Corrigan et al., 2018). The metasedimentary bedrock of the Rachel-Laporte domain is highly metamorphosed, with well-developed folding structures. Less structural control on the topography is evident in the quartzofeldspathic gneiss lithologies in the north of the Rachel-Laporte domain ( Figure 2A). The Lac Tudor shear zone defines the eastern border of the Rachel-Laporte domain with the western orthogneiss and George River lithologies. The De Pas batholith is a felsic intrusive upland within the George River domain that strikes north/south through the middle of the study area (Sanborn-Barrie, 2016). This resistant lithology forms a highland region that roughly bisects the study area ( Figure 2B). The George River valley within this domain has the lowest elevations within the study area (290m; Figure 2C). The De Pas batholith is flanked on both sides by~12km-wide bands of Archean orthogneiss. East of the De Pas batholith and its associated eastern orthogneiss lies the high-grade magmatic Mistinibi-Raude domain, which is composed of paragneiss, FIGURE 2. (A) Simplified regional bedrock geology of northern Quebec-Labrador (modified from James et al., 2003). (Inset) Simplified local bedrock geology of the study area (modified after Wardle et al., 1997;Sanborn-Barrie, 2016;Corrigan et al., 2018). (B) High-resolution residual magnetic field data for the study area, highlighting bedrock structures (MERN, 2019). (C) Digital elevation model of the study area with major rivers and lakes labelled (data source: USGS). (D) Simplified surficial geology of the study area (modified from Rice et al., 2017a,b). [Colour figure can be viewed at wileyonlinelibrary.com] 3157 GIS-BASED MULTI-PROXY ANALYSIS OF THE EVOLUTION OF QUEBEC-LABRADOR-LIS granitoids, diorites, and mafic/ultramafic plutonic rocks (van der Leeden et al., 1990). The Mistinibi-Raude domain is within the larger Core Zone and marks its eastern boundary.

Surficial geology and ice-flow history
The study area is characterized by bedrock uplands that are glacially sculpted and flanked by discontinuous till veneers, with thicker till deposits in the valleys ( Figure 2D; Rice et al., 2017a,b). Glaciofluvial erosional and depositional features (e.g. meltwater channels, eskers, kames, and kettles) are common within the central and eastern parts of the study area, with abundant meltwater channels cut into the bedrock of the central highland. Eskers occur in the eastern part of the study area and are radially oriented towards the east-northeast. Two glacial lakes briefly inundated the study area during ice margin retreat -Lake Naskaupi in the east and Lake McLean in the west (Ives et al., 1976;Jansson, 2003;Dubé-Loubert et al., 2018). Shoreline and littoral features represent the main evidence for the existence of these short-lived lakes as fine-grained ice distal lake deposits were not observed in the study area.
Within the study area, four ice-flow phases were identified ( Figure 3; Rice et al., 2019). The earliest ice-flow phase (Flow 1) was to the northeast and covered most of the study area; it has been documented in other regions of Quebec and Labrador (Klassen and Thompson, 1993;Veillette et al., 1999). Flow 1 was followed by the development of an ice divide within the eastern portion of the study area, which resulted in a northwest-trending ice-flow phase (Flow 2) in the western sector of the study area. This ice-flow phase was likely caused by a shift in ice dynamics associated with ice streaming occurring in Ungava Bay (Clark et al., 2000;Jansson et al., 2003;Margold et al., 2015). A later westward shift in the dynamics of the Ungava Bay ice streams (UBIS)  led to a westward migration of the dispersal divide, commonly referred to as the Ancestral Labrador ice divide (Vincent, 1989), and a concomitant eastward shift of ice flow. This westward shift is recorded by younger east-trending ice-flow features (Flow 3) across the central portion of the study area ( Figure 3; Rice et al., 2019). Flow 3 likely occurred sometime during deglaciation as the eskers in the study area are parallel to the trend of this ice flow phase. Paleo-ice streams also developed during this phase, such as the Kogaluk River ice stream (Margold et al., 2015;Paulen et al., 2015Paulen et al., , 2017a and the Cabot Lake ice stream Rice et al., 2020), but their remaining footprints are mostly located outside of the study area, with the exception of the northeast corner (Lake Mistinibi region; Figure 2C), where subglacial lineations appear to be part of the converging catchment of the Kogaluk River ice stream. Following Flow 3, there was a minor late deglacial ice-flow phase (Flow 4), with ice flow mainly trending to the west, likely sourced from a smaller localized ice centre located on the central De Pas upland . Evidence of Flow 4 is limited to striations, which suggests it was short-lived and associated with thinner, less erosive ice. In contrast to areas further north (Dubé-Loubert and Roy, 2017), frost-shattered bedrock or regions of oxidized glacial deposits were not identified during fieldwork within the study area .
Previous studies of complex glacial landscapes have identified spatially distinct subglacial zones (Stea and Finck, 2001), glaciodynamic zones (Tremblay et al., 2013), or GTZs (Trommelen et al., 2012) that provide important insights into changing subglacial dynamics across glaciated terrains. The separation of different regions based on their surficial characteristics was conducted southwest of Hudson Bay, where ice-sheet reconfiguration and other factors led to multiple ice-flow shifts and transient subglacial conditions which together resulted in spatially coherent, yet distinct, GTZs (Trommelen et al., 2012;Gauthier et al., 2019). Using these principles, four GTZs were identified in the study area ( Figure 3; Rice et al., 2019).
GTZ1 contains numerous northeast-trending ice-flow indicators associated with Flow 1 (average orientation 52°). Streamlined landforms oriented to the northeast in that area typically consist of large rock drumlins and crag-and-tails; this GTZ is interpreted as a relict flow set preserved under later cold-based conditions based on a crosscutting relationship with indicators associated with younger phases . In the northeast corner of the study area, around the Lake Mistinibi region ( Figure 2C), numerous smaller and more elongated drumlins occur within a slightly oblique orientation relative to Flow 1 features; they trend towards 61°. This area was initially included within GTZ1 because of the northeastern orientation of the landforms , however, we have refined this and these smaller oblique landforms are now correlated to the converging catchment of the Kogaluk paleo-ice stream (#187, Kogaluk River; Margold et al., 2015) located further to the northeast (Figure 1). Therefore, the northeast corner of the study area has been amended to include GTZ3b, as the degree of overprinting by late-stage ice flow appears to be stronger in that area than elsewhere within GTZ1. Elsewhere, however, GTZ1 remains the same as in Rice et al. (2019), as it is clearly associated with the oldest regional northeast ice-flow phase. It is characterized by abundant streamlined landforms and striations of similar northeast orientation with very limited overprinting.
This 'relict' landscape is also locally crosscut at a high angle (~45°) by eskers interpreted to have formed during deglaciation. These eskers could have formed as efficient drainage corridors lowering water pressure, raising effective pressure, and inhibiting ice flow (Boulton et al., 2007), or as time-transgressive short segments close behind a retreating ice mass of limited erosive power (Livingston et al., 2015;Lewington et al., 2020). It is important to note that GTZ1, as delineated within our study area, lies approximately between two late-stage paleo-ice streams: the Kogaluk ice stream (#187) to the northeast (Margold et al., 2015;Paulen et al., 2017a) and the Cabot Lake ice stream to the south Rice et al., 2020). Cold-based inter-ice stream conditions could thus have contributed to preserving the relict (Flow 1) features of GTZ1. GTZ2 is dominated by landforms and striations associated with Flow 2 to the northwest, with some striation evidence from Flow 3. GTZ3a is characterized by smaller, east-trending subglacial lineations (rock drumlins, crag-and-tails, till drumlins) overprinted on the larger northeast-trending features associated with Flow 1, as well as by short segments of generally east-trending eskers. GTZ3a is thus defined by relatively strong overprinting of Flow 1 landscape by Flow 3 features, although a palimpsest Flow 1 landscape is still evident in the outcrop and landform records. GTZ4 is the largest zone within the study area and has multiple outcrops with crosscutting striations, double stoss-and-lee outcrops, and few landforms. The internal features of these GTZs show variability that has yet to be characterized. Determining the variation of subglacial proxies within each GTZ, and the contrasting assemblages between them, will offer insights into how subglacial dynamics varied during the different relative ice-flow phases.

Timing of ice-flow phases
The age of the ice-flow phases identified within the study area is unknown. The use of radiocarbon dating can only be useful for deglacial phases associated with ice margin retreat (e.g. Dalton et al., 2020), which means it is of limited use in the study area. Flow 1 ( Figure 4A) likely occurred at a time when the LIS was thick over the study area, as the ice-flow features associated with this ice-flow phase are remarkably parallel and appear to be unaffected by topography. In addition, ice is a good insulator and thick ice can thus trap geothermal heat and promote warm-based conditions in the inner core regions of a thick ice sheet (e.g. Marshall and Clark, 2002). Nonetheless, ice-sheet models indicate low probabilities of warm-based conditions for the study area at any time during the last glaciation but it is during glacial stadials, such as the LGM, that the extent of cold-based zones is at a minimum (e.g. Tarasov and Peltier, 2007).
Recent work by Stokes et al. (2016) and Margold et al. (2018) brought insights into the timing of paleo-ice streams across the LIS. This work has indicated that the earliest paleo-ice stream (#17, Ungava Bay fans 1; Margold et al., 2018) to penetrate the QLD was the Ungava Bay fans 1 that formed~11.5calka BP ( Figure 4B). The Ungava Bay fans 1 correlates to our Flow 2 and indicates that paleo-ice streams propagated a significant distance up ice, close to the hypothesized ice-divide location. Similar paleo-ice stream propagation has been documented in the Keewatin sector during the LGM (Hodder et al., 2016) and in northern Hudson Bay (Ross et al., 2011). Following the establishment of the Ungava Bay fan 1 (#17, Kogaluk River; Margold et al., 2018), a shift in ice-sheet dynamics resulted in the westward migration of the position of the Ungava Bay fan 2 catchment (#16, Ungava Bay fans 2; Margold et al., 2018), creating an entirely different paleo-ice stream flow set that crosscuts the older one. This second phase is suggested to have occurred around 10.1calka BP ( Figure 4C; Margold et al., 2018). This westward shift of paleo-ice streams in Ungava Bay is consistent with the westward migration of the ice divide and the development of east-trending ice flow (Flow 3 of this study) in the study area, and thus correlates to transient events.
Following the paleo-ice streaming events in Ungava Bay, multiple ice streams developed on the eastern margin of the QLD at around 8.9calka BP (#186, Happy Valley-Goose Bay; Figure 4C; Margold et al., 2018). Paleo-ice streams not identified by Margold et al. (2018) but that have been identified through more recent mapping are along the same longitude (Smallwood ice stream, Paulen et al., 2017b; Cabot Lake ice stream, Paulen et al., 2019;Rice et al., 2020). These paleo-ice streams are of a similar orientation to Margold et al.'s (2018) ice streams #186 and #187, suggesting that they also formed during this time. The relatively short time between these paleo-ice streams (~2.6calka BP between Ungava Bay fan 1 and Happy Valley-Goose Bay) suggests relatively rapid changes to the overall organization of the QLD, prior to deglaciation ( Figure 4D).

Methods
For this investigation, multiple analytical methods were collectively used to analyse and map, in an integrated way, the imprint left by basal ice/sediment movement across the study area. Additionally, we interpret warm-based conditions to represent active ice/bed mobility associated with a combination of erosion, till production, and till deposition, whereas cold-based conditions indicate an immobile ice/bed interface with minimal erosion, till production, and deposition. Specifically, we considered features associated with subglacial erosion (e.g. bedrock depressions/lakes, roches moutonnées, rock drumlins/drumlinoids), as well as features that are related to basal ice flow but may reflect more complex erosion/deposition interplay (e.g. drumlins on till blanket). Using a grid overlay, these proxies are combined into one index to classify the study area along a scale from warm-based to cold-based conditions, based on the relative differences in the proxies. Additionally, the CIA from collected till samples and cosmogenic nuclide abundances from bedrock and till samples was examined to further constrain the interpretation. Each of these analytical methods is explained in more detail below.
Surficial signature of the study area Net subglacial erosion of actively flowing basal ice has been inferred from the spatial density of bedrock-controlled lakes across the Canadian Shield (Sugden, 1978;Andrews et al., 1985), on Baffin Island (Briner et al., 2008;Ebert et al., 2015), and over Iceland (Principato and Johnson, 2009). This approach analyses the abundance of lakes within terrains characterized by thin discontinuous veneer and bedrock outcrops, and thus provides insights into variable spatial erosion from basal sliding across a hard bed. The abundance of striations and other outcrop-scale glacial erosional indicators  supports the glacial erosion formation of lake-filled depressions within these terrains. In these regions, spatial lake density can be used to assess the degree and pattern of spatial erosion by abrasion and quarrying from basal sliding across a hard bed. However, within our study area, surficial maps indicate that bedrock and till veneer covers only about 46% of the study area ( Figure 2D). The remaining 54% of the study area is covered with till blanket and other surficial sediments that mask the underlying scoured bedrock surface. Lakes at the surface in these areas of thicker sediments occur in depressions that formed through a variety of processes (e.g. kettle lakes), and are thus not reliable indicators of subglacial erosion related to basal flowing ice. The till in the study area has many characteristics associated with subglacial traction tills (Rice et al., 2017a, b), and thus till blankets (till >2m thickness) are interpreted as representative of areas of actively flowing warm-based ice, even though it may have been deposited later through basal meltout under stagnant ice. Proxies other than bedrock erosional proxies therefore need to be used across these areas of thicker subglacial sediments. Cells containing >50% till blanket, but lacking landforms, were thus automatically categorized as till blanket and used as a proxy for warm-based conditions during cell classification (see below).
Elongated subglacial landforms have been used as a proxy for basal ice-flow velocity through measurement of landform elongation (Stokes and Clark, 2002;Bradwell et al., 2007;Stokes et al., 2007Stokes et al., , 2013Clark et al., 2009;King et al., 2009;Ross et al., 2009Ross et al., , 2011, and landform spatial density (Hart and Smith, 1997;Stokes and Clark, 2001;Briner, 2007;Principato and Johnson, 2009;Dowdswell et al., 2010). In our study area, elongated features of glacial origin occur over bedrock, till veneer, as well as parts of till blanket areas. Therefore, an analysis that considers not only bedrock-controlled lakes, but also the elongated bedrock and bedrock-cored features (e.g. whalebacks and rock drumlins), as well as till-cored drumlinoid landforms and till blanket areas (28% of study area), is likely to provide more complete and thus more useful information about subglacial conditions and overall dynamics across a larger proportion of the study area.
To conduct a uniform analysis across the study area, 2.5×2.5 km cells were created in ArcGIS (v10.3.1) to create a grid overlay that would allow for a comparison between proxies. This grid size has proven effective in other studies with similar geospatial resolution (e.g. Ebert, 2015). It is important to note that about 26% of the study area needed to be excluded from our analysis because it dominantly consists of proglacial and postglacial surficial sediments (Rice et al., 2017a,b) that do not reflect subglacial processes. Grid cells containing >50% of proglacial and postglacial surficial sediments were 3160 J. M. RICE ET AL.
thus masked and not considered further in the analysis. The remaining cells were then categorized based on lake spatial density, lake area, streamlined landform spatial density, streamlined landform elongation, and till blanket. The ggplot2 package (Wickham, 2009) within the R statistical environment (R Core Team, 2015)which is commonly used for geoscience data analysis (Grunsky, 2002) was also used to plot results for further analysis and interpretation of the data. The following subsections describe in more detail how these different proxies of subglacial dynamics were mapped and analysed.

GIS analysis
Lake analysis Waterbodies data (1:50000 scale) were obtained from the Canadian Geogratis website (www.geogratis.gc.ca). Lakes within identified masked regions were removed from the analysis, as their formation could not be associated with subglacial scouring/bed mobility (Sugden, 1978;Andrews et al., 1985;Ebert, 2015). Within the~7000km 2 study area, there are 27 169 lakes; 6738 lakes (24.8%) were within masked units and not used, leaving 20431 lakes for this analysis. Lake spatial density was measured by first computing the centroid of each lake, using the point density function within ArcGIS. Then the number of centroids per grid cell was calculated to derive a lake spatial density per grid cell. This technique has been used in other glaciated regions (Principato and Johnson, 2009;Ebert, 2015). A potential problem with this technique is that the results are sensitive to the extent of the lake relative to the grid cell size. Within the study area, there are numerous small lakes, as well as a few large lakes ( Figures 5A and B). The mean lake surface area is 0.18km 2 , but 16 large lakes exceed the regular cell size of 6.25km 2 and have a combined surface area of 743.6km 2 . The point density function provides only a single centroid for a lake that covers several cells. To address this, a second technique was used, which consists of calculating the area of each cell that is covered by a lake polygon. With this technique, numerous small lakes or few larger lakes will both give a high value. However, the two techniques were retained for this analysis, instead of choosing just one, because they may reveal patterns associated with different types of subglacial bedrock erosion, or different bedrock lithology and structural effects. For instance, characteristics of the underlying bedrock (i.e. lithology and/or structure) and irregular depth of preglacial weathering profiles have modulated the style of bedrock erosion and final shape and size of lakes in regions of areal scouring resulting from mobile basal ice/sediment (e.g. Krabbendam and Bradwell, 2014). During the onset of deglaciation (e.g. Flow 2), low effective pressure and related fast-flowing basal ice (warm-based conditions) favoured the development of longer lee-side cavities down-ice of bedrock obstructions. Frequent water pressure fluctuations within these cavities can lead to large stress variations in the down-ice edge of the bedrock obstruction, creating conditions that promote quarrying and greater subglacial erosion, especially in the presence of pre-existing subvertical fractures in bedrock (Iverson, 2012;Alley et al., 2019). Quarrying can excavate greater volumes of bedrock than abrasion, but requires a relatively rough bedrock surface to be effective (Cook et al., 2020). The up-ice side of the bedrock obstacles is dominated by abrasion due to the sliding and pressing of debris-rich ice along the stoss side. Therefore, the two processes tend to operate together in areas of irregular hard, but jointed bedrock. This can be preserved and form an irregular bedrock landscape with numerous small rock knobs and depressions, which have since ponded with water, creating small lakes; a characteristic of areal scouring landscapes (Sugden, 1978;Andrews et al., 1985). However, in the presence of a more spatially focused network of steeply dipping faults and fractures (weathered bedrock), with some extensive lineaments oriented parallel (anisotropic bedrock) to ice flow, these glacial erosion processes of abrasion and quarrying can be more effective at excavating large volumes of rock along specific narrow corridors. Under the right glaciological conditions, these areas will thus be affected by higher overall erosion rates, which will produce more elongated and deeper bedrock depressions. Therefore, by analysing both lake density and lake surface areas, we can better capture the spatial distribution of these different situations (i.e. numerous small lakes or few large and elongated lakes oriented with paleo-ice flow). This should give insights into the overall style and intensity of bedrock erosion, as opposed to focusing on exclusively mapping the spatially well-distributed abrasion-quarrying process of the areal scouring landscape. Additionally, erosion rates from modern glacial environments are positively correlated to sliding velocities, indicating that faster ice flow is more erosive and providing further evidence that landscape features indicative of high erosion rates correlate to warm-based highly dynamic ice (Cook et al., 2020).
The lake spatial density and the lake area results were grouped into five classes using Jenks natural breaks. The average elevation for each cell was also calculated and plotted against the lake area within the same cell. This was done specifically to evaluate whether the lakes are organized in a pattern consistent with or independent of topographic controls, which has been shown to have a correlation in other ice-divide regions (e.g. Ebert, 2015).
Landform analysis Glacial landform data were obtained from the Geological Survey of Canada surficial geology maps (1:100000 scale; Rice et al., 2017a,b) that were produced by aerial photograph interpretation and contain detailed identification of glacial landforms. To supplement the surficial geology maps, Shuttle Radar Topography Mission (SRTM) DEM (30m vertical and horizontal resolution; United States Geological Survey, 2014) data were coupled with Landsat 7 Enhanced Thematic Mapper+ satellite imagery (Panchromatic band 8) with a resolution of 15m, which was obtained from the United States Geological Survey's Earth Explorer website (earthexplorer.usgs.gov). The DEM and satellite imagery were used to confirm the correct location and extent of the landforms. For this analysis, all types of elongated oriented subglacial bedforms were analysed (Figures 5C-E).
Subglacial landforms have been identified under active warm-based paleo-ice streams resulting from rapid rates of erosion and deposition (e.g. King et al., 2009King et al., , 2016. Although some debate regarding the formational process of subglacial landforms persists, they appear to be the result of a complex interplay of several subglacial processes associated with subglacial conditions (Hindmarsh et al., 1989;Hindmarsh, 1999;Fowler, 2009Fowler, , 2010Stokes et al., 2013;Fowler and Chapwanya, 2014). Therefore, subglacial landform distribution and elongation data offer significant insights into subglacial conditions that would have persisted at least during their formation.
Landform density was calculated using the same regular 2.5 km raster grid used for lake analysis, to allow for an equal comparison between the methods. The represented values correspond to the number of landforms per cell (i.e. landforms per 6.25km 2 ). To obtain landform elongation measurements, polygons first needed to be created. Polygons were outlined manually using a multi-directional hillshade with a vertical exaggeration of 4× from the source DEM. Although automated estimation of glacial landforms' length, width, and height has been used successfully with high-resolution LiDAR data (e.g. Yu et al., 2015;Sookhan et al., 2016Sookhan et al., , 2018, the lower resolution of any DEM currently available for the study area made the automation process ineffective. Jorge and Brennand (2017) performed an analysis of five different methodologies for estimating landform elongation and concluded that drumlin elongation was best estimated using the minimum bounding rectangle. Therefore, ArcGIS (v10.3.1)'s 'minimum bounding geometry' tool was used to fit each polygon with a minimum bounding rectangle to calculate the elongation of the landform (elongation = length of minimum bounding rectangle divided by width of minimum bounding rectangle). The average elongation value was then calculated for each cell of the same raster grid used for landform density. The landform spatial density and elongation results were grouped into five classes using Jenks natural breaks to allow for a similar comparison to lake spatial density distribution.

Classification of cells
Results from all the proxies described above (i.e. lake spatial density and surface area, landform spatial density and elongation, and till blanket) were then combined and interpreted in terms of subglacial regime using a numerical scale representing warm-to cold-based conditions ( Figure 6). Due to the complexity of the surficial signature (i.e. mix of bedrock, till veneer, and till blanket), a weighted index for classifying cells was impractical. Additionally, as each proxy could indicate warm-based conditions (i.e. high lake area will have no landforms, but still be indicative of warm-based conditions), each proxy was given equal weight during cell classification [i.e. highest cells from each classification (red cells in Figures 7-10) were placed on top of second highest cells (orange), which overlaid the third highest (yellow), fourth highest (green), and lowest (blue) cells]. These layers were merged, keeping this hierarchy ( Figure 6). Each proxy was also evaluated independently (Figures 7-10) in order to assess if a single proxy had the greatest effect on cell classification. Because both coldand warm-based conditions are effectively endmembers of the subglacial thermomechanical spectrum (e.g. Tremblay et al., 2013), with complex transitional zones between them, transitional classes were created to represent the subglacial conditions between these two endmember environments (Figure 6). Therefore, five categories were used for this classification, ranging from warm-based (4) to cold-based (0). 'Warmbased cells' (values of 4) indicate cells with high proxy values (i.e. high density of landforms or lakes). Cells with moderate proxy values (values of 3) are classified as 'moderately warmbased'; values of 2 are classified as 'warm-based undefined' and indicate cells with till blanket but otherwise low proxy surficial values, requiring warm-based conditions for till production and deposition, but low dynamics following till deposition, which did not create streamlined subglacial 3162 J. M. RICE ET AL. landforms on the till surface. 'Sluggish ice cells' (values of 1) are cells with low proxy values closest to cold-based conditions. The last category, 'cold-based' (values of 0) was added to include cells with no proxies and no till blanket, indicative of sustained cold-based conditions. Endmember cold-based zones characterized by weathered nonglacial regolith were not observed during fieldwork, but portions of the landscape could still be classified as such in this GIS analysis. In addition, landscapes typically attributed to long-term cold-based conditions were documented north of the study area by Dubé-Loubert and Roy (2017).

Till weathering and cosmogenic exposure
In principle, cells with higher values (4 and 3) should be characterized by relatively freshly exposed bedrock and by till containing similarly freshly exposed material (i.e. not containing material eroded from interglacial exposed surfaces). In contrast, lower-value cells (2 and 1) would be expected to show a comparable, higher degree of surface weathering and cosmogenic exposure inheritance (from previous interglaciations) due to lower overall erosion of these zones during the last glaciation (e.g. Refsnider and Miller, 2010). Within inner-core regions, there is also the additional possibility that relict glacial sediments and surfaces from a previous glaciation were preserved during the last glaciation under a sustained 'frozen patch' (Kleman and Glasser, 2007). For example, areas of relict Flow 1 glacial landscape (GTZ1), with minimal overprint by younger phases could, in theory, have a higher weathering signature and cosmogenic exposure inheritance if Flow 1 was produced during a previous glaciation (Veillette et al., 1999), as opposed to the most recent glacial event. In summary, an area of limited net glacial erosion would have a higher degree of weathering and higher cosmogenic isotope abundance than the warm-based zones of higher erosion (Staiger et al., 2006;Refsnider and Miller, 2010;Corbett et al., 2016), with the possible exception of zones corresponding to an old glacial bed or 'frozen patch' that would have been exposed to cosmic radiation and weathering during the last interglacial. We thus examine bedrock and till across the different zones using a weathering index and cosmogenic inheritance as proxies for the erosional vigour exerted on the landscape during the last glaciation, and to test for the possibility that 'relict' zones (e.g. GTZ1) have higher weathering and exposure inheritance than expected based on their glacial landform record (i.e. they are preserved from a previous glaciation).
Chemical index of alteration The CIA of a till/regolith matrix has been used on Baffin Island to characterize areas affected by sustained cold-based conditions (Refsnider and Miller, 2010). When employed in concert with other proxies, regions of high CIA values (>70) on Baffin Island correlate to regions of low lake spatial density and high cosmogenic exposure inheritance (elevated 10 Be abundances beyond postglacial exposure) in regolith, and these were interpreted as zones of low erosion or inherited CIA values associated with sustained cold-based conditions (Refsnider and Miller, 2010). In contrast, low CIA values of till similar to the local fresh bedrock reflect widespread glacial erosion of fresh bedrock and deposition of unweathered glacial debris during the last glaciation. In essence, the CIA values calculated from till matrix geochemistry reflect the glacial erosion during the last glaciation and, as previously stated, the residence time of till on the landscape. High CIA values thus indicate immature till/regolith produced from old weathered bedrock or a relict till surface from a previous glaciation that was preserved under a 'frozen patch' throughout the last glaciation.
The CIA is a measure of the proportion of primary minerals to chemically weathered minerals, as defined by Nesbitt and Young (1982): Minerals used in this calculation have been reported according to their molar proportions (weight%/molecular mass). CaO* is assumed to be derived from silicate-bearing minerals. It is important to note, however, that the CIA formula is based on the removal of Na and Ca during the chemical breakdown of feldspar from felsic rocks to clay minerals, which results in the enrichment of Al and the depletion of Ca, Na, and K (Fedo et al., 1995). Fresh felsic rocks typically have CIA ranging from 40 to 55, whereas weathered felsic rocks and derived sediments have higher index values >70 (Nesbitt and Young, 1982). For this study, we are confident that the CaO* is derived from silicate-bearing minerals, as there are no known carbonate lithologies in the region (Sanborn-Barrie, 2016). However, the southern Rachel-Laporte domain is characterized by a medium-grade metamorphic quartz-biotite metasedimentary sequence with small bands of amphibolite and ultramafic rocks containing fewer feldspars (0-40%) (Sanborn-Barrie, 2016). The metasedimentary rocks of the Rachel-Laporte domain within our study area could thus exhibit higher 'fresh rock' CIA values or greater CIA variability, which would limit the use of the CIA index in that portion of the study area (e.g. Refsnider and Miller, 2010), or would require a higher threshold for fresh vs. weathered bedrock CIA in that area. However, the bedrock geology underlying the remainder of the study area dominantly contains felsic rocks that have an average fresh bedrock CIA of 51 (Énergie et Ressources naturelles Québec Report, 2010).
CIA calculations were conducted on till-matrix geochemistry results from 74 till samples collected across the study area from hand-dug pits (depths ranging from 0.2 to 0.7m) into active frost boils or the naturally developed soil horizon to collect unoxidized till. At each sample site, a~3kg sample was collected for geochemical analysis (for detailed collection and analytical methodology, see Rice et al., 2017c). Standard reference material (CANMET standard Till 4, n =12) and silica blanks (n =13) were added to the collected till samples to assess the accuracy of the analysis and possible cross-contamination between samples, following procedures of McClenaghan et al. (2013). The <0.063mm fraction was submitted to Bureau Veritas Commodities Canada Limited, Vancouver for analysis. Samples underwent lithium metaborate/tetraborate fusion and nitric acid total digestion followed by ICP-MS (inductively coupled plasma mass spectrometry) determination (BV LF200 package) on 0.2g aliquots. Quality assurance and quality control (QA/QC) procedures outlined by Piercey (2014) were applied to the geochemical data before CIA calculations. Specifically, precision was evaluated using the relative standard deviation method (Jenner, 1996) for reference material using the major oxides used in the CIA formula; a value of 0.47% was obtained, which is 'excellent'. The precision of the analytical method was also evaluated by plotting results of duplicates on scatterplots and Thompson-Howarth plots, which yielded results with better than 5% precision. QA/QC methods yielded acceptable quality data, confirming that results are suitable for further analysis and interpretation. The full results of the QA/QC can be found in Supplementary Appendix S1. Detailed CIA calculations are reported in Supplementary Table T1.
Cosmogenic 10 Be abundances 10 Be abundances have been successfully used as a proxy for estimating subglacial erosion on Baffin Island Briner et al., 2014), in the Cordillera (Stroeven et al., 2010), Fennoscandia (Fable et al., 2002;Linge et al., 2006), and Greenland (Beel et al., 2016;Corbett et al., 2016). 10 Be accumulates within quartz at the surface of the Earth during periods of exposure to cosmic radiation, primarily through a spallation reaction with oxygen within the quartz (Lal, 1991;Gosse and Phillips, 2001). The rate at which 10 Be accumulates (the production rate) largely depends on latitude, elevation, shielding from cosmic rays, and radioactive decay of 10 Be (Dunai, 2010). During glaciation, subglacial bedrock surfaces are shielded from cosmic radiation and are exposed during interglacial or nonglacial periods. Therefore, without any type of shielding (e.g. ice, snow, or sediment), shallow (<2m) quartz grains within a glaciated outcropping bedrock surface are exposed to cosmic radiation and 10 Be is produced. If glacial erosion were limited throughout the last glaciation and thus did not remove preglacial 10 Be, the bedrock surface would contain an amount of 10 Be inherited from a previous interglacial in addition to the recent postglacial accumulation. Therefore, if all other production factors remain constant, the only major effect on 10 Be abundance will be the erosion of the bedrock surface during glaciation, which provides a proxy to assess the relative vigour of glacial erosion, and thus the net effect of the subglacial thermomechanical regime throughout glaciation.
Available 10 Be concentrations from bedrock and glacial erratics were previously used to constrain the timing of ice-margin retreat across the study area along a sampling transect perpendicular to inferred ice-margin retreat . The sample transects also crossed the identified GTZs, which we hypothesized to have discernible differences in 10 Be abundance (e.g. Briner et al., 2006Briner et al., , 2008Staiger et al., 2006). For our analysis herein, 10 Be abundances from bedrock were used in combination with 10  till samples collected within a few metres of the bedrock samples. The analysis of bedrock and till pairs is useful in this case, as the till has been sourced from a wider area than is represented by a bedrock outcrop; therefore, the bedrock gives local information whereas the till should provide insights into inheritance from the surrounding region (Staiger et al., 2006). Areas that have experienced less overall erosion and whose till is immature (short transport) should thus have a high inheritance. Areas of local low bedrock erosion but with relatively long transport till that was deposited by meltout should show more contrast in 10 Be abundance between bedrock and till samples. Finally, the hypothetical situation of a relict glacial bed from a previous glaciation (e.g. possibly GTZ1) should yield higher 10 Be abundances than expected from a glacially scoured landscape.
The till samples were collected from hand-dug pits that exposed the entire till profile from the surface to bedrock to a maximum depth of 0.35m; we sampled the entire profile as a single bulk sample. This method of sample collection ensures homogenous mixing of quartz, given that cryoturbation affects the surficial till in the region. The eight bedrock samples were collected from windswept bedrock outcrops that had high quartz content (>35%) and were free from present-day shielding such as boulders, surrounding topography, or vegetation. 10 Be abundances were normalized based on production rates and scaling [calculated using Cronus online calculator (V.3); https://hess.ess.washington.edu/], to limit any bias in elevation difference between the samples. Elevation values were corrected for isostatic rebound using data from ICE-6G_C database (Argus et al., 2014;Peltier et al., 2015). Details of the methodology used for 10 Be abundance determination are described in Supplementary Appendix S2.

GIS analysis
Lake spatial density and lake area Lake spatial density based on the centroid method shows a low density of lakes in the central and western part of the study area, with a higher density of lakes in the east and northwest ( Figure 7A). The highest density of lakes occurs within the most frequent elevations (between~400 and 500ma.s.l.; Figure 7B), but moderate-to-high densities at less common elevations also occur, particularly at lower elevations within the Mistinibi-Raude and also at higher elevations within the De Pas batholith ( Figure 7C). Lake area coverage per cell shows a different spatial distribution than the lake spatial density map ( Figure 8A). The cells with the greatest lake area are most common at the most frequent elevations (between 400 and 500m; Figure 8B). In addition, lake area is greatest in the northwest (quartzofeldspathic gneiss-Unit 3), in the northeast, and the southeast (paragneiss-Unit 7; Diorite/Gabbro-Unit 8; Figure 8C).
Landform spatial density and elongation A total of 683 streamlined glacial landforms were identified and outlined in the study area (Table 1; Figure 9A). Cells containing only a single landform are the most common ( Figure 9B), likely due to landform size and spacing relative to cell size. However, similar patterns were obtained using other grids (e.g. cell area of 1and 25km 2 ). The maximum number of landforms within a cell is 12. Few landforms occur in the southwest and west-central part of the study area. Conversely, the highest density of landforms occurs in the northwest (quartzofeldspathic gneiss-Unit 3), east-central (orthogneiss-Unit 5), and northeast (paragneiss-Unit 7) part of the study area ( Figure 9C). Landform elongation is, by definition, limited to cells containing landforms ( Figure 10A). Landforms with elongation ratios between 2.0 and 6.0 are most common ( Figure 10B), with the highest elongation ratios in the northeast lithologies ( Figure 10C). However, plotting of the length and width of the landforms against the elongation value shows no clear difference between the six newly (see discussion) defined GTZs (Figures 10D-F).

Classification of cells
A total of 1300 cells were created as an overlay for the assessment of the outlined proxies with 155 (12%) cells having been masked. This left a total of 1145 unmasked cells that were analysed based on all the GIS-based proxies and surficial maps. From this overlay and the combined results of the proxies, 311 (27%) cells were classified as warm-based, 143 (13%) as moderately warm-based, 474 (41%) as warm-based undefined, and 217 (19%) as sluggish ice (Figure 11). No cells were classified as sustained cold-based (i.e. no till blanket, few lakes, no landforms).

Additional constraints
Chemical index of alteration The CIA results are limited to the 74 cells (6.5% of the total grid) in which till samples were collected for geochemical analysis ( Figure 12A). Values range from 46.4 to 64.7, with an average CIA of 50.2 (σ = 3.49). The highest CIA values are from samples in the southwestern part of the study area within the Rachel-Laporte domain (Bedrock Unit 1; Figure 2), where they range from 49.1 to 64.7 ( Figure 12B) with an average CIA value of 56.1. A scatterplot of CIA values and elevation shows no clear correlation but does show clustering of samples from the Rachel-Laporte domain with the highest CIA values, even when compared to samples from similar elevations ( Figure 12C). 10 Be abundances Table 2 shows 10 Be abundances in bedrock (samples with B suffix) and till samples (samples with T suffix), along with elevation corrections and normalized values. The normalized values are used during the discussion of 10 Be abundances. 10 Be abundances fluctuate across the east to west transect of collected till and bedrock samples ( Figure 11A). 10 Be values from bedrock samples are lower in the eastern and western ends of the transect (4.0×10 4 to 6.1×10 4 atomsg À1 ), except for one sample (15-PTA-081B) which has a relatively elevated concentration of 10 Be (9.1×10 4 atomsg À1 ). Bedrock samples collected from the centre of the study area have the highest abundances of 10 Be (7.4×10 4 to 1.9×10 5 atomsg À1 ; Figure 11B). The abundance of 10 Be in till has similar spatial patterns to bedrock results, with lower abundances in the east and west (3.8×10 4 to 7.9×10 4 atomsg À1 ) and higher abundances in the central part of the study area (8.2 to 9.0×10 4 atomsg À1 ). When these abundances are plotted against the longitudinal position from which the sample was collected a general arc in abundance can be seen, with higher abundances in the middle and lower abundances in the east and west ( Figure 11A). Although there is some overlap between the centre, eastern, and western portions, this overall pattern can still be observed ( Figure 11B).

Preservation of old ice-flow indicators
Preserved outcrop-scale ice-flow indicators were identified at 15 locations within the study area ( Figure 11A). Outcrops with multiple crosscutting older ice-flow indicators are limited to the central part of the study area. These outcrops are located mainly within sluggish ice cells (7 of 15), over the most elevated portion of the region (Figure 2). More uniform ice-flow indicators of consistent direction and relative age are observed in the eastern and western parts of the study area. Flow 4 (green striae, Figure 11) was not considered in this assessment due to its lack of inferred erosional vigour within the study area.

Discussion
Variation of proxies across the study area The combination of proxies indicative of warm-based subglacial conditionssuch as streamlined landform density and elongation ratio, lake area and lake density, and surficial sediment classification, using a grid-cell analysisshows great variation in the strength of the subglacial imprint across the whole study area (Figure 11). There is also variation within each individual GTZ, associated with a specific ice-flow phase (e.g. GTZ1 or GTZ2) or unique crosscutting or overprinting (GTZ3a and GTZ3b). Areas of higher landform density tend to cluster together in patches that occur in the centre of GTZ2, GTZ3a, and GTZ3b, indicating that conditions for landform formation (warm-based) are optimal in the middle of the GTZs, with less favourable conditions towards the exterior, possibly due to slower-moving ice towards the lateral margins of the confined ice flow where the margins impose side drag (Stokes and Clark, 2002). Interestingly, GTZ1 has a more evenly distributed high density of landforms with respect to other GTZs, suggesting more uniform subglacial conditions. Elongation values vary throughout GTZ1, GTZ2, GTZ3a, GTZ3b, and GTZ4, with the scaled elongation values indicating that GTZ2 has lowest elongation values ( Figure 13A). However, overall the variation between landform elongations is low between all of the GTZs (Table 1, Figure 10). This could indicate that similar processes produced most of the landforms within the study area (e.g. Ely et al., 2016). Although elongation ratios in other ice-stream corridors are highest in the centre (e.g. Stokes and Clark, 2002;Briner, 2007), elongation ratios across our GTZs show no clear organization. This is likely due to slight variations in subglacial conditions throughout each ice-flow phase that resulted from ice flowing at slower speeds than those required for such streamlined landform organization (i.e. warm-based, but no higher-velocity ice streaming), or a limitation of our data resolution. More likely, the study area does not capture the full footprint of the ice streams and therefore does not have the same pattern (i.e. our study area only captures the catchment zone of the ice stream where typical ice-stream organization has not occurred).
Spatial density of subglacially scoured lakes provides interesting insights. The highest lake counts are found at the most frequent elevations, which are mostly within GTZ1 ( Figures 7A and B). However, there is a higher density of lakes at less frequent elevations as well, notably at lower elevations and also at higher elevations (atop the De Pas batholith; Figures 7A and B) compared to the most frequent elevations where there is a higher proportion of cells with lower densities ( Figure 7B). Higher rates of glacial erosion at lower elevations in glaciated terrains have been documented in other studies (e.g. Ebert et al., 2015), but the abundance of lakes at the highest elevations within the study area suggests that areal scouring also occurred across parts of the De Pas batholith. An example of enhanced small-lake density on elevated surfaces of the batholith is shown in Figure 5A. In contrast, lake surface area appears to follow elevation frequency distribution more closely ( Figure 8B); the larger lakes are thus found in almost every GTZ (Figure 3), except GTZ4 where the highest elevations occur. The larger lakes are particularly concentrated in areas related to paleo-ice stream catchments (GTZ2 and GTZ3b). The larger lakes within GTZ2 may also indicate greater bedrock control of the Rachel Laporte and Lac Tudor shear zone on overall lake size (compare Figure 8A to Figure 2B). Taken together, these results suggest that at higher elevations, moderate/high lake densities may indicate areal scouring during more than one ice-flow phase, as suggested by the multiple sites with crosscutting striations from multiple flow directions at higher elevations ( Figure 11). It is important to note that bedrock hardness over the upland (De Pas batholith) is highest in the study area. Bedrock erosion was therefore limited, hence the abundant striae crosscutting was preserved, but enough subglacial erosion occurred to produce moderate/high lake densities in several places. This result lends additional support that cold-based areas on elevated surfaces were either restricted in time or spatially limited during the last glaciation . Sites where multiple crosscutting outcrop-scale ice-flow indicators have been documented are concentrated within the central portion of the study area. The majority (8 of 15) are within GTZ4, with six sets near the edges of GTZ2, and a single set in GTZ5 (Figure 11). Although striation measurements were largely restricted to regions of till veneer or bedrock, the partial overprinting of older ice-flow indicators by younger indicators at these locations suggests that these areas experienced enough subglacial erosion to sculpt, striate, and/or polish the bedrock surface during multiple ice-flow phases, but the younger phases did not completely remove evidence of older ice flows. The abundance of these preserved striations suggests limited overall subglacial erosion, possibly due to low basal ice velocity and high bedrock hardness. Because the De Pas batholith is a local topographic high, this region is more likely to have acted as an area of basal drag and this would have slowed down the basal ice of that area. The other sites of crosscutting ice-flow indicators are located around the edge of GTZ2 and could therefore be near the onset zone of the UBIS catchment where erosion was more limited than down-ice, in the zone of basal ice-flow acceleration (towards the northwest of the study area within the centre of GTZ2).

Relationship to ice-flow history and glacial terrain zones
We can thus now analyse the results within the context of the ice-flow history and GTZ association. The subglacial terrain zone that contains abundant indicators associated with the oldest northeast ice-flow phase  is GTZ1, in the eastern part of the study area ( Figure 11). This inferred 'relict' glacial terrain is characterized by abundant proxies of warm-based conditions that are well distributed across the zone (Figure 11). Within GTZ1, cells classified as warm-based and moderately warm-based represent 66% of the area, with only 10% of cells classified as sluggish ice. This signature appears to be controlled in part by lake density. It is useful to note here that GTZ1 has very limited recognized overprinting by younger ice-flow phases. This imprint is thus considered to largely reflect the subglacial conditions related to Flow 1. The interpretation is that Flow 1 was laterally extensive and created a broad landscape of areal scouring that created numerous small bedrock lakes (high density of lakes at low elevation in Figure 7B) and few large glacial landforms (mostly elongated glacially sculpted bedrock ridges; Figures 13B and C).
The subglacial imprint from Flow 2 is restricted to GTZ2. Flow 2 indicators only crosscut Flow 1 features locally, as it seems to have completely overprinted most Flow 1 evidence elsewhere in GTZ2 (Figure 11). The GIS results within GTZ2 could thus reflect the cumulative effect of two ice-flow phases, although most elongated glacial landforms within GTZ2 are associated with Flow 2. These results suggest that this ice-flow phase was characterized by warm-based and moderately warm-based conditions, but only in a relatively constrained area within GTZ2. The underlying bedrock structure within GTZ2, in the northern portion of the Rachel-Laporte domain (Figure 2), has a similar orientation to Flow 2's ice-flow direction (northwest), which suggests that the bedrock structure may have played a role in the subglacial erosion and geometry of the upstream portion of the UBIS (e.g. Phillips et al., 2010;Krabbendam and Glasser, 2011). Elsewhere in GTZ2, large areas appear to contain featureless or undulating till blanket. These areas were previously included within GTZ4 (Figure 3; Rice et al., 2019), but are now in GTZ2 to include all the northwest-trending ice-flow indicators on bedrock outcrops. The GIS analysis, however, clearly shows with the abundant warm-based undefined cells (green)that most features associated with these bedrock ice-flow indicators were buried under till that was deposited perhaps under different glacial dynamics conditions (i.e. meltout deposits). This is consistent with the interpretation of a change in dynamics of the UBIS catchment, which led to the stoppage of northwest-trending fast ice flow in the northwest portion of the study area, and a western shift of active paleo-ice stream catchment Margold et al., 2018). The stoppage of fast flow over GTZ2 and shift to reduced warm-based conditions is required to preserve the Flow 2 landscape, which is only weakly overprinted by Flow 3 features along the eastern edge of GTZ2. This reconstruction is further supported by the occurrence of long strips of ribbed moraines that have developed onto and deformed northwest-trending drumlins northwest of our study area (Jansson, 2005). The westward shift of the UBIS catchment is thought to have also caused the ice divide (over GTZ1 and GTZ4) to migrate westward . The ice divide migrated to, or just beyond, the western edge of the study area, resulting in ice flowing to the east (Flow 3) over part of the study area. The westward migration of the divide likely caused a transient shift from more warm-based conditions to more sluggish ice across the western portion of the study area. Till blanket could have been deposited without the development of streamlined landforms (hence the warm-based undefined cells) in some parts of GTZ2 as a result of a change from relatively fast ice flow, when the area was part of the UBIS catchment, to sluggish ice when the area was under the ice divide.
After the ice divide migrated to the west, Flow 3 overprinted portions of GTZ1, creating the GTZ3a and GTZ3b landscapes. In GTZ3a, a strong overprinting of Flow 1 is recognized, whereby large subglacial landforms with a more northeastern orientation are surrounded by smaller, more easterly landforms we associate with Flow 3. GTZ3a is characterized by several outcrop-scale ice-flow indicators associated with the northeast (Flow 1) and by numerous eastward-trending landforms associated with Flow 3. This is thus a zone of mixed landscape inheritance (Flow 1) and overprinting (Flow 3). GTZ3a does contain a small band of warm-based cells (red) that have abundant small, but elongate landforms and, interestingly, no striation evidence for Flow 3 was observed at the accessible outcrops in GTZ3a (n =12). Striation evidence of Flow 3 is limited to the central highlands ( Figure 11). This could be due to a site access bias, as only the more elevated outcrops were accessible in this area. The east-trending landforms at lower elevation could have formed through deformation of the till (e.g. Menzies et al., 2016), with the highlands acting as 'sticky spots' preserving Flow 1 striations (Stokes et al., 2007). Down-ice (east) from the cluster of warm-based cells in GTZ3a (Figure 14), a lack of landforms or striations associated with Flow 3 (marking the eastern boundary of GTZ3a) suggests an abrupt change in subglacial conditions.
The Cabot Lake ice stream, just south of the study area, also trends in an eastward direction, which suggests it likely formed contemporaneously to our Flow 3 Rice et al., 2020). Similarly, GTZ3b is thought to have occurred during Flow 3, as it contains landforms of a similar size and orientation and represents the catchment region of the Kogaluk ice stream to the northeast (Margold et al., 2015;Paulen et al., 2017a). Cell classification within GTZ3b indicates warm-based conditions transitioning to more sluggish conditions along its southern limit, indicative of paleo-ice stream catchment regions and similar to GTZ2. Flow 3 is also associated with deglaciation, as evident from the near parallel eskers to the east of Flow 3 landforms. 3171 GIS-BASED MULTI-PROXY ANALYSIS OF THE EVOLUTION OF QUEBEC-LABRADOR-LIS GTZ4 is largely limited to the central upland area characterized by numerous bedrock outcrops with crosscutting ice-flow indicators ( Figure 11). This zone thus has a mixed record of ice-flow indicators both inherited from Flow 1, as well as the younger phases Flow 3 and Flow 4. The GIS analysis of this zone suggests that the cumulative effect of these ice-flow phases led to a sluggish imprint ( Figure 13D) over most of the GTZ (53%). This suggests sluggish subglacial conditions, although it is important to note that bedrock hardness over the upland (De Pas batholith) is also considered to be the highest of the study area. In addition, there are some places within that zone where bedrock lake density is relatively high, suggesting more areal scouring (indicative of the increase in lake density at the higher elevations in Figure 7B). Nonetheless, it is quite clear that the upstream portion of the UBIS did not extend into GTZ4. It appears that the western edge of the central upland area formed a limit to the extension of the ice-stream catchment zone into the study area. Whether this was a topographic barrier, or a thermal boundary (e.g. Dyke et al., 1992;Kleman, 2005, 2008) is somewhat uncertain, although we favour a topographic control in this case as no cells were classified as cold-based (purple) within the upland.
Finally, GTZ5 is an area that is mostly characterized by till blanket with few other glacial features (96% warm-based undefined; green cells in Figure 13E). This is an area that is interpreted to have experienced a shift from actively flowing ice (few striated sites and subglacial traction till), that later switched to more sluggish ice conditions, perhaps between competing ice-stream catchment zones. The till at the surface could thus have been produced by actively flowing ice, but it was deposited through slow/sluggish ice by basal meltout processes.
Overall, the migration of the ice divide across our study area and influence from multiple paleo-ice stream catchmentscaused transient polythermal conditions, which ultimately created a fragmented mosaic of subglacial features. Some cold-based regions may have been overprinted after a switch to warm-based conditions, but extensive preservation of relict landforms in the east and northwest of the study area requires the opposite situation (i.e. a switch from warm-based conditions to much more sluggish subglacial conditions). This is exemplified by the landforms associated with Flow 1, which formed under warm-based conditions that were later preserved under either sluggish or cold-based conditions following the formation and migration of the ice divide ( Figure 14). This work supports the observations of Clarhäll and Jansson (2003) in the southern Lac aux Goélands area (Figure 2), who proposed a shift from warm to cold-based conditions, preserving older glacial landscapes and subglacial landforms in restricted subglacial zones during restricted time periods. These transient subglacial conditions and the resulting mosaic of subglacial thermomechanical zones created a fragmented landscape similar in many ways to the outer zone of the Keewatin dome, such as in the Canadian Shield terrain of northeastern Manitoba (Trommelen et al., 2012).

Evaluation of study area within a larger glaciological context
The pattern of subglacial conditions identified in this study does not match the symmetrical subglacial thermal regimes observed in other ice-divide regions (e.g. Kleman and Stroeven, 1997;De Angelis and Kleman, 2005;Davis et al., 2006;Staiger et al., 2006;Briner et al., 2008;Refsnider and Miller, 2010), which were affected by sustained cold-based subglacial conditions near an inferred stable ice-divide centre. In studies that examine 10 Be abundances from different subglacial regimes, results with 10 times the abundance of 10 Be (when compared to nearby samples over glacially eroded terrain) are typically interpreted to reflect regions of sustained cold-based conditions (Marquette et al., 2004;Staiger et al., 2006;Briner et al., 2008;Fu et al., 2018). Other studies utilizing the same methodology reported samples as low as 4 times the abundance of 10 Be to be indicative of sustained cold-based regimes (Fable et al., 2002;Briner et al., 2014). Intermediate or transition zones, between warm and cold-based endmember categories, yield 10 Be abundances 1.4 to 4.2 times higher than surrounding warm-based terrains (Li et al., 2005;Staiger et al., 2006;Fu et al., 2018).
The highest 10 Be abundances from our study are 2.1 to 2.6 times greater than what would be expected from postglacial exposure, an order of magnitude lower than those used to infer sustained cold-based conditions in other studies (see above). The highest 10 Be abundances within our study area are thus comparable to those interpreted to reflect intermediate basalthermal conditions in other studies. Similarly, the CIA values from our study (excluding samples from the Rachel-Laporte domain) ranged from 46.4 to 53.2, significantly lower than the 70.0-90.0 reported by Refsnider and Miller (2010) from a region with high 10 Be abundances and inferred cold-based landscape. The fresh bedrock of the Rachel-Laporte domain has a higher CIA (59.2), which leads to a higher threshold in that area for fresh vs. weathered bedrock/till. This means there is no weathering signal (i.e. evidence of highly weathered sediments) from the CIA results in the study area.
Taken together, we interpret our results to reflect transient subglacial conditions associated with the ice-divide migration across the study area with no zone indicative of sustained, stable, cold-based conditions for the full Laurentide glaciation. This is best explained by the migration of the ice-divide eroding enough bedrock material to remove inherited CIA values and remove some, but not all, of the preglacial 10 Be build-up over the central upland area where cold-based conditions were likely the most pervasive. Shattered bedrock and felsenmeer were observed by Dubé-Loubert and Roy (2017) approximately 140km to the north of our study area, but occur over small upland areas only. Those sustained cold-based conditions either did not extend into our study area or were eroded during warm-based conditions during subsequent ice-flow phases ( Figure 14). These results suggest that our study area was more sensitive or susceptible to larger ice-sheet dynamics, possibly due to its proximity to paleo-ice stream catchments.
Interestingly, GTZ1 has a high abundance of warm-based cells (47.7%) whose features were largely (if not solely) created during Flow 1. Ice-flow phases with similar orientations to Flow 1 have been identified across southern and central Quebec and in eastern Labrador (Klassen and Thompson, 1993;Veillette et al., 1999;Parent et al., 2004). The exact age of this extensive, yet discontinuous and 'old' flow set is unknown (Figure 4). However, the CIA and the more limited cosmogenic results (in GTZ1) do not show evidence that would support a pre-Late Pleistocene age for GTZ1 landscape, due to the low inheritance in both bedrock and till (6.5×10 4 and 7.1×10 4 atomsg À1 , respectively; Figure 11). Although uncertainties persist, as Flow 3although its imprint lies outside of GTZ1could have sufficiently trimmed the Flow 1 landscape surface to remove most of the cosmogenic inheritance. We did not find clear erosional evidence related to Flow 3 across GTZ1, but we cannot rule out the possibility of some Flow 3 erosion. In addition, it is unclear whether weathering during the last interglacial would have been enough to create a distinctly higher CIA across GTZ1. Nonetheless, the interpretation, as based on the currently available data, suggests the regional northeast Flow 1 formed at an early stage of the last glaciation rather than a previous glaciation (cf. Veillette et al., 1999). If correct, this would indicate that during the earliest glacial phase (Flow 1), the bed was more broadly warm-based than during subsequent ice-flow phases, with later ice-flow events being more localized and highly influenced by ice streaming events. Ground temperature reconstruction from several borehole temperature profiles collected across the Canadian Shield indicates that the coldest ground surface temperatures (near the pressure melting point for most boreholes) were not reached before the LGM, suggesting widespread warm-based conditions (Pickler et al., 2016); the closest deep borehole site to our study area is in Sept-Iles (~500km south), where minimum reconstructed ground temperatures reached À1.4°C shortly after LGM (Pickler et al., 2016). This may not be representative of the conditions within our study area, but it does provide some evidence that some core regions of the LIS on the Canadian Shield were mostly warm-based throughout the last glaciation, except perhaps for a relatively short time around LGM and shortly thereafter. Ice-sheet models also provide useful insights into this question, and low 'warm-based' probabilities for the entire last glaciation are only predicted for a few small regions (e.g. Tarasov and Peltier, 2007). Interestingly, one of these probable 'sustained cold-based' zones is located within our study area (Tarasov and Peltier, 2007;their figure 11). Progress in numerical modelling, integration of more detailed bedrock thermal properties, and the future development of geothermal energy in remote northern communities, which necessitates deep boreholes, may soon increase the resolution and accuracy of basal LIS temperature reconstructions for the last glaciation.

Conclusions
This research spatially analysed a glacial landscape within a core region of the LIS that represents the net effect of multiple ice-flow phases with complex interactions between paleo-ice stream catchment dynamics and local ice-divide migration that had previously been separated into distinct GTZs . Five different glacial landscape features considered to be proxies of subglacial dynamics were used, namely bedrock-controlled lake density and surface area, glacial landform density and elongation, and till blanket. The results were classified into different relative subglacial dynamic intensities (4-0; warm-based to cold-based). Whereby the warm-based cells represent areas indicative of actively flowing basal ice, moderate warm-based cells represent regions of less erosive warm-based conditions, warm-based undifferentiated cells represent regions of thick till, and thus warm-based conditions at some given time, but with no streamlined landforms indicative of subglacial deformation or erosion. Sluggish ice cells represent regions of few proxies and no till blankets, and cold-based cells represent regions of sustained cold-based conditions. Sustained cold-based conditions did not occur within the study area based on our analysis, but rather the study area was influenced by polythermal subglacial evolution associated with a transient ice-divide migration. Results show distinctly different characteristics between the newly refined GTZs, but the grid-cell classification does identify regions of erosive warm-based conditions in all GTZs. Interestingly, the most relict GTZ1 is associated with the earliest ice-flow phase identified within the study area and contains over 48% warm-based cells and is remarkably well preserved. GTZ2 is associated with the onset of the UBIS and contains warm-based cells over 24% of its area. GTZ3a is associated with ice-streaming events during deglaciation and it has warm-based cells over 31% of its area. The related GTZ3b, which represents a portion of the catchment zone of the Kogaluk ice stream, has warm-based cells over 59% of its area. Local eskers abruptly end at the eastern edge of GTZ3a, and the change in elevation associated with the western boundary of GTZ3a and the formation of eskers may be a partial control on the confinement of GTZ3a and onset of esker formation. GTZ4 is the largest GTZ and centres the study area with warm-based cells over only 3% of its area. GTZ4 represents an area that evolved under less erosive ice conditions that removed some preglacial weathering but did not erode deeply into the bedrock, as indicated by 10 Be abundances, abundant crosscutting striations, lake density, and CIA results. CIA and 10 Be abundance analysis also indicates that the oldest ice-flow phase (Flow 1 to the northeast) most likely occurred during the last glaciation as opposed to a relict flow from a previous glaciation preserved throughout interglaciation and subsequent reglaciation (cf. Veillette et al., 1999). GTZ5 contains warm-based undifferentiated cells over 96% of its area, indicating warm-based conditions at some point during glaciation, which would be required to deposit till blankets, but without the subglacial conditions necessary to produce any measurable proxies during subsequent conditions (e.g. streamlined landforms) or that were masked by deposited meltout till.
Overall, this work has shown that the subglacial dynamics of northeastern Quebec and adjacent Newfoundland and Labrador, a core region of the northeastern sector of the LIS, varied throughout the last glaciation and was subject to time-transgressive changes associated with ice-divide migration. How rapidly these changes occurred remains uncertain, but they seem to have all taken place during the last glaciation. Although more information may be required, specifically concerning the physical properties of the bedrock, such as hardness and thermal conductivities, the findings of this research provide useful constraints or a testable reconstruction for future numerical ice-sheet modelling efforts. Understanding how the reconstructed subglacial conditions relate to sediment erosion and transport would further improve our understanding of how these complex regions of the LIS evolved through the last glacial cycle.