Late history of glacial Lake Agassiz in northwestern Ontario, Canada: a case study in the Sandy Lake basin Open Access

The Sandy Lake basin in northwestern Ontario (Fig. 1) is a potentially key area in discussion of the late history of glacial Lake Agassiz because of its extensive glaciolacustrine deposits and well-preserved shoreline features of this geological episode (e.g., Elson 1967; Thorleifson 1996; Teller and Leverington 2004). However, aside from the general discussions, only cursory mapping has been conducted in this area primarily because of the lack of road access (Hurst 1929; Derry and MacKenzie 1931; Satterly 1939; Prest 1963). These previous surveys were confined to the shore of Sandy Lake and its immediate vicinity, whereas the Opasquia end moraine and the prominent shoreline features on it on the northern rim of the basin have received little study. A brief survey was conducted around this moraine for the creation of the Opasquia Provincial Park in the 1970s, but the full technical report from this survey has never been published (Harvey 1979). Consequently, the Sandy Lake basin has remained largely unknown as to its deglaciation history and records of Lake Agassiz until recently (Gao et al. 2017,2018).

It has been a continued interest in the geoscience community to gain insights into the final drainage of the coalesced Lake Agassiz–Ojibway (Lake Agassiz hereafter for simplicity) into Hudson Bay and its potential causal link to the 8.2 ka cooling event as recorded in the Greenland ice core (e.g., Barber et al. 1999; Teller et al. 2002; Roy et al. 2011; Jennings et al. 2015; Gauthier et al. 2020; Brouard et al. 2021). The Sandy Lake basin remained connected as an arm to Lake Agassiz until the lake drained into Hudson Bay (Fig. 1) (Leverington et al. 2002; Clarke et al. 2004; Teller and Leverington 2004). However, the strandlines of the lake prior to its final drainage, the Ponton and Kinojévis shorelines to the west and east, respectively, have never been adequately dated (e.g., Veillette 1994; Thorleifson 1996; Teller and Leverington 2004; Fisher 2020; Godbout et al. 2020). The only radiocarbon date available was on a freshwater molluscan shell from the Ponton beach near Lake Winnipeg, Manitoba (8310 ± 180 ¹⁴C years BP, GSC-1679), but the true age of this fossil remains unknown because of the potential hardwater effects that are difficult to determine (e.g., Lowdon et al. 1977; McMartin 2000; McNeely and Brennan 2005; Gauthier 2022). Instead, these beaches are inferred at 8.5 ka (Teller and Leverington 2004; Teller and Owen 2021), based on radiocarbon dates on marine shells associated with the Tyrrell Sea incursion in Hudson Bay Lowlands during the final drainage of the lake (Barber et al. 1999). Conversely, recent dating on marine and terrestrial fossil remains indicates the marine incursion occurred as young as at 8.16 cal ka (e.g., Roy et al. 2011; Jennings et al. 2015; Brouard et al. 2021). Further complications arise from the presence of the lower Fiddler and La Sarre shorelines in the Agassiz and Ojibway basins, respectively (Klassen 1986; Thorleifson 1996; Roy et al. 2015; Godbout et al. 2020). These lower strandlines, although not directly dated, are suggested to represent the lake level of the last discharge of the lake responsible for the 8.2 ka cooling event and the Ponton– Kinojévis shorelines represent an early partial drainage probably through subglacial discharge (Clarke et al. 2004; Breckenridge et al. 2012; Daubois et al. 2015; Roy et al. 2015; Godbout et al. 2019; Gauthier et al. 2020; Brouard et al. 2021). The lack of independent age constraints of these shorelines causes confusions about the timing and role of these two lake-level phases in lake drainage and ocean freshening that probably caused a reduction in Atlantic meridional overturning circulation and subsequent cooling.

Recent Quaternary geology mapping by the Ontario Geological Survey in the Sandy Lake basin enabled the writer to sample and date the glaciolacustrine deposits and abandoned shoreline features, including the Ponton shoreline (Gao et al. 2017, 2018). This paper presents the field observations and radiometric dating results from this mapping, which provides insights into the deglacial history of this basin and development of Lake Agassiz as to its inception and final withdrawal from this region. In the discussion below, all ages cited are calendar thousand years ka and those directly from radiocarbon dating are expressed as cal ka or cal years BP.

The Sandy Lake basin is located on the Precambrian Shield consisting of bedrock hills and ridges with intervening lakes and fen–bog wetlands in boreal forest, with a mean annual, July, and January temperature at −1.2, 17 and −23 °C, respectively, and an annual precipitation at 520 mm (1950–1980) (https://climatedata.ca). Sandy Lake is the lowest point of the basin (274 m asl), which drains eastwards via a narrow outlet into the Severn River that flows northeastwards across Hudson Bay Lowlands into Hudson Bay (Fig. 2). A distinct land mark in this area is the west to east aligned Opasquia end moraine that forms the northern rim of this basin (Fig. 2). The moraine, which rises at 409 m asl and has a local relief up to 100 m above the immediately surrounding area, is a local drainage divide separating the Sandy Lake basin from those to the north, e.g., the Island Lake basin to the northwest in Manitoba (Fig. 2). The low-relief Hudwin end moraine in Manitoba is its probable western extension (Fig. 2) (Klassen 1986; Clarke 1989; Gauthier et al. 2022). Adjacent to the Opasquia moraine to the east is the north to south aligned Sachigo interlobate moraine, which lies mostly outside the Sandy Lake basin to the north (Fig. 2). To the south, the basin rises gradually on Precambrian bedrock uplands where the Quaternary deposits are thin and discontinuous (Fig. 2). To the west, it is demarked by Hudwin Lake in Manitoba that drains on thick glaciolacustrine argillaceous sediments into Sandy Lake via the Cobham and upper Severn rivers (Derry and MacKenzie 1931; Clarke 1989; Gao et al. 2017,2018).

Glaciolacustrine deposits consisting predominantly of varved clay are common in the basin, particularly, in its center around Sandy Lake where such clay is extensively exposed along lakeshores. Presently, Sandy Lake undergoes intense shore erosion by lake waves, and consequently, it has a light mocha colour in summer because of the high concentration of suspension material and a cuspate shoreline typified with receding, steep clay banks or bluffs up to 6–8 m high between bedrock ridges and knolls (Fig. A1). There are several small lakes adjacent to Sandy Lake (Fig. 2). Among them, the circular-shaped Niska Lake with a diameter of 6 km is peculiar. This lake does not have regular incoming rivers and is recharged through groundwater and random precipitation runoff. Despite its higher elevation at 275 m asl and distance less than 1 km to Sandy Lake, Niska Lake does not directly connect to Sandy Lake. Instead, it drains via an eastern outlet along a small creek into the Dawson River that flows northwards into Sandy Lake (Fig. 2).

Access in the study area was via helicopter and a power boat. Field methods included section logging, photographing, and sampling in natural exposures along lakeshores and in man-made sections through hand augering and trenching. The surface elevations were estimated from Google Earth™ mapping service for lakes, abandoned shorelines, moraine ridges, and other geological features. In the study area, the Shuttle Radar Topographic Mission digital surface model used by Google Earth™ has a resolution of 30 m/pixel and a root mean square accuracy of 3.64 m in elevation (Gesch et al. 2014).

The surface peat was sampled for pollen analysis and radiocarbon dating. The resultant pollen diagram was included as an Appendix A figure (Fig. A2). The radiocarbon dates were calibrated to calendar years and expressed as cal years BP or cal ka before 1950 using IntCal20 calibration curve at OxCal v4.2.4 (Reimer et al. 2020; Bronk Ramsey 2009). Glaciolacustrine sand was sampled using copper tubes for optically stimulated luminescence (OSL) dating. OSL dating measures the time elapsed since the sediments were last exposed to sunlight and hence was taken to represent the burial age. The infrared stimulated luminescence (IRSL) method was used to measure quartz sand in the samples. The IRSL measurements were carried out with a single-aliquot regenerative dose protocol, and the central age model by Galbraith et al. (1999) was applied to obtain a weighted average age in all calculations and reported as ka. The detailed method was included as Supplementary material S1.

The isobases proposed by Teller and Thorleifson (1983) and Teller and Leverington (2004) and the generalized configuration of Lake Agassiz shorelines summarized by Thorleifson (1996) and refined later by McMartin (2000) were used for projected water planes and shoreline correlations. Recently, Lewis et al. (2022) compiled and re-analyzed some of the known beaches in the Agassiz basin (the Campbell) and concluded that, except for the northwesternmost part of the basin (in Saskatchewan) where the isobases bend toward the west (Rayburn and Teller 2007), the glacial isostatic tilts are very similar to the originally proposed for the Agassiz basin. McMartin’s (2000) work, which was built on Klassen’s (1986) and other previous studies on the late-phase shorelines north and northwest of Lake Winnipeg, has close relevance to the current study because of the presence of correlative shorelines in the Sandy Lake basin. The shorelines she studied are within isobases 7–9 (6.75–9.25). Because the Sandy Lake basin is located farther to the northeast on isobases, the strandlines on her paleo-elevation diagram were manually extrapolated to include isobase 10.5 (Fig. A3). The ages assigned to these late-phase shorelines as summarized by Teller and Leverington (2004) were used in the previous studies (Gao et al. 2017, 2018). Recently, Teller and Owen (2021) updated the assigned ages based on the new OSL dates they obtained (on the Gimli beach). Consequently, these updated ages were adopted in the discussion below.

Recently, Godbout et al. (2023) reconstructed the paleo-topography of North America for various time slices after the last glacial maximum based on the glacial isostatic deformation models (ICE–xG models) (e.g., Peltier et al. 2015; Roy and Peltier 2017). On the paleo-topography, the shoreline of a lake level formed at a particular time, e.g., 8.5 ka, is expected to have the same paleo-elevations (Godbout et al. 2020) and the release of these paleo-topographic maps would greatly aid regional delineation of the glacial lakes and correlation of the former shorelines. Because of the well-defined Ponton and Kinojévis shorelines as reported (e.g., Bell 1978; Klassen 1986; McMartin 2000; Leverington et al. 2002; Roy et al. 2015; Godbout et al. 2020), the Ponton shoreline features in the Sandy Lake basin were examined against the paleo-topography (8.5 ka) and the paleo-elevations extracted were compared with those on the eastern side of Lake Winnipeg in Manitoba and the Kinojévis-type shorelines in the La Sarre area in Quebec (see Fig. 1 for locations). While correlative among these sites, the paleo-elevations along a more extensive Ponton strandline that runs along the western shore of Lake Winnipeg and further to the north were inconsistent, which, against expectations, decrease systematically to the north (see Fig. 1 for location). Included in Appendix A is a 71 km segment of this south to north aligned Ponton strandline (Leverington et al. 2002; Gauthier and Keller 2020) where the paleo-elevations were extracted (Fig. A4). Furthermore, the two paleo-topography models derived from ICE6G and ICE7G (Godbout et al. 2023) were briefly compared based on the known shorelines. The results (not included) showed that the pale-elevations extracted from the former appear relatively more consistent than those from the latter on these Ponton shorelines. Specifically, on the paleo-topography model derived from ICE7G, the Ponton shorelines around Lake Winnipeg are substantially higher in elevation than the Ponton and Kinojévis-type shorelines in the Sandy Lake basin and the La Sarre area. These inconsistencies call for refinement of the models. Until these models are fully tested and adjusted, the isobases and the generalized configuration of Lake Agassiz shorelines that are based on field data (e.g., Teller and Thorleifson 1983; Thorleifson 1996; McMartin 2000) will remain as a useful method for clues in regional shoreline correlations in the Agassiz basin (e.g., Breckenridge 2015; Fisher and Breckenridge 2022).

The Quaternary deposits in the Sandy Lake basin consist of a Late Wisconsin till and overlying glaciofluvial and glaciolacustrine sediments below the surface peat (Gao et al. 2017, 2018). The till, which is sporadically preserved in bedrock depressions, is a noncalcareous, sandy, and stone-rich diamicton with gravel clasts dominated by granitic and metavolcanic lithologies of local origin, including rare, far-travelled carbonate, and greywacke clasts (“omars” pebbles) from the Hudson Bay region. The glaciofluvial deposits consist mainly of sand and gravels in moraines, kames, and eskers, which are commonly modified by lake waves, veneered with glaciolacustrine sand, and strewn with lag boulders (Figs. A5 and A6). The Opasquia and Sachio moraines contain the largest glaciofluvial sand and gravel deposits in this region.

The glaciolacustrine deposits consist of nearshore and beach or littoral sand and gravel and basinal silt and clay. The littoral deposits commonly occur as sand covering the ground surface, in abandoned beach ridges, and as lag boulders on moraines (Figs. A5 and A6) and on bedrock knolls or ridges resulting from lake wave erosion and winnowing of the former glacial and glaciofluvial gravelly deposits. Littoral deposits probably also occur below the basinal deposits as beach or nearshore sand and pebbles, with oversize lag cobbles and boulders draped by basinal fine sand and silt (varved clay) (Figs. 3A, A7, and A8) (Gao et al. 2017,2018).

The fine-grained basinal deposits are extensive and usually varved, referred to informally here as varved clay. The varved clay, up to 8 m thick, is well exposed along the shore of Sandy Lake (Satterly 1939; Prest 1963; Gao et al. 2017, 2018). Locally, it may reach to 12.5 m thick as indicted by the water well records at the Sandy Lake and Keewaywin communities (Gao et al. 2017). In a varve or couplet in the deposit, the summer layer is composed of light grey to pale calcareous silt and fine sand, whereas the overlying winter layer of dark brown clay, which is usually non-calcareous without reaction to 10% HCl solution (Fig. 3B) (Gao et al. 2017, 2018). The couplets, which are tens of centimeters thick at the base of the deposit, become thinner up section where they typically attain cm-scale thickness (Fig. 3B), interpreted as a result of ice recession and reduced sediment supply to the glacial lake (e.g., Palmer et al. 2019).

Abandoned lake shorelines are common on raised grounds, e.g., moraine and bedrock ridges, as well as in low-lying areas. They are in the form of beach ridges and wave-cut scarps or bluffs and commonly seen on the Opasquia and Sachigo moraines and in lowlands adjacent to Sandy Lake.

Several geological sections were sampled for radiometric dating, which include two sections in varved clay for OSL dating and three sections in surface peat for radiocarbon dating (see Fig. 2 for locations). Section 17-CG-014 consists of 1 m of varved clay underlain by 0.1 m of pebbly sand on 0.5 m of till in a bedrock depression exposed in a quarry on a bedrock ridge (286 m asl) at the Sandy Lake Airport (Fig. 3A). As shown in this section, the till has been eroded from the bedrock ridges and is only preserved in a small bedrock depression (Fig. 3A). During the erosion, a thin layer of pebbly sand developed on the erosion surface, interpreted initially as a beach or nearshore deposit in shallow water conditions, consistent with the field observations that indicate extensive erosion of the till by lake waves prior to the accumulation of the varved clay across the basin (Figs. A7 and A8) (Gao et al. 2017, 2018). However, as discussed below, this sand can probably also be explained as being of subaqueous outwash origin. This basal sand was sampled for OSL dating (Fig. 3A).

Section 17-CG-051 is exposed along the lakeshore on a peninsula in Sandy Lake and consists of 7 m of varved clay directly resting on bedrock (Fig. 3B). The varved clay contains about 550 couplets with thickness ranging from 20 cm thick in the lower part of the section to 1 cm or so in the upper (Fig. 3B). The number of couplets recorded here is consistent with that counted by Prest (1963) to the south on the Rathouse Bay, an elongated southwestern arm of Sandy Lake (Fig. 2). A notable feature seen in this section and elsewhere in this area is that, in the couplets, the winter layers are consistent in texture and thickness at 0.5–1 cm, whereas the summer layers become sandier and much thicker at the base of the deposit (Fig. 3B). These thick or giant summer layers also show well-developed horizontal lamination (Fig. 3B). At other localities, they even contain ripple crossbedding (Fig. A9) (Gao et al. 2017, 2018). A pair of samples about 0.5 m apart was collected for OSL dating from one of the large sandy summer layers near the base of the varved clay in this section (Fig. 3B).

Surface peat is ubiquitous in the study area. It was sampled for radiocarbon dating from two lakeshore exposures at 18-CG-41 and 18-CG-113 (see Figs. 2 and 4 for locations). The former consists, on glaciolacustrine silt and clay, of 2.1 m of peat in an abandoned strait that, as discussed later, used to connect Sandy and Niska lakes (Figs. 4, 5A, and 5B), whereas the latter contains 1 m of peat extensively exposed on the flat, former floor of the glacial lake on varved clay (Fig. 5C). As well, to the south, near the southern margin of the Sandy Lake basin, 2.3 m of peat on lacustrine mud was sampled using a hand-held Dutch soil auger in a low-lying peatland between bedrock ridges (18-CG-186) (Fig. 2).

Wave-cut shore scarps or bluffs were examined and sampled on the western end of the Opasquia moraine where recent wildfires had created a large clearance for helicopter landing and mapping (Figs. 2 and 6). The major shore bluffs there have elevations on their shoulders at 360, 345, 325, 315, 295, 285, and 276 m asl (Table 1; Fig. 6). The flat to gentle-rolling, boulder-strewn lake floors or plains on the shore bluffs are veneered with 1–2 m of littoral fine to coarse sand of nearshore origin (Figs. A5 and A6). Minor beach ridges and small wave-cut scarps exist on some of the plains (Figs. 6 and A10). Some of the minor shoreline features on the lake plains are eroded and crosscut by large shore bluffs, e.g., those at 345 and 315 m asl, owing to intense lake wave erosion and shore retreat (Figs. 6 and A10). The littoral fine to coarse sand with a massive to crudely bedded structure was sampled for OSL dating from the former lake plains near the shore bluff shoulders and bases (Fig. 6; Table 2). To the east, the moraine rises in elevation (Fig. 2). Where it reaches 400 m asl and higher at its easternmost end about 83 km east of the locality of the current study, Harvey (1979) delineated a major shoreline feature at 400 m. Somewhere in the eastern part of the moraine, a shoreline feature at 380 m asl was also reported (Harvey 1979). However, its exact location remains unknown.

In the low-lying area around Sandy Lake, a well-developed strandline of shore bluffs at 280 m asl, with their base at 278 m asl, occurs around Niska Lake (Fig. 4). Between this lake and Sandy Lake, this strandline defines a northeast-aligned, shallow corridor interpreted as a former strait that used to connect these lakes (Fig. 4). The strait is about 2.5 m deep and filled with 2.1 m of peat as seen in its northern exit exposed in the shore of Sandy Lake (section 18-CG-41) (Figs. 5A and 5B). The infill peat was, as discussed earlier, sampled for radiocarbon dating to determine when the strait and associated shoreline were abandoned.

The 3 OSL dates from the sections in the basin center on Sandy Lake are in good agreement, which are 11.2 ± 0.6 ka from the basal sand below the varved clay in the section on the bedrock ridge at the Sandy Lake Airport (17-CG-014) and 11.2 ± 0.6 and 11.8 ± 0.6 ka from the paired samples from the basal sandy summer layer of the varved clay exposed in the shore of Sandy Lake (17-CG-051) (Fig. 3; Table 2). These dates together have a mean at 11.4 ± 0.9 ka (mean ± (standard deviation + 0.6)).

The OSL dates from the shore bluffs on the Opasquia moraine are listed in Table 2. The shore bluffs on the moraine resulted from lake wave erosion associated with stepwise drop in lake level and their ages are expected to be younger for those at lower elevations. The OSL dates on the shore bluffs, however, are inconsistent with this interpretation, e.g., the shore bluff at 315 m asl having the oldest date at 15.9 ka, whereas those above it having younger dates at 9.8–12.1 ka for the shore bluffs at 325 and 360 m asl (Fig. 6; Table 2). The inconsistency can probably be ascribed to the presence of poorly bleached quartz sand grains in the sediments.

Radiocarbon dating was conducted on 20 wood and peat samples at various depths in the three peat sections described earlier, with consistent results in chronology (Table 3; see Supplementary material S2 for a full list of the dates). Radiocarbon dates on the basal wood samples indicate the onset of peat accumulation at 8.9 ± 0.1 cal ka (UOC-7893) in the southern part of the basin and 8.3 ± 0.1–8.1 ± 0.1 cal ka (UOC-7883, UOC-7894) in areas adjacent to Sandy Lake in the basin center (Figs. 2 and 5; Table 3).

The consistent OSL dates from the similar stratigraphical horizons below and at the base of the varved clay indicate the withdrawal of the ice and subsequent lake development at the center of the Sandy Lake basin at about 11.4 ± 0.9 ka (Fig. 3; Table 2). The dates are consistent with the reconstructed ice margins that indicate ice withdrawal from the western part of the Sandy Lake basin at 11.5 ka (Fig. 1) (Dyke 2004; Dalton et al. 2020). Around this time, as shown in Fig. 1, the ice margin in northwestern Ontario was located along the northwest-aligned Hartman moraine as determined by multiple radiocarbon dates (e.g., Lowell et al. 2009; Breckenridge et al. 2021); further to the north, it bent northwards and then northeastwards into the Sandy Lake basin. Conversely, recent ¹⁰Be dating work in northwestern Ontario by Lowell et al. (2021) indicates that the Hartman moraine and hence this ice margin could be more than a millennium older.

It is acknowledged that the summer layers of the varve clay may contain partially bleached quartz grains because of turbid water common in glacial lakes and the paired OSL dates from the summer layer could be too old. The same can probably also be said to the OSL date from the basal sand interpreted initially as a beach deposit (Gao et al. 2017, 2018). The projected water plane for this time indicated a deep lake with a depth estimated at over 300 m (above Sandy Lake) for the study area (Breckenridge 2015; Breckenridge 2023, personal communication). This means that the said beach sand is not possible. An alternative interpretation is that the sand is of ice-marginal, subaqueous outwash, or density flow origin, which may be difficult to distinguish from true beach deposits on their appearance alone (Smith and Ashley 1985; Teller 2005). Because sand grains from proglacial, subaqueous environments have a higher potential for partial bleaching, OSL dates from such deposits can be overestimated (e.g., Alexanderson and Murray 2012; Rittenour et al. 2015). As such, although the OSL dates are in agreement with the ice margin configuration interpreted by Dyke (2004) and Dalton et al. (2020), they are treated as the maximum age for the deglaciation of the Sandy Lake basin in the discussion below.

With the recession of the ice, varved clay up to 12.5 m thick accumulated in deep water in the Sandy Lake basin (Fig. 3B). The retreating ice then slowed and reached a standstill along the northern side of the basin, releasing a large amount of sand and gravel through meltwater to construct the Opasquia moraine (Gao et al. 2017, 2018). The moraine was smoothed and bevelled along its upper surface by lake waves, suggesting its submergence in Lake Agassiz that must then have had a lake level above this moraine, with a depth exceeding 100 m. Later, the lake lowered in lake level due to further recession of the ice and opening of lower elevation outlets, forming many strandlines typified by wave-cut large shore bluffs on the moraine (Fig. 6; Table 1).

The terrace-like, wave-cut large shore bluffs on the Opasquia moraine indicate intervals of relatively stable lake levels during the lowering of the lake, which can be correlated with those in the main basin of Lake Agassiz to the west based on projected water planes (Table 1; Fig. A3) (Thorleifson 1996; McMartin 2000). The past lake levels were estimated from the base of the shore bluffs and beach ridges (Table 1). The major strandlines as indicated by the shore bluffs at 360, 345, 315, 295, and 276 m asl on the moraine were correlated with the The Pas, Gimli, Grand Rapids, Drunken Point, and Ponton beaches in the main basin, respectively (Table 1). Those at 325 and 285 m on the moraine were unnamed, which probably developed during some minor intervals with stable lake levels (Table 1). On isobases, the shoreline at 400 m asl reported by Harvey (1979) at the easternmost part of the moraine is located 50 km northeast of the locality of the current study (Fig. 2) and is also likely correlative to the The Pas (Fig. A3).

The OSL dates on the shore bluffs on the Opasquia moraine are inconsistent in chronological succession because of the likely presence of incompletely bleached quartz sand. This can probably be ascribed to the dynamic sedimentary process where the sand grains released from the shore bluffs through lake wave erosion may have been quickly buried without full exposure to sunlight for bleaching. Among the dating results, the date at 9.8 ± 0.6 ka from the shore bluff at 345 m asl is the only one that matches the correlative Gimli beach in the main basin (10 ka) (Table 1). However, the overall chronological inconsistency makes this date dubious. As such, the OSL dates from the moraine were not discussed further. Instead, the recently updated ages of the strandlines in the main basin were adopted in the discussion below (Table 1). Thus, according to Teller and Owen (2021), the The Pas, Gimli, Grand Rapids, and Ponton strandlines occurred at 10.1, 10, 8.6, and 8.5 ka, respectively (Table 1). It is noteworthy that, even in the main basin, direct dating of the beaches has been proven challenging. There, the Gimli beach is among the few, directly dated post-Campbell beaches and its age at 10 ka is a mean of five OSL dates (10.5–9.7 ka) after removal of several outliers (13.9–11.3 ka) (Teller et al. 2018; Teller and Owen 2021). In any event, additional dating is key for better age constraints of these shorelines.

The Opasquia moraine formed after the initial development of the glacial lake in the basin center somewhere after 11.4 ± 0.9 ka but before the development on it of the first major strandline, the The Pas, which Teller and Owen (2021) assigned at 10.1 ka. This chorological framework for the moraine appears in agreement with the reconstructed ice margin configuration where the ice retreated to the northern rim of the basin at 10.9 ka and probably, after a minor re-advance, halted there until its further recession to the north at 10.3 ka (Dyke 2004; Dalton et al. 2020) or slightly later at 10.2 ka (Fisher and Breckenridge 2022, their fig. 2D). The neighbouring low-relief Hudwin moraine in Manitoba is often considered being the western extension of the Opasquia moraine (Klassen 1986; Clarke 1989; Gao et al. 2017; Gauthier et al. 2022). However, an earlier age has been also suggested for the Hudwin moraine (e.g., Dredge and Cowan 1989; Dalton et al. 2020). Additionally, some have suggested that these moraines developed together with the Agutua and Nakina end moraines located to the southeast in a single end moraine system (e.g., Thorleifson 1996; Gauthier et al. 2022), and accordingly, the Opasquia moraine could have formed as young as 9.1 ka, an age inferred for the Nakina moraine based on varve chronologies obtained from the Lake Superior basin (Breckenridge 2007; Breckenridge et al. 2012; Gauthier et al. 2022). While this alternative age assignment is equally possible given the lack of independent age constraints of these moraines, the notion for a single moraine system appears inconsistent with several other observations that the Nakina developed after the Opasquia in a different ice front setting (e.g., Zoltai 1965; Prest 1970; Dredge and Cowan 1989; Barnett 1992). According to the reconstructed ice margins (Dalton et al. 2020), the Agutua–Nakina moraines developed after 10.3 ka but before 9.6 ka, and during its further northward retreat, the ice halted and constructed the Big Beaverhouse moraine at 9 ka (see Figs. 1 and 2 for location). However, it is cautioned that the scant, direct radiometric dates available on these moraines mean that the reconstructed ice margin configuration requires confirmation through further dating (e.g., Lowell et al. 2021).

The lowering of Lake Agassiz in lake level allowed vegetation to colonize and form peat deposits on the exposed lake floor. On the bedrock uplands in the south, the radiocarbon date from the base of the peatland about 45 km south of the Sandy Lake indicates the withdrawal of the lake from this locality by 8.9 ± 0.1 cal ka (Fig. 2; Table 3). Radiocarbon dating on the peat exposed along the shore of Sandy Lake indicates peat accumulation on the emerged lake floor around 8.3 ± 0.1 cal ka (Figs. 2 and 5B; Table 3), suggesting that, by this time, Lake Agassiz had withdrawn from the Sandy Lake basin.

The last presence of Lake Agassiz in the Sandy Lake basin was indicated by the shore bluffs correlative to the Ponton shoreline around Niska Lake (Fig. 4). This strandline defines the former strait that used to link Niska Lake to Sandy Lake (Fig. 4). At that time, both lakes, or more appropriately, their prototypes, shared a same water plane at 278 m as measured from the base of the shore bluffs, and connected with each other through this strait (Figs. 4, 5A, and 5B). The strait and strandline were later abandoned, causing the separation of the two lakes after the withdrawal of Lake Agassiz and subsequent lowering of lake level. To the northwest, the Ponton shoreline is also present as a shallow shore bluff with a base at 274 m asl on the northern foothill of the Opasquia moraine (Fig. 6A; Table 1). Geographically, the strandline at this locality is on the southern margin of the Island Lake basin and its occurrence suggests the southern extent of Lake Agassiz in this basin by that time (Fig. 2). However, no samples were collected at this site for radiometric dating.

At Niska Lake, the radiocarbon date on the basal wood sample in the infill peat in the former strait indicates the disuse of this strait and abandonment of the shore bluffs of Ponton strandline by 8.1 ± 0.1 cal ka (8.2–8.0 cal ka) (Table 3). However, it is unclear whether the infill peat accumulated immediately after the abandonment of the strait. In the nearby lakeshore section 18-CG-113 to the east (Fig. 4), the basal wood remains of the peat on the flat, former lake floor were dated at 8.3 ± 0.1 cal ka (8.4–8.2 cal ka) (Figs. 4, 5C, and 5D; Table 3). Located about 6 km farther to the northeast from Niska Lake on isobases (Figs. 2 and 4), this site although similar in elevation today was likely below the water plane defined by the strait and Ponton strandline. On the reconstructed paleo-topography derived from ICE-6G for 8.5 cal ka (Godbout et al. 2023), it stood at 144 m asl, about 10 m below the strait (155 m asl). A lower paleo-elevation means that the drainage of the lake from this site should have postdated the abandonment of the strait and Ponton shoreline. Therefore, the older age (8.3 ± 0.1 cal ka) at this site suggested a hiatus between the drainage of the lake and accumulation of peat in the strait and that the radiocarbon date at 8.1 ± 0.1 cal ka from the basal wood of the infill peat in the strait is an underestimate. The true age for the abandonment of the strait and Ponton shoreline should have been no later than 8.3 ± 0.1 cal ka. The reason for the hiatus in the strait after its abandonment may have been related to continued rain-wash erosion because of the relatively steep slope gradient, especially, at its exits where the sample was collected, causing delayed peat accumulation there. By this time, the boreal forest had already been established in this area as indicated by the wood of spruce or tamarack and pollen obtained from these sections (Figs. 5D and A2). Normally, as seen in the boreal forest in this region today, it takes no more than a few decades for vegetation to colonize the clearance by wildfires, with woodlands fully established in 50 years (Bergeron et al. 1998; Anyomi et al. 2022). Applying this duration, the age would be slightly older (8.35 ± 0.15 cal ka).

The abandonment of the Ponton shoreline with a minimum-limiting age at 8.3 ± 0.1 cal ka (8.4–8.2 cal ka) is comparable with the previously assigned 8.5 cal ka for this beach (e.g., Teller and Leverington 2004; Teller and Owen 2021). However, Teller and others assigned this age based on the radiocarbon dates acquired by Barber et al. (1999) on marine shells deposited at the incursion of the Tyrrell Sea in Hudson Bay Lowlands during the final drainage of Lake Agassiz. Since then, there have been advances in methods for calibration of radiocarbon dates and estimation of marine reservoir age adjustments (ΔR) for high-latitude oceans, including the former Tyrrell Sea in the Hudson Bay region (e.g., Heaton et al. 2020, 2023; Reimer et al. 2020). Recent dating on marine and terrestrial fossils suggests that the marine incursion occurred as young as at 8.16 cal kaand this age is used to constrain a set of lower shorelines named the La Sarre, thought to represent the lake level of the final drainage instead (e.g., Roy et al. 2011; Jennings et al. 2015; Godbout et al. 2020; Brouard et al. 2021).

In the Minago River area north of Lake Winnipeg, a freshwater bivalve (Lampsilis radiata) from the Ponton beach was radiocarbon dated at 9.7 ± 0.4 cal ka (8310 ± 180 ¹⁴C BP, GSC–1679) (Lowdon et al. 1977; Klassen 1986; McMartin 2000; McNeely and Brennan 2005; Gauthier et al. 2022). This date appears too old for the true age of the molluscan shell due to the potential hardwater effect that can reach 880 and 4000 ¹⁴C years in the Agassiz and Ojibway basins, respectively (Daubois et al. 2015; Gauthier 2022). When compared with the radiocarbon dating result from the current study, the date on this freshwater shell would have a hardwater effect at 825 ¹⁴C years, which is at the high end of its range estimated for Manitoba (40–880 ¹⁴C years) (Gauthier 2022).

Recent mapping by Godbout et al. (2020) in the La Sarre area north of Lake Abitibi in Quebec indicates that the Kinojévis-type shorelines are likely 10 m higher than previously thought (Vincent and Hardy 1979; Veillette 1994; Thorleifson 1996; Roy et al. 2015). This elevation difference is consistent with McMartin’s (2000) work that also indicates a higher Ponton shoreline in similar magnitude (see Table 1 and Fig. 3A). The Ponton strandline in the Sandy Lake basin has a paleo-elevation at 155 m asl as measured on the reconstructed paleo-topography derived from ICE-6G deformation model for 8.5 ka (145 m asl on ICE-7G model) (Godbout et al. 2023). This paleo-elevation is within the cluster of the Kinojevis-type shoreline features ranging from 140 to 165 m asl (Godbout et al. 2020). It is also comparable with those from the nearest known Ponton beaches located about 215 km to the west on the eastern side of Lake Winnipeg in Manitoba where they stood at 145–157 m asl in paleo-elevation (Fig. A4; see Fig. 1 for location). The agreement in paleo-elevation among these sites suggests the prevalence of the Ponton–Kinojévis lake level across the Agassiz and Ojibway basins as previously proposed (Fig. 1) (e.g., Veillette 1994; Thorleifson 1996; Leverington et al. 2002; Clarke et al. 2004; Teller and Leverington 2004). However, as mentioned earlier, the inconsistent paleo-elevations on a more extensive Ponton strandline along the western shore of Lake Winnipeg (Fig. A4) call for further refinement and adjustment of the paleo-topography reconstructions.

Not occurring in the Sandy Lake basin but mapped in northern Manitoba and in northern Quebec and northeastern Ontario are lower shorelines named the Fiddler and La Sarre, respectively (e.g., Bell 1978; Klassen 1986; Thorleifson 1996; McMartin 2000; Roy et al. 2015; Godbout et al. 2020). Their presence indicates that Lake Agassiz drained into Hudson Bay in two separate stages, with the Ponton–Kinojévis shorelines being associated with an early drawdown of the lake and the Fiddle and La Sarre representing the last drainage (Roy et al. 2015; Daubois et al. 2015; Godbout et al. 2019, 2020). The age of the La Sarre shorelines (8.16 cal ka) was inferred primarily from the radiocarbon dates on basal plant remains of surface peat from two localities in the Ojibway basin (Kettles et al. 2000; Packalen et al. 2014; Brouard et al. 2021).

The final drainage of Lake Agassiz is thought to have freshened the North Atlantic and slowed the ocean current circulation, causing the 8.2 ka cooling event as recorded in the Greenland ice core (e.g., Barber et al. 1999; Teller et al. 2002; Kobashi et al. 2007; Kleiven et al. 2008; Roy et al. 2011; Jennings et al. 2015; Brouard et al. 2021). In the past, the lack of independent age constraints of the shorelines has led to confusions about the age and role of the two lake levels or phases in ocean freshening and resultant climate cooling (e.g., Clarke et al. 2004; Teller and Leverington 2004;  Roy et al. 2015; Godbout et al. 2020; Brouard et al. 2021; Teller and Owen 2021). Radiocarbon dating results from the current study indicate the abandonment of the Ponton–Kinojévis shorelines by 8.3 ± 0.1 cal ka, about 1.5 centuries older than the Fiddler–La Sarre shorelines inferred at 8.16 cal ka. Although closely spaced between the two sets of strandlines in chronology, many have, as discussed above, suggested that the Fiddler–La Sarre shorelines represented the lake level of the final lake drainage because their age, although inferred, is more consistent with the timing of the cooling. However, as a cautionary note, the chronologic consistency with the 8.2 ka cooling does not necessarily mean that the final drainage of the lake triggered the climate anomaly, especially, when a smaller lake size of the Fiddler–La Sarre lake level is considered for the volume of water available to freshen the ocean. Alternatively, some have invoked the collapse of the ice saddle over Hudson Bay as the primary cause because of the recent model simulations that often failed producing this cooling that lasted for 1.5 centuries (Gregoire et al. 2012; Matero et al. 2017; Lochte et al. 2019; Gauthier et al. 2020). To test these competing hypotheses, it is important to examine the stratigraphical records, especially, in the former Tyrrell Sea basin in Hudson Bay Lowlands where abundant sedimentary records exist, but only limited studies have been conducted (e.g., Skinner 1973; Roy et al. 2011).

The deposition of extensive varved clay indicates the withdrawal of the ice sheet from the center of the Sandy Lake basin with an OSL-dated maximum age at 11.4 ± 0.9 ka. After its further retreat, the ice halted on the northern side of the basin and generated abundant glaciofluvial sand and gravel to construct the Opasquia moraine sometime before the development on the moraine of the first major strandline of Lake Agassiz, the The Pas inferred at 10.1 ka from previous work. Subsequent stepwise lowering in lake level resulted in many prominent strandlines typified by wave-cut large shore bluffs on this moraine. The major strandlines at 360, 345, 315, 295, and 276 m asl on the moraine were correlated, respectively, with the The Pas, Gimli, Grand Rapids, Drunken Point, and Ponton beaches in the main basin of Lake Agassiz. The last presence of Lake Agassiz in the Sandy Lake basin was indicated by the shore bluffs of the Ponton shoreline around Niska Lake. At that time, Niska and Sandy lakes were connected via a shallow strait, with a lake level at 278 m asl. The strait and associated Ponton strandline were abandoned after the withdrawal of Lake Agassiz from this basin. Radiocarbon date at 8.1 ± 0.1 cal ka on a basal wood sample in the infill peat in the strait is probably an underestimate because of the continued rain-wash erosion in the strait. Instead, the radiocarbon date at 8.3 ± 0.1 cal ka on the basal wood remains of the surface peat at a nearby site on the flat, former lake floor was taken as a better approximation for the age of the shoreline abandonment. This date, although a minimum-limiting age, provides the hitherto best possible age constraint of the Ponton–Kinojévis shorelines, which many hypothesize represent one of the major lake levels during the final drainage of Lake Agassiz into Hudson Bay but have never been adequately dated before.

This work was supported by the Ontario Geological Survey. The author wishes to thank the First Nation communities for allowing mapping on their traditional lands and Blayne Anishinabie, Jason Dyer, Kyle Dzuirban, Terrance Fiddler, Douglas Kakegamic, Nicholas Szumylo, and Mia Tullio for field assistance. Kei Yeung aided in manipulation of remote sensing imagery and GIS applications on many occasions. Charles Turton did pollen analysis. Sebastien Huot at Illinois State Geological Survey conducted OSL dating, and W.E. Kieser at A.E. Lalonde AMS Laboratory, University of Ottawa, undertook radiocarbon dating. The manuscript was improved by discussions with James Teller, David Leverington, Michelle Gauthier, Pierre-Marc Godbout, and Etienne Brouard, and the comments provided by Andy Breckenridge, an anomalous reviewer, and Roger Paulen, the Associate Editor. The manuscript was submitted for publication with the permission of the Director of the Ontario Geological Survey.

Data generated or analyzed during this study are provided in full within the published article and its supplementary material. The surficial geology maps shown in Fig. 4 are from the Ontario Geological Survey’s Preliminary Maps P.3678 and P.3691 at 1:100 000 scale, which, once released (soon), will be available for open access via GeologyOntario at https://www.hub.geologyontario.mines.gov.on.ca/.

Conceptualization: CG

Data curation: CG

Formal analysis: CG

Funding acquisition: CG

Investigation: CG

Methodology: CG

Project administration: CG

Resources: CG

Software: CG

Supervision: CG

Validation: CG

Visualization: CG

Writing – original draft: CG

Writing – review & editing: CG

This research was supported by the Ontario Geological Survey via the Far North terrain mapping project (project No. FN-17-002).

Supplementary data are available with the article at https://doi.org/10.1139/cjes-2023-0014.