In situ cosmogenic ¹⁰Be and ²⁶Al reveal the complex exposure and erosion history of the landscape once covered by the Quebec-Labrador Ice Dome

The rate at which ice sheets erode rock and produce sediment is poorly known. Here, we use paired cosmogenic nuclides in both deglacial and modern sediment to understand better the efficacy with which the Quebec-Labrador Ice Dome (QLID) of the Laurentide Ice Sheet eroded bedrock and generated sand and boulders across the landscape of eastern Canada. We sampled deglacial sediment (esker and delta sand, n= 10), sediment from modern streams (n= 11), one bedrock outcrop, and a bedrock depth profile (n= 7), measuring concentrations of ¹⁰Be and ²⁶Al in quartz isolated from all samples. We also collated published cosmogenic nuclide measurements of boulders and bedrock from eastern Canada (n= 237 samples), and using independent estimates of deglaciation timing, calculated initial nuclide concentrations when the material was exposed by the most recent deglaciation, between 6.3 to 15.2 ka.

At the time of deposition, all 10 deglacial sand samples contained ¹⁰Be and ²⁶Al, on average equivalent to several thousand years of surface exposure. The ubiquitous presence of ¹⁰Be and ²⁶Al in eastern Quebec deglacial sediment is consistent with older-than-expected exposure ages for bedrock outcrops (n= 26 of 46 samples) and boulders (n= 65 of 192 samples) once covered by the QLID. Error-weighted averages of ²⁶Al $/$ ¹⁰Be ratios for both deglacial (6.1 ± 0.3, all uncertainties 1 SD) and modern sediment samples (6.6 ± 0.5) are lower than the measured production ratio at high latitudes (Greenland, 7.3 ± 0.3), suggesting cumulative burial of at least some sediment grains for at least hundreds of thousands of years.

This burial history suggests that ice at the middle of the QLID either survived some interglacials and/or that the average sediment residence time on the landscape is several times longer than a 100 kyr glacial cycle, allowing storage and burial of sediment over multiple glacial cycles, either under ice and/or in thick deposits such as deltas and moraines. Modern river sand contains on average only slightly higher nuclide concentrations than deglacial sediment, suggesting that contemporary river sand is predominately recycled from glacial deposits. Together, the new sediment data (which amalgamate across large areas of the landscape), and our compilation of bedrock and boulder point data, suggest that the average depth of bedrock erosion by ice and the speed of glacial sediment transport in eastern Canada were insufficient to remove material containing cosmogenic nuclides produced during prior interglacial(s).

1 Introduction

Ice sheets are important geomorphic agents of high-latitude landscape change (Patton et al., 2024; Sugden, 1978; Wilner et al., 2024) although the rates at which they erode landscapes and transport sediment are not well constrained. Available data suggest that both erosion and sediment transport rates are spatially and temporally variable (e.g., Alley et al., 2019; Cowton et al., 2012; Hallet et al., 1996; Patton et al., 2022). Rates of glacial erosion appear to be sensitive to the time-scale, making it difficult to accurately compare short-term, contemporary rates of glacial erosion and sediment export to those determined over longer time scales (Ganti et al., 2016; Wilner et al., 2024). Modelling sub-glacial processes, including rates of erosion and sediment transport, is not straightforward and results do not always match field evidence (e.g., Hildes et al., 2004).

During the Last Glacial Maximum (LGM), between 20–25 ka, the Laurentide Ice Sheet (LIS) covered most of Canada and the northern United States (Dalton et al., 2022). It was the largest body of ice in the Northern Hemisphere holding about 70 meters of sea level equivalent (Simms et al., 2019). Its ablation during the latest Pleistocene and early Holocene (characterized by collapse of northern Canadian ice domes) exposed a landscape that was once the ice sheet's bed including rocky outcrops shaped by abrasion and plucking, and deposits of glacial and post-glacial sediment (Occhietti et al., 2011). This relict landscape is an archive of past glacial and interglacial landscape behavior that can be used to untangle the history of landscape change (e.g., Bierman et al., 1999; Booth, 1994; Collins et al., 2025; Schaefer et al., 2016).

Quantitative geochronology techniques, primarily radiocarbon dating of organic material and cosmogenic nuclide exposure dating of bedrock and boulders, has established broad outlines of the post-LGM retreat history of the LIS (Dalton et al., 2023). Yet, uncertainties remain because dates from different chronometers do not always match and the concentrations of cosmogenic nuclides at the time of boulder and bedrock exposure and sediment deposition remain poorly constrained (Bierman et al., 1999, 2015; Briner et al., 2016; Halsted et al., 2024; Staiger et al., 2005). Paired nuclide cosmogenic burial dating has provided ages for pre-LGM LIS advances (Balco, 2005; Balco and Rovey, 2008, 2010).

Our study seeks to understand better the erosivity of LIS ice in eastern Canada by measuring and interpreting the concentrations of the cosmogenic nuclides ²⁶Al and ¹⁰Be in sediment transported by the LIS. To accomplish this goal, we sampled along two transects stretching outward from the former middle of the Quebec-Labrador Ice Dome (QLID) at Labrador City: one eastward toward the Atlantic coast at Goose Bay, Labrador, and another southward to the St. Lawrence River near Quebec City (Figs. 1 and 2). To provide context for our samples, we compiled published cosmogenic nuclide analyses of bedrock and glacially transported boulders from the QLID and used independent age estimates of deglaciation (Dalton et al., 2020) to estimate nuclide concentration at the time these sample sites were first exposed by LIS retreat.

Figure 1 Map of field area with sample sites for new data included in this paper. Marine sediment cores with sample analyses in LeBlanc et al. (2023) shown by open circles. Terrestrial samples with new data provided in this paper shown with symbol color and shape corresponding to sample type. Inset indicates location of field area. Grey lines indicate direction of ice flow (Margold et al., 2018). Blue deglacial isochrons from Dalton et al. (2020).

Figure 2 Photographs of representative field sites. (A) Modern river sediment at MC-03. (B) Deglacial sediment at SS-05. (C) Bedrock outcrop at GB-06. (D) Deglacial sediment at MC-01. (E) Bedrock depth profile (LBD 19 46-52) exposed in roadcut after 1985. Depths in cm shown on outcrop. Note sample at 10 cm depth is offset to left of profile.

The new data, and those we compiled, together allow us to constrain better the accuracy of cosmogenic exposure ages calculated for bedrock surfaces and erratic boulders in areas once covered by the QLID. Our work provides additional context for interpreting cosmogenic nuclide measurements made in sand-sized sediment of glacial and interglacial age present in marine sediment cores collected offshore of former and current ice sheets (e.g, Bierman et al., 2016; Christ et al., 2020; LeBlanc et al., 2023) as well sub-ice bedrock and sediment (Bierman et al., 2024a; Christ et al., 2021, 2023; Schaefer et al., 2016).

Study Site

The QLID occupied the eastern subarctic Canadian Shield. The area we focus on is underlain mostly by Proterozoic quartzofeldspathic gneisses and granites (Hynes and Rivers, 2010). Regolith and soil are commonly thin and patchy (< 2.5 m; Pelletier et al., 2016) and punctuated by prominent bedrock outcrops and glacial erratics (Ullman et al., 2016). Central and southern Quebec-Labrador also includes multiple moraine systems that mark and track the final deglaciation of the ice dome into the early Holocene (Couette et al., 2023; Occhietti et al., 2011; Ullman et al., 2016). There are eskers, and ice contact deltas where the ice front, during deglaciation, met a body of standing water (Liverman, 1997; Storrar et al., 2013). The paraglacial landscape we consider here is still experiencing post-glacial rebound, with the greatest amount having occurred near James Bay and southern Hudson Bay (250–270 m; Lajeunesse and Allard, 2005; Lavoie et al., 2012), compared to ∼ 100 m along coastal eastern Labrador (Andrews and Tyler, 2011).

Eastern Canada is dominated by boreal spruce forests, sedges, and muskegs (shallow bogs covered in moss), a sub-arctic ecosystem prone to burning during dry periods in the summer, with a recorded fire history stretching back to the 1950's (Payette et al., 1989). Northern Quebec and Labrador are classified under the Dfc climate zone (cool continental climate/subarctic) according to the Köppen climate classification system (Amani et al., 2019; Beck et al., 2018). During winter, ground-based measurements record a mean of ∼ 158 mm of snow water equivalent (SWE) for eastern Canadian boreal forests (Larue et al., 2017). The St. Lawrence River flows from southwest to northeast and is located southeast of the middle of the QLID. The Churchill River flows east from the former middle of the QLID, draining into the Atlantic Ocean.

2 Background and Unanswered Questions

2.1 Cosmogenic Nuclides as Recorders of Glacial History and Process

Since the late 1980s, measurement of cosmogenic nuclides, most commonly ¹⁰Be and ²⁶Al extracted from quartz, has helped elucidate the Quaternary history of glaciated landscapes and the rate of glacial processes. These nuclides are produced by cosmic-ray interaction with target atoms at and near Earth's surface. Near the surface, neutron spallation production is dominant, but deeper in a profile (> 2 m), muons dominate nuclide production at rates much lower than near the surface (cf. Braucher et al., 2013).

¹⁰Be and ²⁶Al are produced in both rock and glacially deposited sediment. Analysis of their concentrations can provide information about the timing of deglaciation (e.g., Corbett et al., 2011), the depth of glacial erosion (e.g., Briner et al., 2006; Miller et al., 2006; Stroeven et al., 2002), the persistence of glacial sediment on the landscape (e.g., Nelson et al., 2014), the thermal conditions at the bed (e.g., Staiger et al., 2006), and in the right setting, the limits on the extent and timing of sediment and rock burial by ice sheets (e.g., Bierman et al., 1999; Briner et al., 2006; Christ et al., 2021; Corbett et al., 2016b, c; Marsella et al., 2000; Miller et al., 2006; Schaefer et al., 2016). Measuring the concentrations of ¹⁰Be and ²⁶Al in the same sample can reveal periods of burial during or after near-surface exposure (Klein et al., 1986) because the nuclides decay at different rates; their half-lives differ by a factor of two (∼ 1.36 and ∼ 0.73 Myr, respectively).

Cosmogenic nuclide measurements made in detrital sediment from contemporary paraglacial channels (Nelson et al., 2014) and glacial deposits (e.g., Bierman et al., 2016; Goehring et al., 2010), have been used to infer the concentration of nuclides inherited from prior periods of exposure, to interpret ice sheet basal processes (Corbett et al., 2021), and to infer sediment sourcing (e.g., Balco, 2005; Bierman et al., 2016; Fame et al., 2018; Nelson et al., 2014). Cosmogenic nuclide measurements in marine sediments and in material collected from the base of ice cores allow inference of glacial process, ice sheet history, and interglacial landscape exposure (Bierman et al., 2016, 2024a; Blard et al., 2023; Christ et al., 2020, 2021, 2023; Schaefer et al., 2016; Shakun et al., 2018).

When a landscape is covered by a thick layer of ice, such as the LIS, or if sediment is stored in deposits thick enough to largely attenuate cosmic rays, production of in situ ²⁶Al and ¹⁰Be is diminished or stops. As nuclides produced during initial exposure decay, the ²⁶Al $/$ ¹⁰Be ratio decreases (Bierman et al., 1999), a change that is detectable after several hundred thousand years of burial (Granger, 2014; Granger and Muzikar, 2001). Re-exposure to cosmic rays at or near Earth's surface increases the ²⁶Al $/$ ¹⁰Be ratio; given sufficient time, it approaches the production ratio. However, if sediment remains buried meters below the surface during interglacial periods (such as in deglacial deltas, till sheets, and eskers or under ice), ²⁶Al and ¹⁰Be concentrations and the ²⁶Al $/$ ¹⁰Be at the time of deposition are preserved in whole (deep burial) or in part (shallow burial) (Bierman et al., 2016; Goehring et al., 2010; Nelson et al., 2014).

Concentrations of ¹⁰Be and ²⁶Al in glacial sediment therefore reflect the history of that sediment and of the ice sheet that eroded and transported it over time (e.g., Bierman et al., 2016; Briner et al., 2016; Corbett et al., 2021; Goehring et al., 2010; Nelson et al., 2014). Analysis of sediment deposited by retreating ice sheets allows inferences about erosion and transport efficiency at and near the former or present ice sheet bed. Long interglacial exposures, thin sediment cover, and bedrock that is resistant to erosion will allow high concentrations of nuclides to accumulate – for example, in central North America (Balco, 2005). On the other hand, persistent ice cover, high rates and depths of glacial erosion, and efficient sediment transport by ice will lower nuclide concentrations in glacially derived sediment, such as in parts of Greenland (Corbett et al., 2021; Goehring et al., 2010; Nelson et al., 2014).

2.2 Uncertain Laurentide Ice Sheet History in Eastern Canada

For most of the last glaciation, the LIS was characterized by three major ice domes, regions of especially thick (∼ 4 km in some places), outflowing ice: Foxe Baffin Dome, Keewatin Dome, and the QLID considered here (Stokes et al., 2012). Changes in LIS size are thought to have broadly tracked the marine oxygen isotope record, though uncertainties in ice volume and extent through time reflect the paucity of geologic constraints (Batchelor et al., 2019).

The extent and duration of ice cover in eastern Canada during the Pleistocene remains uncertain. For example, the majority of Canada may have been ice-covered, with the LIS reaching 70 % of its LGM extent as early as Marine Isotope Stage (MIS) 5d (∼ 110 ka) or not until MIS 4 (∼ 60 ka) (Dalton et al., 2022). During subsequent retreat into MIS 3 (∼ 45 ka), a combination of luminescence, ¹⁴C, and cosmogenic nuclide (¹⁰Be and ²⁶Al) dating, along with evidence of a marine incursion into Hudson Bay, was used to argue for deglaciation of the middle of the LIS (Dalton et al., 2019; Pico et al., 2018). Others disagreed, interpretating these ages differently and citing the timing of carbonate-rich Heinrich events H5 and H4, which indicate that an intact, marine-terminating Hudson Strait ice stream existed during MIS 3 (Gauthier et al., 2024; Hodder et al., 2023, 2024; Miller and Andrews, 2019). Recent measurements of cosmogenic nuclide concentrations in quartz-bearing ice rafted debris (IRD), sourced at least in part from LIS-generated icebergs, also suggest ice may have lingered in eastern Canada through many interglacial periods of the last My (LeBlanc et al., 2023). Evidence includes ²⁶Al $/$ ¹⁰Be ratios (4.7 ± 0.8) measured in glacial-marine sediment well below the high-latitude (Greenland) production value of 7.3 ± 0.3 (Corbett et al., 2017). Such depressed ratios require long (> 1 Myr) periods of burial (by ice, sediment, or water) throughout the Pleistocene.

Retreat of the LIS in Quebec and Labrador following the LGM is constrained by radiocarbon ages on post-glacial marine shells and lake organic material (Dalton et al., 2020; Dyke, 2004) as well as cosmogenic ages that are incorporated in Dalton et al. (2023), which we chose not to use here because of potential circularity because their retreat model includes these exposure ages. Together, these data suggest that the initial retreat was gradual, with the eastern margin pulling back slowly from the continental shelf and not reaching the modern coastline until ∼ 11.5 ka (Fig. 1). Retreat rate increased when Hudson Bay opened and the QLID separated from other sectors of the LIS after ∼ 8.5 ka. The margin of the ice sheet then retreated rapidly toward the middle of the landmass, disappearing entirely by ∼ 6 ka.

In situ produced cosmogenic nuclides have been used extensively to date retreat of the QLID after the LGM in Labrador and Quebec. Clark et al. (2003) report 27 measurements (including 2 replicates) of ¹⁰Be on 21 boulders and 4 bedrock surfaces, concluding that ice left valleys in the Torngat Mountains of northeastern Labrador between 12.0 and 13.4 ka. Further south, Ullman et al. (2016) used ¹⁰Be concentrations in 65 boulder samples to demonstrate that ice left eastern Labrador at 10.4 ka and the St. Lawrence coastline in southern Quebec at 9.2 ka, reaching the middle of the Quebec-Labrador landmass at 6.7 ka. Couette et al. (2023) interpreted 37 ¹⁰Be exposure ages on boulders to represent five still-stands or re-advances of the eastern LIS margin as it retreated inland from the coast of eastern Labrador (∼ 12.9, ∼ 11.5, ∼ 10.4, ∼ 9.3, and ∼ 8.4–8.2 ka). Carlson et al. (2007), using 13 boulder measurements, determined that deglaciation across central Quebec occurred rapidly between 8.0 and 6.8 ka. Along western Ungava Bay, Lefebvre-Fortier et al. (2024) obtained ¹⁰Be ages on 14 boulders from five marine strandlines, which suggest north-to-south ice margin retreat from 8.9 to 7.9 ka. Dubé-Loubert et al. (2018, 2021) placed deglaciation between 8.6 and 8.0 ka south of Ungava Bay based on 30 boulder ages from glaciomarine and ice-dammed lake delta surfaces. In the same region of northeastern Labrador, Rice et al. (2019) ¹⁰Be dated 8 bedrock and 2 boulder samples. Rejecting two much older ages, they report a range of ages from 9.8 to 6.2 ka.

Taken together, the cosmogenic data generally agree with Dalton et al.'s (2020) radiocarbon-based chronology of deglaciation in eastern Canada. The ¹⁰Be ages show a pattern of radial ice margin retreat toward the middle of Quebec-Labrador. However, many of the cosmogenic nuclide studies in this region of eastern Canada (Clark et al., 2003; Couette et al., 2023; Dubé-Loubert et al., 2018, 2021; Marquette et al., 2004; Rice et al., 2019) identify some older-than-expected ¹⁰Be ages, which they attribute to inheritance of cosmogenic nuclides from periods of prior exposure. Ullman et al. (2016) note that radiocarbon ages are less by centuries than the ¹⁰Be ages along that study's transects across the western, eastern, and southern part of the former QLID. This discrepancy could represent the pervasive presence of a low concentration of nuclides inherited from an earlier near-surface exposure, a lag in revegetation after deglaciation (cf., Halsted et al., 2024), and/or production rate calibration/scaling inaccuracies. We consider such inheritance in detail in this paper.

2.3 Assessing the Erosivity of Ice Sheets

In some places, ice sheets erode their bed (Alley et al., 2019; Cowton et al., 2012) and in others, they preserve bed materials (Bierman et al., 2024b). Erosion is considered typical of warm-based ice, whereas bed preservation is thought to be common where ice is cold-based and frozen to the bed (e.g., Corbett et al., 2016b; Sugden, 1978). Data generally support these assumptions (Koppes et al., 2015), although stratigraphic observations clearly show that ice, both warm and cold-based, can leave older sediment deposits intact. For example, field surveys in areas overrun by the LIS reveal weathered regolith (Goldthwait and Kruger, 1938) including the fossiliferous Eocene Brandon Lignite in Vermont (Tiffney, 1994). Much of this weathered material is in valleys where ice was warm-based, meaning that pre-glacial, poorly consolidated material has survived repeated LIS advances in New England, south of the field area investigated here. Weathered saprolite is similarly preserved in southeastern Quebec (LaSalle et al., 1985). On Baffin Island, pre-LGM lake sediments remain in basins overrun by ice as does heavily weathered bedrock on which sit fresh, late-glacial erratic boulders (Bierman et al., 2001; Miller et al., 2022).

Other studies have found varying inherited concentrations of cosmogenic nuclides in boulders and bedrock sampled near and inside the former LGM margin of the LIS. They attributed higher-than-expected nuclide concentrations to boulder reworking and insufficiently deep surface erosion during the brief time the ice occupied areas near the ice margin (Balco et al., 2002; Balter-Kennedy et al., 2024; Barth et al., 2019; Colgan et al., 2002; Halsted et al., 2024). In the northeastern United States, Halsted et al. (2024) suggested that boulders on LIS terminal moraines contained concentrations of inherited ¹⁰Be equivalent to 2–6 ka of surface exposure. Mountain summits in New England, which till deposits indicate were covered by LGM ice, also have samples with much higher than expected cosmogenic exposure ages (Bierman et al., 2015; Corbett et al., 2019; Davis et al., 2015; Koester et al., 2021). In south-central Wisconsin, near the LGM ice margin, three out of five bedrock outcrops sampled had concentrations of ¹⁰Be and ²⁶Al many times higher than predicted based on radiocarbon dating (Colgan et al., 2002). In the Torngat Mountains of northern Labrador, measurements of ²⁶Al and ¹⁰Be on bedrock sites and erratic boulders at mountain summits provide evidence of minimal erosion (< 1.4 m Ma⁻¹) and thus large concentrations of inherited nuclides where cold-based ice was predominant before deglaciation (Staiger et al., 2005).

Significant inheritance of cosmogenic nuclides is observed in areas once occupied by other Northern Hemisphere ice sheets as well. Towards the center of the Scandinavian ice sheet (northeastern Sweden) there is evidence that bedrock outcrops and boulder fields have been preserved through many glacial cycles since the late Cenozoic (Stroeven et al., 2002). A compilation of 953 cosmogenic nuclide measurements made on samples collected within the former extent of the Fennoscandian Ice Sheet, showed that nuclide inheritance was common in both bedrock outcrops and erratic boulders (Jansen et al., 2019).

There is also evidence of minimal erosion near the margin of the Cordilleran Ice Sheet, with 8 out of 23 bedrock samples on Whitbey Island having ³⁶Cl $/$ Cl ratios suggesting inheritance of nuclides produced from prior interglacials (Briner and Swanson, 1998). In northwest Greenland, 8 of 28 sampled boulders had high concentrations of ¹⁰Be and ²⁶Al along with low ²⁶Al $/$ ¹⁰Be ratios indicative of burial, indicating minimal subglacial erosion over multiple interglacial and glacial periods when the ice was presumably cold-based (Corbett et al., 2016b). Along Greenland's western ice margin, most subglacial cobbles (72 out of 86) sampled by Corbett et al. (2021) had a low concentration of ¹⁰Be (median = 1.0 × 10³ atoms g⁻¹), indicative of deep subglacial erosion and/or minimal prior surface exposure. Nevertheless, a subset of samples had higher ¹⁰Be concentrations (> 3 × 10³ atoms g⁻¹, n= 14), suggesting sourcing from minimally eroded (cold-based?) regions of the ice sheet where bedrock and sediment retained some ¹⁰Be produced during prior ice-free periods.

Previously-published cosmogenic nuclide concentrations measured in samples of bedrock, boulders, and sediment from the area covered by the QLID suggest a varied pattern of erosion. Some areas were deeply eroded while others show evidence for inherited ¹⁰Be indicating insufficient subglacial erosion to eliminate near-surface-produced cosmogenic nuclides. For example, of the five boulders sampled by Couette et al. (2023) from the early Holocene Paradise Moraine in eastern Labrador, two have implausible ¹⁰Be exposure ages > 20 ka and the other three implausibly date the Paradise Moraine as older than a margin further from the middle of the ice dome. Ullman et al. (2016) likewise found anomalously high ¹⁰Be concentrations in 10 out of 65 boulders along transects stretching eastward and southward from the middle of the QLID to the coast. Marquette et al. (2004) sampled upland bedrock, felsenmeer, and erratics in the Torngat Mountains and analyzed both ¹⁰Be and ²⁶Al, the concentrations of which suggested that non-erosive ice covered the highlands while erosive ice occupied the valleys, retreating between 11.6 and 13.6 ka. Rice et al. (2019) inferred variable erosion conditions across an ∼ 100 km portion of northeastern Quebec based on factor-of-two variations in ¹⁰Be abundances of bedrock and till samples. Only Carlson et al. (2007), who collected a 13-sample transect of boulders along a 650 km transect westward from Labrador to James Bay, found no significant inheritance.

2.4 Cosmogenic nuclides as sediment source tracers

Cosmogenic nuclides have been used to identify sediment sources for both modern and paleo ice sheets. For example, Nelson et al. (2014) sampled rivers in the deglaciated areas of coastal Greenland. They found ¹⁰Be concentrations in sediment sourced from the ice sheet (6500 ± 4100 atoms g⁻¹) were significantly lower than sediment sourced from deglaciated terrain (14 900 ± 8600 atoms g⁻¹). This difference can be explained by contrasting exposure histories. Outboard of the ice margin, ¹⁰Be concentrations in sediment increased because of exposure to cosmic radiation after Holocene deglaciation while concentrations remained low under the ice sheet where production of ¹⁰Be was minimal.

Nelson et al. (2014) found that sediment sourced from a mix of deglaciated and glaciated terrain had ¹⁰Be concentrations much closer to those of the glacial than deglacial endmember. These results suggest that most sediment moving through river systems in Greenland's paraglacial landscape comes from under the ice sheet rather than the adjacent deglaciated terrain. In southwest Minnesota and eastern South Dakota, Balco (2005) used a similar approach to determine sediment sourcing in a part of the midwestern United States previously covered by the LIS. Rivers there carried sediment with ¹⁰Be and ²⁶Al concentrations (∼ 60 000 and 270 000 atoms g⁻¹, respectively) like that of the glacial sediment they sampled. Both materials had lower ²⁶Al $/$ ¹⁰Be ratios (4.70 ± 0.29) than at production (∼ 6.8).

Table 1 Sample Location and Type.

^* Location and elevation were measured in the field using a Garmin eTrex 20 GPS.

Table 2 Calculations for Holocene Exposure-Corrected Concentrations.

^a Deglaciation ages estimated based on proximity to deglacial isochrons in Dalton et al. (2020). ^b Sample depth estimated in the field. Depth uncertainty estimated from photos and considered 95 % confidence interval. ^cMuon and spallation surface production rates estimated from the online calculator formerly known as CRONUS v2.3 (Balco et al., 2008). ^d Total nuclide production rate at sample depth includes production from both muons and spallation, and assumes an attenuation length of 165 g cm⁻² for spallation, 1400 g cm⁻² for muons, and a sediment density of 1.7 g cm⁻³. We applied no snow shielding correction.

3 Methods

3.1 Field Methods

To constrain nuclide concentrations in materials deposited by the LIS, we sampled deglacial sediment deposits (n= 10) including deltas and eskers (Table 1). We collected modern river sediment samples (n= 11) from trunk streams (St. Lawrence and Churchill River) as well as smaller tributaries. We collected deglacial sediment from clean faces in gravel pits or along river bluffs from 2 to 30 m below the upper land surface to reduce the effect of Holocene nuclide production following deglaciation (Table 2). We used shovels to dig ∼ 0.3 m into the side of the landform before collecting ∼ 500 g of sand. Modern river sediment samples were collected along channels away from areas of recent disturbance.

We collected outcrop samples using a rock saw and chisels at two sites, one interior and one coastal. The interior sample is from an exposed outcrop on a glacially molded bedrock knob in the middle of our eastern transect. This sample is close to Ullman et al. (2016) site CL3, in which two of five boulders they sampled contained what the authors interpreted as inherited nuclides from period(s) of prior near-surface exposure. The coastal location is a roadcut 4 km inboard of the Belles Amours Moraine at the northern tip of the Gulf of St. Lawrence, where we obtained 7 samples along a 1.9 m bedrock depth profile; the surface sample (LBD19-46) was reported previously by Couette et al. (2023).

3.2 Laboratory Methods

3.2.1 University of Vermont

To isolate and purify quartz for cosmogenic nuclide analysis at the University of Vermont (all sediment samples and the surface bedrock sample) we used a series of physical and chemical processes (Kohl and Nishiizumi, 1992) after sieving material to between 250 and 850 µm. We performed two 24 h 6 N hydrochloric acid (HCl) etches in heated ultrasonic baths to remove grain coatings. We then used dilute (1 %) hydrofluoric (HF) and nitric acid (HNO₃) etches for three 24-hour periods after which we sonicated samples in 0.5 % HF and HNO₃ for a minimum of two weeks. We evaluated the purity of etched samples using inductively coupled plasma optical emission spectrometry (ICP-OES) after which impure samples were re-etched until each was sufficiently pure.

Table 3 Blank-Corrected Isotopic Data.

Nuclide concentrations and uncertainties incorporate blank corrections. Uncertainties are 1σ. Depths of samples LBD19-46 to LBD19-52 in bedrock depth profile are given after sample name.

We extracted Be and Al from the purified quartz samples (17.3–22.2 g, n= 22) in the National Science Foundation/University of Vermont Community Cosmogenic Facility using methods described in Corbett et al. (2016a). Samples were prepared in two separate batches, each of which included a fully processed blank. We spiked the samples with ∼ 250 µg ⁹Be using a beryl carrier made in the Community Cosmogenic Facility with a Be concentration of 348 µg mL⁻¹ (Table 3). We spiked samples with SPEX ICP Al standard (1000 ppm) as needed based on their quantity of native Al, ensuring at least 1500 µg of total Al in every sample. We quantified total Al in the samples by ICP-OES immediately following sample digestion. Following standard procedures (Corbett et al., 2016a), we removed replicate aliquots from the samples by mass (representing ∼ 2 % and ∼ 4 % of the sample, respectively), added 25 µL H₂SO₄ to each, evaporated the HF, then diluted the residual H₂SO₄ droplets by mass with a 0.25 % H₂SO₄ solution spiked with Y as an internal standard.

Purdue Rare Isotope Measurement Laboratory performed accelerator mass spectrometry analysis (AMS). For ¹⁰Be $/$ ⁹Be, measured ratios were normalized to primary standard 07KNSTD3110 with a ratio of 2.850 × 10⁻¹² (Nishiizumi et al., 2007). For ²⁶Al $/$ ²⁷Al, analyses were normalized to primary standard KNSTD with a ratio of 1.818 × 10⁻¹² (Nishiizumi, 2004).

We used the known concentration of ⁹Be added as carrier, along with the measured isotopic ratio and quartz mass, to calculate the concentration of ¹⁰Be in each sample. Because of the native ²⁷Al within the samples, the concentration of ²⁷Al measured using ICP-OES after quartz dissolution was used to calculate the concentration of ²⁶Al. We subtracted the mean extraction process blank ratios of ¹⁰Be $/$ ⁹Be ((7.41 ± 2.81) × 10⁻¹⁶; n= 2) and ²⁶Al $/$ ²⁷Al ((5.39 ± 0.71) × 10⁻¹⁶; n= 2) from the measured ratios and propagated the uncertainty in quadrature (Table 3).

3.2.2 University of Strasbourg

To isolate and purify quartz for cosmogenic nuclide analysis at the University of Strasbourg, Institut Terre et Environnement de Strasbourg (7 bedrock profile samples, LBD19-46 to 52) we used a series of physical and chemical processes (Bierman et al., 2002; Kohl and Nishiizumi, 1992). We sieved material to between 250 and 1000 µm. For each sample, we performed a 12 h 6 N HCl leaching to eliminate oxides and organic material. Magnetic minerals were then separated with a Frantz magnetic separator. We used dilute (1 %) HF and HNO₃ acid etches in ultrasonic baths for five 12 h periods. We evaluated the purity of etched samples using inductively coupled plasma optical emission spectroscopy (ICP-OES) after which impure samples were re-etched until each was sufficiently pure.

We extracted Be and Al from the purified quartz samples (27.02–56.99 g, n= 7) in the cosmogenic nuclide laboratory at the University of Strasbourg following a protocol modified from Kohl and Nishiizumi (1992). Samples were prepared together with a fully processed blank. We spiked the samples with ∼ 250 or ∼ 500 µg ⁹Be using a Scharlab ICP 1000 mg L⁻¹ Be standard. We spiked only the blank with a Scharlab ICP 1000 mg L⁻¹ Al standard as the quantity of native Al was sufficient in every sample. We quantified total Al in the samples by ICP-OES following sample digestion and three successive fuming steps with 2 to 4 mL of perchloric acid (HClO₄). Be and Al were isolated through anion and cation exchange columns and precipitated as hydroxides. They were then dried and calcinated to BeO and AlO at 750 °C.

AMS analysis was performed at ASTER, the French National Laboratory of Cosmogenic Nuclides in CEREGE (Aix-en-Provence). For ¹⁰Be $/$ ⁹Be, measured ratios were normalized to in-house CEREGE standard STD-11 with an assumed ratio of (1.191 ± 0.013) × 10⁻¹¹ (Braucher et al., 2015). For ²⁶Al $/$ ²⁷Al, analyses were normalized to in-house CEREGE standard SM-AL-11 with a ratio of (7.401 ± 0.064) × 10⁻¹² (Merchel and Bremser, 2004).

We used the known concentration of ⁹Be added as carrier, along with the measured isotopic ratio and quartz mass to calculate the concentration of ¹⁰Be in each sample. Because of the native ²⁷Al within the samples, the concentration of ²⁷Al measured using ICP-OES after quartz dissolution was used to calculate the concentration of ²⁶Al. We subtracted the extraction process blank ratios of ¹⁰Be $/$ ⁹Be ((5.21 ± 0.42) × 10⁻¹⁵) and ²⁶Al $/$ ²⁷Al ((6.25 ± 6.25) × 10⁻¹⁶) from the measured ratios and propagated the uncertainty in quadrature.

3.3 Statistical Analysis

We use the same value (α= 0.05) for all statistical tests we performed. We used both Wilcoxon rank-sum tests (due to the non-normal distribution of nuclide concentrations in modern sediment samples) and Tukey HSD tests to investigate significant differences between sample groups. We used a 2-standard deviation threshold for detectability; that is, if twice the analytical uncertainty exceeded the measured ratio, then we considered the sample to be below detection limits. This provides a 95 % confidence that the isotopic ratios and concentrations we report are finite.

3.4 Data Compilation

To expand our assessment of cosmogenic nuclide concentrations in materials influenced by the QLID, we estimated the concentration of inherited nuclides present in published ¹⁰Be data generated from samples of boulders and bedrock collected in the region by others. We obtained data from 200 samples (42 bedrock, 158 boulders) from Dalton et al.'s (2023) North American ¹⁰Be compilation and recalculated all ages with the online exposure calculator formerly known as CRONUS (Balco et al., 2008) (constant production rate model, version 3.0.2, constants 26 August 2024) using the global production rate and LSDn scaling (Lifton et al., 2014), accounting for differences in standardizations. We also included 37 ¹⁰Be ages (3 bedrock, 34 boulders) from Couette et al. (2023) that were not in the Dalton et al. (2023) compilation and our one surface bedrock sample, following the same protocol.

3.5 Holocene Exposure Correction

For deglacial sediment samples and the single inland bedrock sample, we calculated the concentration of nuclides attributable to Holocene exposure. We used the online exposure age calculator (constant global production rate model, version 2.3) to determine the surface production rate (atoms g⁻¹ yr⁻¹) of ¹⁰Be and ²⁶Al for both muons (P_μ) and spallation (P_s) at each sample site (Balco et al., 2008). The production rate at depth was then estimated using an attenuation length of 165 g cm⁻² for spallation (Λ_s) and 1400 g cm⁻² for muons (Λ_μ). We validated the apparent muonogenic attenuation factor using the CRONUS implementation of Heisinger et al. (2002) (code P_mu_total.m) for sea level–yielding < 1 % difference from our original calculations. We estimated uncertainties in the concentration of nuclides produced during the Holocene due to uncertainties in our sample depths and combined these in quadrature with AMS measurement uncertainties (Table 2). For sediment, we assumed a density (ρ) of 1.7 g cm⁻³.

After establishing the nuclide production rate at depth, we calculated (using Eq. 1) the production since deglaciation (A in yr) of ¹⁰Be and ²⁶Al (H in atoms g⁻¹) at the depth we collected each sample (D in cm), taking deglaciation age estimates for each sample site from Dalton et al. (2020) using radiocarbon-based LIS margin isochrons. For the rock sample, we used half the sample thickness (thickness = 2 cm) as the depth. Dalton 2020 gives isochrons in our study area for 6.3, 6.8, 7.3, 7.9, 8.1, 8.5, 8.7, 8.8, 9, 9.6, 10.3, 10.9, 11, 11.5, 11.8, 12.1, 12.8, 13.5, 14.2, and 14.9 cal ka. Thus, the resolution of the isochrons ranges from 100 to 700 years and averages 500 years. For samples located between isochrons, we estimated their deglaciation ages by linearly interpolating between the isochrons. Thus, typical uncertainties on our estimated deglaciation ages are likely a few centuries.

$$\begin{array}{}\text{(1)}& H=P_{\mathrm{s}}\cdot \mathrm{1}/\exp \left({\displaystyle \frac{D\cdot \mathit{\rho }}{{\mathrm{\Lambda }}_{\mathrm{s}}}}\right)\cdot A+P_{\mathit{\mu }}\cdot \mathrm{1}/\exp \left({\displaystyle \frac{D\cdot \mathit{\rho }}{{\mathrm{\Lambda }}_{\mathit{\mu }}}}\right)\cdot A\end{array}$$

We subtracted the concentration of nuclides produced after deglaciation from the sample's measured ¹⁰Be concentration to calculate the concentration of nuclides inherited at the time of deglaciation. We did not use the more recent reconstruction of LIS retreat from Dalton et al. (2023) to determine deglaciation ages because it incorporates updates based on the ¹⁰Be data and thus risks circularity. We estimated the amount of inheritance in our coastal bedrock depth profile by fitting a production profile to the data representing 11.5 kyr of exposure (the age of the adjacent Belles Amours Moraine; Couette et al., 2023) with varying concentrations of uniformly distributed inherited nuclides presumably from deep muonogenic production during past interglacials.

4 Results

¹⁰Be and ²⁶Al were above detection limit in all 8 bedrock samples and 21 of 22 sediment samples we analyzed, the one exception being the ²⁶Al measurement for sample GB-04, a modern stream sample (Table 3, Fig. 3). For both nuclides, the average concentration is higher in the modern sediment (3.31 ± 1.57 × 10^4 10Be atoms g⁻¹, 2.12 ± 1.18 × 10^{5 26}Al atoms g⁻¹) than in the deglacial sediment (2.25 ± 1.30 × 10^{4 10}Be atoms g⁻¹, 1.47 ± 0.94 × 10^{5 26}Al atoms g⁻¹) but these differences are not statistically significant.

Figure 3 Maps showing new isotopic data excluding bedrock depth profile. (A) ²⁶Al $/$ ¹⁰Be ratio (n= 21). (B) Measured concentration of ¹⁰Be (n= 22). (C) Measured concentration of ²⁶Al (n= 21). Data in Table 3. Sample names and key as in Fig. 1.

The bedrock sample (GB-06) had the highest concentration of ¹⁰Be (7.33 ± 0.39 × 10⁴ atoms g⁻¹) and ²⁶Al (5.91 ± 0.29 × 10⁵ atoms g⁻¹) that we measured. The coastal bedrock depth profile had a surface concentration of 5.02 ± 0.19 × 10^{4 10}Be atoms g⁻¹ and 3.30 ± 0.14 × 10^{5 26}Al atoms g⁻¹ (Fig. 4). Concentrations of both nuclides decrease exponentially with depth in the bedrock profile, with e-folding lengths of approximately 175 to 190 g cm⁻², longer than would be expected from spallation production alone (165 g cm⁻²). We note that the sample at depth 10 cm does not fit well on the profile. It was taken to the left side of the other samples (Fig. 2E) and may have had a different cover history.

Figure 4 Depth profile in bedrock (sample set LBD-19, 46-52) with modelled fit for 11.5 kyr of surface exposure and inheritance of 0 and 2000 atoms g⁻¹. Symbols extend to 1 standard deviation analytic uncertainty. Solid line is fit with no inherited nuclides. Dashed line assumes uniform 2000 atoms g⁻¹ inheritance for ¹⁰Be. ²⁶Al was scaled to equivalent ¹⁰Be concentration assuming ²⁶Al $/$ ¹⁰Be ratio of 7.3. Shaded sample collected to the left of the profile (Fig. 2E) and not considered in fit.

4.1 ¹⁰Be and ²⁶Al Concentrations in Deglacial Samples Corrected for Holocene Exposure

After correcting for Holocene exposure, there is significant difference between the concentration of ¹⁰Be in deglacial (mean, 1SD, n= 10; 1.95 ± 1.38 × 10⁴ atoms g⁻¹) and modern (n= 11, 3.31 ± 1.57 × 10⁴ atoms g⁻¹) sediment samples (p= 0.02) (Table 4, Fig. 5A). For ²⁶Al, the concentration of deglacial (1.25 ± 1.02 × 10⁵ atoms g⁻¹) and modern (2.12 ± 1.18 × 10⁵ atoms g⁻¹) samples are also significantly different (p= 0.04) (Fig. 5B). The modern samples are more variable in ¹⁰Be concentration, with an interquartile range (IQR) of 2.76 × 10⁴ atoms g⁻¹ compared to the IQR of 7.85 × 10³ atoms g⁻¹ for Holocene exposure-corrected deglacial samples.

Figure 5 Comparison of modern and Holocene-exposure-corrected deglacial samples. (A) ¹⁰Be concentrations, (B) ²⁶Al concentrations, and (C) ²⁶Al $/$ ¹⁰Be ratios for deglacial (Holocene-corrected) and modern samples. The dashed line in panel (C) indicates the empirically determined production ratio (7.3) at high latitudes from Corbett et al. (2017). Points represent individual samples with 1σ propagated errors for modern (yellow) and deglacial (red) samples. Box and whisker plots are shown for each dataset with whiskers going to the smallest and highest values. The box extends from the 25th to the 75th percentiles. The line in the middle of each box is the median.

Both the bedrock sample (GB-06) and the bedrock depth profile (LBD19-46-52) contain nuclides inherited from prior exposure. Inferring an exposure age of 8.1 kyr from maps of Dalton et al. (2020), and correcting for Holocene exposure with no erosion (the outcrop is smooth and rounded, Fig. 2C), GB-06 contains 1.79 × 10⁴ atoms g⁻¹ of inherited ¹⁰Be and 2.15 × 10⁵ atoms g⁻¹ of inherited ²⁶Al. Nuclide profiles, modelled to reflect Holocene exposure, best fit the bedrock profile data when they include ∼ 2 × 10³ atoms g⁻¹ of inherited ¹⁰Be and ∼ 1.4 × 10⁴ atoms g⁻¹ of inherited ²⁶Al (Fig. 4).

Our analysis of previously published ¹⁰Be data suggests that inheritance of cosmogenic nuclides from prior period(s) of exposure is common in Quebec and Labrador. Over half of bedrock samples (n= 26 of 46, including our new sample GB-06) and a third of boulder samples (n= 65 of 192) contain significant concentrations of inherited nuclides, which we consider > 3 × 10³ atoms g⁻¹ because that is the median analytical uncertainty in the dataset. The mean and range of inherited nuclide concentrations are larger for samples of bedrock (mean = 1.21 × 10⁵ atoms g⁻¹; IQR = 8.04 × 10⁴ atoms g⁻¹) than boulders (mean = 1.40 × 10⁴ atoms g⁻¹; IQR = 1.28 × 10⁴ atoms g⁻¹).

4.2 ²⁶Al $/$ ¹⁰Be Ratios

The error-weighted mean of ²⁶Al $/$ ¹⁰Be ratios for deglacial sand (correcting for Holocene nuclide production) and modern sediment samples (6.1 ± 0.3 and 6.6 ± 0.5, respectively, one standard deviation) are lower than the measured Greenland production ratio of 7.3 ± 0.3 but more similar to that at mid-latitudes of ∼ 6.8 (Fig. 5C, p= 0.01, Corbett et al., 2017). The ²⁶Al $/$ ¹⁰Be ratios for deglacial samples are more variable (IQR = 2.24) than ratios measured in modern stream samples (IQR = 0.78, Fig. 5C).

The ²⁶Al $/$ ¹⁰Be of deglacial samples exhibit a significant, positive linear trend with increasing distance from the middle of the QLID (r² = 0.45, p= 0.03) (Fig. 6). The five samples closest to the middle of the ice dome (5.3 ± 0.4) have lower error-weighted average ²⁶Al $/$ ¹⁰Be ratios than five samples farther away (7.0 ± 0.4, Fig. 6). Modern samples, in contrast, exhibit no spatial trend in ²⁶Al $/$ ¹⁰Be ratios. Tukey HSD tests show a significant difference between ²⁶Al $/$ ¹⁰Be ratios in deglacial samples closest to and further from the middle of the ice dome (both tests: p < 0.01).

Table 4 Holocene Exposure-Corrected Isotopic Data for Deglacial Samples.

^* Uncertainties account for both analytical and sample depth uncertainties. The depth uncertainties in combination with nonlinear production rate changes with depth yield asymmetrical nuclide uncertainties except for the surface bedrock sample, GB-06. Uncertainties are 1σ.

Figure 6 ²⁶Al $/$ ¹⁰Be ratios versus distance from middle of the QLID, which we have designated as Labrador City. Deglacial data are corrected for Holocene nuclide production and fit with a trendline. Error bars are 1σ. Dashed line shows production ratio at high latitudes (7.3, Corbett et al., 2017).

5 Discussion

Measurements of in situ-produced ¹⁰Be and ²⁶Al, presented here and in previously published papers, indicate that many samples of bedrock (57 %), erratic boulders (34 %), and deglacial sediment (100 %), collected from terrain once covered by the QLID of the LIS, contain cosmogenic nuclides accumulated during periods of exposure prior to the last deglaciation. Our findings for bedrock and boulders are remarkably similar to those of a compilation of nearly 1000 samples collected from areas overrun by the Fennoscandian Ice Sheet and analyzed for cosmogenic nuclides (Jansen et al., 2019). Those authors found that 27 % of erratic boulders contain inherited nuclides along with 64 % of sampled bedrock outcrops (Jansen et al., 2019).

These data, along with many other data sets, show convincingly that in areas once covered by temperate, mid-latitude ice sheets, nuclide inheritance from periods of prior exposure is geographically widespread, commonplace, and spatially variable in both detrital and bedrock samples. Our analysis of glacial sediment presented in this paper provides a different way of examining the magnitude of inheritance because it averages nuclide concentrations retained in earth materials from periods of exposure prior to the LGM on a landscape-scale. Such averaging avoids potential spatial bias of point measurements on outcrops and erratic boulders and provides a more generalized view of glacial erosion, sediment residence and transport times, and cosmogenic nuclide inheritance.

5.1 Differences in ¹⁰Be inheritance between sample types

Compilation of extant cosmogenic data shows that samples collected from erratic boulders contain lower average concentrations of inherited ¹⁰Be (14 000 atoms g⁻¹) than those collected from bedrock outcrops (120 700 atoms g⁻¹) and deglacial sediment (19 500 atoms g⁻¹, Figs. 7 and 8). Of the 192 reported boulder analyses, we calculate that 65 boulder samples (34 % of total) contain ¹⁰Be inherited from pre-LGM periods of exposure on the basis of having > 3000 atoms g⁻¹ of ¹⁰Be (the median analytical uncertainty) above that produced during deglacial exposure as determined from Dalton et al.'s (2020) isochrons. The average boulder inheritance is equal to ∼ 1800 years of surface exposure, a value skewed high by relatively few boulders with large concentrations of inherited nuclides. In contrast, the median boulder age is about 650 years younger than Dalton et al.'s (2020) radiocarbon-based chronology. Such a discrepancy between median boulder exposure ages and radiocarbon ages of organic material likely reflects a combination of boulder surface erosion, boulder burial by snow or till (Schildgen et al., 2005), ¹⁰Be production rate uncertainty (Balco et al., 2009), the effect of isostatic uplift on nuclide production rates, and the incorporation of old carbon in radiocarbon dated samples (see references in Halsted et al., 2024).

Figure 7 Spatial distribution of inheritance interpreted from published data corrected using deglacial age model of Dalton et al. (2020). (A) Boulder data, atoms g⁻¹ inherited ¹⁰Be. (B) Boulder data, inherited years equivalent surface exposure. (C) Bedrock data, atoms g⁻¹ inherited ¹⁰Be. (D) Bedrock data, inherited years equivalent surface exposure. Red circles are deglacial sediment data (this paper). Inheritance scaled to circle diameter.

We find that 26 of 46 sampled bedrock outcrops (57 %) contain detectable ¹⁰Be (mean = 120 700 atoms g⁻¹; median = 21 400 atoms g⁻¹) inherited from pre-LGM periods of exposure. Considered as exposure time, this is an average inheritance for bedrock samples of 14 600 years. Because the distribution is skewed, the median inheritance, when considered as a surface exposure time, is less, 3200 years (Fig. 7). Lower mean and median concentrations of inherited nuclides in erratic boulders than in exposed bedrock from beneath the former QLID suggest that the boulders were quarried from below the present-day bedrock surface, are preferentially derived from more erosive areas of the ice sheet, and/or were reworked from deposits deep enough that interglacial exposure to cosmic radiation was insufficient to create high concentrations of ¹⁰Be.

All 10 deglacial sand samples (100 %) contain inherited nuclides (median = 18 400 atoms g⁻¹; mean = 19 500 atoms g⁻¹), the equivalent of 2700 and 3500 years of surface exposure (Fig. 7). Unlike the distribution of boulder and bedrock sample ¹⁰Be concentrations, that of deglacial sand samples is not skewed as indicated by the similarity of the mean and median. This similarity likely reflects the physical mixing of thousands of sand grains included in the sediment samples, each with its own idiosyncratic burial and exposure history. We recognize that there is likely a sampling bias in the location of bedrock and boulder samples – at least in part dictated by limited road access. Sediment samples, because they integrate across much larger areas, are less subject to such a bias.

We conducted a sensitivity analysis, which shows that our finding of nuclide inheritance is robust to uncertainty in the time/retreat distance model that we used to correct for post-depositional nuclide production. For surface boulder and bedrock samples, if Dalton et al.'s (2020) retreat chronology were 10 % too old, then the number of such samples with inheritance above analytical uncertainty (3000 atoms g⁻¹) would increase from 90 to 120. Conversely, if the Dalton et al. (2020) model reported deglacial ages that were 10 % too young, then the number of samples with inheritance would decrease from 90 to 75. For buried deglacial sand samples, the influence of the age model is much less than for boulders and bedrock collected at the surface because the cosmic-flux is strongly attenuated with depth.

Because isostatic uplift in response to deglaciation has raised sample sites to higher elevations than they were at the time of deglaciation, nuclide production rates at the sampling locations have increased over time (Jones et al., 2019). Using the uplift-corrected data provided in Dalton et al. (2023), we calculated that accounting for such post-glacial uplift would increase the concentration of nuclides inherited at initial exposure by 3400 ± 2400 atoms g⁻¹. While such a correction is uncertain, it suggests that our estimates of nuclides inherited from periods of prior exposure are likely conservative.

5.2 Spatial variability in ice sheet erosion depth

Our compilation of data from 248 samples (deglacial sediment and surface bedrock samples, this paper, and prior published bedrock and boulder data, Supplement) analyzed for cosmogenic nuclides indicates that the QLID of the LIS did not deeply or uniformly erode its bed, nor did it effectively remove cosmogenic nuclides produced in sediment and rock exposed during prior interglacial periods (Figs. 7 and 8). The broad distribution of inherited ¹⁰Be in bedrock is consistent with spatially variable ice sheet erosivity and/or varying thickness of interglacial sediment cover shielding the underlying rock from cosmic rays when the ice was gone and has been noted by others, particularly at high elevations in the eastern part of the study area (Marquette et al., 2004; Staiger et al., 2005).

Figure 8 Histograms (logarithmic) of inherited ¹⁰Be in deglacial sediment (this study), glaciated bedrock (this study and others), and erratic boulders (other studies) sampled from eastern Canadian landscapes once covered by the QLID of the LIS. (A) Inheritance (atoms g⁻¹). (B) Inheritance expressed as years of surface exposure calculated using site-specific ¹⁰Be production rates to remove elevation dependence of nuclide production rates. Sample count, average, and median values indicated with arrows. Values < 3000 atoms g⁻¹, the median uncertainty in measurements, are considered to have no detectable inheritance and are counted as zero for statistical analysis.

Subglacial process modeling over eastern North America (Melanson et al., 2013) suggests a variably erosive LIS in Quebec and Labrador during the last glacial cycle consistent with our findings. Melanson et al. (2013) suggest that total erosion per glacial cycle, predicted using empirical abrasion and quarrying laws, ranges from near zero under the middle of the Quebec-Labrador dome to ≥ 10 m along parts of its Atlantic and St. Lawrence margin. Melanson et al.'s (2013) numerical modeling of the Quebec-Labrador region exhibits similar minima for both basal sliding speed and total ice movement integrated over the last glacial cycle – both variables related to the efficacy of glacial erosion. Simulated ice sliding distances (the integrated basal velocity over the last glacial cycle in millions of meters, Mm) are near zero in the middle of our study area, 1 Mm near Goose Bay, and 3 Mm along the St. Lawrence estuary – an order of magnitude less than for the Hudson Strait ice stream and southern LIS lobes (Melanson et al., 2013). Landscape analysis using remotely sensed data also suggests variable bed conditions and erosion efficacy of the LIS QLID over time and space (Rice et al., 2020).

The average inherited nuclide concentration we measured in LIS deglacial sediment (19 500 ¹⁰Be atoms g⁻¹) is several times higher than the average found in sediment shed by the Greenland Ice Sheet today (∼ 6500 ¹⁰Be atoms g⁻¹, Nelson et al., 2014), but much lower than in sediment deposited by the LIS in the midwestern United States (∼ 60 000 ¹⁰Be atoms g⁻¹, Balco, 2005). The low concentrations of ¹⁰Be in sediment issuing from the Greenland Ice Sheet (as well as low ²⁶Al $/$ ¹⁰Be ratios ∼ 4.5, Bierman et al., 2016) reflect continuous ice cover through many Pleistocene interglacials and erosive warm-based ice in areas from which sediment is sourced. In contrast, the high concentrations of ¹⁰Be in deglacial sediment originating from the southern margin of the Keewatin Ice Dome in Minnesota (Balco, 2005) suggest that ice there was frequently absent and when present, only weakly erosive. Some data suggest that in central North America the LIS was at least in part, cold-based near the margin (cf., Colgan et al., 2002).

5.3 Implications for sample collection and surface exposure dating

Both our compilation (LIS) and that for the Fennoscandian Ice Sheet of Jansen et al. (2019) differ from that of Heyman et al. (2011) who, 15 years prior to our work, suggested that “few boulders from the paleo-ice sheets have exposure ages significantly older than independently known deglaciation ages, indicating that prior exposure is of limited significance.” This difference likely reflects numerous additional cosmogenic nuclide measurements that have been made since the publication of Heyman et al. (2011), a time over which AMS measurement precision has greatly improved and the number of other pertinent dates using different chronometers to constrain the time of deglaciation has increased. As Heyman et al. (2011) clearly point out, such inheritance, if not detected and corrected for, has the potential to increase cosmogenic exposure ages, biasing them high, including those from exposed rock and glacial deposits.

The new and compiled data have implications for sample collection and interpretation in areas once covered by now-vanished ice sheets, specifically the QLID of the LIS. The compiled data show clearly that sampling boulders rather than bedrock is most likely to minimize the influence on calculated landform age of nuclides inherited from exposure prior to the last deglaciation. Because the distributions of inherited ¹⁰Be concentrations in both the bedrock and boulder populations are skewed (Fig. 7), extreme outliers are potentially detectable. Indeed, Ullman et al. (2016) excluded 10 of 65 boulder samples because of their unusually high concentration of ¹⁰Be. Couette et al. (2023) similarly excluded 5 of 37 samples in eastern Labrador because of high ¹⁰Be concentrations.

Lower concentrations of inherited nuclides, likely the result of deeply penetrating muons during interglacial times without ice cover, are more difficult to detect and correct for. This inheritance could skew slightly the average age of an exposure age data set. In that regard, the depth profile data (samples LBD19-46-52) are illustrative because they are consistent with several thousand atoms g⁻¹ of inherited nuclides at depth of several meters (Fig. 4). Such inheritance would increase exposure ages of plucked boulders and glacially eroded outcrops by hundreds of years (a phenomenon noted by Ullman et al. (2016) who reported that ¹⁰Be-based deglacial chronology was older than radiocarbon-based estimates of local deglaciation by centuries, their Fig. 7).

While radiocarbon age lags have been attributed to delayed colonization of the landscape by vegetation following deglaciation (e.g., Halsted et al., 2024; Peteet et al., 2012), the nuclide inheritance we find in Quebec and Labrador suggests that in this region, many ¹⁰Be ages from boulder and bedrock samples will be too old by at least hundreds of years. The overall agreement of the Ullman et al. (2016) cosmogenic ages with those of Dalton et al. (2022) and Dyke (2004) suggests that in our study area a population of boulder dates, screened for outliers, is likely to be representative of deglaciation age at ky resolution.

5.4 Interglacial presence of QLID and ²⁶Al $/$ ¹⁰Be ratio of IRD

Finding ²⁶Al $/$ ¹⁰Be ratios below the nominal Arctic surface production ratio of 7.3 (Corbett et al., 2017) in quartz IRD from North Atlantic Heinrich layers, LeBlanc et al. (2023) suggested that ice sheet remnants may have lingered in parts of eastern Canada for the majority of Pleistocene interglacials burying the land surface and allowing for the ²⁶Al $/$ ¹⁰Be ratio to lower by radioactive decay (Fig. 9). While Heinrich layer sediment was predominantly delivered to the ocean by the Hudson Strait ice stream (Bond et al., 1992), the quartz IRD that LeBlanc et al. (2023) analyzed most likely came from interior areas of the LIS (where crystalline rocks crop out) as opposed to Hudson Strait, which is primarily underlain by carbonate rocks. Burial by till or in outwash deposits could have the same effect – but with the paucity of glacial sediment over most of the QLID landscape (Pelletier et al., 2016), burial of a sufficient volume of sediment at depths of many meters seems implausible.

Figure 9 Summed probability plots of ²⁶Al $/$ ¹⁰Be ratios for Arctic samples. (A) Samples with simple exposure histories in Greenland (Corbett et al., 2017). (B) Modern stream sediment (this study). (C) Deglacial samples corrected for Holocene exposure; black 5 samples most proximal to Labrador City (approximate middle of the dome); grey 5 samples are most distal (this study). (D) IRD samples from North Atlantic (LeBlanc et al., 2023). Analytical error-weighted mean and 1 SD uncertainty above each plot.

Depressed ²⁶Al $/$ ¹⁰Be ratios in the Holocene exposure-corrected deglacial sediment that we sampled are consistent with ice lingering near Labrador City (the middle of the ice dome) for at least some interglacial periods, shielding rock and sediment below from cosmic rays (Fig. 6). In contrast, the higher ²⁶Al $/$ ¹⁰Be ratios in deglacial sediment toward the edge of the ice dome are consistent with repeated, extended exposure during interglacials. Thus, it appears that the quartz analyzed by LeBlanc et al. (2023) is unlikely to have been sourced from much of the area we sampled in southeastern Canada.

It is possible that incorporation of the IRD sampled by LeBlanc et al. (2023) occurred primarily toward the middle of the QLID and not closer to the margins. It is also possible that the IRD sampled by LeBlanc et al. (2023) was sourced from a wider area (the Keewatin Dome and/or Baffin Island) where samples of sediment (Balco, 2005) and those taken from outcrops (Corbett et al., 2016c; Marsella et al., 2000) have sufficiently low ²⁶Al $/$ ¹⁰Be ratios to match those measured by LeBlanc et al. (2023). Low ratios would also result if IRD were stored in Hudson Bay for ∼ 1 Ma before transport to the deep sea, during which time ²⁶Al and ¹⁰Be would decay and the ²⁶Al $/$ ¹⁰Be ratio would decrease. More extensive sampling of eastern Canada, including Quebec-Labrador, and farther north near Hudson Bay and Baffin Island, would provide further evidence on how persistent different sectors of the eastern LIS were during Pleistocene interglacials.

5.5 Sediment Sourcing in Modern Rivers

We estimate the percentage of sediment derived from erosion of deglacial materials using a two-component, linear-mixing model based on the measured average concentrations of ¹⁰Be in both river and deglacial sediment and assumptions about nuclide production since deglaciation. One component is deglacial deposits which, based on our sampling, today contain an average of 22.5 × 10^{3 10}Be atoms g⁻¹ and enter rivers by bank incision (Fig. 2A). The second component is surficial materials which, when eroded, enter the drainage network.

Since we did not sample these surficial materials directly, we calculate their ¹⁰Be concentration by assuming that, at the time of deglaciation (the beginning of exposure), surface sediment contained the average Holocene exposure-corrected concentration of ¹⁰Be for deglacial sediment (19.5 × 10^{3 10}Be atoms g⁻¹). We then use the average deglacial age of 9.32 ka for our field area and the average ¹⁰Be surface production rate (7.82 atoms yr⁻¹ g⁻¹) to calculate that surface material gained about 7.3 × 10^{4 10}Be atoms g⁻¹ since deglaciation. The surface-exposed end member therefore contains 92.3 × 10^{3 10}Be atoms g⁻¹.

Knowing that modern sediment contains on average 33.1 × 10^{3 10}Be atoms g⁻¹ (Table 3), the two-component mixing model suggests that about 85 % of sediment in eastern Quebec and Labrador rivers today is derived from incision of glacial deposits and the remaining about 15 % comes from erosion of surficial sediment and bedrock. Our findings for eastern Quebec and Labrador are similar to those of Balco (2005) in Minnesota. In that previously glaciated region, they suggest most sediment carried by contemporary rivers comes from incision of glacial deposits. Our finding is also consistent with the low sediment yield of forested upland terrains in other glaciated areas of the LIS (e.g., Dethier et al., 2018).

6 Conclusions

Analysis of cosmogenic ¹⁰Be and ²⁶Al in deglacial (n= 10) sediment samples, in a bedrock depth profile (n= 7), in our new as well as previously-published boulders and bedrock sample data (n= 238), and in one bedrock sample we collected, indicates that nuclides produced during prior interglacial exposures are commonplace in samples collected from the landscape beneath the now-vanished QLID of the LIS. These isotopic data indicate erosion depths in bedrock are no more than a few meters per glacial cycle and/or that regolith generated and exposed during previous glacial and interglacial cycles was not efficiently removed from the terrestrial landscape during the last glaciation. Such inheritance explains numerous outliers in previously published data sets and suggests pervasive inheritance of low concentrations of in situ produced cosmogenic nuclides. Boulders in the area covered by the QLID appear to carry, on average, lower concentrations of inherited nuclides than bedrock outcrops suggesting that dating studies with limited resources may benefit from sampling boulders rather than bedrock. On average, nuclide concentrations in modern fluvial sediment are only slightly higher than those in deglacial sediment, implying most sand, transported by rivers today in landscapes once covered by the QLID, is sourced from eroding banks composed of glacial deposits rather than surrounding land surfaces.