New age constraints for glacial terminations IV, III, and III.a based on western Mediterranean speleothem records

The full understanding of climate feedbacks responsible for the amplification of deglaciations requires robust chronologies for these climate transitions, but in the case of marine records, radiocarbon chronologies are possible only for the last glacial termination. Although the assumed relationships between the marine isotopic record and the orbital parameters provide a first-order chronology for glacial terminations, an independent chronological control allows the relationships between orbital forcing and the climate response to be evaluated over multiple previous terminations. To assess this, we present geochemical records from the western Mediterranean, including two speleothems and one marine sediment core. The most notable speleothem, the so-called RAT, established a new long terrestrial climate record for this region, spanning Marine Isotope Stages from MIS 11 to MIS 7. Its absolute U $/$ Th dates provide an exceptional chronology for the glacial terminations IV, III, and III.a. The onset of these three glacial terminations was marked by rapid δ¹⁸O depletions, reflecting ocean freshening by ice melting, thus providing an excellent tie point for regional marine records also sensitive to such freshening. This is exemplified by new δ¹⁸O data of the Ocean Drilling Program (ODP) 977 site from the Alboran Sea, where the speleothem chronology was employed to adjust its age model. The new chronologies reveal an earlier onset of the deglacial melting for the TIV and TIII.a that is in contrast to the generally accepted marine chronologies and indicates that the duration of these deglaciations was variable, with TIV being particularly longer (∼ 20 kyr). This study also supports that the onset of deglacial melting always occurred during a declining precession index, while a nonunique relation occurred with the obliquity parameter.

1 Introduction

Glacial–interglacial transitions exemplify the highest amplitude of climate changes along the Late Quaternary and involve intense and fast processes of ice sheet melting that led to rapid sea level rise (Past Interglacials Working Group of PAGES, 2016). Although the increasing Northern Hemisphere summer insolation (NHSI) is considered the main driver for deglaciations, internal feedbacks from the Earth's climate system were crucial to trigger glacial terminations (Cheng et al., 2009; Drysdale et al., 2009; Tzedakis et al., 2012; Barker and Knorr, 2021). The consequences of the insolation rise depended on the actual extent of the ice sheet build-up during the preceding glacial period, determining the main features of terminations (Tzedakis et al., 2017). Ice sheet instability and intense events of massive icebergs released into the North Atlantic, named Heinrich-like events, are considered characteristic features of deglaciation (Barker et al., 2015, 2019; Barker and Knorr, 2021; Hodell et al., 2013, 2023) and have been well-documented through ice-rafting debris (IRD) records in the North Atlantic (Hodell et al., 2008; Barker et al., 2011). The consequent North Atlantic freshening caused by ice melting became the trigger for powerful feedback processes inducing glacial terminations. The freshwater discharge allowed surface water stratification and provoked the displacement of the deep-water formation, disturbing the strength of the Atlantic Meridional Overturning Circulation (AMOC) (Rahmstorf, 2002). This, in turn, changed the interhemispheric heat transport, inducing a rapid cooling in the North Atlantic region, whereas Antarctica gradually warmed. The asymmetry in the hemispheric heat transport (the bipolar seesaw) and its influence on the atmospheric CO₂ rise was also a fundamental requirement for deglaciation to occur (Cheng et al., 2009; Barker et al., 2011, 2019; Barker and Knorr, 2021; Tzedakis et al., 2012). Thus, variations in the AMOC and greenhouse gases explain much of the variability in the global climate during deglaciations (Clark et al., 2012). In parallel to the weakening of the AMOC, albedo changes due to the reduction in the continental ice sheet extent operated as secondary internal feedback or amplifier. At the end of the deglaciation arises the characteristic abrupt warming that leads to the onset of interglacial periods, but the sea level highstand marking the interglacial conditions occurred systematically a few thousand years after the maximum in the NH summer insolation due to the slow response of the ice sheets to insolation (Past Interglacials Working Group of PAGES, 2016). Several warm marine isotopic substages are considered interglacial acmes or peaks of optimum climate conditions (i.e., the Marine Isotopic Stages (MISs) 7a–c, 7e, 9e, and 11c within them). However, temporal and geographical asynchronicities in the precise definition of their boundaries between marine and terrestrial records or even different proxies within the same archives pose a challenge to resolving the feedbacks acting in these transitional periods (Shackleton, 1969, 2003; Sánchez-Goñi et al., 1999). These feedbacks responsible for past deglacial warming can also play a role in the ongoing situation of climate change. The increase in atmospheric CO₂ due to anthropogenic emissions, the reduction in Earth albedo resulting from ice melting, and even the weakening of the AMOC are currently in operation (IPCC, 2023; Boers, 2021). Studying past glacial terminations provides unique examples for elucidating the role that these internal components of Earth's climate system had during past periods of global climate warming, providing a testing ground for better contextualizing current climate change.

Astronomically tuned marine sediment records provided fundamental information to understand the evolution of the continental ice sheets in the past. Globally averaged δ¹⁸O stacks from benthic foraminifera, reflecting changes in ice sheet build-up and deep-ocean temperature, have provided a global chronological framework for correlating marine sediment sequences, and the marker is its global signal driven by changes in the ice sheet build-up (Imbrie et al., 1984; Martinson et al., 1987; Lisiecki and Raymo, 2005). However, the chronology applied to such stacks is based on the assumption of a stationary relationship between insolation and ice volume changes and does not account for suborbital climate influence in the time sequences of each glacial termination (Hodell et al., 2015). While also reflecting ice volume cycles, the seawater oxygen isotopic signal (δ¹⁸O_sw) from planktonic records can be overprinted by regional temperature or hydrographic changes, thus challenging its use as an isotopic chronostratigraphic tool (Rohling et al., 2015). Hence, orbitally dependent marine isotopic chronologies are limited in resolving precise age constriction of key climate periods, such as glacial–interglacial transitions, complicating any interpretation in terms of orbital forcing.

A growing number of studies uses absolute dating from speleothems to anchor the marine records and achieve astronomically independent age models (e.g., Grant et al., 2012; Cheng et al., 2016). Speleothems with a relatively unambiguous climatic interpretation of their proxies, continuously spanning several glacial–interglacial cycles, are restricted to a few studies. For instance, speleothem records from North and South America support orbital parameters as the predominant driver of glacial terminations, although they are limited to the two most recent cycles (Cruz et al., 2005, 2009; Moseley et al., 2016). Another notable example is the case of Chinese speleothems that provided a great paleoclimatic assembled record of the Asian monsoon lasting 640 kyr and encompassing various deglaciations (Yuan et al., 2004; Cheng et al., 2006, 2009, 2016). They even reported first insights into the existence of the extra glacial termination TIII.a, which reproduces patterns of events equivalent to other glacial terminations despite being amid MIS 7 (Cheng et al., 2009). Superimposed on a strong precessional δ¹⁸O signal forced by regional insolation and sea surface temperature gradients, these precisely dated speleothem records from Sanbao cave in China show positive δ¹⁸O anomalies interpreted to reflect changes in the position and timing of monsoon rainfall coincident with the weakening of the AMOC in the North Atlantic (Wang, 2001; Cheng et al., 2009; Cai et al., 2015; Barker and Knorr, 2021; Jian et al., 2022). If AMOC weakening is linked to glacial meltwater release in the North Atlantic, these records may, therefore, provide indirect indicators of the chronology of glacial terminations. However, the distinctive freshening signature of the glacial terminations and therefore the specific time sequence of its events cannot be recorded in these speleothems located far from the North Atlantic water sources. Hence, well-constrained speleothem records capable of recording the North Atlantic deglacial features through its isotopic signature would be advantageous to establish independent chronology on North Atlantic marine sequences and clarify the forcings of deglaciations. In this respect, the first high-resolution speleothem records successfully used to tune marine records were collected in the NW Iberia peninsula and documented rapid North Atlantic ice sheet melting spanning the two ultimate glacial terminations (TI and II) at millennial–centennial timescales (Stoll et al., 2022).

Figure 1 Map showing the location of this study records. (1) Murada cave, on Menorca, where the RAT speleothem was recovered. (2) Vallgornera cave on Mallorca, where the VALL 2 speleothem was recovered. (3) The marine sediment ODP 977 site. Other sites discussed in the text are also shown as follows. (4) Campanet cave on Mallorca (Dumitru et al., 2018). (5) North Iberian Speleothem Archive (NISA) and Garth speleothem records from north Iberia caves (Stoll et al., 2015, 2022). (6) Corchia cave (Drysdale et al., 2004; Tzedakis et al., 2018). (7) Ejulve cave (Pérez-Mejías et al., 2017). (8) The Tenaghi Philippon arboreal pollen record (Tzedakis et al., 2006). (9) Sanbao speleothem record (Cheng et al., 2009, 2016). (10) The ODP 983 site (Barker et al., 2015, 2019; Barker and Knorr, 2021). (11) The MD01-2444 marine core (Tzedakis et al., 2018; Hodell et al., 2023). The map shows the dominant atmospheric circulation that brings moisture to the Balearic islands (via Mediterranean back-door cold fronts and winter storm tracks). Moisture in Asia is related to the dynamics of the East Asian summer monsoon (EASM) and the Indian summer monsoon (ISM), which depends on the Intertropical Convergence Zone (ITCZ). The map shows the modern average position of the ITCZ in summer.

Regarding the Mediterranean region, available speleothem records that can be potentially used for marine alignments are mostly restricted to the two most recent deglaciations (e.g., Bar-Matthews et al., 1999, 2000, 2003; Drysdale et al., 2009; Tzedakis et al., 2018). In the western Mediterranean basin, only the Ejulve speleothem record from the eastern Iberian Peninsula detected millennial-scale events related to AMOC slowdowns at TIII (Pérez-Mejías et al., 2017), while speleothems spanning TIV and TV remain elusive in the Mediterranean context (Kaushal et al., 2024). This study presents two new speleothem records from the Balearic islands located in the western Mediterranean region (Fig. 1). The VALL 2 speleothem was recovered from Mallorca and encompassed the TV, whereas the RAT speleothem, recovered from Menorca, encompassed two glacial–interglacial cycles with the chronologies for the glacial terminations TIV, TIII, and TIII.a being well-constrained. Hence, the RAT speleothem offers a unique continuous record not provided by any existing stalagmite archive within this region. A geochemical multiproxy approach has been applied in these two speleothems to investigate climate variability, with a special focus on their glacial terminations. The studied deglaciations present different intensities and durations and open an opportunity to investigate the relationship between orbital parameters and paleoclimatic records. We apply the RAT chronology to newly acquired planktic isotope measurements from the Ocean Drilling Program (ODP) 977 site in the Alboran Sea, highlighting the potential of this speleothem record as a referential base for marine chronologies. This approach aims to avoid the circularity of trying to understand the drivers and thresholds of climate variability using marine chronologies dependent on orbital tuning.

2 Material and methods

2.1 Cave settings and regional hydroclimate

The speleothem called Ratpenat (RAT) was collected from inside Murada cave in Menorca (Fig. 1). The Murada cave is located close to the coast and at a low altitude, 80 m a.s.l., in the southwest area of the island (39°57^′58^′′ N, 3°57^′53^′′ E). This cave is the largest in the Barranc d'Algendar, which is a canyon formed by a ravine cutting Miocene calcarenite rocks, which is where a karst was developed. Nowadays, this cave represents a fossil level of the karst drainage conduit. The RAT speleothem was found broken in two pieces, which were matched in the cave. The VALL 2 speleothem grew inside the d'Es Pas de Vallgornera cave, which is located in the southern part of Mallorca (39°22^′00^′′ N, 2°52^′25^′′ E) (Fig. 1). This cave was also developed in upper Miocene rocks and consists of more than 67 km of conduits and chambers, some of them flooded. Sea level reconstruction by means of overgrowths on speleothems has been previously studied in this cave (Dorale et al., 2010).

According to the Köppen climate classification, the Balearic islands present a temperate climate (Csa) with seasonal rainfall, intense storms in autumn, and dry/hot summers that cause the lack of permanent watercourses. The two main origins of the precipitation on the western Mediterranean Sea are the Mediterranean Sea itself and the tropical–subtropical North Atlantic corridor via the dominant eastward atmospheric circulation, leading storm tracks reaching the Mediterranean (Gimeno et al., 2010; Nieto et al., 2010; Krklec and Domínguez-Villar, 2014; Dumitru et al., 2017; Moreno et al., 2021). However, the major source of moisture for precipitation in the Mediterranean surrounding land area has been identified as the basin itself (Schicker et al., 2010). Particularly in the location of the Balearic islands, the moisture is provided by evaporation over the Mediterranean because of the occurrence of back-door cold fronts related to easterly advection leading to heavy precipitation in autumn (Millán et al., 2005; Moreno et al., 2021). The Mediterranean region has important climate contrasts among different nearby areas because it is a transitional zone between subtropical and midlatitude temperate regimens (Lionello, 2012). While the northern part of the Mediterranean region is strongly linked to midlatitude variability, mostly characterized by the North Atlantic Oscillation, the southern and eastern parts are exposed indirectly to the influences of the Intertropical Convergence Zone (ITCZ) dynamics.

2.2 Geochemical analyses of speleothems

2.2.1 Geochronology

In order to determine the absolute ages, both speleothems were longitudinally cut with a saw, taking into account the morphology and the growth axis. The sampling for measuring uranium and thorium isotopes was designed considering the layering to detect possible hiatuses and changes in growth rates. RAT is a 50 cm long speleothem broken at the top; a 2 cm calcite piece was missing when the sample was found in the cave. It does not reveal important changes in the orientation of the growth axis but presents a prominent orange–brown layer at 14 cm from the tip and white layers at 19 and 33 cm (Fig. 2). The VALL 2 speleothem is a wide and dark brown speleothem that is 22 cm long (Fig. A1 in the Appendix). The color at the bottom of the speleothem is light brown, with a porous texture that, towards the top, becomes darker and less porous. The axis presents several changes in the growth direction and possible discontinuities are characterized by light layers. After the last main discontinuity, the speleothem presents two independent growths. The axis selected to accomplish the analyses follows the most developed axis of the last two growths.

Figure 2 On the right, there is an image of the RAT speleothem showing the positions of the radiometric dates sampled for creating the age model (green) and the reversal not used (red). The growth axis path used for the geochemical analyses is highlighted with a white line. At the top left, the age model acquired by StalAge (dark blue line) and the original U $/$ Th dates for the RAT speleothem (dark blue dots) are plotted with the associated errors (error bars and green area). RAT growth rates are shown in dark green. In the lower panel, the results of δ¹⁸O, δ¹³C, and Mg $/$ Ca ratios are accompanied by uncertainties represented as a shaded area. Note that the three parameters are plotted with a reversed axis. The stars highlight negative δ¹⁸O anomalies mentioned in the main text. The pink bars highlight slow growth rate periods, and the dashed grey line points to major color changes in the speleothem layers, also noted as dashed grey arrows in the image.

Samples for U $/$ Th analyses were milled with a 2 mm diameter tungsten carbide micro-drill under a laminar airflow cabinet to prevent contamination. About 100–200 mg of carbonate powder was obtained for each sample. The chemical procedure for isolating uranium and thorium elements was performed in the ultra-clean laboratory at the University of Minnesota and Xi'an Jiaotong University, applying the methodology previously described by Shen et al. (2002) and Cheng et al. (2009) (Appendix A). Several blanks were performed routinely for each laboratory sampling set. The radiogenic isotopes of U and Th were analyzed using a parallel ion-counting multi-collector inductively coupled plasma mass spectrometry (MC-ICP-MS) with a Finnigan Neptune model at both universities. The results of these analyses and the absolute dates are detailed in Table A1. The speleothem depth–age models were performed with R software and the StalAge package (Scholz and Hoffmann, 2011), using 39 radiometric dates for RAT and 12 for VALL 2. The RAT speleothem provided a consistent chronology, with a few single dates exceeding the final modeled dating uncertainty. It presents a remarkably precise age model with absolute uncertainties being even smaller in the youngest part of the speleothem (Fig. 2). In contrast, the dates obtained from the VALL 2 speleothem are approaching the limits of the U-series dating technique, thereby resulting in significant uncertainties that are even larger in the oldest section of the speleothem, where the growth rate is higher (Fig. A1). The growth rates range between 0.02–1.2 cm kyr⁻¹ for VALL 2 and 0.07–1.7 cm kyr⁻¹ for RAT.

2.2.2 Geochemical measurements

Stable isotopes and trace elements have been measured on aliquots of the same powder sample in both speleothems. The δ¹⁸O, δ¹³C, and Mg $/$ Ca ratios have been compared in detail to evaluate their paleoclimatic signals. Samples were drilled using a 0.5 mm diameter tungsten carbide dental burr along the speleothem growth axis. The sample resolution is 2 mm in VALL 2 and 2.5 mm in RAT. However, samples are spaced every 0.5 mm during glacial terminations in the RAT speleothem to increase the temporal resolution due to low growth rates, ensuring that events are captured by more than one measurement.

The δ¹⁸O and δ¹³C measurements were performed with 50 µg of carbonate powder and were accomplished on a Finnigan MAT252 mass spectrometer (isotope ratio mass spectrometer, IRMS) at the Scientific and Technological Centres at the University of Barcelona (CCiTUB). This mass spectrometer is coupled in a single acid bath CarboKiel-III carbonate preparation device to convert solid carbonate samples into simple gas (CO₂) before entering into the ion source of the IRMS. Analytical uncertainties were obtained utilizing two in-house carbonate standards that were calibrated to NBS19 international standard (Coplen, 1996). The uncertainties were ≤ 0.04 ‰ Vienna Pee Dee Belemnite (VPDB) for δ¹³C and 0.08 ‰ VPDB for δ¹⁸O.

Around 3 mg of carbonate powder was used to measure Ca, Mg, and other trace elements such as Ba and Sr to obtain element $/$ calcium ratios. The powder was directly dissolved with 2 mL of ultrapure 2 % HNO₃ (with Rh and Sc as internal standards) and centrifuged to prevent possible lithic particles or impurities in the solution before the analyses. These samples were analyzed employing inductively coupled plasma mass spectrometry (ICP-MS), with an Agilent 7500CE model, at the CCiTUB. Possible contamination was controlled by means of chemistry blank analyses on random days. An in-house high-purity standard solution was measured routinely every four samples to control the accuracy of the measurements. According to the in-house standard analysis, the external reproducibility obtained for the speleothem analyses was 0.86 % (%, 2σ).

2.3 δ¹⁸O Analyses in planktic foraminifera from ODP 977 site

ODP site 977 was retrieved during leg 161 of the Ocean Drilling Program on board the research vessel (R/V) JOIDES Resolution. It was located halfway between the Spanish and Algerian coasts, in the eastern sub-basin of the Alboran Sea, at 36°01.92^′ N and 1°57.32^′ W and at a water depth of 1984 m b.s.l. (meters below sea level; Fig. 1). This study extends the oxygen isotope data obtained in previous studies of this sediment core (Martrat et al., 2004, 2007; Gonzalez-Mora et al., 2008; Torner et al., 2019; Azibeiro et al., 2021). The ODP 977 age models of all these different studies are mostly based on the original age model from Martrat et al. (2004), which is established on the correlation between its Globigerina bulloides δ¹⁸O curve and the SPECMAP stacked curve (Martinson et al., 1987). However, in Torner et al. (2019) some additional tie points were used to align the TII with other Mediterranean cores radiometrically tuned with speleothems. Azibeiro et al. (2021) adjusted the age model along the TV by aligning its alkenone-based sea surface temperature (SST) record to North Atlantic midlatitude temperature records.

Stable isotopes on planktic foraminifers were performed on handpicked specimens of G. bulloides. Shells were crushed between two glass plates to open the chambers and favor the following clay removal. Afterward, samples were cleaned with methanol and sonicated for a few seconds. The supernatant solution was removed with a micropipette, and finally, the samples were dried. The measurements were accomplished by isotope ratio mass spectrometry at the CCiTUB using the same analytical method that was previously explained for speleothems analyses.

3 Results and discussion

3.1 The climate signal in the speleothems

Mg $/$ Ca ratios and stable isotopes in both speleothems reveal climate variability in the western Mediterranean region. According to the obtained depth–age models, VALL 2 spans from 530 to 345 ka (MIS 13–MIS 11), encompassing the glacial termination V, while RAT spans from 387 to 209 ka (MIS 11–MIS 7), continuously covering an exceptionally long time interval (178 kyr) that includes three glacial terminations (TIV, TIII, and TIII.a). Glacial conditions are characterized by slower growth rates and high δ¹⁸O, δ¹³C, and Mg $/$ Ca, while the warm MISs feature low isotopes and Mg $/$ Ca ratios (Fig. 2 and Appendix A1). The RAT record presents an overall variation of 3.2 ‰, 6.3 ‰ and 24.6 mmol mol⁻¹ for δ¹⁸O, δ¹³C, and Mg $/$ Ca, respectively. And the overall variation is 3.4 ‰, 3.8 ‰, and 15.7 mmol mol⁻¹ in the VALL 2 speleothem. Given the insufficient robustness of the VALL 2 age model for the precise evaluation of climatic events, the discussion is primarily focused on the RAT results.

The Mg $/$ Ca ratio in speleothems is sensitive to variations in the precipitation/evaporation balance, reflecting changes in the amount of water infiltrating into the karst system. Frequently, arid conditions, characterized by reduced water infiltration, may increase the groundwater residence time in the karst and even reduce drip rates in the caves (Fairchild et al., 2000). This scenario provides an opportunity to enhance the degassing of groundwater, thereby facilitating carbonate precipitation prior to the arrival of drip water at the surface of the speleothem. This process, known as prior calcite precipitation (PCP), can also be favored by reductions in the partial pressure of CO₂ (pCO₂) within the cave atmosphere. Ventilation typically reduces the CO₂ levels in the cave and occurs when external cooling provokes an exchange of air between the internal and external atmosphere. This mechanism typically occurs on a seasonal scale, and, presuming that the cave structure has remained constant over time, this influence would persist throughout glacial and interglacial periods, exerting a minor impact on the millennial timescale. Moreover, cold and arid conditions cause a decline in vegetation activity, leading to lower respiration rates and diminished CO₂ liberation. Briefly, processes that reduce pCO₂ in the karstic system support PCP. This PCP can preferentially remove Ca from the drip waters, consequently elevating the trace element $/$ Ca ratio of all incompatible elements in the drip waters, such as Mg $/$ Ca, Sr $/$ Ca, and Ba $/$ Ca. When the partitioning coefficients of Mg, Sr, and Ba in stalagmite calcite are relatively stable, the speleothem Mg $/$ Ca, Sr $/$ Ca, and Ba $/$ Ca reflect this drip water PCP signature. The RAT Sr $/$ Ca and Ba $/$ Ca records show similar temporal variations to the Mg $/$ Ca variability (Fig. A2), suggesting that all three elements are dominated by a common process of PCP. The magnitude of change in Mg $/$ Ca is much larger than what would arise due to temperature-driven changes in Mg partitioning, which, in any case, would elevate Mg incorporation during warm interglacials and is in contrast to the observed pattern. PCP can also be enhanced by higher drip water saturation states, even with constant drip rates (Stoll et al., 2012), but the observation of maximum PCP during cold glacial periods characterized by lower soil CO₂ due to diminished microbial activity and vegetation respiration rates and likely lower initial drip water saturation suggests that elevated glacial PCP is caused by slower drip rates and is driven by lower infiltration and lower precipitation–evapotranspiration.

Speleothem δ¹³C can also reflect hydroclimate variations. The initial δ¹³C of the drip water reflects a mixing signal derived from the carbon isotopic composition signatures of limestone, atmospheric CO₂, and the biogenic soil gas derived from both the microbial activity and the vegetation respiration (Hellstrom et al., 1998). Soil autotrophic and heterotrophic respiration typically have lower δ¹³C values. Hence, cold and dry conditions can depress soil respiration rates and contribute to an increase in the δ¹³C. Likewise, the concentration of atmospheric CO₂ (with a higher δ¹³C signal) would be proportionally higher than the relatively low biogenic carbon signature in the soil gas, elevating the δ¹³C values even more. It is for that reason that many midlatitude speleothem studies interpret high δ¹³C as reflecting reduced soil/vegetation activity (Genty et al., 2013; Moreno et al., 2010; Bartolomé et al., 2015; Pérez-Mejías et al., 2017). Additionally, slowed drip rates during dry periods may contribute to greater PCP, which also causes drip water to evolve to higher δ¹³C due to the gradual loss of the light isotope. Thus, via these two processes (initial drip water changes and PCP), dry conditions may contribute to higher speleothem δ¹³C. In the RAT record, while cold temperatures likely depressed soil respiration during glacials and contributed to the high δ¹³C, the general co-variation between Mg $/$ Ca and δ¹³C suggests that the high glacial δ¹³C also reflects drier conditions contributing to low soil CO₂ and slower drip rates favoring more extensive PCP. In another speleothem record from Mallorca, high δ¹³C was similarly interpreted to reflect low biogenic CO₂ productivity in soil and/or a short water residence time, reflecting dry conditions (Dumitru et al., 2018). The RAT δ¹³C signal shows a smoothed signal with respect to the Mg $/$ Ca during glacial terminations, suggesting a slower response of the vegetation compared with precipitation rises (Fig. 2). MIS 11, MIS 9, and MIS 7 were characterized by enhanced precipitation conditions that led to higher soil respiration activity under more extensive vegetation.

The δ¹⁸O RAT record is marked by higher δ¹⁸O during glacials and lower δ¹⁸O during interglacials (Fig. 2). On a global scale, the retreat of the large continental Northern Hemisphere ice volume and release of meltwater lowers the δ¹⁸O_sw in the surface ocean. This lower δ¹⁸O_sw causes a corresponding decrease in the δ¹⁸O of rainfall and the drip waters from which the speleothem is formed, lowering the speleothem δ¹⁸O. Over the most recent glacial cycles, the North Atlantic experiences a greater amplitude change in δ¹⁸O_sw than the globally averaged surface ocean (Waelbroeck et al., 2014). In some regions, major changes in atmospheric circulation patterns cause additional larger changes in rainfall δ¹⁸O which significantly exceed the effect of the surface δ¹⁸O_sw. In coastal regions of the North Atlantic, the change in the surface ocean isotopic composition is a dominant component of speleothem δ¹⁸O variations over glacial terminations (Stoll et al., 2022). A comparison of independently dated marine δ¹⁸O_sw with speleothem δ¹⁸O records from NW Iberia showed that the speleothem δ¹⁸O matched the timing and amplitude of the δ¹⁸O_sw shift in the eastern North Atlantic during termination I (Stoll et al., 2022). Furthermore, several studies demonstrate that the North Atlantic δ¹⁸O_sw signal, changed by those freshwater inputs, was transferred into the Mediterranean Sea (Cacho et al., 1999; Sierro et al., 2005; Jiménez-Amat and Zahn, 2015; Marino et al., 2015). Over termination I, the North Atlantic δ¹⁸O_sw decreased rapidly at around 16 ka and between 12 and 10 ka, and these two shifts are also recorded in the Alboran Sea, as well as around Menorca (Fig. A3). Consequently, we propose that the RAT δ¹⁸O records the change in δ¹⁸O_sw in the moisture source area for the precipitation, encompassing global changes in the δ¹⁸O_sw at a large scale and superimposing the freshening signal transmitted to the Mediterranean Sea during deglaciations. This is also supported by the evaluation of the isotopic composition of rainfall waters and cave drip waters in Mallorca, which points to the western Mediterranean as the dominant source for moisture uptake in coastal caves (Dumitru et al., 2017). The cave location of this study, close to the coast and at a low altitude above sea level, should minimize distillation processes, and therefore, RAT δ¹⁸O captures the signal of the moisture source, which is surface water mostly from the western Mediterranean.

3.2 Mediterranean precipitation modulated by precession

The RAT Mg $/$ Ca record shows a clear precession pattern. Low Mg $/$ Ca, indicating enhanced precipitation, coincides with precession minima, which are times of the high Northern Hemisphere summer insolation (NHSI) (Fig. 3; pink bars), and vice versa for drier periods. The most humid periods are recorded during the strong precession minima of the marine isotopic sub-stages of MIS 9e, MIS 9a, and MIS 7e. MIS 11a shows similar humidity despite a moderate precessional minimum. The low resolution of the Mg $/$ Ca record during the MIS 9a precludes us from evaluating its duration and structure. In the better-resolved interglacial acmes of MIS 9e and MIS 7e, we estimate a duration of ∼ 10 kyr. These orbitally paced changes in Mediterranean precipitation are explained because minima precession led to higher seasonality (high NHSI – low Northern Hemisphere winter insolation, NHWI), which intensified wintertime storm tracks and increased the Mediterranean precipitation. Additionally, higher air–sea temperature gradients enhanced evaporation from the Mediterranean basin itself, leading to greater local convective precipitation (Kutzbach et al., 2014; Wagner et al., 2019; García-Alix et al., 2021).

Figure 3 The RAT speleothem U $/$ Th ages with the 2σ error bars are found at the top of the figure, together with the speleothem dates from the Ejulve and Corchia caves. (a) Precession curve with vertical pink bars highlighting minima precession (Laskar et al., 2004). (b) RAT Mg $/$ Ca record with periods of enhanced precipitation marked with horizontal black lines which correspond to the interglacial/interstadial marine isotopic stages (MISs). (c) The δ¹⁸O_speleo Corchia record (Drysdale et al., 2004; Tzedakis et al., 2018) and the δ¹³C_speleo record from the Ejulve cave (Pérez-Mejías et al., 2017). (d) The Tenaghi Philippon arboreal pollen record from Greece (Tzedakis et al., 2006). (e) The δ¹⁸O_speleo record from the Sanbao cave in China with the U $/$ Th ages and 2σ error bars (Cheng et al., 2016).

This new and uninterrupted record covers enough time to confirm the precession imprint for the first time in a Mediterranean speleothem. Other European speleothems have also recorded enhanced precipitation during precession minima of the studied interglacials; however, they provide shorter records. The iconic δ¹⁸O record of the Antro del Corchia cave, from Italy, shows periods of greater humidity that are in line with those recorded in RAT during the MIS 9 (Drysdale et al., 2004; Tzedakis et al., 2018). However, higher U $/$ Th date uncertainties could account for some chronological discrepancies (Fig. 3). The increment in moisture availability in western Italy was attributed to regional warm SST that led to enhanced evaporation (Drysdale et al., 2004). Carbon isotope records from other European speleothems are more complex to interpret without accompanying trace elements, but the trends are consistent with warmth and/or greater humidity during MIS 7e. The δ¹³C speleothem record from the Ejulve cave, in the northeastern Iberian Peninsula, shows lower values during MIS 7e (Fig. 3) and is interpreted as a proxy for vegetation productivity (Pérez-Mejías et al., 2017). Considering that the typical dating error in the RAT record along the MIS 7e averages ∼ 1.2 kyr and that the Ejulve error record averages ∼ 2 kyr, both records strengthen the enhanced humidity along this period. Other speleothem records from the Austrian Alps (∼ 0.3 kyr average errors) and Sardinia (∼ 1.5 kyr average errors) also reflect warming or more humid conditions along the MIS 7e responding to changes in the North Atlantic realm (Spötl et al., 2008; Columbu et al., 2019; Wendt et al., 2021).

Similarly, long-term pollen records of the Mediterranean region show expansions of the arboreal vegetation coincident with precession minima because of warmer climate and greater water availability (Sánchez-Goñi et al., 2008; Tzedakis et al., 2006, 2009; Wagner et al., 2019). In particular, the arboreal pollen percentages in the Tenaghi Philippon pollen record from Greece, in the eastern Mediterranean region, show an extraordinary resemblance with the precipitation variability recorded in the RAT Mg $/$ Ca, where high arboreal pollen percentages are in line with enhanced precipitation (Fig. 3) (Tzedakis, 2003; Tzedakis et al., 2006). However, certain time discrepancies can be observed mostly during the end of the glacial MIS 10 and for the TIV and TIII.a. The Tenaghi Philippon age model relies on orbital tuning, assuming a constant lag of 1000 years in the orbital alignment, a circumstance that may account for certain discrepancies with the absolutely dated speleothem record. Aside from the age biases during deglaciations, both records exhibit proxy fluctuations associated with glacial/interglacial cycles, as well as stadial/interstadial structures with a strong precession signal.

Precession has a major effect on the Mediterranean Sea due to its ability to disrupt the evaporation-precipitation balance of the sea basin. This influence has been documented through several proxy records that demonstrate modifications in oceanographic and sedimentary processes. For instance, marine sediment records close to the Strait of Gibraltar recorded changes in the Mediterranean Outflow Water (MOW) strength, mainly driven by precession-controlled fluctuations in the Mediterranean hydrologic budget (Nichols et al., 2020; Sierro et al., 2020), and marine sediments from the Balearic Sea documented carbonate cyclicities due to the control of the western Mediterranean cyclogenetic mechanisms (Ochoa et al., 2018). The precession control on the African monsoon is well known, which indirectly also affected the eastern Mediterranean water properties through river runoff along the north African coast, drastically changing the δ¹⁸O_sw composition and allowing sapropel formations (Rohling et al., 2014; Grant et al., 2016). Several studies highlight that the timing of enhanced summer monsoon corresponds with maximum Mediterranean autumn/winter storm track precipitation (Kutzbach et al., 2014; Toucanne et al., 2015; Wagner et al., 2019). The Mg $/$ Ca RAT record also reveals this concordance, where western Mediterranean enhanced precipitation, during interglacial/interstadial periods of precession minima, is consistently in line with higher summer monsoon intensities according to the Chinese speleothem records from Sanbao cave (Fig. 3) (Cheng et al., 2016).

3.3 Onset of glacial terminations

The newly obtained RAT records provide new constraints on the absolute age for glacial terminations TIV, TIII, and TIII.a based on the identification of meltwater release reflected in negative δ¹⁸O excursions (Fig. 4a). Accordingly, the onset of TIV shows a first δ¹⁸O depletion at 351 ± 5 ka, followed by a larger depletion at 343 ± 3 ka, which likely reflects the maximum rate of ice melting. This is consistent with a study based on volcanic interbed deposits of near-coastal aggradational successions in Italy (Fig. 4a; horizontal purple bars), which provided independent ages for sea level rise by means of ⁴⁰Ar $/$ ³⁹Ar absolute dating (Marra et al., 2016, 2019). The onset of the TIII in the RAT speleothem is also in concordance with the sea level rise observed in Italy, and the principal δ¹⁸O negative anomaly occurred at 245 ± 2 ka. This anomaly is comparable in magnitude to the TIV freshening. The TIII.a deglaciation is marked by the negative δ¹⁸O anomaly centered at 227 ± 1 ka (Fig. 4). This characteristic δ¹⁸O feature of deglaciations is also recognized in the Mallorca speleothem (VALL 2) along the TV, although the high age uncertainties associated with the chronological limit of the U–Th technique preclude a precise date for the onset of this termination (Fig. A1).

Figure 4 Comparison between glacial termination TIV (right panel) and glacial terminations TIII and TIII.a (left panel). On the top, the RAT U $/$ Th ages with the 2σ error bars. The vertical blue bars highlight the duration of glacial terminations according to the negative δ¹⁸O structures of the RAT record, while the internal and thinner dark blue bars emphasize the maximum rate of freshening discussed in the main text. (a) The δ¹⁸O_speleo record from RAT speleothem. The horizontal purple lines denote periods of sea level rise according to Marra et al. (2016, 2019). (b) Orbital parameters (Laskar et al., 2004) and the caloric summer insolation curve with maximums are highlighted with red triangles (Tzedakis et al., 2006). (c) The LR04 stack of δ¹⁸O benthic foraminifera (Lisiecki and Raymo, 2005) and the North Atlantic IRD record of the ODP 983 site (Barker et al., 2015, 2019; Barker and Knorr, 2021). (d) The δ¹⁸O_speleo record of the Sanbao cave in China (Cheng et al., 2012, 2016). (e) δ¹⁸O_planktic record from both the ODP 977 site and the MD01-2444 sediment core plotted using the speleo-tuned marine chronologies. The original chronology of the MD01-2444 core is in light green (Hodell et al., 2013).

The peaks of maximum freshening for the three studied glacial terminations in RAT occurred during the falling precession index and are, therefore, consistent with the beginning of the summer insolation increases (Fig. 4b). In contrast, the obliquity phase is very variable among the three terminations. The TIV and TIII.a occurred during the obliquity increase, while this parameter was decreased through TIII. Hence, falling precession is the common orbital feature along the studied terminations, supporting the idea that precession appears more important than obliquity in predicting the onset of Late Pleistocene glacial terminations (Hobart et al., 2023). Tzedakis et al. (2017) proposed that the energy needed to force the glacial terminations also depends on the length and the intensity of the previous glacial period. The evaluation of the glacial conditions also reveals differences among the three studied terminations. Robust glacial conditions properly preconditioned the continental ice sheet volume before the TIV deglaciation. Lower precession amplitude, together with minimum obliquity through the full glacial conditions of the MIS 10, led to less caloric summer insolation before the TIV onset, strengthening the glacial maximum conditions. The caloric summer insolation is controlled by the combination of precession and obliquity and represents the amount of energy integrated over the caloric summer half of the year, certainly playing a causal role in the achievement of interglacial acmes (Milankovitch, 1941; Tzedakis et al., 2006). But it should also influence glacial attributes. MIS 10 culminated with a maximum in the ice sheet extent according to the benthic δ¹⁸O record from the LR04 stack (Lisiecki and Raymo, 2005), but this was also observable in the more extreme positive glacial δ¹⁸O in the δ¹⁸O RAT record. Furthermore, IRD was released continuously in the North Atlantic at the end of MIS 10, highlighting freshwater inputs (Fig. 4c) (Barker et al., 2015, 2019; Barker and Knorr, 2021). These intense pre-deglacial conditions facilitated the early freshening phase at the onset of TIV and may have reinforced the duration of the subsequent pronounced freshening and the relaxation time observed in the δ¹⁸O RAT record. The maximum glacial conditions with noticeable IRD release did not present such a prominent ice volume maximum at the end of the MIS 8 in the benthic δ¹⁸O or in the RAT records. This glacial situation favored the shorter and less noticeable freshening of the TIII (Fig. 4c). The δ¹⁸O RAT shows a similar freshening structure between TIII and TIII.a, despite TIII.a being an exceptional deglaciation because it does not proceed from full glacial conditions. Short and cold stadial conditions, with less IRD release and accumulated ice volume, allowed a strong response of the ice sheet dynamics to the insolation rise; this was probably due to the pronounced orbital changes through this deglaciation. Interestingly, the meltwater release of TIII.a culminated significantly before the summer caloric insolation maxima, reinforcing the exceptional character of this extra termination (Fig. 4c), while TIV and TIII are closer to their respective summer caloric insolation maxima. The extent of the studied freshening is in line with AMOC weakening, according to the Chinese Sanbao speleothems (Fig. 4d).

Figure 5 Comparison of the δ¹⁸O_speleo records from the RAT speleothem and the north Iberian speleothems (Stoll et al., 2022) along the last five glacial terminations. In the left panel, the speleothem records are plotted against the obliquity and precession parameters (Laskar et al., 2004). The onset of each interglacial acme is highlighted with a vertical pink line according to the Sanbao speleothem record (Cheng et al., 2012, 2016), which considered the recovery of the East Asian summer monsoon intensities as the end of the deglaciation. The right panel shows the duration and the magnitude of the δ¹⁸O freshening structures in Iberia using differential values on the x and y axes. The zero value corresponds to the first and most negative δ¹⁸O value of the main freshening structure of each termination. The black arrow underlines the freshening relaxation of the negative isotopic anomaly after the maximum freshening rate. The maximum summer caloric insolation after each glacial termination is also highlighted with a red circle. It is notable how different the positions of the summer caloric insolation maximum and precession minima are with respect to the freshening structure and the interglacial acme among the compared terminations.

A comparison among the last five glacial terminations is presented in Fig. 5. The studied terminations (TIV, TIII, and TIII.a) are directly compared with terminations TII and TI documented in δ¹⁸O records of Cantabrian speleothems from the north Iberia Peninsula caves (Stoll et al., 2015, 2022). The TII and TI also occurred in line with the falling precession index. In the right panel, all the records were centered at the first and most prominent δ¹⁸O negative anomaly structure in order to compare their duration. The same exercise was performed with the isotopic measurements to observe the freshening magnitudes. This comparison underlines the fact that TIII and TIII.a were the ones with less freshening magnitude, while TIV was prominent. Unlike the benthic δ¹⁸O records, which feature a freshening and plateau at interglacial values, in the speleothem records the peak freshening marked by the negative δ¹⁸O relaxes back to intermediate values. During TII, the period of the most extreme negative δ¹⁸O in north Iberia speleothems indicates fresher surface ocean δ¹⁸O_sw, highlighting the concentration of meltwater in the North Atlantic, and the relaxation has been interpreted to reflect the mixing of this meltwater anomaly throughout the global ocean as the rate of meltwater addition slows and the rate of overturning circulation strengthens (Stoll et al., 2022). The RAT record suggests differences in the time required from the deglacial onset until the relaxation of these negative δ¹⁸O anomalies from very rapid (∼ 5 kyr) during TIII.a, similar to the TI relaxation before the Bölling–Alleröd event; and intermediate (∼ 10 kyr) for TIII, similar to TII; to much slower (∼ 20 kyr) for TIV (Fig. 5).

3.4 Implications for marine chronologies

The RAT speleothem provides absolute chronologies for the deglacial time sequences in the western Mediterranean, which can be transferred to marine sequences that have proxies which are also sensitive to deglacial freshening. This is the case of the δ¹⁸O_planktic data measured in G. bulloides from the ODP 977 site in the Alboran Sea, where their combined use with Mg $/$ Ca-derived temperatures has provided δ¹⁸O_sw reconstructions that prove the arrival of North Atlantic freshening associated with the major melting phases of mostly the Fennoscandian ice sheets during TV, TIII.a, TII, and TI (Azibeiro et al., 2021; Gonzalez-Mora et al., 2008; Jiménez-Amat and Zahn, 2015; Torner et al., 2019; Stoll et al., 2022). These studies conducted at the ODP 977 site indicate that the timing and magnitude of the main negative δ¹⁸O_sw shift are synchronous with δ¹⁸O_planktic changes at glacial terminations, suggesting that the freshwater signal may have a dominant influence on the δ¹⁸O_planktic variability (Fig. A4). Moreover, the amplitude between the deglacial anomalies of the δ¹⁸O_planktic record in the ODP 977 and the δ¹⁸O RAT records is similar for the three studied terminations. Hence, the new δ¹⁸O_planktic data for TIV, TIII, and TIII.a provide an exceptional base to tune the millennial-scale features of the marine record to our dated δ¹⁸O RAT record, achieving an orbitally independent chronology. Therefore, the more negative δ¹⁸O_planktic structures have been used to align the chronology with the δ¹⁸O negative anomalies of RAT during the glacial terminations. The resulting ODP 977 age model is shown in Fig. A5, where the tie points used by the alignment can be observed. Freshening due to ice melting is notable in the δ¹⁸O_planktic ODP 977 signature, and the resemblance with the δ¹⁸O RAT record is remarkable during glacial terminations (Fig. 4e). This new ODP 977 chronology is further applied to records from the Iberian Atlantic Margin (MD01-2444), which also contains deglacial negative δ¹⁸O_planktic anomalies (Fig. 4e) (Tzedakis et al., 2018; Hodell et al., 2013). Several tie points have also been used to tune the MD01-2444 in the base of the most relevant structures (Fig. A5). All the correlations among the tuned records were estimated by linear interpolation using the AnalySeries Software package (Paillard et al., 1996). Particularly, the adjusted TIV exhibits a temporal discrepancy of ∼ 7.4 kyr compared to the established chronology, and the mean age uncertainty in the RAT speleothem throughout this termination is less than this difference, which is approximately ±3.4 kyr. In contrast, the TIII aligns well with previously published chronologies due to a minor time discrepancy, which is approximately ∼ 3 kyr and similar to the RAT average age uncertainty of ±1.5 kyr. In the case of TIII.a, the chronological discrepancy is even more significant, reaching ∼ 11 kyr, with average errors in the RAT estimated to be around ±1 kyr. Hence, this exercise reveals significant time discrepancies, particularly for TIV and TIII.a, and exposes that marine chronologies remain unconstrained – mostly for glacial termination onsets. This study supports the idea that millennial-scale variability could strongly influence climate variation along deglaciations (Barker and Knorr, 2021) and highlights that the distinctive characteristic of each glacial termination should be seriously considered to tune the marine chronologies. In addition, it exposes the power of the RAT speleothem as an excellent example to evaluate the time sequences of TIV and TIII.a and encourages the revision of the marine chronologies to clarify glacial termination dynamics independently of the orbital bond.

4 Conclusions

Geochemical proxies in the speleothem RAT reveal alternating dry glacial conditions and wet interglacial conditions in the western Mediterranean region from MIS 11 to MIS 7. Deglaciations were systematically initiated with major ocean freshening events in the North Atlantic and recorded as negative δ¹⁸O anomalies in the speleothem in parallel with extreme arid conditions in the western Mediterranean region according to the Mg $/$ Ca. This study provides further evidence that falling precession was the major orbital trigger for glacial terminations and constrained the deglacial time sequences of each studied glacial termination. Despite the common features, the three deglaciations present differences in terms of magnitude, intensity, and duration. The dissimilarities along the preceding glacial conditions also contributed to these distinctive features. TIV was slower, longer, and more intense than TIII and TIII.a. The RAT record provides an absolute chronology for the intense meltwater release events that initiated the glacial terminations but also for the maximum freshening of each deglaciation, namely TIV onset at 351 ± 5 ka and the maximum freshening peak at 343 ± 3 ka, the TIII onset at 251 ± 3 ka and the maximum freshening peak at 245 ± 2 ka, and finally the TIII.a onset at 229 ± 2 ka and the maximum freshening peak at 227 ± 1 ka. The newly achieved chronology indicates that revision is needed to the most common marine chronologies, mostly for TIV and TIII.a, which are based on orbital tuning. These discrepancies and the role of millennial-scale variability associated with the onset of glacial terminations suggest that a standard orbital tuning of marine sequences will often not yield a robust chronology. The RAT speleothem provides for the very first time the possibility of validating and improving the accuracy of North Atlantic and western Mediterranean marine chronologies by independent absolute dating.

Appendix A: U $/$ Th laboratory procedure

Speleothem carbonate sample powders were dissolved in 7N HNO₃ and spiked with an in-house solution with known ²³³U, ²³⁶U, and ²²⁹Th concentrations. After that, a few drops of concentrated HClO₄ were added to attack and remove organic matter. After 1 h of refluxing, the samples were completely dried and re-dissolved in 2N HCl. The U and Th were isolated from Ca and several trace elements by means of an iron co-precipitation. Then, one–two drops of FeCl₂ were added, and several NH₄OH drops were added later to induce the Fe precipitation, together with U and Th. Afterward, samples were centrifuged and rinsed with super-clean water three times, while the overlying liquid was removed. In order to separate Fe from U and Th, iron samples were dissolved and loaded in anion exchange resin columns. While Fe was removed with 7N HNO₃, Th was collected with 6N HCl and separated from U, which was recovered ultimately with super-clean water. Each pair of samples was dried and dissolved in a low-concentrated nitric acid and analyzed by means of multi-collector inductively coupled plasma mass spectrometry at the University of Minnesota and Xi'an Jiaotong University. Results and dating errors in the RAT and VALL 2 speleothems are found in Table A1.

Table A1 U $/$ Th dates from RAT and VALL 2 speleothems.​​​​​​​

U decay constants: ​​​​​​​λ238 = 1.55125 × 10⁻¹⁰ (Jaffey et al., 1971) and λ234 = 2.82206 × 10⁻⁶ (Cheng et al., 2013). Th decay constant: λ230 = 9.1705 × 10⁻⁶ (Cheng et al., 2013). ^a ​​​​​​​δ²³⁴U = ([²³⁴U $/$ ²³⁸U]_activity − 1) × 1000. ^b ​​​​​​​δ²³⁴U_initial was calculated based on ²³⁰Th age (T), i.e., ​​​​​​​δ²³⁴U_initial = ​​​​​​​δ²³⁴U_measured × e^λ234×T​​​​​​​. Corrected ²³⁰Th ages assume the initial ²³⁰Th $/$ ²³²Th atomic ratio of 4.4 ± 2.2 × 10⁻⁶. Those are the values for a material at secular equilibrium, with the bulk Earth ²³²Th $/$ ²³⁸U value of 3.8. The errors are arbitrarily assumed to be 50 %. ^c BP stands for “Before Present”, where the “Present” is defined as the year 1950 AD.

Figure A1 VALL 2 speleothem age model acquired by StalAge (dark blue line) and the original U $/$ Th dates (dark blue dots), plotted with the associated errors (error bars and green area). The VALL 2 growth rates are shown in dark green. On the right, there is the speleothem image with the position of the sampled radiometric dates used for the age model and the growth axis path used for the geochemical analyses (white line). In the lower panel the δ¹⁸O, δ¹³C, and Mg $/$ Ca results with reversal axes. The pink bar highlights a slow growth rate and the dashed grey line points to major color and growth axis direction changes in the speleothem.

Figure A2 RAT Mg $/$ Ca ratios compared with Sr $/$ Ca and Ba $/$ Ca ratios in counts per second (right axes) and the calculated mmol mol⁻¹ and µmol mol⁻¹ values (left axes).

Figure A3 Transmission of the seawater isotopic signal from the North Atlantic to the Mediterranean Sea along the TI modified from Stoll et al. (2022). The δ¹⁸O_sw values from paired Mg $/$ Ca and δ¹⁸O_plank decline along the deglaciation, with two well-recognized major freshening steps, where one is abrupt at around 16 ka, and the other more gradual from about 12 to 10 ka (vertical blue bars). These freshening events were recorded synchronously and with a similar magnitude in different sediment records from three different locations. The core MD99-2334K from the western Iberian Margin in the North Atlantic (Skinner and Shackleton, 2006) and two sediment cores from the Alboran Sea in the western Mediterranean (Català et al., 2019), with both close to the ODP 977 site, as well as the MD99-2343 core from the Balearic Sea (Català et al., 2019). Sea level rise from Lambeck et al. (2014). In the bottom panel is the gradient of the isotopic signal between the Alboran Sea or the Balearic Sea against the west. The Iberian record shows the characteristic enrichment of the seawater throughout the Mediterranean Sea, with probably higher evaporative processes in the Menorca rise leading to even heavier values in this basin.

Figure A4 Comparison between the G. bulloides δ¹⁸O_plank of the ODP 977 site and the δ¹⁸O_sw signal obtained by removing the temperature effect using the derived SST Mg $/$ Ca also measured in G. bulloides of the same sediment core along glacial termination II (Torner et al., 2019). The maximum rate of freshening along deglaciation coincides in both records as rapid negative δ¹⁸O anomalies and is highlighted by a blue bar.

Figure A5 Marine sediment tuning strategy using the RAT age model. At the top, blue circles and diamonds indicate the tie points used for chronological alignment with the speleothem errors associated according to the RAT δ¹⁸O record. The ODP 977 δ¹⁸O record on the RAT speleothem chronology is in the middle. At the bottom are the MD01-2444 planktic and benthic δ¹⁸O records on the RAT speleothem chronology and with the original chronology in light color (Tzedakis et al., 2018).