Deglacial to Holocene variability in surface water characteristics and major floods in the Beaufort Sea

Surface water characteristics of the Beaufort Sea have global climate implications during the last deglaciation and the Holocene, as (1) sea ice is a critical component of the climate system and (2) Laurentide Ice Sheet meltwater discharges via the Mackenzie River to the Arctic Ocean and further, to its outflow near the deep-water source area of the Atlantic Meridional Overturning Circulation. Here we present high-resolution biomarker records from the southern Beaufort Sea. Multi-proxy biomarker reconstruction suggests that the southern Beaufort Sea was nearly ice-free during the deglacial to Holocene transition, and a seasonal sea-ice cover developed during the mid-late Holocene. Superimposed on the long-term change, two events of high sediment flux were documented at ca. 13 and 11 kyr BP, respectively. The first event can be attributed to the Younger Dryas flood and the second event is likely related to a second flood and/or coastal erosion. The Beaufort Sea was nearly ice-free during the transition from the last deglacial to the Holocene, a period in which two episodes of high sediment flux suggest major glacial flood events, according to high-resolution multi-proxy biomarker records.

A rctic sea ice, characterized by strong seasonal variations, is a critical component in the global climate system as sea ice has direct impact on climate change by regulating energy exchanges between the atmosphere and ocean. The variability of sea ice contributes to changes in the albedo effect, the freshwater system, and the surface energy budget in the Arctic Ocean [1][2][3][4] . In a positive ice-albedo feedback, increases of absorbed energy result in further sea-ice melting 5,6 . Increased melting and sea-ice export can slow down the Atlantic Meridional Overturning Circulation (AMOC) as well as the poleward heat transport 7 . Due to these complex feedback processes, the Arctic is both a contributor to climate change and a region that is most strongly affected by global warming 3,4,8,9 . Arctic sea ice has reduced drastically over the past few decades, recognized in both satellite observations and model simulations 8,10,11 . Since 1980, the sea-ice cover has had a mean annual areal reduction of~20% and an even stronger decrease of 30% in September 7 . This decrease in sea ice is also suggested to be the driver of frequent cold extremes over Eurasia in the past two decades 12 . Historical observations spanning the last few decades are deficient in length to decipher the processes of accelerated sea-ice retreat, thus longer-term and high-resolution proxy records of paleoenvironmental changes are needed. During transition from the last deglacial to the Holocene, the climate system underwent numerous abrupt changes, particularly the Bølling/Allerød (B/A) interstadial and the Younger Dryas (YD) stadial. The Holocene, not concluded yet, has experienced a significantly warmer period during its early stage (Early Holocene Thermal Maximum) 13 . Causes of the abrupt climate change are not fully understood, hence special emphasis on the transition from the deglacial into the Holocene is of significance to understand the forcing mechanisms of the climate system. In this context, the Arctic marginal seas characterized by strong seasonal variability in sea-ice cover, primary productivity, and terrigenous (riverine) input are very sensitive to environmental changes and thus of major importance for paleoclimate reconstructions. One of these marginal seas is the Beaufort Sea ( Fig. 1), which we focus on in this paper.
Today, the Mackenzie River is the largest fluvial source of water and sediment into the Beaufort Sea, characterized by water discharge of 316 km 3 yr −1 and sediment flux of 124-128 Mt yr −1 (Ref. 14 ; Figs. 1 and 2). Since the drainage system was established during the late Wisconsinan (i.e., Marine Isotope Stage 2) as a result of glacial erosion, it was an outlet of glacial meltwater during the last deglaciation 15 . During this period, a cold episode known as YD coincides with a significant reduction of AMOC caused by increased freshwater flux into the North Atlantic deepwater formation region 16 . The freshwater discharge has been linked to an outburst from the Lake Agassiz, and there is an ongoing debate about the pathway of this freshwater discharge, i.e., whether it drained first into the Arctic Ocean and then to the Atlantic Ocean or whether it drained directly into the Atlantic Ocean [17][18][19][20] (Fig. 3). More recently, field data indicate the freshwater of Lake Agassiz could not drain northward during the YD 21 , while model simulations suggest the meltwater drainage via Mackenzie River (northwestern outlet) as a competitive candidate triggering the YD cooling 17,22 . Following the YD cold period, field data from the Fort McMurray and the Mackenzie Delta point to a post-YD flood from proglacial lakes, which drained through the Mackenzie River to the Arctic Ocean during the early Holocene 23,24 . The timing of the flood coincides with the Preboreal Oscillation (PBO) and thus it has been proposed as the trigger of PBO 25 . Above-mentioned events are of significance to climate change, but corresponding marine records are scarce. Therefore, high-resolution records are needed to identify these events and the timing of their inception.
Here, we apply a biomarker approach on a well-dated sediment core from the Beaufort Sea directly off the Mackenzie River ( Fig. 1), covering the time interval of the last deglaciation and the Holocene. Multiple biomarker proxies are used to reconstruct the changes in sea ice, sea surface temperature (SST), primary productivity, and terrigenous input. Based on our records, we demonstrate that (1) the Beaufort Sea region was nearly ice free with variable SSTs during the last deglaciation and the early Holocene, (2) sea-ice cover developed during the mid-late Holocene, coinciding with a drop in terrigenous sediment flux, SST, and primary production, and (3) two major events, characterized by prominent maxima in sediment flux, occurred at 12.83 ± 0.15 and ca. 11 kyr BP. Whereas the former is related to the YD flood event, the origin of the second event is related to a post-YD flood and/or coastal erosion.
The chronology of Core ARA04C/37 was constrained by accelerator mass spectrometry (AMS) 14 C dating on calcareous foraminifera (planktic and mixed species; Fig. 4a and Supplementary Table 1, see "Methods") and, in the uppermost centimeters, excess 210 Pb. The age-depth model is established by a combination of AMS 14 C dates of cores ARA04C/37 and JPC15 ( Fig. 4) 19 , based on the good correlation of AMS 14 C dates and magnetic susceptibility (Figs. 4a and 5a, b) from both cores (see "Methods"). Ages in this paper are provided as calendar years BP. The age-depth relation reveals a mean sedimentation rate of 27 and 100 cm kyr −1 in the upper and lower 300 cm, respectively, with high sedimentation rates in two layers near 13 kyr BP (450-500 centimeters below seafloor (cmbsf)) and 11     transition in high resolution, including the B/A and YD intervals. This enables reconstructions of sea ice, primary productivity, SST, and freshwater discharge in the Beaufort Sea, with special emphasis on the flood events during the transition, in great detail. The sea-ice biomarker IP 25 28,29 in combination with openwater phytoplankton biomarkers (e.g., dinosterol), the so-called "PIP 25 " index, allows semi-quantitative reconstructions of sea-ice concentrations (e.g., Refs. 30,31 ). The newly developed ring index (RI-OH′) of hydroxylated isoprenoid glycerol dialkyl glycerol tetraethers (OH-GDGTs) has been used as a potential tool to reconstruct SST in cold environments (ca. <15°C) 32 . The biomarkers brassicasterol and dinosterol are indicators for marine phytoplankton, whereas ß-sitosterol and campesterol are indicators for terrigenous input [33][34][35] . The brassicasterol in some settings has significant sources from freshwater 33,36 , which is the case in our study (discussed below), and hence we focus on dinosterol to reconstruct marine primary production ( Supplementary Fig. 1a, b). The branched GDGTs (b-GDGT) are indicative for terrestrial input [37][38][39][40] , and the F C32 1,15 index (fractional abundance of C 32 1,15-diol) is indicative for aquatic riverine input, as the C 32 1,15diol is found particularly abundant in freshwater systems [41][42][43][44] . More details can be found in the method section.
Deglacial-Holocene changes in sea ice and primary production. During the deglaciation and the early Holocene (~14-8 kyr BP), minimum concentrations of IP 25 (close to the detection limit) and low PIP 25 values point to dominantly ice-free conditions (Figs. 6a and 7c and Supplementary Fig. 1c). In such a situation, open-water conditions would be expected to promote primary production. However, concentrations of dinosterol (~1-5 µg g TOC −1 ) were low to moderate, implying low primary production during the B/A interstadial (Fig. 6c). This may be explained by the fact that the Beaufort Sea itself is an oligotrophic system 6 , and the Bering Strait was still closed at that time (>13 kyr BP), hence preventing the inflow of nutrient-rich Pacific Water 45,46 . Primary productivity shortly increased during the early YD (Figs. 6c and 7e), probably promoted by higher nutrient levels supplied by the YD flood (see below). In the early Holocene, the primary productivity experienced a rapid increase at ca. 11 kyr BP as well as subsequent moderate increases during the early Holocene (Fig. 7e), possibly caused by high terrigenous input (including nutrients) and the total inundation of Bering Strait, respectively 45 . Furthermore, maximum summer insolation 47 may have caused shorter sea-ice seasons at the core site. More open waters and longer duration of ice-free conditions were favorable for dinosterol production in spring and summer.
During the mid-late Holocene (8-0 kyr BP), sea ice expanded, as evidenced by elevated IP 25 and PIP 25 values (Figs. 6a and 7c), which reveals an evolution from dominantly ice-free conditions (PIP 25 < 0.2) to more extended ice-cover conditions (PIP 25 > 0.8) due to decreasing summer insolation. In turn, expanded sea-ice cover and longer sea-ice duration reduced primary production since the middle Holocene (ca. 8-2.5 kyr BP) as clearly reflected in decreased dinosterol accumulation rates (Fig. 7e). The total organic carbon (TOC)-normalized dinosterol concentrations differ from dinosterol accumulation rates during the middle Holocene, increasing between ca. 8-5 kyr BP (Fig. 6c). This is attributed to simultaneous decreases in marine primary production and terrestrial input, due to sea-ice formation and decrease in river discharge, respectively. During this phase, the reduction of terrestrial organic carbon input was probably larger than that of marine organic carbon input, resulting in increased TOCnormalized dinosterol concentrations. During the late Holocene (ca. 2.5-0 kyr BP), an extended sea-ice cover even hindered primary production by limiting light penetration and nutrient release from sea-ice melt (Fig. 7c, e).

Reconstruction of marginal ice zones (MIZs).
Abundant phytoplankton production is commonly found in MIZs. In these regions, melting sea ice releases nutrients and less-saline waters into the surface layer, forming a stable and nutrient-rich environment above the pycnocline, i.e., favorable conditions for algae blooms [48][49][50][51] .
Dinosterol and HBI-III (Z-isomer) (further referred to as "HBI-III"; see "Methods") have both been proposed to indicate pelagic algal production. Studies of surface sediments from different polar regions show significant enhancement of HBI-III in the MIZs 52-55 , while no significant difference in abundances or distributions of phytoplankton sterols was found between the MIZs and open waters 53,54 . Thus, the HBI-III has been proposed as an MIZ indicator by these authors. Cross-plots of IP 25 vs. HBI-III show a better correlation than IP 25 vs. dinosterol (Supplementary Fig. 1c, d), supporting that HBI-III is closely associated with sea-ice environments. Hence, HBI-III is more or less absent when sea-ice cover is strongly reduced or absent (i.e., between~9 and 14 kyr BP; Fig. 7b, c, e). The IP 25 and HBI-III accumulation rates peaked at ca. 5.6 and 3.5 kyr BP (Fig. 7d, e), probably  representing two periods of MIZ situations (see details in "Discussion").
Deglacial-Holocene changes in freshwater discharge. Synchronous changes (Fig. 6d, e) as well as the positive correlation ( Supplementary Fig. 1a) between terrigenous biomarkers (campesterol and ß-sitosterol) and brassicasterol suggest that brassicasterol in our record has a mixed source from freshwater and marine algae. During the last deglaciation and the early Holocene (~14-8 kyr BP), very negative δ 13 C org values of −28 to −27‰ (Fig. 6f), relatively high terrigenous biomarker accumulation rates of 5-20 µg cm −2 kyr −1 (Fig. 7f), and high F C32 1,15 values up to 0.5 ( Fig. 7g) all point to high terrigenous input. Relatively warm and wet environments promoted terrigenous matter transport, and low sea-level conditions resulted in a closer distance of the core site to the river mouth receiving more terrigenous matter. In the mid-late Holocene (8-0 kyr BP), δ 13 C org values increased (>−27‰) and terrigenous biomarker accumulation rates decreased to 1 µg cm −2 kyr −1 (Figs. 6f and 7f). Similarly, F C32 1,15 values decreased to 0.3 (Fig. 7g), reflecting a descending river discharge and an increasing distance of the core site to the river mouth (cf. Ref. 56 ).    Sea surface temperature variations. During the deglaciation (~14−11.7 kyr BP), SST values were variable and higher in comparison to the mid-late Holocene SSTs (Fig. 7a), indicating relatively warm and open-water conditions. Open-water conditions with less sea ice allowed sufficient heat exchange between ocean and atmosphere. Hence the SST variations followed the global climate history recorded in the Greenland ice core 57 (Fig. 7a, j), with typical higher RI-OH′-SST values in the B/A interstadial and lower SST values during the YD stadial (Fig. 7a).
In the early Holocene (11.7-8 kyr BP), RI-OH′-SST values rose again (Fig. 7a), most likely attributed to the summer insolation maximum. Additionally, the total inundation of the Bering Strait occurred at ca. 11 kyr BP and the Pacific Water inflow may have contributed warm water to this site 45 . In the mid-late Holocene (8-0 kyr BP), RI-OH′-SST values were quite low and stable (Fig. 7a), reflecting expanded/seasonal sea-ice conditions in the mid-late Holocene (Fig. 7c).
Major sediment flux events at 13 and 11 kyr BP. Two distinct peaks at ca. 13 and 11 kyr BP are found in the accumulation rates of terrigenous biomarkers, dinosterol, and bulk sediment (Fig. 7e, f and Supplementary Fig. 2). The first peak may reflect the YD flood event (cf. Ref. 19 ) characterized by large amounts of terrigenous material and nutrients into the Beaufort Sea. Accumulation rates of terrigenous biomarkers ß-sitosterol and campesterol as well as b-GDGTs started to increase at 12.83 ± 0.15 kyr BP (Fig. 7f, h), close to the onset of Beaufort Sea freshening (12.94 ± 0.15 kyr BP) estimated by Keigwin et al. 19,58 . The YD flood transported large amounts of heterogeneous sediments into the Beaufort Sea and resulted in increases in magnetic susceptibility, density, and (probably detrital 59,60 ) CaCO 3 contents (Fig. 5b-d).
These sediments were most probably non-biogenic/organic and thus have caused a significant decrease of TOC content (Figs. 5e and 6g and Supplementary Fig. 2). Such changes in physical and chemical properties occurred slightly before the YD flood, probably attributed to culminated meltwater production in the northwestern outlet area at the end of B/A. In the earliest Holocene, an even stronger peak in terrigenous input occurred (Fig. 7f) at ca. 11 kyr BP. High heterogeneous sediment flux resulted in similar characteristics as caused by the YD event (Fig. 5b-e). Although the second event is characterized by high sediment flux and coincides with the post-YD flood, its origin could not be unambiguously proven yet. The second event may have partially different sediment sources from the YD event (see "Discussion").

Discussion
Core site ARA04C/37 is particularly sensitive to changes in paleo sea-ice cover, because it is located in the area of modern seasonal sea-ice cover (Fig. 2). In contrast to a foraminifera-based study on the nearby Core 750 that proposed a permanent ice cover during the last deglaciation (14-11.5 kyr BP) 61 , the records of PIP 25 and RI-OH′-SST in Core ARA04C/37 clearly point to open-water conditions that allowed heat exchange with the atmosphere (Fig. 7a, c). Additionally, the minor differences of 14 C ages between planktic and benthic foraminifera (Supplementary Table 2) indicate that the surface and bottom water masses were relatively well ventilated, further corroborating the scenario of ice-free conditions during the last deglaciation. During the early Holocene, on the other hand, both the foraminifera study 61 and our biomarker records (Fig. 7c) display more open-water conditions in comparison to the modern Beaufort Sea.
The Beaufort Sea itself is an oligotrophic system 6 and the nutrient-rich Pacific Water inflow is an important nutrient source stimulating primary productivity. Studies (foraminifera,  j Summer insolation (gray) 47 and δ 18 O values from NGRIP ice core (red) 57 . Triangles refer to AMS 14 C dates used in the age-depth model (see Fig. 5).
geochemical, and geophysical data) from sediment cores suggest an initial opening of Bering Strait at about 11.5 kyr BP 46,62,63 , whereas other evidence based on marine species dispersal suggests an earlier connection at about 13.3 kyr BP 64-66 . These contradictory evidences are recently reconciled by a gravitationally self-consistent sea-level simulation 45 . The authors predict the first opening of Bering Strait with shallow inundation at ca. 13 kyr BP. Substantial melting of the Cordilleran Ice Sheet and western Laurentide Ice Sheet (LIS) may have resulted in a local sea-level stillstand during 13-11.5 kyr BP, and then the total inundation of the Bering Strait occurred until ca. 11.5 kyr BP 45 . Primary production in our records remained relatively low during the deglacial to Holocene transition and generally increased from the early Holocene, except for two rapid increases triggered by high terrigenous input (Fig. 7e, f). It implies that even if the Bering Strait opened at 13 kyr BP, probably because of its shallow inundation, the Pacific Water influence on primary production in the southern Beaufort Sea was small. Upon the total inundation of the Bering Strait at the early Holocene, primary production started to increase, reflecting more influence of Pacific Water inflow.
HBI-III has been proposed as an indicator for MIZ [52][53][54][55] . Among these studies, a study of surface sediments from the Barents Sea found that HBI-III maxima are correlated with the winter ice edge and thus proposed HBI-III as an indicator for winter MIZ 53 . This is challenged by a study from the East Greenland shelf where an enhancement of HBI-III was observed near the mid-July ice edge in an East Greenland fjord 55 . The difference in seasonality is controlled by the general sea-ice situations in both areas. Major parts of the Barents Sea are more or less ice free during late spring, summer, and autumn, and thus the ice edge is restricted to the cold/winter season. Along the East Greenland continental margin controlled by the East Greenland Current, on the other hand, sea ice is much more extended throughout the year and thus, the ice edge is more restricted to the summer season.
In our case, HBI-III production was enhanced twice during the Holocene (Fig. 7e). In the early Holocene, i.e., prior to the first HBI-III peak, it was nearly ice free throughout the year, reflected in very low PIP 25 values < 0.2 and variable SSTs (Fig. 7a, c). Then, a general trend toward an increasing sea-ice cover is coinciding with weakened insolation (Fig. 7c, j). Due to sea-ice expansion, the core site was first proximal to an MIZ situation during cold seasons, when HBI-III peaked at ca. 5.6 kyr BP (Fig. 7e). The interval between the two HBI-III peaks (ca. 5-4 kyr BP) characterized by a clear decrease in HBI-III probably represents seasonal sea-ice conditions. As sea ice further expanded, the second period of MIZ may have existed during warm seasons at ca. 3.5 kyr BP, shown by peaked HBI-III (Fig. 7e). During the last 2.5 kyrs, reduced HBI-III and dinosterol production (Figs. 6b, c and 7e) implies a more extended ice cover. Sea-ice variability is illustrated in Fig. 8.
A novel proxy HBI TR 25 shows strong associations with the satellite-derived spring chlorophyll a concentration in the Barents Sea, and thus has been proposed as a spring phytoplankton bloom proxy 67 . Although the proxy is restricted to a regional area, the HBI TR 25 determined in Core ARA04C/37 seems to support our reconstruction of the MIZ history and by this the applicability of the proxy. Higher values of HBI TR 25 may point to enhanced occurrence of spring blooms between 8 and 5 kyr BP, indicating that the core site was relatively ice free in spring. Hence, a winter MIZ may have existed at certain times within this period. On the other hand, lower values of HBI TR 25 may suggest reduced numbers of spring blooms from 5.6 kyr BP to present. One of the explanations is the existence of a spring sea-ice cover, which may have limited the spring blooms. This coincides with the expanding sea-ice cover and longer sea-ice duration during this time interval. To some extent, a summer MIZ may have existed occasionally within this period. These scenarios are in agreement with the HBI-III record, suggesting that the core site was probably close to a winter MIZ at ca. 5.6 kyr BP and to a summer MIZ at ca.

kyr BP.
The changes in RI-OH′-SST values seem to correlate with changes in sea-ice concentrations. That means quite high but variable temperatures occurred during dominantly ice-free periods, while low but stable temperatures during ice-covered periods (Fig. 7a, c), suggesting the proxy's potential as a tool for temperature reconstruction in the Arctic Ocean. To our knowledge, the RI-OH′-SST proxy was first studied on surface sediments from the Chinese coastal seas and the Yangtze Estuary, and then a global surface sediment compilation was assessed 32 . Samples from northern high latitudes are mainly collected from the Fram Strait 68 . In these samples, residual SSTs derived from RI-OH′-SST minus remote sensing SST are relatively small but generally above zero. Therefore, the absolute values should still be interpreted with caution. Nevertheless, based on the close association between the RI-OH′-SSTs and the sea-ice concentrations, we propose that the proxy most likely reflects winter/spring SSTs in the southern Beaufort Sea. Low and stable SSTs during the seasonal sea-ice periods (mid-late Holocene) (Fig. 7a, c) suggest that RI-OH′ may record SST in the ice-covered seasons, likely in winter. As insolation declined in the mid-late Holocene, icecovered seasons might have prolonged and persisted into early spring.
Superimposed on the long-term changes in sea ice, primary production, terrestrial input, MIZ and SST are two extremely high sediment flux events during the deglacial to Holocene transition (at ca. 13 and 11 kyr BP), which are related to the YD flood and the post-YD flood events.
The YD cooling event (12.9-11.7 kyr BP) was the longest interruption to the gradual climate warming since the severe Last Glacial Maximum (Fig. 7j). Meltwater of LIS draining through the Arctic Ocean to the Nordic seas has been suggested to have profound impacts on slowing down the AMOC and triggering the YD cooling 17,22 . Glacial system modeling, and seismic reflection data as well as physical properties from the Beaufort slope, shows massive floods possibly via a northwestern outlet into the Arctic Ocean at the onset of the YD (Fig. 3) 22,69 . Recently, Keigwin et al. 19 provided initial freshening evidence of the southern Beaufort Sea. Core JPC15 recorded the freshening in a decrease of δ 18 O values of the planktic foraminifera Neogloboquadrina pachyderma by at least 1‰ at the start of the YD (Fig. 7i). In a pilot work, paleosalinity reconstruction from δ 2 H of palmitic acid in the same region (Core PS72/291-2; Fig. 1) shows a significant reduction of salinity directly prior to the YD 70 . Other evidence, e.g., decreased δ 18 O values of planktic foraminifera found in the Laptev Sea, Yermak Plateau, Amundsen Basin, and Mendeleev Ridge 71-74 , a prominent increase in sea-ice cover in the Laptev Sea continental margin 35 , and a pulse in detrital dolomitic limestones (Arctic-Canadian sediment source) in the central Arctic Ocean 59 , point to a YD-age Arctic freshening sourced from the LIS (Fig. 3) 75 .
Meltwater production from the northwestern sector of the LIS might be a component causing the Beaufort Sea freshening. Model simulation suggests that runoff to the Arctic Ocean from melting Keewatin Dome alone was sufficient to slow down the AMOC during YD 22 . Thus, it is possible that such melting culminated at the end of the B/A and strengthened freshening. This is supported by changes in physical properties of Core ARA04C/ 37 that slightly predated YD (Fig. 5). Another possible freshwater component is the outburst from proglacial/subglacial lakes, such as Lake Agassiz. Our multiple biomarker records also document this YD event and suggest the freshening was more related to a flood event as outlined in the following (Figs. 5 and 7).
B-GDGTs are synthesized by bacteria and can be found in numerous terrestrial settings, e.g., soils, lakes, and rivers [37][38][39][40] . As freshly deglaciated surfaces probably have negligible soil organic matter, the increases of b-GDGTs at the onset of the YD does likely not indicate input from soils. On the other hand, increase of meltwater in the rivers would not increase b-GDGTs production. Hence the main contribution of b-GDGTs increases seems not from the rivers. Another possible source could be from the melting ice sheet, i.e., b-GDGTs might have been entrained into the base of the ice sheet when it glaciated and were released when it deglaciated. However, such release and accumulation of b-GDGTs should be continuous processes. The sharp increase and peak accumulation of b-GDGTs at the start of the YD (Fig. 7h) point to lakes as the most probable source. In such lakes b-GDGTs could have been accumulated in a relatively stable environment for a long time and then been delivered to the Beaufort Sea within a short period. This interpretation is supported by the riverine input proxy F C32 1,15 (fractional abundance of C 32 1,15-diol). The C 32 1,15diol is found particularly abundant in freshwater systems 41 , and its production has been related to flow regimes 76 . A study from the Amazon River shows that the highest abundance of C 32 1,15diol exists in water bodies with low flow velocity and low turbidity, while the lowest abundance of C 32 1,15-diol is found in conditions with high turbidity and high flow velocity 76 . This suggests that a high-velocity flow regime like a flood producing higher turbidity would limit the C 32 1,15-diol production. Thus, the significant drop in F C32 1,15 values during the YD (Fig. 7g) characterizes a flood event of high erosive power and high flow velocity in turbid waters.
Assuming that the high-resolution terrestrial records of Core ARA04C/37 represent the YD flood event, a duration was estimated based on records of terrigenous sterols and b-GDGTs (Figs. 6d and 7h). The peak discharge started at 12.78 ± 0.15 kyr BP and ended up at 12.63 ± 0.13 kyr BP, resulting in a flood duration of about 150 years. Such estimate is close to the duration of 130 years estimated by Keigwin et al. 19 .
Although our terrestrial biomarker records in combination with previous paleosalinity reconstructions 19,70 indicate a strong YD flood, the flood water source still remains unclear. Coinciding with the YD event, a large lake-level drop of Lake Agassiz (with an estimated freshwater release of~17,000 km 3 ; Ref. 77 ) has been proposed as the freshwater source 23 . According to more recent studies, however, no clear field evidence (ice margins or shorelines) supports that Lake Agassiz drained to any northwestern outlet before 10.8-10.1 kyr BP 21,24 . It cannot be excluded that the northwestern outlet of Lake Agassiz opened at an earlier phase and its evidence has been removed by a readvance of LIS. Otherwise the flood may have other water sources, e.g., outburst from further northern proglacial/subglacial lakes. If so, the freshwater discharge from these lakes alone was too small to trigger the YD cooling. Besides freshening in the Beaufort Sea, input of freshwater during the YD was also found in the Northeast Pacific sea (Ref. 78 and references therein). A recent model simulation suggests that if the Bering Strait partially opened at the start of YD 45 , some of the freshwater from the Northeast Pacific could be transported through the Arctic Ocean to the Nordic seas, contributing to a collapse of the AMOC 78 .
In the earliest Holocene, following the YD flood, a post-YD flood at ca. 11.3 kyr BP has been identified by field data from the Fort McMurray area, and the freshwater source is suggested to derive from proglacial lakes McMurray, Meadow, and Churchill 24 . Evidence in the Fort McMurray region is in line with field data from the Mackenzie Delta supporting that a second highenergy fluvial episode occurred in the delta between 11.7 and 9.3 kyr BP 23 . More recently, seismic and geophysical data from the Beaufort margin also indicate that the second event probably initiated at ca. 11.3 kyr BP 69 . The timing of the second flood coincides with the PBO, thus it has been proposed that the second flood may have triggered the PBO 25,69 . As a large amount of freshwater from the Baltic Ice Lake entered the Arctic marginally predating the PBO 79 , Klotsko et al. 69 regarded the second flood injection as a tipping point for sufficient weakening of the AMOC during the PBO.
Despite the field data indicate a second flood, direct evidence for freshwater input is still missing from the marine records. The δ 18 O record of Core JPC15 shows minimum values only during the YD flood, whereas values were even higher than the 2‰ baseline during the second event (Fig. 7i) 19 . A similar case was described for a piston core from the Mackenzie Trough, where a coarse-layered unit was dated to 11.5-11.3 kyr BP 80 , while the very positive δ 18 O values (~2.7-3.5‰) argued against the hypothesis of an outburst from proglacial lakes. So far, the only freshwater evidence is found in Core P189AR-P45 (see Fig. 1 for location), showing a drop of δ 18 O values at ca. 11.5 kyr BP 81 . However, since the age-depth model of this core is not well constrained 81 , the timing of the freshwater signal might be shifted to the YD period when improving the age model. In comparison to the YD flood, biomarker records also show slight differences during the second event (Fig. 7g, h). Moderate values of F C32 1,15 during the second event indicate relatively stable flow regimes and smaller discharge than the YD flood (Fig. 7g). Otherwise extremely high flow velocity would limit C 32 1,15-diol production. Accumulation of b-GDGTs was also much smaller than during the YD (Fig. 7h), and seems to have different sources. This evidence gives rise to the question whether the high sediment flux was caused by a flood or whether it was controlled by other processes.
In this context, one should keep in mind that the second event also coincided with global meltwater pulse 1B (MWP 1B). Recent studies in the Okhotsk Sea and Pacific Beringia provided direct evidence of coastal erosion induced by the rapid sea-level rise 82,83 . In these records, two distinct maxima of terrigenous material were observed within the rapid sea-level rise at ca. 14 and 11 kyr BP (MWP 1A and 1B). In the Kara Sea, high terrigenous sediment fluxes characterized by significantly increased C/N ratios and refractory organic matter were recorded at about 11 kyr BP and also related to increased coastal erosion due to the postglacial flooding of the shelf sea 84,85 . In the Labrador Sea, a high sedimentation layer containing increased detrital carbonates was dated to ca. 11.5-11.3 kyr BP 86 . Although the authors related it to "Heinrich event 0," almost simultaneous signals widely found in the Arctic may support the interpretation of sea-level rise induced coastal erosion. It is therefore possible that the high sediment flux during the second event was partially caused by coastal erosion and subsequent deposition of the eroded terrigenous material.
Although field data indicate a second flood during the early Holocene, marine evidence suggests that the high sediment flux could be caused by a meltwater flood and/or shelf flooding induced coastal erosion. Further work, e.g., comprehensive proxy reconstruction of freshwater discharge and paleosalinity (cf. Ref. 70 ) as well as carbon source constraint (cf. Ref. 82 ), is hence needed to distinguish different processes.

Methods
Chronology. Radiocarbon dates of Core ARA04C/37 and JPC15 are consistent during the last deglaciation (Fig. 4). A red vertical line is placed close to the onset of the magnetic susceptibility rise (Fig. 5a, b), and the almost same 14 C ages from both cores at the vertical line suggest a good age constraint of the onset of magnetic susceptibility change. 14 C dates from the two cores are highly consistent below 400 cmbsf (Fig. 4a). However, age offsets are enlarged from near 350 cmbsf toward the core top (Fig. 4a), indicating a weakened correlation between Core ARA04C/37 and JPC15. Therefore, we established the age-depth model for core ARA04C/37 by including six AMS 14 C dates from Core JPC15 below 400 cmbsf (Supplementary Table 1). Paired dating shows a mean difference of around 200 years between planktic and benthic foraminifera (Supplementary Table 2). Considerable ventilation allows mixed species (planktic and benthic) being used for age constraints. As discussed in detail by Keigwin et al. 19 , ΔR was defined as 200 ± 100 years in YD and 0 ± 100 in the Holocene and B/A in the age model. The age-depth model was developed using the "Bacon" software of Blaauw and Christen 87 . Bulk sediments of the top cm of the core were analyzed by gamma spectrometry (HPGe planar detector). Based on excess 210 Pb in the uppermost centimeters and detectable anthropogenic 137 Cs in the top sample, surface sediments were identified to be modern, therefore, the core top was fixed to 0 kyr BP. Between 450 and 500 cmbsf (unit 1), 14 C dates are constant (mean 14 C age of 11,160 years) and are indicative for event 1 (Fig. 4a). Although the onset of event 1 is well constrained, because there is no clear termination of event 1 shown in magnetic susceptibility, the range of unit 1 is less constrained. This might lead to an over-or underestimate of sedimentation rate. Similarly, despite that magnetic susceptibility in both cores shows generally synchronous changes between 375 and 420 cmbsf (Fig. 5a, b), enlarged age offsets indicate a weaker correlation between the two cores ( Fig. 4a and Supplementary Table 1). Therefore, we define the unit 2 (event 2) mainly based on four constant 14 C dates (mean 14 C age of 9750 years) between 295 and 375 cmbsf, with poor constraints in both onset and termination of the event. The second high sedimentation unit (295-375 cmbsf) may be thicker than shown by the four dates, therefore we kept the four dates for more precise age model constraints, even though the last date of the plateau is not strictly increasing (Fig. 4). A "boundary" function was applied in the two units, which may have experienced distinctively high sedimentation rates (cf. Ref. 88 ). All the 14 C ages were transformed into calendar ages following Reimer et al. 89 via the Marine13 curve.
Bulk parameters. For organic geochemical analyses, sediment samples were freeze-dried, grounded, and homogenized. TOC contents were determined by a Carbon-Sulfur Analyser (CS-125, Leco) after decarbonization with hydrochloric acids. Total carbon (TC) contents were analyzed by a Carbon-Nitrogen-Sulfur Analyser (Elementar III, Vario), and used for calculation of carbonate contents (CaCO 3 = (TC − TOC) × 8.333). Carbon isotope composition of organic matter (δ 13 C org ) was analyzed by means of a Thermo Delta V Isotope Ratio Mass Spectrometer connected to a Thermo Flash 2000 CHNS/OH Elemental Analyzer employing continuous flow, performed at the Korea Polar Research Institute. δ 13 C org values are given in per mil notation relative to Vienna Pee Dee Belemnite International Standard. Reference gases were calibrated relative to Indiana University Acetanilide#1, USGS40, Urea, and Thermo Soil standards. For a randomly selected set of samples, duplicate analyses were carried out to determine the analytical error. Analytical precision is within 0.2‰.
The tri-unsaturated HBI biomarkers (triene, Z-and E-isomer) are commonly found in marine sediment. Based on studies in the Barents Sea, HBI-III (Z) producers favor nutrient-rich and stratified upper water column at the ice edge 53 . For simplification, the term "HBI-III" is used in the main text, representing the "HBI-III (Z)." The relationship between the relative proportions of HBI-III (Z) and HBI-III (E) and spring chl a in surface sediment from the western Barents Sea suggests the potential of a novel HBI proxy (HBI TR 25 ) indicating the spring phytoplankton bloom (the following equation) 67 : For GDGT analyses, 5 g sediment was ultrasonically extracted with DCM: MeOH (2:1 v/v), and internal standard (C 46 -GDGT, 1 μg per sample) was added prior to analytical treatment. The GDGT fraction was separated from other fractions via open column chromatography and eluted with 5 ml DCM:MeOH (1:1 v/v), dried and re-dissolved in hexane:isopropanol (99:1 v/v), and then filtered via a polytetrafluoroethylene filter with a pore size of 0.45 μm. Compound identification and quantification were carried out by a high-performance liquid chromatography/ atmospheric pressure chemical ionization-MS according to the method described by Meyer et al. 83 RI À OH 0 ¼ 0:0382 SST þ 0:1: Accumulation rates. Accumulation rates of bulk sediment, bulk organic carbon, CaCO 3 , and biomarkers are calculated by using the following equations (cf. Ref. 94 ): BulkAR ¼ SR WBD À 1:026 PO ð Þ ; ð6Þ BulkAR = bulk sediment accumulation rate (g cm −2 kyr −1 ); SR = sedimentation rate (cm kyr −1 ); WBD = wet bulk density (g cm −3 ); PO = porosity (%); TOCAR = total organic carbon accumulation rate (g cm −2 kyr −1 ); CaCO 3 AR = carbonate accumulation rate (g cm −2 kyr −1 ); BM = biomarker concentration (μg g −1 Sed); and BMAR = biomarker accumulation rate (μg cm −2 kyr −1 ).

Data availability
The data sets generated during and analyzed during the current study are available in PANGAEA repository (https://doi.org/10.1594/PANGAEA.915048) 95 .