Submerged speleothems and sea level reconstructions: a global overview and new results from the Mediterranean Sea

: This study presents a global overview of the submerged speleothems used to reconstruct paleo sea level and reports new results from two stalactites collected in the Mediterranean Sea. Coastal cave deposits significantly contributed to the understanding of the global and regional sea level variations during the Middle and Late Quaternary. The studied speleothems cover the last 1.4 Myr and are focused mainly on Marine Isotope Stages (MIS) 1, 2, 3, 5.1, 5.3, 5.5, 7.1, 7.2, 7.3 and 7.5. Results indicate that submerged speleothems represent extraordinary archives providing detailed information on former sea level changes. The two stalactites collected in the central Mediterranean Sea, at Favignana and Ustica islands (Sicily, Italy) are both characterized by continental, phreatic or marine layers. The U-Th and 14 C ages of the new speleothems provide results of great interest for relative sea level changes over the last 1000 years. GIA-driven RSL variability the interstadials and glacials. arguet that the geoidal deformations,

The study of submerged speleothems in coastal caves significantly contributed to constrain past sea level variations for the last 1.4 Myr, especially for the climaticallyimportant Marine Isotope Stages (MIS) 1, 2, 3, 5.1, 5.3, 5.5, 6.5, 7.1, 7.2, 7.3 and 7.5. Submerged speleothems can be dated using the U-series disequilibrium methods that provide precise and accurate age determination [e.g. [1][2][3][4][5]. 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 seawater invaded the cave and reached the speleothem.
In addition to speleothems, other biological-physical markers have been used to reconstruct past sea level variations, including shallow-water corals, lagoonal and marine cores, tidal notches and archaeological remains that are well connected with sea level fluctuations. However, most of these indicators suffer of large uncertainties in the age determination due to sample alteration and/or in absolute elevations (e.g. paleowater depth for corals). Submerged speleothems have been studied so far: 1) continental (vadose) carbonates (in few cases containing hiatuses), 2) speleothems with overgrowth(s) of marine organisms (e.g. the marine polychaete worm Serpula massiliensis, bryozoans and sponges), and 3) speleothems with phreatic overgrowth(s). The general features of these carbonate cave deposits are described in [6]. In the Mediterranean Sea, more than 300 submerged speleothems were sampled in 32 caves and more than half were analysed for paleo sea level reconstructions [6] (Figure 1). Other speleothem-based sea level records come from few sites in the Caribbean region (Yucatán Peninsula and Cuba), Bahamas, Bermuda and Japan (Minami Daitō island) (Figure 2 and 3), which show very different geological and stratigraphic structures. Here we provide a comprehensive review of published sea level reconstructions based on submerged speleothems and present new data from two stalactites collected in the Ustica and Favignana Islands in the central Mediterranean Sea. Published data for the Mediterranean Sea have been recently summarized in [6]. Among the most important findings of speleothems retrieved from flooded coastal caves worldwide and often showing hiatuses or marine overgrowths, a special attention deserves the submerged cave deposits from Bahamas [3,4,7], Bermuda [8,9], Yucatan Peninsula [10,11], Cuba [12], and Japan [13]. [14,15] provided an exhaustive review on the use of sea-level indicators preserved in coastal caves, including speleothems and cave sedimentary deposits.   Table 1.

Western North Atlantic-Caribbean region and Northern Philippine Sea
Sea level reconstructions based on vadose speleothems and marine overgrowths from the Yucatán Peninsula have been obtained by [10,11], through the analysis of samples collected at elevations between +1.5 and -15.1 m. The Yucatán Peninsula is located in the southeastern Mexico and separates the Caribbean Sea from the Gulf of Mexico. It is almost entirely composed of carbonate and soluble rocks, being mostly limestone although dolomite and evaporites are also present at various depths [16]. The whole Yucatán Peninsula is an unconfined flat lying karst landscape. Sinkholes, locally known as cenotes, are widespread in the northern lowlands [17]. The samples studied by [10,11] formed within Plio-Pleistocene reef and back-reef limestones in a tectonically stable area with little deformation [18].
Phreatic overgrowths and speleothems from Santa Catalina Cave on the northern coast of Cuba were investigated by [12,19]. This cave is carved into Late Pliocene-Early Pleistocene reef limestones [20] and hosts a variety of normal and exotic speleothems, including the mushroom-shaped speleothems, that precisely mark the water level.
Several speleothems and flowstones were studied by [3,4] in the Bahamas islands for past sea level reconstruction. The Bahamas islands include aeolian sands and limestone built on the basement rock of the Florida-Bahamas Platform and are considered a quasi stable area with a low subsidence rate (5 m/100 kyr) [3,21].
Speleothems and marine aragonite overgrowths were analysed by [8,9] in Bermuda. This volcanic seamount is located in the central-western Atlantic and has a carbonate cap with eolian limestone interbedded with paleosol layers. The island's volcanic basement rock is relatively shallow (75 m bsl) and includes 700 m of tholeiitic lavas aged 33 million years. The oceanic crust around the island is about 120 million years old. Bermuda is considered a tectonically stable island but a "near field" with respect to the forbulge ice [22].
Phreatic overgrowths on speleothems (POS) were studied by [13] in the Minami Daito island in Japan (Northern Philippine Sea) to reconstruct Holocene sea level variations. The island is formed as an atoll in the mid-Eocene (ca. 42 Ma ago) on a subsiding volcanic edifice [23] and experienced an average uplift rate of ca. 5 cm/kyr for the last 125 kyr. Based on fossil tidal notches and dated corals the island is considerd stable [24]. Table 1 reports the basic information of the cave deposits from the Western North Atlantic-Caribbean region and Minami Daito island the discussed in the present review. For further information regarding the cave morphology, sampling method, date of sampling and the error associated to the radiocarbon and U-Th ages, the reader should refer to the original publications [3][4][5][10][11][12][13]19,[25][26][27][28][29][30][31][32][33][34][35][36][37][38].

Mediterranean region
The coastline in the Mediterranean Sea (MS) extends for about 46,000 km, with approximately half of it being rocky and characterized by plunging cliffs, sloping coasts, screes and shore platforms [39]. The MS experiences warm and dry summers and mild and rainy winters owing to its geographic position at the boundary between the subtropical and mid-latitudes zones [40]. It is a highly evaporative basin that is connected with the Atlantic Ocean through the Strait of Gibraltar (sill depth ~ 300 m) and with the Black Sea through the Strait of Dardanelles (sill depth ~100 m) and the Bosphorus Strait (sill depth ~65 m). The coastal geometry and bathymetry of the MS, together with the weather conditions, strongly control the tidal amplitude that varies from place to place along the Mediterranean coasts. The average tidal amplitude is about 40 cm, with the exception of the large tides in the Gulf of Gabes (Tunisia) and in some areas of the North Adriatic Sea where they may reach amplitudes up to 1.80 m. However, in general, the Mediterranean tides have lower amplitudes compared to the other oceanic basins Figure 1.

Ustica Island
The Ustica Island (central MS, Figure 1) represents the emerged part of a submarine volcano that rises more than 2200 m from the bottom of the southern Tyrrhenian Sea and lies on a deformed continental crust (15-20 km; [41]) along the Africa-Europe plate boundary [42]. The island is located ca. 60 km north of Palermo (northwestern Sicily coast) and covers an area of 8,6 km 2 , reaching a maximum elevation of 248 m a.s.l. The structural pattern of Ustica Island and its offshore is characterized by ENE-WSW and N-S-directed extensional and compressional fault systems [42][43][44]. Ustica island is made up of alkalibasaltic to trachytic composition lava flow and pyroclastic products [45][46][47][48][49] locally covered by marine and continental sedimentary deposits [50]. The present morphology of the Ustica Island and its offshore is the result of the prolonged interplay between volcanic, tectonic and eustatic processes [42,51] active since the Pleistocene and which contributed to the birth of the island occurred between 750 and 130 ky as dimostrated by stratigraphic and geochronological studies corried out on the exposed rocks [49,52].
During the Ustica Island volcanic history, there have been several overlapping cycles of marine ingression and regression, resulting from Middle-Upper Pleistocene glacioeustatic movements. The sea-level stands related to these oscillations generated typical sedimentary terraces and submerged plains [42,51] which were later uplifted by tectonics at variable height from their original position. Five orders of level surfaces have been identified at heights ranging from about 100 m asl for the oldest and highest terrace, to 5 m asl for the most recent and lowest one [49,53]. The occurrence of marine terraces of variable age at decreasing elevation from the oldest to the youngest, helps in defining the amount of uplift recorded by the island, which has to have been no more than 120 m, at least since the formation of the oldest level surface, that is 350 ka [49,53].
Regarding the geological long term vertical movements, the maximum elevation of MIS 5.5 is 33 m [52] with an uplift rate of 0.21 mm/yr. The recent discovery of fossil barnacle (Chtamalus) bands radiocarbon dated to 1823 ± 104 yr AD and covered by a flowstone dated to 1894 ± 24 yr AD in the Grotta Segreta [51], reveals a cosismic event. An uplift of at least 40 cm is suggested that might be related to the strong earthquake in 1906, which obliged all the people to leave the island. This example advises against the use of instrumental data alone for the long-term projections of future sea level rise by 2100 (e.g. GPS data do not reveal vertical movements at Ustica).

Favignana Island
Favignana (Egadi archipelago, Figure 1) is the closest island to Sicily. The archipelago are formed by carbonate and terrigenous Meso-Cenozoic deposits pertaining to different palaeogeographical domains, covered by unconformable Plio-Quaternary pelagites and clastics. MIS 5.5 calcarenite and biorudite outcrops [54,55] are widely distributed in the eastern sector of Favignana, from 5 up to 10 m [56]. Favignana and the Egadi Figure 4 The stalactite sampled at Ustica. a) longitudinal section with sampling area of U\Th ages, refer to Table 2. b) horizontal section with radiocarbon ages, samples US1, US2. US3, refer to Table 3. c) the stalagmite sampled.
Islands may be considered as a stable or slightly uplifted area at least for the last 125 kyr [57][58][59]. Many outcrops of polygenic conglomerates contain fossil shells of Senegalese

Speleothem sampling in Favignana and Ustica
A partly submerged stalactite was collected in 2016 in the Stalattiti cave (Figure 1, 4) on the Ustica island (about 1.4 km from the harbor Figure 5) using a hammer. At the time of the sampling, the deepest portion of the speleothem was 40 cm below the sea level. Elevation measurements were corrected for tide and pressure in order to relate them to the mean sea level (MSL), by using data from the tide gauges of Palermo (www.wxtide32.com) and from the meteorological website www.wunderground.com. The correction was 7 cm for a 1019 hPa pressure. The uncertainty on the elevation measurements of the partly submerged speleothem relative to the sea level, which is on the order of few centimeters, can be considered negligible, even though large uncertainties may arise from the assumed sea level position of markers (Table 1).
A partly submerged stalactite was retrieved in 2014 in the Stalattiti cave (Figure 1, 6) on the NE coast of Favignana island using a hammer. At the time of the sampling, the deepest portion of the speleothem was 30 cm below the sea level. Elevation measurements were corrected for tide and pressure in order to relate them to the mean sea level (MSL), by using data from the tide gauges of Marsala (www.wxtide32.com) and from the meteorological website www.wunderground.com. The correction was 2 cm for a pressure of 1015 hPa.
The uncertainty on the elevation measurements of the partly submerged speleothem relative to the sea level, which is on the order of few centimeters, can be considered negligible, even though large uncertainties may arise from the assumed sea level position of markers (Table 1).

U-Th and 14 C dating
Three fragments of spelean calcite were collected within the vadose portion of the stalactite from Ustica (Figures 4,5) and analysed through the U-series disequilibrium method to determine their absolute age. The fragments were carefully cleaned using a small diamond blade to remove any visible contamination and dissolved with diluted HCl before being mixed with known concentrations of 229 Th, 233 U and 236 U, calibrated against a Harwell Uraninite solution (HU-1) assumed to be at secular equilibrium. The uranium and thorium fractions were separated using UTEVA resin (Eichrom Technologies, USA) and analysed by standard-sample bracketing using a Thermo Scientific Neptune Plus multi-collector inductively coupled plasma-mass spectrometer at the Laboratoire des Sciences du Climat et de l'Environnement (LSCE) (Gif-sur-Yvette) following the protocol developed at LSCE [60]. The 230 Th/U ages were calculated from measured atomic ratios through iterative age estimation [61], using the 230 Th, 234 U and 238 U decay constants of [62] and [63]. A correction was applied for the non-radiogenic (detrital) 230 Th fraction using an initial 230 Th/ 232 Th activity ratio of 1.5 ± 0.75. However, due to the low 232 Th concentration (< 0.05 ng/g) and high 230 Th/ 232 Th ratio (> 700), the detrital correction has a negligible effect (   Radiocarbon dating for the Ustica and Favignana samples was carried out at CEDAD-Centre for Applied Physics, Dating and Diagnostics at the University of Salento, [67]. The selected samples were carefully removed from the matrix by using a diamond blade. The extracted portion of the sample were then analyzed at the optical microscope in order to identify and mechanically remove adering particles or unwanted farction of the matrix. The sample were then attacked with H2O2 in order to remove the uppermost layer and then hydrolised under vacuum to CO2 by usign HCl. Extracted CO2 was then catalitically reduced to graphite which was then used to measure the radiocarbon age by AMS (Accelerator Mass Spectrometry). The radiocarbon dating beamline was used, based on a 3MV Tandetron accelerator (Mod. HVEE 4130HC) to measure the radiocarbon content. After correction for machine and processing background the conventional radiocarbon age of the samples was calculated. For continental samples conventional radiocarbon ages were calibrated by using the INTCAL20 calibration curve while marine samples were calibrated by using the MARINE20 [68] curve and a local ΔR=-88±50 yrs [69]

Numerical modelling of glacial-and hydro-isostatic adjustment and sea-level change
Land-based sea-level indicators such as speleothems, keep track of past sea-level changes that are referenced to the surface of the solid Earth. Because the latter can be affected by vertical motions, the recorded sea-level fluctuations in speleothems are defined as local Relative Sea Level changes (RSL).
In particular, ice-driven sea-level changes stem for the contribution of ocean mass variations (in response to continental ice mass variations) as well as solid Earth deformations and vertical geoid perturbations [70][71][72]. In fact, when water mass is transferred from ice sheets to oceans, and viceversa, solid Earth deforms in order to restore the isostatic equilibrium under a different surface loading setting. Solid Earth deformations and ice masses, then, behave as density anomalies, thus affecting the vertical position of the mean sea surface, which is an equipotential surface of gravity (geoid). This process, which is known as glacial-and hydro-isostatic adjustment (GIA), is responsible for regionally-varying RSL changes that are modulated in time the the viscous flow of the Mantle material.
Because of mass conservation, the ocean averaged RSL change corresponds to the global mean, or eustatic, sea-level change. Based on the geographical distance from the ice sheets, the local RSL change can be very closte to the eustatic (ice distal o far-field regions), or significantly different, if not opposed to the eustatic (ice proximal areas). The largest deviations from the eustatic signal are found where ice sheets grow and shrink. In particular, during glacial perionds and underneath the growing ice sheets, the crust subsides and forces the upper mantle material to flow outwards, thus causing the uplift of the ice-proximal peripheral area that forms the so-called forebulge. Between the forebulge and the ice margins, the local RSL rises due to the crustal subsidence and the geoid uplift that is caused by the gravitational pull that the ice masses exert on the ocean water. On the other hand, on top of the forebulge, the local RSL drop is much larger than the eustatic during glacial periods because of the combination of crustal uplift and geoid subsidence. When the ice sheets melt, the ice-covered areas undergo uplift, while the forebulge ares experience a strong subsidence that, combined to the eustatic signal, results in a much higher-than-average RSL rise. The vertical motions that affect the forebulge areas are not just relevant for the ice proximal location, but also for the farfield equatiorial and continental coastal areas. Here, in fact, despite the distance from the ice sheets, the RSL changes that follow major deglaciations may strongly deviate from the eustatic and show up with an opposite sign. In fact, because of mass conservation, ocean water flows towards the subsiding peripheral forebulges, and locally results in a RSL drop that follow an initial highstand. This process is known as ocean syphoning.
Interestingly, some of the Caribbean and Yucatan speleothems, are sitting in areas where the compex forebulge ocean syphoning dynamics operate.
Therefore, RSL indicators provide constraints on past ice-sheets chronologies and on the relevant Earth rheological parameters that regulate deformation of the planet in response to surface ice and water loading.
To compute the local ice-driven relative sea-level changes, the gravitationally selfconsistent sea level equation (SLE; 1-4) must be solved for a prescribed ice sheets chronology and Earth rheological model. 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 of the solid surface of the Earth, the difference of the two being the relative sea-level (RSL) change. We employ SELEN [70][71][72] to solve the SLE and retrieve the local GIA-modulated RSL curves at each site considered in this paper. The latest, up-to-date version of open source model SELEN includes the rotational feedbacks as well as the form treatment of the time-dependent ocean function (i.e. variable coastlines).
We assume a self-gravitating, rotating, spherically symmetric, radially stratified, deformable but not compressible Earth model. The latter is 1-dimensional and all the relevant rheological parameters depend on the Earth's radius only. 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 [73] 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 5). We combine this with a lithosphere thickness of 90 km. We use this as a reference model to be combined with three ice sheets chronologies.

VM2
LT (  We employ the following ice sheets models: Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 6 May 2021 doi:10.20944/preprints202105.0100.v1 • ICE-5G [73]: this global model describes the ice-sheet thickness variations over North America, Eurasia, Greenland and Antarctica during the last 123 ka. Geological and modern geodetical observations are used to constrain the ice thickness chronology between 26 ka and present-day [6]. Through an iterative procedure, the SLE is solved for an a priori ice sheets configuration and a prescribed fixed solid Earth model (mantle viscosity profile and Lithosphere thickness). The predicted RSL elevations and modern vertical rates are compared to, respectively, paleo RSL indicators and measured crustal and geoidal vertical velocities. The differences between predictions and observations are then used to to modify the ice sheets models until a satisfactory solution is found. The ice sheets volumetric evolution between 123 and 26 ka is tuned to the  18 O curve [74] and loosely constrained by surface geological evidences that define the ice margins.. Here we combine two ICE-5G chronologies in time in order to cover the time span under consideration and to evaluate the GIA-modulated RSL response to the melting of the MIS 6 ice sheets during the MIS 5e. The latter is characterized by a ~0.9 m eustatic sea-level highstand above present-day msl between 129.5 ka and 123 ka, mostly coming from a reduced (with respect to present-day) Greenland Ice Sheet.
• ICE-6G [75]: this model follows from an improvement of the ICE-5G. Here we apply the same repetition in series of ICE-5G. According to the ICE-6G, the maximum MIS 5.5 eustatic highstand between 129 ka and 123 ka is ~3.1 m above present-day msl as a consequence of reduced Greenland and West Antarctic Ice Sheets • ANICE-SELEN [76,77]: this global chronology model is the result of an inverse forward modelling procedure where the  18 O stack [74] is decoupled into global ice sheets 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 feedbacks. ANICE-SELEN is not yet constrained by geological or instrumental data. Here we use the original ANICE-SELEN chronology, which is characterized by ~2.5 m eustatic highstand during MIS 5 e (125-117 ka), in response to a reduced Greenland Ice Sheet volume.
The temporal volume variations of three ice sheets models between are therefore dependent, to a certain degree, on the benthic oxygen curve [74].

Results
We divived the paragraph in two parts, the first one regarding published data in the Western North Atlantic-Caribbean region and Japan (par. 3.1), while the second one regarding new data in the MS (par. 3.2).

Bermuda
The Western North Atlantic-Caribbean region, and Bermuda in particular, is strongly affected by the vertical motions of the periglacial forebulge that forms as a result of upper Mantle flow in response to the expansion of the Laurentide ice-sheet during glacial periods. U-Th age uncertainties range between ±7 ka and ± 1 ka. These data support previous controversial indications that Bermudian RSL was significantly higher than RSL at other locations during MIS5.5 and MIS 5.3 Figure 2, Table 1.    Figure 6).

Santa Catalina Cave (Cuba)
The Santa Caterina cave is carved in a marine terrace, which occurs along the Cuban coast. [19] suggested a general slowly uplifting trend during the Pleistocene, even if U-Th datings of these terraces are difficult to achieve considering local conditions. The  [11] found and dated 4 marine levels between 8.9, 8.6, 8.4 6 ka cal BP collected on a spectacular stalactite sampled at -6 meters in a Cenote in the Yucatan Peninsula. He interpreted dated levels as sea level fluctuations: "therefore, it seems highly likely that four relative sea-level 'Oscillation Events are responsible for submergence and re-emergence.

Yucatan
These oscillations are presumably due to GIA movements since all the area is considered tectonically stable". Authors also discussed the results obtained by [78] on the U-Th ages from the Acropora palmata fossil corals and suggested that the error related to the paleoposition of the sea level can reach 10 m, while a submerged speleothem with marine layers show a centimetric error bars, also considering the local microtidal amplitude of 0.3 m (Table1, Figure 2).
The results summarized in Table 1 Table 4).
The speleothem shows an emerged part made of brown coloured subaerial carbonatic facies, and a 6 cm white carbonate facies suggesting a phreatic genesis. Between the aforementioned facies we also noted and analized (Figure 4) a marine overgrowth made up of serpulids. Radiocarbon ages of subaerial carbonatic overgrowth provided an age of 23.9 kyrs cal BP, while the marine overgrowth of 1.3 kyrs cal BP and the last presently submerged phreatic layer provided an age of 0.7 kyrs cal BP.

Ustica
The submerged cave (location and sketches in Figure 5) has a small entrance at -13 m bsl. A 10 m x 10 m room with open-air at the sea level was found inward. The cave roof was estimated to be about 2 m in height asl. Hundreds of stalactites hang from the vault of the cave, all formed only of continental carbonate. There are alsosome plant roots. In the submerged portion, both on the cave walls and the submerged stalagmites and stalac tites, Polychaeta overgrowth can be observed. In some cases, they appear to be very large, up to 7-8 mm in diameter) Figure 4. The stalactite sampled around the present-day sea level is sprinkled with living Polychaeta. Once dissected, it showed (see Figure 4) a dark brown central part due to continental/phreatic deposition occurred in subaerial environment. On the contrary, microscope analysis show that the external part of the speleothem, made oflighter tones, presents a mix of both marine, with polychaetes and continental origin ( Figure 5). The age of this stalactite resulted to be about 1000 years (Table 2). Radiocarbon datings carried out on the outermost part of the speleothem (Table  3), show an older age while maintaining a stratigraphically correct sequence.

GIA numerical predictions
For the three ice sheets models, the predicted RSL curves computed at each location and for each of the four mantle profiles, tend to be distributed above the eustatic curves, with Japan being the only exception as it is almost continuously below the eustatic Figure 9 Cumulative predicted RSL curves for ICE-5G (a,d), ICE-6G (b,e) and ANICE-SELEN (c,f) ice sheets models in combination with different mantle viscosity profiles (see legend and Table 5 and 6 for details). The black curves represent the eustatic sea level variation for each ice sheets model. The colored curves represent the RSL change for each site under consideration. throughout the last 250 ka (Figure 9 a-c). The Mediterranean sites (Ustica and Favignana) are very close to the eustatic component and do not show significant variability as a function of mantle viscosity profiles. The same holds for Japan, where local RSL variability is larger for ICE-5g and ICE-6g and negligible for ANICE-SELEN. In the Caribbean and Yucatan, instead, the local RSL variability is of the order of 20-30 m during glacials and interstadials, and decreases during the MIS 5.5. This implies that Bermuda and Bahamas, being closer to the north American ice sheets margins, are more sensitive to GIA-modulated vertical motions, which can be emphasized by specific viscosity profiles.
RSL predictions vs data (MIS 5.1 -5.3) Overall, the predicted RSL peaks are generally lower than the observed terrestrial limiting elevations (deposition) as well as the hiatuses and the marine transgressions inferred from the speleothems (Figure 10). The largest RSL variability is predicted at Bermuda (Figure 9, 10c) because of the relative proximity to the ice sheets margins and peripheral forebulge. Here, though, the maximum predicted highstands at MIS 5.1 and MIS 5.3 are, respectively, ~20 and ~15 m lower than the speleothem elevation. Furthermore, the predicted highstands at MIS 5.1 and 5.3 are, respectively, ~5 and ~13 m higher than those at Yucatan (Figure 10 b). These differences are in line with the results of [9], but significantly smaller. We argue that this might be related to the very different ice sheets models (ICE-1 and ICE-2) used by [9]. Also, our predictions at Bermuda are only slightly higher than at Bahamas (Figure 10d), implying that both sites share the same GIA response for the three ice sheets models and four mantle viscosity profiles.

Discussions
Considering data reviewed in the paper and in [6], the use of submerged speleothems provided remarkable results for the global sea level change history during the Holocene and for the long-term sea level change reconstruction, such as during the Late to Middle Pleistocene, in particular for MIS 1, 3, 5.1, 5.3, 7.1, 7.2, 7.3 7.5, MIS 9 highstands (Tables 1, 2 ,3, 4). In particular, figures 11 and 12 depict all the results for the Caribbean and Japan submerged speleothems and some of the main results of the aforementioned review for the Mediterranean area [6] have been compared with the global curve of [79]. In the table 1 presented by [6], a speleothem sampled in Italy (Custonaci) show marine and continental levels and seems to be the older: aged at MIS 25-37 (1.4 MA).

Favignana stalactite
This speleothem provides an interesting result concerning sea level change for the last 1300 years, due to the fact that sea level data are very scarce in tectonically stable areas of the Mediterranean sea [80,81]. The submerged portion of the stalactite is truncated and bored by Lithophaga lithophaga, The cave ( Figure 8) could be recently greatly enlarged by marine processes, such as storm waves, also considering the fact that the almost all the other stalactite near the one we have studied are truncated. Therefore we hypothesize the evolution following described: 1) During LGM the sea level was far from here, at about -129 m bsl [82] and the speleothem presumably formed a column ( Figure 11) aged 23.9 ka BP (Table 3).
2) sea level rose up until 1395 years cal BP. The sea level stillstand at -0.32 m shows a stillstand which promote the growth of colonies of serpulids and sponges that formed on the broken surface of the speleothem, as reported at the Argentarola cave by [38,83].
3) the following millimetric layer of phreatic carbonate (Figures 11) shows that the aforementioned sea level stillstand was up to -0.35 m, considering the deepr phreatic level,. Radiocarbon datings of the most recent level is 705 kyrs cal BP (Table 3) . 4) the phreatic level was then interrupted by marine ingression, so now the submerged part of the stalactite is bored by Lithophaga sp., while the external carbonatic surface of the submerged part of the speleothem is colonized by small Serpulids, tetracorals, limpets and sponges. The levels of phreatic carbonates have been used in the Mediterranean as sea level markers by [28,30]. In our case, the entire phreatic formation has not been preserved, but only some non-eroded fragments were observed. For this reason we can hypothesize (as happens for the cave from Palma de Mallorca) a thickness of the phreatic deposit (in the shape of a ball) on average as high as the average tide of Favignana (annual mean: 32 cm, OSU tidal prediction model, http://volkov.oce.orst.edu/tides/otps.html), it could also be assumed a greater thickness in depth of the phreatic.
Furthermore, the age correction due to dead carbon is not difficult in this case. It showed that the value of the term DCP is known. In the case of the Favignana stalactite, the age of 806 BP (calibrated 705 cal BP), that is, of the sample S3 would become 478 BP (calibrated 518 cal BP) with a dead carbon percentage of 4%. Even a percentage of DCP of 10% would be enough to make the sample modern. In other words, if the DCP percentage of "dead" carbon in the cave were 10%, the sample would correspond to the present. This is obviously not reasonable, also because the stalactite is submerged. This suggests that the effect of carbon dead is very low (or negligible) in the studied cave. After all, the dating of 705 BP, without considering dead carbon effects, fits perfectly with the curve and this cannot be considered a coincidence. We took into account this question, but the analysis of datings leads us to think that this effect is completely negligible.
In Figure 12, the stalactite age with different ice models for Favignana (see caption figure 10 and Table 3) Vermetid reef from S.Vito lo Capo (see Table 4). The sea level change (blue) curve seem to intercept the Warm Medioeval Period and the Little Ice Age. and other geomorphological and geoarchaeological markers of the same time range are compared. The Serpulid overgrowth dated back to 1395 cal BP yields quite above almost any curve. The phreatic speleothem dated back to 705 years is instead in agreement with the average of the curves. Then, we added, as shown in Figure 16, some archaeological and biological markers included in the 2 ka time range. For example a fishtank in Salerno at -1.4 m bsl aged back to 1.9 kyrs from [84]. At about -0.40 m we also Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 6 May 2021 doi:10.20944/preprints202105.0100.v1 found a less precise marker (a walking surface of the basilica of San Nicola [81] at -0.32 m in Middle Ages. We also added, some radiocarbon ages on Vermetid reefs [65] at -0.3 m dated back to 169, 163, 189 years cal BP and another one at -0.4 cm dated back to 279 years cal BP. Even if we advance these hypotheses with few data, it would seem to notice a peak of the sl curve corresponding with the Medieval Warm Period between 900 and 1200 years AD and a cooling corresponding with the Little Ice Age between 1550 and 1900 years AD). This trend is in agreemen with what was also found in Israel by [80].

Ustica stalactite
The stalactite sampled at the present-day sea level is about 0.8 cm long, with a diameter ranging from 4 to 8 cm (Figure 4 c). As shown in Tables 2 and 3, this speleothem also covers a very narrow time range of about 900 to 1000 years (Figure 4 a). There is a high concretion rate which could surprise, given that we are in a volcanic area. But the Ustica lava contains many fragments of carbonate platform erupted together with the lava, and often these carbonate ouotcrops werefound within lava outcrops up to few dozen meters. In some cases oalso with a modern tidal notch [51]. The most recent portion of the speleothem is formed by large living serpulid tubes. In the portion, between the continental part and the living polychaete worms, a continental sediment showing a different color was observed (Figure 4 b) together with porosity and with polychaete tubes, sponges incorporated in the sediment, is clearly visible.
The mixed deposition of continental carbonate with serpulid tubes (Figure 5b) costituites a clear signal of positive and negative vertical coseismic movements. The radiocarbon age of this outer portion, shows older values (albeit in perfect stratigraphic order) with respect to the central part (Figure 5a) of the satalactite aged with U-Th. The last part is younger. Our explanation is that the dead carbon of the very ancient carbonate inclusions in the Ustica lava have altered the date but not the stratigraphic order. The average tide amplitude in Ustica of 36 cm (https: / /www.tide-forecast.com/) and the length of the sampled speleothem of about 85 cm (Figure 4 c), add a proof that over the last 1 ka the island of Ustica was affected by positive and negative cosismic events. This consideration is confirmed by the very recent coseismic movement described for the near, Grotta Segreta, with a coseismic uplift of 0.4 m due to the 1906 earthquake [51]. The studied stalactite indicates both positive and negative vertical movements, with a positive average component (uplift) that can be deduced from the altitude of the MIS 5.5 deposits containing Persistrombus latus at 35 m altitude [52] and by the presence of evident marine terraces raised on the lava, with an uplift rate ranging from 0.13 to 0.24 mm/yr [53] (see also Figure 5).

Caribbean and Japan speleothems
U-Th datings on the speleothems sampled at Bermuda, Bahamas and Yukatan caves ( Figure 2, Table 1) provided detailed information on the timing, duration, and altitude of sea level highstands during late and middle Pleistocene. [11] analyzed a peculiar speleothem that preserves 4 marine hiati between 10.9 e 58 kyrs cal BP. The region under consideration is considered to be tectonically stable and it is very unikely that GIA could have been responsible for these anomalous crustal and/or sea surface fluctuations.

MIS 2 -Holocene
[13] studies a cave on the island of Minami Daito in Japan, which was made of a beautiful phreatic overgrowth. The authors investigated the relationship between the highest POS level, which represents the maximum elevation reached by the palaeogroundwater level and the observed modern highest peak. Regardless of a relevant local mean tide fluctuation of 1.5 m, the authors found that between 5.4 and 4.6 cal BP, the local sl reached a maxmimum highstand of ~ 35 cm during the mid-Holocene. This can be explained by the GIA-drivem ocean syphoning.  Figures 13 and 14 show that the GB-89-25-5c speleothem that was sampled at an elevation of -18.1 m bsl in Bahamas [4] preserves a hiatus between 63.7 ± 1.7 kyrs and 58.5 ± 0.7 kyrs (MIS 3.5, Table 1 Figure 6). The lack of RSL data related to MIS 3.3 (and MIS 3 in general) make this a very unique and important indicator from a nearly tectonically stable area. Figure 13 shows that the global sea-level curve from [79] is characterized by a much lower highstand of -48 m bsl.

MIS 3
Our GIA predictions for the three ice sheets models and mantle viscosity profiles are characterized by lower RSL peaks at MIS 3 (Figure 9 a-c). Because the three ice sheets models are all somehow tuned to the δ 18 O stack curve, we argue that their volumes might be biased towards larger values, while the speleothems here described might require a volume reduction throughout the MIS 3. The latter would, in fact, result in a much higher eustatic level. Interestingly, this is supported by the laters ice sheets model reconstruction from [85] that, for the MIS 3, provides a eustatic peak of ~28 m bsl. The authors used ice proximal glaciological and geological evidences to contrain the maximum extent of ice sheets through time and reconstruct the ice thickness by modelling the 2D ice profiles based on viscoplastic ice rheological flow law and 2D bottom shear stress. The model, therefore, is fully independent from far-field RSl observations and deep sea sea-level proxies.

MIS 5.1
The Yucatan speleothems do not record any hiatuses during the MIS 5.1 highstand, which therefore did not exceed the elevation of -14.6 m bsl [10]. Similarly, the DWBAH flowstone from Bahamas (corrected for vertical tectonics), does not show any evidence of local sea level higher than -13.5 m bsl ( Figure 13). Hiatuses are also lacking in the GB-89-25-5c speleothem (Bahamas; [4]), which therefore limits the maximum highstand to -18.1 m bsl. Yucatan and Bahamas are in an intermediate-to-far-field position with respect to the ice sheet, where ice-loading related deformations and ocean syphoning might prevail ( Figure  15). Our GIA-modulated RSL curves for Yucatan and Bahamas (see Figure 10 b,d) are all lower than the measured elevations (also corrected for tectonics). [9] explain the elevation of the MIS 5.1-5.3 hiatuses at Bermuda, sampled at +1.5 m (Table 1), as a consequence of the local GIA. The island, in fact, is considered tectonically stable but very  Our highest MIS 5.1 sea-level prediction in Bermuda reach -18 m, while lower elevations of -25 and -20 m are predicted, respectively, at Yucatan and Bahamas (Figure 10 b,d). This confirms, indeed, that Bermuda is more affected by the GIA-related forebulge processes than Yucatan and Bahamas ( Figure 15) The elevations inferred from the Caribbean speleothems of Bahamas and Yucatan (-13.5 to -18.1 m) are somehow comparable that observed in the Mediterranean site of the Argentarola cave (-18.8 m) [6,86]. We argue that an elevation of -19 m could be representative of the intermediate-field areas such as the Caribbean and the Mediterranean Sea, in agreement with Linberg (1999) and [85]. Therefore, by removing the expected GIA contributions, we also hypothesize an eustatic elevation for the MIS 5.1 of approximately −17 m, obtained as average value between −18.8 and −13.5 m.

MIS 5.3
The speleothem from Bermuda which lyes at+1.5 m asl preserves several hiatuses and is likely affected by the forebulge processes (Figure 9 a-c, 10c and 16).
The DBWAH flowstone from Bahamas at -13.0 m bsl shows a MIS 5.3 hiatus that, again, cannot be explained with our GIA results.
A hiatus is recorded in Yucatan, where the sea level might have been higher than -11m after 109 ka, but lower than -5 m.
As it occurs with the MIS 5.1 elevations, our GIA predictions for the MIS 5.3 highstand cannot reach the observed elevations (Figure 9 a-c, 10 ab-d). Again, we argue that a reduction of ice sheets volumes is required to shift the eustatic value upwards. We suggest a eustatic highstand between -12 and -6 m.

MIS 5.5
All the submerged speleothems from the Caribbean seashow hiati between -6 m bsl at 117 kyrs (Yucatan), -18.1 m bsl at 148 kyrs (Bahamas), at -14 at 134 kyrs of the flowstone DWBAH. The only difference is the timing of the transgression. In Bermuda the hiatus of MIS 5.5 is between 114 kyrs and 137 kyrs.

Middle Pleistocene
The only age oldest than MIS 5.5 was provided by the DWBAH flowstone [3,5] which contains 5 marine hiati (MIS 9, 7.5, 7.3, 5.5, 5.3) We can provide some considerations on differences in tide amplitude between the Carribean sea and MS. Table 1 shows the maximum tide amplitudes in Bermuda and Bahamas of about 1.3 m and lower than 0.5 m in the other sites in the MS. Some corals used by some key studies, such as [87,88], used Acropora palmata with a depth range from -1 m bsl to -20 m bsl (as measured by F. Antonioli in Barbados). In our opinion the precision in elevation of the submerged speleothems is higher than corals. This is due to the fact that speleothems are precise indicators with respect to the mean sea level for tides lower than 0.5 m, (Table1).

Conclusions
Submerged speleothems can grow for significantly long temporal spans and, therefore, can record multiple intervals during which deposition stops and/or erosion takes over in response to sea-level highstands. The alternance of continental depostional layers and hiatuses or marine layers records, therefore, the long-term climat-related global ice-sheets fluctuations. In this review, we show that the use of submerged speleothems significantly contributed to constrain the short-term (Holocene) and long-term (Middle and mostly Late Pleistocene) sea level changes in the Caribbean, Japan and Mediterranean Sea.
The longest time span is recorded in the DWBHA flowstone [3] at Bahamas and covers the last 326 ka with 5 hiatuses recorded. The Bahamas speleothem GB-89-25-5 covers the time span between 80 and 24 ka and contains one hiatus. At Bermuda, the QB speleothem covers the time span between 70 and 146 ka and records 4 hiati. Several speleothems from Yucatan [10] cover a time span ranging from 61 to 115 ka, with one evident hiatus. In the Mediterranean Sea, the ASI speleothem from Argentarola [86] covers the last 210 ka and records 4 hiatuses. The K-18 stalactite from Croatia [6,26] spans from 55 to 131 ka ka and contains 1 hiatus.
Overall, the discussed speleothems reflect the effects of the regionally-varying GIA signal, which shows up with its RSL fingerprints that stem from ice-sheets configurations and solid Earth response.

Holocene:
Continental, POS and marine overgrowth (serpulids) portions from submerged speleothems in Favignana (Mediterranean Sea) allow the construction of a sea level curve for the last 2000 years ( Figure 12).
The Yucatan speleothems record 4 hiati between 10.8 e 4,6 ka, which suggesthighfrequency sea-level fluctuations that might have been overimposed to the alocal falling sea-level in response to the GIA-driven ocean syphoning and continental levering (Mosely et al 2015).
The presence and/or absence of hiatuses in the Bahamas and Mediterranean speleothems indicate a sea level between -13.5 and -19 m. The higher elevation at Bermuda is related to the GIA.

MIS 5.3:
The presence and/or absence of hiatuses suggest a eustatic sea level between -12 and -6 m. At Bermuda, the a MIS 5.3 exceeds the +1.5 m of elevation because of local GIA.

MIS 5.5:
All the speleothemes record hiatuses between 117 and 148 ka. In particular, the -6m elevation of the Yucatan speleothem, constrains the end of the MIS 5e local highstand. Similarly, the elevations of -18.1 and -14 m constrain, respectively at 148 and 134 ka, the local RSL rise at Bahamas before the actual MIS 5e highstand.
The maximum vertical differences between the GIA-modulated MIS 5.5 highstands from the different Caribbean sites, have been quantified on the order of ~8 m, typical of intermediate-field regions. In the Mediterranean Sea, instead, the vertical differences are less than 2.5 m, thus showing that this is an intermediate-to-far-field region [59]. The Bahamas flowstone at -14 m does not preserve any hiatuses. The Argentarola stalagmite shows marine conditions at -18.5 m, but the local sea leve did not exceed -18.0 m -18. Therefore, we argue that the global sea level was between 14 and 18 m b.m.s.l.

MIS 7.5:
The observed highstand in the DWBAH flowstone is certainly higher than -18 m and presumably lower than -12 m, thus in agreement again with the global curve.
We show that the GIA-driven RSL variability at the sites under consideration increases during the interstadials and glacials. We arguet that this stems from the contribution for the ice-loading-induced crustal and geoidal deformations, which maintain strong spatial RSL gradients.