Regional subsurface mapping of glacial depositional systems preserved in buried bedrock paleovalleys, and quantitative analysis of new LiDAR imagery of surface glacial landforms using machine learning techniques, when combined, are powerful tools for assessing the dynamics of the Laurentide Ice Sheet (LIS) during the last (Wisconsinan) glaciation in southern Ontario. While age dating of deposits preserved below Last Glacial Maximum tills (LGM: marine isotope stage (MIS) 2 < c.24 000 years B.P. (ybp)) is still sparse, newly available sedimentological data derived by cored drilling, combined with legacy outcrop data, identify thick (100 m+) successions of glaciolacustrine sediments and a lack of till(s), indicating that the ice sheet margin did not extend beyond the Niagara Escarpment at the western end of Lake Ontario, during the earliest phases of the glaciation (MIS 4) or the ensuing mid-Wisconsinan (MIS 3). Ice was able to extend into New York State blocking the Rome outlet to the Hudson Valley ponding deep proglacial lakes in the glacio-isostatically depressed Huron–Ontario–Erie basins recorded by thick glaciolacustrine sediments in paleovalleys. These were cannibalized by an expanding Late Wisconsinan ice sheet after ∼24 000 ybp recorded by extensive till sheets resting on a marked erosional unconformity, with drumlinized surfaces. Analysis and visualization of LiDAR data identifies discrete statistically validated flow sets of highly elongated streamlined bedforms (mega-scale glacial lineations (MSGLs)). These provide key evidence of a major reorganization of the ice sheet margin during deglaciation into lobate paleo ice streams shortly after 17 400 ybp. MSGLs are cut across earlier LGM drumlinized tills creating widespread “palimpsest” surfaces. At least two principal phases of fast ice flow can be identified, marked by large fluxes of sediment and the rapid building of large gravel and sand-dominated moraine complexes within interlobate depocentres, the largest glacial landforms in southern Ontario. Analysis of LiDAR data further reveals the common presence of DeGeer moraines where ice margins retreated in water, and iceberg scours. Future work using LiDAR mapping has the objective of fully documenting the number, extent, and timing of ice streams to enhance glaciological modelling when the ice sheet rapidly lost mass.

Glaciological modelling of ancient ice sheets is fundamentally reliant on temporal, stratigraphic, and geomorphological data that constrain the timing, extent, volume, and changing dynamics of the ice sheet. In this regard, the growth and evolution of the last (Wisconsinan) Laurentide Ice Sheet (LIS) in central Canada, at the conclusion of the last interglacial and its initial expansion during the early Wisconsinan (marine isotope stage; MIS 4) through succeeding middle (MIS 3) and late Wisconsinan (MIS 2), is still poorly known (see Dalton et al. 2020, 2022 and references therein). This paper adds to the understanding of this interval by describing recent developments in subsurface stratigraphic analysis and surface geomorphic mapping that together provide new perspectives on the extent of the ice sheet during the last glacial cycle, during its early and middle phases, and especially, the dynamic reconfiguration of this part of the LIS during its deglaciation.

This study draws on the results of deep (+200 m depth) subsurface sediment drilling by the Ontario Geological Survey and subsurface geological and geophysical data collected from subway, road, and landfill construction projects in Canada’s most rapidly urbanizing region (e.g., Interim Waste Authority 1992; Boone 1996; Boyce and Koseoglu 1996; Bajc and Shirota 2007; Burt and Dodge 2011, 2016; Bajc et al. 2012, 2019; Burt 2013, 2017; Burt and Webb 2013; Burt and Chartrand 2014; Mulligan 2014, 2017; Campbell and Burt 2015; Crow et al. 2017a, 2017b; Sharpe et al. 2018). This analysis provides insights into the subsurface stratigraphic and paleoenvironmental record preferentially preserved in the principal buried bedrock paleovalleys of southern Ontario. This record contrasts with sparse pre-Last Glacial Maximum (LGM) tills sediments in the glacially scoured Huron, Ontario, and Erie basins, which are locally overdeepened well below modern mean sea level (Fig. 1).

We summarize recently acquired subsurface data collected from selected paleovalleys using the depositional systems approach, widely used in the petroleum industry, for synthesizing complex deposits of the last glaciation (Wisconsinan) across southern Ontario. The resulting stratigraphic framework clarifies the extent of the LIS after the last interglacial through to the LGM and deglaciation. In addition, the interpretation of the uppermost parts of stratigraphic sequences has benefitted from quantitative analysis and visualization of high-resolution (≤1 m) LiDAR imagery that reveals new geomorphic detail of glacial landforms. This permits a refined understanding of the dynamics of the LIS, especially during its final phases of deglaciation when ice was thinning, and highly dynamic. As we shall show, it provides key insights into the mode of construction of large moraine complexes, such as the Oak Ridges Moraine (ORM), deposited in inter-ice stream depocentres (Sookhan et al. 2018). Analysis of LiDAR imagery reveals that many glacial landscapes are palimpsest in nature recording partial overprinting of older LGM landforms by ice streams during deglaciation.

Southern Ontario is essentially a narrow peninsula underlain by a southwest-dipping succession of Paleozoic sedimentary strata (predominantly limestones and shales) of Ordovician–Devonian age (Fig. 2; see Armstrong and Carter 2010; Armstrong and Dodge 2010, for a full description of stratigraphic nomenclature). The succession “youngs” to the southwest with successively younger Paleozoic strata occurring downdip to the southwest (Fig. 2). The peninsula is flanked by glacially overdeepened lake basins (Huron–Georgian Bay, Erie, and Ontario; Fig. 1). Long-term Phanerozoic differential weathering of alternating shales and carbonates of the peninsula has created a subdued cuesta topography. This is now largely buried by varying thicknesses of glacial sediments, locally exceeding 200 m thickness in former pre-Pleistocene fluvial valleys and (or) at surface as large moraine complexes (Figs. 1 and 3).

The southwest–northeast alignment of the peninsula reflects draping of the Paleozoic cover over the axis of a prominent arch in the underlying Precambrian basement peneplain (Findlay–Algonquin Arch; Fig. 3). This arch separates the Michigan and Appalachian foreland sedimentary basins to the northwest and southeast respectively, which formed in the interior of eastern North America accompanying the assembly of Pangea (see Brunton 2009 and refs. therein). The Algonquin Arch gives rise to a prominent bulge in surface topography, known informally as the “Ontario Island” in early investigations in the region (Leverett and Taylor 1915; Chapman and Putnam 1966; Chapman 1985). Flanked by overdeepened Great Lake basins, the topographic high was a major influence on regional ice flow during successive phases of ice advance during the Late Wisconsinan, and during subsequent deglaciation when the ice sheet margin was much thinner and markedly lobate. At the eastern end of Lake Ontario, another basement arch (the northwest–southeast-trending Frontenac Arch; Fig. 2) exposes Precambrian crystalline strata of the Frontenac–Adirondack Belt (Carr et al. 2000) forming the higher standing bedrock sill separating Lake Ontario from the St. Lawrence River.

Topography in the northern part of the study area is dominated by the rugged high relief terrain of the Algonquin Highlands (Fig. 1) underlain by resistant Precambrian crystalline rocks of the Central Gneiss Belt (Carr et al. 2000), which reaches a maximum elevation of 586 m above mean sea level (m.a.s.l.). This high standing massif, together with adjacent uplands such as the Algoma Highlands, the Adirondack Mountains in the USA, and the Laurentian Highlands of Quebec (Fig. 1) functioned as so-called “gatekeepers” during successive Pleistocene glaciations forcing ice to flow through low-lying “saddles” such as the Nipissing–French River lowlands, the Ottawa–Bonnechere Graben, and St. Lawrence Valley (Fig. 1). At a glance, a qualitative visual assessment identifies a close geographic relationship between the location of such structurally controlled low-lying saddles and overdeepened sections of the Huron, Ontario basins and to some degree that of the Lake Erie basin. This suggests that enhanced erosion of offlapping Paleozoic strata has been the product of preferred, long-term topographic steering of ice through low-lying saddles by basement highs in the Shield. For example, the location of the overdeepened Rochester Basin in Lake Ontario (244 m below mean sea level; m b.m.s.l. Figs. 13) likely reflects ice from Quebec and Labrador flowing around the eastern flanks of the Algonquin Highlands. This overdeepened basin likely played a major role in guiding ice into the Finger Lakes of upper New York State, where the bedrock floor of Cayuga Lake has been eroded to 306 m b.m.s.l. (Mullins and Eyles 1996; Sookhan et al. 2021). In the same manner, glacial erosion of the Lake Huron and Georgian Bay basins may have been enhanced by ice forced southwestward around the Algonquin Highlands through the Ottawa–Bonnechere Graben and Nipissing–French River lowlands. In this respect, the deepest part of Lake Huron (229 m b.m.s.l.; Figs. 13) occurs downflow of the Nipissing–French River lowlands close to where the has been breached by glacial erosion (the Main Channel; Fig. 1).

The spatial distribution of glacial sediments across southern Ontario is remarkably disjunct (Fig. 4). Little or no sediment occurs in the north across the Shield other than glaciolacustrine deposits, gravelly valley fills, and eskers. The presence of rare, isolated till-cored bedforms may suggest a more extensive sediment cover in the past that has been largely stripped. In contrast, a thick glacial sediment cover begins abruptly along the Shield–Paleozoic boundary, suggesting a fundamental and immediate change in glacial erosional processes (and sediment production) once ice moved off the Shield onto softer, more easily eroded Paleozoic sedimentary substrates (Fig. 4). Bukhari et al. (2021) identified the importance of enhanced subglacial quarrying of Paleozoic sedimentary strata, especially jointed, well-bedded limestones, which produced large volumes of carbonate debris that was reworked southwards below the ice sheet. This contrasts with lower rates of erosion across crystalline strata of the Shield in response to areal scouring and localized quarrying. Regional-scale quarrying of the northern limit of Paleozoic strata is reflected in large spreads of hummocky debris formerly mapped in south-central Ontario as the Dummer Moraine (Fig. 4). Exposed carbonate plains immediately south of the Shield, extending from Manitoulin Island in Lake Huron in the west, to the eastern end of Lake Ontario on the flanks of the Frontenac Arch, are extensively grooved by subglacial abrasion. North-facing escarpments have been heavily modified by glacial erosion (Fig. 4; see Discussion below).

Counterintuitively, the oldest Pleistocene sediments in southern Ontario are selectively preserved not along the floors of the Erie, Huron, and Ontario basins but on their landward margins within paleovalleys incised into softer Paleozoic bedrock (Fig. 3). The consensus view is that these valleys are the product of differential erosion of sedimentary strata of varying lithology, which largely accomplished during the pre-Pleistocene when southern Ontario underwent regional uplift (e.g., Davis et al. 2020). These valleys formed part of a larger preglacial mid-continental drainage system that underwent modification and enhancement during subsequent glacial–interglacial cycles (Mulligan and Bajc 2018; Sharpe et al. 2018; Naylor et al. 2021; Sears and Bacanek 2022). This has given rise to a series of escarpments and recessive strike valleys underlain by weaker bedrock, now partially or fully filled by glacial sediments (Figs. 3, 5A, and 5B).

The largest bedrock paleovalley, the Laurentian paleovalley, occurs where subcropping shales of the Ordovician Blue Mountain Formation form the floor of a broad (30–50 km wide) elongate basin linking Georgian Bay with the Lake Ontario basin (Fig. 3). Similarly, a sediment-filled bedrock depression (Ipperwash paleovalley; Fig. 3) connects the Erie and Huron basins cut into soft Hamilton Group shales (Fig. 3). Easily eroded interbedded evaporites, dolostone, and shale of the Silurian Salina Group form the floor of a prominent bedrock depression west and south of the Niagara Escarpment that can be traced northwestwards into the Walkerton Trough. Its southernmost margin is defined by the Onondaga Escarpment marking the strike of resistant Devonian carbonates (Fig. 3).

The Niagara Escarpment is the most prominent bedrock topographic feature (Figs. 1 and 3). It is the surficial expression of highly resistant Silurian reef dolostones, emerging from under the thick Pleistocene sediment cover preserved within the Laurentian paleovalley (Fig. 5A). The escarpment has been severely modified by glacial erosion and shows several narrow V-shaped “re-entrant” valleys, some now partially filled by sediment such as those along the Niagara Peninsula near Hamilton (principally the Dundas Valley and the Erigan Channel near St. Catharines; Fig. 3), which again, likely have pre-glacial antecedents as part of a mid-continent drainage system (e.g., Spencer 1890). The strike of the Niagara Escarpment north of the Main Channel is noticeably offset (to the northeast) from that part south of the Grenville Front Tectonic Zone (Fig. 2) indicating the role of basement reactivation, fracturing of overlying Paleozoic strata, and locally enhanced glacial erosion of fractured rock by ice directed into the Huron Basin through the French River–Lake Nipissing lowlands, as related above.

Fracturing of Paleozoic strata above reactivated basement faults and terrane boundaries is recognized as a key control on the location of linear lake basins in south-central Ontario (e.g., namely Balsam Lake and Pigeon, Chemong, Clear lakes, and Lake Simcoe; see Sanford et al. 1985; Fig. 1 and figs. 2A and 2B in Bukhari et al. 2021). These finger-like basins show broad similarities with the much better-known Finger Lakes of New York State and lie along the trend of underlying Precambrian terrane boundaries and shear zones resulting from accretion of terranes during the ∼1.5–1.0 Ga Grenville Orogeny (Fig. 2; Eyles et al. 1993; Boyce et al. 2002). The most prominent basement structure is the Central Metasedimentary Belt Boundary Zone (CMBBZ: Fig. 2), which now forms a prominent south–southwest-trending structure that influences surface topography, forming the drainage divide between Lake Huron and Lake Ontario. Terrane boundaries are exposed to the north across the Algonquin Highlands and their buried extensions to the south under Paleozoic cover strata marked by strike–slip faults with significant offsets due to structural reactivation of subcropping Precambrian structures during the mid-Paleozoic Appalachian Orogeny (Sanford et al. 1985).

Existing frameworks for the Late Pleistocene stratigraphy of southern Ontario are complex, reflecting a long history of investigations employing different methodologies, temporal frameworks, and conceptualizations of glacial sedimentation, climate change, and ice sheet structure and behaviour (e.g., Taylor 1913; Chapman and Putnam 1966; Karrow 1987; Barnett 1992). This complexity is also due, in no small part, to a lack of subsurface data (especially cored boreholes). Many investigations have been dictated by the practical need to complete surface mapping of the region at 1:50 000 scale, accompanied by simplified summaries of stratigraphic successions, to identify and inventory resources such as aggregates, and to determine the thickness of the drift cover (Karrow 1963, 1967, 1968, 2004a, 2005). This approach has given rise to a plethora of local stratigraphic names from one area to another requiring careful navigation during paleoenvironmental, hydrogeological, and geoengineering investigations. A major complication has been the adoption of a simple “layer cake” lithostratigraphic approach based on the North American Stratigraphic Code (1983) originally adopted for simple flat-lying Paleozoic strata the in mid-continent. In this scheme, formally defined lithostratigraphic units are organized hierarchically into members, formations, and groups, related simply to synchronous fluctuations in water depths resulting from global eustatic variations ( and many references therein). However, the succession of lithostratigraphic units has been shown to reflect more complex plate tectonic processes and basin formation mechanisms (e.g., Miall 2016 and many references therein). In similar fashion, early lithostratigraphic work on Quaternary deposits in southern Ontario adopted the same approach and assumed a simple regional layer cake stratigraphy controlled by regional climates and synchronous fluctuations in ice flowing through widely separated lake basins (Dreimanis and Karrow 1972; Karrow et al. 2000). Tills were assigned to formally recognized cold climate stades, and glaciolacustrine or glaciofluvial sediments, to warm interstadials. This simple approach does not lend itself to detailed investigations of complex glacial deposits and environments characterized by markedly diachronous ice advances and retreats and marked spatial variations in depositional settings and ensuing depositional facies.

A depositional systems approach

Stratigraphic data, whether derived from the study of outcrops, drilling, and seismic reflection and geophysical downhole logging (along with the corresponding legacy of nomenclature), can be simplified and rationalized into a succession of stacked sediment “packages” (“depositional systems”; Galloway 1998; Figs. 5 and 6B). These packages are defined by Miall (2016) as ‘a three-dimensional array of sediments or lithofacies that fill a basin. Depositional systems vary according to the types of sediments available for deposition as well as the depositional processes and environments in which they are deposited. Each system is composed of a wide range of genetically related three-dimensional lithosomes composed of variable facies types, and are defined and bounded top and bottom, by major region-wide unconformities. Each system recognized in the present study essentially constitutes part of a regional “event stratigraphy” across southern Ontario that records distinct phases of sedimentation and erosion brought about either by extrabasinal changes in climate or elevation brought about glacio-isostatic movements, or more local intrabasinal controls brought about by changing positions of ice margins and the size and depth of ice frontal lakes etc. The subsurface distribution and stacking of genetically related depositional systems define the broader “depositional architecture” of any one sedimentary basin (see Galloway and Hobday 1996; Miall 2016).

The depositional systems approach lends itself to the analysis of complex glacial stratigraphic records where extrabasinal controls such as glacio-isostatic crustal movements, diachronous fluctuations in the extent of ice margins, and (especially relevant in the study area) changes in lake level and lake extent. As a consequence, individual lithostratigraphic units (such as tills) are spatially discontinuous or evolve spatially into genetically related facies, defying simple “like-for-like” correlations implicit in more traditional stratigraphic schemes. Instead, individual units are enveloped within much broader sedimentary packages (depositional systems) bounded by unconformities that record large-scale events across glaciated basins. The reader is referred to detailed discussions of this so-called “allostratigraphic” approach, which underpins the use of depositional systems in Frazier (1974), Bhattacharya and Posamentier (1994), Vanderburgh and Roberts (1996), Galloway (1998), Walker and Wiseman (1995), Miall (2016), and Eyles et al. (2022). This approach is widely used in reservoir delineation in oil and gas exploration and is readily adapted to groundwater investigations in Pleistocene glaciated terrains requiring mapping of complex subsurface aquifers and aquitards (e.g., Burt and Dodge 2011; Burt 2018).

In the first application of this approach to Late Pleistocene sediments, Eyles (1987) outlined a series of interglacial and glacial depositional systems for the Greater Toronto Area (GTA). This approach used then available subsurface data from subway construction projects and from outcrops along the Lake Ontario shoreline at Scarborough Bluffs that expose sediments of the Laurentian paleovalley (Fig. 6). This scheme can now be refined and extended regionally in light of new geophysical and drilling data generated during the search for new waste management facilities for the GTA (Boyce et al. 1995; Boyce and Eyles 2000) and by deep stratigraphic drilling data collected from the infills of the principal paleovalleys by the Ontario Geological Survey and Geological Survey of Canada (see reviews in Bajc et al. 2018; Burt 2018, 2020; Mulligan and Bajc 2018; Sharpe et al. 2018; Fig. 5C).

The Pleistocene stratigraphy of southern Ontario can be broadly subdivided into seven principal depositional systems (DS I to VII; Figs. 5 and 6), which are described in detail below. In brief, each system records a distinct phase in the glacial geological evolution of southern Ontario extending from the Illinoian glaciation (MIS 6), through the ensuing Sangamonian interglacial (MIS 5) and Wisconsinan glaciation (MIS 4, 3, 2), to final deglaciation (Figs. 5 and 6).

It is emphasized that for simplicity, the subsurface portion of this paper is focussed primarily on the Laurentian paleovalley where there is a rich subsurface database and where the depositional system approach was originally applied (Eyles 1987). Assessment of the wider Ontario Geological Survey subsurface database is ongoing, and preliminary findings indicate that it applies equally well to paleovalleys outside of the focus area. It is anticipated that the depositional system approach will eventually prove to be an effective way to consider Pleistocene sedimentary systems outside of Ontario.

The principal depositional systems are described in the following section in order of age from oldest to youngest. In brief, the oldest glacial deposits so far recovered (principally tills of Illinoian glaciation age; MIS 6: DS I; Figs. 5 and 6) rest on weathered or glacially scoured Paleozoic bedrock (Burt 2018; Mulligan and Bajc 2018). In turn, Illinoian tills are succeeded by warm climate fluvial–lacustrine Sangamon age deposits (DS II; Figs. 5 and 6). These pass upwards conformably into cold climate deeper water fine-grained glaciolacustrine and deltaic facies that record damming of the St. Lawrence Valley during early expansion of the LIS (DS III) during MIS 4 (Eyles and Williams 1992; Williams and Eyles 1995). The ice margin subsequently expanded into southern Ontario and ice-contact glaciolacustrine conditions, and fluctuating water depths are recorded in paleovalleys during the Early- and Mid-Wisconsinan (MIS 4-3: DS IV, see below).

Glaciolacustrine deposits of DS III-IV were truncated unconformably (and removed across large areas) when ice advanced over the entire region during the LGM. This event left widespread and thick till deposits (DS V; Figs. 5 and 6) whose surface is extensively drumlinized. In contrast, the two overlying depositional systems are of similar sedimentological and geomorphic characteristics (DS VI and DS VII) but are of different ages being the product of two episodes of ice streaming when the highly lobate and dynamic ice sheet margin, now melting rapidly, retreated eastward from the Erie Basin into the Ontario basin and northwards out of the Huron Basin. The beds of these ice stream lobes have been mapped and visualized in detail using analysis of LiDAR data. This reveals statistically validated flow sets of highly elongated streamlined subglacial bedforms called mega-scale glacial lineations (MSGLs; Clark 1993) that completely or partially reshaped older LGM till surfaces. Most significantly, ice streams and their meltwater systems delivered substantial volumes of glaciofluvial sediment to large, high-volume interlobate “kame” moraine systems that dominate the glacial geomorphology of southern Ontario and are of strategic hydrogeological significance (e.g., Waterloo and Oak Ridges moraines; Fig. 1).

A final phase of lateglacial sedimentation across southern Ontario saw the formation of large proglacial lakes (such as Whittlesey, Algonquin, and Iroquois in the Erie, Huron, and Ontario basins, respectively; see Barnett 1979, 1992). These left recessional moraines and locally, an uppermost cap of fine-grained glaciolacustrine deposits such as low relief clay and sand plains, iceberg scours and uplifted beach systems, terraces, shoreline bluffs, and wind-blown dune fields.

Illinoian and Sangamonian (DS I–DS II)

Illinoian till of MIS 6 (York Till in the GTA and the Bradtville Till in the Lake Erie basin; Dreimanis and Karrow 1972) is the oldest recognized subglacial depositional system in southern Ontario (DS I; Fig. 6). Another possible candidate is the Canning Till, which underlies the Waterloo area west of the Niagara Escarpment (Bajc et al. 2014). At Toronto, a bouldery lag deposit resting directly on shale bedrock, possibly York Till reworked in a shallow lacustrine environment, was formerly exposed in a now closed brick pit near Toronto in the Don Valley and in the Rouge River to the east. Stratigraphically equivalent till is locally encountered on bedrock along the flanks of the Laurentian paleovalley at surface (Mulligan 2017) and buried beneath younger sediments within the paleovalley (Mulligan et al. 2018; Sharpe et al. 2018; Bajc et al. 2019). In the central and northern parts of the Laurentian paleovalley east of the Niagara Escarpment, DS I contains two till units; a coarse-grained lowermost till with abundant Shield clasts, and an upper fine-grained till enriched in Paleozoic shale and carbonate lithologies, locally separated by glaciolacustrine deposits (Mulligan 2017; Bajc et al. 2019). Illinoian tills are also reported west of the Niagara Escarpment (Bajc and Dodge 2011; Bajc et al. 2015; Burt 2018) and at the base of rapidly eroding coastal bluffs on the north shore of Lake Erie (Dreimanis 1992). DS I is characterized by an upper bounding surface that is deeply weathered (Dreimanis 1992; Burt 2018; Mulligan and Bajc 2018).

Succeeding warm-to-cool climate interglacial lacustrine sediments (Don Formation: DS II; Fig. 6) are preserved along the Laurentian paleovalley. Study of exposures near Toronto at the southern end of the paleovalley shows that sediments accumulated in shallow water (<10 m) on the wave-influenced shoreface of an expanded ancestral Lake Ontario, named Lake Coleman (see Occhietti et al. 2016 for history of studies). The inferred water plane is at least 10 m above the modern level (75 m a.s.l.) possibly recording complete postglacial glacio-isostatic recovery of the bedrock sill at the eastern end of an interglacial lake near Kingston (the Duck–Galloo Ridge; Fig. 3). Conversely, subaerial weathering of the upper surface of DS II deposits in the northern Laurentian paleovalley, and where exposed along the north shore of Lake Erie, suggests a complex history of water level fluctuations. The uppermost part of the Don Formation at Toronto passes into deeper water basinal facies with a diverse pollen, plant macrofossil, and insect assemblage that identifies regional climate cooling and a corresponding increase in water depths in the ancestral Ontario basin (Lake Scarborough; Williams and Eyles 1995). This rise in lake level is consistent with the onset of regional glacio-isostatic depression and (or) blockage of the St. Lawrence River and the Rome outlet into the Hudson Valley by an expanding LIS (see Occhietti et al. 2015, 2016 for a discussion of regional correlations in southern Quebec). Water depths then rose to at least 40 m above the modern level of Lake Ontario recorded by deposition of an extensive boreal deltaic depositional system along the Laurentian paleovalley (Scarborough Formation; Early Wisconsinan: DS III: Fig. 6; see below). Berger and Eyles (1994) determined thermoluminescence ages of 80 ± 19 ka for the basal Don Formation and 67 ± 9 ka for its upper part just below the overlying cold climate Scarborough Formation. Occhietti et al. (2015, 2016) provided revised ages that are considerably older (>100 ka) and suggested that final Don Formation sedimentation marks regional cooling at the end of MIS 5e. At this time, the early LIS was large enough to block the St. Lawrence River, exert significant regional glacio-isostatic depression, and impound a deep lake(s) to the west in the Ontario and Erie basins.

Early Wisconsinan DS III

Outcrops of the cold climate deep water Early Wisconsinan Scarborough Formation (DS III; Fig. 6) were formerly well exposed along Scarborough Bluffs on the Lake Ontario shoreline east of Toronto. These are now largely obscured by engineering works designed to prevent erosion of the heavily urbanized coastline. The Scarborough Formation is ∼ 50 m thick and broadly divisible into a lowermost fine-grained pro-delta sequence dominated by varve-like couplets of silt and clay interbedded with muddy debris flows (Scarborough Clays), overlain by a delta-front and delta-top fluvial sequence (Scarborough Sands) in which subaerially formed ice-wedge casts indicate a local maximum water plane elevation of at least 40 m higher than the modern lake with mean annual temperatures of −6 °C or lower (Kelly and Martini 1986; Westgate et al. 1999). The overall paleogeographic setting likely consisted of a large glacier-fed sandy braided river system flowing south along the Laurentian paleovalley from Georgian Bay, discharging into an ice-dammed ancestral Lake Ontario. Significant glacio-isostatic depression is indicated, suggesting that LIS was sufficiently extensive to block the Rome outlet of the Ontario basin that drains to the Hudson Valley (see Fig. 1). The distinctive coarsening-upward character of the deltaic Scarborough Formation is recognizable inland on subsurface downhole geophysical data (Fligg 1983; Bajc et al. 2019) where continuously cored boreholes record an evolving fluvial system (Mulligan and Bajc 2018).

Early–Middle Wisconsinan DS IV

Sandy braid plain fluvial deposits of the uppermost Scarborough Formation (Scarborough Sands: Fig. 4) are blanketed by a unit of fine-grained predominantly deep water glaciolacustrine pebbly muds (Sunnybrook Drift) and the overlying glaciolacustrine Thorncliffe Formation (DS IV). This requires a substantial increase in water depths in the Ontario basin indicating damming of the St. Lawrence (Occhietti et al. 2011). Deep drilling data from the Erie, Huron, and Ontario basins show that DS IV is not restricted to the Ontario basin (e.g., Tyrconnell Formation in the Lake Erie basin; Dreimanis 1992; Fig. 5). Possible water depth changes may be recorded at the Scarborough Bluffs by the thin gravelly Pottery Road Formation present in the base of a large channel cut into the Scarborough Formation. Its wider regional paleoenvironmental significance remains unclear.

DS IV is found in all paleovalleys testifying to long-term influxes of fine sediment to a deep ice frontal lake(s) that spanned multiple Great Lake basins (Figs. 5 and 6). Much debate has focussed hitherto on the origin of fine-grained glaciolacustrine pebbly muds (diamicts) exposed at Scarborough such as the Sunnybrook Drift and other such facies (the overlying so-called Meadowcliffe and Seminary “tills” within the Thorncliffe Formation; Eyles and Westgate 1987; Hicock and Dreimanis 1992; Eden and Eyles 2001, 2002; Eyles et al. 2005). Diamicts are crudely bedded massive and laminated pebbly (and also stone-free) silty-clay (mud) facies typical of subaqueously deposited “rain-out” and debris flow deposits and contain in situ ostracod faunas (see below). These facies are conformably interbedded with laminated silty clays and ripple cross-laminated shoreface sands. Prior to modern sedimentological approaches to describing glacial sediments, pebbly mud diamict facies were interpreted as subglacial tills interbedded with “recessional” glaciolacustrine facies. Considered together they were seen as a record of successive advances and retreats of the LIS across the Great Lakes and well into mid-continent during MIS 4 and MIS 3 (Goldthwait and Dreimanis 1973). They were once employed as stratotypes defining a succession of cold climate stades and warmer interstadial climates named after Toronto-area localities (e.g., Sunnybrook Stade and Guildwood Stade; Dreimanis and Karrow 1972; Karrow et al. 2000).

Detailed facies and basin analysis studies of DS IV beginning in the early 1980s revealed a more complex paleoenvironmental setting for so-called subglacial “tills”. These identified deposition in deep water adjacent to a grounded LIS ice margin to the east involving the rainout of mud from suspension and ice rafting of clasts to produce pebbly muds, with water level changes, large storms, and slumping and mass flow of muddy rain-out facies, and scouring of bottom sediments by icebergs or partially floating ice (see Hicock and Dreimanis 1992; Eyles et al. 2005 and many references therein). A low diversity in situ benthic ostracod fauna occurs within the Sunnybrook and Thorncliffe Formation and is dominated by deep water species Candona subtriangulata and Candona caudata (Rutka and Eyles 1989). Schwarcz and Eyles (1991) determined the oxygen and carbon isotopic compositions of these two species ranging from the top of the Scarborough Formation through overlying DS IV sediments of the Thorncliffe Formation including pebbly mud diamicts (Sunnybrook, Meadowcliffe) to the base of the overlying LGM till (DS V; Fig. 6C). This yields a picture of a complex evolving ice-dammed waterbody. There is a clear upward decrease in δ18O in both species that is more marked in C. caudata. The shift in δ18O of ostracode calcite from interglacial values of −6 to −8‰ to values of −17 to −21‰ immediately below the LGM till (Fig. 6C) is the result of lowering of δ18O of lake waters by pulses of isotopically light meltwater from an advancing and thickening LIS, combined with lowering of δ18O of local precipitation due to falling temperatures. In the uppermost glaciolacustrine deposits immediately below the overlying LGM till (DS V), a δ18O value −17 to −21‰, indicates an ice-contact lake composed of between 35% and 55% glacial meltwater. Differences in δ13C values between C. caudata and C. subtriangulata increase upward within DS IV at Scarborough recording enhanced partitioning of carbon isotopes between surface and bottom waters as a result of increased water depths, photosynthesis by algae, or changes in the input of dissolved organic carbon from the oxidation of organic matter (Fig. 6C). Collectively, isotopic data provide valuable information regarding the evolving isotopic composition of lake waters during the early (MIS 4) and middle (MIS 3) phases of the last glacial cycle recorded at Scarborough. When combined with subsurface data from paleovalleys showing the regional extent of DS IV, there is strong evidence that southern Ontario was only partially glaciated at this time (Fig. 7). Deep drilling of paleovalleys across southern Ontario indicates the presence of glaciolacustrine DS IV beyond the Laurentian paleovalley (Mulligan 2017; Mulligan and Bajc 2018; Sharpe et al. 2018: Figs. 5 and 6) extending west of the Niagara Escarpment into the Erie and Huron basins.

The margin of the LIS in the Ontario basin when DS IV was accumulating in large ice frontal lake(s) is not precisely known. Correlative subglacially deposited tills have been reported from exposures along coastal cliffs 80 km to the east of Scarborough, near Port Hope (Martini and Acton 1975; Brookfield et al. 1982; Martini and Brookfield 1995). Karig and Miller (2013) also report subglacial tills inferred to be Middle Wisconsinan age from west central New York State. Kozlowski et al. (2014, 2019, 2021) and Porreca et al. 2018) describe multiple tills separated by glaciolacustrine and fluvial facies thought to span MIS 3-4 in the Finger Lakes. This suggests ice did not extend beyond the western end of the Ontario basin, possibly as a result of topographic confinement by the Niagara Escarpment but was able to expand southwards into upper New York State, where it was sufficiently thick to block the Rome outlet of the basin ponding deep water across southern Ontario (Fig. 7).

The absence of a regionally extensive subglacially deposited till(s) or a marked regional unconformity cut by glacial erosion, between DS II (the last interglacial) and DS IV (the Late Wisconsinan; Figs. 5 and 6) in the Laurentian paleovalley and in other paleovalleys is highly significant. This points to a prolonged phase of Early–Middle Wisconsinan glaciolacustrine sedimentation (Figs. 5 and 6) in a regionally extensive and likely complex glaciolacustrine proglacial water body. At some stage in this lengthy interval of regional flooding at the margin of the LIS, it can be speculated that water levels may have been as much as several hundred metres above the current level of Lake Ontario (c. 75 m a.s.l.). This suggestion rests on estimates of water depths from the highest elevation reached by lake floor deposits (521 m a.s.l.; Figs. 5 and 6) in addition to the substantial water depths indicated by the bathymetric preferences of ostracod faunas in these sediments (Rutka and Eyles 1989). Certainly, there were substantial fluctuations in water depths and thus lake size; at Scarborough, for example, DS IV contains intervals of storm-influenced shoreface sands characterized by hummocky and swaley cross stratification, with deformation structures and ice-rafted large boulders attributed to grounding icebergs (Eyles and Clark 1986; Eden and Eyles 2001,2002).

The record of limited ice extent and extensive glaciolacustrine conditions in southern Ontario during MIS 4 to 3 agrees with data elsewhere in the Great Lakes basin where the last few decades have seen a major re-evaluation of the pre-LGM extent of the LIS in mid-continent (e.g., Szabo 1997; Dyke et al. 2002; Szabo and Chanda 2004). In this respect, Mickelson and Colgan (2003) concluded that tills “formerly thought to be early Wisconsinan have recently been reinterpreted as Illinoian in many key locations in Illinois, Indiana, and Ohio”. Szabo (1997) similarly stated that “the literature suggests that nonglacial surficial processes may be inferred from sites in Ohio of early through middle Wisconsinan age”. In southern Ontario, the Canning Till was originally interpreted as early Wisconsinan, (e.g., Dreimanis and Karrow 1972), but recent drilling and sedimentological study suggest that the till may be of Illinoian age (DS I; Bajc and Dodge 2011; Bajc et al. 2015). In a review of the findings of deep drilling across southern Ontario, Burt (2018) concluded that “widespread nonglacial conditions spanned the Early and Middle Wisconsin”.

Last Glacial Maximum (DS V)

The LIS had sufficient volume to expand south and southwestward across the entire Great Lakes region reaching its maximum extent by 24 000 years B.P. (ybp) (Heath et al. 2018; Figs. 6 and 8). In southern Ontario, a prominent region-wide high relief unconformity was cut subglacially across older glaciolacustrine sediments of DS IV and covered by thick coarse-grained and frequently bouldery tills of DS V (variably named Northern, Newmarket, Bowmanville, or Catfish Creek tills; Figs. 5 and 6). Deep drilling around the Toronto area for a long-term waste disposal site shows these commonly contain glaciotectonically derived rafts of underlying glaciolacustrine sediment of DS IV near their base, testifying to widespread erosion and cannibalization of older glaciolacustrine deposits (Northern Till; Boyce et al. 1995). Tills were formerly regarded as homogenous massive units in early geoengineering and hydrogeological studies (see Gerber 1998). In contrast, a wealth of drilling and downhole geophysics data has allowed a detailed “architectural element analysis” (e.g., Miall 1985) of the three-dimensional structure of these tills, which reveals a more complex composite internal stratigraphy. This consists of stacked till “elements” up to several metres in thickness, recording punctuated aggradation by deforming subglacial debris being moved below the ice sheet base (Boyce and Eyles 2000; Meriano and Eyles 2009). Individual till elements are bounded by disconformities with boulder pavements and thin (<1 m) interbeds of sorted sediment, likely marking short-lived local phases of erosion across the ice sheet bed during overall subglacial aggradation (Fig. 6B). The orientation of drumlins on the surface of LGM tills indicate predominantly south to southwest-directed regional flow across southern Ontario (Fig. 8). Additional examination of the internal structure and lithology of LGM tills may shed light on evolving regional ice flow directions immediately prior to the LGM.

Deglaciation after 17 400 ybp: (DS VI)

During final phases of Late Wisconsinan glaciation after ∼ 17 400 ybp, the ice sheet margin had thinned significantly and become highly lobate as a consequence of strong topographic control on ice flow by the Niagara Escarpment, and Huron and Erie–Ontario basins. Much prior work has focussed on the evolution of large complex ice-contact glacial lakes that formed as lobes retreated (e.g., Barnett 1992).

Access to newly available high resolution LiDAR topographic digital datasets (Ontario Ministry of Natural Resources and Forestry Ontario 2018) and novel processing and visualization techniques developed at the University of Toronto (summarized in Fig. 9) allow quantitative analysis of the geomorphological record of the final stages of the LIS in southern Ontario. This is especially the case in mapping drumlins and their morphological variants such as narrow, highly elongate bedforms as MSGLs (Clark 1993) that are produced under ice streams where surface ice velocities were likely in excess of 1 km year−1. MSGLs are subtle, low relief landforms (previously mapped by many studies as “flutes”) and their full extent in general across southern Ontario and upper New York State is only now becoming evident from quantitative analysis of LiDAR data. This reveals statistically validated “flow sets” of MSGLs defining individual paleo ice streams that were part of broader ice lobes flowing along each basin (Fig. 10).

Two such episodes of fast flow can be broadly recognized across southern Ontario; what have been named the Port Bruce (Figs. 11 and 12) and Port Huron stages (Karrow et al. 2000; Fig. 13). During both these intervals, substantial volumes of sediment and rock were eroded and entrained by energetic fast-flowing ice streams (and their associated meltwater systems; Figs. 14 and 15) and deposited as large interlobate moraine complexes (e.g., Waterloo, Orangeville, and Oak Ridges; Russell et al. 2004; Sharpe et al. 2007; Bajc et al. 2014; Burt and Dodge 2016; Burt 2018).

Port Bruce stage of fast flow

Rapid construction of the interlobate Waterloo Moraine occurred during the Port Bruce phase just after ∼17 400 ybp (Weaver and Arnaud 2011) when ice lobes converged from the Georgian Bay and Ontario basins (Figs. 11 and 16). Their streamlined beds are defined by longitudinal corridors (flow sets) of MSGLs that are underlain by thin “upper” tills of highly varying texture. These flow sets terminate downglacier at large interlobate moraines or prominent till-cored moraines (Figs. 11 and 12). By ∼16 900 ybp, the southern LIS margin may have been located at the Valley Heads Moraine in New York State, at the southern end of the glacially overdeepened Finger Lake basins (Fig. 10; Bloom 2018; Young et al. 2020). Sookhan et al. (2021) documented the effects of large deep ice-contact lakes in New York State, triggering late-stage fast flow when LGM drumlins were partially eroded and overprinted by flow sets of highly elongated MSGLs resulting in palimpsest landforms.

Dalton et al. (2023) noted a lack of control on the position of the retreating LIS margin in the eastern Great Lakes between the LGM and c. 17 500 ybp. A detailed analysis of successive ice margin positions marked by moraines of various types, and thus inter-regional correlation of ice stream margins, is beyond the scope of the present paper and awaits further geomorphic mapping using LiDAR. In brief, at this point in time, the Port Bruce stage is generally correlated with the Lake Escarpment Moraine in Pennsylvania (Braun 2004) and the Defiance Moraine in Ohio (Eschman 1985). Other possible correlative positions of the LIS margin at this time have been suggested by Karrow (2004), Ridge et al. (2012), Franzi et al. (2016), Porreca et al. (2018), Kozlowski et al. (2018, 2019, 2021), and Young et al. (2020). Significantly, Ross et al. (2006) recognized a synchronous change in ice sheet structure at about 17 000 ybp along the upper St. Lawrence Valley. This was interpreted as marking the onset of a large St. Lawrence Valley Ice Stream flowing east to the Gulf of St. Lawrence (e.g., Parent and Occhietti 1999; Occhietti et al. 2001, 2011). This event, and the succeeding Port Huron event of fast flow (see below) can be tentatively related to the abrupt onset of the short-lived Bølling–Allerød warming event in Europe sometime shortly after 17 000 ybp (Rasmussen et al. 2006), which is associated with Meltwater Pulse 1A when global sea level rose ∼20 m in less than 500 years (Deschamps et al. 2012). While we recognize increasing evidence for diachronous regional climatic shifts at this time (see Lowe et al. 2008), we note that the model of Tarasov et al. (2012) for the deglaciation of the LIS indicated that the interval of fast flow recognized in southern Ontario and adjacent New York State coincides with maximum ice loss in North America (see also Norris et al. 2022). The term “ice sheet collapse” is appropriate (Mullins and Eyles 1996; Deschamps et al. 2012).

Port Huron stage of fast flow

About 1000 years after the Port Bruce stage, a second (and final) phase of fast flow (the Port Huron stage) affected southern Ontario. Large volumes of sediment were eroded and remobilized by ice streams, which resulted in formation of the largest glacial landform in southern Ontario, the ORM. This feature forms a major topographic high extending some 160 km east of the Niagara Escarpment. As the last interlobate moraine complex to be deposited during regional deglaciation of southern Ontario, it is a befitting finale to the Wisconsinan glaciation. The short-lived Port Huron phase ended some several hundred years (or less) later, with final ice retreat and the flooding of the newly deglaciated Ontario basin by glacial Lake Iroquois, and its equivalent waterbody glacial Lake Algonquin in the Huron basin sometime after c. 14 500 ybp (Dalton et al. 2020, 2022,2023).

The subsurface geology of the ORM is dominated by subaqueously deposited sand and gravels that accumulated between converging ice streams. Flow sets of MSGLs underlain by thin (<5 m) ice stream tills (e.g., Halton Till) cut across much thicker underlying LGM till, converging on the ORM from north and south (Sookhan et al. 2018; Eyles et al. 2018; Figs 13 and 14). Recent LiDAR mapping of flow sets of MSGLs using the processing methodology outlined in Fig. 9 now permits expansion of the number and dimensions of converging ice streams (Fig. 10). The floors of these ice streams were lowered by subglacial erosion under fast-flowing ice and MSGLs cut into older drumlinized tills. High standing drumlinized LGM till surfaces define inter-ice stream areas. LiDAR images capture a range of “palimpsest” subglacial bedforms recording modification of the original drumlins form (Fig. 17). For example, LGM drumlins on the Northern Till on the south flank of the ORM were partially modified by later northwest ice flow of the Halton Ice Stream producing distinctive modified drumlin forms that were briefly noted by White (1985) and Karrow (1967). LiDAR data from these areas reveal subtle flow sets of MSGLs imprinted on LGM drumlins, oriented oblique to the original long axis drumlin trend (Fig. 17).

The lateral margins of the south-flowing Simcoe, Scugog, and Rice Lake ice streams are marked by prominent triangular shaped beads (the so-called “wedges” of Barnett et al. 1998) where the width of the ORM abruptly increases. These beads are the morphological expressions of large fan-delta bodies deposited in a complex interlobate water body trapped between converging ice streams (Figs. 10 and 13). Studies of modern ice streams have shown that the location of subglacial tunnels (and thus eskers) is structurally controlled; tunnels can only be kept open under more slowly flowing inter-ice stream areas and are ephemeral features below highly dynamic ice streams (e.g., Dowdeswell et al. 2015). A full glacio-sedimentological treatment of the ORM combining the results of new LiDAR mapping with subsurface data is being presented in full elsewhere. In brief, each fan delta body shows a broadly consistent tripartite internal structure of a lowermost core of coarse gravel resting on drumlinized LGM till (Northern Till) and likely deposited in subglacial conduits as eskers (Barnett et al. 1998; Russell et al. 2004; Sharpe et al. 2004, 2018; Mulligan et al. 2018; Sharpe 2022). These conduits fed a composite west–east-oriented subglacial esker system when ice streams were fully confluent (Fig. 13). At this stage, subglacial meltwaters flowed westward across the Niagara Escarpment into the Erie Basin. In turn, basal gravels are overlain by large discrete fan delta bodies that accumulated sequentially and diachronously from west to east, as the margins of the northern ice streams began to progressively “unzip” and withdraw from the Halton Ice Stream, which remained in the Ontario basin ponding a complex water body between the now divergent lobes.

LiDAR-based geomorphic mapping is ongoing across the ORM to identify successive ice marginal positions during final deglaciation, focussing on push moraine ridges noted by Barnett (1996) and Mulligan et al. (2019) and the orientation of the many small moraine ridges, newly identified on LiDAR data from both flanks of the ORM interpreted as De Geer moraines (Fig. 17). These indicate that fast ice flow had likely ended and the ice margin was retreating incrementally year by year in shallow water; the spacing of ridges suggests annual retreat rates of between a few metres to ∼100 m (see also Zilliacus 1989; Ottesen and Dowdeswell 2006; Ojala et al. 2015; Bukhari et al. 2021). In contrast, a calving bay margin is recorded by iceberg scours in deeper water of glacial Lake Iroquois, which flooded the Ontario basin as ice retreated eastward toward the St. Lawrence Valley (Fig. 18D).

The St. Lawrence Ice Stream had collapsed and fully retreated by c.13 500–12 800 ybp (Parent and Occhietti 1999; Occhietti et al. 2001; Dalton et al. 2020). This marks the end of ice streaming in southern Ontario, the final drainage of glacial Lake Iroquois and deglaciation of the study area (Dyke 2004; Donnelly et al. 2005).

Several themes emerge from this paper that we briefly explore here. In regard to describing the complex and long-studied late Pleistocene glacial stratigraphy of southern Ontario, a depositional system approach brings a regional coherency to existing stratigraphic nomenclature. This approach also furthers understanding of the paleoenvironmental history and paleoglaciology of the LIS, especially during the MIS 3-4 interval. It provides a firmer basis for applied geotechnical and hydrogeological investigations reliant on subsurface data. Individual depositional systems can be readily identified by nonglacial specialists by reference to the principal facies types seen in drill core, in outcrops or inferred from geophysical studies, and the broader stratigraphic “context” provided by overlying and underlying sediments. In particular, coarse-grained, heavily overconsolidated, bouldery LGM tills (DS V) resting on a marked regional-scale unconformity, with drumlinized surfaces, form a regionally important marker horizon. Equally identifiable is the coarsening-upward lithology of the deltaic Scarborough Formation (DS III) along the Laurentian paleovalley, and the overlying thick fine-grained glaciolacustrine package of the Thorncliffe Formation with its distinctive pebbly mud facies (DS IV) preserved in the Laurentian paleovalley and now being traced to other basins (Fig. 6).

Dalton et al. (2020) suggested that ice sheet reconstructions for the last glaciation would benefit from detailed studies in southern Ontario. The current paper finds that the LIS expanded during the LGM after 24 000 ybp across a thick and regionally extensive succession of glaciolacustrine deposits. This indicates that ice did not extend westward beyond the Lake Ontario basin and the Niagara Escarpment before that time. Consequently, there is no compelling geological evidence across most of southern Ontario outside the Ontario basin for pre-LGM tills between the Illinoian glaciation (MIS 6) to the LGM. Based on outcrop and subsurface data showing the wide extent of glaciolacustrine deposits of DS IV (Figs. 5 and 6), it appears that for much of the early and middle phases of the Wisconsinan glaciation, the margin of the LIS terminated in large deep ice-contact lakes, at times contiguous across parts of the Ontario, Erie, and Huron lake basins, other low-lying portions of southern Ontario and into upper New York State (Fig. 7). This large water body likely experienced marked changes in water depths and abrupt inputs of isotopically light meltwaters indicating a dynamic ice sheet (Fig. 6).

The history of the LIS in southern Ontario outlined above is consistent with the “rapid late growth” model proposed by Weisenberg et al. (2017) and others (e.g., Carlson et al. 2018). This argues that the ice sheet was of restricted volume and configuration during its early and middle phases before growing very rapidly over the course of a 15 000 year interval, immediately prior to the LGM. Extensive ice-frontal lakes several hundred metres deep within the glacio-isostatically lowered Ontario, Erie, and Huron basins may have acted as a large moat, presenting a significant obstacle to further ice expansion during the earlier and middle phases of the ice sheet’s history. The role of the overdeepened Rochester Basin (Fig. 3) in steering ice flow into New York State while parts of southern Ontario west of the Niagara Escarpment remained ice free should be examined further. As a final comment, we stress that while subsurface and outcrop geological data provide a consistent picture of the stratigraphic record of an expanding ice sheet following its inception, it has to be recognized that age dating of pre LGM events in southern Ontario still remains poorly constrained (e.g., Berger and Eyles 1994; Occhietti et al. 2015,2016).

The record of the final deglacial phases of sedimentation in southern Ontario is recorded geomorphologically by well-defined flow sets of MSGLs and thin tills (Fig. 10). These record two phases of ice streaming (Port Bruce and Port Huron phases around 18 000–16 980 and 16 000–14 300 ybp, respectively). Overall, these phases were very short-lived events within the broader history of Wisconsinan glaciation in southern Ontario but are remarkable for the very large volumes of sediment and bedrock debris that was eroded and moved by fast-flowing ice (e.g., Waterloo and Oak Ridges moraines). Eyles (2012), Krabbendam et al. (2016), Bukhari et al. (2021), and Mulligan et al. (2023) all highlighted rapid quarrying and shaping of up-ice facing escarpments of Paleozoic carbonates by fast-flowing ice. This has resulted in a wide extent of grooved bedrock surfaces (Figs. 18A8C) and extensive spreads of hummocky and bouldery carbonate-rich tills (e.g., Dummer Moraine; Fig. 4). Bouldery carbonate-rich debris covers a large area of the dip slope of the Niagara Escarpment (Cowan 1978; see fig. 13 in Karrow 1987) and is present in the cores of several moraine ridges (Paris, Galt, and Moffat moraines; Arnaud et al. 2018) and multiple moraines along the Niagara Peninsula (Crystal Beach, Fort Erie, Niagara Falls, Vinemount, and Waterdown Moraines; Burwasser 1979; Maclachlan and Eyles 2013; Burt and Mulligan 2017; Barnett and Karrow 2018; Burt 2020). This highlights, once again, the profound geomorphic influence of short-lived phases of ice streaming and fast flow in southern Ontario.

A key emerging element in the study of glaciated terrains is the availability of LiDAR-derived topographic data, which allows detailed visualization of landscapes and bedforms. The imagery used in this study relies on meticulous processing steps (Fig. 9) to achieve the optimum level of visual detail (Figs. 14 and 18). The results of machine learning analysis of LiDAR data reveal statistically unique assemblages of subglacial bedforms based on orientation data (flow sets; Fig. 13A) in southern Ontario. This identifies that much of the geomorphology of southern Ontario owes its origin in substantial part, to short-lived phases of ice stream flow during final deglaciation. In this regard, these glacial landscapes can be classified as “ice stream landsystems” (e.g., Stokes 2011; Evans et al. 2014; Klages et al. 2016) where drumlinized and usually heavily overconsolidated LGM tills were reconfigured by later fast flow (see also Hättestrand et al. 1999; Kleman et al. 1999; Putkinen et al. 2017; Sookhan et al. 2022; Eyles et al. 2016, 2023).

A Chautauqua Ice Stream has very recently been identified as part of ongoing LiDAR mapping during the course of this study (Figs. 10 and 16). This flowed south from the Erie basin, most likely during the Port Bruce Stadial immediately after 17 000 ybp. This discovery underscores the need for a comprehensive LiDAR-based mapping of the entire Great Lakes basin. This exercise will ultimately generate a comprehensive understanding of the distribution of fast-flowing lobes and their constituent ice streams, their timing and likely triggering mechanisms in regard to evolving ice frontal lakes and underlying topography (e.g., Sookhan et al. 2022). There is already (Fig. 10) a much more complex picture of the extent of ice streaming within the eastern Great Lakes margin of the LIS during deglaciation than previously documented (Tarasov et al. 2012; Margold et al. 2015, 2018; Stokes et al. 2016). Such a study would address the lack of understanding of the retreat phases of the LIS after LGM highlighted by Dalton et al. (2023).

We are confident that by detailed mapping of subtle topographic features, such as recessional moraines and newly recognized De Geer moraines (Figs. 15, 17, and 19), a clear picture will emerge of the stages of ice marginal retreat that will be well beyond the current resolving power of radiocarbon dating. Moreover, using statistical analysis of landform elongation data as a proxy for ice flow velocity (e.g., Eyles et al. 2023; Fig. 15), we will be able to draw inferences about the depositional conditions under which LGM tills were deposited, as well as the evolution of widespread drumlin bedforms, including their later modification by fast flow (Fig. 17). Superposition of MSGLs on older bedforms is typically associated with so-called “flow switching” where thin, fast-moving, and highly dynamic ice streams abruptly change their direction and partially reshape their beds (see Dowdeswell et al. 2006; Greenwood et al. 2012; Winsborrow et al. 2012). LiDAR mapping of swaths of De Geer moraines on both flanks of the ORM explains enigmatic subtle changes in relief noted previously from analysis of air photographs (e.g., Karrow 1991). These ridges present a significant opportunity to reconstruct in detail the pattern and rate of ice recession and grounding line dynamics of thin (<several hundred metres?) ice streams (e.g., Lindén and Möller 2005).

Modern ice streams in Antarctica are, in some cases, topographically controlled by basement structures in the underlying Shield (see Bell et al. 1998; Campo et al. 2017). In this regard, we note a possible geographic relationship between ice streams north of the ORM and lake basins in Paleozoic bedrock, with underlying Precambrian terrane boundaries (e.g., Figs. 2 and 13). This suggests the subtle structural–topographic steering of ice flow by basement structures such as identified in eastern Ontario by Gadd (1980) where ice flow was guided into the Ottawa Valley and thence into the Ontario basin, by the prominent Madawaska Fault on the eastern flank of the Algonquin Highlands. Further mapping of topography and structure using LiDAR can fully resolve this possible relationship.

A depositional system approach rationalizes the complex subsurface glacial stratigraphy of southern Ontario, preferentially preserved in several large paleovalleys that were protected in part, from later LGM erosion. The approach not only aids reconstruction of paleoenvironments but also provides a simple, practical, and above all regionally coherent framework for applied geoengineering and hydrogeological investigations.

Documentation of depositional systems addresses significant gaps in knowledge of the paleoenvironmental history of the LIS in southern Ontario. New subsurface data from paleovalleys confirm previous interpretations based largely on outcrop studies, of restricted LIS extent in the Huron, Erie, and western part of the Ontario basins during the early and middle phases of the last glaciation (MIS 4-3). Ice was present in New York State of sufficient extent to block the Rome outlet of the Ontario basin, creating large ice-dammed lakes across the Ontario–Erie and Huron basins. LIS expanded rapidly westward across the Niagara Escarpment to its LGM limits by 24 000 ybp (MIS 2). Its deglacial phases after 17 400 ybp were dominated by paleo ice streams within lobes steered by the Ontario–Erie and Huron basins.

Analysis of LIDAR data provides crucial information regarding the presence and origins of palimpsest streamlined surfaces that traditionally have been poorly documented. It underscores the urgent need for expanded national LiDAR coverage to fully resolve the geomorphic complexities of complex glacial landscapes.

We are very grateful to Kirsten Kennedy, Shane Sookhan, Naki Akçar, and Mike Doughty for helpful discussions and assistance in the field and laboratory. We thank Andrew Koslowski and Etienne Brouard for their supportive critical comments on an earlier version of this paper. We are also particularly grateful to Jessy Rice for editorial comments on the final manuscript.

Data used in this study are openly available online through Ontario GeoHub and OGSEarth libraries.

Conceptualization: SB, NE, RM, AB, CE, RP, MR, NP

Data curation: SB, NE, RM, AB, CE, RP, MR, NP

Formal analysis: SB, NE, RM, AB, CE, RP, MR, NP

Funding acquisition: NE, RP

Investigation: SB, NE, RM, AB, CE, RP, MR, NP

Methodology: SB, NE, RM, AB, CE, RP, MR, NP

Project administration: SB, NE

Resources: SB, NE, RM, AB, CE, RP, MR, NP

Software: SB, NE

Supervision: SB, NE

Validation: SB, NE, RM, AB, CE, RP, MR, NP

Visualization: SB, NE, CE

Writing – original draft: SB, NE, RM, AB, CE, RP, MR, NP

Writing – review & editing: SB, NE, RM, AB, CE, RP, MR, NP

This investigation and compilation was funded by a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada to NE. Additional support for this research was provided by the Geological Survey of Canada (GSC) Geo-Mapping for Energy and Minerals Program (GEM-GeoNorth) and is currently funded by the GSC Critical Mineral Geoscience and Data (CGMD) Program. This is Geological Survey of Canada Contribution # 20230122