Preservation of probable MIS 7 deglacial and nonglacial deposits near the edge of the Hudson Bay Lowland in Manitoba, Canada

It has long been thought that glacial sediments at the surface belong to the youngest glacial episodes, and that preservation of older sediments deposited during previous glacial and interglacial episodes is rare in northern America (Hildes et al. 2004; Melanson et al. 2013; Batchelor et al. 2019). And yet, Marine Isotope Stage (MIS) 5- and 7-aged nonglacial sediments were documented at the base of lakes in the high arctic (Briner et al. 2007, MIS 5, 130–71 ka; MIS 7, 243–191 ka; Lisiecki and Raymo 2005; Miller et al. 2022), and in a few stratigraphic sections from eastern James Bay (Dube-Loubert et al. 2013). Are these examples extremely rare, or has a lack of absolute dating influenced the paradigm? It is now accepted that fragments of older landscapes are found preserved (inherited) beneath both the former Laurentide (Ross et al. 2009; Trommelen et al. 2012, 2013, 2014; Trommelen and Ross 2014; McMartin et al. 2021), Scandinavian (Kleman 1992; Kleman and Stroeven 1997; Kleman and Hättestrand 1999; Kleman et al. 1999; Tylmann et al. 2012; Ebert et al. 2015; Helmens et al. 2018; Salonen et al. 2018), and British–Irish ice sheets (Finlayson et al. 2010; Rex et al. 2023). Henceforth, there should also be fragments of older landscapes buried in the Quaternary stratigraphic record.

Lithostratigraphy and allostratigraphy can be used to identify disjoint patches of older sediments in 3D (Gauthier et al. 2019; Hodder et al. 2024). Exactly how old those sediments are, however, is an ongoing puzzle. Radiocarbon dating requires the presence of organic material and has an upper age limit of ∼50 ¹⁴C ka. Paleomagnetic geochronology is restricted sediments older than 778 ka, since polarity has been normal between the present and MIS 17 (Barendregt 2011; Andriashek and Barendregt 2017; c.f. Fullerton et al. 2004). In situ cosmogenic nuclide methods cover a larger time frame, but may be complicated by inheritance and overprinting (Bierman et al. 1999; Colgan et al. 2002; McMartin et al. 2022). Recently though, optical dating of quartz grains has provided ages consistent with MIS 5 and 7 on sorted sand deposits in the western Hudson Bay Lowland (HBL, Hodder et al. 2023; Hodder et al.,2024). Based on extensive fieldwork and regional till characterization (Gauthier et al. 2022b), we identified two gravel pits near Gillam, Manitoba, where the stratigraphy suggests near-surface (<4.5 m depth) subtill sediments were deposited during an old (pre-MIS 5) interglacial period. Pollen and macrofossil data are used to reconstruct the paleo-environment during deposition of a nonglacial stratigraphic bed, and optical dating of quartz grains is used to constrain the age of the underlying stratigraphic bed.

The HBL is a critical location for studying the ice-sheet history of North America because of its location near the geographic center of the ice-sheet complex, as well as the extensive (yet fragmentary) preservation of glacial and nonglacial stratigraphy. The western HBL was glaciated by ice flowing from two major dispersal centres (Keewatin and Québec–Labrador), leading to a dynamic and complex paleoglaciology over at least three glacial cycles (Dredge and Nixon 1992; Dredge and McMartin 2011; Gauthier et al. 2019). The river banks in the western HBL expose up to 60 m of stratigraphy, containing sparse but widespread subtill organic-bearing sediments (Tyrrell 1916; Craig et al. 1967; Netterville 1974; Nielsen and Dredge 1982; Klassen 1986; Nielsen et al. 1986; Dredge et al. 1990; Roy 1998; Hodder et al. 2017). This initial work led to the identification and naming of four till units, a nonglacial unit, and a second older interglacial period as evidenced by a paleosol and possibly marine sediments (Nielsen et al. 1986; Dredge et al. 1990). Previous work typically interpreted all organic-bearing sorted-sediment units as interglacial and correlated them as one unit referred to as the “Nelson River Sediments”, often at regional- to reconnaissance-scales of study (Nielsen et al. 1986; Dredge et al. 1990; Roy 1998). To the east in Ontario, all nonglacial sediments were similarly lumped into the Missinaibi Formation (Terasmae and Hughes 1960; Skinner 1973). Various attempts have been made to date these subtill nonglacial deposits over time, through amino acid racemization in marine shells (Andrews et al. 1983), radiocarbon dating (McDonald 1968; Dredge et al. 1990), and thermoluminescence and optical dating (Berger and Nielsen 1990; Roy 1998). These methods have met with mixed success; namely, identifying ice-free conditions during some part of MIS 5 and 7 (Allard et al. 2012; Gao et al. 2020). Some authors proposed additional ice-free conditions during MIS 3 (Dalton et al. 2019, 57–29 ka, Lisiecki and Raymo 2005; Dalton et al. 2016), but the feasibility of ice-free conditions during that time, along with the validity of the geochronological dataset, is debated (Miller and Andrews 2019; Kerr et al. 2021; Gauthier et al.,2024). Thorough reanalysis of western HBL stratigraphy, focusing on till, has now led to the recognition of multiple till units that are spatially discontinuous yet indicative of a long record of changing ice-flow orientations (Hodder et al. 2017, 2020, 2024). Using quartz optical dating methods together with till stratigraphy has led Hodder et al. (2023) to determine that 6 out of 7 of those potential western HBL MIS 3-aged sites were more likely deposited during or prior to MIS 5 (Fig. 1). The question then remains—was the first occurrence of subtill nonglacial sediments in the stratigraphic record deposited during MIS 5 or has correlation been oversimplified?

Two gravel pits were identified for further study (Pit 1 and 2, Fig. 2), as part of a larger western HBL stratigraphy and mapping program (e.g., Gauthier et al. 2019, 2022b). These pits, situated along a road that leads to Gillam Manitoba, provide excellent 3D access to the local stratigraphy. In contrast, regional mapping sites (Fig. 2) were examined through a mixture of Dutch auger holes (1.2 m deep), hand-dug holes, shoreline sections, road cuts, and shallow construction borrow pits.

Glacial Lake Agassiz inundated most of the region during the most recent deglaciation (Klassen 1983). There are iceberg scours east of Pit 1, between 175 and 180 m above sea level (asl), meaning Lake Agassiz was likely shallow and possibly short-lived (Gauthier et al. 2022b). Very little fieldwork has been completed north of the study area, and 1:250 000 regional-scale mapping was completed without ground-truthing (Klassen and Netterville 1980). Though mapping is incomplete, glaciolacustrine sediments were not mapped within 8 km of Pit 2 ¹. The postglacial Tyrrell Sea transgressed to a maximum of ∼135 m asl (Gauthier et al. 2020b), approximately 29 km east of Pit 1¹.

The study area is less than 50 km west of the boundary between the Precambrian shield and Paleozoic Hudson Bay Basin (HBB) (Fig. 1). The Precambrian shield is composed mainly of granite and gneiss (Corkery et al. 1992; Bohm et al. 2006). The Paleoproterozoic Fox River belt also crosses the area, consisting of metasedimentary and metavolcanic rocks (Rinne 2020). The flat-lying HBB onlaps both and is composed of calcareous sandstone, limestone, and dolomite (Nicolas and Armstrong 2017). Depth to bedrock is largely unknown in this understudied area (Hodder and Gauthier 2021). In the study area, bedrock crops out at ∼140 m asl along the shores of Stephens Lake and the Limestone River (Fig. 2). Interestingly, pits 1 and 2 are dug into a thick-sediment interfluve that is situated between these two regions of outcropping bedrock. Quaternary sediments thicken towards Pit 2 in the east.

Clear descriptions of sediment exposures are necessary to help establish the Quaternary stratigraphy in an area that is largely devoid of radiometric ages and subtill organic units. The sedimentology, composition, organic matter, structure, and stratigraphic contacts of each exposure was described in detail. To aid field descriptions, eight diamict samples were taken at the gravel pits to document moist matrix colour (Munsell Color-X-Rite Inc. 2015), matrix texture (<2 mm size fraction), total carbonate content (<63 µm), and clast lithology (2–8 mm size-fraction). An additional gravel sample from Pit 1 was collected to compare the clast lithology with diamict samples. Units were later defined using a lithostratigraphic approach.

For each sample, the lithological composition of ∼400 clasts (2–8 mm) per sample were identified under an optical microscope. Clast lithologies were initially separated into 18 classes and later grouped into four categories, HBB, undifferentiated granitoids–gneissic (GR), identified greywacke and greenstones (IGG), and undifferentiated greywacke and greenstones (UGG). HBB clasts, sourced from the Paleozoic platform (Nicolas and Armstrong 2017), include calcareous clasts (grey, tan, white, or pink) and rarer (<1 ct.%) Kenogami Formation calcareous sandstones and hard red exotic sediment. IGG are metavolcanic and metasedimentary rocks that are likely sourced from the underlying greenstone belt. UGG include all “unidentifiable” small grey/black/green clasts that could be sourced from either the underlying greenstone belt, the Omarolluk Formation of the Belcher Islands in eastern Hudson Bay (Ricketts and Donaldson 1981; Prest et al. 2000), or the Sutton Inlier (Stott et al. 2010).

Clast fabric measurements were conducted within the diamict in order to determine the strain direction, as a proxy for ice-flow orientation (Holmes 1941; Andrews and Smith 1970). Sites were chosen based on uniformity of diamict, where no sand lenses or discontinuous bedding was present. At each site, at least 30 clasts were carefully excavated and measured from within a “box” consisting of three vertical faces of different orientations, over a maximum distance of 30 × 30 × 30 cm. Accepted clasts included in these analyses were (1) free to rotate in the matrix at the time of deposition (not clast-supported or close to much larger clasts), (2) rod, tabular-rectangle or wedge-shaped (ratio of the A:B axis was 1.5 or greater), (3) had a plunge of the a-axis less than 60° (averaged 18° plunge with a standard deviation of 11°), and (4) a plunge of the b-axis less than 60° (c.f. Holmes 1941). At some sample sites, parallel striations on the upper surfaces of in situ, lodged cobble- to boulder-sized clasts were measured instead of till fabrics.

Interpreting the ice-flow orientation from a clast fabric within a diamict interpreted as till can be challenging. Clast fabric data were represented graphically using Rockware Stereostat V1.6.1. Ice-flow interpretations were assigned to the data based on stereonet patterns, rose-diagram patterns, and regional ice-flow history (Gauthier et al. 2019), in addition to the eigenvector of the largest eigenvalue or V₁ (Mark 1974). Multimodal till fabrics were not used to interpret ice-flow direction (c.f. Hicock et al. 1996)

Pollen assemblages were analyzed from 14 organic-rich sediment samples at the University of Toronto following standard palynological techniques, namely, KOH to remove humic acid, HF to digest the sediment matrix, and acetolysis to clean pollen grains (Faegri and Iversen 1975). Sodium polydentate, a heavy liquid, was used to remove any sediments that remained after the aforementioned processing (Zabenskie 2006; Campbell et al. 2016). Ceramic palynospheres were used to estimate pollen concentration at each interval (Kitaba and Nakagawa 2017). At least 150 grains of herb, arboreal, and shrub species were counted in each interval, using the identification key of McAndrews et al. (1973) and a compound microscope at 400×. Intervals with low pollen concentration (less than 5000 pollen grains/cm³⁾ or those containing many broken pollen grains (>10% of counted grains) were excluded from further consideration. Stratigraphic plots of pollen data were created using R package “rioja” (Juggins 2015).

We used the pollen-based modern analogue technique (Overpeck et al. 1985) to reconstruct total annual precipitation and mean summer temperature (June, July, and August) during the fossil interval. This involved the statistical comparison of the herb, arboreal, and shrub fossil pollen data to an extensive calibration set (4882 sites spanning North America; Whitmore et al. 2005; Dalton et al. 2017). Only those sites from the calibration set with a squared chord distance dissimilarity of <0.25 were considered appropriate analogues. The three closest analogues to each fossil interval were used to reconstruct paleoenvironmental variables using cross-validation (n = 500 bootstrap iterations) and the R package “analogue” (Simpson 2007; Simpson and Oksanen 2014).

Plant macrofossils were also examined at the University of Helsinki for seven intervals that had the highest amount of organic content. For each interval, a sample measuring 20 cm³ was washed using water over a 100 µm mesh sieve. Plant macrofossils were examined using a stereomicroscope. Finally, the organic content for each examined interval was estimated using loss-on-ignition methods (Heiri et al. 2001).

Chronological constraints were obtained using both radiocarbon and optical dating. A flattened piece of wood from section Pit 1-southeast was submitted for accelerator mass spectrometry ¹⁴C dating. Standard acid–alkali–acid pre-treatment was applied to the wood (Crann et al. 2017). For optical dating, quartz grains extracted from waterlaid sands at section Pit 1-southeast and section Pit 2-west were prepared at University of the Fraser Valley ². Quartz extracts of sediment were mounted on aluminum discs, each aliquot containing 50–100 grains. Ideally, aliquots containing a single grain should be measured so that various age populations, if they exist, can be resolved. But this is often not practical for samples where few grains are suitable for dating, and for “old” samples where long laboratory irradiation times are required; in each case, an inordinate amount of machine time would be required. For these samples, measuring aliquots with 50–100 grains assured that a sufficient number of aliquots could be accepted for dating. Equivalent dose (D_e) values were measured using a single-aliquot regenerative-dose (SAR) protocol (Murray and Wintle 2000, 2003) following the same laboratory protocols described in Hodder et al. (2023),². Aliquots were accepted for SAR analysis if their natural signals showed a clear “fast” component (e.g., Fig. 3) and if all the SAR quality control criteria were met (c.f. Hodder et al. 2023). Finally, the efficacy of the chosen SAR protocol was tested by means of dose recovery tests (Wintle and Murray 2006) as described in Hodder et al. (2023). Optical dating of feldspar (i.e., infrared-stimulated luminescence, IRSL, dating), typically K-feldspar, is often used for samples for which the quartz signal is not suitable for dating, or in cases where it is saturated. This is because the IRSL signal from K-feldspar saturates at higher doses. However, K-feldspars suffer from anomalous fading (AF), which results in calculated ages being underestimated. AF can be corrected for, but only for samples that have a linear dose response (i.e., “young” samples, c. 20–50 ka, c.f. Huntley and Lamothe 2001), and this is not the case for the samples studied here. Other options could have been using post IR–IR signals from feldspar, for which AF is typically lower, and sometimes (rarely) insignificant, or the thermal-transferred signal from quartz (TT-OSL dating), which can yield dose responses that are linear to much higher doses. However, both the post IR–IR and TT-OSL signals need prolonged exposure to direct sunlight to be reset, and the depositional environments sampled for this study are not conducive to that. Environmental dose rates were determined by drying and milling a representative fraction of the bulk sediment used for dating and determining U, Th, and K concentrations by neutron activation analysis (Table 1). As our samples consisted of well-drained sand, the as-collected water content was used in each case in the dose rate calculations, with a ±10% uncertainty that is intended to account for moisture fluctuation over time. The contribution from cosmic ray radiation was estimated using present burial depths and the relationship described in Prescott and Hutton (1994).

This study has identified five informal lithostratigraphic units at Pit 1, recognized based on standard lithologic criteria, including unit thickness, texture, sedimentary structure, colour, clast lithology, fossil content, stratigraphic position, and nature of contacts. The establishment of the five lithostratigraphic units are based on observations made at three exposures: Pit 1-west, Pit 1-southeast, and Pit 1-northwest ³ (Fig. 4). Simplified stratigraphy is shown in Fig. 5, while each unit is described below.

Unit A is gravel with granule- to cobble-sized clasts that is at least 8 m thick and extends to the base of Pit 1 (Fig. 5). The upper contact is at ∼161–165 m asl throughout the pit, suggesting a tabular depositional geometry or architecture (Fig. 4c). The gravel is poorly sorted, and matrix- to clast-supported (Fig. 6). The larger clasts (framework grains) are dominantly rounded to subrounded, and the gravel contains rare, abraded mm-sized marine shell fragments. The clast lithology, from one sample, is dominantly carbonate (62 count (ct.)% HBB, 19 ct.% GR, 13 ct.% IGG, 6 ct.% UGG). There is planar bedding and weak cross-stratification, defined by beds of cobbles (Fig. 6a) and (or) changes in clast size (Fig. 6c). At the Pit 1-northwest section, cross-bedded gravels dip at 15° towards 020°³.

Unit A is overlain unconformably by unit B-1, both subhorizontally (Fig. 6d) and undulatory as infill of a depression (Fig. 6e) or by units D (Fig. 6f) or E. Fractured limestone clasts are present in the gravels near the upper contact with both units (Figs. 6b and 6f). The thickness, bedload, sorting, architecture, and matrix- to clast-supported nature of these gravels suggest deposition by relatively high-energy currents in a proximal glaciofluvial environment.

Unit B is a sandy to pebbly sand that is between 0 and 2.5 m thick (Fig. 5). At Pit 1-southeast, unit B consists of three lithofacies: a basal 0.5 m thick, well-sorted, fine-grained sand (Figs. 7a and 7b); a middle 1.0 m thick, moderately sorted, fine- to very-coarse grained sand with rare pebbly beds (Figs. 7b and 7c); and an upper 0.5 m thick, moderately sorted, gravelly sand (Fig. 7d). The basal lithofacies (B-1) is laminated and cross-bedded, the middle lithofacies (B-2) has laminations outlined by heavy minerals that dip 45° toward 055° (northeast), and the upper lithofacies (B-3) is massive to stratified. There are 1–2 mm sized marine shell fragments present in lithofacies B-2 and B-3. All three lithofacies within unit B at Pit 1-southeast have normal faults and collapse structures (Figs. 6e, 7b, and 7c). Faulting does not appear to extend downward into Unit A (Figs. 6d, 6e and 7a). The upper contact of unit B has been bulldozed along most of Pit 1 (Fig. 4), except at the Pit 1-west section where it has a sharp erosional contact with Unit D (Fig. 6d). The sand-size grains and moderate sorting suggest deposition by relatively low-energy currents in a fluvial or lacustrine environment, while their presence at the height of land suggests a glacial source.

Unit C consists mainly of sorted silt and clay that is at least 3 m thick at Pit 1-southeast (Fig. 5)⁴. The top of the exposure was significantly bulldozed during gravel extraction, and no upper contact could be confidently established. The lower contact of unit C was not exposed, since gravel pit excavation had created a horizontal bench at this location (Fig. 4c). However, unit C is directly above unit B (within 0.2 m of lithofacies B-3 using a soil probe) and laterally equivalent in elevation to unit D. Unit C consists of three conformable lithofacies. Lithofacies C-1 is a 1.0 m thick grey to black laminated silt that is horizontally bedded with oxidized sand (Fig. 8a). The fine-grained sand beds are 0.1–10 cm thick. The silt beds are 3–10 cm thick and contain organics including wood (Figs. 8b and 8c). Lithofacies C-2 is a 0.2 m thick interbedded grey silt and coarse-grained sand to pebbly gravel (Fig. 8d). The sand or gravel beds are 0.01–0.08 m thick and poorly sorted. The thickest beds coarsen upward from coarse-grained sand to pebbly sand. Lithofacies C-3 is a 2.0 m thick dense, laminated light brown silt (Figs. 8a and 8e). The uppermost 1.3 m is heterogeneous and chaotically bedded, which likely means they are anthropogenically disturbed. All three lithofacies have undeformed subhorizontal bedding planes. Unit C is interpreted to have been deposited in a nonglacial quiet-water pond or floodplain environment based on the fine-grained texture, laminated structure, and organic content of the sediment.

Pollen assemblages were analyzed from samples taken from a 1.45 m thick section of unit C (through all three lithofacies, from C-1 through to C-3, Fig. 5). Organic content ranged from 2% to 8% in the examined interval. Seven of the 14 intervals were excluded based on poor preservation (Dalton et al. 2022). Intervals with sufficiently preserved pollen (Fig. 9) were characterized by typical boreal forest groups: typically 60%–80% Picea with less than 20% Pinus and Betula, along with Salix and Alnus in lesser amounts. An exception to this assemblage was noted at 0.95 m, which contained >80% Betula grains. Herbaceous groups generally consisted of <10% Poaceae and Asteraceae. In terms of wetland indicators, very few were noted in these samples, with Cyperaceae and Sphagnum comprising generally <10% of the herb, arboreal, and shrub count. Aquatic indicator Pediastrum was present in three intervals.

Squared chord distance dissimilarity analysis suggests the modern analogue technique is an appropriate tool for reconstructing paleoenvironmental variables (Dalton et al. 2022). The examined fossil interval contained pollen assemblages similar to those in the present-day boreal forest, peatland, and tundra regions of North America. During the fossil interval, mean summer temperature was similar to present day, ranging from 9.4 to 15.2 °C (present day is 12.7 °C). Reconstructed total annual precipitation was also similar to present day (482 mm) with some intervals earlier in the fossil record suggesting slightly higher, up to 920 mm. The sample with an anomalously high Betula count was excluded from paleoclimate reconstruction.

Plant macrofossils are present in all examined intervals, but diversity is low (Fig. 10). The most observed macrofossils are conifer remains, bark, and needles, most of which were intact and sometimes large (Fig. 8c). Many intervals also contained Najas flexilis seeds. Charcoal is also present in many intervals, along with Cristatella statoblast. Rarely, Potamogeton and Carex seeds are preserved. There is no obvious succession in the fossil intervals, although charcoal is slightly more abundant at the bottom of the section.

Unit D is a brown (Munsell colour 10YR 5/3) to olive-brown (Munsell colour 2.5Y 4/3), matrix-supported, massive, clayey sandy silt diamict that is 0–2.2 m thick (Fig. 5). The diamict is highly overconsolidated, jointed, and contains 10%–15% clasts (Fig. 11). The clast lithology, from two samples, is dominantly carbonate (60 and 70 ct.% HBB, 20 and 16 ct.% GR, 1 ct.% IGG, 19 and 13 ct.% UGG). Two clast fabric measurements were conducted within unit D at Pit 1-west (Fig. 5b). The lower clast fabric is unimodal (c.f. Hicock et al. 1996) and suggests shear stress was to the west (∼268°, S₁ = 0.72, Table 2). The upper clast fabric is weaker (S₁ = 0.54, Table 2), multimodal, and cannot be used to interpret shear stress. Unit D diamict unconformably overlies units A or B and is situated at the same elevation as unit C at Pit 1-southeast (Fig. 4c). The massive, highly overconsolidated nature, texture, clast shape (bullet-shaped and striated), clast fabric, and lateral continuity indicate that unit D is a subglacial traction till (Evans et al. 2006).

Unit E consists of sorted silt and clay, at the surface, that overlies unit D or unit B (Fig. 4, Fig. 5). At Pit 1-west, 1.0 m of dark brown, horizontally bedded, silty clay (lithofacies E-1, Fig. 12a) has a transitional lower contact with diamict (unit D), which consists of 0.1–0.3 m of silt, clay, and sand that is mixed or interbedded with diamict. A 2022 pit expansion (pit 1-northwest, Fig. 4a) exposed 0.4 m of olive-brown (Munsell colour 2.5Y 4/3), matrix-supported, massive sandy silt diamict (lithofacies E-2, Fig. 5) sharply overlying gravel³. The diamict is consolidated, not jointed or stained, and contains 15% clasts. The clast lithology, from one sample, is dominantly carbonate (63 ct.% HBB, 22 ct.% GR, 10 ct.% IGG, 4 ct.% UGG). A weak bimodal fabric measured in this diamict suggests shear stress was to the west (∼278°, S₁ = 0.57, Table 2). Stratified diamict and interbedded clay, silt, and diamict (0.4 m thick, lithofacies E-3, Fig. 5)³ overlie the massive diamict, which is then overlain by horizontal rhythmically bedded silt and clay (2 m thick, lithofacies E-4)³. To the west of this site, horizontal rhythmically bedded silt and clay (1–2 m) overlies 0.5–3.5 m of folded rhythmically bedded silt and clay (Fig. 12c and 12d) that contains <1% pebbles. These folded beds appear to infill a local depression. The northern Pit-1 wall is mostly colluvium, but at some point, a diamict pinches out to the east, and up to 2.0 m of rhythmically bedded silt and clay (Fig. 12b) sharply overlie unit A gravel.

Regionally, massive to rhythmically bedded silt and clay was deposited at surface throughout most of the study area, and generally varies in thickness from 0.1 to 1.0 m (averaging ∼0.4 m, Gauthier et al. 2022b)¹. The fine texture of these sediments, together with the presence of dropstones and rhythmic bedding, suggests that these sediments were deposited in a glacial lake; either glacial Lake Agassiz or a smaller regional proglacial lake at the edge of the ice margin (Klassen 1983; Thorleifson 1996; Gauthier et al. 2020b). The folds observed in the lower bed of lithofacies E-4 are interpreted to be glaciotectonic structures. The folded beds suggest that after deposition of the lower rhythmites, but before deposition of the upper undeformed beds, a re-advance of an ice margin occurred in the area (e.g., Stone and Koteff 1979; Stephens Lake re-advance of Gauthier et al. 2022a).

An additional gravel pit was available for study, 46.5 km east of Pit 1 (Fig. 2). This pit is inactive and is situated on the edge of a meltwater corridor (Fig. 13a). Observations at Pit 2 are based on two exposures near the top of the pit, as the remainder of the pit is covered by extensive colluvium (Pit 2-west and Pit 2-east, Fig. 13)⁵,⁶. The stratigraphy exposed consists of three units (Fig. 14), though only the upper unit was described in detail. The lowermost unit (P-A) is gravel that is well to poorly sorted, massive to horizontally bedded, and contains granule to small pebble-sized clasts (Fig. 15a)⁵. The beds vary in clast concentration from clast-supported gravels to matrix-supported gravelly sand. A soft-bodied fossil within a limestone cobble was collected from the surface of the lower gravels and was identified as Devonian in age (pers. comm. Dave Rudkin 2021)⁵. This rounded cobble would have travelled at least 180 km towards the southwest from its source. Unit P-A is at least 5.2 m thick, and the lower contact is obscured. The upper contact with a sandy to pebbly sand unit (P-B) is conformable and is situated at ∼137 m asl. Unit P-B is approximately 2.5 m thick, and consists of massive well-sorted sand (0.6 m thick), massive to cross-bedded sand and gravel with 30% granule to small pebble-sized clasts (0.9 m thick), horizontally bedded well-sorted sand (0.7 m thick), and moderately sorted sand with 5% pebbles (0.2–0.4 m thick; Figs. 14 and 15)⁵. No faulting was observed, and some beds contain rare, abraded mm-sized marine shell fragments. Unit P-D sharply overlies unit P-B (Fig. 13), and consists of two lithofacies^5,6. Lithofacies P-D-1 is a brown to olive-brown, matrix-supported, massive, clayey sandy silt diamict that is 2.8–4.3 m thick (Fig. 14). Lithofacies P-D-1 contains 10%–15% clasts that are dominantly carbonate (59–68 ct.% HBB, 9–14 ct.% GR, 1–16 ct.% IGG, 12–21 ct.% UGG; n = 5). Consolidation increases with depth, jointing is present throughout, and there is minor oxidation staining along joints (Figs 15c–5e). Three spread unimodal clast fabric measurements were collected within unit P-D-1 at Pit 2-east, in addition to measurements of striated and lodged cobbles at Pit 2-west (Fig. 14). The lower clast fabric suggests shear stress was to the west (∼271°, S₁ = 0.70, Table 2). The weaker middle clast fabric suggests shear stress was to the northwest or southeast (140–320°, S₁ = 0.55, Table 2). The upper clast fabric suggests shear stress was to the south–southwest (∼208°, S₁ = 0.79, Table 2). The second, lithofacies P-D-2, is a light yellowish-brown, matrix-supported, massive to stratified, silty sand diamict that is 0–1.8 m thick (Figs. 13 and 14). Lithofacies P-D-2 contains 15%–20% clasts that are dominantly carbonate (64 ct.% HBB, 22 ct.% GR, 1 ct.% IGG, 13 ct.% UGG; n = 1). Lithofacies P-D-2 is soft, has minor jointing structure, with no oxidation on joint surfaces, and is relatively less consolidated than the underlying diamict (Figs. 13b and 15f). The clast fabric completed within unit P-D-2 is unimodal and suggests shear stress was to the south–southwest (V₁ = 225, S₁ = 0.81, Table 2). The massive and consolidated nature, texture, clast-shape (bullet-shaped and striated), clast fabrics, mixed clast lithology, and lateral continuity indicate that units P-D-1 and P-D-2 are subglacial traction tills (Evans et al. 2006).

The units at pit 2 are likely correlative to those observed in pit 1, though specific lithofacies differ. At each pit, gravel (unit A) or sandy gravel (unit P-A) is overlain by sands and pebbly sands (unit B and P-B) and a dark coloured, highly overconsolidated, clayey diamict (unit D, P-D-1; Figs. 5 and 14). Units C and E are missing from pit 2, while an additional light coloured, consolidated, sandy diamict overlies the darker diamict at Pit 2-west (unit P-D-2, Fig. 14a).

A non-finite radiocarbon age was obtained from a 4 by 14 cm flattened wood piece, collected from within lithofacies C-1 at Pit 1 (>50 ka ¹⁴C BP, UOC-13391, Fig. 8b; Table 3).

Two samples were collected for optical age determination. At Pit 1, the laminated, moderately sorted, fine- to very-coarse grained sand (lithofacies B-2) was sampled (sample 115-20-200-OSL011) at 7 m below the ground surface at Pit 1-southeast (Figs. 5 and Fig. 16a). At Pit 2, the horizontally bedded, well-sorted medium- to coarse-grained sand (unit P-B) was sampled (sample HBL22-1) at 3.5 m below the ground surface at Pit 2-west (Figs. 14 and Fig. 16c). The depositional environment for both sand samples are interpreted as distal outwash or glaciolacustrine.

The luminescence signals recorded from all aliquots were bright and showed the initial rapidly decaying signal indicative of the desired thermally stable “fast component”. For each aliquot, dose response data were fitted with a saturating exponential plus linear function (e.g., Fig. 3). D_e was determined using 24 (sample 115-20-200-OSL011) and 44 (sample HBL22-1) aliquots; 21 and 25 aliquots were accepted, respectively (Table 4). Some aliquots were rejected due to interpolation of the natural signal beyond the last measured dose point on their respective dose response curves (nine aliquots for sample HBL22-1 and five for 115-20-200-OSL01), or due to the dose response signal being saturated (two aliquots for sample HBL22-1 and three for sample 115-20-200-OSL011). A few aliquots (two or three) were rejected from each sample because recycling ratios and recuperation values were too high (cf. Murray and Wintle 2000). A dose recovery test on sample HBL22-1 showed that a given dose like the sample's D_e can be adequately recovered (Table 4).

Ages were calculated using both a central age model (CAM) and minimum age model (MAM) using accepted aliquot D_e values (Galbraith et al. 1999; Galbraith and Roberts 2012). (Fig. 16; Table 4). For sample 115-20-200-OSL011, which had a relatively low overdispersion (OD) value of 18 ± 4%, the CAM and the MAM gave ages that are consistent with each other at 1σ (222 ± 14 and 214 ± 19 ka, respectively). For sample HBL22-1, which has an OD value of 29 ± 5%, the ages are 286 ± 24 ka (CAM) and 214 ± 22 ka (MAM), which are also consistent with each other at 1σ. However, when the depositional environments (distal glacial outwash or glaciolacustrine sediments), we consider the MAM ages to be a better approximation of the true depositional age. This is because heterogeneous bleaching of the quartz grains before burial in the inferred depositional environments is expected. For this reason, even the ages derived from the MAM may overestimate the time of deposition. But since the ages derived from the two sites are consistent with each other, we assume that any overestimation is small.

The lower stratigraphy at the gravel pits record a deglaciation of the western HBL. Thick gravels deposited by relatively high energy currents suggest deposition in a proximal glaciofluvial environment (units A, P-A), either within a braided outwash plain or as a subaqueous fan(s). Sands and gravelly sands (units B, P-B), sometimes faulted, overlie the gravels, and indicate deposition within a lower energy environment. While the non-finite ¹⁴C marine shell fragments in unit B could indicate deposition within an older ocean that had a higher marine limit than the Holocene (∼30 m higher), the dip of bedding in lithofacies B-2 indicates drainage towards Hudson Bay. Thus, while the bay may have been ice-free, it is more likely that the shell fragments were eroded from pre-existing sediment. Similar shell fragments occur in the local tills (Roy 1998; Dredge and McMartin 2011), indicating that glacial transport of shells outside of the marine basin is common. The normal faulting, collapse structures, and sedimentology of the sands and gravelly sands (unit B, Pit 1 only) may be indicative of distal outwash, or perhaps melting ice within glaciolacustrine deposition (Mcdonald and Shilts 1975). The optical age estimates of the time of deposition of these sands, which is 214 ± 19 and 214 ± 22 ka (units B-1 and P-B; Table 4) using the MAM. These ages correspond to MIS 7, an interval that is generally believed to represent ice-free conditions on the continent according to the δ¹⁸O record (243–191 ka; Lisiecki and Raymo 2005). While it may seem problematic to interpret melting of ice in glaciofluvial sediments during a supposedly ice-free interval, units A and B could have been deposited during the cool MIS 7-d stage (within 1σ range) or during the end of MIS-8 (within 2σ range, 300–243 ka, Lisiecki and Raymo 2005). Moreover, terrestrial glaciation is not always in sync with marine isotope stages (Gibbard and Hughes 2021). Near the end of glaciations, large ice sheets may have acted independent of global climate (e.g., Jackson et al. 2017), and this provides additional uncertainty when correlating MIS intervals with our terrestrial optical ages. Moreover, some caution needs to be exercised with the optical ages. Although exponential plus linear dose responses have been observed for quartz elsewhere, in some cases, it seems apparent that the ages may underestimate the true depositional age of deposits when compared to other age information. But there have been other cases where ages are in good agreement (see the discussion in Hodder et al. 2023). More work needs to be done to fully understand the utility of optically dating quartz from the HBL that has exponential plus linear dose responses.

The Pit 1 organic-rich silts of unit C are interpreted as nonglacial low-energy floodplain or pond sediments. Wood from this unit is beyond the age of radiocarbon dating (>50 ka ¹⁴C BP, UOC-13391) and the stratigraphic position, directly above the MIS 7 optical age determined from unit B, means that unit C was deposited during either MIS 7 or MIS 5. The unfortunate human excavation of the lateral and upper contacts means that we are unable to determine an accurate relationship with the laterally equivalent till (unit D); however, it most likely correlates to nonglacial deposition after the unit that was optically dated to MIS 7.

The pollen assemblage from the organic-rich silts of Unit C is similar to other subtill fossil sites in the western HBL, with a predominance of boreal and herbaceous groups, along with wetland indicators (Netterville 1974; Nielsen et al. 1986; Dredge et al. 1990). A dominance of well-preserved conifer remains (leaves/bark/needles) confirms conifers were present locally, while Carex seeds suggest sedges were also present. Charcoal along with charred plant remains suggests relatively local fires, similar to present day. The rare presence of Pediastrum suggests a local aquatic depositional environment, which is supported by plant macrofossils Najas flexilis seeds and Potamogeton, both aquatic indicators of a limnic environment. These two plant macrofossils also require specific temperatures to grow: Najas flexilis is a warm climate indicator, and the current distribution suggests July temperatures of 16–17 °C, whereas Potamogeton requires July temps of no cooler than 13–14 °C (Hultén 1950; Lampinen and Lahti 2013; Väliranta et al. 2015). The relatively low organic content suggests high minerogenic input during the fossil interval, and hence a limnic depositional environment with minerogenic components imported from the catchment.

Quantitative paleoclimate reconstruction shows that summer temperatures were statistically indistinguishable from those of present day, whereas total annual precipitation was similar or possible higher than present-day values. Subtill fossil pollen records yielding similar-to-present temperatures and higher precipitation were noted in the eastern Hudson Bay Lowlands (Dalton et al. 2018). The aforementioned sites were constrained via optical dating of quartz grains to broadly the MIS 4 interval (Dalton et al. 2018), but this age is tenuous and difficult to reconcile with the presumed fall in global sea level during that time (documented in records of relative sea level as well as the O¹⁸ record, Sasaki et al. 2004; Lisiecki and Raymo 2005), which would have most likely corresponded to a fully glaciated HBL. A second subtill site in western Hudson Bay contained fossil pollen records yielding similar-to-present day temperatures and precipitation (Black Duck River, Dalton et al. 2016, 2022). This second site was constrained via optical dating of quartz grains to MIS 5e (Hodder et al. 2023).

The deglacial to nonglacial sediments are capped by three different tills; a darker coloured, highly overconsolidated, clayier till (1–4 m thick, unit D or P-D-1, Fig. 11), a darker coloured, relatively less consolidated, silty till (0.4 m thick, unit E-2, Fig. 5), and a lighter coloured, also less consolidated, sandier till (0–2 m thick, unit P-D-2, Figs. 13b and Fig. 15f). Lithofacies E-2 is conformable with, and transitions up into, glaciolacustrine rhythmites, indicating that the till must have been deposited near or at the ice margin during a re-advance. It is likely that the rhythmites were formed during a penultimate Lake Agassiz drainage event (c.f. Bennett et al. 2000; Gauthier et al. 2020b), before the west-trending advance of the Stephens Lake sublobe (deglacial phase of Gauthier et al. 2019). The upper till at Pit 2-west (P-D-2), interpreted to have been formed by south–southwest-trending ice flow, is parallel to the Little Limestone esker (Fig. 2) and correlates with the last major regionally extensive MIS 2 ice-flow phase (220°–250°, phase 8 of Gauthier et al. 2019). In addition, lithofacies P-D-2 is the dominant surface diamict in the study area, and unconformably overlies lithofacies P-D-1 or bedrock (Fig. 2, Gauthier et al. 2022b). As such, we can be reasonably sure that the upper lighter coloured, less consolidated, sandier till was deposited during MIS 2. The lower darker coloured, more consolidated, and more clay-rich till, which is the only till present at Pit 1 (unit D, Fig. 5), directly overlies sand with optical age estimates of 214 ± 19 ka (1σ MAM (115-20-200-OSL011)) and 214 ± 22 ka (1σ MAM (HBL22-1)). At both pits, the lowest till samples are interpreted to have formed by west-trending ice flow (Table 2). At Pit 2-east, the middle till is interpreted to have formed by northwest or southeast ice flow, and the upper till by south–southwest-trending ice flow. It is not clear where the southwest-trending striated boulder pavement at Pit 2-west fits into that record. Regionally, tills similar to lithofacies P-D-1 contain a record of complex ice flow orientations, indicating that the lithofacies records paleoglaciological changes over a long time period (Gauthier et al. 2022b). As such, we tentatively assign the lower till (unit D, P-D-1) to the MIS 6 glaciation (130–191 ka, Lisiecki and Raymo 2005).

The till-clast composition is similar for all tills, with some variations up-section (Fig. 17). HBB clast content varies from 59–70 ct.%, Precambrian undifferentiated granitoid clasts from 9%–22% and UGG from 14%–32%. Unlike most regional till samples (Dredge and McMartin 2011; Gauthier et al. 2022b), there were no obvious indicator clasts of eastern or northern provenance within these diamict samples (i.e., Kipalu Formation, Omarolluk Formation, marine shell fragments, nor Dubawnt Supergroup). The only systematic variation is an up-section increase in UGG in unit P-D-1, likely indicative of increased uptake of local bedrock. The underlying gravel at Pit 1 (unit A) is also similar in composition to the tills (Fig. 17), meaning that clast composition is unable to differentiate glacial and deglacial events.

This study examined the highest-exposed nonglacial sediments at two sites near the edge of the western HBL. In contrast to the three regional MIS 5- or 6-aged subtill nonglacial sediments at similar depths (Fig. 1, Hodder et al. 2023), these new sites have sands dated to 214 ± 19 ka (1σ MAM (115-20-200-OSL011)) and 214 ± 22 ka (1σ MAM (HBL22-1))—and hence are interpreted to have been deposited during MIS 7. Importantly, the highest-exposed nonglacial sediments at another section in far northeast Manitoba also have an MIS 7 optical age determination (191 ± 29 ka (1σ MAM, 112-19-630-OSL025), (Fig. 1, “Echoing River” section of Dredge et al. 1990). Despite being the first occurrence of a nonglacial bed at this section, the regional stratigraphic framework indicates that it belongs to a second, older, ice-free period (Hodder et al.,2024), highlighting the fragmented nature of the stratigraphic record. Hence, for at least three western HBL sites, nonglacial sediments near the surface belong to an older interglacial (MIS 7), and there is no preservation of the penultimate interglacial sediments (MIS 5) at these investigated sites. It is clear, then, that revisions to the existing HBL Quaternary stratigraphy (Skinner 1973; Nielsen et al. 1986; Dredge et al. 1990; Thorleifson et al. 1993; Dredge and McMartin 2011) are necessary. Our new data confirm the need for detailed multiparameter till-stratigraphy, geochronology, and paleoenvironmental studies (c.f. Dube-Loubert et al. 2013; Hodder et al.,2024) to avoid erroneously correlating all nonglacial sediments to the same interglacial. Preservation of these old sediments is best explained as the result of spatio-temporal variability in erosion and deposition at the ice-bed interface—termed the “subglacial bed mosaic” by Piotrowski et al. (2004). This paradigm shift, from equivalent-chronologic units (Dredge et al. 1990; Dredge and McMartin 2011) to patchy units of different chronologies (c.f. Gauthier et al. 2019; Hodder et al. 2,024), now allows us to re-examine the Quaternary record in more detail.

The new ages indicate that between MIS 6 and 2, the net amount of diamict deposition at the gravel pits was between 2.5 and 4.5 m. This is valuable information, as the few attempts to model glacial erosion beneath the Laurentide Ice Sheet are typically validated using present-day sediment thickness and palimpsest debris trains (e.g., Hudson Bay Paleozoic carbonates, Hildes et al. 2004; Melanson et al. 2013). As noted by Larson and Mooers (2005), carbonate-bearing tills are not continuous sheets (e.g., Dredge and Cowan 1989) but rather palimpsest within a subglacial bed mosaic. The net carbonate-dispersal pattern is complex in Manitoba, and surface tills locally contain a range of carbonate concentrations that relate to overprinting and inheritance during till transportation and deposition (McMartin et al. 2016; Gauthier et al. 2019; Gauthier 2023). We now show that all “surface tills” were not formed during the same event, requiring further caution before using surface dispersal patterns in the interior region of the Laurentide Ice Sheet. Ongoing work on the till-stratigraphic composition at depth (e.g., Hodder et al.,2024) will provide additional debris transport and remobilization parameters to further constrain new subglacial sediment transport models.

The studied gravel pits, at the edge of the western HBL, provide a new local Quaternary stratigraphic framework. They expose ∼6 m of deglacial gravels and sands, the top of which was deposited around 214 ± 22 ka (1σ MAM (samples 115-20-200-OSL011, HBL22-1)) using optical dating of quartz grains from sand. These sediments were likely deposited either during the latter end of the cool MIS 7-d stage (within 1σ range) or near the end of MIS-8 (within 2σ range). The preserved deglacial sediments are overlain by ∼4 m of interglacial silts with organic horizons at one site; deposited in a low-energy floodplain or pond environment, within a climate similar to, or warmer than, present day but with higher precipitation. These interglacial sediments are preserved near the surface in a flat area; human excavation precludes their upper stratigraphic context. The deglacial gravels and sand are overlain by a darker coloured, highly overconsolidated, clayier till that was deposited during west, then northwest or southeast ice flow, followed by southwest-trending ice flow; possibly during MIS 6. A lighter coloured, relatively less consolidated, sandier till was later deposited during southwest-trending ice flow; likely during MIS 2.

If the optical ages approximate time of deposition, which we think they do, then the chronostratigraphy at these pits extends over multiple glaciations, despite being an 11 m thick stratigraphic package. As such, revisions must be made to the existing HBL Quaternary stratigraphy, given that previous workers correlated all organic-bearing sorted-sediment units to MIS 5. Moving forward, it is important to note that nonglacial “marker beds” may be locally or regionally absent and correlative to either MIS 5 or MIS 7, or older. Other beds, such as tills, will need to be used to develop local-scale stratigraphic frameworks. Providing a robust absolute chronology is crucial for correlating the complex and patchy sediments across the HBL. Optical dating provides this opportunity, but its practicability is limited by the long machine time required for samples of this age (weeks per sample) to perform adequate preliminary tests and analyses, and then measure an adequate number aliquots per sample from, which a representative D_e can be determined for the age calculation. Machine time would increase by one or two orders of magnitude if single grains are measured, as it is expected that only a small percentage of grains provide the signal measured from multigrain aliquots. But it is single grain work that would likely provide the most robust age information for HBL stratigraphy, so it should be attempted in the future for at lease some key samples.

More generally, this study provides chronologic evidence that nonglacial and penultimate glacial sediments are preserved at or near surface beneath parts of the former polythermal-based Laurentide Ice Sheet in the HBL. At these two sites, nonglacial sediments near the surface belong to an older interglacial (MIS 7), and there is little to no preservation of the penultimate interglacial sediments (MIS 5). Preservation of sediments of this age is best explained as the result of spatio-temporal variability in erosion and deposition within the subglacial bed mosaic. Future research on paleo-ice sheet regions must consider that the landscape (2D and 3D) can consist of patchy “units” of different chronologies.

Assistance in the Luminescence laboratory was provided by Nicola Ferguson, Justine Stoeckly, and Alison Goeres; Vanessa Brewer assisted in the field. Two anonymous reviewers are thanked for their helpful reviews.

Data generated or analyzed during this study are provided in full within the published article, its Supplementary material, and Appendix A.

Conceptualization: MSG, TJH

Funding acquisition: MSG

Investigation: MSG, TJH, ASD, OL, MS, MV

Methodology: MSG, TJH

Project administration: MSG

Resources: OL, SAF

Supervision: MSG, SAF

Visualization: MSG, ASD

Writing – original draft: MSG, TJH, ASD, OL

Writing – review & editing: MSG, TJH, ASD, MR

This research was funded by the Manitoba Geological Survey and the Geological Survey of Canada's Geo-mapping for Energy and Minerals program (GEM-GeoNorth). The Luminescence Dating Laboratory at University of the Fraser Valley (UFV) is further supported in part by Discovery grants and Research tools and Instruments grants from the Natural Science and Engineering Research Council (NSERC) of Canada, and from internal UFV grants.

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