Was there a nonglacial episode in the western Hudson Bay Lowland during Marine Isotope Stage 3?

Abstract

Establishing the timing of glacial and nonglacial intervals in the core regions of the Laurentide Ice Sheet (LIS) is essential to constrain ice-sheet configuration at times of globally reduced ice volume, such as during Marine Isotope Stage (MIS) 3 (~57–29 ka). Hudson Bay Lowland (HBL) deglaciation, at the centre of the LIS, has been inferred at MIS 3 based on near-infinite wood radiocarbon ages and limited luminescence ages. To better constrain the age of the penultimate deglaciation of the western HBL, this study initially identified the youngest intertill nonglacial sediments, based on extensive fieldwork and till characterization. Next, vetted radiocarbon ages were combined with revised stratigraphy to show that five previously identified “MIS 3” sites were likely deposited during an earlier ice-free period. Finally, new optical ages targeting the youngest intertill nonglacial bed at three localities yielded ages ranging from 166 to 146 ka; all older than MIS 3. These ages indicate that the penultimate deglaciation of the western HBL likely occurred during MIS 5e. This interpretation better explains accompanying paleobotanic data sets that indicate vegetation similar to vegetation existing under present interglacial conditions. Currently no firm evidence exists in the terrestrial stratigraphic record for the deglaciation of the western HBL during MIS 3.

INTRODUCTION

The Laurentide Ice Sheet (LIS) repeatedly waxed and waned over the North American continent during the Quaternary period. Considerable progress has been made in understanding ice-sheet configuration and evolution during and after the last glacial maximum (LGM; e.g., Dalton et al., 2020). Yet large uncertainties persist, especially for earlier stages of the late Pleistocene glaciation. Knowing the extent of former ice sheets at different times is critical for understanding long-term behavior of ice sheets, and this requires evidence from the geological record. For instance, mapping of ice-margin positions provides critical domain boundaries in ice-sheet models, which are used to get insights into past ice-sheet evolution that in turn affect the global ice volume budget used in modelling exercises (e.g., Pico et al., 2018; Gowan et al., 2021). One key time period that requires robust field-based evidence to reconstruct ice margins is the Marine Isotope Stage (MIS) 3 interstadial, spanning the time interval of ~57–29 ka (Lisiecki and Raymo, 2005). During this interstadial, the extent and volume of the LIS was restricted compared with the MIS 4 (71–57 ka) and MIS 2 (29–14 ka) stadials (Stokes et al., 2012; Dalton et al., 2019, 2022b; Gowan et al., 2021), but the exact configuration of the LIS geometry during this time period is uncertain, with large differences between estimated minimal and maximal extents (Dredge and Thorleifson, 1987; Batchelor et al., 2019; Gowan et al., 2021). Hence, there is a need for robust geochronological constraints that are supported by detailed lithostratigraphic and palaeoecological studies to elucidate the LIS geometry during MIS 3, as well as other periods of relatively reduced ice volume.

Data from the Hudson Bay Lowland (HBL), which is situated approximately in the geographic centre of the LIS between two major ice dispersal centres (Fig. 1), have long been used to argue for or against MIS 3 deglaciation in central Canada (Andrews et al., 1983; Dredge and Thorleifson, 1987; Dalton et al., 2019; Miller and Andrews, 2019). Initial interpretations, using relative dating methods, assumed that a mapped subtill nonglacial stratigraphic bed (Bell Sea sediments; Shilts, 1982; Andrews et al., 1983) was deposited during the last interglacial period around MIS 5e (~130–115 ka, peak highstand at 123 ka; Lisiecki and Raymo, 2005). The first direct dating of intertill sorted sediments in the western HBL using thermoluminescence (TL) age determination resulted in four MIS 3 ages and one MIS 5 age (Berger and Nielsen, 1991); however, the MIS 3 ages conflicted with a nonfinite radiocarbon age from wood within the same stratigraphic bed (Dredge et al., 1990; Fig. 1, site 2). More recently, a deglaciated HBL during MIS 3 has been interpreted using newer optical and near-finite radiocarbon ages (Dalton et al., 2016, 2019), but this was contested by Miller and Andrews (2019). The latter argued that only two radiocarbon ages from the HBL provided in Dalton et al. (2019) could plausibly be considered as finite radiocarbon ages compared to a significant amount of nonfinite radiocarbon ages. Furthermore, they argued that the timing of Heinrich events H4 (~38 ka) and H5 (~46 ka) provide firm evidence for an extensive LIS during this time interval, which is further supported in the terrestrial record of northeastern Canada (Edwards et al., 2022). Recent studies of the marine record in the Labrador Sea provide further support for an extensive LIS during MIS 3 (Parker et al., 2022).

Figure 1. Geochronological data for Quaternary intertill nonglacial deposits in the western Hudson Bay Lowland region of central Canada. New ages from this study are green, bolded, and italicized. Sites 1–7 were interpreted as possible Marine Isotope Stage (MIS) 3-aged deposits (Dalton et al., 2019), and readers are referred to the Supplementary Material 1 for additional details about site-specific geochronological data and references. Some sites are shown without data, due to vetting of potentially unreliable radiocarbon ages. The inset map shows the location of the study area in central Canada, relative to deglacial position of the Keewatin (K) and Quebec-Labrador (Q-L) Domes of the Laurentide Ice Sheet (LIS) and the Batchelor et al. (2019) minimum and maximum reconstruction limits for the LIS during MIS 3 at 45 ka. IRSL, infrared-stimulated luminescence; OSL, optically stimulated luminescence; TL, thermoluminescence.

The western HBL is a key region in the heart of the LIS, with an extensive pre-LGM stratigraphic record (Nielsen et al., 1986; Dredge et al., 1990; Dalton et al., 2018). The stratigraphic record is fragmented (Gauthier et al., 2019) and contains organic-bearing sediment deposited during at least two different ice-free interstadials and/or interglacial periods (Nielsen et al., 1986; Dredge and McMartin, 2011). As such, detailed stratigraphic studies are necessary to identify the youngest intertill nonglacial sorted sediment bed(s) in the region. Following extensive stratigraphic work in the western HBL (Gauthier et al. 2019 and references therein, 2021, 2022; Hodder and Gauthier, 2019, 2021), we document a series of new ages, from both optical and radiocarbon dating methods, on the uppermost intertill organic-bearing and sorted sediments in the western HBL. We show that key stratigraphic beds, including some that were previously assigned to MIS 3, can now be more confidently assigned to MIS 5e.

METHODS

Field data

Sediment exposures investigated were first cleared of slumped detritus, exposing a continuous section with no gap, and then described in detail. This includes describing the grain size, structure, colour, thickness, lateral variations, composition (clast lithology and geochemistry), and ice-flow data of identified stratigraphic beds. Ice-flow data were obtained from studied sections by measuring the long-axis orientations, or fabric, of clasts within till. Elongate clasts, defined by a minimum 1.5 ratio of the a-axis (longest) to the b-axis (middle), will, in most cases, rotate within the till matrix and orient parallel to the direction of shear that the overriding glacier exerts on the till (Holmes, 1941; Hicock et al., 1996). A minimum of 30 elongated clasts were measured at each of the nine clast-fabric sites, and the data set is provided in Supplementary Material 2.

This work highlights three new sites where the stratigraphically youngest intertill nonglacial bed was identified through regional mapping and the sediment was suitable for optical dating (Fig. 1, sites 15, 20, and 23). It also includes revision of two known sites, with the support of new stratigraphic correlations and radiocarbon ages from wood, where nonglacial beds were previously assigned to MIS 3 (Fig. 1, sites 3 and 4).

Optical dating

At sites 15, 20, and 23 (Fig. 1), the stratigraphically youngest intertill sorted sediments were sampled for optical age determination. Optical dating samples were collected by inserting opaque plastic tubes (~20 cm long and 5 cm in diameter) into the sediment faces, extracting them, and then sealing each end to preserve water content. Sample tubes were opened under appropriate laboratory lighting, and the light-exposed material at the ends of the tubes was removed and discarded. Bulk sediment extracted from the middle of the tubes was treated with HCl acid to remove any carbonates, and then by concentrated H₂O₂ to remove any disseminated organic matter. Sediment in the 180–250-μm-diameter size range was then extracted by wet sieving, and the quartz fraction was separated by density separation in lithium metatungstate. Quartz grains were then etched with concentrated HF acid to remove the outer alpha affected surfaces of the grains and to dissolve any remaining feldspar (e.g., Wintle, 1997). Aliquots consisted of 50–100 quartz grains mounted on aluminum disks using silicon oil as an adhesive.

Equivalent dose (D _e) values were found using the single-aliquot regenerative-dose (SAR) method (Murray and Wintle 2000, 2003). Irradiations and luminescence measurements were made using a Risø TL/OSL DA-20 reader. Laboratory irradiations were applied using a calibrated ⁹⁰Sr/⁹⁰Y beta-source on the reader that delivered ~4.7 Gy/min to the sample grains. Aliquots were stimulated by exposing them to 56 mW/cm² of blue (470 ± 20 nm) light for 100 s while being held at 125°C (Murray and Wintle, 2003). Ultraviolet emissions (~350 nm) were detected by an Electron Tubes Ltd. 9235QB photomultiplier tube placed behind a 7.5-mm-thick Hoya U-340 optical filter. An infrared “wash” (exposure to ~130 mW/cm² of 870 ± 40 nm light for 100 s while at 50°C) was administered before stimulation with blue light to reduce or eliminate the signal from any contaminating feldspar (Olley et al., 2004; Wintle and Murray, 2006). Aliquots were accepted for SAR analysis if their natural signals showed a clear “fast” component. Dose–response data were generated by giving the aliquots various laboratory doses, up to ~400 Gy, and one of the doses was repeated to test the effectiveness of the chosen test dose to correct for sensitivity change incurred during the SAR cycles. A preheat temperature and duration of 220°C for 10 s was chosen, which is about midway in the range of typical values (160–300°C; Wintle and Murray, 2006) used by most practitioners for quartz. The preheat treatment administered after the test dose involved heating the aliquot to 160°C at a rate of 5°C/s and then turning off the heat. A “zero-dose” measurement was also included to measure the degree of thermal transfer (recuperation). For each aliquot, the luminescence measured over the last 20 s of 100 s of stimulation was subtracted from that recorded over the first 0.4 s of stimulation, and this was divided by that measured after the test dose. These data were plotted as a function of laboratory dose and were best fit with saturating exponential+linear curves. The luminescence intensity measured from a preheated aliquot of prepared grains that had received no laboratory dose (the “natural”) was interpolated onto its dose–response curve, and a D _e value was measured from the dose axis; a representative example is show in Figure 2. Standard SAR quality control tests were applied, and these included assuring that the intensities of the recuperation signals were less than 5% of those of the natural signals and that recycling ratios were no more than ±10% of unity (Wintle and Murray 2006). Finally, the efficacy of the chosen SAR protocol was tested for each sample by means of dose-recovery tests (Wintle and Murray, 2006). For these, 15 aliquots were bleached by natural sunlight for 30 min, then given a laboratory dose approximately equal to the samples’ D _e, and then the same SAR protocol was applied to see if that dose could be recovered.

Figure 2. Luminescence decay curve (main graph) showing the “natural” signal and a dose–response curve (inset) for sample 112-19-605-OSL-004. These curves are typical for the three samples analyzed in this study. Note that the initial part of the luminescence decay is dominated by the desired “fast” signal component. In each case, the final 20 s of the signal was subtracted from the initial 0.4 s, and this value (L_i) was divided by that measured from a subsequent test dose (T_i) to produce the normalized (i.e., sensitivity-corrected) signal, which is plotted on the vertical axis of the dose–response curve graph. Note that the dose response is best fit by an exponential+linear function. The equivalent dose (D _e) is estimated by interpolation of the natural signal onto the dose–response curve as shown.

Environmental dose rates (Table 1) were determined by drying and milling a representative fraction of the bulk sediment used for dating and analyzing for U, Th, and K by neutron activation analysis (Table 1). As all our samples consisted of well-drained sand, the as-collected water content was used in the dose-rate calculations, with a ±10% uncertainty (1σ) that is intended to account for moisture fluctuation over time. Dose rates were calculated using standard formulae (Aitken, 1985; Berger, 1988; Lian et al., 1995) and the dose-rate scaling and conversion factors of Guérin et al. (2011) and Brennan (2003), respectively. The contribution from cosmic ray radiation (Table 1) was estimated using present burial depths and the relationship of Prescott and Hutton (1994).

Table 1. Optical age sample water content, radioisotope concentrations, sample depths, and calculated dose rates.

^a Δ^w = water content (mass water/mass dry sediment). For the dose-rate calculations, each Δ^w value included an uncertainty of ±10% (1σ) to account for fluctuations in water content over time.

^b d = sample depth beneath the ground surface.

^c Ḋ _C = cosmic ray dose rate; Ḋ _T = total dose rate (that due to cosmic rays plus that due to β and γ radiation).

Radiocarbon dating

Radiocarbon ages determined by accelerator mass spectrometry (AMS) were obtained from detrital wood (n = 8) and carbonate shell material (n = 5). Wood specimens were physically pretreated before analysis at the A.E. Lalonde AMS Laboratory or Keck Carbon Cycle AMS Laboratory at the University of California by removing the outer layer of the specimen and then chemically pretreated with an acid–alkali–acid wash. Carbonate shell specimens were physically pretreated before analysis at the A.E. Lalonde AMS Laboratory by removing porous or recrystallized areas with a Dremel bit and then chemically pretreated by etching with 0.2 N HCl to remove the outer 20–30% (Crann et al., 2017). Radiocarbon ages are interpreted as minimum limiting ages for the stratigraphic beds from which they were sampled.

Vetting of radiocarbon ages

Given the geochronological inconsistencies in the western HBL, we first vetted all new and previous radiocarbon ages for intertill sediments across the study area (Fig. 1, Supplementary Material 1), building upon Dalton et al. (2019). This vetting needs to be completed before interpretation to remove unreliable ages. Wood is the best material to attempt radiocarbon age determination when the material is near the upper limit of this geochronometer (e.g., Dalton et al., 2016; Bayliss and Marshall, 2019), and all vetted radiocarbon constraints for the western HBL are from wood (Fig. 1). Like others (e.g., Young et al. 2021), we removed all ages from bulk organic samples (n = 6), because this material can more easily be contaminated by inclusion of older carbon (old detrital fragments) and/or by younger carbon (modern-day rootlets). For example, site 6 (Fig. 1) was included as a potential MIS 3 site in Dalton et al. (2019) based on an untreated bulk organic age of 37 ± 1.6 ¹⁴C ka BP (WAT-1378). Yet there are another four nonfinite ages from pretreated bulk organics collected during this initial study and a nonfinite AMS age from wood (>50.8 ¹⁴C ka BP, ISGS A1658) at the same site (Supplementary Material 1). Given that only one of the six ages is finite, and the source material is the least reliable of the group, this site should not be used to support MIS 3 deglaciation. We also removed the finite shell ages that are >40 ka (n = 4), because shell material is susceptible to some isotopic exchange with the surrounding environment (e.g., Ramsey, 2008). This issue is not always minimized with physical and acid-etching pretreatment (method used in this study), and more advanced pretreatment procedures are necessary to gain confidence in shell radiocarbon ages (e.g., Douka et al., 2010).

RESULTS

Assessing the geochronological constraints

Ideally, multiple lines of geochronological evidence should be used to support any MIS 3 (considered to span ~57–29 ka) age determination without dispute. The vetted compilation of new and existing radiocarbon ages on sub-till wood shows only 1 out of 14 sites has finite radiocarbon ages; that one site also contains wood with a nonfinite radiocarbon age, indicating that these sediments are potentially older than MIS 3 (Fig. 1, site 5, Supplementary Material 1). Sediments with MIS 3 TL ages at sites 2 and 3 also contained wood with conflicting nonfinite radiocarbon ages (Fig. 1, Supplementary Material 1). Finding nonfinite radiocarbon age wood at two different sites decreases the likelihood of “transported” wood (i.e., reworked older wood). As pointed out by Roy (1998), the TL ages at sites 1–4 are likely underestimations due to the now better understood effect of anomalous fading on the luminescence signal of feldspars. Berger and Nielsen (1991) used long-term storage of their laboratory irradiated samples, some at an elevated temperature, in an attempt to reduce the effect of anomalous fading of the TL signal, but it has been found that this treatment is, at least for some samples, not sufficient (cf. Spooner 1992, 1993). Roy (1998) also questions these TL ages on the basis of the range of glow curve temperatures from which the TL signal used to build dose–response curves was determined and notes that inclusion of the “unstable” 270°C TL peak could have resulted in the calculated TL age being too young. This might be the case, as others have found that the infrared-stimulated luminescence (IRSL) signal associated with the 270°C TL peak is thermally unstable but are not certain that this is also the case for the TL signal from this peak (Duller and Bøtter-Jensen, 1993). Another point of concern are the dose rates reported by Berger and Nielsen (1991) for their samples, which are up to a factor of two higher than those reported by Roy (1998) for sediments extracted from the same stratigraphic bed. Indeed, one of our samples (not part of this study) collected from the same stratigraphic bed as the apparent MIS 3 age of Berger and Nielsen (1991) has radioisotope concentrations close to those of Roy (1998). These values are much lower than those presented in Berger and Nielsen (1991). If a dose rate that is the average of those of Roy (1998) is applied to the equivalent dose values of Berger and Nielsen (1991), the ages will almost double, and they would increase even more if corrected for anomalous fading. For these reasons, we believe that the previously reported TL ages of Berger and Nielsen (1991) are likely to be underestimations of the true depositional age, and we suggest that the ages should be removed from the MIS 3 age data set.

New optical ages

Three samples were collected for optical age determination (Fig. 1, sites 15, 20, and 23). For each sample, all the aliquots passed the SAR quality-control tests, but for some (5 to 10 aliquots per sample) their natural signals plotted beyond the last data point of their respective dose–response curves, and these were excluded from further analysis. The D _e values from accepted aliquots were plotted on radial plots (Galbraith, 2010; Fig. 3), and overdispersion (OD) values were calculated (Table 2). Overdispersion is the variation between individual aliquot D _e values above and beyond that associated with analytical uncertainties. It has been demonstrated that for samples consisting of quartz grains that were mostly likely fully bleached, OD values from small multigrain aliquots (fewer than 100 grains per aliquot) are, on average, about 20% (Galbraith et al., 1999; Lian and Roberts, 2006; Jacobs and Roberts, 2007), although for some of these samples, OD values are higher than 30% (Arnold and Roberts, 2009). The convention therefore has been to use a weighted mean D _e found using the central age model (CAM) (Galbraith et al., 1999) for the age calculation for samples for which OD is 20% or less (Jacobs and Roberts, 2007). For samples with OD values higher than 20%, which are assumed to consist of grains that have been heterogeneously bleached before deposition, the minimum age model (MAM) is appropriate. In that case the MAM D _e would give an approximation of true depositional age of a sample if a significant number of grains had been fully bleached, or it would give an age closer to the true age than that derived using the CAM if this is not the case. Experiments have shown that for some samples, OD measured from D _e values derived from large aliquots (~80 grains) is less than that derived from small aliquots (<10 grains), and much less than that derived from aliquots that contain single grains (Galbraith et al., 2005; Arnold and Roberts, 2009). These variations are likely due to averaging of various single-grain signals in multigrain aliquots. These observations, together with the fact that the sedimentology at our sample sites indicates that heterogeneous bleaching was highly probable before deposition and burial, suggest that the variation in D _e values of the individual grains in our samples could be higher that that suggested by our multigrain OD values. As such, we have provided ages using both CAM and MAM (Table 2), but we suggest that the MAM ages should be used for our stratigraphic interpretations due to the nature of the depositional environments. The ages derived using the MAM are 146 ± 20 ka (site 15, sample 112-19-629-024), 166 ± 12 ka (site 20, sample 112-16-407-019), and 150 ± 22 ka (site 23, sample 112-19-605-004). Results of the dose-recovery experiments (Table 2) suggest that the chosen SAR protocol was sufficient for producing accurate D _e and age values. In consideration of the inferred depositional environments for the units from which our optical dating samples were collected, these MAM ages should be considered to be upper limits. But because all three ages are statistically consistent (within 1σ of each other), the true age probably lies somewhere within the range indicated by their analytical uncertainties.

Figure 3. Radial plots (Galbraith 2010) showing the distribution of equivalent dose (D _e) values for each sample and estimations of representative D _e values (weighted mean) used for dating, using either the central age model (CAM) or the minimum age model (MAM; solid line that intersects the curved axis). Values that plot in the shaded zones fall within 2σ of the weighted mean value.

Table 2. Number of aliquots measured, overdispersion values (OD), equivalent dose values (D _e), optical ages, and dose-recovery (DR) ratios.

Note: OD and D _e values calculated using the software of Liang and Forman (2019).

^a N = number of aliquots accepted/number of aliquots measured.

^b C determined using the central age model (CAM).

^c M determined using the three-parameter minimum age model (MAM).

^dDR = dose-recovery ratio (given dose/recovered dose).

Stratigraphy and age constraints

We highlight herein three sections that contain sediment suitable for optical dating and that our till-stratigraphic framework has identified as lying directly under an advance-phase till from the last glaciation. Two other sections then show the impact till stratigraphy can have on understanding the relative age of intertill nonglacial sediments.

Site 15

The exposed stratigraphy at site 15 consists of five beds (Fig. 4). The youngest, uppermost bed (Bed A) consists of horizontally bedded sand and gravel with minor clay (Fig. 4). Shell fragments were observed throughout the upper sediments and paired whole valves of Hiatella arctica occur between 4 and 5 m depths, indicative of deposition within a marine environment. The horizontal bedding and grain-size variations in Bed A are interpreted to represent postglacial marine transgression at this site. A 7.8-m-thick layer of brown, massive, matrix-supported diamicton (Fig. 4, Bed B) sharply underlies Bed A. The diamicton contains 10–15% clasts and has a sandy-silt matrix (27% sand, 56% silt, and 13% clay). The diamicton has a crumbly appearance with minor oxidation rind developed on joint surfaces. A 0.6-m-thick layer of matrix-supported, massive, sandy gravel (Bed C) sharply underlies Bed B. The gravel is poorly sorted and contains 40% clasts that are granule- to large pebble-sized. The lower 0.2 m of this gravel contains rip-up clasts of diamicton. This gravel is interpreted as a glaciofluvial or fluvial deposit based on the poorly sorted nature of the sediments. A 1.6-m-thick layer of well-sorted sand sharply underlies the gravel (Bed D). The sand is horizontally bedded, with bedding defined by textural differences between alternating beds of very fine sand and beds of medium- to coarse-grained sand (Bed D). Shell fragments were observed in the upper 0.1 m of the bedded sands. The sands were likely deposited in a shallow-marine environment, based on the rhythmic bedding, the low elevation (41 m asl), the position beneath postglacial Tyrrell Sea sediments, and the proximity to the modern shoreline (~90 km). More than 6.0 m of dark greyish-brown, massive, matrix-supported diamicton (Fig. 4, Bed E) sharply underlies Bed D sand. The diamicton contains 15% clasts and has a clayey sandy-silt matrix (25% sand, 54% silt, and 21% clay). The diamicton has a blocky appearance with oxidation rind on the joint surfaces and is relatively denser compared with the upper diamicton (Fig. 4). The massive structure, texture, typical glacigenic shape of the clasts (bullet-shaped, faceted, and striated), the modality and strength of the pebble fabric, and lateral continuity indicate that Bed B and E diamicton is subglacial traction till (Evans et al., 2006). The stratigraphy at site 15 is interpreted to consist of postglacial marine sediments that are underlain by two till beds separated by intertill sorted sediments. The interpretation of a marine depositional environment for Bed D requires that Hudson Bay was ice-free and that relative sea level was higher than the present day. The age of the bedded, well-sorted sand (Bed D) was determined to be 146 ± 20 ka (Figs. 3 and 4, Table 2).

Figure 4. Quaternary stratigraphy and sediment photos of site 15. Till fabric data are plotted on equal area, lower hemisphere projection stereonets, with the principal eigenvalue (S₁) and eigenvector (V₁) values provided. Contours were generated using the Kamb method (Vollmer 1995). The sample location for optical age determination is shown by the black circle with an “L.”

Site 20

The exposed stratigraphy at site 20 consists of four beds (Fig. 5). The uppermost stratigraphic bed (Bed A) consists of 1.0 m of massive, moderately to well-sorted sand that is interpreted to have been deposited in a postglacial subaqueous environment, which was likely Lake Agassiz, as this site is above the local Holocene marine limit (~145 m asl; Gauthier et al., 2020) and the region was inundated by Lake Agassiz following deglaciation. Bed A sharply overlies 9.0 m of dark yellowish-brown, massive, matrix-supported diamicton (Bed B). The diamicton contains 10% clasts and has a sandy clayey-silt matrix (21% sand, 51% silt, and 28% clay). Bed B is in sharp contact with 1.4 m of grey, massive, matrix-supported diamicton that contains rip-up mud clasts (Bed C). The diamicton contains 10% clasts and has a silty-clay matrix (5% sand, 35% silt, and 60% clay). The massive structure, texture, and glacigenic shape of the clasts (bullet-shaped, faceted, and striated), the modality and strength of the pebble fabric, and lateral continuity indicate that Bed B and C diamicton is subglacial traction till (Evans et al., 2006). Considering the fine-grained texture of Bed C and mud rip-up clasts, it is interpreted that the ice that deposited this bed likely advanced over fine-grained sediments. Tills B and C were likely deposited during the same ice-flow event, based on the similar ice-flow direction interpreted from till fabrics (Fig. 5). Underlying these two tills is 3.8+ m of well-sorted sand (Bed D). Within the upper 1.0 m of this bed, there is convolute bedding and dewatering structures at places laterally along the contact (Fig. 5). These structures could have formed syndepositional in a subaqueous environment or from postdepositional dewatering. Beyond 1.0 m depth into this bed, the sand is cross-bedded to ripple cross-laminated sand. This thick sequence of cross-bedded sands indicates unidirectional bedload transport within a low-flow regime. No evidence was observed for settling from a suspension cloud. The stratigraphy at site 20 is interpreted to consist of well-sorted, cross-bedded sand deposited during a nonglacial period that is overlain by glacial and postglacial sediments. The optical age of the rippled sands (Bed D) was determined to be 166 ± 12 ka (Figs. 3 and 5, Table 2). This optical age had a low overdispersion (15 ± 3%) and a unimodal aliquot D _e population, which supports that the quartz grains were likely well bleached before burial. However, it is nevertheless possible that more than one age population exists and is masked, as our aliquots consist of 50–100 grains.

Figure 5. Quaternary stratigraphy and sediment photos of site 20. Till fabric data are plotted on equal area, lower hemisphere projection stereonets, with the principal eigenvalue (S₁) and eigenvector (V₁) values provided. Contours were generated using the Kamb method (Vollmer 1995). The sample location for optical age determination is shown by the black circle with an “L.”

Site 23

The exposed stratigraphy at site 23 consists of four beds (Fig. 6). The youngest, uppermost bed (Bed A) consists of 1.5 m of horizontally bedded sand and gravel. Bed A is interpreted to have been deposited in a postglacial marine environment based on the sediment texture and regional setting, as this site is situated below the local Holocene marine limit (~145 m asl; Gauthier et al., 2020) and the top of the section is 26 m above the modern sea-level elevation. Bed A sharply overlies 1.4 m of brown to dark greyish-brown, massive, matrix-supported diamicton (Bed B). The diamicton contains 5% clasts and has a sandy clayey-silt matrix (24% sand, 51% silt, and 25% clay). The diamicton has a minor blocky structure with oxidation rind on the joint surfaces. This diamicton, like the others described at sites 15 and 20, has all the main characteristics of a subglacial traction till. Bed B sharply overlies massive silt and massive to horizontally bedded fine sand to silty fine sand that contains disseminated organic material (Bed C). The presence of organic material requires terrestrial input into the depositional environment. Thus, the sediments that were targeted for optical dating were likely transported within a fluvial environment before deposition, which should be favourable to sufficient light exposure before burial. Bed C sharply overlies 3.3+ m of poorly sorted, matrix-supported sandy gravel (Bed D). Bed D gravel contains granule- to small cobble-sized clasts and rare shell fragments and is interpreted to have been deposited in a subaqueous environment. The subtill sediments (Beds C and D) indicate an overall fining upward sequence, possibly deposited in a fluvial environment. The stratigraphy at site 23 is interpreted to consist of nonglacial sediments that are overlain by glacial and postglacial sediments. The optical age of the organic-bearing fine sand (Bed C) was determined to be 150 ± 22 ka (Figs. 3 and 6, Table 2).

Figure 6. Quaternary stratigraphy and sediment photos of site 23. Till fabric data are plotted on equal area, lower hemisphere projection stereonets, with the principal eigenvalue (S₁) and eigenvector (V₁) values provided. Contours were generated using the Kamb method (Vollmer 1995). The sample location for optical age determination is shown by the black circle with an “L.”

Till stratigraphy as an important tool

Sites 3 and 4 (Fig. 7) are located ~10–12 km apart, and the stratigraphy at both sites was previously interpreted to contain one correlative intertill nonglacial bed (Dredge and Nielsen, 1985; Dredge et al., 1990; Dredge and McMartin, 2011). At site 3, the Bed C gravel is discontinuous across the exposure and pinches out to a till–till contact between Beds B and D in places. Bed C gravel contains an in situ peat bed with large (>10 cm) wood fragments that yielded a nonfinite radiocarbon age (>55 ¹⁴C ka BP, UOC-10974). Similarly, a detrital wood fragment from the gravel underlying the peat bed yielded a nonfinite radiocarbon age (>55 ¹⁴C ka BP, UOC-10975). Paleobotanic data from these sediments indicate a northern boreal forest to forest–tundra environment, similar to the present day for that area (Dredge et al., 1990). Hence, these uppermost intertill nonglacial sediments are interpreted to have been deposited during the last interglacial period (MIS 5e; Dredge et al., 1990), leading us to reject the MIS 3 TL age at site 3 (46 ± 4 ka; Berger and Nielsen, 1991) when conditions must have been colder based on the global sea-level record (Lisiecki and Raymo, 2005).

Figure 7. Quaternary stratigraphy at sites 3 and 4 in the on the Nelson River at Hudson Bay. The stratigraphic beds identified have been interpreted based on qualitative and quantitative observations. The thermoluminescence (TL) age estimations for stratigraphic beds (Berger and Nielsen, 1991) are interpreted to be age underestimations.

Bed C is missing between beds B and D at site 4 (Fig. 7). However, site 4 does expose nonglacial sediments below Bed D that have a TL age of 38 ± 3 ka (Berger and Nielsen, 1991; Fig. 1). Bed E consists of a clay-rich diamicton that contains whole valves of Macoma calcarea (Fig. 7), interpreted as a glaciomarine deposit. Underlying Bed E is massive, blue clay (Bed F-1), interpreted as a marine deposit. Owing to the principles of superposition and relative age, and perhaps ignoring the differences between till Beds B and D, Bed C at site 3 was originally correlated to Bed F-1 at site 4 (Dredge et al., 1990; Dredge and McMartin, 2011).

At site 3, we documented a newly discovered second, lower, intertill nonglacial bed, which indicates that there are actually two intertill nonglacial beds exposed in the stratigraphic record at site 3 (Fig. 7, Bed F-2). This lowermost intertill nonglacial bed consists of well-sorted, horizontally bedded sands (Bed F-2). To support correlations between sites 3 and 4, we examined the tills in addition to the sorted sediments. The upper till (Bed B) is significantly less dense, is sandier, and has a higher proportion of carbonate detritus compared with the middle till (Bed D), enabling stratigraphic correlation between these two sites (Supplementary Material 3). These differences, together with the knowledge that Bed C was deposited between the two tills, shows that the subtill marine sediments at site 4 are correlated to the second, older nonglacial unit at site 3 (Fig. 7, Bed F). Thus, while Bed C was deposited during MIS 5e, based on nonfinite radiocarbon ages, regional stratigraphy, and pollen data, the Bed F intertill sediments are even older (pre–MIS 5).

The reassessment of the stratigraphy at sites 3 and 4 indicates that the previously reported MIS 3 TL age from Bed F at site 4 (38 ± 3 ka), like the other TL ages in the region, is underestimated (Fig. 1). Furthermore, the stratigraphic evidence confirms that the near-finite radiocarbon age (49.9 ± 1.0 ka ¹⁴C yr BP, UOC-13465) from a Macoma calcarea whole valve in Bed E is not accurate and supports our approach of removing near-finite shell ¹⁴C ages from the data set. The stratigraphy here also highlights the importance of studying till beds, because the first bed of nonglacial sorted sediments encountered below till may not relate to the most recent ice-free period (i.e., site 4).

DISCUSSION

Optical dating

Doses of up to 400 Gy needed to be applied in the SAR sequences to be able to interpolate the natural quartz luminescence signals from a sufficient number of aliquots onto their dose–response curves. The high doses resolved a dose response that was best fit with a saturating exponential that includes a linear component (Fig. 2). An exponential+linear dose response has been observed for quartz in many other studies, but the reason for its existence and its efficacy for dating remains uncertain. Some workers have found that for samples expected to have ages >100 ka that exhibit exponential+linear dose responses, standard SAR protocols perform well and yield ages consistent with independent information (e.g., Murray et al., 2008; Pawley et al., 2008, 2010), while there have been other studies that show that exponential+linear dose responses produce ages that are too young when compared with independent information (see review by Lowick and Preusser, 2011). For many of these studies, however, reliable independent age information is lacking. A detailed study that examined the effectiveness of the SAR protocol for dating fine-grained quartz samples with exponential+linear dose responses has found that it performs well (e.g., Lowick and Preusser, 2011); it is therefore possible that the usefulness of the exponential+linear signal in quartz for dating is sample dependent. Nonetheless, more work clearly needs to be undertaken to better understand the underlying physical processes; until that is done, ages found from samples with exponential+linear dose responses should be considered with some caution (e.g., Li et al., 2015).

This study targeted intertill sorted sediment beds that were identified using a regional till stratigraphic framework to have been deposited during the penultimate deglaciation of the region. At face value, two of our optical ages are within 1σ of MIS 5e, which is an expected ice-free period in the HBL considering this was a period of time with significantly reduced global ice volume (Lisiecki and Raymo, 2005). The third age (166 ± 12 ka) is not within 1σ of MIS 5e, but based on the stratigraphy at this site, these sediments could have been deposited during the deglaciation phase of MIS 6. Furthermore, terrestrial glaciation is not always synchronous with marine isotope stages (Gibbard and Hughes, 2021), and it is therefore possible that the timing of glaciation in the HBL was slightly out of phase with the timing of climate shifts recorded in the marine record. It is likely that any overestimation of our optical ages is due to the inclusion of a sufficient number of grains whose signals were not adequately reset before burial.

K-feldspar is also routinely used for optical dating, and its upper age limit is usually higher than that for quartz but it was not used in this study, because its dose response was expected to be nonlinear in the region of our samples’ natural luminescence signals, and it therefore would not be possible to correct them for anomalous fading (Huntley and Lamothe, 2001). The post-IR IR (pIR-IR) signal from K-feldspar, which usually fades less (in some cases negligibly), was also not investigated, as its signal resets (bleaches) relatively slowly, and it was therefore anticipated that this would be a problem for the fluvial and lakeshore sedimentary environments sampled here. Another option for dating quartz in this region would be to investigate the thermally transferred optically stimulated luminescence (TT-OSL) signal, which has been observed to maintain a dose response defined by a single saturating exponential, even at high doses (e.g., Neudorf et al. 2019), but like the pIR-IR signal from feldspar, the TT-OSL signal from quartz bleaches much more slowly than the standard quartz signal used for dating (Jacobs et al. 2011; Neudorf et al. 2019), and it would therefore be expected to be of limited use for the sedimentary environments studied here; aeolian units would need to be identified. Dating aliquots that consist of single grains of quartz would be an advantage, as D _e populations would be easier to identify using the MAM, or the finite mixture model, which is suited to analysis of D _e values derived from single-grain aliquots, but can present serious problems when applied to multigrain aliquots (Arnold and Roberts 2009). However, routine SAR dating of single grains from the samples studied here is not practical, as the sensitivity of the quartz in the region is variable, and would require thousands of grains being analyzed. Moreover, the depositional age of the samples requires high-dose irradiations and therefore exceedingly long machine time. Nevertheless, single-grain dating of a few key strata in the region, including some where reliable independent age information is available, should be considered in future research.

Resolving the age of the Black Duck River nonglacial sediments

The organic-bearing nonglacial sediments at site 5 (Fig. 1) are described to be near the surface and overlain by ~1.5 m of gravel that is interpreted to have been deposited in a postglacial marine environment (Dalton et al., 2016). Three of the four radiocarbon ages on wood from the organic-bearing nonglacial sediments are near-infinite (45.7 ± 1.3 ¹⁴C ka BP, ISGS-A1656; 49.6 ± 1.0 ¹⁴C ka BP, UOC-0587; 50.1 ± 3.3 ¹⁴C ka BP, ISGS-A1995), which led to the interpretation that these sediments were deposited during MIS 3 (Dalton et al., 2016, 2019). Dalton et al. (2016) considered these ages reliable, based on the repeatability of the radiocarbon ages from different laboratories, and stratigraphic-age agreement. Importantly though, the paleobotanic data from the site 5 organic-bearing bed indicate a paleoenvironment dominated by boreal and peatland taxa, similar to vegetation currently present in the region (Dalton et al., 2022a). Quantitative paleoenvironmental reconstruction of these sediments indicates the mean summer temperature varied between 9.4 and 14.5°C compared with the present-day 11°C, and total precipitation was similar to the present-day precipitation (467 mm; Dalton et al., 2022a). This evidence of climatic conditions similar to those of the present day supports deposition during an interglacial period, such as MIS 5e, as it is difficult to reconcile a climate reconstruction similar to the present day with the LIS still covering large areas of North America in the most restrictive reconstructions for MIS 3 (e.g., Batchelor et al., 2019; Gowan et al., 2021). As such, it seems likely that the nonfinite radiocarbon age is more reliable than those near-finite from the same strata, and that deposition at site 5 occurred during MIS 5e.

This study attempted to obtain an optical age from the nonglacial sediments at the Black Duck River site 5 to test the accuracy of the near-finite radiocarbon ages from wood; however, the exact site could not be accurately located in the field based on published location information. While the stratigraphy at site 23 on the Black Duck River (Fig. 6) in this study is different from site 5, considering that both of these sections contain near-surface interglacial deposits and are situated within ~1 km of each other, we interpret that the organic-bearing Bed C at site 23 is likely correlated to the organic-bearing nonglacial bed at site 5 studied by Dalton et al. (2016). Owing to the susceptibility of near-infinite radiocarbon ages to modern contamination (e.g., Reyes et al., 2020), the presence of one nonfinite radiocarbon age, and paleoenvironmental information consistent with a climate warmer than that expected during MIS 3, we consider our site 23 to have more reliable stratigraphy and chronological data than site 5 (Dalton et al., 2016). Thus, we suggest that the optical age at site 23 (150 ± 22 ka), within 1σ error of MIS 5e, more accurately dates the deposition of interglacial sediments situated on the Black Duck River and is in better agreement with pollen profile from this bed at site 5.

HBL evidence for MIS 3 deglaciation?

Vetting of previous radiocarbon ages, together with new optical and radiocarbon ages, has led to removal of sites 1–6 from the MIS 3 age data set in the western HBL (Fig. 1, inset). That leaves site 7, where two MIS 3 optical ages were obtained from quartz grains (42.2 ± 4 and 52.5 ± 5 ka; Dalton et al., 2016), as well as an older TL age of 73 ± 10 ka (Forman et al., 1987) from subtill sediments interpreted to be marine. A third MIS 3 optical age provided in Dalton et al. (2016) is from the Drowning River near James Bay (~540 km to the southeast of site 7). At this site, the quartz optical age of 42.9 ± 3.7 ka is incompatible at 1σ error with two nonfinite AMS radiocarbon ages on wood (>48.4 ¹⁴C ka BP, ISGS 2463; >51.7 ¹⁴C ka BP, ISGS |A2425; Dalton et al., 2016). The inconsistency between the optical ages and nonfinite radiocarbon ages at this site raises questions about the accuracy of the optical ages in Dalton et al. (2016). Furthermore, Gao et al. (2020) report new chronological data from a site just 16 km away from the Dalton et al. (2016) Drowning River site, with an optical age from feldspar that is much older than MIS 3 (118 ± 13 ka); correlation of the two subtill units is unclear. The optical age provided by Gao et al. (2020) is likely a lower limit, as an age of this magnitude is expected to have a dose response that is nonlinear, and such ages cannot be corrected for anomalous fading using the method of Huntley and Lamothe (2001). However, Gao et al. (2020) do not provide representative dose–response curves necessary to fully evaluate the efficacy of their fading correction procedure. As such, site 7 will require additional investigation to confirm the sediments there are indeed MIS 3. For now, it remains the only site with an MIS 3 age in western HBL, in contrast to the great number of sites with nonfinite radiocarbon ages or older optical ages (Fig. 1).

Our new minimum limiting radiocarbon and optical ages for the stratigraphically youngest intertill nonglacial sediments in the western HBL indicate that these sediments were likely deposited during MIS 5e to late MIS 6. Hence, we find there is insufficient stratigraphic and geochronological evidence for models and reconstructions (e.g., Carlson et al., 2018; Gowan et al., 2021) to consider the western HBL region ice-free during MIS 3. Additional work is needed to pinpoint the southern margin of the LIS in central Canada during this time period.

As pointed out by Miller and Andrews (2019), the timing and provenance of large iceberg discharge events into the North Atlantic, which are recorded by Heinrich layers (e.g., Hemming, 2004), can provide an important proxy for ice extent and need to be reconciled with HBL age constraints. Heinrich events H4 (~38 ka) and H5 (~45 ka), which occurred during MIS 3, have a provenance signature similar to that of Heinrich event H2 (~24 ka), which has long been associated with the Hudson Strait Ice Stream, and a catchment zone extending deep into Hudson Bay and perhaps beyond (e.g., Hemming, 2004; Ross et al., 2011; Miller and Andrews, 2019). As pointed out by Roy et al. (2009), some uncertainties persist about the exact provenance of the sediment, but an ice sheet over Hudson Bay during MIS 3 feeding a long, topographically controlled Hudson Strait Ice Stream remains the most plausible model for the source of Heinrich layers H4 and H5. Configuration of a “big” ice sheet is also in agreement with the interpretation that the Hudson Strait Ice Stream was stable at the ice shelf during H4 and H5 (Rashid et al., 2019). Heinrich events thus provide additional independent support for a glaciation of Hudson Bay and the western HBL during MIS 3.

CONCLUSIONS

New age determinations from the uppermost intertill nonglacial sediments, combined with consideration of the paleobotanic data sets in the western HBL stratigraphic record, suggest that the region was last deglaciated during MIS 5e. New radiocarbon ages and vetting of the radiocarbon age data set, together with the new optical ages and detailed stratigraphic correlations herein, have removed six sites from the HBL “MIS 3” data set, leaving just one potential MIS 3 site that has conflicting age assignments.

The Quaternary stratigraphy of the HBL is highly fragmented, which hampers correlations and also indicates that previous investigations of the region have provided an oversimplified stratigraphic framework. This complex stratigraphic record, including detailed lithostratigraphy of enclosing glacial units, needs to be considered alongside geochronology results to provide robust reconstructions of Late Quaternary nonglacial events. This study suggests that the LIS southern margin did not retreat north or northeast of the western HBL and that Hudson Bay remained glaciated during MIS 3. Future studies to determine the maximal or minimal extent of the LIS during MIS 3 should continue reconstruction of the stratigraphic record south of the western HBL.