Introduction

The global threat of rising sea levels has triggered a large and diverse research effort to improve the understanding of the long-term stability and transient dynamics of ice sheets and their capacity to deliver vast fluxes of meltwater to the oceans. The knowledge on the non-linear interplay between oceanic and atmospheric circulation and ice sheets is central to climate system dynamics and in identifying key thresholds, tipping points and abrupt shifts. An issue repeatedly highlighted in this research has been to investigate how factors, such as air and ocean temperatures, mass balance and subglacial conditions, ice sheet rheology, and calving dynamics control the stability of the marine sectors of the Greenland and Antarctic ice sheets; a knowledge which is fundamental for reliable prediction of future sea level rise1. A major challenge in these endeavours is the relatively short timeframes of direct observations documenting ice sheet retreats and associated environmental forcing across glaciated regions. Investigation of former ice sheets offers an alternative perspective and can contribute with important knowledge on ice sheet behaviour by identifying the key controls on the pace and pattern of deglaciation over centennial to millennial timescales. A limitation in this respect is that while the retreat rates and patterns of the Northern Hemisphere former ice sheets (including the Eurasian and Laurentide ice sheets) have been fairly well constrained for their terrestrial parts, the marine parts have been less documented2,3,4.

In this study, we focus on the last deglaciation of the marine sectors of the Eurasian Ice Sheet (EIS) and its disintegration into three separate ice masses: the British-Irish (BIIS), the Fennoscandian (FIS) and the Kara-Barents Sea-Svalbard Ice Sheets (KBSIS)(Figs. 1–3). Previous reconstructions of the last deglaciation of the marine sectors of EIS have mostly been based on radiocarbon dates from sediment cores and extrapolation from adjacent land areas2,5,6. In this study we are taking advantage of new knowledge on the last deglaciation of the northwestern European shelf areas obtained the last decade, especially through large research efforts such as the BritIce-Chrono and GLANAM (Glaciated North Atlantic Margins) programmes7,8. The compilation presented (Fig. 2) includes results from recent papers with a large number of new dates and information of the changing geometry of the EIS, including the location of ice streams, deduced from ice marginal and ice directional features9,10,11,12,13,14. Also, the increasing availability of high-resolution seismic data has greatly helped both in interpretation of sediments (tills vs. glacial marine) and in providing a better understanding of the deglaciation patterns. From a compilation and recalibration of existing radiocarbon dates and of new and published mapping of glacial landforms (see Methods section, Fig. 1 and Supplementary Note 1), we present a refined reconstruction of the EIS disintegration in 1 ka time slices from the North Sea to the Kara Sea from 20 to 14 ka (Fig. 2). How the new reconstructions of the ice sheet extent for the marine sectors deviate from those of DATED-12, the last major time slice reconstruction carried out for the entire EIS, is presented in Supplementary Fig. 1. Our reconstructions are melded with those of Hughes et al.2 for the western seaboard of the Britain Isles and with those of Stroeven et al.3 across the terrestrial sectors of Fennoscandia. These studies are also used to extend the series to cover 25 to 10 ka. Following this, a steady-state ice sheet model is applied to each 1 ka time slice from 25 to 10 ka to calculate net ice-volume changes through the glaciation. Furthermore, we have compiled summer sea surface temperature (SSST) data spanning 25 to 10 ka from two well-dated cores situated under the trajectory of present inflowing Atlantic water in the eastern Nordic Seas (Figs. 2 and 4e). Using these refined empirical and modelled ice sheet reconstructions in combination with oceanic temperature records, we then discuss the key factors controlling the retreat dynamics of the marine sectors of the EIS, including the key role of the two major EIS ice streams, the Norwegian Channel Ice Stream (NCIS) and the Bear Island Ice Stream (BeIS) (Figs. 2 and 3 and Supplementary Figs. 2–4).

Fig. 1: Maps of investigated areas and mapped glacial landforms.
figure 1

a Overview map (bathymetry from www.gebco.net) with Last Glacial Maximum (20 ka) extent (black line) of the Eurasian Ice Sheet and boundaries (broken black lines) between the Kara-Barents Sea-Svalbard Ice Sheet (KBSIS), Fennoscandian Ice Sheet (FIS) and the British-Irish Ice Sheet (BIIS) indicated. Location of deep sea cores (LINK17 and JM373; white squares) and dated samples (Supplementary Note 1 and Supplementary Data 1–3) (black dots) are indicated. b, c Ice marginal features (red lines) and glacial lineations (black lines) mapped on the Mid Norwegian margin and compiled for the North Sea, largely based on refs. 9,10,12,22. d Ice marginal features (red lines) and glacial lineations (black lines) in the Barents Sea-Svalbard region, compiled using datasets from refs. 32,61,62,67,82,83. NC Norwegian Channel, FG Fladen Ground, DB Dogger Bank, SuT Suladjupet Trough, SkT Sklinnadjupet Trough, TT Trænadjupet Trough, BiT Bjørnøya Trough, ST Storfjorden Trough, StAT St Anna Trough, SI Shetland Islands, EA East Anglia.

Full size image
Fig. 2: Reconstructed EIS ice margins for every 1 ka for the period 20-14 ka.
figure 2

Ice margins from the terrestrial areas and marine areas west of the UK are from Hughes et al.2 and terrestrial Fennoscandian region from Stroeven et al.3. Extent of ice dammed lake in the southern North Sea21 and location of the Norwegian Channel (NCIS) and Bear Island (BeIS) ice streams are indicated. Bathymetry from www.gebco.net.

Full size image
Fig. 3: Modelled ice thickness and isostatically adjusted topography.
figure 3

Eurasian Ice Sheet configuration immediately before and after break-up of the major EIS saddles in the south at 19–18 ka (a, b) and in the north at 16–15 ka (c, d). Bathymetry and topography from www.gebco.net.

Full size image
Fig. 4: Time series of Eurasian Ice Sheet development and possible controlling factors for the period 25 to 10 ka.
figure 4

a Ice volume of the EIS (black line), model run with shear stress values by ±20% (broken black lines) and volume change in thousand km3 (broken brown line). b Ice volume of the FIS (green line), KBSIS (red line), and BIIS (yellow line) in km3. c Length of marine margin in km (colour coding as in a, b). d Summer insolation at 65°N in Wm-284 (orange line), EIS volume in thousand km3 (black line) and global sea level85 (blue line). e Summer sea surface temperature estimates based on planktonic foraminifera assemblages from two deep sea cores, Link 17 and JM373. f Modelling of eastern Nordic Seas ocean temperatures36. Grey columns indicate the millennium for break-up of the Eurasian Ice Sheet (EIS) into the Kara-Barents Sea Ice Sheet (KBSIS), Fennoscandian Ice Sheet (FIS) and British-Irish Ice Sheet (BIIS). Red column indicate the Bølling-Allerød (B/A) time period and Meltwater Pulse-1A (MWP-1A) are indicated with a blue box.

Full size image

The methods utilised and background data are further presented in Methods and in Supplementary Information.

Results

EIS deglaciation

From 25 to 20 ka, commonly referred to as the Last Glacial Maximum (LGM), the continental shelves of NW Europe were covered by a contiguous ice mass that spanned from offshore Ireland, across the North Sea, to the Kara Sea in the high Russian arctic. Except for a somewhat larger uncertainty in the northeastern sector, there is a general agreement on the maximum extent of the EIS during this period2. It is also generally accepted that the maximum extent of the EIS was asynchronous. For example, in northwestern Russia it has been suggested that the maximum was reached as late as c. 17 ka15, whereas across the Irish shelf, off the Shetland Islands and along the western Svalbard shelf there is evidence that the maximum ice extent close to the shelf edge was attained before c. 24 ka9,16,17,18. Around 19 ka, the western EIS margin had largely started to retreat from the areas close to the shelf edge (Fig. 2). The timing of this retreat has been determined from dates of ice-fed glacigenic debris flows (GDFs) found on the upper continental slope, indicative of shelf-edge glaciation19, and of deglacial sediments on the shelf. This first phase of retreat was most rapid within in the glacial troughs, which hosted fast flowing ice streams, located along the seaward margin from the North Sea to the Arctic Ocean.

For the North Sea, our reconstruction deviates both in the timing and pattern of deglaciation compared to Hughes et al.2 (Supplementary Figure 1). In the eastern North Sea the retreat of glacial ice from the shelf edge develops into a rapid retreat of the Norwegian Channel Ice Stream (NCIS) resulting in a deglaciation of the Norwegian Channel all the way to the Skagerrak/Swedish west coast by c. 17 ka20 (Fig. 2). At c. 18.7 ka the FIS and BIIS separated along a north-south line east of the Shetland Islands and west of the Norwegian Channel (Fig. 3), leading to drainage of an ice-dammed lake south of the Dogger Bank12,21. The drainage of the ice-dammed lake is evidenced by laminated plume deposits at the upper south-eastern Norwegian Sea continental slope dated to c. 18.7 ka, flood deposits and fluvial break-through channels in the Dogger Bank region21. For the western part of the North Sea we largely follow the reconstructions presented by Evans et al.10 and Roberts et al.22 in the southern part, by Sejrup et al.23 in the Fladen area and by Bradwell et al.24 around the Shetland Islands. The BIIS, including the North Sea Lobe along the eastern UK coast, retreated northwards from its maximum limit in East Anglia and the Dogger Bank area, and the Fladen Ground area and shelf off the Shetland Islands were mostly deglaciated between 16 and 17 ka9,10,22. Evans et al. (2021) suggested, with question marks, north-south trending ice margins of the BIIS in the central North Sea at 19 and 20 ka. However, these suggestions are not compatible with the dates from the eastern part of the basin and the relatively well-dated evidence for the separation of the BIIS and FIS in the central North Sea at 18.7 ka12,20,21.

Our reconstruction suggests that on the mid Norwegian margin grounded glacial ice remained at the shelf off Møre in the south, and west of Lofoten, in the north, as late as 16 ka (Fig. 1a and 2). In the areas between, the ice sheet had probably receded to a mid-shelf position at 16 ka, an exception being the Trænadjupet Trough, which was more or less ice free at this time (e.g., Laberg et al.25 and Nygård et al.26). At c.15 ka the ice margin on the mid Norwegian margin was close to the coastline.

Shelf-edge glaciation along the northern and western margins of the Barents Sea during maximum extent is widely evidenced27,28. In contrast, considerable uncertainty remains over the maximum extent, ice configuration and deglaciation pattern along the eastern margin of the KBSIS, reflecting sparse geological and chronological datasets in the Russian sector. One uncertainty is related to whether the ice sheet extended onto the Taymyr Peninsula (Fig. 1) during the LGM as it was tentatively suggested based on radiocarbon dates on molluscs within melt-out tills29,30,31. However, such a marginal lobe is unlikely, given the uncertainty of the age of these tills and the unlikely ice margin configuration such an extension implies, and given that Severnaya Zemlya is confidently asserted to be ice-free at this time31.

Our proposed deglaciation of the eastern KBSIS represents a major deviation from the DATED-1 time slices presented by Hughes et al.2. Whereas previously the Saint Anna Trough and Novaya Zemlya (Fig. 1 and Supplementary Fig. 1) were proposed ice free by 18 ka, we suggest that the main trunk of the Saint Anna Trough and the deepest parts of the neighbouring Franz Victoria and Kvitøya troughs were not ice-free until 16 ka (consistent with Hogan et al.32). Our reconstruction also indicates that partial deglaciation of Novaya Zemlya first commenced at 15 ka, with deglaciation of the northern island not completed until 14 ka, 4 ka later than the most credible estimate of Hughes et al.2. As Hughes et al.2 clearly acknowledge, the few available age constraints from Novaya Zemlya do not provide a consistent or sufficient account of the timing or pattern of deglaciation, with some suggesting a late33 and others an early deglaciation34. We note that recent transient ice-model simulations support a late deglaciation with ice retreating towards these topographic highs, rather than away from them35,36.

Along the western Barents Sea margin, retreat from the shelf edge is first noted in the Storfjorden Trough at 20 ka37, followed c. 1 ka later by a retreat along the rest of the western margin (Fig. 2). In the Bjørnøya Trough (Fig. 1), ice stream retreat is marked by a series of major grounding zone wedges. Only the outermost wedge is reliably dated and indicates a significant readvance at c. 17 ka38, whilst the morphology and spacing of other grounding zone wedges suggest large fluctuations in the retreat rate and multiple readvances/stillstands. At 15 ka the FIS and the KBSIS had separated in the Barents Sea, first through a corridor along the Norway-Kola Peninsula coast39 (Fig. 3). The separation of the KBSIS and FIS and subsequent collapse of the marine-based ice sheet seem to coincide with Melt Water Pulse (MWP)-1A40,41 and is reflected in sediment cores from the upper continental slope along the western parts of the Barents Sea (Fig. 4). Here evidenced by rapidly deposited laminated sediments/plume deposits topped by an IRD layer dated to between 15.1 and 14.3 ka28. By 14 ka most of the KBSIS is gone, but some glacial ice still remained across the archipelagos adjacent to the Arctic Ocean.

As seen from Figs. 1a and 2, the number of available dates and ice marginal features mapped are not evenly spread over the marine areas investigated. Therefore, our reconstructions of ice extent, for the time slices focused on, have larger uncertainties in some regions. Especially, the eastern parts of the Barents and Kara Seas, and also the south-eastern parts of the North Sea, are hampered by limited amount of data.

Ice sheet modelling

A steady-state ice sheet model (ICESHEET v1.042) assuming perfectly plastic flow43 is applied to our reconstructed time slices of the EIS from 25 to 10 ka (Supplementary Figs. 2 and 3 and Supplementary Data 4–6). Input boundary conditions include the bed topography, isostatically adjusted through the glaciation according to a first-pass ice-thickness solution, and basal shear stress values (Supplementary Fig. 5), parameterised on the basis of topography and sediment cover44. These values were perturbed ±20% to envelope a reasonable uncertainty range on potential ice volume changes between each time slice (see Methods section) (Fig. 4a, b).

At its peak volume, the EIS locked up a maximum of 7.9 (7.1–8.6) × 106 km3 of ice at 19 ka, equivalent to 20.04 (17.98–21.90) m of sea-level rise. Net ice losses increase steadily from this time slice—the period of maximum extent into eastern sectors—culminating at c. 15 ka through the collapse of the marine-based KBSIS and the retreat of the FIS from its marine sectors.

The partitioning of modelled ice volume reconstructions between the three major nucleation centres of the ice complex reveals a nuanced, asynchronous response of the EIS marine sectors during deglaciation. At 20 ka BP, approximately 40% of the total EIS is grounded in marine sectors, and diminishes at an accelerating rate—particularly after 16 ka - to only 14% by 14 ka (Supplementary Figs. 2 and 3). The increasing significance of marine sectors on deglaciation is demonstrated by the change in the length of the ice margin exposed to a marine setting (Fig. 4C). At 20 ka, ~38% of the total EIS margin is marine terminating, representing a total distance of 5055 km. By 16-17 ka, this proportion has risen to 47% (c. 7700 km). However, the timing of this destabilisation of marine sectors exhibits a strong latitudinal trend, with an increasing persistence of marine-based ice northwards (Fig. 4C). Between 20 and 14 ka the EIS lost a volume of 12.57 m (11.24–13.78 m) sea-level rise equivalent. Peak contributions to eustatic sea-level rise of 3.95 m occurs at c. 15 ka, partitioned roughly equal between Fennoscandian and Barents Sea sources. While around half of this contribution is from ice grounded in terrestrial sectors, the signal of the KBSIS collapse is clear, accounting for 87% (1.75 m) of the marine-based component (Supplementary Fig. 3). This peak follows the dominant influence of marine retreat across Fennoscandia and the British Isles by c. 3 ka, when these two ice sheets contributed to 56% (0.46 m) of the total marine-based sea-level rise contribution from the EIS (Supplementary Fig. 3).

Discussion

Generally, the overall ice loss of the EIS follows the summer insolation at 65°N (Fig. 4d), demonstrating its primary influence in triggering the onset of Northern Hemisphere deglaciation through increased summer ablation. The notable deviation of the ice-sheet loss curve from the insolation and sea level curves between 16 and 14 ka can largely be explained by the rapid loss of the mostly marine-based KBSIS. We further note that when most of the North Sea is deglaciated close to 17 ka, the Barents Sea is still largely glaciated (Figs. 2 and 3). This asynchronous retreat pattern of the EIS demonstrates a strong latitudinal difference in the ice sheets sensitivity to controlling factors.

Over the last decade numerous modelling studies have focussed on the interplay between ice sheet disintegration and stability of the Atlantic Meridional Oceanic Circulation (AMOC)45,46,47,48,49,50,51. It has been suggested that the combined effect of deglacial meltwater sourced from EIS, the Greenland Ice Sheet and the Laurentide ice sheets into the North Atlantic/Nordic Seas following the LGM led to a severe weakening of the AMOC which resulted in surface water cooling (shallower than 400 m) and warming of the deeper sub-surface water in the Nordic Seas45. Based on ocean and ice sheet modelling it has also been proposed36,49 that ice retreat in the Barents Sea between 20 and 15 ka could largely be explained by gradual warming of the sub-surface water in the eastern Nordic Seas (Fig. 4e). This gradual warming could result in increased sub-marine melting and grounding line instability, at a time when atmospheric temperature conditions were otherwise too cold in the Barents Sea region to impart significant surface melt. The modelled change in surface water in the eastern Nordic Seas is supported by planktonic foraminiferal evidence37,52 which reveal a general cooling in near-surface (c. 10–200 m) summer temperatures from the LGM and culminating close to the Bølling/Allerød transition at c. 15 ka (Fig. 4e). The fact that the EIS has its LGM maximum when the eastern Nordic Seas surface waters (Fig. 4e) were warmer than during most of the deglaciation suggests that the oceanic influence on atmospheric temperature over the ice sheet was not critical. However, the warm surface water during LGM must have contributed significantly to ice sheet growth through enhanced precipitation. Numerical modelling shows, especially in the north, that the decrease in precipitation following the decrease in surface water temperatures after LGM would contribute to a negative Surface Mass Balance (SMB) during the deglaciation phase53.

Our reconstructions and modelling show that the Bjørnøya Trough and the Norwegian Channel—host to the two largest ice streams of the EIS—are deep enough, when isotatically adjusted, to continuously expose the ice margin to deeper and warmer water throughout the deglaciation (Supplementary Figure 4). In addition, recent palaeo-evidence from benthic foraminiferal assemblages show sub-surface water warming along the western Barents Sea margin close to the BeIS during deglaciation37,54. In the North Sea the bulk of the deglaciation appears to take place well before (19–17 ka) the sub-surface water is considered to start to warm up (Fig. 3f). Planktic foraminiferal evidence from the Norwegian Channel55,56 and from the SE Nordic Seas (Fig. 3e) indicate that the near-surface water is still cooling during this time, presumably in response to AMOC slow-down subsequent to LGM deglaciation as suggested by climate models36. Model sensitivity experiments of the NCIS deglaciation also emphasise the significant vulnerability of this ice stream to melt-driven feedbacks that enhance surface lowering and instabilities associated with retreat along its retrograde bed57.

Based on the reasoning above we infer that the disintegration of the NCIS and the following deglaciation of the marine-based ice sheet covering the North Sea was not strongly influenced by sub-surface ocean warming as has been suggested to be the case for the retreat of the BeIS and the Barents Sea deglaciation. This is in line with evidence showing relatively cold North Atlantic surface waters during periods with large meltwater output from the southern sectors of the FIS58. Also, it indicates that the early melting of glacial ice in the North Sea could be attributed to increased surface melting and changes in SMB. Ice dynamics related to subglacial stress dependent on geomorphology, sediment cover and subglacial meltwater may have also played a role. The large number of channels in the North Sea, suggested to be the result of erosion by subglacial meltwater (tunnel valleys)59,60, as opposed to the sparsity of such features beyond the central sector of the Barents Sea, suggest different subglacial conditions between the two regions. The abundant evidence of subglacial water in the North Sea may also be an indication that ablation controlled by air temperatures are more important here than in the Barents Sea during the deglaciation.

The increase in meltwater delivery to the Nordic Seas between 16 and 14 ka, results largely from the rapid disintegration of the KBSIS (Figs. 2 and 4a). A possible extreme event shortly after 15 ka that could have contributed to the MWP-1A with as much as 3-6 m of sea level equivalents in <500 years has been proposed by Brendryen et al.41. Abundant meltwater features from the central Barents Sea associated with collapsing ice streams, including bedrock-incised tunnel valleys and subglacial lake basins61,62, infer that the KBSIS was becoming increasingly susceptible to surface melt processes through the latter stage of its collapse. The period with most meltwater delivery from the EIS takes place between 16 and 14 ka and predominantly associated with the collapse of the KBSIS (Fig. 4b). This coincides with a period of marked water column stratification in the western Barents Sea during HS1, characterised by a strong contrast between warm bottom water conditions and cold surface water63,64. In this period a peak in warm planktonic foraminiferal faunas seen in both the northern and southern core suggest a brief strengthening of the AMOC (Fig. 4e). This may imply that meltwater events sourced from the EIS did not have a profound and long-lasting impact on the North Atlantic oceanic circulation. In addition, modelling of the effects of different and rapid disintegration scenarios of the marine-based sectors of the EIS on ocean circulation suggest a similarly minor influence65. The period after 14 ka is characterised by a strongly reduced total meltwater delivery from the EIS, and is associated with observations indicative for a strengthened AMOC (Fig. 4e). This demonstrates that the AMOC strength in eastern Nordic Seas during deglaciation is strongly influenced by factors external to the regional domain, including meltwater delivery also from the Laurentide and Greenland ice sheets66.

Conclusions

On the basis of our reconstructed 1 ka time slices, it is apparent that the factors controlling the initiation of the last EIS deglaciation in marine sectors are complex and varied. The EIS collapse was not homogeneous, with the tipping points for the collapse characterised by regional oceanic and climatic differences. The North Sea deglaciation occurred relatively early in such a development and is therefore more likely to have been influenced by direct insolation and SMB changes than through the impact of ocean circulation changes. Across the Barents Sea, we find that a delayed heating of the subsurface waters strengthened by a post-LGM gradual decrease in precipitation related to the cooling of surface waters, coincides with the acceleration in retreat between 16 and 15 ka. A strong decrease in meltwater delivery once the majority of the marine sectors of the EIS had disintegrated at c. 15 ka could have provided an important contribution for AMOC strengthening. The interplay of factors controlling EIS deglaciation highlighted here, and the importance of the marine sectors, provides an important benchmark for fully coupled ice sheet modelling that must account for ice-ocean-climate feedbacks in a holistic approach. Yet, addressing this challenge will also require more accurately constrained empirical data. Lessons learned from the demise of Pleistocene ice sheets may be an imperative context for efforts in forecasting the fate of our still existing large ice sheets.

Material and methods

Mapping of landforms

Results from offshore mapping of glacial landforms along the NW European seaboard have become increasingly available over the last decade8,14. The mapping results presented here (Fig. 2) are largely based on the Olex bathymetric database (www.olex.no), but also other types of seabed imagery/bathymetric data of different resolution and quality, as well as shallow acoustics, have been utilised in the cited studies. In the present compilation we have focused on ice marginal features (including terminal moraines and grounding zone wedges) and ice directional features (features identified as glacial lineations) identified in the North Sea, on the Mid Norwegian margin, the Svalbard continental shelves and in the Barents Sea (Fig. 1). For the North Sea and Barents Sea the compilation is largely based on published data (e.g. Sejrup et al.12 and Winsborrow et al.67) whereas for the Mid Norwegian margin an extensive new mapping effort has been done for this study.

The ArcMAP (www.esri.com) and QGIS (www.qgis.org) softwares were used to generate the map figures (Figs. 1–3, Supplementary Figs. 1, 4 and 5)

Geochronology

In the present study we are presenting only radiocarbon dates as the reconstructions for the marine domain are mostly based on the interpretation of such dates. We acknowledge the promising work by Roberts et al.68 in utilising OSL dates in the southern North Sea and that these dates generally support the placing of the LGM in this region close to the Dogger Bank. However, we find that the large reported uncertainty and additional uncertainties related to water content, bleaching and context, do not warrant to put too much weight on these dates in constructing our 1 ka time slices from 20 to 14 ka. Available radiocarbon dates from the mapped regions have been compiled (Fig. 1a and Supplementary Data 1–3). This compilation is largely based on the database published by Hughes et al.2 and with dates published after 2015 added. Only dates on marine carbonate (molluscs and foraminifera) from marine sites are included. Dates considered to be of less value concerning the dating of the last deglaciation have been omitted. All the dates have been recalibrated utilising the Marine20 curve69 and the OxCal programme. A major problem in radiocarbon dating of marine carbonate is the uncertainty in Delta R (marine reservoir effect offset) which will change both with time and space. A number of papers have addressed this problem partly by the identification of terrestrial dated Icelandic tephras and also by comparing marine timeseries with well-dated records such as Greenland ice cores or Asian cave stratigraphies, assuming that pattern seen in timeseries of isotopes, sediment properties or fossil content can be directly correlated18,41,52,70,71. Because of the challenges with such correlations and the possibly of spatial and temporal variability of Delta R in the region, we have chosen the conservative approach and set Delta R to zero. In evaluating the weight put on individual dates we have also considered the stratigraphical context, possibility for mixing in foraminiferal assemblages and if the foraminiferal and mollusc dates can be assumed to date the sediments in which they are found. We acknowledge the uncertainty in the dating, however, this uncertainty will not jeopardize the broader picture of the deglaciation and the conclusions in the present work. Original references, as well as information about location, type of material dated and context is given in Supplementary Note 1 and Supplementary Data 1–3.

Ice sheet modelling

Using the new time slices, a 3D reconstruction of the EIS was created using an ice sheet model with the aim of producing a transient history of ice volumes and meltwater delivery through the last deglaciation. The results of this modelling are presented in Figs. 3 and 4 and in Supplementary Data 4–6. Inverse-type, steady-state reconstructions of the evolving ice sheet geometry were created using the ICESHEET 1.0 programme42, based on the assumption of a perfectly plastic ice rheology72,73. Although it is highly unlikely the EIS was ever in a steady state, this method otherwise provides an effective first-order approximation for the 3D ice sheet architecture when climate/ocean forcings and ice dynamics are otherwise poorly constrained. Minimal input parameters are required, and include the ice margin, basal shear stress (Supplementary Fig. 5) and basal topography.

Basal shear stress values are parameterised on the basis of topography and sediment cover, with maximum values applied from the PaleoMIST 1.0 dataset44 and originally guided to the ice sheet evolution according to DATED-12. In lieu of tuning the basal shear stress values to relative sea level records from across the domain, additional experiments perturbing these shear stress values by ±20% were run (Fig. 4a and Supplementary Data 4–6). These PaleoMIST shear stress values were originally guided to the ice sheet evolution according to DATED-12. These sensitivity experiments therefore envelope a probable uncertainty range for the evolving ice-sheet volume, in particular where our marine-limits differ most from that of DATED-1 e.g., in the North Sea and eastern Barents Sea. The model outputs incorporating ±20% and median basal shear stress values are given in Supplementary Data 4–6 and presented for the total volume estimates in Fig. 4a. Otherwise, the median experiments are used as the reference.

The GEBCO2020 grid was used as a base topography. Flexure of the Earth’s crust as a result of ice-sheet loading significantly affects ice volume calculations. Glacial isostatic adjustment (GIA) was therefore calculated by passing an initial iteration of the steady-state ice thickness through SELEN v474, which solves the sea level equation, including accounting for shoreline migration, adjustments for grounding line position, and rotational feedback. Since there is a viscous component to the ice loading response, the ice loading history was supplemented with reconstructed ice thicknesses derived from the DATED-1 margin dataset from 25 to 21 ka BP and the UiT model4,35 from 30 to 26 ka BP. Far field effects from the additional global ice sheets were included using ICE-6G_C outputs. The VM5i Earth rheological profile was used—an ad hoc realisation of the ICE-6G_C (VM5a) GIA model by Peltier et al.75—applied using a Tegmark grid resolution parameter of R = 100 (c. 20 km). The resulting Earth deformation was added to the modern topography to produce a time-transgressive paleo topography. The ice sheet reconstruction was again calculated, using the adjusted bed on a projected coordinate system at 2 km resolution. One iteration of GIA increased total ice volume by ~5%. The global eustatic sea level curve by Waelbroeck et al.76 was used to partition marine- and terrestrial-based sectors of the modelled ice sheet. Contributions to eustatic sea level rise (in metres) were calculated based on the volume of grounded ice above freeboard (Supplementary Fig. 3) and assuming a constant global ocean coverage of 3.618E8 km2.

Ocean temperature timeseries

Two previously published records of planktic foraminiferal species distributions from high-resolution marine piston core records, LINK17 from the North Sea margin52 and JM03-373PC from the western Svalbard margin37 (Fig. 1a), are used. Sea surface temperatures (Fig. 4e) was calculated by transfer functions applying the C2 programme77 and the Weighed Average Partial Least-Squares (WAPLS) method using one component. The calculations were based on the percentage data of species counted in the size-fraction >100 µm and the modern 100-µm database of ref. 78 with the same modifications as in ref. 79. The modern sea surface temperatures were from the World Ocean Atlas 200980. The Root Mean Square Error of Prediction (RMSEP) is 1.9329. This value is caused by the spread of data in the training set; the more outliers removed, the lower the RMSEP. It should be emphasised that due to low or mono specific assemblages of planktonic foraminifera in the arctic, SST estimates with this method has a lower limit close to 4 °C. We acknowledge the simplicity in this approach by using only this proxy as a multitude of different oceanic scenarios including variability in sea ice cover, length of season, seasonality in sea surface conditions, presence of polynyas etc. are not reflected. However, we are confident that by selecting the planktonic foraminifera proxy record from two offshore sites we are identifying the periods with warmest conditions during summer (melting season) in the region. The fact that the same pattern is seen in both cores, from more than 1400 km apart, suggest that they represent a regional pattern. The chronology of the timeseries is based on radiocarbon dates calibrated after the same protocol as the deglaciation dates.