Regional subglacial quarrying and abrasion below hard‐bedded palaeo‐ice streams crossing the Shield–Palaeozoic boundary of central Canada: the importance of substrate control

Three‐dimensional surface visualization models derived from high‐resolution LiDAR data provide new information about the type and scale of erosional processes below Late Wisconsin palaeo‐ice streams traversing the boundary between Canadian Shield crystalline rocks with offlapping Palaeozoic limestones in central Ontario. The hard bed is directly analogous to that found below ice streams in East Antarctica and East Greenland and provides insight into the effects of abrupt changes in substrate type on subglacial processes. Erosion of hard crystalline Canadian Shield rock was largely ineffectual consisting of areal abrasion of rounded whalebacks and local lee side plucking. In contrast, fast flow over the strike of gently dipping well‐bedded and jointed Palaeozoic limestones cut large flow‐parallel grooves and ridges akin to mega‐scale glacial lineations reflecting intense abrasion below narrow streams of subglacial debris dominated by hard crystalline Shield clasts (erodents). Regionally extensive plucking of structurally weak, well‐jointed and bedded limestone produced large volumes of rubbly carbonate debris leaving a 25‐km‐wide belt of uncontrolled hummocky rubble terrain (long known as the Dummer Moraine in Southern Ontario) some 350 km long and locally as much as 10 m thick. Subglacial plucking and abrasion under fast flowing ice were highly effective in stripping limestone cover rocks from Precambrian basement, and over many glacial cycles, may have played a role in the location and excavation of numerous large and deep lake basins around the Shield–Palaeozoic boundary zone in North America.

The geomorphological imprint of fast flowing ice streams on soft (sediment/soft rock) beds is relatively well studied (e.g. Stokes & Clark 1999Evans et al. 2008Evans et al. , 2014Stokes 2011;Spagnolo et al. 2014;Sookhan et al. 2018a, b). This contrastswith the situation where ice flowed over hard (rock) beds and 'mixed' beds in areas of non-existent or thin sediment cover (Eyles 2012;Eyles & Putkinen 2014;Evans et al. 2015;Krabbendam et al. 2016;Veillette et al. 2017) and where the imprint of fast flow is more difficult to decipher (e.g. Margold et al. 2015Margold et al. , 2018. In their review of modern and ancient hard bedded Antarctic ice streams Livingstone et al. (2012: p. 98) stated that it is 'surprising' that so few studies have taken advantage of exposed palaeo-ice stream beds for comparison with modern ice streams whose beds cannot be directly observed (e.g. Schroeder et al. 2014;Brisbourne et al. 2017).
This paper addresses this gap in knowledge of subglacial processes by mapping the well-exposed hard bed of Late Wisconsin palaeo-ice streams in Southern Ontario in central Canada composed of Precambrian crystalline and Palaeozoic sedimentary strata, using recently acquired high-resolution (0.5 m) LiDAR data (Ontario Ministry of Natural Resources and Forestry Ontario, 2018;Fig. 1A). The geology of the study area is directly comparable to the beds of modern ice streams in Antarctica and Greenland where up-ice areas are commonly underlain by resistant crystalline Precambrian Shield rocks with abrupt down-ice boundaries with sedimentary strata and/or soft sediments (e.g. O Cofaigh et al. 2002;Livingstone et al. 2012).
The study area (4000 km 2 ) was covered by the Laurentide Ice Sheet (LIS) during the Late Wisconsin glacial maximum and extends some 350 km along the southern periphery of the Canadian Shield and the northern limit of overlying Palaeozoic platformal sedimentary strata. This sharp boundary extends from southern Georgian Bay in the west to upper New York State east of Lake Ontario and the St. Lawrence Valley (Fig. 1). This is representative of a large part of the circumference of the Canadian Shield of North America (e.g. Stott & Aitken 1993) and is also seen in Scandinavia around the margins of the Fennoscandian Shield (e.g. zoic Boundary zone (S-PBz) to emphasize its unique geology. During regional deglaciation in Ontario after c. 14.5 ka, a large trunk ice stream (Ontario Ice Stream, OIS; Fig. 1B) was steered into the Ontario Basin from the Quebec-Labrador centre of the LIS, and was confluent with narrower, as yet unnamed ice streams flowing southward off the higher parts of the Shield (Fig. 1B). The primary objective of the present study is to employ high-resolution LiDAR mapping, combined with field investigations, to determine the role of changing substrate on subglacial processes and landforms under these fast flowing ice masses. In turn, data from Ontario throw light on broader continent-scale erosional processes operating below ice streams within the LIS that flowed from the Canadian Shield onto surrounding sedimentary rocks. This zone is characterized by numerous large lake basins whose origins are still not well understood.
Physical setting and glacial history of the study area Southern Ontario is essentially a large inter-lake peninsula (140 000 km 2 ) underlain by platformal Palaeozoic strata extending southwestward from the edge of the Canadian Shield, and bounded by the bedrock basins of lakes Erie, Huron and Ontario. The region was completely overrun by the LIS just after 23.0 ka Mulligan & Bajc 2018), when regional flow was predominantly to the southwest. Deglaciation about 14.5 ka was characterized by strong topographic control on the thinning and melting ice sheet resulting in fingerlike ice lobes flowing along the axis of each lake basin. Their existence was first recognized by Taylor (1913) using early aerial photographs to map their bounding moraine ridges, and he referred to them as 'ice streams'. Recent work identifies their palaeoglaciological role as fast flowing ice streams based on LiDAR mapping of large swaths of mega-scale glacial lineations (MSGLs) on their soft beds (Ross et al. 2006;Margold et al. 2015Margold et al. , 2018Eyles & Doughty 2016;Sookhan et al. 2018a). Sookhan et al. (2018b) extended LiDAR mapping from Southern Ontario into New York State and presented geomorphological evidence of late-stage fast flow by ice from the Ontario Basin southward into the Finger Lake basins as the palaeo-Seneca-Cayuga Ice Stream. This deposited the Valley Heads Moraine at the southern end of the Finger Lakes basin at about 14.5 ka. This event broadly coincides with the Bølling-Allerød warm phase when the eastern sector of the LIS experienced largescale re-organization and ice streams were switched on possibly because of widespread thinning and melting (e.g. Margold et al. 2018;Sookhan et al. 2018b). Another factor was the presence of large deep water bodies along the front of the ice sheet and resulting destabilization. OIS for example, terminated in an extensive ice-contact lake (glacial Lake Iroquois; Fig. 2) that stood as much as 45 m above the modern level of Lake Ontario (Pair & Rodrigues 1993) resulting in water depths that approached 300 m in the glacially overdeepened Rochester Basin. This large ice-frontal lake was short lived (>500 years; Donnelly et al. 2004;Lewis & Anderson 2020) and left a well-defined bluff with boulder beaches and thin (<3 m) offshore clays that drape much of the southern portion of the hard bed described herein. Lake Iroquois drained southwards through the newly ice-free Hudson River Valley at about~13.0 ka (Rayburn et al. 2011;Bird & Kozlowski 2016). The abrupt transition from an ice-contact water body standing at least 50 m above the level of modern Lake Ontario to a much smaller early postglacial water body is recorded in highresolution seismic reflection records of bottom sediments of the Finger Lakes and Oneida Lake in upper New York State (e.g. Mullins & Eyles 1996;Zaremba & Scholz 2019). Dyke (2004) reports a date of~12.9 ka for early Lake Ontario immediately following glacial lake drainage. Lewis & Anderson (2020) report that following glacial lake drainage early Lake Ontario was confluent with the Champlain Sea by about 12.8 ka indicating that ice had fully retreated from the hard bed and the upper St. Lawrence Valley.
Ongoing LiDAR-based work in Southern Ontario has the long-term objective of identifying the number of palaeo-ice streams and their interactions, and their role in building large inter-ice stream morainal systems, which are major sources of groundwater in a rapidly urbanizing area (e.g. Mulligan 2017;Mulligan et al. 2018;Bajc et al. 2019). In this regard, just after 13.5 ka, the northern flankof the OIS underwent a shortlived phase of flow switching and flowed northwest as the Halton Ice Stream to impinge against the Simcoe Ice Stream to build the Oak Ridges Moraine by rapid dumping of subaqueous sediment in an interlobate water body akin to the Valley Heads Moraine (Sookhan et al. 2018b).
The present paper focusses on those parts of the beds of the OIS and confluent ice streams where they flowed across Precambrian crystalline rocks of the Canadian Shield and onto Palaeozoic sedimentary rocks (Figs 1-3). The sharp boundary between these two contrasting geological terrains is marked by an abrupt change in the Fig. 1. A. Bedrock geology and topography of study area in central Canada. The hard bed of the Ontario Ice Stream (B) is composed of Precambrian crystalline strata of the Canadian Shield and offlapping Palaeozoic strata outcropping as extensive limestone plains. Outline of rubble moraine (Dummer Moraine) is also shown. White arrows show generalized ice-flow directions of fast flowing ice draining off the Shield. The downice soft bed is dominated by thick (<200 m) glacial sediment cover and extensive streamlined till plains (e.g. Peterborough Drumlin Field). B. Generalized schematic view of glacial conditions in the Ontario Basin following deposition of the interlobate Oak Ridges Moraine sometime after 13.3 cal. ka BP and retreat of the Ontario Ice Stream by calving in early glacial Lake Iroquois. Palaeo-water depths along the axis of the basin exceeded 300 m in the Rochester Basin; iceberg plough marks are depicted in Fig. 13B. [Colour figure can be viewed at www.boreas.dk] thickness and extent of glacial sediments, from sporadic thin covers on the Shield to Palaeozoic carbonates heavily obscured by coarse-grained carbonate rubble that is expressed geomorphologically as a distinctive belt of hummocky moraine traditionally named the Dummer Moraine (Chapman & Putnam 1966;Fig. 1). The distinctive geology and topography of the Dummer Moraine has been attributed to plucking of well-jointed limestones (see Leyland & Russell 1984;Mihychuk 1984;Shulmeister 1989;Marich 2016a, b, c) but the source(s) of carbonate debris and its precise origin(s) remain unknown because of a heavy cover of vegetation. The  Fig. 1A for location) showing lakes in structurally controlled valleys (e.g. Pigeon Lake, Clear Lake etc.). Large Palaeozoic limestone fault blocks (see B) whose planforms often resemble cast-iron flat irons ('drumlinoids' in text) with more rounded up-ice margins and abrupt fault-controlled leesides, have been extensively plucked to create carbonate rubble terrain of the Dummer Moraine (Figs 1, 3-9). Limestone blocks north of Buckhorn Lake and Stony Lake are outliers of limestone resting north of the Shield-Palaeozoic boundary; those to south form the serrated northern edge of glacially scoured limestone plains that pass under sediment comprising the drumlinized soft bed (Fig. 1). x-y = line of structural section shown in (B). Prominent eskers are shown by red lines. B. West-east oriented structural cross-section through drumlinoid limestone blocks and intervening structurally controlled valleys resulting from reactivation of underlying lineaments in Precambrian basement rocks. Modified from Sanford (1993). Fault blocks have commonly been glacially streamlined into large north-facing 'drumlinoids'. [Colour figure can be viewed at www.boreas.dk] major contrast in sediment cover from Shield to Palaeozoics prompted the present study to determine how subglacial processes changed across this lithological boundary.

Canadian Shield
Resistant Precambrian crystalline rocks of the Canadian Shield form a high-standing regional dome (Algonquin Highlands; Fig. 1) reflecting the presence of a northeast trending arch on the surface of the Shield having a maximum height of~500 m above sea level (m a.s.l.). The Highlands are dominated by broad belts of high-grade metamorphic rocks such as gneisses, and mylonites and large intrusive granite plutons, which are all resistant to glacial abrasion. Narrow linear zones of more easily eroded metasedimentary rocks such as marbles and metaturbidites are more easily eroded and create a corrugated surface to the Shield (Fig. 3). All these lithologies collectively belong to the Grenville Province, which records crustal accretion to ancestral North America during the assembly of Rodinia between c. 1.5 and 1 billion years ago (Rivers 1997). Northwestdirected accretion and compression is reflected in the regional structure of the Grenville Province dominated by closely spaced southwest-northeast trending shear zones and terrane boundaries (Figs 2, 3) that include faults and fractures that lie oblique to regional strike. A mountainous relief developed during the Grenville Orogeny was levelled to an undulating peneplain by the Fig. 3. Typical corrugated topography of Canadian Shield (see Fig. 2A for location) reflecting selective areal glacial scour of lineated highly strained crystalline and metamorphic rocks that form narrow north-south ridges, and softer recessive metasedimentary rocks such as metaturbidites and marbles. Note numerous oblique-to-strike structural lineaments that underwent reactivation during the Late Palaeozoic resulting in block-faulting of overlying Palaeozoic cover rocks (Fig. 2B). White arrow shows ice-flowdirection, and on later figures. [Colour figure can beviewed at www.boreas.dk] beginning of the Cambrian and prior to deposition of Palaeozoic limestone. The present Shield surface is the result of a combination of deep weathering during warm climates of the Mesozoic and Tertiary (e.g. Bouchard & Jolicoeur 2000) followed by glacial stripping of clayey regolith after 2.5 Ma (Clark & Pollard 1998) creating what has been described as a 'glacially scoured etch plain' (Twidale & Romani 2005). Areal scouring of this preglacial surface has left large areas relatively unscathed (e.g. Lindstr€ om 1988; Lidmar-Bergstr€ om, 1997; Krabbendam & Bradwell 2014; see Discussion below) and has emphasized the structure of the Shield by selective glacial erosion along lineaments and fractures, and softer highly strained lithologies such as metaturbidites and marbles ( Fig. 4; Krabbendam et al. 2016).
There is widespread agreement that Precambrian Shields in North America and Scandinavia were not major sources of sediment during times of Pleistocene ice cover (e.g. Sugden 1976Sugden , 1978Kaszycki & Shilts 1979;Hay et al. 1989), other than the removal of clay-rich preglacial regolith. The local-scale topography of the Shield surface in Ontario is little different from that preserved below Palaeozoic strata to the south (e.g. Ambrose 1964) and is dominated by rounded monadnocks and roche moutonn ees (whalebacks) of gneiss and granite with plucked lee side margins. The Shield is largely free of any significant sediment cover apart from narrow ribbon-like bodies of glacifluvial outwash, including eskers deposited along structurally controlled valleys that guided subglacial meltwaters. Coarsegrained often boulder-rich sediment is locally preserved in the lee of bedrock highs forming large crag-and-tails or bedrock-drift complexes. Any direct geomorphic imprint of fast ice flow is missing from the Shield, most likely as a result of its geological complexity and its resistance to abrasion. However, the imprint of fast flow is clearly expressed on Palaeozoic limestone surfaces lying immediately south of the Shield indicating fast ice flow over crystalline rocks of the Algonquin Highlands.

Palaeozoic carbonates
To the south of the S-PBz, the surface of the Canadian Shield is overlain by gently dipping (3-5°; 3-5 m km À1 ) Palaeozoic platformal carbonates (Figs 1, 2) of the Upper-Ordovician Simcoe Group (c. 540-450 Ma) of the Sauk megasequence. These were deposited in shallow inland seas that once covered most of North America (Laurentia) and were subsequently stripped back exposing basement crystalline rocks as the Canadian Shield. In the study area Palaeozoic carbonates form a west-east trending belt of plains, underlain by limestones and locally more dolomitic facies, whose surface is largely free of any sediment cover. Plains extend some 350 km from Georgian Bay and Lake Huron in the west to Kingston at the eastern end of Lake Ontario (e.g. Carson 1981a, b;Sanford 1993a;Armstrong & Carter 2010;Figs 1, 2). Gently dipping platformal carbonates were block faulted during the Devonian when underlying Precambrian structures were re-activated during successive Palaeozoic Appalachian orogenies during the formation of Pangea Sanford 1993a, b; Fig. 2A, B). Long term selective erosion along faults results in an irregular and strongly serrated Palaeozoic boundary with the Shield characterized by large commonly tilted blocks of carbonate typically as much as 2 km in width transverse to ice flow, and 4 km long separated by 're-entrant' basement fault-controlled valleys commonly occupied by lakes (e.g. Kawartha Lakes; Fig. 4) or eskers ( Fig. 2A). In many instances, the overall planform of these blocks resembles vintage 'flat irons' with prows rounded off by glacial scour with straight lateral margins and wide, fault controlled, commonly scalloped, lee-side margins (Fig. 2B). Livingstone et al. (2012) employed the term 'drumlinoid' for similar bedrock features on the floor of Antarctic palaeoice streams and this is employed herein. Isolated carbonate blocks and streamlined drumlinoids also occur as outliers on the Shield north of the Shield-Palaeozoic boundary varying in size from a few hundred square metres to almost 15 km 2 (Carson 1981a, b; Fig. 2A). The southern margin of carbonate plains passes under a cover of drumlinized and megalineated Late Wisconsin till resting unconformably on a thick succession of penultimate glacial (Illinoian) and younger deposits . The southern margin of the hard bed is very irregular with isolated bedrock drumlinoids and streamlined bedrock escarpments poking through the cover of glacial sediment as noted by Gravenor (1957) creating a 'mixed bed' along the southern fringes of the hard bed. Detailed geomorphic mapping of the hard bed is now possible using LiDAR data and the following section describes the analytical methods used to process raw data and generate high-resolution topographic images. Eyles & Doughty (2016) described the general topographic characteristics of the hard bed but their study and its findings were constrained by using low-resolution 5-and 10-m digital topographic data compared to the 0.5-m resolution data now available.

Methodology for geomorphological mapping
Georeferenced LiDAR data were obtained from the Ontario Ministry of Natural Resources and Forestry across an area of 4245 km 2 and consist of a series of 1km 2 non-overlapping tiles covering the northern shore of Lake Ontario with avertical accuracy of 15 cm. The tiles are composed of disk image files that are geographically referenced using the equal arc-second raster map system. Lidar Eastern Acquisition Project (LEAP; Ontario Ministry of Natural Resources and Forestry Ontario, 2018) and Ontario Ministry of Agriculture, Food and Rural Affairs (OMAFRA; Ontario Ministry of Natural Resources and Forestry Ontario, 2018) data were projected in the Universal Transverse Mercator zone 18 horizontal coordinate system with the Canadian Geodetic Vertical Datum of 1928 used for the vertical datum. The horizontal and vertical unit of measure (coordinate system axis units) for all raster grid cells is 1 metre (m). The same methodology was employed to generate all LiDAR figures used in this study. Non-overlapping 1km 2 data tiles were assembled into a seamless raster model. Each package (obtained in ZIP format) was decompressed, and its contents separated into groups of 100 files. These were imported into QGIS Desktop 3.2.2 and the mosaic raster tool was used to create a series of 100-km 2 seamless GeoTIFFs. A final LiDAR mosaic representing the entire area was rendered by repeating this mosaic process.
Topographic data were visualized using a traditional hillshade function that employs a single light source to create a shaded relief map or a multidirectional hillshade function combining light from multiple sources, more applicable to areas with marked relief where a single light source creates large shadows. These two methods were used throughout the study where best suited. The hillshade functions do not provide accurate elevation data but are instrumental in qualitatively assessing topography and were completed using the image analysis tool in ArcMap 10.7.1. The 'z factor' is also used in ArcMap as a conversion factor that accounts for the variation between units of measure in elevation vs. units of measure on a horizontal plane. It is also employed as a vertical exaggeration factor for visualization purposes but requires frequent adjustment to avoid flattening or exaggerating topography. Accordingly, to create a final image that accurately represents topography (e.g. Fig. 3) it is necessary to overlay the shaded relief map on the original raster. This was accomplished by setting the transparency of the original raster to 55% allowing the concurrent visualization of significant topography and accurate elevation on a two-dimensional plane. The original raster was then colourized to demonstrate relative elevation changes in the area using the commonly used 'stretched' method where raster cells are displayed as a continuous ramp of colours. Supporting data such as geology and surface hydrology were obtained from online databases. Once a usable image was created, the entire data set was inspected at a fixed scale of 1:25 000, in order to identify areas of interest, with a new shapefile created to serve as a boundary for each area.
As a final step, rasters were converted into highly detailed 3D models using ArcScene 10.7.1 to identify smaller, more subtle variations in topography. Clipped data were imported into ArcScene and the vertical exaggeration was set to 5 to accentuate slight changes in topography. A detailed image was created using rendering properties to shade the raster relative to the light position. The colour ramp was adjusted to reflect relative changes in elevation and the clipped raster was used as the custom elevation surface to provide the base heights and render the 3D features. This process uses the elevation data within the raster to create a 3D image from the initial 2D plane. The X and Y cell sizes are altered to match the cell size of the original raster. Using the 'illumination' settings, the azimuth and altitude were adjustedwhere the light source is projected to traverse the long axis of the landforms. The light source was adjusted to the visually optimal setting for identifying small-scale geomorphic features and is typically within the range of 25-35 degrees.

Geomorphology of the hard bed
As already related, the up-ice part of the hard bed comprising crystalline rocks of the Canadian Shield was subject to regional areal abrasion. Its geomorphology is distinctly corrugated with rounded monadnocks and elongate whalebacks that are the expression of a structurally controlled surface and selective glacial erosion of softer metasedimentary strata within resistant gneisses and granites (Figs 2A, 3). The geomorphology of the hard bed changes abruptly down-ice with the change to Palaeozoic sedimentary strata.

Rubble terrain
The dominant geomorphic feature of the surface of Palaeozoic carbonates down-ice of the Shield is a broad belt of coarse carbonate rubble (Fig. 4A-F) mapped and named previously as the Dummer Moraine (Gravenor 1957;Barnett et al. 1991;Chapman & Putnam 2007) given its resemblance in planform to a large end moraine (Fig. 1A). In detail however, LiDAR data show it to be a discontinuous blanket of hummocky, chaotically bedded carbonate rubble (Fig. 5) with a maximum north-south extent of 25 km (Shulmeister 1989;Barnett 1992;Marich 2016a, b, c) and an estimated volume of some 4.3 km 3 . The term 'rubble terrain' seems appropriate given prior usage in Alberta for broad expanses of highly irregular hummocky terrain underlain by glaciotectonized rafts of Mesozoic bedrock (see Fenton et al. 1993;Atkinson et al. 2018;Evans et al. 2020).
LiDAR mapping identifies large areas of chaotic rubble terrain made up of low-relief mounds of debris up to 2 m high, with steep-sided cone-like hummocks as high as 10 m with steep side slopes standing at, or near, the angle of repose (Fig. 5A, B). These surfaces are litteredwith limestone slabs and angular blocks (Fig. 4E, F). No preferred organization or distribution of hummocks can be identified and the term 'uncontrolled' is appropriate. In areas of discontinuous rubble terrain, underlying limestone surfaces have been extensively scoured and grooved (see below). Rubble terrain is underlain by clast-supported disorganized carbonate debris lacking any internal bedding or systematic sorting. Grain size varies from silt to boulder-sized jointbounded blocks up to detached bedding plane slabs as much as 10 m in length and up to 1.5 m thick (Fig. 4B,  D). Boulders of Precambrian gneiss and granite are conspicuous but overall, finer grained crystallinederived debris is rare having been diluted by abundant locally derived carbonate. Locally, areas of rubble terrain have a more organized 'controlled' topography in the form of a corrugated washboard topography of closely spaced (50-100 m) ridges of rubble lying transverse to ice flow (Fig. 6). The largest of these ridges is 1.4 km long, as much as 33 m in height and is a composite landform consisting of multiple superposed ridges. Locally, rubble debris has been streamlined into flutes and drumlins (Fig. 6).

Interpretation of rubble terrain
High-resolution LiDAR mapping shows that carbonate rubble is derived from the cliffed lee sides of large Palaeozoic bedrock drumlinoids with glacially scoured  (Figs 7-9). Detached bedding plane slabs and boulders of distinctive carbonate lithologies readily distinguished by colour, body and trace fossil assemblages and/or degree of dolomitization, can be 'walked back' up-ice and matched with in situ source beds in cliffs on the lee-side margins of carbonate blocks. Aprons of rubble extend down-ice from the lee sides of drumlinoids and in gross form at least, resemble large rockfalls with their headwall scarps and hummocky runout masses of rubble (e.g. Eisbacher 1979;Figs 7B, 9B). A postglacial landslide origin for headwall scarps and rubble aprons on the lee-side margins of drumlinoids can clearly be ruled out given the restricted height of the drumlinoids above the surrounding surface (<10 m), the low overall dip of beds and especially by the large extent of debris aprons down-ice of the S-PBz. Nonetheless, there is a strong possibility that the plucking process also involved the collapse of unstable lee side cliffs in open subglacial cavities (e.g. Boulton 1974) made larger by faster ice flow.
Extensive plucking of well-bedded carbonate along the S-PBz was facilitated by bedrock jointing, which results in strongly anisotropic rock mass properties effectively reducing the extensional bulk rock strength to zero (e.g. NWMO 2011) and rendering such strata highly susceptible to plucking. Bulk rock mass strength appears to have been further lowered by bed-parallel shear planes within soft shale interbeds; their role was likely analogous to that played by weakened shale beds in facilitating bedding plane failures and resulting slides in steeply dipping sedimentary rocks (e.g. Hart 2000). Hydraulic jacking by high-pressure subglacial meltwaters being forced along joints and bedding planes apparently lifted off entire bedding planes (Fig. 4D) while injecting massive, structureless and unconsolidated sand along bedding planes and joints (see Kyrke-Smith et al. 2014;Hall et al. 2020). Hydraulic jacking may reflect the flow of high-pressure subglacial waters from relatively impermeable Precambrian crystalline strata to the north onto well-jointed Fig. 6. LiDAR image of controlled rubble terrain with a large composite moraine ridge up to 33 m high composed of multiple closely spaced minor ridges formed by repeated stillstands and minor advances of the ice stream margin. In the process, rubble was also overridden and streamlined (at left) into larger drumlin-like forms and more elongated flutes, which in turn, are covered by closely spaced minor moraines of bulldozed rubble creating a corrugated washboard appearance. See Fig. 2A for location. [Colour figure can be viewed at www.boreas.dk] Fig. 7. A. Planform LiDAR topographic image of stepped lee-side slopes of limestone drumlinoid (see Fig. 2A for location) with hummocky mantle of carbonate rubble resulting from subglacial plucking. The uneven irregular surface of the drumlinoid with small steps reflects large-scale spalling and detachment of limestone beds as a result of shear through shale interbeds; hydraulic jacking is also indicated (Fig. 4D). A large Shieldderived crystalline erratic boulder is circled. B. Three-dimensional viewof limestone drumlinoid (see Fig. 2A  limestones with a much greater bulk hydraulic permeability. In this regard, it can be noted that freezing on of basal waters to the ice base has been suggested to occur when ice streams either decelerate, or shut down, or when basal meltwater experiences abrupt drops in hydraulic pressure (Bougamont et al. 2003;Christoffersen et al. 2006) such as in central Ontario when subglacial meltwater was forced through well-jointed platformal carbonates. Subglacially plucked limestone debris was incorporated in the ice base either by freezing on of these meltwaters, or possibly by the excavation of rubble from lee-side cavities undergoing closure in response to fluctuating ice-flow velocities. During down-ice transport, bedding plane slabs were broken into joint-bounded blocks and granular carbonate debris was generated by frictional contact and crushing of adjacent blocks. Crushing likely resulted in an ice-rich melange of angular blocks in a coarse gritty matrix of comminuted carbonate till (Fig. 4C). Basalenglacial carbonate debris was deposited on top of the hard bed of glacially scoured limestone by in-situ meltout as the ice margin retreated, generating an irregular, regional-scale belt of uncontrolled hummocky rubble terrain (Fig. 5A). Patches of better organized 'controlled' morainal topography (Fig. 6) with crudely streamlined ridges of rubble oriented parallel to ice flow, suggest episodic re-advances of the ice margin when chaotic rubble terrain was bulldozed and overrun to form a washboard pattern of moraine ridges, and drumlin-like bedforms.

Bedrock cut grooves and ridges
The upper surfaces of limestone drumlinoids and carbonate plains show the effects of intense glacial abrasion in the form of ice flow-parallel sets of grooves and remnant ridges (Figs 10-12). Grooves are as much as 4 km in length and are initiated at their up-ice ends on the shoulders of small step-like escarpments up to 10 m in height. These mark the north facing strikes of more resistant dolostone beds within thick gently southwarddipping carbonate strata. The planform of these source escarpments is distinctly 'digitate', resembling rounded fingertips (Fig. 11). The width and depth dimensions of grooves comprise a continuum from striations a few millimetres deep and centimetres wide but of many tens of metres in length, to shallow grooves less than 1 m wide and 10 cm deep (Fig. 12A) to large grooves as much as 3 m deep and 5 m wide and 4 km in length (Fig. 12B). These are widely expressed as strips of vegetation on the surface of plains reflecting the selective trapping of sediment and soil in grooves. Groove orientation is parallel to ice-flow direction and independent of regional joint sets, antecedent topography, or structural control (e.g. Lo 1978;Kingston et al. 1985;Andjelkovic et al. 1998;Lam & Usher, 2011;NWMO 2011) indicating the Variation in limestone lithology gives rise to changes in the cross-sectional morphology of bedrock cut grooves. Those on the backs of large drumlinoid blocks along the Shield-Palaeozoic boundary and on smaller limestone outliers to the north are cut into resistant bioclastic or nodular and commonly dolomitic carbonate facies (Lower Gull River Formation). In this situation, grooves have a more irregular, predominantly rectangular crosssectional form characterized by less well-defined interfluve ridges (Fig. 8) when compared to the smooth-walled 'half-pipe' grooves cut on finer-grained and homogeneous micritic 'lithographic' facies of the younger Upper Gull River Formation to the south (Figs 11, 12B). Those cut on relatively resistant dolomitic facies are also characterized by slightly divergent trends in the orientation of long axes, which is attributed to divergent flow pathways of ice moving across and around high standing drumlinoids, likely reflecting variations in ice velocity or ice thickness. Their much more irregular cross-sectional form is attributed to the greater role of lateral spalling and plucking of jointed and massively bedded sidewalls (Krabbendam & Bradwell 2011) compared to the smoothly abraded floors and sidewalls of those on more homogenous finer-grained carbonate facies. The same variation was noted by Smith (1948) in a study of subglacially cut bedrock grooves on limestones near Norman Wells in the Northwest Territories.

Interpretation of bedrock grooves
The grooved bedrock surfaces identified here from LiDAR imaging are part of a much broader regional assemblage of grooves and ridges found on Palaeozoic rocks immediately south of the Shield. Grooves are typical of the hard bed of the Huron Ice Stream where ice overrode resistant Silurian dolomitic carbonates of the Niagara Escarpment and Bruce Peninsula (Eyles 2012) and, also in upper New York State (Leverett 1902;Sookhan et al. 2018b). Bedrock grooving on this scale is now widely associated with fast flow (Krabbendam et al. 2016) and is a very common feature on the beds of Antarctic ice streams (Livingstone et al. 2012). Grooves near Kingston (Fig. 14) were previously attributed to erosion by catastrophic discharges of subglacial meltwater (Shaw 1988(Shaw , 1994Shaw & Gilbert 1990;Brennand & Shaw 1994) but their flow parallel orientation and linearity across the entire study area including those of New York State, are clearly the product of direct glacial abrasion very unlike the highly sinuous bedrock grooves produced by sediment-laden meltwaters (e.g. Shaw 1988;Snow et al. 1991;see discussion in Pair 1997;Krabbendam et al. 2016;Eyles et al. 2018;Sookhan et al. 2018, Sookhan et al. 2019. Significantly, the grooved limestone surfaces described here all share a common palaeoglaciological setting of having been cut into relatively soft Palaeozoic sedimentary strata immediately down-ice of resistant crystalline rocks of the Canadian Shield. This setting likely reflects enhanced abrasion of relatively soft sedimentary strata by coarse-grained Shield-derived subglacial debris dominated by hard lithologies. Large crystalline erratic boulders up to 4 m in diameter occur on the grooved backs of drumlinoids near or on the Shield and these likely acted as hard 'erodent' particles when swept across limestones under fast-flowing ice. Grooves were likely initiated as striations and expanded A B  Fig. 1A for location). Note prominent lateral berm of displaced clay (black arrow) and common overprinting of marks by later scours. This area would have been a shallow (<10 m) bank in glacial Lake Iroquois (Fig. 2) Smith 1948;Snow et al. 1991;Krabbendam et al. 2016). Flow-parallel grooves are essentially subglacial tool marks (in effect 'megastriations') cut into flat, carbonate bedding planes, and separated by remnant ridges. These are the hardbedded equivalents of mega-scale glacial lineations present on soft sediment or soft rock beds. Whilst subglacial quarrying along the S-PBz was highly effective at generating coarse rubble debris, the intense scouring of other parts of the hard bed also functioned as an areally extensive 'cement factory' where subglacial abrasion produced large amounts of fine-grained reactive carbonate silt and carbonate-saturated meltwaters (see Fairchild et al. 1993 for review).

Ice stream retreat
LiDAR mapping reveals much new information on the style of retreat of fast flowing ice from the hard bed especially where the margin of the OIS along the axis of Lake Ontario retreated eastward in standing glacial Lake Iroquois (Fig. 1B). Previous work has assumed on a priori theoretical grounds that the Ontario lobe (now ice stream) retreated eastward from the Ontario Basin as a calving bay margin in the deeper offshore waters of glacial Lake Iroquois (Thomas 1977;Pair & Rodrigues 1993; Fig. 2).

Crevasse squeeze ridges
As has already been related, moraine ridges composed of bulldozed rubble (Fig. 8) record episodic stillstands and small re-advances of retreating ice streams. This incremental style of retreat is also indicated by large numbers of small moraine ridges on that part of the hard bed flooded by glacial Lake Iroquois during deglaciation (Fig. 1B). Vreeken (1994) mapped 'ridged clay terrain' of small, closely spaced linear moraine ridges near Kingston (Fig. 13A). These were interpreted as De Geer ridges formed under a partially buoyant ice margin that retreated by the calving of tabular icebergs each summer, followed by pushing forward of the grounding line in winter. LiDAR mapping reveals a more complex form in agreement with earlier descriptions by Løken & Leahy (1964). Ridges are highly sinuous, up to 6 m high (with a mean of 3 m) with a mean spacing between 50 and 60 m. They are commonly flat-topped and often connected by low indistinct ridges creating a honeycomb pattern (Fig. 13A). Regionally, this distinctive landform element occurs along the floors of valleys incised into carbonate plains where water depths in glacial Lake Iroquois were as much as 100 m. Ridges are composed of deformed (glaciotectonized) glaciolacustrine rhythmites that Løken & Leahy (1964) argued resulted from subglacial pressing of soft clay into transverse basal crevasses ('subglacial crevasse-squeeze ridges'; see also Zilliacus 1989;Beaudry & Prichonnet 1995). Their overall form is, to some degree, comparable to crevasse squeeze ridges produced during late-stage stagnation and downwasting of the Maskwa Ice Stream resting on soft clayey sediment in western Canada (Evans et al. 2016). However, the regular spacing of ridges in Ontario is not consistent with a pressing model where some degree of irregularity in crevasse width and distribution could be expected. In some respects, the planform shape, spacing and dimensions of the ridges described here (Fig. 13A) also resemble the 'rung ridges' transverse to ice flow between MSGL ridges described by Dowdeswell et al. (2020) and so named because of the overall topographic resemblance to a ladder. Rung ridges were produced by extrusion of soft sediment from the grounding line of the retreating Larsen Ice Shelf in Antarctica as it rose and fell in response to tidal cycles. While a tidal mechanism clearly does not apply in Ontario, the broadly analogous size and form of the Canadian ridges supports a mechanism involving squeezing of clays at the front of an incrementally retreating ice stream terminating in deep water of glacial Lake Iroquois. This supports the model of Vreeken (1994) that these ridges are De Geer moraine ridges (e.g. Larsen et al. 1991;Kotilainen et al. 2012;Ojala et al. 2015;Høgaas & Longva 2018;Sinclair et al. 2018) related to repeated calving in deep water and suggests an annual retreat rate of up to 60 m. Confirmation of a calving ice stream margin is provided by geomorphological evidence in the form of iceberg scours.

Iceberg scours
LiDAR images reveal hundreds of linear to gently curved, predominantly west-east oriented grooves often with prominent lateral berms of deformed Iroquois clay (Fig. 13B). Using low-resolution black and white air photographs, Gilbert et al. (1992) and Gilbert (1994) mapped 164 predominantly west-east oriented shallow curvilinear scours across a~93-km 2 portion of the former floor of Lake Iroquois approximately 25 m above the current level of Lake Ontario between Trenton and Kingston. Scours are shallow with mean depths of <1 m because of the thin sediment cover on underlying limestones. The mean azimuthal orientation of all scours is 261°, the longest is 3.6 km and the maximum recorded width is 174 m. Bedrock is commonly exposed along their floors often accompanied by distinct lateral berms of displaced clayey sediment (Fig. 13B). LiDAR mapping greatly expands the number and density of scours and identifies many generations of scouring events where individual scours are crossed by younger forms. Reconstructed water depths in areas of mapped scours using palaeo-shoreline elevation data indicate water depths between 9 m on bank tops and as deep as 90 m (e.g. Pair & Rodrigues 1993) suggesting the presence of large icebergs. Gilbert et al. (1992) favoured scouring by masses of lake ice being blown westward by strong katabatic winds away from the ice front but an icebergrelated origin is consistent with calving in deep water of glacial Lake Iroquois (Thomas 1977;Pair & Rodrigues 1993; Fig. 2).

Discussion
The effect of changing substrate on subglacial processes This study provides a good illustration of the importance of substrate control on the geomorphic imprint of fast flow. Reconstructions of palaeo-ice streams flowing over hard beds under the LIS are hampered by the absence of definitive geomorphic criteria for recognizing fast flow over hard crystalline substrates (Stokes & Clark 1999;Stokes 2011Stokes , 2018Margold et al. 2015Margold et al. , 2018. The data presented here show that fast ice flow on the Shield in central Ontario occurred across a corrugated surface of crystalline rocks consisting of smooth whalebacks likely aided by high-pressure water at the ice base. A broad zone of up-ice areal abrasion and limited plucking on the Shield is fringed by a narrow linear zone transverse to ice flow of regional-scale plucking of the trailing edges of faulted Palaeozoic carbonate blocks (drumlinoids). Strongly focussed subglacial abrasion by debris streams cut flow-parallel grooves whose morphology was controlled by facies variations within underlying carbonates. The export of large crystalline 'erodent' clasts from the Shield likely enhanced subglacial erosion of much softer limestones down flow, and likely produced large volumes of fine carbonate silt. The geology of the hard bed described here from Ontario and the range of subglacial processes are fundamentally the same as those found under many ice streams flowing off the Antarctic Craton of East Antarctica and that of eastern and northern Greenland (an extension of the Canadian Shield), which all share a common geological heritage of having been conjoined within Rodinia in the Meso-Neoproterozoic and later, in Pangea (Elliot et al. 2015;Brisbourne et al. 2017). Ice streams flow over large drumlinoids of sedimentary strata along the outer boundary of the East Antarctic Craton, which are also associated with deep bedrock-cut grooves (see O Cofaigh et al. 2002;Livingstone et al. 2012: table 3 and references therein). It is possible that these forms record abrasion by hard crystalline-dominated debris being rapidly swept across softer strata. Similarly, the glacially scoured and grooved bedrock surfaces along the S-PBz in Southern Ontario are most likely long-lived surfaces that have evolved over repeated glacial cycles. In this respect, Livingstone et al. (2012) suggested that 'heavily eroded bedrock exhibiting grooving and streamlining could potentially indicate a legacy of repeated ice streaming over several glacial cycles'. Field evidence suggests this is indeed the case in Ontario because the limestone hard bed passes southward under till of the soft bed (Figs 1, 3). This indicates that parts of the carbonate hard bedwere buried andwere 'softened' by a cover of glacial sediment. In turn, as the soft bed itself undergoes erosion along its northernup-ice boundary the underlying hard bed is re-exposed to glacial erosion. The repeated alternation of softening by burial and exhumation by erosion has likely played out on numerous occasions from one glaciation to another or even within individual glaciations as ice evolved from early steady-state flows to late-stage streaming. In a similar fashion, large ice-marginal lakes and their finegrained deposits were likely instrumental in softening large portions of the Shield possibly triggering and promoting fast flow (e.g. Stokes & Clark 2003, 2004. In this fashion, what constitutes a 'hard' bed and a 'soft' bed under large ice sheets evolves through time.

Debris fluxes from hard beds and the role of subglacial plucking
In central Ontario, extensive subglacial plucking of jointed limestones occurred along some 350 km of the S-PBz leaving a broad belt of rubble terrain as much as 25 km wide (Fig. 1A). Plucking is a widely reported subglacial erosional process (e.g. Evans et al. 1998;Glasser et al. 2020) and commonly observed in subglacial cavities below valley glaciers (Boulton 1974;R€ othlisberger & Iken 1981;Anderson et al. 1982;Lliboutry 1994;Rea & Whalley 1994. The results reported herein represent a striking departure from previous considerations of the importance of this process by recognizing that plucking can occur on a much larger regional scale. Additional evidence of the importance of spatially extensive subglacial quarrying has also recently been described from the Fennoscandian Shield in eastern Sweden (Hall et al. 2020) where it was referred to as 'glacial ripping'. This process involves hydraulic jacking along subhorizontal fractures in basement rocks and produces extensive boulder spreads as relatively narrow belts parallel to ice flow. More examples of the effects of large-scale quarrying both parallel to, and transverse to ice flow, are likely to emerge as the availability of LiDAR data allows detailed mapping of former hard beds.
The Shield-Palaeozoic boundary zone in Ontario likely functioned as a long-term source of detrital carbonate. Glacial deposits in the region are dominated by carbonate clasts and fine reactive silt, which has commonly resulted in postglacial cementation of tills and glacifluvial deposits. Angular limestone rubble clasts derived from the northern exposures of Palaeozoic limestones are a conspicuous visual component of tills in Southern Ontario but there has been little consideration of where (and how) such large volumes of carbonate were eroded and the resulting landscape effects in source areas. The generation of such large volumes of carbonate debris can be considered as an initial step in the till-forming process under the LIS. In this respect, Marich (2016a) mapped the down-glacier evolution from primary carbonate rubble to immature coarse-grained clast-rich till (mapped as 'Dummer till') evolving down-glacier to matrix-supported till facies resulting from the progressive down-ice comminution of carbonate debris and its mixing with overridden glacifluvial and glaciolacustrine sediments (e.g. Northern Till; Boyce & Eyles 2000;Mulligan & Bajc 2018).

Broader relevance to modelling erosion around the margins of the Canadian Shield
The results presented here may have a wider significance for the study of glacial landscapes across North America. The large volume of rubble (<4 km 3 ) produced by plucking where ice traversed the S-PBz is exceptional and not considered by current, largely theoretical models of glacial landscape evolution (e.g. Iverson 1991Iverson , 2012Hildes et al. 2004;D€ uhnforth et al. 2010;Ugelvig et al. 2016). However, the S-PBz in Ontario is representative of some 7000 km of terrain along the circumference of the Canadian Shield (Fig. 14) where offlapping Cambrian-Devonian carbonate rocks form gently dipping homoclinal successions where more resistant facies give rise to a gentle 'cuesta and dip slope' topography with numerous up-ice facing escarpments (Stott & Aitken 1993). This zone is notably also characterized by many hundreds of bedrock-cut lake basins including the five Great Lake basins, Great Bear, Great Slave, Athabasca, Winnipeg, and Lake Huron whose age and origins are still not unequivocally established and surprisingly little studied in the modern era. White (1972White ( , 1988 proposed that the Canadian Shield experienced deep glacial erosion during the Pleistocene but this was rejected by Gravenor (1975) who argued instead, for minimal glacial erosion of the Shield in agreement with earlier views that it had been 'negligible' (Ambrose 1964: p. 851); the current view is that apart from deep fjord incisions around part of its margin, the mean depth of erosion of crystalline rock has likely been less than 80 m (Sugden 1976(Sugden , 1978Hay et al. 1989). What was lost in the discussion of White's thesis, however, was his assertion that glacial erosion was an effective stripper of Palaeozoic rocks from the surface of the Shield. He suggested that the northern limit of Palaeozoic carbonates was peeled back downdip during successive glaciations creating what he called an 'arc of exhumation' marked by large lake basins (Fig. 14). Indeed, the geographical scale and intensity of subglacial plucking along the S-PBz in Ontario now identifies a highly effective mechanism for rapid downdip stripping of jointed Palaeozoic cover strata as suggested by White (1972White ( , 1988. It can be noted for example, that the Shield-Palaeozoic boundary on the floor of Lake Winnipeg is marked by a prominent limestone escarpment (Todd et al. 1998); similarly, Hall et al. (2020) Fig. 14) and its description closely resembles the Dummer rubble terrain described herein. Lithological data in support of the thesis of enhanced stripping of Palaeozoic carbonates from the Shield were presented by Hicock (1988) who identified carbonate-rich till north of Lake Superior, which he suggested implied enhanced quarrying of carbonate debris as an ice stream moved across Palaeozoic carbonates within James Bay. Notably, Smith (1948) described extensive glacially cut grooves on limestones in the Northwest Territories (C in Fig. 14) no different from those described from the part of the S-PBz in Ontario that we interpret as rock-cut mega-scale glacial lineations. Widespread glacial erosion of carbonate strata is also suggested by dispersal fans of carbonaterich till that extend southwestward onto the Shield from areas underlain by Palaeozoic limestones south of Hudson Bay (Karrow 1992) and by other prominent rubble-rich dispersal fans of detrital carbonate in the Canadian Arctic such as on Melville Island, Prince of Wales Island, King William Island, and across southern Baffin Island (Andrews & Miller 1979;Dredge 2008;D in Fig. 14). The same Palaeozoic carbonate strata in the Hudson Bay Basin are the primary source of detrital carbonate in Heinrich Layers in Atlantic Ocean sediment cores linked to subglacial erosion below a large ice stream along Hudson Strait (e.g. Andrews & Tedesco 1992;Andrews et al. 1995;Hemming et al. 2000;Hodell & Curtis 2008;Naafs et al. 2011;Margold et al. 2018).
In summary, the literature cited above collectively support a hypothesis of enhanced glacial erosion of Palaeozoic limestones around the margins of the Shield in North America. Plucking was likely focussed in those areas of the S-PBz where underlying basement structures had been reactivated resulting in block faulting of overlying strata thereby increasing their susceptibility to plucking and the formation of lake basins. The latter, in turn, likely triggered instability of ice margins and topographically constrained fast flow off the Shield in successive glacial cycles. Additional detailed assessment of this hypothesis is now warranted using modern highresolution topographic data sets of the geomorphology of the S-PBz and lake bathymetry around the circumference of the Canadian Shield where it abuts carbonate platforms.

Conclusions
High-resolution topographic images generated by processing of newly available LiDAR data reveal new details of the hard crystalline and sedimentary rock bed of fast flowing ice within the Late Wisconsin Laurentide Ice Sheet in central Ontario. The imprint of fast flow across Shield crystalline lithologies is minimal and fast flow may have been aided by smoothed crystalline bedrock whalebacks; in this case subglacial erosion was dominated by areal scour with selective erosion along softer lithologies and fractures. Subglacial conditions changed dramatically across the Shield-Palaeozoic boundary zone where enhanced plucking and crushing of jointed and bedded limestones by fast ice flowoccurred on the lee sides of large limestone blocks (drumlinoids). Plucking was in turn, succeeded down-ice by focussed glacial scouring of limestones surfaces marked by selective linear erosion possibly by hard Shield-derived crystalline boulders (erodents) to form a continuum of tool marks from striations to large grooves ('megagrooves') identified as bedrock-cut mega-scale glacial lineations. The likely role of enhanced bulk hydraulic conductivity of Palaeozoic carbonates in enhancing hydraulic jacking and freezing on of debris to the ice base can also be identified.
The Shield-Palaeozoic boundary zone in Canada is a well-exposed analogue for the hard beds of modern ice streams in Antarctica and Greenland composed of Precambrian crystalline cratonic rocks and offlapping softer sedimentary strata, and thus provides insights into the range of associated subglacial processes occurring today. The same plucking/scouring processes were also likely widespread during successive Pleistocene glaciations in North America and may have resulted in rapid stripping of Palaeozoic cover rocks from Precambrian basement, and the excavation of large lake basins. In turn, the topography of these basins likely controlled the location of ice streams within successive ice sheets and underwent further deepening. More broadly, this study highlights the importance of geoinformatic and associated analytical techniques in advancing the understanding of ice-sheet beds and reiterates the need for further research into the application of more advanced 3D computer vision techniques to glacial geomorphological systems. manuscript by J. Campbell, S. Livingstone, M. Margold and J. A. Piotrowski.
Author contributions. -All co-authors shared equal responsibility for concept development, fieldwork, analysis, andwriting. High-resolution LiDAR data were processed by SB.