Challenges of initial Thorium and Approaches to Robust Speleothem Age Models: A case study from the Yucatán peninsula, Mexico

Speleothems, such as stalagmites and flowstones, are invaluable archives of past environmental and climatic conditions due to their layered growth and suitability for precise ²³⁰$\mathrm{Th}/\mathrm{U}$ dating. In tropical karst settings, however, elevated and variable detrital thorium contamination can severely compromise age accuracy, especially for young speleothems. Here, a suite of stalagmites from Áaktun Kóopo' Cave from the Yucatán Peninsula is used as a case study to develop and evaluate a multi-method framework for constraining initial thorium and constructing robust age–depth models. We combine extensive ²³⁰$\mathrm{Th}/\mathrm{U}$ dating with local isochron analysis, stratigraphic approaches, and annual $\mathrm{Sr}/\mathrm{Ca}$ layer counting to constrain elevated and initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios and resolve pronounced chronological inversions. High uranium concentrations in the speleothems (average ca. 1 ppm) allow precise measurements, yet we infer a high and largely unsystematic variability of high initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios in space and time, with values spanning between 4 and 68, thus far exceeding standard bulk earth values. This variability demonstrates that elevated and temporally variable initial Th may be common in tropical karst settings and that multiple, independent constraints on initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios are essential for reliable chronologies in such environments. Applying this framework to Áaktun Kóopo' Cave stalagmites yields internally consistent stalagmite age–depth models spanning the past 2.7 kyr, as well as evidence for earlier glacial and interglacial growth phases. These chronologies provide the basis for high-resolution, multi-proxy reconstructions of climatic and environmental changes in the northeastern Yucatán Peninsula during the entire era of Maya cultural evolution, and they illustrate how robust speleothem chronologies can be obtained even in cave systems affected by elevated and highly variable initial thorium.

1 Introduction

Reconstructing past environmental and climatic variability relies on archives that combine continuity, sensitivity to climate, and robust dating methods (Fairchild and Baker, 2012; Koltai et al., 2026; Medici et al., 2026). Speleothems fulfil these criteria, as their layered carbonate growth can be precisely dated using U-series disequilibrium dating (Wendt et al., 2021; Scholz and Hoffmann, 2008; Cheng et al., 2013). Constructing accurate age models based on a large number of individual ²³⁰$\mathrm{Th}/\mathrm{U}$ ages is fundamental for interpreting speleothem climate proxy records, as chronological uncertainties directly affect the absolute timing of events, their duration, and partially the magnitude (smoothing) of inferred climate variability (Comas-Bru et al., 2020; Moseley et al., 2016). However, ²³⁰$\mathrm{Th}/\mathrm{U}$ dating assumes a closed system for the exchange of uranium and its decay products upon speleothem formation, and it further presumes absence of initial ²³⁰Th. In such an ideal theoretical scenario the $\mathrm{Th}/\mathrm{U}$ age precision and accuracy will solely depend on the analytical precision and accuracy of isotope ratio measurements and isotope half-lives (Cheng et al., 2013). Any deviation of the carbonate precipitating environment from this ideal scenario may influence the age accuracy, either through U-series system opening via dissolution and reprecipitation of carbonate, or via a variable degree of initial Th contamination. To date, U-series system opening cannot be corrected for cave deposits due to the lack of constant boundary conditions such as a time-invariant U-isotope composition of drip water. Such U-series open system models only exist for well-preserved marine carbonates and assume a near constancy of seawater U-isotope ratios as additional constraint (Frank et al., 2006). For the presence of initial ²³⁰Th, however, several strategies exist to correct for its influence. In most cases ²³⁰Th is related to the detrital origin of ²³²Th assuming a constant detrital activity ratio for linear correction models applied to the measured ${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$ activity ratio. Appropriate correction techniques have to be applied, which are, however, associated with relatively large uncertainties. Consequently, high Th contamination and young ages of speleothems result in ²³⁰$\mathrm{Th}/\mathrm{U}$-age errors exceeding analytical precision by orders of magnitude (Ludwig and Titterington, 1994). This issue is particularly problematic for speleothems with low uranium concentrations and visible non-carbonate traces, since even small amounts of initial Th can lead to erroneously old age estimates (e.g., Faraji et al., 2023; Kerber et al., 2025; Fensterer et al., 2010).

A primary challenge with initial Th correction arises from the difficulty in accurately estimating the initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio of detrital material (further on referred to as R02i), which can vary depending on the characteristics of the host rock, soil, and chemical complexing agents during aqueous transport of Th (Ivanovich and Harmon, 1992). For geological material in ²³⁸U–²³⁴U–²³⁰Th secular equilibrium, and assuming a crustal ${}^{\mathrm{232}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$ activity ratio of 3.8 (Taylor and McLennan, 1985), the resulting atomic ${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$ ratio is $\mathrm{4.4}\pm \mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{6}}$ (with an assumed error of 50 %), equivalent to an activity ratio of ∼0.8. This ratio, widely utilized for detrital ²³⁰Th correction (Spötl et al., 2002), is commonly applied to estimate initial ²³⁰Th and correct ²³⁰$\mathrm{Th}/\mathrm{U}$ ages. This is particularly effective when absolute ²³²Th concentrations are low (< a few ng) and speleothems are older than several thousand years. However, tropical karst systems frequently exhibit elevated and variable R02i in drip waters, modern carbonates, and speleothems, sometimes up to 80 times higher than the commonly assumed value of ∼0.8 for the upper continental crust. For example, speleothem studies from Mexico report a wide range of R02i between 3.5 and 32 (Moseley et al., 2015; Kennett et al., 2022; Stinnesbeck et al., 2020). R02i ranging from 8 to 41 have been used to correct a Belizean stalagmite (Ridley et al., 2015), and values between 11 and 61 were applied to Puerto Rican stalagmites (Kerber et al., 2025; Vieten et al., 2024a, 2024b; Warken et al., 2020; Rivera-Collazo et al., 2015). Also at other sites including the Bahamas (Beck et al., 2001; Arienzo et al., 2015; Hoffmann et al., 2010; Richards and Dorale, 2003), Cuba (Fensterer et al., 2010), or in the tropical Pacific (Faraji et al., 2021; Partin et al., 2007; Carolin et al., 2013) elevated R02i were found. These observations thus demonstrate a vast range of initial Th variability, also within single cave sites. However, large (and poorly constrained) variations in R02i directly translate into additional, often dominant, age uncertainties, particularly for young and/or strongly contaminated speleothems. Classical correction schemes based on a single bulk earth detrital ratio can therefore produce stratigraphic age inversions and biased chronologies when applied in such settings.

In general, there are two main approaches to assess R02i in speleothems: (1) direct measurement of the contaminating phase, e.g., from cave drip water, soil, or zero-age carbonate (ZAC) (Kerber et al., 2025; Hu et al., 2008; Li et al., 2022; Wortham et al., 2022); and (2) estimation of R02i using linear binary mixing models or a priori assumptions to obtain continuous chronologies. In cases where direct measurement is not possible, the isochron method can serve as a first order alternative (Kaufman, 1993; Richards and Dorale, 2003; Frank et al., 2002; Ludwig and Titterington, 1994). However, for this technique to be effective, multiple sub-samples from the same lamina (or group of laminae) must exhibit varying ${}^{\mathrm{232}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$ ratios and follow first-order binary linear mixing behaviour (Wenz et al., 2016; Ludwig and Titterington, 1994). Another approach is the stratigraphic constraint method, where R02i is estimated by ensuring that the corrected ages increase steadily with distance from the top (dft) of the stalagmite (Roy-Barman and Pons-Branchu, 2016; Hellstrom, 2006). While monitoring detrital activity ratios in ZAC and drip waters provides useful information about initial activity ratios, it offers little insight into whether and how these ratios change over time. Overall, most studies still rely on only one of these methods, and the resulting uncertainty in initial Th variability is often poorly quantified. This highlights the need for more systematic assessments of initial Th when correcting ²³⁰$\mathrm{Th}/\mathrm{U}$ ages from such environments.

High-resolution stalagmite records from the Yucatán Peninsula (YP) have been analysed to explore hydroclimate variability and its intersection with Maya sociopolitical evolution (e.g., Kennett et al., 2012; Medina-Elizalde et al., 2010; James et al., 2026). However, in the northeastern YP, high-resolution records are few, with only one stalagmite record available, spanning 1037 BCE to 397 AD (Medina-Elizalde et al., 2016b). Critically, YP records often face dating challenges due to low uranium concentrations or high detrital contamination (e.g., Frappier et al., 2014; Webster et al., 2007; Akers et al., 2016, 2019; James et al., 2025; Medina-Elizalde et al., 2010; Kennett et al., 2022; James et al., 2026). These complexities make it difficult to reconstruct past climate and environmental changes with precision, obscuring the evaluation of their effects on Maya cultural evolution, particularly during the Terminal Classic Period (TCP; ∼800–1000 AD), a time of profound social and political upheaval marked by widespread population decline and the abandonment of many urban centres (e.g., Douglas et al., 2015).

This study focuses on constructing precise ²³⁰$\mathrm{Th}/\mathrm{U}$ chronologies for stalagmites from Áaktun Kóopo' Cave, located in the northeastern YP. Through a detailed study of (i) the speleothems' U and Th isotope compositions, (ii) their micro-facies and geochemistry, as well as (iii) their stable isotope compositions we establish multiple age-depth models that demonstrate active speleothem deposition over the past 2.7 kyr. These chronologies will allow for filling a notable gap in high-resolution paleoclimate records for the region and provide a foundation for future investigations into hydroclimatic and environmental changes, particularly during periods of profound Maya societal changes. Crucially, our study reveals highly variable initial ²³⁰Th coprecipitation within the karst system that impacts, and in some layers inhibits, high-precision ²³⁰$\mathrm{Th}/\mathrm{U}$ age determination, thereby limiting the accuracy of speleothem chronologies. The chemical transport processes of Th isotopes in meteoric water, however, remain poorly understood.

2 Material and Methods

2.1 Study site

Áaktun Kóopo' Cave is situated in La Estrella, a small Maya community in the northeastern YP (Fig. 1a), approximately 30 km northeast of Chemax and 7 km southwest of Xcan, at 20°48^′26.7" N, 87°44^′15.8" W. This region experiences a warm, humid tropical climate with highly seasonal rainfall (wet season June–November), high year-round humidity and frequent tropical storms and hurricanes (Romero and Alfaro, 2024; Medina-Elizalde et al., 2016a) and the soils are typically thin and variable, with pockets rarely reaching 0.5 m between areas of exposed bedrock (Lases-Hernandez et al., 2019).

Áaktun Kóopo' Cave is developed within the low-relief Yucatán carbonate platform, which consists of thick Cretaceous and Tertiary shallow-marine limestones, forming a great limestone plain with a shallow water table and almost no surface rivers (Ward et al., 1995; Perry et al., 2009; Smart et al., 2006). The cave comprises more than 2 km of horizontal passages and at least nine water bodies (Fig. 1b), with seasonally fluctuating water levels that can locally flood parts of the cave (Thomas, 2011; Tec, 2021). Cave sediments are dominated by red clay (kankab) that ranges from profound to muddy conditions. Archaeological evidence inside the cave includes Maya paintings, polychrome pottery, and skeletal remains of 23 individuals (juveniles and adults), dating to the Late Preclassic (−250 to 300 AD) and Early Classic (300–600 AD) periods (Gómez, 2020). The Maya also modified parts of the cave by constructing walls using limestone blocks, speleothems, and soil to create chambers (Tec, 2021). The presence of imported materials such as flint and obsidian (Gómez, 2020) suggests that members of the local elite performed mortuary and ritual practices in the cave. In addition to these archaeological features, Áaktun Kóopo' Cave is decorated with numerous speleothem features, including stalactites, stalagmites, flowstones, and draperies, making it an important site for both archaeological and paleoclimate research.

Figure 1 (a) Location of the Áaktun Kóopo' Cave (yellow asterisk) in the northeastern Yucatán Peninsula, Mexico, relative to other Maya archaeological settlements in the area. Orange diamonds represent large settlements while large (small) triangles represent intermediate (smallest) ones, based on Witschey and Brown (2010). Well-known settlements are e.g., Coba, Punta Laguna, El Naranjal, and Ek Balam. (b) Map of Áaktun Kóopo' Cave with sample locations. Red stars show locations of stalagmites that have been extracted from the cave, while brown stars represent locations where samples were drilled or scraped from stalagmites. Blue symbols mark sites of drip water collection. Map taken from Thomas (2011) and modified after Gómez (2020) and Tec (2021). Green lines indicate locations of Maya walls.

2.2 Sample material

During three cave visits in 2018, 2022, and 2023, samples for ²³⁰$\mathrm{Th}/\mathrm{U}$ dating were collected from 15 different speleothems inside Áaktun Kóopo' Cave, including stalagmites, flowstones, and modern carbonate deposits. In addition to in situ drill cores from various speleothems, ten stalagmites, ranging from 6.5 to 34 cm in length, were extracted completely for more detailed analysis (Fig. S1 in the Supplement). For six of these stalagmites (E0, E1, E2, E7, E8, E23-3) their drip sites were still actively dripping at the time of collection. Stalagmites E0, E1, E7, E8, and E23-3 were all collected from the same cave chamber in the northwestern part of the cave, with E0, E1, E7 and E8 even from the same two square metre area (Fig. 1b). Interestingly, both E0 and E7 exhibited an L-shaped, or “double” formation (Fig. S1). This morphology occurs when a stalagmite falls over, and continued carbonate precipitation causes growth to proceed at a perpendicular angle to the original axis, resulting in the observed L-shape. Notably, E0 was not calcified to the ground but rather embedded in soft sediment, which may have contributed to its collapse. During the extraction, a third stalagmite (E0-A) was uncovered beneath 10 cm of sediment. Stalagmite E2 was collected from the central part of the cave, while stalagmite E4 was retrieved from the northeastern section, closer to the cave's entrance. Stalagmite E3 was collected from the southern portion, located in front of a large Maya wall, from which a small stalagmite (EMW5) detached and was collected as well. For the follow-up analyses, the stalagmites were cut along their presumed growth axis into slices.

Furthermore, we collected three drip water samples (EDW-1, EDW-2 and EDW-3) from different parts of the cave in 2022 (Fig. 1b) to evaluate the R02i in drip water. EDW-3 consists of an integrated water sample collected from different drips (soda straws) with low drip rates within an area of about 1 m² within the chamber in the western part of the cave where most of the analysed stalagmites were extracted (Fig. 1b). After drip water collection, the plastic bottles (125 mL) were immediately sealed with Parafilm for transportation and storage upon measurement. From the floor below the drip site EDW-1, we collected a small sample of what appeared to be modern carbonate precipitates.

2.3 Mineralogical characterization of stalagmites

2.3.1 XRD measurements

Slices of stalagmites from Áaktun Kóopo' Cave revealed several changes in colour throughout their growth (Fig. S1), which may indicate changes in mineralogy or crystal fabrics. To investigate these variations and determine the mineralogical composition of the stalagmites (e.g., calcite vs. aragonite), nine samples were analysed from visually distinct sections of stalagmite E0-C for XRD analysis (Fig. S2a). Additionally, a sample was collected from the bottom portion of stalagmite E1, where the fan-like crystal texture and silky gloss suggested aragonite as the dominant polymorph (Figs. S1 and S2). Each sample consisted of approximately 1 g of carbonate powder that was drilled using a Proxxon machine equipped with a diamond-coated stainless-steel drill bit. XRD measurements were performed at the Institute of Earth Sciences, Heidelberg University, using a Bruker D8 ADVANCE Eco diffractometer. The measurements were conducted with the following parameters: Cu Kα radiation source, 30 kV voltage, 33 mA current, Ni filter, SSD160 detector. Samples were analysed in rotating sample holders over an angular range of 2θ from 5 to 70°, with increments of 0.019° and a counting time of 1 s per step. Peak positions and intensities were processed using Diffrac.Suite EVA software (Bruker).

2.3.2 Thin section inspection

Thin sections have been used to study the crystal fabrics in further detail, i.e. the individual expression of mineral phases, direction, non-carbonate contamination and/or signs of dissolution (e.g., Frisia et al., 2002; Perrin et al., 2014). Thin sections were prepared from stalagmites E0-B, E0-C, E4, E8, and the upper section of E1. These sections were analysed under a polarized transmitted-light microscope (KEYENCE VHX-6000) at the Institute of Earth Sciences, Heidelberg University.

2.4 Geochemistry

2.4.1 In situ strontium-calcium analysis

Laser ablation ICP-MS analyses of selected stalagmites and for $\mathrm{Sr}/\mathrm{Ca}$ ratios followed established protocols (Warken et al., 2021; Jochum et al., 2012; Schorndorf, 2024) and were performed using a 193 nm ArF excimer laser (NWR193UC by New Wave Research) coupled to an inductively coupled plasma quadrupole mass spectrometer (Thermo Fisher iCAP-Q) at the Institute of Environmental Physics, Heidelberg University. Line scans were conducted along the growth axis of the stalagmite slices using a rectangular spot size of 25×100 µm. The repetition rate was set to 20 Hz, with scan speeds of 5 µm s⁻¹ for stalagmites E1 and E8, and 10 µm s⁻¹ for stalagmite E23-3. To eliminate potential surface contamination, a pre-ablation step was performed along each scan path using a circular spot size of 150 µm at a scan speed of 100 µm s⁻¹ and a repetition rate of 10 Hz. Background counts were measured with the laser in the off mode and subtracted from the raw data. Data processing involved an outlier correction using a floating 1.5 interquartile range (IQR). To account for matrix effects, blank-corrected count rates of ⁴³Ca and ⁸⁸Sr were normalized to the ⁴⁴Ca signal, which served as an internal standard for the total ablation efficiency. External calibration was performed using silicate glass NIST SRM 612, with reference values provided by Jochum et al. (2011). Drift corrections were applied through linear interpolation between two standard measurements taken before and after each ablation path. The resulting $\mathrm{Sr}/\mathrm{Ca}$ elemental ratios are presented as mass concentration ratios.

2.4.2 Stable isotope analysis

For stalagmites E0-C, E1, and E23-3, approximately 50–90 µg of powdered samples were drilled using a Micromill (2010 A Sherline) equipped with a micro milling cutter (TiSiN) of 500 µm in diameter. Most samples were analyzed for their stable isotope composition of oxygen (δ¹⁸O) and carbon (δ¹³C) at the Institute of Earth Sciences, Heidelberg University, using a ThermoFinnigan MAT253Plus gas source mass spectrometer equipped with a Thermo Fisher Scientific Kiel IV carbonate device. Quality control is based on the analysis of an in-house standard (Solnhofen limestone, δ¹³C${}_{\mathrm{VPDB}}=+\mathrm{1.38}\pm \mathrm{0.03}$ ‰ and δ¹⁸O${}_{\mathrm{VPDB}}=-\mathrm{4.59}\pm \mathrm{0.06}$ ‰) calibrated to the reference material IAEA-603 (calcite; δ¹³C${}_{\mathrm{VPDB}}=+\mathrm{2.46}\pm \mathrm{0.01}$ ‰ and δ¹⁸O${}_{\mathrm{VPDB}}=-\mathrm{2.37}\pm \mathrm{0.04}$ ‰). A part of the samples of E0-C was analysed at Elemtex Ltd., Cornwall, UK, using an isotope ratio mass spectrometer (IRMS) (Thermo Scientific Delta V) coupled with a Gasbench II. Here, quality control is performed using a Carrara marble standard (δ¹³C =2.10 ‰, δ¹⁸O $=-\mathrm{2.01}$ ‰) and a second in-house calcite standard (δ¹³C =2.89 ‰, δ¹⁸O $=-\mathrm{6.15}$ ‰), calibrated directly against NBS18 and NBS19.

External reproducibility at both labs is better than 0.04 ‰ for δ¹³C_VPDB and 0.07 ‰ for δ¹⁸O_VPDB at the 1σ confidence level, and all δ¹³C and δ¹⁸O values are reported relative to Vienna Pee Dee Belemnite (VPDB) standard.

2.5 Chronology

2.5.1 ²³⁰$\mathrm{Th}/\mathrm{U}$ dating

A total of 157 samples from Áaktun Kóopo' speleothems as well as the three drip water samples were analysed using high-precision ²³⁰$\mathrm{Th}/\mathrm{U}$ dating at the Institute of Environmental Physics, Heidelberg University. For the analysis, small sample plates weighing approximately 100 mg were cut perpendicular to their growth direction using a diamond wire saw and subsequently leached in weak acid, weighed, dissolved in ultra-clean nitric acid and spiked with our in-house triple-spike solution (“TriSpike”, ²³³U, ²³⁶U and ²²⁹Th) (Kerber et al., 2023). For drip water analysis, about 15 to 30 mL of sample was weighed and spiked with “TriSpike”, acidified with ultra-clean nitric acid, subsequently dried down and re-dissolved. Overall, the preparation involving sample leaching, dissolution, spiking, and wet column extraction chromatography using Eichrom UTEVA resin followed established protocols (Wefing et al., 2017; Matos et al., 2015; Douville et al., 2010). Chemical blanks are regularly monitored and commonly below 0.4 fg for ²³⁴U and 0.04 fg for ²³⁰Th (Kerber et al., 2023, 2025).

Isotope measurements were performed with a multi-collector inductively coupled plasma source mass spectrometer (MC-ICP-MS, Thermo Fisher Neptune plus) using a semi-static multi-cup setting as described by Kerber et al. (2023). Mass spectrometry techniques and data treatment procedures are detailed in Kerber et al. (2023, 2025). Data analysis was performed using a Python-based GUI application (Kerber et al., 2025). Ages were calculated using the half-lives provided by Cheng et al. (2013) for ²³⁰Th and ²³⁴U as well as by Jaffey et al. (1971) for ²³⁸U. All uncertainties are reported at the 2σ level, excluding errors associated with half-life values.

2.5.2 Radiocarbon dating

To constrain recent carbonate deposition, the tops of presumably modern speleothems were analysed using radiocarbon dating. Small calcite samples (10–20 mg) were cut from the tops of stalagmites E0-C, E1, and E8 using a diamond wire saw at the Institute of Environmental Physics, Heidelberg University. To avoid contamination by ambient air, the samples were leached in 4 % hydrochloric acid prior to hydrolysis and graphitization. A detailed description of the sample preparation routine is provided in Therre et al. (2021). The resulting graphite–iron compounds were measured using a 200 kV tandem mini carbon dating system (MICADAS) accelerator mass spectrometer (AMS) at the Curt-Engelhorn Center for Archaeometry in Mannheim, Germany (Kromer et al., 2013; Synal et al., 2007). Details regarding long-term blank values and external standard reproducibility can be found in Beisel et al. (2025) and Therre et al. (2021).

2.5.3 Speleothem age-depth models

The age-depth models for selected stalagmites were constructed using a suite of established approaches: linear interpolation, linear regression, Bchron (Haslett and Parnell, 2008), Bacon (Blaauw and Christen, 2011; Blaauw et al., 2021), COPRA (Breitenbach et al., 2012), and StalAge (Scholz and Hoffmann, 2011). Model construction followed and extended the protocols of the Speleothem Isotopes Synthesis and Analysis (SISAL) working group (Roesch and Rehfeld, 2020), which is provided as a supplement to Comas-Bru et al. (2020); SISAL.AM. In case of identified hiatuses, the age-depth models were split at the depth of the growth interruption into separate sections. This approach improves model convergence and allows more methods to return successful chronologies, rather than including the hiatuses in the models. For each model, the median and the 2.5 %–97.5 % quantiles across the ensemble were extracted to define the central age estimate and its 95 % uncertainty range at each depth. To derive the final chronology, all models that returned successful chronologies were included. At each depth, the composite chronology was calculated as the arithmetic mean of the selected models (equal weights):

$${\mathit{\mu }}_{\mathrm{final}}\left(d\right)={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{M}_i\left(d\right)$$

where M_i(d) is the age of model i at depth d, and N is the number of selected models. The corresponding 95 % uncertainty was derived by combining (i) the between-model spread (standard error of model ages) and (ii) the mean within-model uncertainty (average 95 % half-widths) in quadrature:

$${\mathit{\sigma }}_{\mathrm{final},\mathrm{95}}=\sqrt{{\mathit{\sigma }}_{\mathrm{between}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{{\mathit{\sigma }}_{\mathrm{within}}}}^{\mathrm{2}}}$$

Where geochemically defined annual cycles from $\mathrm{Sr}/\mathrm{Ca}$ variations were identified, the resulting layer-count chronologies were treated as an independent constraint rather than being incorporated into the composite chronology. These were converted to absolute ages by anchoring the layer-count chronology either (i) to the year of stalagmite extraction, if continuous growth until collection can be assumed, or (ii) to the U-Th age with the smallest analytical uncertainty within the counted interval. For intervals in which seasonal variations were less pronounced or individual cycles were indistinct, a conservative uncertainty of ±0.5 layers per 10 counted layers was assigned, reflecting the potential for missed or ambiguous cycles.

3 Results

3.1 Speleothem petrography

One XRD analysis of the bottom portion of stalagmite E1 verified the presence of aragonite, as initially suggested by its fan-like crystal texture and silky gloss (Fig. S2b). In contrast, all nine XRD measurements from different sections of stalagmite E0-C indicated the presence of pure calcite (Fig. S2a). Thin section analyses of stalagmites E0-B, E0-C, E1, E4 and E8 revealed only calcite phases. These primarily consisted of two fabrics: compact- and open-elongated calcite, with crystals exhibiting a length-to-width ratio greater than 6:1 (Frisia, 2015) (Fig. S3). Macroscopically, the compact elongated calcite phase appeared more translucent, while the open elongated fabric displayed a whitish appearance (Fig. S1). Growth interruptions (hiatuses) were macroscopically visible either as major flooding layers (Fig. S1) or as smaller detrital layers, typically near the tops of the stalagmites but occasionally occurring throughout their growth (Fig. S3). Under the microscope, these interruptions were characterized by clear rhombohedral crystal terminations and the presence of detrital minerals such as clay (Fig. S3). No evidence of aragonite-to-calcite recrystallization, such as preserved aragonite textures or “ghosts” of aragonite (Frisia et al., 2002; Perrin et al., 2014), was observed in any of the thin sections. Therefore, it is assumed that stalagmites E0-B, E0-C, E4, E8, and the upper portion of stalagmite E1 consist of primary calcite, predominantly in the form of compact or slightly open elongated fabrics.

3.2 Radiocarbon dating

A total of four ¹⁴C measurements were performed on three different stalagmite tips (E0-C, E1, and E8), with one duplicate measurement taken from the E0-C sample (Table S1). The measured ¹⁴C activities (a¹⁴C_meas) ranged from 97.54±0.20 to 94.95±0.19 pmC (percent modern Carbon). Assuming even only a minor contribution of “dead carbon” (<10 %) to the carbonate in the stalagmite would result in a ¹⁴C_meas values >100 pmC, which suggests modern growth (post-1950 AD) for the stalagmite tips.

3.3 ²³⁰$\mathrm{Th}/\mathrm{U}$ dating

²³⁰$\mathrm{Th}/\mathrm{U}$ dating was performed on 157 samples from the Áaktun Kóopo' Cave speleothems (Table S3). The uranium content of the analysed speleothems spans from 106 to 2237 ng g⁻¹. Measured δ²³⁴U values time-independently average 14±8 ‰, which is close to secular radioactive equilibrium with minimal variation, mirroring the low activity ratios of the drip (Table S3) as well as local ground water (Schorndorf et al., 2023). ²³²Th concentrations are generally low with values <0.1 ng g⁻¹, resulting in (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios >300, indicating minor residual initial Th influence. Overall, ²³⁰$\mathrm{Th}/\mathrm{U}$ dating suggests late Holocene speleothem growth for most studied stalagmites from Áaktun Kóopo' Cave, since most uncorrected ages cluster during the past ca. 2.7 kyr (Table S3). Notable exceptions are stalagmites E6 or E23-1 where uncorrected ages suggest deposition between 333±13 to 264.1±5.7 kyr and 74.90±0.20 to 16.735±0.061 kyr, respectively.

However, several samples show elevated ²³²Th (up to 26.5 ng g⁻¹), young ages, lower uranium concentrations, or a combination of these factors. Consequently, measured ²³⁰$\mathrm{Th}/\mathrm{U}$ ages do not align stratigraphically, and many of the inversions are not resolved when using a bulk earth detrital (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio of ∼0.8, typical for the upper continental crust (Taylor and McLennan, 1985). Notably, the modern carbonate sample yielded an elevated (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio of 6.57±0.03, while (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios in three drip water samples from the cave ranged from 1.659±0.046 to 8.9±1.0 (Table S2), suggesting high levels of initial ²³⁰Th in the system. Since using the drip water or modern carbonate (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio for age correction does not yield satisfying results either, we infer that a more detailed inspection of the elevated initial ²³⁰Th systematics is mandatory to obtain reliable chronologies.

3.4 $\mathrm{Sr}/\mathrm{Ca}$ ratios

For stalagmites E1, E8, and E23-3, LA-ICP-MS analyses revealed pronounced oscillations in $\mathrm{Sr}/\mathrm{Ca}$ ratios from 0.1 to 0.66 mg g⁻¹, permitting identification of strong geochemical cycles (Fig. S4). These cycles serve as a potential proxy for seasonal abundance changes of Sr relative to Ca and thus annual geochemical laminae (Warken et al., 2018; Huang et al., 2001). The interval between two successive $\mathrm{Sr}/\mathrm{Ca}$ minima is here interpreted as representing one year. Annual cycles of strongly variable amplitude and distance were identified by plotting $\mathrm{Sr}/\mathrm{Ca}$ ratios against the dft and manually counting the distance between successive $\mathrm{Sr}/\mathrm{Ca}$ minima. The minima appeared as the best choice since those are sharper and more distinct than the broader maxima. This approach provides an independent means to refine age-depth models derived from ²³⁰$\mathrm{Th}/\mathrm{U}$ dating (Smith et al., 2009) and may ultimately allow for annual relative age precision of climate proxy data.

3.5 Stable isotope values

The stable isotopes of oxygen and carbon analyzed in stalagmites E0-C (Table S4), E1 (Table S5), and E23-3 (Table S6) yielded δ¹³C values ranging between ca. −7.9 ‰ and ca. −12.4 ‰, while the range of δ¹⁸O values is between ca. −7.1 ‰ and ca. −3.8 ‰. Notably, the lowest and highest δ¹³C and δ¹⁸O values in all three stalagmites are very similar, respectively.

4 Discussion

4.1 Determination of initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) ratios

²³⁰$\mathrm{Th}/\mathrm{U}$ dating of Áaktun Kóopo' Cave speleothems reveals a vast spread, including significant age inversions. Petrographic observations suggest no obvious recrystallization or diagenetic alteration; the stalagmites show well-preserved primary fabrics and display largely coherent U concentrations and δ²³⁴U values along their growth axes. Therefore, we consider wide-spread post-depositional U remobilization to play a subordinate role for our dated material. There is, however, strong evidence for elevated initial ²³⁰Th in the Áaktun Kóopo' Cave system. Most age inversions are associated with elevated ²³²Th concentrations and a relationship between (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$) and (${}^{\mathrm{232}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$) activity ratios, consistent with significant Th transport through the karst water, either associated to particles, colloids or chemical complexing agents (Ludwig, 2003). The resulting relationship follows a first order binary mixing, when activity ratios are plotted as Osmond type II isochrons (Fig. 2a). Modern carbonate and drip water samples yielded variable R02i of up to 8.9±1.0, but even with these elevated ratios applied to correct for detrital Th contamination, many speleothem age inversions persist.

Figure 2 Activity ratios from various stalagmites in Áaktun Kóopo' Cave. (a) Relationship between (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$) and (${}^{\mathrm{232}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$) activity ratios. The observed linear trend indicates significant detrital contamination. Different symbols and colors correspond to individual stalagmites. Both axes are plotted on a logarithmic scale. (b) (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios. Samples with (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios <200 (n=68, red squares) were excluded from age-depth models due to significant initial Th contamination. Remaining samples (n=89, black squares) were included in age-depth models and further analysis. No samples with (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) ratios >200 were observed within 7.5 mm of the stalagmite tops, due to the young age of these carbonates and limited ²³⁰Th ingrowth.

Consequently, higher R02i are necessary in cases to meet the stratigraphic constraint of increasing ages from top to bottom (Hellstrom, 2006). To estimate the range of R02i, different approaches are tested, including (i) isochrons and (ii) stratigraphic constraint. In addition, (iii) geochemical $\mathrm{Sr}/\mathrm{Ca}$ cycles were used as annual layer counts to estimate relative age differences between $\mathrm{Th}/\mathrm{U}$ ages and to independently test the other approaches (Domínguez-Villar et al., 2012; Nagra et al., 2017). Given the difficulty in accurately constraining R02i, all samples with (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) ratios <200 (n=68) were excluded from the age–depth modelling and from further chronological interpretation (Fig. 2b), because for such highly Th-contaminated samples the propagated age uncertainty is dominated by the initial Th correction and does not add meaningful information to the chronologies. This criterion also applies to all samples within 7.5 mm of the stalagmite tops, which are generally highly contaminated and have limited ²³⁰Th ingrowth.

4.1.1 Comparison of different approaches

Dense sampling in some stalagmites allows for first-order linear two-endmember mixing regressions (Osmond type 1 – like “isochrons”). In this study, such “isochrons” are applied only to closely spaced, co-genetic subsamples and are interpreted as local, first-order estimates of R02i, similar to the approaches of Warken et al. (2020) or Stinnesbeck et al. (2020). The “isochron” approach using ISOPLOT (Ludwig and Titterington, 1994; Ludwig, 2008; Vermeesch, 2018), produced five successful “isochrons” from three stalagmites (E0-C, E4, E8), which yielded R02i ranging from 6.4±2.4 to as large as 19±10 (Fig. S5, Sect. S1.1). To apply stratigraphic constraints (e.g., Roy-Barman and Pons-Branchu, 2016; Hellstrom, 2006), we used samples with no detectable ²³²Th as stratigraphic anchors, as their ages are unaffected by initial Th contamination and because those provide a steady growth pattern. In addition, we adopt the assumption that stalagmite tops with elevated ²³²Th are younger than uncorrected ages below them, but older than the year of collection. From Stalagmite E0-C, 34 of 51 samples could be placed in stratigraphic order using individual R02i ranging from 11.5 to 49 (Fig. 4b). In contrast, for other stalagmites, such as E0-B or E1, only single ages can be corrected using a stratigraphically derived ratio. Lastly, we exploit distinct cyclic variations in $\mathrm{Sr}/\mathrm{Ca}$ ratios from LA-ICP-MS measurements (Fig. S4) to estimate relative age differences between ²³⁰$\mathrm{Th}/\mathrm{U}$ measurements, thus allowing to correct the ages and estimate R02i (James et al., 2025). For example, in stalagmite E23-3, in total 331±17 $\mathrm{Sr}/\mathrm{Ca}$ minima are counted from 0 to 84 mm dft. Assuming the topmost layer to correspond to 2022 AD (the year of collection), the R02i is estimated to 26 ± 7 to align the ²³⁰$\mathrm{Th}/\mathrm{U}$ ages with the layer-counted age model presuming $\mathrm{Sr}/\mathrm{Ca}$ seasonal cycles within uncertainties of both methods. In contrast, in the upper section of stalagmite E8 (8.5–26 mm dft), ages could only be corrected with individual R02i ranging from 7 to 35.5 to match the $\mathrm{Sr}/\mathrm{Ca}$-based age-depth model. Overall, to correct stalagmite ages with the different approaches, R02i ranging from 4 to 68 are found (Table S3). Notably, only one sample from stalagmite E1 requires a ratio as high as 68, while most other samples need ratios between 30 and 50 or lower. R02i < 20 are only found for three ages via stratigraphic constraint and for nine ages by isochrons. Therefore, samples without specific estimates are corrected by using the mean ratio of 35.5 ± 32.5 to account for the full range of possible initial Th corrections (Table S3). Taken together, these independent approaches indicate that strongly elevated initial Th activity ratios are plausible and necessary at least in parts of the record to obtain stratigraphically consistent age models, even though the exact numerical values of R02i for individual, highly contaminated samples remain poorly constrained.

4.1.2 Elevated and variable initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) ratios in Áaktun Kóopo' Cave

Overall, our analysis shows that significantly elevated and strongly variable R02i (up to 68) are required to correct stalagmite ages from Áaktun Kóopo' Cave for initial Th contamination (Table S3). These values greatly exceed the commonly assumed upper continental crust ratio of (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) ∼0.8 (Taylor and McLennan, 1985). However, they align with findings from other tropical sites (e.g., the Caribbean and Southeast Asia) where similarly high ratios have been reported (e.g., Beck et al., 2001; Richards and Dorale, 2003; Partin et al., 2007; Carolin et al., 2013; Arienzo et al., 2015; Hoffmann et al., 2010; Moseley et al., 2015; Ridley et al., 2015; Stinnesbeck et al., 2020; Vieten et al., 2024a, 2024b; Steidle et al., 2021; Warken et al., 2021; James et al., 2025). Furthermore, our measurements demonstrate that this pronounced variability is largely unsystematic, occurring not only spatially across different stalagmites within the same cave, similar to findings in Larga Cave, Puerto Rico (Vieten et al., 2024b; Warken et al., 2020; Kerber et al., 2025), but also temporally within individual stalagmites.

The vast spread of R02i may be attributed to the varying contributions of multiple Th reservoirs within the karst and epikarst environments (Li et al., 2022; Wortham et al., 2022; Huang et al., 2024). Rather than a single source, Th can be transported as a complex mixture of detrital particles (e.g., clay minerals, aluminosilicates, iron oxyhydroxides), colloidal phases, and organic complexes (Morton et al., 2001; Nascimento et al., 2019; Reiller et al., 2002; Richards and Dorale, 2003). In tropical environments with relatively young reef limestones as host rock, such as the YP, the rapid dissolution and varying mobilization of these different organic and inorganic fractions lead to highly variable initial Th isotopic signatures. This dynamic reflects the observations from various aquatic systems, where the (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio has been shown to vary by orders of magnitude (Scott, 1968; Hubert et al., 2006; Hirose et al., 2012; Töchterle et al., 2022; Moore and Sackett, 1964; Moore, 1981). Consequently, highly elevated and time-variable R02i should be regarded as the rule rather than the exception in such settings, fundamentally questioning the applicability of a static “bulk earth” value for age corrections. In such settings, initial Th corrections are best treated as approximations with explicit uncertainty, constrained where possible by multiple, independent estimates of R02i.

Because we do not find any systematic behaviour that would allow to predict R02i a priori, developing accurate age models remains challenging. Initial Th corrections are particularly critical for younger speleothems (less than a few thousand years old), as even small uncertainties in the R02i can significantly influence age models. While monitoring ZAC and drip water provides a baseline for estimating R02i (Wortham et al., 2022; Li et al., 2022), this approach often fails to capture long-term temporal variability. In Áaktun Kóopo' Cave, for example, the (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio from supposedly modern carbonate and drip-water samples remains far below the highest estimated ratios (up to 68) required to align stalagmite ages with the sedimentary sequence. Therefore, to reduce the uncertainty of corrected ages in similar karst settings, it is strongly recommended to employ a combination of independent validation methods. Integrating ²³⁰$\mathrm{Th}/\mathrm{U}$ dating with radiocarbon or other radiometric dating methods (Akers et al., 2019; Huang et al., 2024; Sengupta et al., 2025) or annual layer counting (e.g., Warken et al., 2025; Domínguez-Villar et al., 2012; James et al., 2025; Faraji et al., 2021; Nagra et al., 2017) provides a more robust framework to constrain these erratic detrital inputs and secure reliable chronologies.

4.2 The Áaktun Kóopo' chronology

The finally corrected ages of speleothems from Áaktun Kóopo' Cave range from as old as 353 ± 7 kyr to as young as 0.13 ± 0.015 kyr (Table S3). Five older samples indicate speleothem growth during earlier glacial and interglacial periods, but most samples (84 out of 89) are younger than 3 kyr (Fig. 3). Even though most ²³⁰$\mathrm{Th}/\mathrm{U}$ samples with low (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios <200 were rejected, radiocarbon measurements on speleothem tops (E1, E0-C, E8) as well as seasonal $\mathrm{Sr}/\mathrm{Ca}$ chronologies (E1, E8, E23-3) suggest carbonate deposition within the last 130 years, likely continuing up to the present (Sect. S1.1). Overall, the stalagmites from Áaktun Kóopo' Cave provide a robust basis for paleoclimatic and paleoenvironmental reconstructions over the late Holocene.

Figure 3 (a) Distribution of speleothem samples according to their corrected ages on a logarithmic scale. Most samples are younger than 2.7 kyr, with some carbonate deposition at ∼75, 260, and 350 kyr. (b) ²³⁰$\mathrm{Th}/\mathrm{U}$ ages plotted against the LR04 benthic δ¹⁸O stack (Lisiecki and Raymo, 2005), highlighting deposition during earlier glacial and interglacial periods, with marine isotope stages shown in blue and red.

Figure 4 Age-depth models for Áaktun Kóopo' Cave stalagmites E0-B (a), E0-C (b), E1 (c), E4 (d), E8 (e) and E23-3 (f), generated using Bacon, Bchron, CopRa, linInt, linReg, and StalAge, respectively. The final composite age-depth model is shown in black with its corresponding 95 % uncertainty band in grey, representing the arithmetic mean of the six model chronologies (equal weights). The ²³⁰$\mathrm{Th}/\mathrm{U}$ ages used in the models are indicated along with identified growth hiatuses and major flooding layers (E0-C, E1, E4, E8, E23-3). A scan of each stalagmite is shown on the right of each plot. Speleothem scans are at the same vertical scale as the y-axis (dft) values.

4.2.1 Late Holocene speleothem growth history

Figure 4 shows the individual age-depth models constructed for stalagmites E0-B, E0-C, E1, E4, E8, and E23-3, as described in detail in Sects. 2 and S1.1. Although some of the analysed stalagmites exhibit decade-long growth interruptions, their chronologies overlap and suggest continuous carbonate deposition in Áaktun Kóopo' Cave over the last 2.7 kyr with mean growth rates of 200–300 µm a⁻¹ (Fig. 5).

Figure 5 Mean age-depth models for stalagmites E0-B, E0-C, E1, E4, E8, and E23-3 from Áaktun Kóopo' Cave. Numbers indicate the average growth rates (in µm a⁻¹) during different growth periods. The coloured horizontal bars at the bottom illustrate the overlapping growth periods of the stalagmites over the last 2.7 kyr, highlighting intervals of simultaneous growth across multiple stalagmites. Lighter shades within the bars reflect the associated age uncertainties. Grey bars indicate periods covered by other published stalagmite records from the northern Yucatán Peninsula (YP). Notably, while most records are from the northwestern YP (Kennett et al., 2022; Medina-Elizalde et al., 2010, 2016a; James et al., 2025, 2026; Frappier et al., 2014), there is presently a marked absence of stalagmite proxy records from the northeastern YP for the last 1600 years, with only one record covering an earlier period between 1037 BC and 397 AD (Medina-Elizalde et al., 2016b).

While E0-B covers the oldest growth period (633 ± 161 BC to 403 ± 54 AD, Fig. 4a), the growth of its “sibling” E0-C can be divided into three periods separated by hiatuses: (1) 538 ± 56 to 811 ± 52 AD; (2) 994 ± 29 to 1311 ± 15 AD; and (3) 1451 ± 22 to 1835 ± 53 AD (Fig. 4b). The last millennium is also covered by several speleothems (Fig. 5), with E1 growing from 1440 ± 39 to 1577 ± 42 AD and between 1635 ± 24 to 1926 ± 20 AD (Fig. 4c). Stalagmite E4 grew from approximately 228 ± 88 to 1611 ± 118 AD, even though high initial Th contamination led to the rejection of 24 of 28 samples (Fig. 4d). Similar challenges arose for the chronology of E8, where growth could be tightly constrained for the latest part (1605 ± 24 to 1733 ± 33 AD, Fig. 4e). Lastly, the growth history of stalagmite E23-3 can be divided into two periods, with the older part covering 1070 ± 7 to 1347 ± 9 AD, and the upper part spanning 1441 ± 5 AD to the year of collection in 2022 AD (Fig. 4f).

Figure 5 compares the late Holocene Áaktun Kóopo' Cave growth history with periods covered by other published stalagmite records from the northern YP. Notably, while most records are from the northwestern or central YP (Kennett et al., 2022; Medina-Elizalde et al., 2010; James et al., 2025, 2026; Frappier et al., 2014; Medina-Elizalde et al., 2016a), there is a marked absence of stalagmite proxy records from the northeastern YP for the last 1600 years, with only one record covering an earlier period between 1037 BC and 397 AD (Medina-Elizalde et al., 2016b). Given the known regional heterogeneity in past hydroclimate reconstructions and projected trends, the Áaktun Kóopo' Cave record thus fills an important gap in the region (Douglas et al., 2016; Obrist‐Farner et al., 2023; Steinman et al., 2022).

4.2.2 Origin of flooding layers and other growth interruptions

While precise ²³⁰$\mathrm{Th}/\mathrm{U}$ dating establishes the absolute chronological framework for these stalagmites, the numerous growth interruptions (hiatuses) observed throughout their growth histories provide an additional, independent means to cross-validate these highly corrected age models. Detectable either as macroscopic millimeter- to centimeter-thick mud layers or via microscopic thin-section analysis, these hiatuses indicate discrete environmental events that simultaneously affected multiple stalagmites within the cave (Fig. 4).

One prominent flooding layer (∼1 cm thick) appears in stalagmites E0-C, E1, and E23-3, which were all collected from the same cave area at a similar elevation. This layer likely represents a single event around 1400 AD. By comparing the ²³⁰$\mathrm{Th}/\mathrm{U}$ dates across the stalagmites, the timing can be constrained between ∼1350 and ∼1430 AD. The latest confirmed carbonate precipitation prior to the event is recorded by stalagmite E23-3, which yielded a ²³⁰$\mathrm{Th}/\mathrm{U}$ age of 1347 ± 9 AD. Following the flood, stalagmite E1's growth resumed directly upon the detrital substrate at 1428 ± 40 AD. The sediment sequence of Sullivan et al. (2025) from Cenote Muyil ca. 80 km south of Áaktun Kóopo' Cave recorded a period of intense hurricane activity between 1285 and 1420 AD in the area, which suggests that a major storm event may have caused the flooding of the cave during that time. Indeed, historical records report a major hurricane (“hurricane of the four winds”) that impacted the area in 1464 AD (Sullivan et al., 2025; Masson and Lope, 2014), which is however not fully consistent with our date, suggesting another extreme event caused the mud layer in Áaktun Kóopo' Cave. Thinner layers (0.5–1 mm) in stalagmites E1 and E8, which also stem from the same area in the cave, correspond to a shorter interruption that may have lasted only months to years. The timing can be precisely dated to 1624 ± 24 AD within the continuous age models of stalagmite E8, perfectly aligning with a brief hiatus observed in E1 between ∼ 1590 and ∼ 1630 AD. In addition to detecting and dating single extreme events that may have impacted the environment and the people living in the area, such markers could serve as key stratigraphic tie-points, analogous to ash layers in marine sediments or ice cores, for anchoring proxy records and refining chronologies (Warken et al., 2025). Fine laminations (alternating carbonate and mud/clay) of stalagmite E4 even suggest recurrent flooding in the part of the cave near the entrance. While the current efforts limit a more precise dating of the individual layers, future studies could deliver a more detailed history of local flooding (Frappier et al., 2014; Denniston and Luetscher, 2017; Gonzalez-Lemos et al., 2015).

Although many hiatuses coincide with clear flooding layers, other interruptions could reflect periods of reduced effective moisture due to reduced precipitation and/or enhanced evapotranspiration, or other processes. For example, the special morphology of Áaktun Kóopo' Cave (Fig. 1b) makes it conceivable that a sufficiently high-water table sealed cave ventilation in some sections of the cave, leading to elevated cave CO₂. Because Áaktun Kóopo' Cave was visited and modified by the Maya, at least during the Late Preclassic (−250 to 300 AD) and Early Classic (300–600 AD) periods (Gómez, 2020), human influence on speleothem development can also not be ruled out. Deforestation and agricultural practices above the cave, for example, may have altered soil cover and infiltration pathways, thereby influencing recharge and ultimately affected speleothem formation. Likewise, cave use and ritual activities might have modified the cave's microclimate (e.g., CO₂ concentration, humidity, temperature) or sedimentological conditions in ways that impacted speleothem formation, although such influences remain speculative. Disentangling natural hydrological variability from potential anthropogenic influences on speleothem formation remains a challenging task and requires a comprehensive multi-proxy approach that explores all aspects of speleothem formation, aimed at separating climatic from human-induced signals.

4.3 Age model uncertainties and algorithm comparison

4.3.1 Evaluation of the chronologies

Evaluating the diverse performance of the applied age-depth modelling algorithms reinforces the necessity of our ensemble approach. Consistent with Comas-Bru et al. (2020), no single algorithm universally outperformed the others across all speleothems, as the underlying statistical assumptions of each model respond differently to variations in growth rates and dating resolution. This methodological variability becomes especially evident when comparing the models against independent layer-counting chronologies: while the Bacon model closely matched the visual layer counts for stalagmite E23-3, Bchron and linear interpolation provided the tightest fits for different sections of E8. Furthermore, the individual algorithms diverged significantly in their chronological treatment of extreme events, such as the major flooding layers and resulting hiatuses. The chosen algorithm significantly affects not only the estimated duration of an event, but also its absolute timing. For instance, the first hiatus in stalagmite E0-C (at 40 mm depth) has an estimated duration ranging from 96 ($+\mathrm{58}/-\mathrm{55}$) years (linear regression) to 174 ($+\mathrm{85}/-\mathrm{88}$) years (Bacon). Consequently, the absolute timing of the growth interruption shifts across models: Bchron dates the onset of the hiatus to 1302 ($+\mathrm{25}/-\mathrm{4}$) AD, whereas linear regression places it almost four decades later at 1341 ($+\mathrm{26}/-\mathrm{25}$) AD. Constructing a composite “mean” age model is therefore justified, as relying on any single algorithm risks either artificially over-smoothing or over-estimating the timing and duration of discrete environmental events.

Figure 6 Stable oxygen (δ¹⁸O) and carbon (δ¹³C) records of speleothems E1, E23-3 and E0-C for the period from 1400 to 2022 AD binned to three-year averages, respectively. In panel (a), the records are plotted on their originalage scale. In addition, the respective ²³⁰$\mathrm{Th}/\mathrm{U}$ ages are shown to indicate the uncertainty in the respective chronologies. Panel (b) shows the δ¹⁸O and δ¹³C records on the composite timescales, demonstrating a remarkable agreement. Note that the overlap between E23-3 and E0-C is too short to calculate meaningful correlation coefficients.

4.3.2 Validation of age models via stable isotope comparison

To further assess the robustness of our age models, we compared the independent stable oxygen (δ¹⁸O) and carbon (δ¹³C) isotope records of stalagmites E0-C, E1, and E23-3 for their overlapping growth periods. The chronology of stalagmite E23-3 is considered the most reliable anchor for the last 400 years. It is constrained by continuous counting of seasonal $\mathrm{Sr}/\mathrm{Ca}$ cycles and anchored to the year of extraction (2022 AD), resulting in a cumulative counting error of ±17 years at 1691 AD. When plotting the stable isotope time series using the individually established ²³⁰$\mathrm{Th}/\mathrm{U}$ age-depth models, the records already show an apparent coherence in both δ¹⁸O and δ¹³C values over the past 600 years (Fig. 6a). Despite being fed by different drip sites within the same cave chamber, the stalagmites capture nearly identical multi-decadal to centennial-scale geochemical patterns. This coherence suggests that the individual proxy records are driven by a common external forcing, such as regional hydroclimate variability affecting cave recharge and prior calcite precipitation (PCP), rather than localized kinetic or drip-specific effects (Skiba and Fohlmeister, 2023; Baker et al., 2019; Treble et al., 2022). However, the individual chronologies appear to be shifted relative to each other, reflected in the weak or insignificant correlations between the stable isotope records. For example, the largest significant coefficient of 0.2 (p<0.01) occurs for the δ¹⁸O records of E23-3 and E1. To precisely align the time series, a wiggle-matching procedure was performed using QAnalySeries software (Kotov and Pälike, 2018; Pälike, 2024), similar to the approach of Huang et al. (2024). The upper section of the E1 δ¹⁸O record (3–61 mm dft) was tuned to the E23-3 reference record using 10 tie points, while the upper section of E0-C (5–38.6 mm dft) was subsequently matched to E1 using 11 tie points. As a result, the maximum shift of E1 with respect to E23-3 is ca. 20 years directly after the thin flooding layer around 1620 AD, while the other sections are only corrected by <7 years. The chronology of E0-C requires larger shifts by around 25–70 years with respect to E23-3. Overall, the resulting adjustments required to align these independent records fall entirely within the calculated 95 % uncertainty ranges of their respective ²³⁰$\mathrm{Th}/\mathrm{U}$ age-depth models (Fig. 6b). After the synchronization, the δ¹⁸O (δ¹³C) record of E1 exhibits a correlation of 0.62 (0.61) with E23-3 and 0.42 (0.51) with E0-C (p<0.01, respectively), demonstrating a striking coherence between the stalagmite records. The successful replication of stable isotope patterns across three distinct stalagmites and the tuning required to synchronize them strongly validate our geochronological approach. It demonstrates that the strategies used to correct for elevated and highly variable initial ²³⁰Th contamination, specifically the integration of ²³⁰$\mathrm{Th}/\mathrm{U}$ dating with seasonal trace-element layer counting and stratigraphic constraints, yield robust and accurate chronologies suitable for high-resolution paleoclimate reconstruction on the YP.

5 Conclusions

Reconstructing reliable paleoclimate records from speleothems requires precise chronologies, which are frequently compromised by variations of initial Th contamination through time in tropical karst environments. While low-Th speleothems with simple initial ratios would be preferable, such material is not always available at sites of high archaeological and paleoclimate interest, so working with variable Th-contaminated samples is sometimes unavoidable if any chronological information is to be obtained.

In this study, we generated robust age-depth models for multiple late Holocene stalagmites from Áaktun Kóopo' Cave, Yucatán, spanning the past 2.7 kyr. ²³⁰$\mathrm{Th}/\mathrm{U}$ dating revealed significant and highly variable initial Th contamination, with initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios ranging unpredictably between 4 and 68. These values vastly exceed standard bulk earth corrections and demonstrate substantial spatial and temporal variability within the cave system. This indicates that static initial Th corrections are insufficient in such settings, and modern drip-water or ZAC monitoring may not capture the full range of initial ²³⁰Th variability. Considering the vast number of other studies suggesting high and variable R02i ratios in (sub-)tropical caves, we expect that this could be also a characteristic of other, comparable cave systems. To determine R02i and resolve chronological inversions, we successfully employed a multi-method approach, combining isochron methods, stratigraphic constraints, and relative seasonal layer counting ($\mathrm{Sr}/\mathrm{Ca}$). We note, however, that R02i constrained by local “isochrons” assuming binary-mixing relationships should be regarded as helpful local, first-order approximations to delimit the plausible range of initial Th activity ratios, rather than as exact descriptions of the underlying Th sources.

Furthermore, our analysis highlights that different age-depth algorithms (e.g., Bacon, Bchron, copRa, StalAge, Linear Regression) may yield different results, particularly when handling growth discontinuities. This can affect both the estimated duration and the absolute timing of hiatus or other events; thus, the integration of multiple constraints can yield more reliable chronologies.

Ultimately, the improved chronologies and accompanying coherent proxy records from Áaktun Kóopo' Cave now provide a robust foundational framework for future high-resolution reconstructions of hydroclimate and environmental change throughout the timeline of Maya cultural evolution on the YP, e.g., how long-term trends and short-term extreme events, such as the identified flood layers, are related and how they may have intersected with phases of cultural change.