Overview of LIG sea-level evidence

The nature of LIG sea-level variability remains strongly debated, with emphasis on two issues. First, near-field sites (close to the ice sheets) in NW Europe suggest LIG sea-level stability, although resolution and age control remain limited and other N European sites might support sea-level fluctuations25. Second, there is a wealth of global sites (mostly in the far field relative to the ice sheets) that implies LIG sea-level variability (Fig. 1), but which also reveals a striking divergence between site-specific signals with respect to both timing and amplitude of variability (Supplementary Note 1). This suggests that individual sites are overprinted by considerable site-specific influences—e.g., prevailing isostatic, tectonic, physical, biological, biophysical, and biochemical characteristics—rather than reflecting only global sea-level changes. Regardless, a more coherent pattern seems to be emerging from the more densely dated and stratigraphically well-constrained sites, which include the Seychelles, Bahamas, and also Western Australia (Supplementary Note 1, synthesis). The Seychelles coral data are radiometrically precisely dated, avoid glacio-isostatic offsets among sites, and include stratigraphic relationships that unambiguously reveal relative event timings3,22. The Bahamas data comprise stratigraphically well-documented and dated evidence of different reef-growth phases23. Nevertheless, the overall coral-based literature suggests at least two plausible types of LIG history (early vs. late highstand solutions) that remain to be reconciled (Supplementary Note 1, synthesis).

Fig. 1 Global summary of stratigraphic evidence for Last Interglacial sea-level instability in coral-reef deposits and coastal-sediment sequences. Blue dot is the location of Hanish Sill, the constraining point for the Red Sea sea-level record. Red squares with white centres are stratigraphically superimposed coral reef or sedimentary archives for sea-level oscillations within the Last Interglacial (LIG). Solid red dots are locations where sea-level oscillations are inferred but where there is no stratigraphic superposition. The underlying map is of the difference between maximum Last Interglacial (LIG) relative sea level (RSL) values for glacio-isostatic adjustment (GIA) modelling results based on two contrasting ice models (ICE-1 and ICE-3) for the penultimate glaciation using Earth model E1 (VM1-like set up). The ICE-1 model is a version of the ICE-5G ice history (LGM-like), whereas ICE-3 has both reduced total ice volume relative to ICE-1, and a different ice-mass distribution (i.e., a smaller North American Ice Sheet complex and larger Eurasian Ice Sheet) that is consistent with glaciological reconstructions of the penultimate glacial period4 Full size image

Updated Red Sea age model

Regarding the Red Sea RSL record, we improve its LIG-end age control10,18 by comparing the entire dataset (the stack) with radiometrically dated coral-data compilations11,26 and Yucatan cave-deposits that indicate when sea level dropped below the cave (i.e., a “ceiling” for sea level)24. This comparison reveals that the 95% probability limit of the Red Sea stack on its latest chronology10,19 dropped too early (123 ka; see Methods and Supplementary Note 2) relative to the well-dated archives (119–118 ka; Fig. 2b, c; Supplementary Figs. 2 and 3). We, therefore, adjust this point to 118.5 ± 1.2 ka (95% uncertainty bounds) (Fig. 2, Supplementary Figs. 2 and 3), and accordingly revise all interpolated LIG ages with fully propagated uncertainties (Supplementary Fig. 2).

Fig. 2 Variability in Last Interglacial sea-level time-series. Yellow bar: time-interval of Heinrich Stadial 11 (HS11)19. Orange bar: approximate interval of temporary sea-level drop in various records. Dashed line: end of main LIG highstand set to 118.5 ka (cross-bar indicates 95% confidence limits of ±1.2 ka), based on compilations in b and the speleothem sea-level “ceiling” (c). a GrIS contributions to sea level from a model-based assessment of Greenland ice-core data (blue)9, and changes in surface sea-water δ18O at Eirik Drift (black; this study) with uncertainties (2σ) determined from underpinning δ18O and Mg/Ca measurement uncertainties and Mg/Ca calibration uncertainties. b Ninety-five per cent probability interval for coral sea-level markers above 0 m11 (brown), and LIG duration from a previous compilation (black)26. c Red Sea RSL stack (red, including KL23) with 1σ error bars. Smoothings are shown to highlight general trends only, and represent simple polynomial regressions with 68% and 95% confidence limits (orange shading and black dashes, respectively). Purple line indicates the sea-level “ceiling” indicated by subaerial speleothem growth (Yucatan)24. d Probability maximum (PM, lines) and its 95% confidence interval for Antarctic temperature changes (red)68, and proxy for eastern Atlantic water temperature (ODP976, grey)69. Blue crosses: composite record of atmospheric CO 2 concentrations from Antarctic ice cores19. e Individual records for Red Sea cores KL11 (blue, dots) and KL23 (red, plusses), with 300-year moving Gaussian smoothings (as used in ref. 1). Also shown is a replication exercise to validate the single-sample earliest-LIG peak in KL23 (grey, filled squares) with 1 standard error intervals (bars, σ/√{N}, based on N = 5, 5, 4, 4, and 5 replications, from youngest to oldest sample, respectively). Separate blue cross indicates typical uncertainties (1σ) in individual KL11 datapoints prior to probabilistic analysis of the record. f Probabilistic analysis of the KL11 Red Sea RSL record, taking into account the strict stratigraphic coherence of this record. Results are reported for the median (50th percentile, dashed yellow), PM (modal value, black), the 95% probability interval of the PM (dark grey shading), and both the 68% and 95% probability intervals for individual datapoints (intermediate and light grey shading, respectively) Full size image

Estimates of Greenland mass loss

Next, we compare the Red Sea sea-level information (Fig. 2b, c, e, f) with estimates of GrIS-derived LIG sea-level contributions from a model-data-assimilation of Greenland ice-core data for summer temperature anomalies, accumulation rates, and elevation changes9 (Fig. 2a). We add independent support for the inferred late GrIS contribution9, based on a newly extended record of sea-water oxygen isotope ratios (δ18O sw ) from a sediment core from Eirik Drift, off southern Greenland. In this location, δ18O sw reflects Greenland meltwater input with a sensitivity of 4 ± 1.2 m global sea-level rise for the −1.3‰ change seen in the δ18O sw record from ~128 to ~118 ka (Fig. 2a) (Methods, Supplementary Note 3). This record suggests (albeit within combined uncertainties) generally lower GrIS contributions than Yau et al.9, which may agree with results from other modelling studies for GrIS14,15. Both the modelling and δ18O sw approaches indicate a late GrIS contribution to LIG sea level, which is further supported by wider N. Atlantic and European palaeoclimate data, which reveal that contributions started after 127 ka, while GrIS started to regain net mass from 121 ka27.

AIS and GrIS distinction

Although GrIS did not affect LIG sea-level change significantly before 126.5–127 ka (Fig. 2a), the Red Sea and coral data compiled here imply that sea level crossed 0 m at 130–129.5 ka, during a rapid rise to a first highstand apex that was reached at ~127 (Fig. 2b, c, e, f). The Seychelles record indicates specifically that sea level reached 5.9 ± 1.7 m by 128.6 ± 0.8 ka3. We infer that both the first LIG rise above 0 m and the subsequent rapid rise between 129.5 and 127 ka resulted from AIS reduction. Similar qualitative inferences about an early-LIG AIS highstand contribution have been made previously3,9,19, including attribution to sustained heat advection to Antarctica during Heinrich Stadial 11 (HS11; 135–130 ka)19, when a northern hemisphere deglaciation pulse (~70 m sea-level rise in 5000 years) caused overturning-circulation shutdown28, a widespread North Atlantic cold event, and southern hemisphere warming (Fig. 2d). Here we present a quantitative AIS and GrIS separation with comprehensively evaluated uncertainties.

First, we determine centennial-scale LIG sea-level variability from the continuous (and contiguous) single-core RSL record of central Red Sea core KL11 on our new Red Sea LIG age model. We validate this record with new data for high-accumulation-rate core KL23 from the northern Red Sea; i.e., from a physically separate setting than KL11 (Methods) (Fig. 2e). Given this validation, we continue with KL11 alone because it remains the most detailed record from the best-constrained (central) location in the Red Sea RSL quantification method, where δ18O is least affected by either Gulf of Aden inflow effects in the south, or northern Red Sea convective overturning and Mediterranean-derived weather systems in the north16,29.

Second, we perform a Monte Carlo (MC)-style probabilistic analysis of the KL11 record (Fig. 2f), which accounts for all uncertainties in individual-sample RSL and age estimates (cf. blue cross in Fig. 2e). This procedure mimics that applied previously to the Red Sea stack10,18, but now contains an additional criterion of strict stratigraphic coherence (Methods). The analysis leads to statistical uncertainty reduction based on datapoint characteristics, density, and stratigraphy. Remaining RSL uncertainties are ±2.0 to 2.5 m for the 95% probability zone of the probability maximum (PM, modal value; Fig. 2f; Methods).

Both PM and median reveal an initial RSL rise from ~129.5 to ~127 ka to a highstand apex centred on ~127 ka, followed by a drop to a lowstand centred on 125–124 ka at a few metres below 0 m, and then a small return to a minor peak above 0 m at ~123 ka (Fig. 2f). To quantify AIS contributions, we apply a first-order glacio-isostatic correction (with uncertainties) to translate the record from RSL to global mean sea level (GMSL) (Supplementary Note 4) (Fig. 3a), and then subtract the GrIS-contribution records (Figs. 2a and 3b). Our results quantify significant asynchrony and amplitude-differences between GrIS and AIS ice-volume changes during the LIG (Fig. 3b, c). A caveat applies in intervals where the reconstructed AIS sea-level record drops below −10 m, because at that stage the maximum AIS growth limit is approximated (AIS growth is limited by Antarctic continental shelf edges). Whenever the reconstructed AIS sea-level record falls below −10 m (notably after ~119 ka), North American and/or Eurasian ice-sheet growth contributions likely became important. This timing agrees with a surface-ocean change south of Iceland from warm to colder conditions27.

Fig. 3 Identification of Greenland Ice Sheet and Antarctic Ice Sheet contributions to Last Interglacial sea-level variations. a Global Mean Sea Level (GMSL) approximation based on the probabilistically assessed KL11 PM (black line) and its 95% probability interval (grey). This record is shown in terms of RSL in Fig. 2f, but here includes the glacio-isostatic correction and its propagated uncertainty. Black triangles identify limits between which sea-level rises R1, R2, and R3 were measured. Rates of rise with 95% bounds: R1 = 2.8 (1.2–3.7) m c−1; R2 = 2.3 (0.9–3.5) m c−1; R3 = 0.6 (0.1–1.3) m c−1. b Blue: GrIS sea-level contribution from the model-data assimilation of ref. 9 (shading represents the 95% probability interval). Grey: GrIS contribution based on Eirik Drift δ18O sw . Uncertainties as in Fig. 2a. Orange: AIS contribution from subtraction of the blue GrIS reconstruction from the record in a. Green: AIS contribution found by subtracting the grey GrIS reconstruction from the record in a. Orange and green AIS reconstructions are shown as medians (lines) and 95% confidence intervals (shading). Reconstructed AIS contributions cross downward through a fine dashed when they fall below –10 m, which indicates a rough maximum AIS growth limit in terms of sea-level lowering (AIS growth is limited by Antarctic continental shelf edges). When the green/orange curves fall below these limits, North American and/or Eurasian ice-sheet growth is likely implied. The key result from the present study lies in identification of GrIS and AIS sea-level contributions above 0 m. c Southern Ocean ODP (Ocean Drilling Program) Site 1094 authigenic uranium mass accumulation rates, on its original, Antarctic Ice Core Chronology (AICC2012) tuned, age model. Dashed lines indicate potential offsets (within uncertainties) between the ODP 1094 AICC2012-based chronology36 and our LIG chronology (see refs. 10,19 and this study) Full size image

Intra-LIG sea-level variability

Red Sea intra-LIG variations are generally consistent (within uncertainties) in timing with apparent sea-level variations in the well-dated and stratigraphically coherent coral data from the Seychelles, and Bahamas3,22,23, but with larger amplitudes. Northwestern Red Sea reef and coastal-sequence architecture reconstructions offer both timing and amplitude agreement (although age control needs refining)30,31 (Supplementary Note 1). The reef-architecture study in particular30 indicates an early-LIG sea-level rise with a post-128-ka culmination at 5–10 m above present, followed by a millennial-scale ~10 m sea-level drop to a lowstand centred on ~124 ka.

In more detail, the probabilistic Red Sea record suggests a statistically robust dual substructure within the initial LIG sea-level rise (Fig. 2f), which is replicated between Red Sea records (Fig. 2e). It is not (yet) supported in wider global evidence (Methods, Supplementary Note 1), but there are indications that certain systems may have recorded it independently. For example, southwestern Red Sea reef-architecture reveals two main reef phases with a superimposed minor patch-reef phase1,32, reaching total thicknesses up to 10 m. But more precise dating and support from other locations are needed to be conclusive. In this context, we calculate with a basic fringing-reef accretion model that the rapid rises and short highstands inferred here (Fig. 2e, f) may have left limited expressions in reef systems, except for rare ones with exceptionally high accretion rates, or where rapid crustal uplift offset some of the rapid sea-level rises (Supplementary Note 5). Hence, we consider wider palaeoceanographic evidence to evaluate the suggested sea-level history.

Palaeoceanographic support

AIS meltwater pulses implied by sea-level rises R1 and R2 (Fig. 2f) should have left detectable signals around Antarctica. The early-LIG AIS sea-level contribution occurred immediately after Heinrich Stadial (HS) 11, when overturning circulation had recovered from a collapsed HS11 state (Figs. 2–4)28. This likely enhanced advection of relatively warm northern-sourced deep water into the Circumpolar Deep Water (CDW), which impinges on the AIS. At the same time, there was a peak in Antarctic surface temperatures (Figs. 2d and 4c) and Southern Ocean sea surface temperatures (ODP Site 1094 TEX 86 L, ODP Site 1089 planktic foraminiferal δ18O) (Fig. 4c–e), and Southern Ocean sea ice was reduced (Fig. 4b). We infer that early-LIG AIS retreat resulted from both atmospheric and (subsurface) oceanic warming, which—together with minimal sea ice (important for shielding Antarctic ice shelves from warm circumpolar waters, e.g., ref. 33)—drove enhanced subglacial melting rates and ice-shelf destabilisation, and thus strong AIS sea-level contributions between 130 and 125 ka.

Fig. 4 Timing of Antarctic Ice Sheet retreat relative to circum-Antarctic climate and ocean warming. LIG records of a. Antarctic ice core composite atmospheric CO 2 (ref. 70), b EPICA Dome C sea-salt Na flux (on a logarithmic scale), which reflects Southern Ocean sea-ice extent71, c Vostok δD (lilac)67,72, d Site 1089 planktic foraminiferal (G. bulloides) δ18O (red)38, e Site 1094 TEX 86 L-based sea surface temperatures (orange)36, f Site 1089 planktic minus benthic foraminiferal δ18O (‰) plotted as 3-point running mean (red) and sample average including combined 1-sigma uncertainty (light red shading)38, g Site 1094 authigenic uranium (aU) accumulation where higher values indicate bottom-water deoxygenation36, h Site 1063 ε Nd (dark blue, measured by MC-ICP-MS; light blue, measured by TIMS)28, and i bottom-water δ13C records from Site 1063 (blue, 3-point running mean, based on benthic foraminifera Cibicidoides wuellerstorfi, Melonis pompilioides, and Oridorsalis)28, MD03–2664 (yellow, 3-point running mean, C. wuellerstorfi)73, Site 1089 (red, C. wuellerstorfi)36, and Site 1094 (orange, C. wuellerstorfi)36. h and i Indicate North Atlantic Deep Water (NADW) influence as denoted. Map inset includes marine core locations, plotted using Ocean Data View (https://odv.awi.de) Full size image

Wider palaeoceanographic evidence can be used to test the concept that major AIS melt will provide freshwater to the ocean surface, which density-stratifies the near-continental Southern ocean, impeding Antarctic Bottom Water (AABW) formation34,35, which in turn will lead to reduced AABW ventilation/oxygenation and an increase in North Atlantic Deep Water (NADW) proportion vs. AABW proportion in the Atlantic Ocean28,36. Thus, we infer strong support for early-LIG AIS melt from palaeoceanographic observations. For example, an anomaly in authigenic uranium mass-accumulation rates (aU MAR) in Southern Ocean ODP Site 1094 has been attributed to bottom-water deoxygenation (AABW reduction/stagnation), due to strong Antarctic meltwater releases and consequent water-column stratification36 (Figs. 3c and 4g). Also, increased bottom-water δ13C, due to expansion of high-δ13C NADW at the expense of low-δ13C AABW, occurred at the end of HS11 in both the abyssal North Atlantic (ODP Site 1063, core MD03-2664) and South Atlantic (Sites 1089 and 1094) (Fig. 4i). Moreover, ε Nd changes in Site 1063 (ref. 28) support the δ13C interpretation (Fig. 4h). Given that intensification of relatively warm NADW likely plays a key role in subglacial melting and resultant AABW source-water freshening33,37, we infer a positive feedback. In this feedback, meltwater-induced AABW reduction warmed CDW through increased admixture of relatively warm NADW, which then caused further subglacial melting and AABW source-water freshening, driving additional AABW decline. Finally, a distinct early-LIG minimum in the Site 1089 planktic–benthic foraminiferal δ18O gradient indicates a persistent surface buoyancy anomaly, which agrees with strong AIS meltwater input38 (Fig. 4c–f). Surface buoyancy/stratification increase would restrict air–sea exchange and subsurface heat loss. Analogous to explanations offered for high melt rates in some regions of Antarctica today and for even higher melt rates in a warmer future climate39, we therefore propose another positive feedback for the LIG, in which melt-stratification led to subsurface ocean warming, which then intensified ice-shelf melting.

Finally, we note that the aU MAR variations in Southern Ocean Site 1094 (ref. 36) also agree in more detail with our inferred dual substructure in the AIS-related early-LIG highstand (Fig. 3b, c). It is not yet possible to eliminate robustly the inferred offsets (which fall within uncertainties) between the ODP 1094 AICC2012-based chronology36 and our LIG chronology (see refs. 10,19 and this study) (Fig. 3b, c), but the offsets may also (partly) arise from time-lags between meltwater input at the surface and oxygenation decline at the sea floor. Given the position of ODP Site 1094 (South Atlantic sector), the aU MAR record may be to some extent site-specific, in which case it suggests a likely meltwater source from the West Antarctic Ice Sheet (WAIS). The lack of later aU MAR spikes for our further inferred AIS contribution may then suggest either that most of WAIS had been lost during the earliest LIG, or that it had at least retreated far enough to stop contributions as is also indicated by ice-sheet studies14,40,41,42,43.