Submerged Speleothems and Sea Level Reconstructions: A Global Overview and New Results from the Mediterranean Sea

Abstract

This study presents a global overview of the submerged speleothems used to reconstruct paleo sea levels and reports new results from two stalactites collected in the Mediterranean Sea. Coastal cave deposits significantly contributed to the understanding of global and regional sea-level variations during the Middle and Late Quaternary. The studied speleothems cover the last 1.4 Myr and 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. The results indicate that submerged speleothems represent extraordinary archives that can provide 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 ¹⁴C ages of the new speleothems provide results of great interest for relative sea-level changes over the last 1000 years.

1. Introduction

The study of submerged speleothems in coastal caves significantly contributed to explaining past sea-level variations for the last 1.4 Myr, especially for the climatically-important 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 currently below sea level and 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 variation including shallow-water corals, lagoonal and marine cores, tidal notches, and archaeological remains that are well connected to sea-level fluctuations. However, most of these indicators exhibit large uncertainties in age determination due to sample alteration and/or in absolute elevation (e.g., paleowater depth for corals). Submerged speleothems that have been studied so far include (1) continental (vadose) carbonates, in a 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 analyzed for paleo sea-level reconstructions (Figure 1) [6].

Other speleothem-based sea-level records come from sites in the Caribbean region (Yucatán Peninsula and Cuba), The Bahamas, Bermuda, and Japan (Minami Daitō island) (Figure 2 and Figure 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 from the Ustica and Favignana islands in the central Mediterranean Sea. Published data for the Mediterranean Sea have recently been summarized in [6]. Some of the most important findings have come from speleothems retrieved from flooded coastal caves worldwide, often exhibiting hiatuses or marine overgrowths, and the submerged cave deposits in The Bahamas [3,4,7], Bermuda [8,9], Yucatán Peninsula [10,11], Cuba [12], and Japan [13] merit special attention. There have been exhaustive reviews of the use of sea-level indicators preserved in coastal caves including speleothems and cave sedimentary deposits [14,15].

2. Materials and Methods

2.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 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, which precisely mark the water level.

Several speleothems and flowstones were studied by [3,4] in The Bahamas 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 of 5 m 100 ky⁻¹ [3,21].

Speleothems and marine aragonite overgrowths were analyzed 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 below sea level, bsl) and includes 700 m of tholeiitic lava 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 forebulge 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 was 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 5 cm ky⁻¹ over the last 125 ky. The island is considered stable based on fossil tidal notches and dated corals [24].

Table 1 reports the basic information of the cave deposits from the western North Atlantic–Caribbean region and Minami Daito island. For further information regarding the cave morphology, sampling method used, date of sampling, and the error associated with the radiocarbon and U-Th ages, please 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].

2.2. Mediterranean Region

The coastline in the Mediterranean Sea (MS) extends for about 46,000 km, with approximately half 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-latitude 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 weather conditions, strongly control the tidal amplitude that varies along the Mediterranean coast. 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 of up to 1.80 m. In general, the Mediterranean tides have lower amplitudes compared to the other oceanic basins.

2.3. Ustica Island

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 about 60 km north of Palermo (northwest Sicily coast) and covers an area of 8.6 km², reaching a maximum elevation of 248 m a.s.l. The structural pattern of Ustica Island and its offshore area is characterized by ENE–WSW and N–S directed extensional and compressional fault systems [42,43,44]. The island is made up of alkali-basaltic to trachytic composition lava flow and pyroclastic products [45,46,47,48,49] that are locally covered by marine and continental sedimentary deposits [50]. The present morphology of Ustica Island and its offshore area results from the prolonged interplay between volcanic, tectonic, and eustatic processes [42,51], which have been active since the Pleistocene and contributed to the birth of the island between 750 and 130 ky ago, as demonstrated by stratigraphic and geochronological studies carried out on exposed rocks [49,52].

In the volcanic history of Ustica Island, there have been several overlapping cycles of marine ingression and regression, resulting from Middle–Upper Pleistocene glacio-eustatic 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 terrace [49,53]. The occurrence of marine terraces of variable age at decreasing elevations from the oldest to the youngest helps to define the amount of uplift recorded by the island, which must have been no more than 120 m, at least since the formation of the oldest level surface, which is 350 ky ago [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 (genus Chthamalus) bands radiocarbon dated to 1823 ± 104 yr Anno Domini AD and covered by a flowstone dated to 1894 ± 24 yr AD in the Grotta Segreta [51] reveals a coseismic event. An uplift of at least 40 cm is suggested that might be related to the strong earthquake in 1906, which forced the human population to evacuate the island. This example advises against the use of instrumental data alone for the long-term projections of future sea-level rise by 2100 (i.e., GPS data do not reveal vertical movements at Ustica).

2.4. Favignana Island

Favignana Island (Egadi Archipelago, Figure 1) is the closest island to Sicily. The archipelago was formed by carbonate and terrigenous Meso-Cenozoic deposits pertaining to different paleogeographic 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 to 10 m [56]. Favignana and Egadi islands are considered to have been stable or slightly uplifted areas for at least the last 125 ky [57,58,59]. Many outcrops of polygenic conglomerates contain fossil shells of Senegalese fauna including Persistrombus latus (previously Strombus bubonius). Uplifted Middle Pleistocene (MIS 9–11) forms and deposits reach an elevation of 40 m [56].

2.5. Speleothem Sampling in Favignana and Ustica Islands

A partly submerged stalactite was collected in 2016 from Stalattiti cave (Figure 1 and Figure 4) on 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 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 the meteorological website www.wunderground.com. The correction used 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 from Stalattiti cave (Figure 1 and Figure 6) on the north-east coast of Favignana Island using a hammer. At the time of the sampling, the deepest portion of the speleothem was 30 cm below 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 of the elevation measurements of the partly submerged speleothem relative to the sea level, which is in 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).

Figure 4. The stalactite sampled at Ustica: (a) Longitudinal section with sampling area of U\Th ages (refer to Table 2); (b) horizontal Scheme 1. US2, and US3, refer to Table 3; (c) the stalagmite sampled.

Figure 4. The stalactite sampled at Ustica: (a) Longitudinal section with sampling area of U\Th ages (refer to Table 2); (b) horizontal Scheme 1. US2, and US3, refer to Table 3; (c) the stalagmite sampled.

Figure 5. (a) Ustica and the location of the Grotta delle Stalattiti and Grotta Segreta caves. (b) Section of the Stalattiti cave. (c) Location of Stalattiti cave, Ustica town, and a well-carved marine terrace aged MIS 9 [53]. (d) Inside the cave at the sampling area, at sea level in an air bubble.

Figure 5. (a) Ustica and the location of the Grotta delle Stalattiti and Grotta Segreta caves. (b) Section of the Stalattiti cave. (c) Location of Stalattiti cave, Ustica town, and a well-carved marine terrace aged MIS 9 [53]. (d) Inside the cave at the sampling area, at sea level in an air bubble.

Figure 6. Section of the stalactite sampled at Favignana Island: (1) Continental portion sedimented during MIS 2 (Table 2); (2) the area highlighted in green contains a narrow (3–4 mm) band of marine deposition (Serpulids); (3) the final portion of the speleothem is composed by a stratified white carbonate, typical of phreatic deposition.

Figure 6. Section of the stalactite sampled at Favignana Island: (1) Continental portion sedimented during MIS 2 (Table 2); (2) the area highlighted in green contains a narrow (3–4 mm) band of marine deposition (Serpulids); (3) the final portion of the speleothem is composed by a stratified white carbonate, typical of phreatic deposition.

2.6. U-Th and ¹⁴C Dating

Three fragments of spelean calcite (~ 300 mg) were collected within the vadose portion of the stalactite from Ustica Island (Figure 4 and Figure 5) and analyzed 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 ²²⁹Th, ²³³U and ²³⁶U, and 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, Lisle, IL, USA) and analyzed 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 ²³⁰Th/U ages were calculated from measured atomic ratios through iterative age estimation [61] using the ²³⁰Th, ²³⁴U and ²³⁸U decay constants of [62] and [63]. A correction was applied for the non-radiogenic (detrital) ²³⁰Th fraction using an initial ²³⁰Th/²³²Th activity ratio of 1.5 ± 0.75. However, due to the low ²³²Th concentration (<0.05 ng/g) and high ²³⁰Th/²³²Th ratio (>700), the detrital correction had a negligible effect (Table 2).

Table 1. Locations and geographic coordinates of the caves and speleothems reported in this study, with information on references, dating method, time range, elevation, GIA, Tectonic corrections and local tide (https://www.tide-forecast.com/).

Site	Coordinates (Lat °N, Long °E)	Sample Name	References	Dating Method	Time Range (kyr BP)	Depth (m)	Tectonic Corrections	GIA Corrections MIS 5.1 m	GIA Corrections MIS 5.3 m	Local Tide m
Yucatán Peninsula	20.215 −87.4094	Y01-7-2(S) Y01-10-2(S)   Y01-22-1(F) Y01-23-2(F)	10	U-Th (MC-ICPMS)	86.6–61.3 109.2–107.7 hiatus 75.8–70.9 95.3–68.9 115.2–93.8	−14.6 −11 −11 −11 −9.9 −4.9	No	−34 m	−32	0.3
Yucatán Peninsula	20.215 −87.4094		11	U-Th (MC-ICPMS)	Marine: 8.9, 8.6, 8.4, 6 Cont.: 10.9, 8.4, 7.2, 5.8	−6	No	-	-	0.3
Bermuda	32.300 −64.7331	QB	9	U-Th (MC-ICPMS)	68–146 79–75 (hiatus) 84–82 (hiatus) 105-92 (hiatus) 137–114 (hiatus)	+1.5	No	−23 m	−19	1.35
Japan (Minami Daito Island)	25.8440 −131.2301	-	13	14C-AMS	5.1–4.6	−0.35	No	-	-	1.5
Cuba, Matanzas (Cueva de Santa Catalina)	23.0859 −81.4184	-	12, 19	U-Th (MC-ICPMS)	MIS 11, MIS 9, MIS 7, MIS 5 75.2–363.5	+16	Yes	−32	−33	1.08
Bahamas Lucayan Caverns	26.13029 −77.17577	Flowstone DWBAH	3, 5	U-Th (TIMS)	38.8–337.5 5 marine hiati	−15, −10	Yes	−27	−24	1.3
Bahamas Grand Bahamas South Andros. Sagittarius flowstone	24.1502 −77.7821	GB-89-25-5	4	U-Th (TIMS)	23.9–79.7; hiatus between 63.7 and 58.5. Tectonic was not considered	−18.1	Yes	−32	−27	1.35

Table 1. Locations and geographic coordinates of the caves and speleothems reported in this study, with information on references, dating method, time range, elevation, GIA, Tectonic corrections and local tide (https://www.tide-forecast.com/).

Site	Coordinates (Lat °N, Long °E)	Sample Name	References	Dating Method	Time Range (kyr BP)	Depth (m)	Tectonic Corrections	GIA Corrections MIS 5.1 m	GIA Corrections MIS 5.3 m	Local Tide m
Yucatán Peninsula	20.215 −87.4094	Y01-7-2(S) Y01-10-2(S)   Y01-22-1(F) Y01-23-2(F)	10	U-Th (MC-ICPMS)	86.6–61.3 109.2–107.7 hiatus 75.8–70.9 95.3–68.9 115.2–93.8	−14.6 −11 −11 −11 −9.9 −4.9	No	−34 m	−32	0.3
Yucatán Peninsula	20.215 −87.4094		11	U-Th (MC-ICPMS)	Marine: 8.9, 8.6, 8.4, 6 Cont.: 10.9, 8.4, 7.2, 5.8	−6	No	-	-	0.3
Bermuda	32.300 −64.7331	QB	9	U-Th (MC-ICPMS)	68–146 79–75 (hiatus) 84–82 (hiatus) 105-92 (hiatus) 137–114 (hiatus)	+1.5	No	−23 m	−19	1.35
Japan (Minami Daito Island)	25.8440 −131.2301	-	13	14C-AMS	5.1–4.6	−0.35	No	-	-	1.5
Cuba, Matanzas (Cueva de Santa Catalina)	23.0859 −81.4184	-	12, 19	U-Th (MC-ICPMS)	MIS 11, MIS 9, MIS 7, MIS 5 75.2–363.5	+16	Yes	−32	−33	1.08
Bahamas Lucayan Caverns	26.13029 −77.17577	Flowstone DWBAH	3, 5	U-Th (TIMS)	38.8–337.5 5 marine hiati	−15, −10	Yes	−27	−24	1.3
Bahamas Grand Bahamas South Andros. Sagittarius flowstone	24.1502 −77.7821	GB-89-25-5	4	U-Th (TIMS)	23.9–79.7; hiatus between 63.7 and 58.5. Tectonic was not considered	−18.1	Yes	−32	−27	1.35

Table 2. Results of the U/Th dating of the vadose speleothem (Ustica 1, 2, 3; Figure 5). ²³⁸U and ²³²Th concentrations were determined using the enriched ²³⁶U and ²²⁹Th isotopes, respectively. δ²³⁴U_m = {[(²³⁴U/²³⁸U)_sample/(²³⁴U/²³⁸U)_eq] − 1} × 1000, where (²³⁴U/²³⁸U)_sample is the measured atomic ratio, (²³⁴U/²³⁸U)_eq is the atomic ratio at secular equilibrium, and δ²³⁴U₍₀₎ corr is the initial value and was calculated by the following equation: δ²³⁴U₍₀₎ = δ²³⁴U_m exp(δ₂₃₄t), where t is the corrected age in years and δ₂₃₄ is the decay constant for ²³⁴U. * Ages corrected for detrital ²³⁰Th using an inital ²³⁰Th/²³²Th activity ratio of 0.85 ± 0.36 [64].

Sample	(238U)	(232Th)	δ234Um	(230Th/238U)act	(230Th/232Th)act	Age	Age corr
	(µg/g)	(ng/g)	(‰)			(kyr)	(kyr BP) *
Ustica-1	1.194 (0.001)	0.0483 (0.0005)	67.4 (0.6)	0.00939 (0.00007)	702.7 (5.4)	0.97 (0.01)	0.89 (0.01)
Ustica-2	1.143 (0.001)	0.0340 (0.0004)	67.5 (0.8)	0.01024 (0.00023)	1041.5 (24.0)	1.05 (0.03)	0.98 (0.03)
Ustica-3	1.543 (0.002)	0.0236 (0.0005)	70.5 (0.9)	0.00989 (0.00021)	1944.3 (42.3)	1.01 (0.02)	0.95 (0.02)

Table 2. Results of the U/Th dating of the vadose speleothem (Ustica 1, 2, 3; Figure 5). ²³⁸U and ²³²Th concentrations were determined using the enriched ²³⁶U and ²²⁹Th isotopes, respectively. δ²³⁴U_m = {[(²³⁴U/²³⁸U)_sample/(²³⁴U/²³⁸U)_eq] − 1} × 1000, where (²³⁴U/²³⁸U)_sample is the measured atomic ratio, (²³⁴U/²³⁸U)_eq is the atomic ratio at secular equilibrium, and δ²³⁴U₍₀₎ corr is the initial value and was calculated by the following equation: δ²³⁴U₍₀₎ = δ²³⁴U_m exp(δ₂₃₄t), where t is the corrected age in years and δ₂₃₄ is the decay constant for ²³⁴U. * Ages corrected for detrital ²³⁰Th using an inital ²³⁰Th/²³²Th activity ratio of 0.85 ± 0.36 [64].

Sample	(238U)	(232Th)	δ234Um	(230Th/238U)act	(230Th/232Th)act	Age	Age corr
	(µg/g)	(ng/g)	(‰)			(kyr)	(kyr BP) *
Ustica-1	1.194 (0.001)	0.0483 (0.0005)	67.4 (0.6)	0.00939 (0.00007)	702.7 (5.4)	0.97 (0.01)	0.89 (0.01)
Ustica-2	1.143 (0.001)	0.0340 (0.0004)	67.5 (0.8)	0.01024 (0.00023)	1041.5 (24.0)	1.05 (0.03)	0.98 (0.03)
Ustica-3	1.543 (0.002)	0.0236 (0.0005)	70.5 (0.9)	0.00989 (0.00021)	1944.3 (42.3)	1.01 (0.02)	0.95 (0.02)

Table 3. Radiocarbon data for Favignana Island stalactite cave. Column 1: * are from [65]. Column 4 * unavailable error bar. The obtained conventional radiocarbon ages were calibrated by using the Calib 8.10 software [66].

Sample	Corrected Altitude (cm)	Radiocarbon Age (BP)	δ13C (‰)	Kind of markers	Cal BP 1 Sigma
LTL15761A	−32 ± 5	1912 ± 45	+2 ± 0.7	S2 serpulids	1485–1304 1395 ± 45
LTL15762A	−35 ± 15	806 ± 45	+3 ± 0.5	S3 phreatic	735–675 705 ± 30
LTL15762B	−32 ± 5	19,886 ± 100	−29.5 ± 0.5	S1 continental	24,013–23,795 23,903 ± 109
R–2741 *	−30 ± 10	531 ± 32	+1.40 *	Cala Mancina St Vito Lo Capo	232–106 169 ± 63
R–2742 *	−30 ± 10	527 ± 36	+1.33 *	Cala Mancina St Vito Lo Capo (Tp)	232–94 163 ± 69
R–2764 *	−30 ± 10	554 ± 35	+1.31 *	Calazza (Pa)	251–127 189 ± 62
R–2580 *	−40 ± 10	784 ± 37	+1.76 *	Barcarello (Pa)	412–146 279 ± 133

Table 3. Radiocarbon data for Favignana Island stalactite cave. Column 1: * are from [65]. Column 4 * unavailable error bar. The obtained conventional radiocarbon ages were calibrated by using the Calib 8.10 software [66].

Sample	Corrected Altitude (cm)	Radiocarbon Age (BP)	δ13C (‰)	Kind of markers	Cal BP 1 Sigma
LTL15761A	−32 ± 5	1912 ± 45	+2 ± 0.7	S2 serpulids	1485–1304 1395 ± 45
LTL15762A	−35 ± 15	806 ± 45	+3 ± 0.5	S3 phreatic	735–675 705 ± 30
LTL15762B	−32 ± 5	19,886 ± 100	−29.5 ± 0.5	S1 continental	24,013–23,795 23,903 ± 109
R–2741 *	−30 ± 10	531 ± 32	+1.40 *	Cala Mancina St Vito Lo Capo	232–106 169 ± 63
R–2742 *	−30 ± 10	527 ± 36	+1.33 *	Cala Mancina St Vito Lo Capo (Tp)	232–94 163 ± 69
R–2764 *	−30 ± 10	554 ± 35	+1.31 *	Calazza (Pa)	251–127 189 ± 62
R–2580 *	−40 ± 10	784 ± 37	+1.76 *	Barcarello (Pa)	412–146 279 ± 133

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 using a diamond blade. The extracted portion of the sample was analyzed using an optical microscope in order to identify and mechanically remove adhering particles or unwanted fractions of the matrix. The sample was then attacked with H₂O₂ in order to remove the uppermost layer and then hydrolyzed under vacuum to CO₂ by using HCl. Extracted CO₂ was then catalytically 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 the machine and processing background, the conventional radiocarbon age of the samples was calculated. For continental samples, conventional radiocarbon ages were calibrated using the INTCAL20 calibration curve, while marine samples were calibrated using the MARINE20 [68] curve and a local ΔR = −88 ± 50 y (

Table 4. Radiocarbon data of the outer overgrowth of Ustica stalactite. The Note column shows the sample name used in Figure 5b.

Sample in Section	Sample	Radiocarbon Age (BP)	δ13C (‰)	Note
1	LTL20043A	2231 ± 45	−18.3 ± 0.6	US1
2	LTL20044A	1954 ± 45	−17.1 ± 0.7	US2
3	LTL20045A	1929 ± 45	−23.7 ± 0.4	US3

, Figure 5) [69].

2.7. 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. As 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 from the contribution of ocean mass variation in response to continental ice mass variation, 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 vice versa, solid earth deforms in order to restore the isostatic equilibrium under a different surface loading setting. Solid earth deformations and ice masses behave as density anomalies, thus affecting the vertical position of the mean sea surface, which is an equipotential surface of gravity (geoid). This process, known as glacial- and hydro-isostatic adjustment (GIA), is responsible for regionally varying RSL changes that are modulated in time with the viscous flow of mantle material.

As a result of mass conservation, the ocean average RSL change corresponds to the global mean, or eustatic, sea-level change. Based on the geographical distance from the ice sheets, local RSL change can either be very close to (ice distal or far-field regions), or significantly different from the eustatic sea-level change (ice proximal areas). The largest deviations from the eustatic signal are found where ice sheets grow and shrink. In particular, during glacial periods and beneath the growing ice sheets, the crust subsides and forces the upper mantle material to flow outwards, thus causing 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 top of the forebulge, the local RSL drop is much larger than the eustatic equivalent 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 areas experience strong subsidence that, combined with 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 far-field equatorial 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 change. Due to mass conservation, ocean water flows towards the subsiding peripheral forebulges, which locally results in an RSL drop that follows an initial highstand. This process is known as ocean syphoning. Interestingly, some of the Caribbean and Yucatán speleothems are in areas where complex forebulge ocean syphoning dynamics operate.

Therefore, RSL indicators provide constraints on the chronologies of past ice sheets and on the relevant earth rheological parameters that regulate the deformation of the planet in response to surface ice and water loading.

To compute local ice-driven relative sea-level change, the gravitationally self-consistent sea-level equation (SLE; 1–4) must be solved for a prescribed ice sheet 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 the solid earth surface. The difference between the two is the relative sea-level (RSL) change. We employed SELEN [70,71,72] to solve the SLE and retrieve the local GIA-modulated RSL curves at each site. The most 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 assumed a self-gravitating, rotating, spherically symmetric, radially stratified, deformable but not compressible Earth model. The latter is one-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 was characterized by linear Maxwell viscoelastic rheology. The mantle can be discretized into a number of layers.

First, we discretized the VM2 mantle viscosity profile [73] over a four-layer model where the upper mantle (UM), lower upper mantle (LUM), transition zone (TZ) and lower mantle (LM) are characterized by a Maxwell viscoelastic rheology and laterally uniform viscosity (

Table 5. Elastic lithosphere thickness and viscosity of the four-layer VM2 profile.

VM2	LT (km)	UM (× 1021 Pa⋅s)	LUM (× 1021 Pa⋅s)	TZ (× 1021 Pa⋅s)	LM (× 1021 Pa⋅s)
	90	0.67	0.44	0.46	2.53

and

Table 6. Elastic lithosphere thickness and viscosity of the three-layer MVP 1–3 profiles.

MVPs	LT (km)	UM (× 1021 Pa⋅s)	TZ (× 1021 Pa⋅s)	LM (× 1021 Pa⋅s)
MVP1	100	1	1	2
MVP2	100	0.5	0.5	5
MVP3	100	0.25	0.5	10

). We combined this with a lithosphere thickness of 90 km. We used this as a reference model to be combined with three ice sheet chronologies.

We employed the following ice-sheet models:

• ICE-5G [73]: This global model describes the ice-sheet thickness variation over North America, Eurasia, Greenland and Antarctica during the last 123 ka. Geological and modern geodetical observations were used to constrain the ice thickness chronology between 26 ka and the present day [6]. Through an iterative procedure, the SLE was solved for an a priori ice-sheet configuration and a prescribed fixed solid earth model (mantle viscosity profile and lithosphere thickness). The predicted RSL elevations and modern vertical rates were compared with paleo RSL indicators and measured crustal and geoidal vertical velocities, respectively. The differences between predictions and observations were then used to modify the ice-sheet models until a satisfactory solution is found. The ice-sheet volumetric evolution between 123 and 26 ka was tuned to the δ¹⁸O curve [74] and loosely constrained by surface geological evidence that defines the ice margins. Here, we combined two ICE-5G chronologies 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 the present-day mean sea level of 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 is a modified version of ICE-5G. Here, we applied the same repetition as in the 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 δ¹⁸O stack [74] is decoupled into 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 to include all the GIA feedbacks. ANICE-SELEN is not yet constrained by geological or instrumental data. Here, we used 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 variation of the three ice-sheet models is dependent, to a certain degree, on the benthic oxygen curve [74].

3. Results

This section is divided into two parts: The first reports published data in the Western North Atlantic–Caribbean region and Japan (Section 3.1), while the second reports new data from the MS (Section 3.2).

3.1. Review Data

3.1.1. Bermuda

The Western North Atlantic–Caribbean region, and Bermuda in particular, is strongly affected by the vertical motions of the periglacial forebulge that formed as a result of upper mantle flow in response to the expansion of the Laurentide ice sheet during glacial periods.

A previous review [9] summarized literature focused on the Bermuda area and the NE coast of the USA. The authors also published the results of a new stalactite sampled in a cave at +1.5 m asl in the age range from 137 ka to 71 ka. Four hiatuses are clearly visible in a continuous deposition on the speleothem (Table 1). A first hiatus (MIS 5.5) occurs from 137.0 ka to 114.7 ka, and a second from 105.9 ka to 92.0 ka corresponds with MIS5.3. The third hiatus ranges from 84.3 ka to 82.0 ka and corresponds with MIS5.1, and the presence of an additional hiatus from 79.9 ka to 75.0 ka corresponds with MIS5.1bis.

U-Th age uncertainties range between ±7 ka and ± 1 ka. These data support previous controversial indications that Bermudian RSL was significantly higher than the RSL at other locations during MIS5.5 and MIS 5.3 (Figure 2, Table 1).

Andros Island (The Bahamas)

A previous study [4] examined 116 speleothems that were collected from 22 caves that are part of a wide cave system in Andros Island (The Bahamas) during several expeditions from 1984 to 1992 (Figure 2 and Figure 7).

Table 1 reports the results of one of the most interesting speleothems collected at −18.1 m and covering almost the entire glacial period. The speleothem was sampled from −53.6 m and −1813.8 m. Regarding the flowstone sequence, which was called GB-89-25-5A/ B/ C, at −18.1 m bsl (Table 1, Figure 1 and Figure 2) the authors stated that other samples grew continuously throughout this period. “We cannot, therefore, support the suggestions for high sea-stand events near present sea level between 93 and 15 ky”. The same speleothem shows a non-erosional hiatus between 63.7 ± 1.7 ky and 58.5 ± 0.7 ky corresponding to MIS3 (Table 1, Figure 6).

3.1.2. Lucayan Cavern DWBAH Flowstone (The Bahamas)

The Lucayan Cavern DWBAH flowstone can be considered a historical discovery. The results of the work were published in 1989 using precise U-Th TIMS dating [3]. After 32 years, it remains one of the most complete sequences. The flowstone started to form about 300 ky ago at −10 m and ceased its growth about 39 ky ago. The flowstone shows five hiatuses in accordance with MIS9, MIS7.5, MIS7.3, MIS5.5, and MIS5.3 [3,5]. The tectonic rate considered for The Bahamas area is 1 m every 50 ky, as suggested by [7] (Table 1, Figure 2 and Figure 7).

3.1.3. Santa Catalina Cave (Cuba)

The Santa Catalina cave is carved into a marine terrace, which occurs along the Cuban coast. It was suggested by [19] that a general slowly uplifting trend occurred during the Pleistocene, even if U-Th dating of these terraces is difficult to achieve considering local conditions. The minimum age of the terraces can be constrained by dating carried out on speleothems from coastal mixing (flank margin) caves formed during past sea-level highstands together with the carving of marine limestones (Table 1, Figure 2). The authors suggested that speleothems collected in the Santa Catalina cave (Matanzas) provided ages >360 ky and showed several cycles of subaerial and subaqueous corrosion and speleothem growth. This suggests that the cave was carved during the MIS 11 sea-level highstand or earlier. Some stalagmites grew from MIS 11 to MIS 8 and were submerged twice, once at the end of MIS 11 and again during MIS 9. Phreatic overgrowths (POS) covering the speleothems suggest the presence of anchialine conditions within the cave during MIS 5.5. The altitude of these speleothems (16 m asl) indicates a late Pleistocene uplift rate of less than 0.1 mm ky⁻¹.

3.1.4. Yucatán

A previous study [11] found and dated four marine levels between 8.9, 8.6, 8.4 6 ka cal BP that were collected on a stalactite sampled at −6 m in a cenote on the Yucatán Peninsula. The authors 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. These oscillations are presumably due to GIA movements since all the area is considered tectonically stable”. The 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 gives error bars of cm, especially considering the local microtidal amplitude of 0.3 m (Table 1, Figure 2).

Another study [10] provided an age for numerous speleothems sampled in different submerged caves and cenotes in the Tulum area (Table 1, Figure 2 and Figure 7) at −14.6, −11, −9.9 and −4.9 m. The results (summarized in Table 1) indicated that between 109 ky and 61 ky ago the highstand of MIS 5.1 did not reach the lowest altitude of the speleothems studied (−14.9 m bsl). With regard to the MIS 5.3 highstand, it did not reach an altitude of −11 m bsl. For MIS 5.5, the precise speleothem age measured at −4.9 m bsl provides a robust chronology for the timing of the fall at the last interglacial termination, or the 5.5–5.4 transition.

3.1.5. Minami Daito Island (Japan)

A previous study [13] researched the partially flooded Akiba Chiteiko cave, in a small atoll in the northern Philippine Sea (Figure 3). The sample shows a phreatic overgrowth on speleothem (POS). Sampling and dating the highest POS level, which corresponds to the high tide mean level (Figure 7), the authors suggest that the sea level reached its Holocene maximum highstand between 5.1 and 4.6 ka cal BP, after which it remained more or less stable until the present day, with a possible minor sea-level drawdown of about 35 cm. The sea level at Minami Daito did not have a distinct highstand in the Holocene, which coincides with the mid-Holocene maximum observed elsewhere in the Indo-Pacific basin.

3.2. New Data

3.2.1. Favignana Island

We sampled and analyzed a speleothem with continental, marine and phreatic layers from a small cave at Favignana Island (Egadi Archipelago, Sicily) (Figure 4 and Figure 5). The cave entrance is narrow, with a length of 14 m, an average width of about 3–4 m, and a depth of about 7 m at the entrance that reduces to about 4 m inside (Figure 8). The studied speleothem, more precisely a stalactite, was blade-shaped with a width of about 1 m and was collected at the bottom of the cave. Today, the cave is flooded and the sea level lies 0.35 m above the speleothem (Table 4).

The speleothem shows an emerged part made of brown colored subaerial carbonate facies, and a 6 cm white carbonate facies that suggests a phreatic genesis. Between the aforementioned facies, we also noted and analyzed (Figure 4) a marine overgrowth made up of serpulids. Radiocarbon ages of subaerial carbonate overgrowth provided an age of 23.9 ky cal BP, while the marine overgrowth of 1.3 ky cal BP and the last presently submerged phreatic layer provided an age of 0.7 ky cal BP.

3.2.2. Ustica

The submerged cave (location and sketches in Figure 5) has a small entrance at −13 m bsl. A 10 m × 10 m room with open-air at 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 also some plant roots. In the submerged portion, both on the cave walls and the submerged stalagmites and stalactites, Polychaeta overgrowths can be observed. In some cases, they appear to be very large (up to 7–8 mm in diameter; see Figure 4). The stalactite sampled around the present-day sea level was covered with living Polychaeta. Once dissected, it showed a dark brown central part due to continental/phreatic deposition that occurred in the subaerial environment (Figure 4). Microscope analysis showed that the external part of the speleothem, which had lighter tones, presents a mix of both marine (with polychaetes) and continental origin (Figure 5). The age of this stalactite was estimated to be about 1000 y (Table 2). Radiocarbon dating carried out on the outermost part of the speleothem showed an older age while maintaining a stratigraphically correct sequence (Table 3).

3.3. GIA Numerical Predictions

For the three ice-sheet 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 always below the eustatic curve throughout the last 250 ka (Figure 9a–c). The Mediterranean sites (Ustica and Favignana islands) 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 Yucatán, the local RSL variability was in the order of 20–30 m during glacials and interstadials, and decreased during the MIS 5.5. This implies that Bermuda and The Bahamas, being closer to the North American ice-sheet margins, are more sensitive to GIA-modulated vertical motions, which can be emphasized by specific viscosity profiles.

RSL Predictions vs. Observed Values (MIS 5.1–5.3)

Overall, the predicted RSL peaks are lower than the observed terrestrial limiting elevations (deposition) and the hiatuses and the marine transgressions inferred from the speleothems (Figure 10). The largest RSL variability was predicted at Bermuda (Figure 9 and Figure 10c) because of the relative proximity to the ice-sheet margins and peripheral forebulge. Here, the maximum predicted highstands at MIS 5.1 and MIS 5.3 are ≈20 and ≈15 m lower than the speleothem elevation, respectively. Furthermore, the predicted highstands at MIS 5.1 and 5.3 were ≈5 and ≈13 m higher than at Yucatán (Figure 10b). These differences are in line with the results of [9], albeit significantly smaller. We argue that this might be related to the very different ice-sheet models (ICE-1 and ICE-2) used by [9]. Moreover, our predictions at Bermuda are only slightly higher than at The Bahamas (Figure 10d), implying that both sites share the same GIA response for the three ice-sheet models and four mantle viscosity profiles.

4. Discussion

Considering the data reviewed in this paper and [6], the use of submerged speleothems appears to provide remarkable results for the global sea-level change history during the Holocene and for long-term sea-level change reconstruction during the Late to Middle Pleistocene, in particular for MIS 1, 3, 5.1, 5.3, 7.1, 7.2, 7.3 7.5, and 9 highstands (Table 1, Table 2, Table 3 and Table 4). Figure 11 and Figure 12 depict the results for submerged speleothems from the Caribbean and Japan, and some of the main results of the aforementioned review for the Mediterranean area [6] are compared with the global curve of [79]. In [6], a speleothem sampled in Italy (Custonaci) exhibited marine and continental levels and seems to be older, aged at MIS 25–37 (1.4 MA).

4.1. Favignana Stalactite

The Favignana speleothem provides an interesting result concerning sea-level change over the last 1300 years, as sea-level data are very scarce in tectonically stable areas of the Mediterranean Sea [80,81]. The submerged portion of the stalactite was truncated and bored by Lithophaga lithophaga. The cave (Figure 8) may have recently been enlarged by marine processes, such as storm waves, a hypothesis that is supported by the fact that almost all the other stalactites near the one we studied were truncated. Based on this, we hypothesize the following evolution:

(1)

During LGM, the sea level was far from Favignana, at about −129 m bsl [82], and the speleothem presumably formed a column (Figure 11) aged 23.9 ka BP (Table 3).

(2)

The sea level rose until 1395 years cal BP. The sea level stillstand at −0.32 m promoted the growth of colonies of serpulids and sponges that formed on the broken surface of the speleothem, as reported from the Argentarola cave by [38,83].

(3)

The millimetric layer of phreatic carbonate (Figure 11) shows that the aforementioned sea-level stillstand was up to −0.35 m, considering the deeper phreatic level. Radiocarbon dating of the most recent level is 705 ky cal BP (Table 3).

(4)

The phreatic level was then interrupted by marine ingression, so that the submerged part of the stalactite was bored by Lithophaga sp., while the external carbonate surface of the submerged part of the speleothem was 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 hypothesize (similar to the cave from Palma de Mallorca) that the thickness of the phreatic deposit (in the shape of a ball) is 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), and that it has a greater thickness in the depth of the phreatic, the speleothem was sampled on 23 September 2014.

In addition, the age correction due to dead carbon was 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 yr BP (calibrated 705 cal yr BP, sample S3) would become 478 yr BP (calibrated to 518 cal yr BP) with a dead carbon percentage of 4%. Even a DCP of 10% would be enough to make the sample modern. In other words, if the DCP percentage of “dead” carbon in the cave was 10%, the sample would correspond to the present day. This is obviously not reasonable, especially because the stalactite was submerged. This suggests that the effect of dead carbon is very low, or negligible, in the studied cave. The dating of 705 yr BP, without considering dead carbon effects, fits perfectly with the curve, and this cannot be considered a coincidence. The dating analysis therefore leads us to believe that this effect is completely negligible. Figure 12 shows the estimated stalactite age according to different ice models for Favignana (see caption from Figure 10 for ice model, and Table 3) and other geomorphological and geoarchaeological markers of the same time range are compared. The serpulids overgrowth dated back to 1395 cal BP and was above almost all other curves. The phreatic speleothem dated back to 705 years and was in agreement with the average of the other curves. We also added some archaeological and biological markers in the 2 ka time range (Figure 12). For example, a fish tank in Salerno at −1.4 m bsl was dated back to 1.9 ky in [84]. At about −0.40 m, we found a less precise marker (a walking surface of the basilica of San Nicola [81]) at −0.32 m in the Middle Ages. We also added some radiocarbon ages from Vermetid reefs [65] at −0.3 m, which dated back to 169, 163, 189 years cal BP, and another at −0.4 cm which dated back to 279 years cal BP. Even if we advance these hypotheses with few data, there appears to be a peak of the sea-level curve that corresponds with the Medieval Warm Period between 900 and 1200 y AD and a cooling corresponding with the Little Ice Age between 1550 and 1900 y AD. This trend is in agreement with what was found by a study in Israel [80].

4.2. Ustica Stalactite

The stalactite sampled at the present-day sea level was about 0.8 cm long, with a diameter ranging from 4 to 8 cm (Figure 4c). As shown in Table 2 and Table 3, this speleothem covers a very narrow time range of about 900 to 1000 years (Figure 4a). There is a high concretion rate, which is somewhat surprising, given that it is located in a volcanic area. However, the Ustica lava contains many fragments of carbonate platform that erupted together with the lava, and these carbonate outcrops are often found within lava outcrops of up to 30 m, in some cases also with a modern tidal notch [51]. The most recent portion of the speleothem was formed by large living serpulid tubes. In the portion between the continental part and the living polychaete worms, continental sediment showing a different color was observed (Figure 4b) together with porosity and polychaete tubes, and sponges incorporated in the sediment were clearly visible.

The mixed deposition of continental carbonate with serpulid tubes (Figure 5b) constitutes 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 final part is younger. Our explanation for this is that the dead carbon of the very ancient carbonate inclusions in the Ustica lava has 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 (on 17 September 2016) speleothem (on of about 85 cm (Figure 4c) provide further evidence that, over the last 1 ka, the island of Ustica was affected by positive and negative coseismic events. This consideration is confirmed by the very recent coseismic movement described at the nearby 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 marine terraces raised on the lava, with an uplift rate ranging from 0.13 to 0.24 mm yr⁻¹ (see also Figure 5) [53].

4.3. Speleothems from the Caribbean and Japan

U-Th dating of the speleothems sampled at Bermuda, The Bahamas and Yucatán caves (Figure 2, Table 1) provided detailed information on the timing, duration, and altitude of sea-level highstands during the late and middle Pleistocene.

4.3.1. MIS 2—Holocene

A previous study analyzed a peculiar speleothem that preserves four marine hiatuses between 10.9 and 58 ky cal BP [11]. The region where this occurred is considered to be tectonically stable and it is unlikely that GIA could have been responsible for these anomalous crustal and/or sea surface fluctuations.

Another study [13] focused on a cave on the island of Minami Daito in Japan, which was made of a phreatic overgrowth. The authors investigated the relationship between the highest POS level, which represents the maximum elevation reached by the paleo-groundwater level and the observed modern highest peak. Regardless of a relevant local mean tide fluctuation of 1.5 m, the authors found that at between 5.4 and 4.6 cal BP, the local sea level reached a maximum highstand of about 35 cm during the mid-Holocene. This can be explained by GIA-driven ocean syphoning.

4.3.2. MIS 3

Figure 13 and Figure 14 show that the GB-89-25-5c speleothem that was sampled at an elevation of −18.1 m bsl in The Bahamas [4] had 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) makes this a unique and important indicator from an almost 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.

Our GIA predictions for the three ice-sheet models and mantle viscosity profiles are characterized by lower RSL peaks at MIS 3 (Figure 9a–c). As the three ice-sheet models are all somehow tuned to the δ¹⁸O stack curve, we suggest that their volumes might be biased towards larger values, while the speleothems described here might require a volume reduction throughout the MIS 3, which would result in a much higher eustatic level. Interestingly, this is supported by the ice-sheet model reconstruction from [85], which, for the MIS 3, provides a eustatic peak of 28 m bsl. The authors used ice proximal glaciological and geological evidence to constrain 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 was therefore fully independent from far-field RSl observations and deep-sea sea-level proxies.

4.3.3. MIS 5.1

The Yucatán speleothems do not show any hiatuses during the MIS 5.1 highstand, indicating it did not exceed the elevation of −14.6 m bsl [10]. Similarly, the DWBAH flowstone from The 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 (The Bahamas, [4]), which therefore limits the maximum highstand to −18.1 m bsl.

Yucatán and The 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 Yucatán and The Bahamas (see Figure 10b,d) were lower than the measured elevations, which were also corrected for tectonics. A previous study [9] reported 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 close to the peripheral forebulge that surrounded the North American ice sheet (see Figure 10c and Figure 15). Bermuda, therefore, is significantly affected by the vertical motions of the forebulge, i.e., uplift during ice-sheet expansion and subsidence during ice-sheet retreat (Figure 15).

Our highest MIS 5.1 sea-level prediction in Bermuda reached −18 m, while lower elevations of −25 and −20 m were predicted at Yucatán and The Bahamas, respectively (Figure 10b,d). This confirms that Bermuda is more affected by the GIA-related forebulge processes than Yucatán and The Bahamas (Figure 15).

The elevations inferred from the Caribbean speleothems of The Bahamas and Yucatán (−13.5 to −18.1 m) are somewhat comparable to that observed at 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 [85].

Therefore, by removing the expected GIA contributions, we also hypothesize a eustatic elevation for the MIS 5.1 of approximately −17 m, obtained as average value between −18.8 and −13.5 m.

4.3.4. MIS 5.3

The speleothem from Bermuda which lies at +1.5 m asl preserves several hiatuses and was probably affected by forebulge processes (Figure 9a–c, Figure 10c, and Figure 16).

The DBWAH flowstone from The Bahamas at −13.0 m bsl shows an MIS 5.3 hiatus that cannot be explained by our GIA results. A hiatus was recorded in Yucatán, where the sea level might have been higher than −11 m after 109 ka, but less 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 9a–c, Figure 10a–d). Again, we argue that a reduction of ice-sheet volumes is required to shift the eustatic value upwards.

We suggest a eustatic highstand of between −12 and −6 m.

4.3.5. MIS 5.5

All the submerged speleothems from the Caribbean Sea show hiatuses at −6 m bsl at 117 ky (Yucatán), −18.1 m bsl at 148 kyrs (The Bahamas), and −14 at 134 kyrs of the flowstone DWBAH. The only difference between these was the timing of the transgression. In Bermuda, the hiatus of MIS 5.5 was between 114 kyrs and 137 kyrs.

4.3.6. Middle Pleistocene

The only date that was older than MIS 5.5 was provided by the DWBAH flowstone [3,5], which contained five marine hiatuses (MIS 9, 7.5, 7.3, 5.5, 5.3).

We can provide some considerations on differences in tide amplitude between the Caribbean Sea and MS. Table 1 shows the maximum tide amplitudes in Bermuda and The 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 (measured on Barbados). We suggest that the precision in elevation of the submerged speleothems is higher than corals, because speleothems are precise indicators with respect to the mean sea level for tides lower than 0.5 m (Table 1).

5. Conclusions

Submerged speleothems can grow for significantly long timespans and, therefore, can record multiple intervals during which deposition stops and/or erosion takes over in response to sea-level highstands. The alternation of continental depositional layers and hiatuses or marine layers records show long-term climate-related global ice-sheet 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 the Mediterranean Sea.

The longest time span was recorded in the DWBHA flowstone [3] at The Bahamas and covered the last 326 ka, with five hiatuses recorded. The Bahamas speleothem GB-89-25-5 covered the time between 80 and 24 ka and contained one hiatus. At Bermuda, the QB speleothem covered the time between 70 and 146 ka and recorded four hiatuses. Several speleothems from Yucatán [10] cover a time span ranging from 61 to 115 ka, with one hiatus. In the Mediterranean Sea, the ASI speleothem from Argentarola [86] covers the last 210 ka and recorded four hiatuses. The K-18 stalactite from Croatia [6,26] spanned 55 to 131 ka and contained one hiatus.

Overall, the discussed speleothems reflect the effects of the regionally varying GIA signal, which is indicated by its RSL fingerprints that stem from ice-sheet configurations and the solid earth response.

Holocene: continental, POS, and marine overgrowth (serpulids) portions from submerged speleothems in Favignana (Mediterranean Sea) facilitate the construction of a sea-level curve for the last 2000 years (Figure 12).

The Yucatán speleothems exhibit four hiatuses between 10.8 and 4.6 ka, which suggests high-frequency sea-level fluctuations that might have been superimposed on the local falling sea level in response to the GIA-driven ocean syphoning and continental levering [11].

Similarly, the phreatic speleothem from Minami Daito Island in Japan records a mid-Holocene sea-level highstand of 0.35 m, which is characteristic of the far-field Asian sites where the ocean syphoning process dominates [13].

MIS 3: the Bahamas elevation was −18.1, which is much higher than the typical global eustatic curves, but also higher than the observed value from the Mediterranean sites (from −29 to −19 m). This might be due to the different GIA responses of the two regions.

MIS 5.1: the presence or absence of hiatuses in The Bahamas and the Mediterranean speleothems indicate a sea level of between −13.5 and −19 m. The higher elevation at Bermuda is related to the GIA.

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

MIS 5.5: all the speleothems exhibited hiatuses between 117 and 148 ka. In particular, the −6 m elevation of the Yucatán speleothem constrains the end of the MIS 5e local highstand. Similarly, the elevations of −18.1 and −14 m at 148 and 134 ka, respectively, constrained the local RSL rise at The Bahamas before the actual MIS 5e highstand.

The maximum vertical differences between the GIA-modulated MIS 5.5 highstands from the different Caribbean sites were quantified in the order of ≈8 m, typical of intermediate-field regions. In the Mediterranean Sea, the vertical differences were less than 2.5 m, indicating that it is an intermediate-to-far-field region [59].

MIS 7.1: the flowstone at −14 m did not record any hiatuses, while the three speleothems from the Argentarola cave (Figure 13) showed a clearly visible serpulid overgrowth.

MIS 7.2: a clear highstand was recorded in the Argentarola cave. The sea level oscillated between −18.5 and −21 m, similar to the global curve that reaches −27 m.

MIS 7.3: the Bahamas flowstone at −14 m did not exhibit any hiatuses. The Argentarola stalagmite shows marine conditions at −18.5 m, but the local sea level did not exceed −18.0 m. 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 was certainly higher than −18 m and presumably lower than −12 m and is therefore in agreement with the global curve.

We show that the GIA-driven RSL variability at the sites under consideration increased during interstadial and glacial periods. We argue that this stems from the contribution of the ice-loading-induced crustal and geoidal deformations, which maintain strong spatial RSL gradients.