Arctic drainage of Laurentide Ice Sheet meltwater throughout the past 14,700 years

Abstract

During the last deglaciation substantial volumes of meltwater from the decaying Laurentide Ice Sheet were supplied to the Arctic, Gulf of Mexico and North Atlantic along different drainage routes, sometimes as catastrophic flood events. These events are suggested to have impacted global climate, for example initiating the Younger Dryas cold period. Here we analyze the authigenic Pb isotopic composition of sediments in front of the Arctic Mackenzie Delta, a sensitive tracer for elevated freshwater runoff of the retreating Laurentide Ice Sheet. Our data reveal continuous meltwater supply to the Arctic along the Mackenzie River since the onset of the Bølling–Allerød. The strongest Lake Agassiz outflow event is observed at the end of the Bølling–Allerød close to the onset of the Younger Dryas. In context of deglacial North American runoff records from the southern and eastern outlets, our findings provide a detailed reconstruction of the deglacial drainage chronology of the disintegrating Laurentide Ice Sheet.

Introduction

Continuing research efforts over the last decades allowed the identification of major drainage routes for meltwater from the decaying Laurentide Ice Sheet (LIS) (South: Mississippi River; East: St. Lawrence River; North: Hudson Bay/Strait¸ Northwest: Mackenzie River) (Literature review in Teller (2013)¹ and Fisher (2020)²). While consensus exists that the drainage to the south via the Mississippi River was the major meltwater drainage route for several millennia after the Last Glacial Maximum (LGM)^3,4, there has been active debate on the dispersion of drainage routes following the onset of the warm Bølling–Allerød^{1,2,5,6,7,8,9,10,11,12,13,14,15}. Proglacial lakes formed and increased in size with continuous deglaciation, with glacial Lake Agassiz as the most prominent and voluminous lake, existing for more than 6000 years. The potential importance of its drainage history was underlined by Broecker et al. (1989)¹⁶ who suggested that a major freshwater outburst from Lake Agassiz via the St. Lawrence River at the onset of the Younger Dryas may have disturbed the Atlantic Meridional Overturning Circulation resulting in the last major deglacial cooling stage. Since then, numerous studies aimed to reconstruct the deglacial evolution of the LIS and the routing of meltwater around the Younger Dryas^6,7,10,11 (review in Teller, 2013¹). While one of the largest lake level drops of ~90 m of Lake Agassiz could be tied to the onset of the Younger Dryas^16,17,18,19, suggesting rapid lake drainage at this time, it is still strongly debated whether the routing along the St. Lawrence River was possible during the Younger Dryas, since terrestrial landforms like elevated unglaciated surfaces or moraines as well as blockage by the LIS potentially prevented drainage of Lake Agassiz into the Atlantic^12,20,21. Additionally, terrestrial evidence of a southern or northwestern outflow at the onset of the Younger Dryas do not provide a uniform and coherent picture². While identification and dating of ancient strandlines and gravel beds represents a major advance in understanding the lake level history, the spatial and temporal evolution of the major drainage routes at the onset and throughout the Younger Dryas remain enigmatic. More recent studies argued that the northwestern outlet around the Fort McMurray area and through the Mackenzie River was the most likely drainage route of Lake Agassiz during that time^2,6,7,11, although the terrestrial evidence is either lacking or debated^1,5. Missing terrestrial evidence in the Fort McMurray area can be ascribed to (1) erosion of terrestrial landforms by the severity of the corresponding flood^2,12, (2) slow re-growth of the vegetation after the flood, due to the harsh climate of the Younger Dryas resulting in a lack of radiocarbon-datable organic material^1,22,23, and (3) as yet undiscovered (missing) field evidence. Since terrestrial records hence provide insufficient hints for the reconstruction of the deglacial evolution of Laurentide meltwater discharge, marine records from the Gulf of Mexico^3,4, Laurentian Fan^5,24,25 and Beaufort Sea^6,26 have been investigated aiming to provide the missing supportive evidence. Compared to terrestrial archives marine proxies are less sensitive to even small changes in freshwater runoff and weathering regime but often preserve continuous and well dated information about the most extreme events, which are often the most important triggers for climatic perturbations.

As early as 1975 Kennett and Shackleton²⁷ identified low δ¹⁸O in planktonic foraminifera in the Gulf of Mexico as indicators of deglacial freshwater input from the melting Laurentide Ice Sheet. With the same method, a major freshwater event associated with the Younger Dryas in front of the Mackenzie Delta was identified⁶. However, planktonic δ¹⁸O records suffer from important limitations: First, if the marginal environment becomes too fresh the planktonic foraminifera population decreases or can entirely disappear. Second, the δ¹⁸O signatures incorporated by planktonic foraminifera may be less sensitive to record (fast) changing environments as tens to hundreds of individuals are needed for analysis. Moreover, a δ¹⁸O signal does not bear any information about the source of the freshwater signal. In contrast, the inorganic runoff signal preserved in the Pb isotopic composition of sedimentary authigenic Fe–Mn–oxyhydroxides has been shown to be more sensitive than δ¹⁸O, while additionally providing information towards the source area of continental runoff supplied to a marginal marine basin²⁵.

In this light, Pb isotopes in the authigenic sedimentary phase are a promising sensitive tracer for changes in the runoff signal, especially on glacial/interglacial timescales²⁸. The Pb isotopic signature is controlled by two major processes: First, the detrital signature of any sediment is typically defined by the isotopic bedrock signature of the sediment source (i.e., its provenance), and this isotopic source signature is often also incorporated in the authigenic phase. Second, Pb isotopes are powerful tracers for prevailing weathering conditions due to their variable enrichment in accessory uranium- and thorium-enriched mineral phases, which commonly yield more radiogenic (i.e. higher) Pb isotope compositions than the bulk rock signature. These accessory minerals have a different resistance and dissolution behavior with respect to chemical weathering processes (temperature and precipitation)^29,30,31. As a result, the Pb isotopic signal released during early chemical weathering can be incongruent, with continental runoff being either more radiogenic (in ²⁰⁶Pb/²⁰⁴Pb) or less radiogenic (in ²⁰⁸Pb/²⁰⁴Pb) than the continental bulk Pb isotope signature. Therefore, the bedrock signature as well as the weathering intensity defines the isotopic signature that is transported to the oceans. In the oceans this isotopic signal can be incorporated into authigenic sedimentary phases during the formation of Fe–Mn–oxyhydroxides. The combined analysis of the detrital and authigenic Pb isotopic signature can hence be used to resolve sediment sourcing as well as the prevailing physico-chemical weathering regime in the continental source area. Differences in the detrital and authigenic phase of the same sediment are a powerful indicator for pronounced incongruent weathering or different sourcing of (i) the sediment (detrital material) and (ii) runoff (archived in the authigenic phase). Employing this method, the St. Lawrence River has previously been identified as a major freshwater route during the late deglacial with a particularly strong multi-centennial increase of freshwater routing into the northwestern Atlantic during and after the Younger Dryas²⁵, although a low δ¹⁸O meltwater pulse reflecting the diversion of water from the Gulf of Mexico has never been clearly identified²⁴.

Here we analyzed the Pb isotopic composition recovered from the authigenic and detrital phase of three locations that were cored on cruise HLY1302 (JPC-9; JPC-15/27; JPC-19; Fig. 1) from the Beaufort Sea margin and Amundsen Gulf. With these high-resolution records covering the period from the late Heinrich Stadial 1 to today, we are able to present a detailed reconstruction of freshwater routing into the Arctic Ocean since the Bølling–Allerød. These results provide updated key constraints on freshwater routing from Lake Agassiz into the Arctic Ocean shortly before the onset of the Younger Dryas which has been postulated as one mechanism capable of causing the Younger Dryas cooling^6,13.

Fig. 1: Overview map of North America and the study area.

a Reconstructed ice margins are indicated for climatic key periods (LGM, onset of the Bølling–Allerød warm period at 14.5 ka BP and onset of the Younger Dryas cold period at 13 ka BP)²¹. Blue shaded area indicates the present-day catchment area of the Mackenzie River. The pale gray area depicts the area Lake Agassiz-Ojibway have occupied in the course of the late Deglacial and early Holocene melting history¹⁹. The red frame highlights the area displayed at expanded scale in panel b. b Detailed map of the Mackenzie Delta region. Core sites investigated in this study are displayed as red circles. Ice margins are the same as in panel a. The blue line indicates the Mackenzie River in its modern bed. The flow direction and position of the Shelfbreak Current and Beaufort Gyre are indicated by arrows (light blue and brown, respectively) .

Results

The authigenic as well as the detrital Pb isotopic composition of sediment core JPC-15/27 is marked by a first shift towards more radiogenic (higher) values in ²⁰⁶Pb/²⁰⁴Pb and ²⁰⁸Pb/²⁰⁴Pb at 1330 cm core depth marking the onset of the Bølling–Allerød, the first deglacial warm period around ~14.7 ka BP (Figs. 2 and 3). The authigenic isotopic composition becomes slightly more radiogenic towards the mid-Bølling–Allerød. Throughout the Bølling–Allerød the authigenic and detrital ²⁰⁶Pb/²⁰⁴Pb signals are very similar. In the same interval, the detrital signature of ²⁰⁸Pb/²⁰⁴Pb traces the authigenic trend but with an offset to more radiogenic values. The authigenic Pb isotopic composition experienced its strongest radiogenic excursion just before the transition into the Younger Dryas at 13.1 ka BP (530 cm core depth). In the course of the Younger Dryas and the early Holocene until 8.4 ka BP (150 cm core depth) the authigenic Pb isotopic signature recorded three more radiogenic excursions (Fig. 4) with decreasing magnitude of change but increasing duration. During this interval the detrital signature of ²⁰⁶Pb/²⁰⁴Pb diverges from the authigenic one to less radiogenic values. For ²⁰⁸Pb/²⁰⁴Pb the authigenic signature instead converges and become more similar to the detrital signature. After 8.4 ka BP the authigenic isotopic signature stabilizes at a ²⁰⁶Pb/²⁰⁴Pb ratio of 19.3 and ²⁰⁸Pb/²⁰⁴Pb ratio of 39.1, which are similar to values seen during the early Bølling–Allerød. During this period, the detrital phases of both isotope ratios diverge again from the authigenic signature to more radiogenic values.

Fig. 2: Pb isotopic signatures of Beaufort Sea sediment cores.

Authigenic (red) and detrital (blue) Pb Isotope records (²⁰⁶Pb/²⁰⁴Pb) for sediment cores JPC-15/27 (a), JPC-19 (b), and JPC-9 (c). Replicated authigenic compositions in (a) are shown as orange diamonds. Arrows (gray and black) indicate radiocarbon ages. Colored areas depicts the timing of major intervals explained in the discussion. The IRD Unit was defined by Klotsko et al. (2019)²⁶.

Fig. 3: Geochemical and bulk sediment properties of sediment core HLY1302 JPC-15/27.

Data are shown against depth (cm) due to the large sedimentation rates >1 m/ka. Dashed lines mark the position of the climatic periods Bølling–Allerød and Younger Dryas. a and b Authigenic (red) as well as detrital (dashed blue) ²⁰⁸Pb/²⁰⁴Pb (a) and ²⁰⁶Pb/²⁰⁴Pb (b) signature. Orange diamonds display results of replicate analyses of the authigenic samples. c Stable oxygen isotope data (δ¹⁸O) in ‰ of planktonic foraminifera N. pachyderma s. and d magnetic susceptibility⁶. e ARM_30mT record of JPC-15/27 complemented with f the ARM_30mT/ARM_0mT ratio of the same samples. g Ca/Sr ratio of the authigenic phase.

Fig. 4: Geochemical and bulk sediment properties of core HLY1302 JPC-15/27 for the late Bølling–Allerød, the Younger Dryas and early Holocene.

Data are shown against age (ka). a and b Authigenic (red) as well as detrital (dashed blue) ²⁰⁸Pb/²⁰⁴Pb (a) and ²⁰⁶Pb/²⁰⁴Pb (b) signature. c Stable oxygen isotope data (δ¹⁸O) in ‰ of planktonic foraminifera N. pachyderma s. and d magnetic susceptibility⁶. e ARM_30mT record of JPC-15/27 complemented with f the ARM_30mT/ARM_0mT ratio of the same samples. g NGRIP δ¹⁸O⁶⁹. Light yellow areas indicate major Lake Agassiz drainage events through the Mackenzie River.

ARM_30mT data of JPC-15/27 mirror the evolution of the magnetic susceptibility data⁶, indicating both signals are primarily driven by concentration of ferrimagnetic minerals (Figs. 3 and 4). Differences between the ARM_30mT and magnetic susceptibility are driven by changes in the relative proportion of fine-grained magnetic minerals, as also tracked by the coercivity ratio ARM_30mT/ARM_0mT. The higher values of the coercivity ratio below 1330 cm core depth (~14.7 ka BP) indicate a magnetically fine mineral assemblage that transitions to a coarser assemblage with less variability since the onset of the Bølling–Allerød. This shift in magnetic coercivity is not accompanied by a persistent shift to coarser physical particle size, but instead mirrors a shift in geochemical provenance indicators (Supplementary Fig. S1). Given that higher coercivity is associated with fine magnetic particles (see the “Methods” section), these are further associated with a less radiogenic source with high Ti/Ca and low Ca/Sr (Fig. 3; Supplementary Fig. S1).

Discussion

Initiation of enhanced Arctic meltwater discharge

The most remarkable characteristics of JPC-15/27 are the major shifts in the Pb isotopic composition at the onset of the Bølling–Allerød as well as during the transition into the Younger Dryas which have also been recorded as low δ¹⁸O in planktonic foraminifera in the same core⁵ (Fig. 3). Both δ¹⁸O excursions have been interpreted as being controlled by major flood events from the retreating LIS. At the onset of the Bølling–Allerød warm period coeval trends in the authigenic and detrital Pb isotopic composition (Fig. 3) support the suggested sharp change in sediment and continental runoff sourcing^6,26. Sediment deposited at Site JPC-15/27 before the Bølling–Allerød consisted mainly of transported IRD predominantly from the Amundsen Gulf^26,32. However, the detrital Pb isotopic composition of JPC-15/27 is markedly different (e.g. lower for ²⁰⁶Pb/²⁰⁴Pb) than detrital signatures found in coeval pre-Bølling–Allerød-aged sediments in the Amundsen Gulf (JPC-19; Supplementary Fig. S2). This isotopic difference suggests that sediments deposited before the Bølling–Allerød at Site JPC-15/27 do not exclusively originate from the Amundsen Gulf. Sea ice, icebergs and associated IRD from the Amundsen Gulf was directed to the central Beaufort Sea before these became transported westward within the clockwise circulation of the Beaufort Gyre³³. Consequently, Site JPC-15/27 was only partially influenced by material arriving from the Amundsen Gulf explaining the different isotopic signatures of Site JPC-15/27 and JPC-19 sediments.

With the onset of the Bølling–Allerød the area between Amundsen Gulf, Mackenzie Delta and the Great Bear Lake became gradually deglaciated, thereby allowing unhindered freshwater drainage into the Arctic Ocean through the Mackenzie valley²¹ (Fig. 1). This substantial change in the hydrological system was initiated by a first strong freshwater pulse that is traceable by low planktonic foraminifera δ¹⁸O⁶ and a first shift to more radiogenic Pb (Fig. 3). With the activation of this freshwater route into the Arctic Ocean the amount of runoff substantially increased and sediments from further inland of northern Canada could now be transported to the Beaufort Sea resulting in an increase in regional sedimentation rates^6,26. Additionally, changes in bulk sedimentary magnetic susceptibility, ARM_30mT and ARM_30mT/ARM_0mT as well as authigenic Ca/Sr trace this provenance change of delivered sediment (Fig. 3). For example, authigenic Ca/Sr is a qualitative indicator for detrital carbonate, which is exposed in the northern Canadian section of the interior platform that is built up of Devonian limestone³⁴ and can now be transported to the Arctic. Furthermore, the shift in ARM_30mT/ARM_0mT is indicative for a change from a fine to relatively coarser magnetic mineral assemblage at the lithologic transition from the physically coarse lower IRD unit to a fine laminated unit above^6,26 (Supplementary Fig. S1) and, in the context of other proxies, is therefore indicative for a provenance change (see the “Methods” section). The opposite variability in physical particle vs. magnetic grain sizes at this transition is consistent with the provenance/source change from near coastal material to material from further inland that is indicated by the detrital Pb isotopic ratios and authigenic Ca/Sr. In contrast, the magnetic susceptibility and ARM_30mT data increase at the boundary between the two lithologic units, coincident with an increase in >425 μm lithic clasts and bulk physical particle size^6,26, tracing a short increase in particle size and IRD around the time of the δ¹⁸O defined freshwater event. This first abrupt freshwater event at the onset of the Bølling–Allerød, which was first observed by Keigwin et al. (2018)⁶, equally marks this abrupt and persistent shift in provenance.

In the course of the Bølling–Allerød the deglaciated area increased exposing larger parts of today’s catchment area of the Mackenzie parts of which were subsequently covered by proglacial lakes^35,36. These lakes continuously supplied freshwater via the Mackenzie valley to the Arctic Ocean alongside the delivery of sediment from the same region as evidenced by the quasi invariant behavior of all sedimentary proxies in the course of the Bølling–Allerød for JPC-15/27 (Fig. 3). Only the Pb isotopic signatures (authigenic and detrital) became successively more radiogenic towards the mid-Bølling–Allerød (1000–800 cm core depth). This gradual change could have been caused by intensifying incipient chemical weathering of proglacial sediments and bedrock releasing predominantly a radiogenic signature recorded in the authigenic phase^30,33. However, the associated trend to more radiogenic ²⁰⁸Pb/²⁰⁴Pb of the detrital and authigenic phases argue against increasing incongruent weathering leading to this more radiogenic Pb isotopic signal and rather supports a change of sediment provenance. The gradual Pb isotopic shifts are accompanied by a slight trend in ARM_30mT/ARM_0mT to lower values during the Bølling–Allerød (~0.05 difference; Supplementary Fig. S3), although the amplitude of this variation is much smaller and not as obvious as the previous provenance shift described above.

The Bølling–Allerød-aged sediments from the Amundsen Gulf (JPC-19) are dominated by a distinct radiogenic excursion in the detrital and strong radiogenic excursion in the authigenic phase lasting for at least one kyr, which in the following will be referred to as Amundsen Gulf event (Fig. 2; Supplementary Fig. S2; first 500 cm). The retreating Amundsen ice stream likely made the Amundsen Gulf a major outflow route of freshwater (in form of icebergs and meltwater) during the Bølling–Allerød^32,37,38. As discussed above, sediment and freshwater transport from the Amundsen Gulf was generally directed to the central Beaufort Sea leaving site JPC-15/27 mainly unaffected. However, the extremely distinct Amundsen Gulf event, observed during the Bølling–Allerød (Fig. 2, Supplementary Fig. S2), has also been identified as a sedimentary unit outside the Amundsen Gulf proximal to site JPC-15/27 (site JPC-25²⁶) by magnetic susceptibility and seismic data. Although this event was not recorded in magnetic susceptibility and seismic data at Site JPC-15/27, the Pb isotopic signatures during the mid-Bølling–Allerød potentially trace this Amundsen Gulf (freshwater-)event influencing the whole region (Fig. 2a, b). If causally related, this would show that Pb isotopic compositions in the detrital and authigenic phases are more sensitive to even minor sedimentary (provenance-) changes which are not recognizable in classical bulk sedimentary proxies at Site JPC-15/27. This sensitivity for shifts in sediment sourcing is due to the diverse geology in crustal sections bordering the Mackenzie delta and Amundsen Gulf^34,39. The influence of Amundsen Gulf sediments on JPC-15/27 decreased again towards the end of the Bølling–Allerød marked by slightly less radiogenic values.

Younger Dryas freshwater pulses and their sources

The end of the Bølling–Allerød warm period is marked by the most radiogenic authigenic ²⁰⁶Pb/²⁰⁴Pb as well as ²⁰⁸Pb/²⁰⁴Pb excursions in JPC-15/27 associated with strong excursions in the magnetic susceptibility, ARM_30mT records (Figs. 3 and 4), and physical particle size (Supplementary Fig. S1)^6,26. It is tempting to interpret this proxy behavior as a major change in the depositional regime and sedimentary source. However, the detrital Pb isotopic signatures do not change alongside the authigenic signal in ²⁰⁶Pb/²⁰⁴Pb and only partially in ²⁰⁸Pb/²⁰⁴Pb throughout the late Bølling–Allerød and the early Younger Dryas. In other words, a clear provenance change as seen at the onset of the Bølling–Allerød equally affecting the authigenic and detrital signature is not observed here. Taken together these lines of evidence argue against a considerable sediment provenance change and are likely driven by the increased sedimentary particle sizes⁶, which do not have to affect the detrital Pb isotope signature, if the delivered sediment is from the same source. Given these constraints, the more radiogenic authigenic Pb isotope excursions which are decoupled from the detrital Pb isotope trends rather represent a runoff signal from another origin than glacial Lakes Mackenzie and McConnell located within the deglaciated Mackenzie catchment. Earlier reconstructions of the LIS suggested the Fort McMurray area as an important spillway for glacial Lake Agassiz to the Arctic^7,8,22,40. However, due to the lack of terrestrial evidence it is still debated if this drainage route was opened at the onset of the Younger Dryas or not^2,7. In this context, the authigenic Pb isotopic excursion in core JPC-15/27 very likely represents the first direct evidence for Lake Agassiz outflow through the Mackenzie valley via the Athabasca spillway in the Fort McMurray area at the transition from the Bølling–Allerød into the Younger Dryas. Interestingly, this outflow event, which defines the most prominent authigenic Pb isotopic excursion seen in the whole deglacial record of the Mackenzie area, predates the Younger Dryas (Fig. 4) and would therefore agree with the hypothesis of an Arctic Ocean freshening as a plausible cause for the Younger Dryas cooling^6,10,11,13. Strikingly, our authigenic Pb isotope peak also predates the δ¹⁸O excursions of JPC-15/27⁶, which was interpreted as an indicator for the same freshwater event. This circumstance reflects the sensitivity of authigenic Pb isotopes to runoff (weathering rate and freshwater supply) changes. The high sensitivity of our Pb isotopic proxy for variable freshwater supply in this geographic setting is further underlined by recurring smaller isotopic excursions during the Younger Dryas and early Holocene until 8.5 ka BP (Fig. 4). While δ¹⁸O does not suggest any further freshening events after the early Younger Dryas, the authigenic Pb isotopic signatures recorded at least two more distinct radiogenic excursions around 12, 11.2 ka BP and a potential minor one starting at 9.7 ka BP (Fig. 4). The two excursions around 12 and 11.2 ka BP are in good temporal accordance with findings from gravel bed deposits in the Mackenzie Delta¹¹. As continuous freshwater outflow persisted from northern glacial lakes (e.g. Lake Mackenzie and Lake McConnell), the radiogenic Pb isotope excursions at 12 and 11.2 ka BP might also be linked to overflow events of Lake Agassiz through the Athabasca spillway¹¹. Discussions about the opening of this spillway and potential ice re-advances in this area, which have blocked this drainage routes from time to time modulating flood like events, are ongoing^2,7,8,11,22. Furthermore, the source of the last overflow event at 9.7 ka BP is rather unclear and difficult to determine. There are suggestions that overflow of Lake Agassiz was exclusively to the east after 10.6 ka BP¹⁹ while other authors favor minor overflows of Lake Agassiz to the north^6,11. However, isostatic and geomorphological reconstructions of Lake Agassiz during the deglacial period strongly suggests that it was impossible for Lake Agassiz to drain to the Northwest and only the eastern outlet was active after 10.6 ka BP¹⁹. Consequently, the excursions observed at 9.7 ka BP was most likely be caused by overflow of smaller proglacial lakes (e.g., glacial lakes Meadow, McMurray, and McConnell) north of Lake Agassiz and the Athabasca spillway⁷. Irrespective of the freshwater source, such extensive runoff fluxes into the Arctic Ocean could not only have triggered the Younger Dryas but may have also been involved in mechanisms triggering the Preboreal oscillation, a period of climatically unstable conditions^18,26,41.

Whether or not the decreasing amplitudes of the Pb isotopic excursions are driven by decreasing outflow intensities or rather reflect different water mass sourcing for each event cannot be answered with the available data. The formation of the runoff Pb isotopic signal is a complex process, which is mainly controlled by the bedrock signature and weathering/erosion rates (intensity/time). Lake Agassiz was located mainly on granitic bedrock of the Superior Province, which is one of the most radiogenic cratons regarding its Pb isotopic composition in North America⁴² (Supplementary Fig. S4). It could therefore easily generate the observed highly radiogenic Pb isotope signatures seen in the authigenic phase of JPC-15/27 at the onset of the Younger Dryas even under congruent release of crustal Pb during early chemical weathering in the catchment area. However, the weathering of these granitic sequences fails to explain the observed pattern of decreasing radiogenic compositions for each further event unless freshwater additions from other catchments modified the ambient Lake Agassiz Pb isotope signature or the net volume of freshwater runoff became successively smaller with each drainage event. The first outflow event of Lake Agassiz to the Arctic likely was the most voluminous as witnessed in the strongest lake level drop of Lake Agassiz during its whole evolution^2,18,19. Each subsequent drainage event of Lake Agassiz towards the Arctic Ocean has to be considered as a unique event with distinct environmental and geographical boundary conditions. The catchment area of Lake Agassiz changed continuously with the retreating LIS^19,43,44 during which different lithologies were exposed to the atmosphere and therefore to weathering. Furthermore, the exposure time of sedimentary substrate (e.g. glacial flour) to weathering and subglacial supply of meltwater has to be taken into account. As a consequence, the Pb isotopic composition of Lake Agassiz lake water must have been modified continuously in the course of the deglaciation but have to be confirmed in future analysis of Lake Agassiz sedimentary sequences. Last but not least, the isotopic signature of each water drainage event from Lake Agassiz is expected to be modified differently during transport toward the Arctic depending on the above mentioned conditions as well as the time and distance of transport.

In contrast to the Pb isotopic signatures of the authigenic phase, the detrital phase shows a different trend in the course of the Lake Agassiz drainage period (13.1–8.5 ka; Fig. 4). Klotsko et al. (2019)²⁶ assumed that the sedimentary unit deposited during the Younger Dryas originated from Lake Agassiz. However, sediment from Lake Agassiz itself should be more radiogenic due to its location within the Superior Province (Supplementary Fig. S4), but such radiogenic detrital isotope compositions are not observed here. Patterns observed during the Bølling–Allerød with very similar authigenic and detrital ²⁰⁶Pb/²⁰⁴Pb signatures and a more radiogenic detrital than authigenic ²⁰⁸Pb/²⁰⁴Pb are in line with findings of identical detrital and authigenic Pb isotopic compositions in post-glacial lake sediments/deposits that were supplied from a common source⁴⁵. Yet the different ²⁰⁶Pb/²⁰⁴Pb and ²⁰⁸Pb/²⁰⁴Pb evolution of the authigenic and detrital signatures during the main drainage period starting at 13.1 ka BP point to independent controls of the respective signatures during transport. The modern Mackenzie catchment mainly includes lithologies of the Interior Platform³⁴ (Supplementary Fig. S4) and large amounts of sediments derived from the Canadian Shield could therefore only be transported along the Mackenzie Valley to the Arctic Ocean by strong flood events originating from, for example, Lake Agassiz. However, while freshwater can flow nearly unhindered from the source to the oceans, Lake Agassiz sediments are much less effectively transported towards the Arctic Ocean. Sediment that was potentially transported northward from Lake Agassiz was reported to be deposited in glacial lakes such as glacial Lake Mackenzie, effectively acting as a sediment trap hindering these sediments from reaching the Arctic Ocean⁴⁶. Therefore, relative source contributions of freshwater (recorded in the authigenic signal) and terrigenous sediment (detrital signature) reaching the Arctic Ocean during this period have to be considered individually. Additionally, the invariant ARM_30mT/ARM_0mT ratio argues for a generally unchanged sediment composition during this period with respect to the magnetic mineral assemblage. Furthermore, high magnetic susceptibility values and ARM_30mT data during the Younger Dryas rather represent the strong variation/increase in particle size^6,26 caused by the higher transport energy of Lake Agassiz outflow waters as equally seen at the onset of the Bølling–Allerød.

The trends in detrital Pb isotopic signatures during the Younger Dryas and early Holocene is at foremost tracing continuous geographic changes in sediment sourcing. The progressing deglaciation increased the catchment of the Mackenzie drainage route towards today’s Mackenzie catchment area (Fig. 1; Supplementary Fig. S4). The reason for the switch back to more radiogenic detrital Pb isotope compositions (most obvious in ²⁰⁶Pb/²⁰⁴Pb) around 8.5 ka BP cannot ultimately be resolved in this study. It is obvious that the timing coincides with the postulated final and total drainage of Lake Agassiz into the Hudson Bay⁴⁷. However, as discussed above, a connection between Lake Agassiz and the Arctic Ocean was rather unlikely at this time. One possible explanation could be a teleconnection between Lake Agassiz and the Mackenzie area due to the instability of the LIS. The abrupt end of Lake Agassiz can be linked to the collapse of the Hudson Bay Ice Saddle and the final outflow of Lake Agassiz releasing some 163,000 km³ of freshwater¹⁷ into the Hudson Bay⁴⁷ that was previously stored along the southern margin of the LIS. The collapse of the Hudson Bay Ice Saddle and the drainage of Lake Agassiz might have caused further instabilities and breakup of the LIS that have caused re-routing of meltwater from northern glacial Lakes, which was previously directed to the Arctic Ocean, now to the former Lake Agassiz area and the Hudson Bay. Such a sequence of events would have the capability of terminating extensive freshwater supply from the southern Mackenzie catchment into the Arctic Ocean. As a result, the early Holocene Mackenzie drainage area was established and sediment supply similar to today towards the Arctic Ocean was manifested within this area. The same event at 8.5 ka BP was also recorded in the authigenic isotopic signature, at this point resulting in a return to remarkably invariant Holocene unradiogenic values and a shift in the ARM_30mT/ARM_0mT ratio. The change towards less radiogenic authigenic Pb isotope compositions not only marks the onset of uniform sedimentation detached from glacial processes inland after 8.5 ka BP, but the more radiogenic Pb isotope compositions before that time bear evidence for the very dynamic runoff environment with an occasionally more or less pronounced contribution from glacial lakes.

Analogously to JPC-15/27, sediment core JPC-9 to the west of the Mackenzie Delta witnessed a similar deglacial trend (from the Younger Dryas until today; Fig. 2; Supplementary Fig. S2). At this site the detrital phase recorded a nearly constant Pb isotope signal pointing to regional steady sediment sourcing while the authigenic phase yielded most radiogenic values during the Younger Dryas and early Holocene, likely equally influenced by freshwater from Mackenzie outflow. The amplitude of these radiogenic excursions at site JPC-9 are less intense as observed at site JPC-15/27, suggesting preferential freshwater transport to the east and central Beaufort Sea by the eastward directed shelfbreak current given that JPC-9 and JPC-15/27 are located in a similar distance to the Mackenzie Delta (Fig. 1).

Deglacial drainage chronology of the LIS since the Bølling–Allerød

Reconstructing the major Younger Dryas drainage route(s) (e.g., to the Arctic Ocean) was previously shown to be difficult, in particular with regards to its timing^1,2,23,48. Consequently, it is even more challenging to provide a complete and consistent drainage chronology of the melting LIS for the entire deglacial period. With our constraints from the northern drainage route we compiled marine records from the Gulf of Mexico (planktonic δ¹⁸O from core EN32-PC6^3,4), Laurentian Fan (authigenic Pb isotopes and planktonic δ¹⁸O from composite core 26GGC/14GGC^24,25), Hudson Strait (planktonic δ¹⁸O from core MSM45-19-2⁴⁷) and Beaufort Sea margin (planktonic δ¹⁸O of core JPC-15/27⁶ as well as authigenic and detrital Pb isotopes from this study) to provide an updated overview of the major meltwater routing events/periods from the decaying LIS from the onset of the Bølling–Allerød (Fig. 5).

Fig. 5: Chronology of active drainage routes of meltwater from the Laurentide Ice Sheet during the Deglaciation and early Holocene.

a Pb isotope (²⁰⁶Pb/²⁰⁴Pb) data of this study are compared to δ¹⁸O⁶ of the same core (JPC-15/27) as well as ²⁰⁶Pb/²⁰⁴Pb²⁵ data planktonic δ¹⁸O²⁴ from the Laurentian Fan (composite core: 26GGC/14GGC), δ¹⁸O from the Gulf of Mexico (EN32-PC6; brown³, black⁴), δ¹⁸O from the Labrador Sea as an indicator for Hudson Strait freshwaters (MSM45-19-2⁴⁷) and NGRIP δ¹⁸O⁶⁹. b Colored arrows mark the main drainage routes. Numbers displayed in white arrows correspond to phases of different drainage route activation independent from prominent climate periods. The pale gray area depicts the area Lake Agassiz-Ojibway have occupied in the course of the Deglacial and early Holocene melting history. Wavy areas correspond to the main areas of Lake Agassiz and Lake Ojibway during the Bølling-Allerød and Younger Dryas before merging¹⁹. c Chart displays the major active drainage routes and their potential freshwater source. Northern glacial lakes (NGL) are Lake Mackenzie and Lake McConnell; Southeastern Glacial Lakes (SEGL) represent former glacial lakes in the Great Lakes area as well as glacial Lake Ojibway.

After the LGM, in the course of Heinrich Stadial 1 the majority of meltwater from the LIS was transported southward via the Mississippi to the Gulf of Mexico. In the ongoing deglaciation before the Bølling–Allerød, proglacial lakes had formed around the LIS and new drainage routes became active. Due to the establishment of multiple drainage routes from the beginning Bølling–Allerød, we distinguished five phases defined by the major drainage route activity (Fig. 5). During the early Bølling–Allerød (14.7–13.8 ka BP; first phase) the freshwater flux into the Gulf of Mexico increased^3,4,49, sourced from Lake Agassiz as well as glacial lakes in the region of today’s Great Lakes^9,19,49. At the same time a first strong freshwater signal was traceable in the Arctic^6,50 (this study), sourced from retreating LIS sequences in the northern Mackenzie and Amundsen Gulf area³⁶ and the drainage of northern glacial lakes²⁶. In the course of the Bølling-Allerød the freshwater flux into the Arctic Ocean increased due to the collapse of the Amundsen Gulf Ice Stream (see discussion above). The second half of the Bølling–Allerød (13.8–13.1) can be regarded as the second phase of deglacial freshwater routing. The southward outflow from Lake Agassiz increased dramatically resulting in very high runoff rates sustained over several hundred years into the Gulf of Mexico^3,4. At the same time a first resolvable runoff signal was recorded in Laurentian Fan sediments²⁵. It is unclear if this meltwater is from a rather local source of accelerated ice margin retreat²⁵ or has been potentially supplied by an overflow of Lake Agassiz into the Great Lake area and further into the North Atlantic⁹. These enhanced meltwater fluxes during the Bølling–Allerød are the result of accelerated melting of the LIS, a contemporaneously growing volume of Lake Agassiz overflow¹⁹ as well as glacial lakes in the Great Lakes area and finally an ice margin retreat opening new drainage and freshwater routes (at this time the eastern additionally to the southern outlet). The latest Bølling–Allerød marks the beginning of the third deglacial drainage phase that extends well into the Younger Dryas (13.1–12.4 ka BP). The retreating ice margin of the LIS opened the Athabasca Spillway in the Fort McMurray area and allowed major drainage to the Arctic (see discussion above). Marine records of the Gulf of Mexico⁴ and the Laurentian Fan²⁵ witnessed a substantial decrease in freshwater input to these sites. With this first overflow of Lake Agassiz into the Arctic the southern drainage to the Gulf of Mexico became minuscule and was not a major outlet for the remainder of the deglaciation². The role of the eastern drainage route along the St. Lawrence valley during the earliest Younger Dryas is still enigmatic. There is evidence of some freshening with a potential Lake Agassiz source^5,40 but also strong evidences of only limited freshwater routing^6,12,25,51. In light of the findings of our study, it is most plausible that in the transition to the Younger Dryas Lake Agassiz almost exclusively drained to the Arctic Ocean and that freshwater directed along the St. Lawrence valley originated from the Great Lakes. In the course of the Younger Dryas the drainage route to the Laurentian Fan was gradually reactivated²⁵ and initiated the fourth deglacial drainage phase (12.4–8.6 ka BP). During this phase, meltwater from the LIS was routed nearly exclusively along the St. Lawrence valley into the North Atlantic as well as along the Mackenzie River into the Arctic Ocean. While Lake Agassiz continuously supplied freshwater through the St. Lawrence valley into the NW Atlantic, two potential overflow events into the Arctic Ocean took place at 12.0 and 11. 2 ka BP.

With further retreat of the LIS to the North and continuous melting, the second largest glacial lake of North America, Lake Ojibway, formed east of Lake Agassiz during the early Holocene, which later merged with Lake Agassiz to form the largest existing glacial lake during the complete melting history of the LIS^2,9,19. The last drainage phase (8.6–8.2 ka BP) of meltwater from the LIS is characterized by the last and most voluminous freshwater outburst of Lake Agassiz-Ojibway. Due to the collapse of the Hudson Bay Ice Saddle Lake Agassiz-Ojibway drained completely into the Hudson Bay between 8.5 and 8.2 ka BP^47,52 and further into the Labrador Sea causing the 8.2 ka cooling event^53,54. However, evidence was recently presented for large volumes of freshwater reaching the Hudson Bay before the collapse of the ice saddle through subglacial flows^41,55. In either case, meltwater drainage to the other outlets of Lake Agassiz-Ojibway was immediately terminated after the collapse of the Hudson Bay Ice Saddle around 8.5 ka BP. After this event, most North American rivers and catchments attained their recent configuration like the Mackenzie River, witnessed by the invariant Pb isotopic signature at site JPC-15/27 (Figs. 3 and 4).

Overall, the Mackenzie Pb isotope record of JPC-15/27 documents continuous northward directed meltwater routing into the Arctic Ocean from the northern side of the LIS since the onset of the Bølling–Allerød, lasting almost until the final disintegration stages of the Laurentide Ice Sheet around 8.5 ka BP. Our data solve the long-standing conundrum surrounding the activity of a northwestern drainage route of Lake Agassiz into the Arctic Ocean at the onset of the Younger Dryas. The authigenic Pb isotope records present direct evidence of freshwater routing from Lake Agassiz along the Mackenzie valley at the very end of the Bølling–Allerød and therefore support the hypothesis of a freshening Arctic Ocean as a potential trigger for the Younger Dryas cooling. Our northwestern drainage chronology of Lake Agassiz since the Younger Dryas provides insights into the melting chronology of the LIS and its runoff along the southern and eastern drainage routes. In this context, the role of the eastern drainage route along the St. Lawrence valley as well as into the Labrador Sea should be investigated in greater detail to better constrain its importance as outlet routes of Lake Agassiz freshwater during the Younger Dryas.

Methods

Sediment cores and age models

Piston cores JPC-9, JPC-15, JPC-19, and JPC-27 were obtained during the USCGC Healy cruise 1302 along the Canadian Arctic continental margin in the Beaufort Sea as well as the Amundsen Gulf. Core JPC-9 (70.58°N, 142.42°W, 394 m water depth; Fig. 1) is located at the westernmost location of the investigated area and reaches back at least until 13 ka BP⁶. This core was previously dated with three ¹⁴C dates (Keigwin et al., 2018). We obtained two additional radiocarbon ages on planktonic foraminifera N. pachyderma s. (800 cm core depth: 10.15 ka cal BP; 1000 cm core depth: 11.32 ka cal BP). The new ¹⁴C dates were measured at the LARA laboratory at the University of Bern⁵⁶ and tied to the Marine20 calibration curve using the CALIB 8.2 online tool⁵⁷ assigning the identical ∆R-values as presented in Keigwin et al. (2018)⁶ that vary within 0 and 200 conventional radiocarbon years for the deglacial interval and Holocene. Cores JPC-15 and JPC-27 are located at the same position and water depth (71.10°N, 135.13°W, 687 m water depth). During coring, JPC-27 over-penetrated by 205 cm below the sediment surface into the sediment but recovered a longer sediment sequence than JPC-15. Therefore, both cores were aligned by magnetic susceptibility and a composite core was formed⁶ (JPC-15/27). This composite core was dated with 14 radiocarbon ages and dates back at least until 15 ka BP⁶. Core JPC-19 (71.29°N, 126.28°W, 442 m water depth) is located in the central Amundsen Gulf and represents the easternmost site of the study area. This core was dated with four new radiocarbon ages on mixed samples of planktonic and benthic foraminifera due to the scarcity of abundant foraminifera in the sediments. These ¹⁴C dates were also measured at the LARA laboratory at the University of Bern⁵⁶ and calibrated using the Calib 8.2 online tool⁵⁷ with ∆R-values of Keigwin et al. (2018)⁶ (Fig. 2; Supplementary Figs. S2 and S5). A detailed lithological description of all investigated cores can be found in Keigwin et al. (2018)⁶ and Klotsko et al. (2019)²⁶.

Isotope and elemental analyses

Lead isotopic compositions were analyzed in the authigenic Fe–Mn–oxyhydroxide fraction as well as in the detrital phase. The authigenic phase was extracted from ~0.5 g weighed out sediments with a weak reductive leaching solution⁵⁸ (0.005 M hydroxylamine hydrochloride, 1.5% acetic acid, 0.003 M Na-EDTA solution buffered to pH 4 with 25% ammonia solution), yet exposing the sediment only ~30 s to the reductive leaching solution on a vortex mixer as applied for Southern Ocean sediments⁵⁹. Before digestion of the detrital sediment the remaining authigenic Fe–Mn–oxyhydroxides were removed with a tenfold more concentrated reductive solution⁶⁰ (0.05 M hydroxylamine hydrochloride, 15% acetic acid, 0.03 M Na-EDTA solution buffered to pH 4 with 25% ammonia solution) applied for at least 12 h. Afterwards sediments were ground and homogenized before a sub-sample of ~50 mg was weighed out prior to dissolution in a microwave system in a mixture of concentrated HNO₃ an HBF₄. The Al/Pb ratios of the detrital and extracted authigenic phase was monitored to assess any detrital contribution to the authigenic signal that should result in elevated Al/Pb (»100)⁶⁰ (Supplementary Fig. S5). After leaching and/or dissolution, aliquots were taken for elemental analyses while Pb was separated by ion chromatography exchange resin AG1-X8⁶¹. Purified samples were analyzed for their Pb isotopic composition using a Thermo Scientific Neptune Plus MC-ICP-MS at GEOMAR Helmholtz Centre for Ocean Research Kiel, Germany. Mass bias correction was performed by Tl-doping^62,63 aiming for a Pb/Tl ratio of ~4. The secondary standard USGS NOD-A-1 were reproduced with an accuracy of 0.0055 for ²⁰⁶Pb/²⁰⁴Pb and 0.015 for ²⁰⁸Pb/²⁰⁴Pb (n = 25). The total procedure blanks were on the order of 500 pg which represents less than 0.1% of the sample signal and are therefore negligible. Detailed information about the measurement procedure can be found in Süfke et al. (2019)⁴⁵ and Huang et al. (2021)⁵⁹.

Finally, elemental concentrations of major elements (Li, Al, P, Ca, Ti, Mn, Fe, Sr, Ce, Nd, Pb, Th and U) were measured for the extracted authigenic phase and the digested detrital samples with an Agilent Series 7500 ICP-MS at GEOMAR Kiel. The relative reproducibility of element measurements is better than 2% (2SD) based on repeated determination of the elemental composition of the secondary standard USGS NOD-A-1.

ARM analyses

Anhysteretic remanent magnetization (ARM) data were generated for u-channel samples from sediment core JPC-15/27. ARMs were applied using a 100 mT peak alternating field (AF) with 0.05 mT bias field on a 2 G Enterprises^TM model 755-1.65UC superconducting rock magnetometer (SRM) with inline AF coils optimized for u-channel samples. The magnetizations were measured for sediments before (ARM_0mT) and after treatment with a 30 mT peak AF (ARM_30mT). Variation in both ARM and magnetic susceptibility data primarily reflect the concentration of magnetic minerals, however, ARM is particularly sensitive to the concentration of fine grained magnetite while magnetic susceptibility can be more sensitive to coarser grained magnetite^64,65. The ratio ARM_30mT/ARM_0mT track the bulk magnetic coercivity. In the case of constant magnetic mineralogy, this ratio is proportional to the grain size of magnetic particles with higher values indicating a higher proportion of finer magnetite particles⁶⁶. However, magnetic coercivity can often be decoupled from the physical particle grain size of the sediments and reflect sediment provenance changes controlled by magnetic mineral assemblage changes and to magnetic particles existing as inclusions in larger silicate minerals⁶⁷. Thus sediment magnetic mineral concentration and coercivity can reflect sediment provenance, transport, and/or diagenesis^67,68, which in turn can be a function of changes in sediment source, erosion/weathering and hydrological regimes. Therefore, the combined consideration of ARM_0mT, ARM_30mT, Magnetic Susceptibility and the ARM_30mT/ARM_0mT ratio can provide additional insights to disentangling the effects of changing particle size or sediment source along with sedimentological and geochemical data.