Collapse of Eurasian ice sheets 14,600 years ago was a major

Rapid sea-level rise caused by the collapse of large ice sheets is a global threat to human societies1. In 10 the last deglacial period, the rate of global sea-level rise peaked at more than 4 cm/yr during Meltwa11 ter Pulse 1a, which coincided with the abrupt Bølling warming event ~14,650 yr ago2–5. However, the 12 sources of the meltwater have proven elusive6,7, and the contribution from Eurasian ice sheets has 13 until now been considered negligible8–10. Here we show that marine-based sectors of the Eurasian 14 ice sheet complex collapsed at the Bølling transition and lost an ice volume of between 4.5 and 7.9 15 m sea level equivalents (95% quantiles) over 500 yr. During peak melting 14,650 14,310 yr ago, 16 Eurasian ice sheets lost between 3.3 and 6.7 m sea level equivalents (95% quantiles), thus contribut17 ing significantly to Meltwater Pulse 1a. A mean meltwater flux of 0.2 Sv over 300 yr was injected 18 into the Norwegian Sea and the Arctic Ocean during a time when proxy evidence suggests vigorous 19 Atlantic meridional overturning circulation11,12. Our reconstruction of the EIS deglaciation shows 20 that a marine-based ice sheet comparable in size to the West Antarctic ice sheet can collapse in as 21 little as 300-500 years. 22

Contour lines represent ice margins at different stages of the deglaciation. Thick lines represent ice margin positions at boundaries between the deglacial phases used in the Bayesian chronology (Supplementary Data Fig. 7 and 8 and Supplementary Data File). Black lines are the inferred ice margin following the late Heinrich Stadial 1 ice advance. Pink lines are the ice margins that followed the separation of the BSIS and FIS. Yellow lines mark ice margins when the BSIS are constrained on the archipelagos and shallow banks in the northern Barents sea. The median age of each margin is indicated. The accompanying transparent fields mark the geographic uncertainties associated with the respective ice margins. Thin lines mark the suggested ice sheet retreat pattern within each phase as synthesized from the literature listed in Methods. The black stippled line marks the separation between the FIS and the BSIS used in the area-volume calculation when they were confluent. Black filled circles mark sites used to constrain the Heinrich Stadial 1 extent of the ice sheet. The positions of the stratigraphic records and dates used to constrain the deglacial phases are marked with gray, pink, yellow and white filled circles. White diamonds mark the position of cores used to reconstruct the Norwegian Sea 14 C reservoir age. White lines indicate ice margins adopted from the Dated-1 reconstruction.   22,27,28 . B, Magnetic susceptibility from two Norwegian sea sediment cores ( Fig.  1), aligned with the speleothem δ 18 O record in (A) (Methods). C, Average δ 18 O record from Greenland summit ice cores (GISP2 and GRIP) on the GICC05 chronology 30 . D, Planktonic foraminifera δ 18 O (Neoglobigerina pachyderma sinistral) from three Norwegian Sea sediment cores [31][32][33] . E, Proxy records of ice rafted detritus from Norwegian Sea cores 31,33 . F, Compiled AMS 14 C ages from Norwegian Sea sediment cores (GS07-148-17GC, this study; GIK23074 31,34 ; MD95-2010 33 ; HM79-6 32 ). Horizontal error bars represent the 68.2% quantiles (equivalent to 1σ) of the GS07-148-17GC deposition model. Gray shading represents ±1σ of the Monte Carlo sampling of the probability density functions of both the stratigraphic and chronological core alignments and the 14 C uncertainty. G, Norwegian Sea 14 C reservoir age, R is calculated as the difference between the conventional 14 C ages (at the median age) and the IntCal13 atmospheric 14 C curve 29 . Vertical error bars are the root sum of squares of the 14 C uncertainties. The average global reservoir age represented by the Marine13 calibration curve 29 is plotted for reference. H, Reconstructed ice volume for the Eurasian Ice Sheet (EIS) complex expressed as m sea level equivalents (SLE; 25 yr running mean of median and 95% quantiles). FIS: Fennoscandian Ice Sheet; BSIS: Barents-Svalbard Ice Sheet. I: median rate of ice volume loss in cm SLE per yr and as meltwater flux (Sv) (colors as in (H)). and 14 kyr BP), and is comparable to the estimated contribution from the much larger North American ice 100 sheet (5-6 m SLE in ref. 39 , 6.4-9 m SLE (interpolated to 340 yr) in ref. 40 , and 4-7 m SLE in ref. 10 ). Al-at this locality. For the high-end local sea level rise estimate of 17.3 m at the Sunda shelf 6 , our mass loss estimates correspond to 20-40% of the local sea level rise. A more accurate estimate of the eustatic sea-117 level contribution from the EIS collapse will require additional constraints on the effect of glacio-isostasy 118 and ice volume below flotation. Nevertheless, our findings provide strong empirical evidence that the EIS 119 was a major source of the MWP-1a. Combined with recent estimates for the North American Ice Sheet 120 MWP-1a contribution 10,40 our EIS mass loss estimates are sufficient for explaining the far-field relative sea 121 level observations without a major Antarctic contribution, consistent with the lack of field evidence for a 122 large retreat of the Antarctic Ice Sheet 45 . 123 Our new account of the EIS collapse is an important step towards solving the mysteries of the Bølling 124 event and the MWP-1a, which also raises a number of research questions pertinent to climate change sce-125 narios for the near future.

126
(1) What triggered the collapse of the marine-based EIS? In addition to the abrupt atmospheric and 127 surface ocean warming at the Bølling transition 32,46,47 , proxy records from core JM02-460 suggest a marked 128 subsurface warming on the Barents Sea continental shelf during the late Heinrich Stadial 1 48 , close to the 129 inferred ice sheet grounding line (Fig. 1). A vast ice-ocean interface rendered marine-based EIS sectors 130 potentially very sensitive to subsurface warming and melting at the grounding line, which is considered to 131 be one of the main drivers of current 49,50 and past 51 mass loss from the Antarctic ice sheets.

132
(2) Which mechanisms drove the rapid EIS retreat? In addition to surface melting and the likely in-133 volvement of mass-balance/elevation feedback 39 , continuity between subglacially carved lineations and 134 iceberg ploughmarks in the Bear Island Trough suggests calving of deep-keeled icebergs at the ice front 52 . 135 These findings are consistent with the operation of the marine ice cliff instability mechanism (MICI) 53,54 136 during the rapid ice sheet retreat. The current water depth in the SW Barents Sea is 400-500 m, less than 137 the~800 m thought to be required by MICI 53 . Isostatic depression by ice sheet loading 55 , however, may 138 have lowered the bed sufficiently for this mechanism to operate. Alternatively, the MICI may operate at 139 shallower depths than currently parameterized in models. Although past Antarctic deglaciation events can 140 be explained without invoking this specific mechanism 56 , the MICI is featured in the model yielding the 141 high-end future rate of ice loss from the Antarctic Ice Sheet 18 .

142
(3) What was the impact of EIS meltwater on ocean circulation? We estimate that a meltwater flux of 143 0.2 Sv over 300 yr was injected into the Norwegian Sea and the Arctic Ocean during the early Bølling, a 144 time period when proxy evidence suggests vigorous Atlantic meridional overturning circulation 11,12,57 . This 145 result implies that the relationship between freshwater injection and North Atlantic deep water formation is 146 not clear-cut, and highlights the need to resolve meltwater routing 58 .

147
Our reconstruction of the EIS deglaciation shows that an ice sheet comparable in size to the West 148 Antarctic ice sheet can collapse in as little as 300-500 years. Ice sheet models used to predict the future of 149 marine-based Antarctic ice sheets differ markedly in their predicted rates of ice loss and in the mechanisms 150 involved 17,18 . We provide new empirical constraints that raise the prospect of using the marine-based EIS 151 collapse as a benchmark for validating such ice sheet models and ultimately improve projections of future 152 sea-level rise. The estimated rates of ice loss from the EIS during the early Bølling (~1.6 cm SLE yr −1 153 averaged over 300 yr, peaking at~2.2 cm SLE yr −1 ) are comparable to high-end values of mass loss 154 projected for the West Antarctic ice sheet in the next centuries.

156
Temporal evolution of the marine radiocarbon reservoir age (R) δ 18 O and δ 13 C of N. pachyderma sinistral ( Supplementary Fig. 4). The alignment to the GS07-148-17GC depth scale was performed with the Oxcal v4.3.2 software 62 , using the P_Sequence sediment deposition 172 model 63 and the variable k option 64 . We assume an uncertainty of ± 2 cm (1σ) for each tie-point.

173
Absolute age control of the core records including 14 C was obtained by event-stratigraphic correlation 174 with the U/Th dated H82 speleothem δ 18 O record from Hulu Cave, China 27 and isotope records from 175 Greenland Summit ice cores 30 ( Supplementary Fig. 1). The rationale for this correlation rests on the close 176 relationship between Greenland temperatures, North Atlantic Ocean temperature and circulation, and the 177 Asian Monsoon on decadal to millennial time scales 22-25 . 178 For the correlation we used the MS record of core GS07-148-17GC determined in 2 mm steps by a

179
Geotek TM multi sensor core logger and a Barlington2 point sensor. MS in Norwegian Sea sediments is 180 considered to be a proxy for the strength of the warm Atlantic Water inflow over the basaltic Iceland Scot-181 land Ridge through ocean current erosion and transport of magnetic mineral grains that are subsequently 182 deposited in the S-Norwegian sea; the Atlantic water inflow is in turn tightly linked to the general North 183 Atlantic climate, including Greenland temperatures 33,65-67 . 184 We used the Hulu cave speleothem H82 δ 18 O record as the Norwegian Sea MS correlation target be-185 cause of its high temporal resolution, and because it contains high-amplitude signals that covary with the 186 MS record. This covariance has been attributed to fast atmospheric teleconnections between ocean circu-187 lation in the North Atlantic and regional Asian monsoon intensity and isotopic fractionation captured in 188 the speleothem δ 18 O 23,68 . The covariation between Greenland ice core δ 18 O and Norwegian sea MS is 189 more subdued, especially during Heinrich Stadial 1 (HS1), which has been attributed to a diminishing ef-

210
To assess the sensitivity of our results to the reconstructed chronology, we explored an alternative depo- independently constrained by paired marine and terrestrial 14 C dates 75 . We then used the Marine13 cali-217 bration curve 29 with a ∆R of 100 ± 50 yr, and the same depositional model as in our preferred chronology, 218 invoking the default general outlier model 76 . Due to a lack of pre-Bølling age constraints, this alterna-219 tive chronology expectedly shows much greater pre-Bølling age uncertainty than our preferred chronology.  Fig. 2). Notably, the alternative chronology yields a drop in 14 C age at the Bølling transition that is 222 steeper than in our preferred chronology, implying an even more abrupt EIS collapse. Hence, we conclude 223 that the inferred drop in R at the Bølling transition is unlikely to be an artefact of the age model, and that 224 our estimates are conservative in terms of the rate of EIS mass loss and its contribution to the MWP-1a. can be used directly in radiocarbon calibration software.

236
Our R record are consistent with R values previously reported from the North Atlantic and the Norwe-237 gian Sea and coast 34,75,[77][78][79] . Although a different approach was used to constrain the calender ages of core 238 GIK23074 34 , we arrive at similar reservoir ages.

239
Tephrochonology 240 Tephra shards were quantified in the >100 µm grain fraction in~20 cm interval of core GS07-148-17GC 241 corresponding to the Younger Dryas chronozone. This interval was chosen with the aim of finding the 242 Vedde Ash tephra that is a key chronostratigraphic marker horizon in the North Atlantic region, and is also 243 found in the Greenland Ice cores 30 and several of the Norwegian Sea cores used in this study 32,33 . Based on 244 their colour and morphological character, tephra particles were grouped into a transparent-white rhyolitic 245 type of tephra and a brown basaltic type of tephra. The total count from each of these tephra types was performed directly on the shards and without any leveling or polishing the beam will hit the surface from 255 different angles. This resulted in that the counting rate of the different elements becomes slightly more 256 scattered than during analysis on a polished thin section. The major element composition is, however, 257 consistent with published major element data from the Vedde Ash ( Supplementary Fig. 3).

258
Ice sheet margin reconstructions 259 We reconstructed the deglaciation of the EIS complex in a Bayesian chronological framework using Oxcal 260 4.2.4 62-64,76 . The prior model was constructed using available chronological, stratigraphical and morpho-261 logical data that were aggregated, independently for the BSIS and the FIS, into a sequence of phases with 262 known relative ages. A phase in this context refers to the retreat (or advance) of the ice sheet in a specific 263 area. 264 We grouped the deglaciation of the FIS ice sheet into two phases: (i) late HS1 advance and (ii) deglacia-265 tion on the continental shelf and outer coasts. Following the deglaciation of the continental shelf, we use the 266 ages and ice sheet geometries provided by the Dated-1 reconstruction 19 in the 14-10 ka interval, as these 267 are predominantly based on terrestrial dates not affected by our recalibration of the marine 14 C dates. The to the Svalbard archipelago. An early deglacial phase was added before the late HS1 advance, without 277 assigning ice sheet margins. At 12-10 ka we used the Dated-1 19 BSIS ice sheet geometries. 278 We adapt a previously proposed ice sheet retreat pattern for the southern Barents Sea, suggesting

290
To account for possible deviations in R from the reconstructed Norwegian Sea 14 C and Marine13, we 291 add a ∆R of 0 ± 50 14 C years (1 σ) to each marine radiocarbon age determination. To calibrate marine 292 conventional 14 C ages younger than 11800 14 C years, we use the Marine13 curve 29 , terrestrial dates are 293 calibrated with the IntCal13 29 .

294
For each phase of the deglaciation we outlined a succession of ice margins ( Fig. 1) based on published 295 sediment core data, geomorphological interpretations and ice sheet reconstructions for the BSIS 19,21,48,42,47,72,73,[80][81][82]111,[123][124][125][126][127][128][129][130][131][132][133][134][135][136] . The available information is, however, too sparse to yield continuous 297 time-synchronous margins and we stress that the reconstructed margins are intended to capture the general 298 pattern of retreat rather that to be accurate representation of the ice sheet at a specific time. To account for 299 uncertainty in the ice sheet geometry, we follow the approach of 19 and construct accompanying maximum 300 and minimum margins (Fig. 1). These are treated as the 95% quantiles. For margins derived from the 301 Dated-1 reconstruction, we use the their max and min margins 19 .

302
Ice sheet volume estimates 303 We converted the reconstructed ice sheet areas to volumes using the approximation proposed by Paterson 37 : where V is volume and S is area. Paterson's formula was determined empirically 305 by regression of measurements on six extant ice sheets and ice caps, the boundary conditions of which 306 are not directly comparable to those of the EIS. To assess the sensitivity of the volume estimates to the 307 regression assumptions, we also used the area-volume relationships from the output of a recent ice-sheet 308 model of the EIS 38 to convert the reconstructed areas volume ( Supplementary Fig. 9). Although the model-309 based regression yields an EIS volume that is 2.7 m SLE smaller than the Paterson approximation at the 310 start of the deglaciation, the difference in the estimated ice loss between 14.7 and 14.4 kyr BP is only~0.2 311 m SLE, which is negligible with respect to our conclusions ( Supplementary Fig. 9). Paleo-depths of the 312 continental shelves on which the EIS was grounded are obscured by an unknown amount of isostatic uplift 313 since deglaciation. Without correcting for ice volume below flotation through glacio-isostatic modelling, 314 which is outside the scope of our study, our estimated volumes cannot be interpreted as eustatic sea-level 315 change. For each ice sheet margin reconstruction and associated uncertainty estimates, we generated a 316 PDF of the volume estimate using Gaussian kernels. The volume-PDF and accompanying age-PDF of each 317 reconstructed ice sheet were resampled using a Monte Carlo technique detailed at https://github. compelled us to discard these 14 C dates from the 14 C reconstruction ( Supplementary Fig. 4).

327
To assess the potential impact of ambient biogenic sediment mixing on the observed decline in R at as possible, we let the hypothetical true decline in R be an instantaneous step change superimposed on the 336 overall linear trend in the observed 14 C record, and we assumed a constant mixed layer depth. Under this 337 scenario, if we invoked a drop in the modeled R record of~1,220 14 C yr from 14,600 to 14,550 calendar yr 338 BP and used a mixed layer depth of 6 cm, then the bioturbated 14 C ages simulated by TURBO2 provided a 339 reasonable fit to the observed 14 C record ( Supplementary Fig. 6). Hence, the effect of bioturbation would 340 be to temporally smear out a more abrupt event in the 14 C record. This smearing effect pushes the recali-341 brated 14 C ages for the start of the deglaciation backwards in time, and attenuates the estimated EIS melt 342 water flux. An upward bias towards older ages affects 14 C dates between~13,200 and 14,000 14 C yr BP in 343 particular, and is important to bear in mind if the 14 C record is to be used as a regional calibration curve.     The dark-and light-colored bands represent the respective 68.2% and 95.4% credible intervals of the model. The model is made by defining tie-points (diamonds and vertical dashes between (B) and (C)) between the magnetic susceptibility record of core GS07-148-17GC (C) and the δ 18 O record from Hulu cave (B) 27 . While the Bølling transition is associated with high sedimentation rates and deposition of plumites closer to the continental shelf edge and the ice sheet grounding line 82,100,109 , core GS07-148-17GC is located in a more distal setting where the direct influence from sediment-laden meltwater plumes is less likely. The interval with high sedimentation rates centered at about 17.5 kyr cal BP is related to the deposition of a plumite sourced from the Norwegian Channel Ice Stream [71][72][73]138,139 . Horizontal error bars in B-C represent the 1σ uncertainty of the Oxcal-generated agemodel for the respective records. (D), The average of the δ 18 O record from the Greenland summit ice cores (GISP2 and GRIP aligned on the GICC05 chronology 30 ), which is plotted for reference. The peak occurrence of the Vedde Ash in core GS07-148-17GC and the Greenland ice cores is indicated by the blue line. Note that the Vedde Ash has not been used to constrain the GS07-148-17GC chronology, yet the difference in the Vedde Ash ages is only 10 years. E, The distribution of tephra shards found in core GS07-148-17GC, including rhyolitic (black) and basaltic (red) shards. Arrows mark levels sampled for geochemical analyses of tephra shards ( Supplementary Fig. 3).  Supplementary Fig. 1 27 . C-D, the MS record of core GS07-148-17GC on the preferred (C, magenta) and alternative (D, blue) deposition model. The horizontal error bars in B, C and D represent the 1σ uncertainty of the Oxcal-generated deposition models for the respective records. E, the average of the δ 18 O record from the Greenland summit ice cores (GISP2 and GRIP aligned on the GICC05 chronology 30 ) plotted for reference. F, the 14 C ages of the Norwegian Sea compilation plotted both on our preferred chronology (magenta) and the alternative chronology (blue), the light pink field is the Norwegian Sea 14 C reconstruction.    Figure 4: Norwegian Sea data records plotted on GS07-148-17GC depth scale. A, Depth models of cores HM79-4, GIK23074-1 and MD95-2010 constructed using the P_Sequence option in OxCal 63 . Light-colored uncertainty envelopes represent the 95.4% quantiles, while darker colored represent the 68.2% quantiles of the depth model PDF. The models are made by defining tie-point between the cores and core GS07-148-17GC using the records of (B) δ 18 O 31-33 , (C) δ 13 C 31-33 , (D) IRD 31,33 , and (E) magnetic susceptibility 33 . F, Compiled AMS 14 C 31-34 . Circles mark the dates that are excluded from further analysis due to distortion of the core stratigraphy from deep burrows ( Supplementary Fig. 5). Horizontal error bars in B-F represent the 1σ uncertainty of the depth model for the respective cores. Supplementary Figure 6: The effect of bioturbation on the 14 C reconstruction at the Bølling transition. To assess the potential impact of bioturbation, we used the TURBO2 model 137 (Methods). As input we used 1,024 simulated abundance vectors (gray; top panel) generated as normally distributed random values centered on the best-fit linear trend and with the standard deviation of the observed abundance of foraminifera in core MD95-2210 33 (top panel). If we assume a constant mixed layer depth of 6 cm, then the observed change in 14 C age can be reproduced with reasonable accuracy in TURBO2 by invoking a hypothetical true 14 C age with an abrupt step change 14.56 kyr ago (lower panel). This result is not an attempt to infer the true 14 C age history, but rather to demonstrate that the effect of bioturbation would be to smear out the true event. As a consequence, our reconstruction is likely to overestimate the time scale of the EIS collapse and underestimate its contribution to the global MWP-1a.  Figure 7: Bayesian deglacial chronology of the Norwegian continental shelf. As prior information, all radiocarbon dates or probability density functions of sediment unit boundaries are grouped into phases according to geographical and/or stratigraphical context. A phase in this context refers to a retreat (or advance) of the ice sheet in a specific area. The phases are ordered in a sequence following the relative chronological order. The PDF's of unmodeled conventional 14 C dates are calibrated using the new Norwegian Sea 14 C age reconstruction (Fig. 2) and is shown as light gray. Dark gray mark the modeled posteriori PDF of the same dates. Red PDF's show the posteriori age probabilities of undated events that corresponds to reconstructed ice margins depicted in Fig. (1).

Sequence boundary
Phase grouping dates from coastal SE Svalbard and south of the the Hinlopen Strait As prior information, all radiocarbon dates or probability density functions of sediment unit boundaries are grouped into phases according to geographical and/or stratigraphical context. A phase in this context refers to a retreat (or advance) of the ice sheet in a specific area. The phases are ordered in a sequence following the relative chronological order. The PDF's of unmodeled conventional 14 C dates are calibrated using the new Norwegian Sea 14 C age reconstruction (Fig. 2) and is shown as light gray. Dark gray mark the modeled posteriori PDF of the same dates. Red PDF's show the posteriori age probabilities of undated events that corresponds to reconstructed ice margins depicted in Fig. (1).