The Use of Submerged Speleothems for Sea Level Studies in the Mediterranean Sea: A New Perspective Using Glacial Isostatic Adjustment (GIA)

: The investigation of submerged speleothems for sea level studies has made significant contributions to the understanding of the global and regional sea level variations during the Middle and Late Quaternary. This has especially been the case for the Mediterranean Sea, where more than 300 submerged speleothems sampled in 32 caves have been analysed so far. Here, we present a comprehensive review of the results obtained from the study of submerged speleothems since 1978. The studied speleothems cover the last 1.4 Myr and are mainly focused on Marine Isotope Stages (MIS) 1, 2, 3, 5.1, 5.3, 5.5, 7.1, 7.2, 7.3, and 7.5. The results reveal that submerged speleothems represent extraordinary archives providing accurate information on former sea level changes. New results from a stalagmite collected at Palinuro (Campania, Italy) and characterized by marine overgrowth are also reported. The measured elevations of speleothems are affected by the local response to glacial and hydro-isostatic adjustment (GIA), and thus might significantly deviate from the global eustatic signal. A comparison of the ages and altitude values of the Mediterranean speleothems and flowstone from the Bahamas with local GIA provides a new scenario for MIS 5 and 7 sea level reconstructions.

The importance of the sea level in controlling groundwater elevation in both inland and coastal karst was first proposed by [7], but it was only in the early 1970s that the potential of using speleothems to reconstruct former sea levels was fully recognized. For example, [8] discusses Late Pleistocene sea level fluctuations based on radiocarbon ( 14 C) and U-Th ages from a submerged stalagmite retrieved from Ben's Hole cave on Grand Bahama Island, and [9] interprets the origin of the encrustations from Cova de Sa Bassa Blanca cave (Mallorca) as carbonate formations related to the Pleistocene sea level variations.
Later on, other studies used submerged speleothems to constrain relative sea level changes, taking advantage of the improvement of the U-series disequilibrium dating method by mass spectrometry for precise and accurate age determination, e.g., [1,10,11].
These carbonate deposits are usually found in caves that are below the sea level and presently flooded and they often contain marine layers/encrustations or growth hiatuses that identify periods of highstands, when the sea level reached the speleothem.
Coastal caves can be divided into two major types: Sea caves and flooded caves [12]. The first type includes caves formed by marine processes, such as wave abrasion, while the second type comprises caves that developed under subaerial processes acting in dry conditions and were then submerged by the Holocene marine transgression. They can be actively forming along present-day coastlines or occur as relict sea caves on former coastlines. The formation of sea caves is mainly controlled by geological weakness, such as bedding planes, joints, and faults along sea cliffs [13]. Ref. [14] recognized a third group, called flank margin caves. These are solutional landforms whose genesis is due to water mixing in sealed chambers. Sea caves can develop in different salinity environments. The most common coastal caves are anchialine caves that are defined as having a fresh to slightly brackish water lens overlying seawater [15].
The investigation of submerged speleothems from coastal caves has provided considerable scientific data that have contributed to precisely constraining the relative sea level variation curves. In addition to speleothems, the main markers used to reconstruct past sea level variations are shallow-water corals, ice, lagoonal and marine cores, and archaeological remains. The Mediterranean Sea is devoid of any reef-building tropical corals from after the Miocene. However, few studies have used the endemic coral species Cladocora caespitosa as a sea level marker, e.g., [16]. This coral species forms 3D structures comparable to tropical reefs [17] and shows a fossil record dating back to the Late Pliocene [18]. Despite the great potential of this temperate coral to reconstruct sea level variations in the Mediterranean Sea, the major drawback is the large depth range where this species lives, ranging from a few meters to around 40 m, which increases the uncertainty in sea level reconstructions.
When a rising sea level inundates a coastal cave, speleothems stop growing once submerged and several events can happen:
Among the most important findings of speleothems retrieved from flooded marine caves worldwide and often showing hiatuses or marine overgrowth, special attention deserves to be paid to the submerged cave deposits from the Bahamas [1,10,30], Bermuda [31], the Yucatan Peninsula [32,33], and Cuba [34]. Ref. [35] provides an exhaustive review on the use of sea level indicators preserved in coastal caves, including speleothems and cave sedimentary deposits.
The measured elevations of speleothems are affected by the local response to glacial and hydro-isostatic adjustment (GIA), and thus might significantly deviate from the global eustatic signal. The contribution of GIA to local relative sea level (RSL) changes is primarily a function of the distance with respect to the ice sheets. The closer one is to the ice sheets, the larger the deviation from the global mean (i.e., the eustatic sea level change). In near-field, ice-proximal settings, the RSL signal is opposite in sign (i.e., trend) and can be 1-2 orders of magnitude larger than the eustatic signal. In these areas, the ice sheet configurations (i.e., areal ice thickness variation) and the spatial patterns of retreat and advance are the dominant factors. The Mediterranean Sea is at an intermediate distance with respect to the ice loading centers/margins, and the northern coastlines are likely to be directly affected by the ice-loading term of the Eurasian ice sheets. However, GIA is strongly modulated by the solid Earth rheological behavior, which can cause large deviations from the eustatic signal, even in far-field regions. Here, in fact, the relevant Earth rheological parameters can enhance the solid Earth response to the water-loading term (counterpart of the ice-loading term). Therefore, the shape and size of the ocean basins become important.
Previous studies have shown that GIA-related vertical motions might have affected the elevations of speleothems [5,36,37]. Different sites (the Bahamas, Bermuda, Mexico, Cuba, and the Mediterranean Sea) show distinct relative sea level responses to glacio-and hydro-isostatic adjustment as a function of the relative distances to the continental ice sheets during MIS 2.
However, the maximum vertical variability of the MIS 5.5 maximum sea level highstand predicted by GIA models in the Mediterranean Sea [36,38] is 2-3 m. Models were tested on a field by MIS 5.5 coastal fossil deposits and fossil tidal notches. This variability largely depends on the mantle viscosity profile and the retreat pattern of the MIS 6 ice sheets. So far, no GIA models have been published for MIS 7 or MIS 5.1 and 5.3 in the Mediterranean Sea. GIA results for MIS 5.1 and 5.3 along a transect between Barbados and Northern Florida have been published by [39]. The Bahamas flowstone is in a far-field position w.r.t. north American ice sheets and should be comparable to the Mediterranean sites (see Figure A2).
The number of speleothem studies for sea level reconstructions significantly increased after the 80s, when the development of mass spectrometric techniques for the measurement of uranium and thorium isotopes led to increases in the typical precision and accuracy of U-series ages compared to alpha-counting measurements and reduced sample size requirements [40,41]. Speleothems can usually provide more reliable records than, for instance, shallow-water corals, given that their dense calcitic or aragonitic structure is less susceptible to post-depositional alteration, significantly reducing the isotope exchange with the surrounding environment.
Isotope-dilution mass spectrometry for submerged speleothem dating was first applied by [1,11] on a flowstone (DWBAH) recovered 15 m below the present sea level of Grand Bahama Island, which contains five hiatuses ( Figure A3), and provided a record of the sea level over the past 280 kyr. These studies reconstructed one of the first sea level change curves which, together with the pioneering works by [42,43] on fossil corals, has considerably contributed to knowledge in the field of sea level change.
Similarly, the use of thermal ionization mass spectrometry (TIMS) and multi-collector inductively coupled plasma mass spectrometry (MC-ICPMS) for the precise U-Th dating of Mediterranean speleothems has greatly improved our understanding of former sea level changes.
In addition to U-Th dating, vadose stalagmites and stalactites and POS can also be dated by 14 C. However, the radiocarbon ages can be affected by the "dead carbon proportion" (DCP), which is the percentage of dead carbon incorporated in the speleothem or POS at the time of formation. The DCP is mainly the fraction of carbon within the calcium carbonate that is derived from equilibration with the soil CO2 and the radiocarbon-free "dead carbon" from the bedrock [44]. High DCP values (e.g., >20%), resulting from a stronger water-soil-rock interaction of the precipitating water, overestimate the radiocarbon ages and cause a large 14 C-U-Th age offset. In a recent study, [45] reported combined 14 C and U-Th measurements on vadose stalactites and POS from Mallorca and showed that the use of radiocarbon dating can be problematic when DCP values are high. Here, we provide a thorough review on the use of Mediterranean submerged speleothems to reconstruct past sea level variations, reporting literature and new data from samples retrieved in the western and central Mediterranean Sea, including the Balearic Islands, Croatia, Malta, and Italy ( Figure 1).

Figure 1.
Overview map of the Mediterranean region showing the locations of the coastal caves discussed in this study. Specific information (latitude, longitude, sampling depths, dating method, and ages) is reported in Table 3 and 4.

General Overview of the Mediterranean Region
The Mediterranean Sea (MS) is a marine basin largely bordered by land that covers an area of about 2.5 million km 2 . Its geographical features have been comprehensively reviewed by [46]. The coastline extends for about 46,000 km and approximately half of this is rocky, with plunging cliffs, sloping coasts, screes, and shore platforms [47]. The Mediterranean Sea can be divided into the Western and Eastern Mediterranean along an imaginary axis between the Strait of Sicily and Tunisia, but a large amount of scientific literature also identifies it as an undefined "central Mediterranean area" [48][49][50].
The MS is situated at the boundary between the subtropical and mid-latitude zones and its climate is characterized by warm and dry summers and mild and rainy winters [51]. It is a semi-enclosed and highly evaporative basin that is connected to the Atlantic Ocean through the Strait of Gibraltar (at a Camarilla sill depth of 284 m) and with the Black Sea through the Strait of Dardanelles (sill depth of ~65 m) and the Bosphorus Strait (sill depth of ~42 m). The Atlantic water that enters the Mediterranean Sea spreads throughout the entire basin and participates in the formation of intermediate and deep waters that contribute to the Mediterranean thermohaline circulation [52].
Tides vary from place to place along the coasts of the Mediterranean, depending on several parameters, such as the coastal geometry and bathymetry. In general, the Mediterranean is characterized by tides with lower amplitudes compared to other oceanic basins. The average tidal amplitude is about 40 cm, with the exception of the remarkably large tides observed in the Gulf of Gabes (Tunisia) and in parts of the North Adriatic Sea, where they may reach amplitudes of up to 1.80 m. In other areas, such as in Greece or Sicily, tides are very small, especially near the amphidromic points, where the tidal range is almost non-existent. At the Strait of Gibraltar, the tide increases to 1.5 m due to the influence of the Atlantic Ocean, but it quickly decreases eastwards. Weather conditions can significantly reduce or amplify the tidal amplitude, with variations of up to 1 m.

Numerical Modeling of Glacial-Isostatic Adjustment (GIA)
We can compute the ice-driven, GIA-modulated RSL curves by solving the gravitationally self-consistent Sea Level Equation, hereafter called SLE [38,53,54]. The SLE is a linear equation that yields, for any variation of land-based (i.e., grounded) ice thickness, the time-and space-dependent vertical deformation of the mean sea surface and solid surface of the Earth, with the difference between the two being the relative sea level (RSL) change.
By neglecting gravity and solid Earth deformations, and by assuming that the continents are bounded by infinite vertical walls, the solution of the SLE for a prescribed ice mass variation yields the eustatic sea level change, i.e., the change of ice sheet volume (or mass) as meters of equivalent sea level. This solution is straightforward and directly depends on the ice mass change and conversely, on the (constant) ocean area. Furthermore, it does not depend on the geographical location, implying that each location in ocean areas experiences the same amount and sign of sea level change. Accordingly, the local RSL changes are equal to the global eustatic sea level change when deformations and gravity are neglected.
By including gravity, the SLE becomes gravitationally self-consistent and its solution is spatially varying, mostly depending on the distance from the changing ice masses [55]. Self-gravitation, in fact, allows for deviations from the eustatic signal that increase towards the changing ice masses, where the sea level change switches sign with respect to the eustatic signal. Therefore, any ice mass increase, or drop, would result, respectively, in a local sea level rise, or drop, in the proximity of an ice sheet.
Any net surface ice mass gain, or loss, is instantaneously accompanied, respectively, by an increase, or decrease, of the mutual gravitational pull between ice and ocean water mass. The latter is free to move and its surface is an equipotential surface of gravity: The gravity vector is always perpendicular to the mean sea surface. Accordingly, an increase in the ice mass acts as a positive density anomaly, which in turn deflects the gravity vector towards the ice mass, causing a local sea level rise. Because mass must be conserved, the sea level drop in the far field is larger than the eustatic signal, and the average over the oceans of the sea level change is equal to the eustatic signal. Therefore, the local RSL changes are equal to the local changes of the mean sea surface.
The complexity and realism of the SLE increase by including solid Earth deformations in response to surface ice-and water-loading changes. Accordingly, the RSL changes also stem from the vertical motions of the Earth's solid surface and therefore depend on the relevant rheological parameters of Earth. Earth exhibits a dual rheological behavior with respect to surface load changes: (i) Solid elastic deformation (solid behavior), which is instantaneous (i.e., time-independent) and a linear function of the surface load, and (ii) viscous flow (fluid behavior), which linearly depends on the surface load, on the viscosity, and therefore on the amount of time since the activation of the load (i.e., load change). Therefore, a commonly used model for GIA simulations is the linear Maxwell viscoelastic model, where an elastic component (spring) is combined in series with a viscous element (dashpot). This allows for the calculation of both the immediate (elastic) and delayed (viscous) deformations that decay exponentially over time as a function of viscosity and provided that gravity is the restoring force that tends to restore the isostatic equilibrium. Solid Earth deformations increase the regional variability of RSL changes, but the average over ocean areas remains equal to the eustatic signal (mass conservation).
Solving the gravitationally self-consistent SLE for a deformable Earth model that is characterized by a linear Maxwell viscoelastic rheology implies that, for each point on the surface of the Earth, the predicted RSL change at a specific time depends on all the changes in ice and water thickness (ice and water loading terms) that have occurred everywhere and since the beginning of the chronology. Because the variations of water thickness are the actual RSL changes that one aims to compute, the SLE is integral in its nature. Therefore, a solution can only be obtained through an iterative procedure. Initially, the SLE is solved under eustatic approximation and the uniform but time-dependent water loading term is computed over time. This is then loaded over the ocean areas to retrieve the combined hydro-and glacio-isostatic contributions to water height variations (RSL changes) over the ocean areas. Again, this is loaded on the ocean areas during the next iterative step. The solution converges rather quickly and three to four iterative steps are usually required for the gravitationally self-consistent solution. The latter complies with the requirements of mass conservation and equal gravitational potential for the sea surface.
Here we use the latest, up-to-date version of the open source model SELEN [38,54], which solves the SLE by means of the pseudo-spectral method. The latter implies a discretization of all the variables through spherical harmonics expansion. We include the rotational feedbacks, as well as the form treatment of the time-dependent ocean function (i.e., variable coastlines).
The first component required by SELEN is the Earth rheological model. We assume a self-gravitating, rotating, spherically symmetric, radially stratified, and deformable but not compressible Earth model. The latter is one-dimensional and all the relevant rheological parameters only depend on the Earth's radius. Accordingly, lateral heterogeneities are not accounted for. The outer shell of the model is perfectly elastic and represents the lithosphere. Between the lithosphere and the inviscid core is the mantle, which is characterized by linear Maxwell viscoelastic rheology. The mantle can be discretized into a number of layers.
First, we discretize the VM2 mantle viscosity profile [56] over a four-layer model, where the upper mantle (UM), the lower upper mantle (LUM), the transition zone (TZ), and the lower mantle (LM) are characterized by a Maxwell viscoelastic rheology and laterally uniform viscosity (see Table 1). We combine this with a lithosphere thickness of 90 km. We use this as a reference model to be combined with three ice sheet chronologies. We also employ a three-layer model (UM, TZ, and LM) and test three mantle viscosity profiles (MVP 1-3), each combined with a lithosphere thickness of 100 km ( Table 2). The viscosity gradient of viscosity as a function of the Earth's radius increases by one order of magnitude from MVP1 to MVP3. Here, we focus on the vertical gradient of viscosity and on the role of the LM viscosity, rather than on the role of the lithosphere thickness. In fact, it was shown by [57] that the LM viscosity is the prime mover of the GIA-modulated RSL changes in the Mediterranean Sea, mostly because of the extent of the water load term. In fact, the latter is much larger than the thickness of the lithosphere, and therefore, is sensitive to the viscosity values at deeper positions in the mantle.
We employ the following ice sheet models:  ICE-5G [56]: This global model describes the ice-sheet thickness variations over North America, Eurasia, Greenland, and Antarctica during the last 123 ka. The ice sheet evolution between 123 and 26 ka is tuned to the  18 O curve by [58] and lacks any glaciological and geological constraints. Between 26 ka and the present day, ICE-5G is constrained by glaciological, geological, and modern geodetical observations from remote sensing [56]. Here, we combine two ICE-5G chronologies over time in order to cover the time span under consideration and evaluate the GIA-modulated RSL response to the melting of the MIS 6 ice sheets during MIS 5e. The latter is characterized by a ~0.9 m eustatic sea level highstand above present-day msl between 129.5 and 123 ka, mostly coming from a reduced (with respect to the present day) Greenland Ice Sheet;  ICE-6G [59]: This model follows from an improvement of the ICE-5G. Here, we apply the same repetition in a series of ICE-5G. According to the ICE-6G, the maximum MIS 5.5 eustatic highstand between 129 and 123 ka is ~3.1 m above present-day msl as a consequence of reduced Greenland and West Antarctic Ice Sheets;  ANICE-SELEN [38,60,61]: This global chronology model is the result of an inverse forward modeling procedure, where the  18 O stack [58] is decoupled into the global ice sheet volume and deep-water temperature. For this purpose, 3D thermomechanical ice-sheet models for North America, Eurasia, Greenland, and Antarctica are dynamically coupled to SELEN, in order to include all the GIA feedback. ANICE-SELEN is not constrained by geological or instrumental data. Here, we use the original ANICE-SELEN chronology, which is characterized by a ~2.5 m eustatic highstand during MIS 5 e (125-117 ka), in response to a reduced Greenland Ice Sheet volume.

Speleothem Samples
The Mediterranean speleothems reported in this study were collected in coastal caves of the Balearic Islands and along the Croatian, Maltese, and Italian coasts ( Figure 1; Table  1). The Italian samples were retrieved between 1990 and 2007, at water depths ranging between −0.3 and −48 m. The depth was determined with a digital depth gauge (typical error of ±0.1 m). The depth has not been corrected for tidal data, which means that the elevation error of the studied speleothems is at least ±30 cm (i.e., the sum of the instrument uncertainty and the Mediterranean tidal mean of ±20 cm [50]). It is worth noting that the error associated with the vertical position of the Mediterranean submerged speleothems is lower than for the samples retrieved from other oceanic basins due to the microtidal regime of the Mediterranean. Samples were measured and collected by sawing through the base of the speleothems. For further information regarding the cave morphology, sampling method, date of sampling, and error associated with the radiocarbon and U-Th ages, the reader should refer to the original articles [1][2][3][4][5][6]11,19,[21][22][23][24][25][26]28,29,[62][63][64][65][66].

New Stalactite from Palinuro, Italy
A ~10 cm stalagmite was collected in a fossil tidal notch along the carbonate coast of the Palinuro promontory in the Campania region (South Italy) at 2.1 m above the sea level. A variety of living organisms, including the crustose coralline algae Lithophyllum, the acorn barnacle Chthamalus, and the limpet Patella, were observed below the tidal notch, exposed to strong wave conditions. The external surface of the speleothem was completely covered by a blackish patina that is likely the result of the precipitation of Fe-Mn oxides that occurred when the rising sea level approached the speleothem (Figure 2d,g). The top portion of the stalagmite is encrusted by a marine overgrowth that consists of several specimens of the barnacle Chthamalus stellatus, which typically lives in the intertidal zone on high energy rocky shores and can be found above the highest tidal level [67].

14 C and U-Series Dating of the Stalactite from Palinuro
Three carbonate fragments were sub-sampled from the upper portion of the stalagmite from Palinuro (Figure 2d,g,h) and prepared for radiometric dating. In particular, subsample 1 (Palinuro-1) was collected within the marine overgrowth (barnacle) and analysed for AMS− 14 C. The fragment was carefully cleaned using a small diamond blade to remove any visible speleothem calcite and only retrieve the encrusted overgrowth. The sample was further processed and analysed for 14 C at the Center of Applied Physics, Dating and Diagnostics (CEDAD) Laboratory of the University of Salento, Italy, using an Accelerator Mass Spectrometer (AMS). The radiocarbon age was converted to calendar years (cal. yr BP, BP = AD 1950) using the IntCal20 calibration curve [68] and the Calib 7.10 program [69]. The ages of the other two sub-samples (Palinuro-2 and Palinuro-3) were determined through the U-series dating method at the Laboratoire des Sciences du Climat et de l'Environnement (LSCE) at Gif-sur-Yvette, France. The samples were first leached with 0.1 M HCl and then dissolved with diluted HCl and equilibrated with a mixed 236 U-233 U-229 Th spike. The uranium and thorium fractions were separated using UTEVA resin (Eichrom Technologies, USA) and analysed using a ThermoScientific Neptune Plus multicollector inductively coupled plasma-mass spectrometer, following the protocol developed at LSCE [70]. The 230 Th/U ages were calculated from measured atomic ratios through iterative age estimation [71], using the 230 Th, 234 U, and 238 U decay constants of [72,73]. The ages were corrected for the non-radiogenic detrital 230 Th fraction using an initial 230 Th/ 232 Th activity ratio of 0.85 ± 0.36, which corresponds to the mean upper crust value [74].
The U-Th ages of the DWBAH flowstone reported in [11] have been re-calculated using the most recent 230 Th, 234 U, and 238 U decay constants of [72,73] and corrected for the detrital 230 Th fraction using an initial 230 Th/ 232 Th activity ratio of 0.8 ± 0.8 [74] (Table A1), as suggested by [75] for the speleothems from the Bahamas. The difference between the original U-Th ages [11] and the re-calculated and corrected ages is minimal (<3.5 kyr for ages younger than 270 kyr).

Results
Tables 3 and 4 report the relevant details of the sampling locations and the published U-Th and 14 C ages of the vadose speleothems, POS, and marine biogenic overgrowths used for the present review. In particular, Table 3 summarizes the main results of all the studies dealing with submerged Mediterranean speleothems, which were published between 1992 and 2020: 33 caves (Figure 1, Table 3) carved on carbonate lithology that preserved vadose speleothems, POS, and marine biogenic overgrowth that were used to reconstruct former sea level variations. An essential condition for the conservation of submerged speleothems is that the cave's entrance must be small in size so that the energy of the waves can be almost completely attenuated.  Table A1 −15 Table 4. Location, depth, sample code, U-Th, and radiocarbon ages of the speleothem studied in Croatia (also see Figure  3). 14 C ages of speleothems L-1 and L-10 measured through liquid scintillation counting and speleothems Z-41, B-38, B-36, B-34, B-28, P-23, and R-21 measured by an alpha proportional counter. The youngest parts are marked with an A and S, while B refers to the oldest parts. 14 C ages are expressed as conventional 14 C corrected for A0 = 85% and measured δ 13 C (-8‰ when not measured), and as calibrated ages. From [23,24]. On the Adriatic coasts of Croatia, speleothems from seven submerged caves distributed along ~300 km of coast (Figure 1 and Figure 3a) were sampled and dated by 14 C and U-Th [3,23,66]. Their ages range between 7 and 217 kyr (MIS 1 and MIS 7), and the results are summarized in [25]. Since 1974, a pool of Spanish and international speleologists and geologists have investigated vadose speleothems and POS collected from nine different caves located along the southern coast of Mallorca in the western Mediterranean ( Figure 4). The phreatic speleothems were retrieved at 1.2 and 5.8 m above the present sea level and their ages vary between 1 and 231.9 kyr [4,5,27,28,45,62,76]. The main results are reported and discussed in two important studies: [4] focuses on POS sampled in five different caves at 1.2 to 1.5 m, which precipitated within a very narrow temporal range, from 81.95 ± 0.56 to 80.43 ± 0.48 kyr, whilst [5] reports the U-Th ages of 11 POS from eight caves that enabled the relative sea level changes between 126.6 ± 0.4 and 116 ± 0.8 kyr to be reconstructed. In Italy, the first speleothems sampled for sea level studies were retrieved from 15 submerged caves in the Tyrrhenian Sea (Figures 1 and 5) between 1992 and 1994 [19,63]. Between 2002 and 2017, a couple of studies on speleothems from the Argentarola and Custonaci caves, characterized by an alternation of continental layers, hiatuses, and encrusting marine organisms allowed the sea level change during MIS 7 [2] and the Middle Pleistocene Transition [6] to be constrained.

Croatia (2005-2010)
The Adriatic coast of Croatia is a karstic coast characterized by more than 230 submarine caves that have been investigated over the last two decades to reconstruct former sea level changes ( Figure 3, Table 4). Slightly more than half of the caves (n = 126) show full marine conditions, 75 are anchialine (mostly) pits, 13 are submarine springs (vruljas), and 21 are not determined by hydrological or environmental conditions. More than 140 caves contain speleothems [77].
In particular, 16 submerged speleothems were collected at depths between −41.4 and −1.5 m in seven caves (out of 235 visited) located along the Croatian coast over a distance of ~300 km. The age model of these samples was constrained by 15 radiocarbon measurements and 36 U-Th analyses. Among these, the U-Th ages obtained from the stalagmites K-14 and K-18 sampled in U vode Pit (Krk Island, Figure 2) at a −14.5 and −18.8 water depth deserve particular attention. Eleven U-Th ages from K-18 span a ~50 kyr time range, from 97 kyr (MIS 5.3) to 50 kyr (MIS 3), and identify two hiatuses, whereas 5 U-Th measurements from K-14 cover a ~10 kyr period, from 91 to 82 kyr (MIS 5.1), with the presence of one hiatus marked by calcite, gypsum, and halite precipitation (Figure 3a,b) The deposition of speleothem K-18 ceased during MIS 3 (58-28 kyr BP) and few data have been published for MIS 2 ( Table 2). The radiocarbon ages of the continental layers of the speleothems L-1-S, L-10-S, and Z-41-S have provided some interesting data on the sea level rise during the Holocene (Table 4): −1.5 m at 3.3 kyr, −10 m at 7.8 kyr, and −41.5 m at 9.2 kyr cal BP. The comparison with the sea level curve predicted by [78] for Croatia reveals a good agreement for the first two points (L-1-S and L-10-S), whereas the value for sample Z-41-S slightly deviates from what is expected.

Italy (1992-2017)
Between 1992 and 1994, [19,63] published the first relative sea level curve for the Tyrrhenian Sea during the Holocene using the radiocarbon dating of samples at the contact between serpulid overgrowth and continental stalactites. The samples were retrieved at a 41 and 6 m water depth at the sites of Palinuro and Argentarola and 14 C ages range between 26 and 3.5 kyr cal BP. The marine overgrowth observed on the Argentarola and Scaletta (Palinuro) speleothems started to encrust the speleothem surface when the rising sea level flooded the caves ( Figure 5). The biogenic overgrowth mostly consists of whiteyellow tubes precipitated by the polychaete worm Serpula massiliensis, and remaining bryozoans and sponges [2]. This encrustation plays a significant role because it protects the continental portion of the speleothem from degradation, bioerosion, or dissolution. The speleothems from the Argentarola and Rumena caves in Italy [2,6] represent the only cases in the world where more than one continuous deposition of continental and marine layers was observed in the same speleothem (for subsequent marine and continental highstands).
Argentarola cave (Grosseto, [2,19,21,22,63]): Submerged speleothems were collected from −21.7 and −3.5 m and dated by the 14 C and the U-series disequilibrium method. The marine overgrowth of the highstand older than 40-50 kyr was not dated given the limit of the radiocarbon method when approaching these ages. The results obtained from speleothems of the Argentarola cave have greatly contributed to understanding of the sea level changes and environmental conditions during MIS 7.5, 7.4, 7.3, 7. Plemmirio cave (Siracusa, [64]): Submerged speleothems were collected from −25.5 and −26.6 m. The radiometric ages obtained from the U-series disequilibrium method of the speleothem calcite and 14 C of the encrusting serpulids range from MIS 1 to 5.1. Nettuno cave (Capo Caccia, [81]): A submerged speleothem was sampled at −52 m, representing the deepest speleothem collected in the MS. Refs. [82,83] reported the age of phreatic speleothems collected from the same cave which, compared with oxidation bands and tidal notches in the external portion of the cave [84], confirms the MIS 5.5 age with a local highstand at 4.2 m (Figure 2c). Other speleothem samples were retrieved from Stalattiti Ubriache grotta Zingaro reserve (Trapani), Grotta di Tragara Capri, Risorgenza di Cala Luna Orosei (Nuoro), and Grotta Azzurra in Ustica [19] and radiocarbon dated. The results from samples at a −9 and −0.6 m depth revealed very recent speleothem formation (between 2 and 4 kyr) (Figure 2i, Table 3). All the sites described in Italy are tectonically stable, with the exception of Plemmirio, which is characterized by an uplift rate of 0.3 mm/yr and Rumena cave, which experienced a continuous tectonic uplift of 0.081 mm/yr for the last 1.4 Myr.

Malta
Ref. [65] studied a submerged speleothem collected at a − hem collected at d spelti-onItems":[{"id":"ITEM-1","itemData":{"DOI":"10.1177/095968(Table 1, Figure 2a). Since the cave was mainly formed in a subaerial karst environment, radiocarbon dating of the speleothem with serpulid encrustations (Figure 2f) enabled the sea level during the mid-Holocene, when the cave was fully flooded, to be reconstructed. In particular, the mean radiocarbon age of 7.6 kyr perfectly aligns with the sea level curve predicted by [78] for Malta. Ref. [65] concluded that the Maltese islands were tectonically stable during the mid-Holocene, and this tectonic behavior still persists today.

Mallorca (1974-2020)
The results of more than 30 years of work carried out in 10 different caves of Mallorca, all located along the southern coast of the island (Figure 1), have been published in numerous papers [4,5,27,28,45,62,76]. POS have been extensively studied and radiometrically dated, including both modern samples formed a few cm above and below the water table and fossil remains collected up to 15 m above the present sea level. The results have significantly contributed to sea level reconstruction of MIS 5.1 and 5.5. In particular, [9] first highlighted the great potential of studying POS in the caves of Mallorca and paved the way for further in-depth investigations of these distinct encrustations. Ref. [28] published the U-Th dating results of 15 POS collected between 1.4 and 5.8 m above the present sea level, spanning ages between 71 and >350 kyr; [76] identified a tilting of the southern coast of Mallorca which reconciles the different heights of speleothems dated to MIS 5.1. Ref. [45] analyzed eight POS samples from Cova de Cala Varques, comparing 14 C with U-Th ages.
The results obtained from the two dating methods generally agree, but the use of radiocarbon dating can be problematic when the dead carbon proportion values are high. Refs. [4,5] reported the results of several U-Th measurements conducted on POS samples located at 1.2 to 1.5 m and at 2.2. to 3.2 m above sea level, respectively. The ages obtained by [4] fall within a narrow range, from 81.95 ± 0.56 to 80.09 ± 0.48 (MIS 5.1), whereas the ages calculated by [5] cover MIS 5.5, from 126.6 ± 0.4 to 116 ± 0.8 kyr. Finally, ref. [62] carried out several U-Th analyses on vadose speleothems collected along the eastern coast of Mallorca at 2 and 5.5 m above sea level. The U-Th ages span a large range, from ~440 to ~1500 kyr (Table 3, Figure 4). As also mentioned by the authors, the oldest ages (>650 kyr), which were obtained using the 234 U/ 238 U disequilibrium method, suffer from large uncertainties, mainly associated with the errors of the initial  234 U estimates needed to calculate  234 U ages.

New Results from the Palinuro Stalactite
The U-Th dating of the vadose portion of the speleothem provided ages between 5.49 ± 1.65 and 4.46 ± 3.74 kyr ( Figure 6, Table 5), whereas the radiocarbon measurement of the fragment of a Chthamalus barnacle gave an age of 1.58 kyr cal BP (median probability). The large uncertainties of the U-Th ages, especially for the sample Palinuro-2, derive from the high detrital Th contamination ( 230 Th/ 232 Th < 2).

Predicted RSL Curves.
The predicted RSL curves for ICE-5G, ICE-6G, and ANICE-SELEN and the VM2 mantle viscosity profile (gray curves in Figure 7c) show an increase of variability in the Mediterranean, with respect to the eustatic signal, when moving from the interglacials (i.e., MIS 5.5) through the interstadials (i.e., MIS 5.1-5.3 and MIS 3) and glacials (e.g., Last Glacial Maximum, about 21 kyr cal BP, LGM).
The relatively limited regional RSL variability during the interglacials is characterized by a maximum deviation from the eustatic signal of ~2.5 m for the three ice sheet models (Figure 7a-c), confirming previous findings [36,38,85]. This also supports the results of [39] for the Bahamas, where the RSL gradient as a function of latitude vanishes during the MIS 5.5. Our predicted RSL curves for the Bahamas (red curves in Figure 7a-c) fall within the Mediterranean response, thus confirming that the Bahamas is an intermediate-field site (w.r.t. north American ice sheets), and similar to the northern Mediterranean sites (w.r.t. Fennoscandia and Eurasian ice sheets).
The predicted maximum highstand during the MIS 5.5 in the Bahamas is delayed with respect to the eustatic signal. For both ICE-5G and ICE-6G, the MIS 5.5 peak occurs at the end of the MIS 5.5 and exceeds the eustatic value. For ANICE-SELEN, instead, the delayed maximum peak is lower than the eustatic value.
The predicted RSL variability during the interstadials shows up as deviation from eustasy with the order of ~15 and ~25 m for MIS 5.1-5.3 and MIS 3, respectively. This confirms previous findings for the Bahamas [37,38], and shows that the ice-loading term also causes strong regional RSL gradients in the Mediterranean Sea.
Additionally, in the Bahamas, the deviation from the eustatic signal increases towards the interstadials and glacials [39], implying its dependency on the fluctuations of the north American ice sheets.
In the Bahamas, the predicted RSL curves are above the eustatic value for both ICE-5G and ICE-6G. On the other hand, ANICE-SELEN results are lower than the eustatic curve. The predictions at the Mediterranean sites are distributed around the eustatic values.
The cumulative predictions in the Mediterranean for three ice sheet models and mantle viscosity profiles (MVP1-3) are characterized by a larger variability (gray curves in Figure 7d-f). Again, the variability increases towards the interstadials and glacials, implying that the viscosity gradient is a significant factor in the regional RSL variability in the Mediterranean Sea. The predicted variability in the Bahamas is significantly larger than in the Mediterranean (red, green, and blue curves in Figure 7d-f) and fully bounds the predictions at the Mediterranean sites. Increasing the viscosity gradient, i.e., moving from MVP 1 to MVP 3, results in an upward shift of the RSL curves between MIS 5.5 and LGM. The predicted total RSL glacial-interglacial excursions (glacial-interglacial) for MVP1 are larger than for MVP3, indicating that a larger viscosity gradient (and a larger viscosity value for the lower mantle) results in a significant delay of the solid Earth response. This can be appreciated during the MIS 5.5, when the maximum peak of transgression appears at the end of the interglacial and during the Holocene (i.e., steeper RSL rise curves Figure  8).

Discussion
The use of submerged speleothems has provided remarkable results for the sea level change history during the Holocene and for long-term sea level change reconstruction of the Late Pleistocene, in particular, for MIS 1, 7.1, 7.2, 7.3, and 7.5 highstands (Argentarola stalagmites), the MIS 5.1 highstand (Croatia and Mallorca), the MIS 5.5 highstand (Mallorca), and the MIS 6.5 highstand (Nettuno cave in Italy) (Tables 3 and 4). Unfortunately, marine overgrowth layers encrusting Italian or Croatian speleothems with ages beyond the radiocarbon dating range (i.e., >40-50 kyr) cannot be analysed by the U-series disequilibrium method due to diagenetic overprints that cause the exchange of uranium [86].

Palinuro Speleothem
For the interpretation of the 14 C and U-Th results obtained from the new Palinuro sample, a study of the tidal area was conducted along the promontory of Palinuro over a distance of ~3 km from the sampling site. The survey revealed the presence of a fossil tidal notch, below which we found specimens of the crustose coralline algae Lithophyllum, the acorn barnacle Chthamalus, and the limpet Patella. Based on [78], the sea level was 4.4, 3.3, and 0.73 m lower that the present level at 5.5, 4.5, and 1.6 kyr, respectively (Table 3, Figure  6). Therefore, the barnacles started to colonize the stalagmite at 1.6 kyr, when the rising sea reached a level that exposed the sample to strong wave conditions, favorable for barnacle growth. However, no living barnacles were observed during the sampling of the stalagmite and this might be due to changes of the environmental conditions since the deposition of the fossil barnacles.

Speleothems from Italy
U-Th dating and geochemical analyses of the speleothem samples from the Argentarola cave provided detailed information on the timing, duration, and position of sea level highstands during MIS 7.1, 7.2, 7.3, 7.4, and 7.5 [2,22]. Based on the  18 O results obtained from the serpulid overgrowth on a stalagmite collected at −18.5 m from the same cave (Figure 5f) and the correlation with foraminifera  18 O values from the Mediterranean ODP site 975 [87], ref. [21] identified MIS 5 as the time of marine overgrowth deposition, without any visible interruption from MIS 5.1 to MIS 5.5 (Figure 9). It is worth mentioning that the timing of the MIS 7.3 and 7.2 highstands identified by U-Th dating on the speleothem agrees with the global sea level curve reconstructed by [88]. In the Nettuno cave, the deepest speleothem sampled in the Mediterranean at −52 m was dated by the alpha-counting U-series method (see Table 1 Figure 4), which provided an age of 165 ± 12 kyr, corresponding to MIS 6.5 and in agreement with [88]. Finally, the speleothem Plemm A (Table 3) was collected in the Plemmirio cave [62] at −20.2 m, which was corrected to −35.3 m when considering the tectonic uplift (0.2 mm/yr, [88])

Speleothems from Croatia
The speleothems retrieved along the Croatian coast and studied by [3,24,25,66,89,90] greatly contributed to the sea level reconstruction from 1.5 to 220 kyr. In particular, highprecision U-Th dating of the speleothems K-14 and K-18 sampled at −18.8 and −14.5 m in the U vode Pit cave (Krk Island) provided significant data to constrain the relative sea level curve during the MIS 5.1 highstand. Regarding the presence of two hiatuses in the K-18 stalagmite, on the basis of the lack of XRD analyzes that invoked the presence of a hiatus between 90.8 and 82.9 kyr, and in correspondence with a change in color, we suppose that only the second hiatus was present between 77.7 and 64.5 kyr (Figure 3).
The ~13-18 m difference between the elevation of the two speleothems (−18.8 and −14.5 m) with the ice-volume-equivalent global sea level curve by [91] has been justified by [25], invoking a long-term regional tectonic uplift with an average rate of 0.15-0.25 mm/yr during the last 75-85 kyr. However, the vertical tectonic movements along the Croatian coast are still a matter of debate. For example, [90] suggested a generalized downlift in Istria and uplift along the southern cost of Croatia. Ref. [12] investigated the carbonate coast stretching from Trieste to southern Croatia and showed that the tidal notch is always submerged between 1.8 and 0.5 m below the mean sea level. Furthermore, archaeological studies along the Croatian coast revealed the presence of Roman remains (2 kyr BP), including a pier and fishtanks, at ~1.5-1.8 m below sea level as a result of tectonic subsidence or coseismic downlift [92,93]. Finally, none of the fossil outcrops related to the MIS 5.5 highstand that are usually found at ~+6 m a.s.l. in tectonically stable regions [36,94] have been observed thus far in the north-eastern Adriatic Sea. Therefore, we hypothesize that most of this area is affected by tectonic subsidence (a minimum of 6 m since MIS 5.5), with possible local uplift movements [25]. Accordingly, we decided to correct the elevation of the stalagmites K-14 and K-18 by 4 m, to −10.4 and −14.7 m, respectively, assuming a constant subsidence of 0.048 mm/yr for the last 80 kyr.

Speleothems from Mallorca
A large number of precise U and Th isotope measurements have been performed on submerged speleothems sampled from 10 coastal caves on the island of Mallorca. Data published by [5] agree with previous studies reporting MIS 5.5 fossil deposits in Mallorca [95,96]. These deposits contain the Tyrrhenian Senegalese faunal assemblage that was analysed by the amino-acid dating technique [95] and remainders of the coral Cladocora caespitosa that were dated by U-Th [96].
The elevation of MIS 5.1 POS samples studied by [4] in Mallorca still represents a scientific debate.
The reconstructions obtained by [88,97] indicate that the sea level during MIS 5.1 was 21.2 and 26.7 m below the present sea level, respectively. However, the global sea level curves do not take into account the local tectonics and GIA. Based on the U-Th results of [96] on fossil corals from Mallorca, we can rule out the presence of MIS 5.1 coral deposits on the Island.
In Sardinia, the presence of fossil tidal notches, attributed to MIS 5.5, was reported by [66] in the Orosei Gulf, which is a tectonically stable coastal area in Italy. The fossil tidal notches in this area show a lateral continuity of tens of kilometers, while there is no evidence of younger fossil tidal notches at a lower elevation.
In other regions worldwide, like in the Caribbean Sea, observational data of MIS 5.1 deposits from submerged speleothems seem to disagree with the finding of POS samples located at ~1.4 m a.s.l. in Mallorca and dated to ~81 kyr BP. The U-Th results for the DWBAH flowstone recovered from 15 m below the present sea level on Grand Bahama Island clearly indicate that, during MIS 5.1, the sea level did not rise as high as −15 to −10 m [11].
Discrepancies in the elevations of different sea level markers, ranging from −20 to + 5 m, are found along a latitudinal transect in tectonically stable areas for the coasts of Haiti, Bermuda, the Bahamas, South Carolina, and Florida. However, these apparent inconsistencies are reconciled when taking into account the isostatic response of the Earth to the North American ice sheets melting during the last glacial cycle [39].
Results from the MIS 5.1 POS samples in Mallorca agree with the finding of some deposits from the same age between + 2 and + 1.5 m near Gibraltar [98], but disagree with the following:


All the eustatic curves (in particular, [88,97]);  The U-Th ages reported by [96] for Mallorca;  The speleothem from Plemmirio located at −20. The morphology of the MIS 5.5 tidal notches [36] and the lack of younger tidal notches;  The mushroom-like landforms with tidal notches at −25 m found on Tavolara island (Sardinia) and interpreted by [36] as being MIS 5.1 in age;  The results of glacial and hydro-isostatic adjustment (GIA) of the Mediterranean Sea (see Figures 7 and 8). We decided to compare the vertical position of most of the submerged speleothems described in this review with the DWBAH flowstone [1,11]. This exceptional flowstone contain five hiatuses identified by the deposition of a thin layer of red or yellow mud ( Figure 10 and A3), and was extensively dated by the U-series disequilibrium method, providing precise ages that unraveled, for the first time, the history of sea level variations from ~300 to ~30 kyr. The hiatuses mark the marine highstands as follows: H1 between 231 and 237 kyr (MIS 7.5) Table A1).
Based on the correction for the subsidence, the DWBAH flowstone started to form ~300 kyr ago at −10 m and ceased growing ~39 kyr ago. The flowstone is presently located at −15 m. Similarly, the elevation of the two speleothems sampled in the Plemmirio cave in Sicily (Table 3) was corrected from 20.2 and 35.3 m b.s.l. to 22.7 and 38.2 m b.s.l., considering an uplift rate of ~0.2 mm/yr, as evaluated by [64].
In Croatia, the elevation of the stalagmites K-14 and K-18 was corrected considering a subsidence of ~4 m in 81 kyr. The correction derives from the fact that MIS 5.5 at 6 m a.s.l., or at higher elevations, has never been found in Croatia. This means that MIS 5.5 is below the present-day sea level and 4 m is extrapolated from a subsidence of the coast of at least 6 m in 125 kyr, corresponding to MIS 5.5. The elevation of the other speleothems considered in the present review (i.e., Argentarola cave, Mallorca, Malta, and Sardinia) was not corrected for uplift or subsidence owing to the tectonic stability of these areas.
The distribution of coastal caves in the Mediterranean Sea necessarily reflects the presence of carbonate rocks thick enough to enable the formation of caves, even when the sea level was 60-70 m lower than the present level. Generally, submerged speleothems are rarely preserved in carbonate formations from regions with a strong uplift or subsidence. In the western Mediterranean, research has largely focused on Mallorca, which is a tectonically-stable island with calcarenites that hosts several POS deposits. In Croatia, where the tectonic activity is limited (e.g., small subsidence) and the lithological features are monotonous, the formation of submerged speleothems is very active and preservation is good. In Italy, sea caves predominantly formed in tectonically-stable areas and carved limestones with different ages, textures, and microstructures, including Mesozoic-Paleogene carbonates and Miocene biocalcarenites. Karst in specific sites along the Italian coasts has developed for millions of years, shaping unique caves (e.g., Grotta Rumena, Palinuro, and Argentarola), where speleothems with multiple marine and continental layers have been preserved.

GIA-Modulated RSL Curves
The GIA-modulated regional RSL variability in the Mediterranean is minimal for the interglacials (~2.5 at MIS 5.5), but increases during the interstadials (~15 m at MIS 5.1-5.3), thus complicating the task of finding the correct eustatic value (gray curves in Figure 7). In fact, the GIA signal is strongly dependent on the solid Earth rheological profile, which is itself an unknown. In particular, the viscosity gradient and the lower mantle viscosity are important parameters in the Mediterranean Sea, as well as in the Bahamas. The predictions for the Bahamas (Figure 8a   When compared to the observations, the predictions are generally significantly lower, with the MVP 3 scenarios reducing the vertical gap between model predictions and observations (Figure 8a-f). In the Bahamas, for each ice sheet model, the MVP 3 profile shifts the RSL curve upwards and closer to the observed transgression at MIS 5.3 ( Figure  8a). This implies that a slightly higher eustatic value (hence smaller ice sheet volume) might be required to explain the observed submersion. The observed +1.5 m sea level in Mallorca at 81 ka, on the other hand, stands out as a clear outlier and cannot be reconciled with the current ice-sheet and Earth rheological models (Figure 8b). Predictions for Croatia, Plemmirio, and Argentarola (Figure 8c-e) are in agreement with the observed vertical limits. However, the models cannot exceed the −20 m limit at 81 ka in Argentarola, thus implying that a further melting of ice sheets is required.
The predicted MIS 5.1 regional RSL variability, according to ANICE-SELEN and MVP 3, is characterized by higher-than-eustatic values towards the East-North East and around the continental margins (Figures 11 and 12). On the other hand, the predicted values are lower than the eustatic value at the center of the basins towards the South East (where the ocean loading term dominates). Interestingly, Plemmirio and Mallorca are very close to the eustatic signal ( Figure 12).
Predictions based on ANICE-SELEN and MVP 1-3 are in satisfactory agreement with the MIS 7.1-7.2 observations at Argentarola (Figure 8f).

Conclusions
The use of submerged speleothems has significantly contributed to reconstructions of the short-term (Holocene) and long-term (Middle and mostly Late Pleistocene) sea level changes in the Mediterranean. Even though there are still some open questions, the different groups that have been working on the speleothems from Mallorca, Croatia, and Italy for many years have carried out extensive and valuable work, which has significantly contributed to improving our knowledge on the sea level changes that occurred during the Mid and Late Pleistocene.
The present review reports the published elevations and radiometric ages ( 14 C and U-Th) of the submerged speleothems in the Mediterranean. Elevation data have been corrected for local tectonics (uplift or subsidence) and compared with results from the DWBAH flowstone on Grand Bahama Island (Figures 13-15).    [22], and stalagmite ASI (−18.5) from Argentarola cave [2]; black line: global sea level curve reconstructed by [88]; dark blue line: sea level drown with observed data.
Below, we summarize the main results obtained, following the key time periods investigated:


Holocene: Most of the data on the Holocene sea level curve in the Mediterranean result from the investigation of continental portions of submerged speleothems, marine overgrowth (serpulids), and Lithophaga boring mussels in samples collected in Italy and Croatia;   [81] suggests that the sea level was lower than the elevation of the submerged sample;  MIS 7.1: Three speleothems from the Argentarola cave ( Figure 13) clearly show a serpulid layer between MIS 6 and MIS 7.2 [2,22]. The sea level was therefore higher than 18 m b.s.l.. Based on the results from the DWBAH flowstone, lacking a hiatus at 197 kyr [11], the sea level was lower than 12 m b.s.l., in agreement with the global sea level [88].  [88]. However, the DWBAH flowstone shows a hiatus during MIS 7.3, which seems to be older than the marine layers from the Argentarola stalagmites. The duration of the sea level highstand was between 217.2 and 201.6 kyr [22];  MIS 7.5: The results from the stalagmites of the Argentarola cave and DWBAH flowstone indicate that the sea level was between 18 and 12 m b.s.l., again in agreement with the global sea level curve [88]. The duration of the sea level highstand was between 248.9 and 231.0 kyr [22].
Some GIA models have been published for the Mediterranean basin [5,36,38]. The latter two models have been tested in the field, with a comparison of coastal fossil deposits dated to MIS 5.5 and fossil tidal notches of the same age. The maximum reported differences for the entire Mediterranean basin for MIS 5.5 are less than 2.5 m in elevation. This value is consistent with the predicted GIA-driven RSL variability within the Mediterranean basin during MIS 5.5.
Here, we show, for the first time, that the regional RSL variability increases during interstadials MIS 5.1-5.3 and glacials (MIS 2), most likely as a consequence of the ice-loading-induced spatial RSL gradients. However, because the GIA signal also depends on Earth rheology, this further complicates the search for eustatic constraints in the Mediterranean Basin.
Overall, we argue that the observations call for a revised eustatic value (i.e., ice sheet volume) during MIS 5.1-5.3 and MIS 7.1-3. In particular, we speculate that a reduction of the ice sheet volume, and thus an increase of the eustatic value, is needed for MIS 5. 1-5.3. This is consistent with the recent findings of [99], which are based on a novel ice sheet modeling technique.
In particular, for MIS 5.1, based on the elevations below sea level observed in the ASI stalagmite (Figures 13 and 15, Table 3) between the well-dated continental layers of MIS 6 and MIS 3 at −18 m [2,21] and the elevation of the Bahamian flowstone [5] at −13.5 m, which is considered mid-field like the Mediterranean (Figures 7 and 8), and therefore comparable with the Mediterranean speleothems described in Figures 13, 14 and 15, it is possible to hypothesize a eustatic altitude of MIS 5.1 of approximately −16 m (average between −18.5 and −13.5 m). According to Figure 12, the Mallorca site is located around this eustatic altitude, while the Croatian site (samples K-14 and K-18, Table 4) would be raised for isostatic reasons by about 7 m, in agreement with [99] and the novel ice sheet modeling technique described.
However, we also argue that the high value observed in Mallorca (1.5 m above m.s.l.) cannot be reconciled with GIA and would require a drastic reduction of the ice sheet volume, which is not required at other sites.    Table A1. Ages of the flowstone DWBAH [11] corrected for the detrital 230 Th component. Age re-calculated represents the age re-calculated using the most recent 230 Th, 234 U, and 238 U decay constants of [72,73]. Age cor. represents the age corrected for the detrital 230 Th fraction using an initial 230 Th/ 232 Th activity ratio of 0.8 ± 0.8, as suggested by [75] for the speleothems from the Bahamas.