Speleothem evidence for Late Miocene extreme Arctic amplification – an analogue for near-future anthropogenic climate change?

The Miocene provides an excellent climatic analogue for near-future runaway anthropogenic warming, with atmospheric CO₂ concentrations and global average temperatures similar to those projected for the coming century under extreme-emissions scenarios. However, the magnitude of Miocene Arctic warming remains unclear due to the scarcity of reliable proxy data. Here we use stable oxygen isotope and trace element analyses, alongside clumped isotope and fluid inclusion palaeothermometry of speleothems to reconstruct palaeo-environmental conditions near the Siberian Arctic coast during the Tortonian (8.68 ± 0.09 Ma). Stable oxygen isotope records suggest warmer-than-present temperatures. This is supported by temperature estimates based on clumped isotopes and fluid inclusions giving mean annual air temperatures between +6.6 and +11.1 °C, compared with −12.3 °C today. Trace elements records reveal a highly seasonal hydrological environment.

Our estimate of > 18 °C of Arctic warming supports the wider consensus of a warmer-than-present Miocene and provides a rare palaeo-analogue for future Arctic amplification under high-emissions scenarios. The reconstructed increase in mean surface temperature far exceeds temperatures projected in fully coupled global climate models, even under extreme-emissions scenarios. Given that climate models have consistently underestimated the extent of recent Arctic amplification, our proxy data suggest Arctic warming may exceed current projections.

1 Introduction

The Arctic and sub-Arctic regions have warmed nearly 4 times faster than the global average since 1979 (Rantanen et al., 2022), and this disproportionate warming is expected to continue over the coming decades (Ma et al., 2022). Climate models have consistently underestimated the magnitude of this so-called Arctic amplification (Chylek et al., 2022; Rantanen et al., 2022) and show considerable discrepancies in its predicted magnitude (Smith et al., 2020; Taylor et al., 2022). Uncertainty in Arctic temperature projections is greater than for any other region on the planet (Taylor et al., 2022).

Improving Arctic climatic projections is imperative for informing local adaptation efforts. Warming is driving permafrost thaw and ice retreat and impacting ecosystems, infrastructure, and access to food, whilst weather-related changes are increasing natural hazards (Ford et al., 2015). Localized warming could have global impacts. The Arctic is home to most of the planet's permafrost, a possible climate tipping element (McKay et al., 2022; Nitzbon et al., 2024). Permafrost is the Earth's largest terrestrial carbon pool (Strauss et al., 2024), and its future thaw will play a major role in potential climate trajectories. Furthermore, modelling studies suggest Greenland Ice Sheet mass loss, driven by Arctic warming, will be a major contributor to 21st century sea level rise, with the potential to significantly weaken the Atlantic Meridional Overturning Circulation (Hofer et al., 2020).

Detailed palaeoclimate records present an excellent opportunity to verify climate models and improve future projections. Whilst past warm intervals represent imperfect analogues for near-future climate because of differences in climate forcing, these periods can provide important examples of planetary response to atmospheric warming. Recently, Steinthorsdottir et al. (2021b) proposed the Miocene (23.03–5.33 Ma) as a suitable palaeo-analogue for future anthropogenic climate change. The distribution of continental landmasses was similar to the present day, and atmospheric pCO₂ concentrations were close to modern levels. Reconstructions generally agree on atmospheric pCO₂ values between 400 and 600 ppm throughout the Miocene (Steinthorsdottir et al., 2021a). Whilst some estimates suggest this may have reached as high as 800–1000 ppm during the Miocene Climatic Optimum (MCO) (Rae et al., 2021), most studies suggest it was considerably lower, likely between 430 and 630 ppm (Sosdian et al., 2018; Super et al., 2018). Subsequent cooling of ca. 6 °C was coincident with a ca. 125 ppm decline in atmospheric CO₂ (Super et al., 2018). By the Late Miocene (ca. 7–5.4 Ma), multiple reconstructions converge on atmospheric carbon dioxide concentrations between 300 and 500 ppm (Rae et al., 2021; Sosdian et al., 2018; Super et al., 2018), close to modern-day levels (ca. 420 ppm), with the upper estimate similar to temperatures projected for the coming decades at current emission rates (IPCC, 2021).

Burls et al. (2021) synthesized available proxy temperature reconstructions to deduce a global mean surface temperature of 21.21 ± 0.56 °C for the Early Miocene. Temperatures declined rapidly following the MCO, particularly in the high latitudes, where cooling up to 8 °C is observed in marine records (Herbert et al., 2016). By the Late Miocene, global mean sea surface temperature (SST) was ca. 6 °C higher than the present day (Herbert et al., 2016), and mean terrestrial surface temperature was ca. 4.5 °C above the modern day (Pound et al., 2011). Proxy reconstructions suggest a much-reduced latitudinal temperature gradient, which saw Arctic mean annual temperature (MAT) 11–19 °C warmer than present during the Middle Miocene (Steinthorsdottir et al., 2021b, and references therein), declining to between 5 and 6 °C above the pre-industrial level by the Late Miocene (Pound et al., 2012). Whilst this period saw the establishment of a stable Antarctic ice sheet (Miller et al., 2020), biomarker evidence suggests that ice-free summers in the Arctic persisted (Stein et al., 2016).

Modelling reconstructions, driven by inferred Miocene atmospheric CO₂ concentrations, have consistently failed to reconstruct the high temperatures and reduced latitudinal temperature gradient seen in the Miocene proxy records (Goldner et al., 2014). It is unclear whether this discrepancy arises due to missing climatic feedbacks within climate models or consistent under-/over-estimation of atmospheric CO₂ concentrations and temperatures in proxy reconstructions, but resolving it is of importance for estimation of future warming. Here we present a multi-proxy analysis of four speleothems (cave carbonate deposits) from the Siberian Arctic to infer environmental conditions during the Tortonian Stage of the Late Miocene Sub-epoch. We use isotope measurements of speleothem fluid inclusions and clumped isotope analyses to estimate multi-annual mean surface air palaeotemperatures and stable oxygen isotope and trace element data to reconstruct palaeo-hydrology and seasonality, providing a new, high-precision dataset to understand Miocene climate sensitivity.

2 Study site and sample material

2.1 Study site

Samples were collected from the Taba Ba'astakh cliffs (72.27° N, 126.94° E), which rise ca. 140 m above the eastern riverbank near the entrance of the Lena River delta, ca. 100 km northwest of Tiksi, Sakha Republic, Russian Federation (Fig. 1). The cliffs are composed of Carboniferous-aged carbonates (Mikhaltsov et al., 2018). Lower strata (up to ca. 50 m above current river level) comprise fine quartzitic sandstone, siltstone, sandy dolomite, and dolomitic limestone with sparse fossil coral and ostracod and foraminifera assemblages broadly linked to the Tournaisian (358.9–346.7 Ma) (Izokh and Yazikov, 2017). Above, ca. 50 m from the riverbank, bioclastic carbonate, calcareous siliciclastics, and carbonates form a subsiding platform. Continuous permafrost, 400–600 m thick, is found in this region, with an active-layer thaw depth between 0.6 and 0.4 m (Boike et al., 2019). All caves encountered during the expedition contained significant ice deposits, rendering modern speleothem formation impossible. Modern tundra vegetation is sparse, consisting mostly of grasses, mosses, and lichens, with a small number of bushes (mainly polar willow). Trees are not found above the caves today.

Figure 1 (a) Location of Taba Ba'astakh sampling site, Samoylov Research Station, and relevant Miocene palaeoclimate reconstruction sites: Temmirdekh-khaj and Omoloy (Popova et al., 2012). The green shading shows forested areas (map created in Esri ArcGIS). (b) Photograph of the Taba Ba'astakh cliffs (height ca. 120–140 m) with cave locations outlined with yellow dashed lines (adapted from Vaks et al., 2025).

Meteorological data recorded at the Samoylov Island Research Station in the central Lena River Delta (72.37° N, 126.48° E), ca. 20 km northwest of our study area, indicate a polar tundra climate (ET according to the Köppen–Geiger classification; Peel et al., 2007), with a mean annual air temperature (MAAT) of −12.3 °C, July average of +9.5 °C, and February average of −32.7 °C (all temperature data averaged between 1998 and 2017) (Boike et al., 2019). The closest available mean annual rainfall estimate is 309 mm (1980–2018), measured at the Tiksi meteorological station 90 km southeast of Taba Ba'astakh (Buchwal et al., 2023).

2.2 Sample description

Fourteen speleothem samples were collected from the base of the Taba Ba'astakh cliffs, on the bank of the Lena River in 2014. The modern caves are ice filled and inaccessible after a few metres, but erosion of the cliff face has exposed the interior of relic caves, with speleothems observed along the cliff walls. Observations of ongoing weathering of cave walls suggest that these samples originated from the cliff face 70–120 m above the sampling site where some in situ speleothems were also found. All 14 samples have been dated to the Tortonian Stage (Vaks et al., 2025) using U–Pb dating following a modified method described in Vaks et al. (2020) and Mason et al. (2022). The ²³⁵U–²⁰⁷Pb–²⁰⁸Pb/²⁰⁷Pb isochron age of all samples yielded an age of 8.68 ± 0.09 Ma, with a mean square weighted deviation of 1.2 (Vaks et al., 2025). Overlapping uncertainties of top and bottom ages only allow the construction of relative age–depth models along the growth profile of individual samples.

Four samples, STBB I – 1, STBB II – 7, STBB 4 – 2, and STBB 4 – 3 (Fig. 2), were selected for this study, as their high-precision dates and large size make them amenable for coupled reconstruction of temperature and climate conditions at their time of deposition. All samples are composed of calcite (Vaks et al., 2025).

Figure 2 Speleothem samples used in this study. Thick translucent tracks, overlain with dashed lines, show the positions of stable isotope and trace element sampling profiles, respectively (STBB I – 1 and STBB II – 7 only). Boxes indicate clumped isotope sampling positions.

3 Methodology

3.1 Clumped isotope geothermometry

Clumped isotope analyses for estimating speleothem formation temperatures were performed on subsamples of all four speleothems. Measurements were conducted at Northumbria University using a Nu Instruments Perspective dual-inlet isotope ratio mass spectrometer coupled with an automated NuCarb preparation system. Samples were drilled from layers with well-constrained ages and 325 ± 25 µg of powder loaded into sample vials, which were then evacuated and reacted with concentrated orthophosphoric acid at 70 °C. The reactant gas was dehydrated and cleaned using standard procedures (Petersen et al., 2015). Briefly, samples were dehydrated using two liquid nitrogen cryotraps, cooled to −80 °C, and scrubbed of contaminants by passing them through a 1 cm static trap filled with a Porapak™ Q (Waters Corporation) absorbent material, cooled to −30 °C. Each replicate digestion was measured 40 times in dual-inlet mode (total measurement time ca. 40 min). No pressure baseline correction was required.

Final clumped isotope (Δ₄₇) values were calculated using the software Easotope (http://www.easotope.org, last access: 15 June 2023; John and Bowen, 2016) and the D47crunch Python package (Daëron, 2021). Both methodologies used the International Union of Applied Chemistry (IUPAC) parameters for ¹⁷O correction and isotopic ratios for Vienna Pee Dee Belemnite (VPDB) and Vienna Standard Mean Ocean Water (VSMOW) (Bernasconi et al., 2018; Brand et al., 2010; Daëron et al., 2016). Application of both methods allowed Δ₄₇ values to be inputted into both our in-house calibration, derived using Easotope, and the composite calibration of Anderson et al. (2021), which is derived using D47crunch.

For the Easotope method, sample Δ₄₇ values were projected onto carbon dioxide equilibrium space (I-CDES-90) following the methodology of Dennis et al. (2011), using standards ETH1, ETH2, and ETH3, and I-CDES Intercarb Δ₄₇ values from Bernasconi et al. (2021). Long-term instrument performance was monitored with an internal standard, Pol-2 (a natural cave pool rim deposit), giving a long-term Δ₄₇ external standard deviation of 0.032 ‰ (Δ₄₇ = 0.656 ± 0.003 ‰ I-CDES-90, $T_{{\mathrm{\Delta }}_{\mathrm{47}}}$ = 2.4 ± 1.0 °C, N = 65).

Replicates were measured across multiple runs over a long period of time, similar to the sliding-window approach employed by, for example, Meinicke et al. (2020) and described by Daëron and Gray (2023). Stable isotope and Δ₄₇ outliers (> $\stackrel{\mathrm{‾}}{x}$ ± 2σ), alongside samples with elevated Δ₄₈, indicative of contamination (Eiler and Schauble, 2004), were discarded from final Δ₄₇ calculations (Table S4 in the Supplement). After pruning, a minimum of 14 replicate analyses were made for each sample such that Δ₄₇ 95 % confidence intervals of ≤ 0.015 ‰ were achieved (i.e. ≤ 0.007 ‰ standard error). All Δ₄₇ uncertainties are quoted as standard errors and 95 % confidence intervals, according to best practices outlined by Fernandez et al. (2017).

For the D47crunch method, contaminated samples were pruned (identified by elevated Δ₄₈ and Δ₄₉), alongside outlier ETH3 replicates (> $\stackrel{\mathrm{‾}}{x}$ ± 2σ for δ₄₅–δ₄₇) from each run. A linear drift correction was calculated for δ₄₅–δ₄₇ using ETH3, which was then applied to all sample and standard data. Outlier measurements of ETH1, ETH2, and ETH3 (> $\stackrel{\mathrm{‾}}{x}$ ± 2σ for δ₄₅–δ₄₇) were then pruned from each run before being input into D47crunch to calculate I-CDES-90 Δ₄₇ values. Stable isotope and Δ₄₇ outliers (> $\stackrel{\mathrm{‾}}{x}$ ± 2σ) were pruned from the replicate data before final Δ₄₇ calculations.

Clumped isotope temperatures ($T_{{\mathrm{\Delta }}_{\mathrm{47}}}$) were calculated by inputting D47crunch-derived Δ₄₇ values into the composite regression of Anderson et al. (2021) and Easotope-derived Δ₄₇ values into an internal laboratory calibration calculated from 17 natural inorganic calcite samples precipitated at known temperatures (Fig. S1 in the Supplement, Tables S1 and S2). The derived regression equation, calculated using the York least-squares method (York et al., 2004) in the R geostats package, is

$$\begin{array}{}\text{(1)}& {\mathrm{\Delta }}_{\mathrm{47}}=\mathrm{0.0372}(\pm \mathrm{0.0008})\cdot {\displaystyle \frac{{\mathrm{10}}^{\mathrm{6}}}{T^{\mathrm{2}}}}+\mathrm{0.166}(\pm \mathrm{0.01}).\end{array}$$

Bracketed numbers denote 1 standard error. This in-house regression has a slightly lower gradient and higher intercept than the composite regression derived by Anderson et al. (2021) (gradient = 0.0391 ± 0.0004, intercept = 0.154 ± 0.004). Measurement and calibration uncertainties were propagated together to calculate temperature uncertainty following the methods of Huntington et al. (2009).

3.2 Fluid inclusion analysis for temperature estimation

Fluid inclusion analyses were conducted on STBB I – 1 and STBB II – 7, where microscopic inspection confirmed the presence of fluid inclusions. Measurements were performed at the University of Innsbruck, Austria, following the method of Dublyansky and Spötl (2009). An aliquot of 1.5 g of carbonate was crushed in a heated crusher under He flow. The released fluid inclusion water was cryo-trapped and then admitted into the TC/EA analyser (Thermo Scientific). After pyrolysis on contact with glassy carbon at 1400 °C, the evolved H₂ and CO were separated in the GC column and admitted to a Delta V Advantage mass spectrometer (Thermo Scientific). Calibration was performed by measuring various amounts of reference waters with isotopic compositions bracketing the expected compositions of the sample. Fluid inclusion isotope data are reported with respect to VSMOW. Accuracy of the measurement is 1.5 ‰ for δ²H and 0.2 ‰ for δ¹⁸O.

Temperatures from fluid inclusions were calculated using the calibrations of Tremaine et al. (2011) and Coplen (2007) for water–calcite oxygen isotope fractionation.

3.3 Stable oxygen isotopes

Samples STBB I – 1 and STBB II – 7 were micromilled at 50 µm resolution along their growth length using a Sherline micromill, following the methodology outlined in Lechleitner et al. (2020), producing a total of 225 and 823 subsamples, respectively.

Stable isotope analysis was conducted at Northumbria University using a method adapted from Spötl and Vennemann (2003). 110 ± 10 µg of the sample was loaded into a 12 mL borosilicate exetainer tube, flushed with helium, and reacted with concentrated orthophosphoric acid at 70 °C. Liberated CO₂ was dried using a Gasbench II and analysed for carbon and oxygen isotope ratios. All samples were measured on a ThermoScientific Delta V Isotope Ratio Mass Spectrometer coupled with a ConFlo IV. We used an in-house laboratory carbonate standard (Plessen), alongside international standards NBS18 and IAEA603, measured every 10 samples to evaluate the runs. An in-house carbonate standard (Pol-2) was used to evaluate long-term external standard deviation, achieving < 0.1 ‰ for δ¹⁸O. All carbonate-derived stable isotope data are reported on the VPDB scale.

3.4 Trace element analysis

Trace element to calcium ($\mathrm{X}/\mathrm{Ca}$) ratios of 24 trace elements (Tables S7 and S8) were measured along the stable isotope sampling profiles of STBB I – 1 and STBB II – 7 at 2.6 µm resolution, using laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) at the University of Western Australia following the methods outlined in Finestone et al. (2022).

We prioritized elements regularly utilized as hydrological proxies including Mg, Sr, and Ba (Stoll et al., 2012; Treble et al., 2015), alongside a vegetation proxy (P), indicators for detrital input (Al, Th) and heavy metals (Fe, Mn) associated with organic matter binding (Hartland et al., 2012). See Tables S7 and S8 for a full list of elements analysed. The stable isotope record was aligned to the trace element record using multiple visual markers such as prominent growth layers and surface blemishes.

Principal component analysis (PCA) was utilized to identify common variability in the Taba Ba'astakh trace element records (Orland et al., 2014). A Gaussian kernel smoothing was applied to trace element records to remove missing data and minimize noise artefacts, before log transformation and z-score normalization to account for the PCA's sensitivity to variable scaling. The PCA was performed using the FactoMineR package in R. We applied spectral analysis to smoothed trace element and isotope records to reveal periodicities and dominant frequencies along the speleothem growth length.

4 Results

4.1 Clumped isotope temperature estimates

Clumped isotope results are summarized in Table 1. The Δ₄₇ values of our samples range from 0.626 to 0.641 for the Easotope method and 0.640 to 0.656 for the D47crunch method, returning clumped isotope temperature ($T_{{\mathrm{\Delta }}_{\mathrm{47}}}$) estimates between 6.6 and 11.1 °C, using our in-house calibration ($T_{{\mathrm{\Delta }}_{\mathrm{47}}\text{Easotope}}$) and 5.9 and 10.5 °C for the D47crunch method and Anderson et al. (2021) regression ($T_{{\mathrm{\Delta }}_{\mathrm{47}}\text{Crunch}}$). Uncertainty in Δ₄₇ is larger using the D47crunch method (standard error between 0.11 ‰ and 0.14 ‰, 95 % confidence limits between 0.022 ‰ and 0.026 ‰) than the Easotope method (standard error between 0.06 ‰ and 0.07 ‰, 95 % confidence limits between 0.012 ‰ and 0.015 ‰). This is expected given that D47crunch additionally propagates the uncertainty associated with the normalization standards (Daëron, 2021). This uncertainty is not accounted for using Easotope. Both Easotope and D47crunch produce temperatures within uncertainty of each other, with higher $T_{{\mathrm{\Delta }}_{\mathrm{47}}}$ uncertainty in the D47crunch method due to larger Δ₄₇ standard errors and the differing methods used to propagate calibration uncertainties (Sect. 3.1).

Table 1 Clumped isotope data and temperature estimates from Taba Ba'astakh samples. N = number of replicate Δ₄₇ measurements included in final Δ₄₇ calculations, with the number of excluded samples in parentheses. The number of sample replicates varies between methods due to different outlier exclusion methodologies and the exclusion of replicates with insufficient standard bracketing from the Easotope method calculations. All Δ₄₇ values are reported on I-CDES90. The fluid inclusion temperature (T_FI) from STBB I – 1 temperature range is presented for comparison with the mean in parentheses. Replicate measurement data are presented in Table S3. n/a: not applicable.

4.2 Fluid inclusion temperature estimate

Sample STBB I – 1 yielded 0.11 µL of fossil water with a δ¹⁸O value of −17.5 ‰ VSMOW and a δD value of −127.7 ‰ VSMOW, giving a deuterium excess value (d-excess) of 12.3 ‰. This value lies within uncertainty of the modern Global Meteoric Water Line (GMWL) (Craig, 1961) and ca. +9 ‰ in δ²H above the modern Local Meteoric Water Line (LMWL) from the Samoylov Island Research Station (Spors, 2018) (Fig. S2). Since large sample volumes are required for fluid inclusion measurements, our δ¹⁸O value is obtained from a speleothem section incorporating the entire growth length, representing an average of the entire speleothem. Thus, we applied a constant δ¹⁸O value of −17.5 ‰ (VSMOW) for dripwater (δ¹⁸O_dw) across the entire calcite δ¹⁸O record to obtain minimum and maximum temperature estimates. Sample STBB II – 7 failed to yield sufficient inclusion water for analysis, and thus no temperature estimate could be derived from this sample. The lower end of the STBB I – 1 fluid-inclusion-derived temperature estimate (T_FI = 9.0–19.2 °C) overlaps with the independently derived $T_{{\mathrm{\Delta }}_{\mathrm{47}}}$ estimate ($T_{{\mathrm{\Delta }}_{\mathrm{47}}}$ =  11.1).

4.3 Stable oxygen isotopes

Speleothem δ¹⁸O ranges between −17.5 ‰ and −14.6 ‰ VPDB, with mean values of −15.6 ‰ and −16.9 ‰ for STBB II – 7 and STBB I – 1, respectively (Fig. 3). STBB I – 1 δ¹⁸O exhibits an initial decline, reaching a minimum around 4 mm where the trend reverses. A sharp increase is observed around 7 mm where the record peaks and remains relatively constant thereafter. The STBB II – 7 δ¹⁸O record exhibits a similar initial decline of ca. 1 ‰ until ca. 20 mm, where this trend reverses. At ca. 37 mm δ¹⁸O declines by ca. 1 ‰ in the final ca. 4 mm of the record.

Figure 3 δ¹⁸O and selected trace element records for STBB I – 1 and STBB II – 7. δ¹⁸O data are shown with a 0.25 mm rolling window smoothing (thick lines). Trace elements are smoothed to 0.01 mm rolling window. $\mathrm{Cu}/\mathrm{Ca}$. $\mathrm{Fe}/\mathrm{Ca}$, and $\mathrm{Mn}/\mathrm{Ca}$ ratios have been mean normalized for easier comparison between the two samples.

4.4 Trace elements

Our analysis focuses on 18 elements, predominantly cation substitutes Ba, Mg, Sr, and U and metals Cu, Co, Fe, and Mn. Concentrations of Y, La, Ce, Nd, Yb, and Th were below the detection limit, and thus these elements were removed from subsequent analysis.

We identify two dominant principal components (PCs) in each sample, accounting for 55.5 % and 70.2 % of the variance in STBB I – 1 and STBB II – 7, respectively (Table 2, Fig. S4). The first PC correlates with Ba, Sr, Mg, and U and accounts for 37 % of the variance in STBB I – 1 and 15.8 % in STBB II – 7. These elements are commonly utilized as hydrological proxies in speleothems (Tooth and Fairchild, 2003; Treble et al., 2015), and we therefore refer to this PC as the hydrological PC. This PC also correlates with Na, P, and S in STBB I – 1. The second PC correlates with B, Fe, Mn, and Pb in both samples (Al, Cu, Co, Li, P, and Rb are also included in STBB II – 7) and accounts for 18.5 % of the variance in STBB I – 1 and 54.4 % in STBB II – 7. These elements are commonly associated with detrital input via flushing (Hartland et al., 2012), and we therefore refer to this PC as the detrital PC herein.

Table 2 Results of multivariate principal component analysis on Taba Ba'astakh speleothem samples. Two major elemental groupings are highlighted: the hydrological principal component, generally defined by elements Ba, Mg, P, Sr, and U, and the detrital principal component, defined by elements Al, B, Cu, Co, Fe, Mn, Rb, and Pb.

Spectral analysis revealed strong cyclicity in the hydrological PC elemental ratios (Ba, Mg, Sr, and U) in both samples, with dominant wavelengths at ∼ 0.3 and 0.5 mm in STBB I – 1 and 0.2 mm in the STBB II – 7 (Fig. S3). We also observe cyclicity in P (∼ 0.3 and 0.5 mm in STBB I – 1, ∼ 0.2–0.25 mm in STBB II – 7) and Cu (∼ 0.2 mm, in STBB I – 1 only). Detrital PC metals remain relatively constant with sporadic spikes (Fig. 3).

5 Discussion

The Late Miocene climate is widely accepted to have been several degrees warmer than today. Global SSTs were ca. 6 °C warmer than present (Herbert et al., 2016), and mean global surface temperature was ca. 4.5 °C above modern (1870–1900) (Pound et al., 2011). Most temperature estimates for the Miocene stem from marine sediments, and terrestrial data from high northern latitudes are exceedingly sparse (e.g. Popova et al., 2012). Our speleothem records give a rare insight into terrestrial high-latitude environmental changes, including temperature and hydrological conditions during this period.

Taba Ba'astakh lies deep within the modern continuous permafrost zone with a MAAT of −12.3 °C. We infer a mean annual ground temperature (MAGT) of −8.4 °C, which we calculate by averaging temperature along a 27 m borehole at the Samoylov Island Research Station (Boike et al., 2013). Since speleothem growth depends on liquid water supply, their presence at Taba Ba'astakh indicates a much warmer climate at the time of their formation (8.68 ± 0.09 Ma).

5.1 Quantitative Tortonian temperature estimates

For comparison we calculated our clumped isotope temperatures using both the Easotope software and our in-house calibration and the D47crunch data reduction algorithm with the composite calibration of Anderson et al. (2021). Both methods produce temperatures within uncertainty, giving confidence in our reconstructions. $T_{{\mathrm{\Delta }}_{\mathrm{47}}}$ uncertainty is greater for the D47crunch method due to propagation of the uncertainty in normalization standards and the different method used to calculate calibration uncertainty (Sect. 3.1). We also provide a direct comparison between clumped isotope and fluid inclusion-based temperature estimates on one sample (STBB I – 1), with the $T_{{\mathrm{\Delta }}_{\mathrm{47}}}$ estimates ($T_{{\mathrm{\Delta }}_{\mathrm{47}}\text{Easotope}}$ = 9.6 ± 6.0 °C and $T_{{\mathrm{\Delta }}_{\mathrm{47}}\text{Crunch}}$ = 11.1 ± 2.1 °C) overlapping with the T_FI estimates of both Tremaine et al. (2011) (9.0–19.2 °C, mean = 11.8 °C) (Table 1) and Coplen (2007) (10.9–21.5 °C, mean = 14.5 °C) calibrations (Table 1). A full discussion of temperature reconstruction uncertainties is given in the Supplement. Herein we limit our discussion to temperatures derived using the Easotope method and our in-house calibration since it is derived solely from subaqueous cave carbonates precipitated at temperatures similar to our samples. Both methods give similar mean temperatures across all four samples (mean $T_{{\mathrm{\Delta }}_{\mathrm{47}}\text{Easotope}}$ = +8.4 °C compared with mean $T_{{\mathrm{\Delta }}_{\mathrm{47}}\text{Crunch}}$ = +8.6 °C).

The large sample sizes required for clumped isotope analysis (ca. 4 mg) make sampling of individual growth layers impossible, and thus our temperature estimates derive multi-annual means. The long-held assumption that cave temperature corresponds to mean surface temperature (e.g. Wigley and Brown, 1976) has been questioned (e.g. Domínguez-Villar et al., 2013). Cave temperature is impacted by numerous factors including ventilation (Pflitsch and Piasecki, 2003), thermal conductivity rates through the overburden (Domínguez-Villar et al., 2013), and vegetation and snow cover (Domínguez-Villar et al., 2013; Töchterle et al., 2024), which can all act to produce offsets of several degrees between the surface and the cave environment.

Since erosion has mostly removed the Miocene overburden and brought our samples to the surface, it is impossible to know the full impact that ventilation and conduction of heat may have had on the Taba Ba'astakh cave temperatures. We suggest the insulating effects of winter snow and shading from summer vegetation likely counteracted each other with minimal overall effect. For instance, in cold regions, snow acts to insulate the ground, reducing heat loss to the atmosphere (Molnar, 2022). This insulating effect has been shown to lead to cave temperatures 5–7 °C higher than surface air temperatures in cold regions with persistent (ca. 233 d yr⁻¹) snow cover (Töchterle et al., 2024). Our $T_{{\mathrm{\Delta }}_{\mathrm{47}}}$ reconstructions between 6.6 and 11.1 °C suggest a mean annual surface temperature between modern-day Stockholm (Moberg, 2021) and London (Met Office, 2024), which experience significantly less than 233 d yr⁻¹ of snow cover. We therefore envisage a limited effect of snow insulation at Taba Ba'astakh. In addition, there is palynological evidence that the Miocene treeline stretched as far north as 80° N (Steinthorsdottir et al., 2021b), and thus it is reasonable to assume a degree of forest cover at Taba Ba'astakh during that time. Monitoring studies in Eagle Cave, Spain, showed that the transition from shrubland to forest resulted in a reduction in cave temperature up to 2 °C due to changes in insolation and modification of soil properties (Domínguez-Villar et al., 2013). Given the higher latitude of Taba Ba'astakh, it would be reasonable to assume a reduced impact from insolation shielding compared with Eagle Cave; however, a small offset (< ± 2 °C) is possible between our cave reconstructions and surface temperatures. Given global average sea levels ca. 10 m higher during the Late Miocene compared with modern levels (Miller et al., 2005), Taba Ba'astakh might have occupied a more coastal position than the present day. In the modern Arctic, lingering summer sea ice can act to reduce coastal air temperatures compared with inland locations at the same elevation (e.g. Tuktoyaktuk and Inuvik in the Canadian Arctic; Hammar, 2022) through increased albedo and latent heat effects (Vihma, 2014). Miocene Arctic Sea ice was quite reduced compared to the modern day (Stein et al., 2016), suggesting this effect may be limited. Nearby cold month temperature reconstructions from the coastal site Temmirdekh-khaj (Fig. 1) of between −2.8 and +1.1 °C (Popova et al., 2012), considerably warmer than modern levels, support this notion, although we cannot rule out the possibility of lingering cold-season sea ice reducing temperatures more than equivalent latitude inland locations.

To our knowledge, the Taba Ba'astakh temperature reconstructions are the most northerly terrestrial Tortonian MAAT estimates to date and support the consensus of a reduced latitudinal temperature gradient during the Late Miocene (Burls et al., 2021; Gaskell et al., 2022). Our temperature estimates agree well with independent high-latitude temperature reconstructions (Fig. 4). Late Miocene fossil pollen assemblages from nearby sites Temmirdekh-khaj (71° N, 132° E) and Omoloy (70° N, 133° E) (Fig. 1) yield MAAT estimates between 7.3 and 16.1 °C (Popova et al., 2012). The same authors infer a Late Miocene MAAT of 9.7 °C for the whole of eastern Siberia. Exceptionally high summer SSTs between 10 and 16 °C have been reconstructed for the Arctic Ocean at the Lomonosov Ridge (88.5° N) (Stein et al., 2016). Our temperature reconstructions provide new evidence of terrestrial MAAT ca. 18 to 23 °C warmer than present in the Siberian Arctic during the Tortonian. This is far in excess of negative temperatures simulated by models of the Alaskan Arctic during the Middle Miocene (Goldner et al., 2014) and the ca. 3.5 °C above present-day zonal mean estimated in general circulation modelling of the Tortonian (Micheels et al., 2007). Thus, our findings confirm the long-established discrepancy between model simulations and proxy reconstructions with the most northerly terrestrial Miocene proxy temperature record to date and suggest that Miocene Arctic amplification is not accurately captured by general circulation models.

Figure 4 Compiled Northern Hemisphere Late Miocene temperature reconstructions. Green filled circles are Tortonian terrestrial mean annual temperatures reconstructions from the Bradshaw et al. (2012) database, with additional sites Temmirdekh-khaj and Omoloy from Popova et al. (2012) shown unfilled. Error bars for terrestrial reconstructions show minimum and maximum estimates. The Late Miocene SST dataset from Burls et al. (2021) is shown as purple diamonds, with error bars showing uncertainty in reconstruction estimate. Additional SST estimates from the ACEX borehole are shown as unfilled purple diamonds (Stein et al., 2016). ACEX borehole temperatures show the average of the four different alkenone-derived temperature calibrations quoted in Stein et al. (2016), with error bars showing the maximum and minimum values. Both terrestrial and SST datasets have been filtered to show only data with age uncertainties overlapping our Taba Ba'astakh estimates. The red circles are the four Taba Ba'astakh clumped isotope temperature reconstructions from this study.

Since a Tortonian global mean surface temperature of 4.5 °C above modern day (Pound et al., 2011) was similar to the 2–5 °C above pre-industrial predicted for the end-of-century projections under high-emissions scenarios (IPCC, 2021), our findings provide estimates for end-of-century Arctic temperature amplification under extreme levels of anthropogenic warming. Our estimates are considerably higher than the ca. 10–12 °C of Arctic warming by the end of the century projected in modelling studies (Nazarenko et al., 2022; Xie et al., 2022) and exceed reconstructions for the Middle Pliocene Warm Period (3.3–3.0 Ma), which suggested Arctic surface air temperature warming of 7.2 °C given global mean temperatures 3.2 °C above the pre-industrial level (de Nooijer et al., 2020).

Our estimates provide useful constraints on near-future Arctic warming, but we emphasize that different planetary-scale boundary conditions in the Miocene compared with today mean our deep-time reconstructions provide imperfect analogues for anthropogenic warming. Whilst the Miocene marked a large-scale expansion of global ice sheets, both Northern and Southern Hemisphere ice sheets were highly dynamic, particularly in Greenland, which was likely only partially glaciated (Steinthorsdottir et al., 2021b). Evidence from ice-wedge pseudomorphs suggests the onset of high-latitude Northern Hemisphere permafrost formation occurred in the Late Pliocene (∼ 3 Ma) (Opel et al., 2025), and it is very likely that the Northern Hemisphere was permafrost-free during the Tortonian (Vaks et al., 2025). These, and other, slow planetary-scale feedbacks will have played a major role in the Miocene global energy budget that are unlikely to be of similar significance in driving near-future temperatures.

5.2 Stable oxygen isotope records

Both Taba Ba'astakh speleothem δ¹⁸O records show very negative values (means: −15.6 ‰ and −16.9 ‰). Modern studies have shown that dripwater δ¹⁸O (δ¹⁸O_dw) predominantly reflects changes in δ¹⁸O of precipitation (δ¹⁸O_p) above a cave (Baker et al., 2019), which, at high latitudes, is largely driven by temperature (Dansgaard, 1964). Monthly mean values of δ¹⁸O_p from both the Samoylov Island Research Station (2013–2017) and Tiksi (2004–2017) are highly correlated to mean monthly air temperatures (Spors, 2018). Thus, we interpret δ¹⁸O_p, and in turn δ¹⁸O_dw and speleothem δ¹⁸O, as indicative of temperature variations at Taba Ba'astakh.

Similarly negative δ¹⁸O values (ca. −16 ‰ to −11 ‰) are found in Greenland speleothems that formed during Marine Isotope Stage (MIS) 15a–14 (Moseley et al., 2021), a Middle Pleistocene period of unusually warm and prolonged interglacials (Rodrigues et al., 2011; Hao et al., 2015). Holocene Arctic speleothem reconstructions are absent from the literature; however, the Taba Ba'astakh isotope records agree well with modern Holocene values found in central Siberia, 15–20° latitude farther south at Botovskaya Cave (55° N), Okhotnichya Cave (52° N) (ca. −18 ‰ to −13 ‰) (Lechleitner et al., 2020), and Kinderlinskaya Cave (54° N) (ca. −14 ‰ to −11 ‰) (Baker et al., 2017). We note that the more negative δ¹⁸O of the Miocene Ocean may have compounded a shift to more negative δ¹⁸O values (Westerhold et al., 2020), but this would have been somewhat offset by reduced continentality at Taba Ba'astakh (generally associated with positive shifts in δ¹⁸O) given the ca. 10 m higher global sea level (Miller et al., 2005). Given these competing influences, we do not propose any firm assertions on the impact of source values of δ¹⁸O compared with the modern day on the δ¹⁸O signal. Thus, we assume that the stable oxygen isotope records from Taba Ba'astakh support the independent temperature estimates of a much higher MAAT during the Tortonian than at present.

5.3 δ¹⁸O signal of palaeo-dripwater

We estimate δ¹⁸O_dw using our clumped-isotope-derived temperatures in the temperature-dependent water–calcite oxygen isotopic fractionation relationship from Tremaine et al. (2011) (Table 3). Given the control of δ¹⁸O_p on δ¹⁸O_dw, the latter provides an estimate for δ¹⁸O_p. Due to the large sample sizes required for clumped isotope analysis, the samples integrate multiple growth layers, and therefore δ¹⁸O_dw values reflect multi-annual means.

Table 3 Dripwater isotopic composition estimates for the four Taba Ba'astakh speleothems. δ¹⁸O_cc is the δ¹⁸O of calcium carbonate. δ¹⁸O_cc uncertainties are standard errors. δ¹⁸O_dw uncertainty is calculated from propagation of $T_{{\mathrm{\Delta }}_{\mathrm{47}}}$ and δ¹⁸O_cc uncertainties.

δ¹⁸O_dw ranges between −17.6 ± 0.4 ‰ and −18.4 ± 0.4 ‰ (VSMOW), within error of the fluid inclusion measurement from STBB I – 1 (−17.5 ± 1.0 ‰). Agreement between these two independent δ¹⁸O_dw derivations provides additional confidence in our estimate. We stress that this assessment does not consider the influence of additional post-precipitation evaporative processes that complicate isolation of a pure precipitation signal from speleothem fluid inclusions (Lachniet, 2009).

Assuming δ¹⁸O_dw reflects δ¹⁸O_p, our reconstructed values for the Tortonian suggest a mean annual precipitation regime between the modern summer (−15.9 ‰) and autumnal regimes (−19.1 ‰) (modern annual mean = −21.6 ‰) (Bonne et al., 2020).

5.4 Seasonal hydrological regime

We propose that the cyclical behaviour of trace elements is driven by a strong seasonal hydrological regime at the time of speleothem deposition. The identified hydrological PC groups in both STBB I – 1 (PC1) and STBB II – 7 (PC2) are associated with element ratios widely utilized as hydrological proxies: most commonly $\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Sr}/\mathrm{Ca}$, and $\mathrm{Ba}/\mathrm{Ca}$ (e.g. Tooth and Fairchild, 2003; Treble et al., 2015). These alkali metals are transported via dripwaters and substituted into the carbonate lattice during speleothem deposition. In periods of low throughflow, Ca is preferentially removed from dripwaters through prior carbonate precipitation (PCP), increasing relative concentrations of Mg, Sr, and Ba.

Correlation of Sr and Mg is often used as an indicator of PCP from infiltration waters before they reach the speleothem formation site (Wassenburg et al., 2020). $\mathrm{Sr}/\mathrm{Ca}$ and $\mathrm{Mg}/\mathrm{Ca}$ exhibit no consistent relationship in either Taba Ba'astakh sample, with periods of high correlation interspersed with periods of no correlation (Fig. S5). Deviation from constant $\mathrm{Mg}/\mathrm{Ca}$ vs. $\mathrm{Sr}/\mathrm{Ca}$ ratios can arise from growth rate variability (Sliwinski et al., 2023) and changing mixing ratios of dripwater solutions derived from multiple geological endmembers (Tremaine and Froelich, 2013). Dating constraints hinder identification of growth rate variability in the Taba Ba'astakh records; however, $\mathrm{Sr}/\mathrm{Ca}$ and $\mathrm{Mg}/\mathrm{Ca}$ decoupling could arise from periods of hydrological rerouting or changing endmember dissolution rates between the carbonate and the overlying dolomitic limestone (Izokh and Yazikov, 2017) altering $\mathrm{Mg}/\mathrm{Ca}$ concentration input.

Tremaine and Froelich (2013) showed that in order to diagnose $\mathrm{Mg}/\mathrm{Ca}$ and $\mathrm{Sr}/\mathrm{Ca}$ as reliable “wet vs. dry” proxies, $\mathrm{Sr}/\mathrm{Mg}$ must remain constant along the entire sampling profile. Whilst this is not the case in the Taba Ba'astakh records, $\mathrm{Sr}/\mathrm{Ca}$ is remarkably highly correlated with $\mathrm{Ba}/\mathrm{Ca}$ (r = 0.98, p value < 0.01 in STBB I – 1 and r = 0.55, p value < 0.01 in STBB II – 7) and $\mathrm{U}/\mathrm{Ca}$ (r = 0.59, p value < 0.01 in STBB II – 7 and r = 0.83, p value < 0.01 in STBB II – 7) along the entire growth length of both Taba Ba'astakh samples (Figs. S3 and S4). $\mathrm{Ba}/\mathrm{Ca}$ has been extensively utilized as a PCP proxy (e.g. Stoll et al., 2012) and might constitute a more robust PCP proxy that is less likely affected by host rock composition at Taba Ba'astakh. Few studies have considered $\mathrm{U}/\mathrm{Ca}$ a potential indicator of infiltration changes. Dripwater U is derived from bedrock dissolution and is readily incorporated into the calcite lattice (Oster et al., 2023). The strong correlation of U with both Sr and Ba in the Taba Ba'astakh records suggests that $\mathrm{U}/\mathrm{Ca}$ may also reflect changes in PCP in the Taba Ba'astakh records. Finally, $\mathrm{P}/\mathrm{Ca}$ is also included in the hydrological principal component in STBB I – 1 and is strongly anticorrelated with Sr, Ba, and U. P is scavenged from infiltration water during calcite precipitation, thus decreasing $\mathrm{P}/\mathrm{Ca}$ in dripwater during periods of elevated PCP (Johnson et al., 2006) – as is observed in the Taba Ba'astakh records (indicated by elevated $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Ba}/\mathrm{Ca}$ and $\mathrm{U}/\mathrm{Ca}$, Fig. 5).

Figure 5 High-resolution trace element concentrations and major principal components along the STBB I – 1 and STBB II – 7 sampling profiles. The $\mathrm{Cu}/\mathrm{Ca}$ ratio has been mean-normalized for easier comparison between the two samples. Grey dashed vertical lines mark the positions of $\mathrm{Sr}/\mathrm{Ca}$ troughs.

The detrital principal component is characterized by metals which are commonly associated with organic matter flux into the cave (Hartland et al., 2012). These metals are poorly soluble and colloidally transported through organic matter binding during soil infiltration (Hartland et al., 2012). Hartland et al. (2012) demonstrated Cu and Co organic binding by showing speleothem $\mathrm{Cu}/\mathrm{Ca}$ vs. $\mathrm{Co}/\mathrm{Ca}$ ratios in line with those predicted based on n1 NICA-Donnan humic and fulvic binding affinity ratios (Milne et al., 2003). We repeat this calculation for our records and find that the measured $\mathrm{Cu}/\mathrm{Ca}$ vs. $\mathrm{Co}/\mathrm{Ca}$ ratios show good agreement with predicted ratios in STBB II – 7 (Fig. S6), confirming these metals are indeed derived from organic binding in this sample. There is poor agreement between predicted and measured ratios in STBB I – 1, suggesting detrital metals in this sample are not associated with organic binding. Phosphorus is also included in the detrital principal component in STBB II – 7. P supply has been linked with vegetative leaching during storm events (Pionke et al., 1997), with maximum $\mathrm{P}/\mathrm{Ca}$ linked to elevated infiltration rates (Borsato et al., 2007). Thus, the strong positive correlation observed between P and organically bound metals provides further evidence that these element variations reflect infiltration rates in STBB II – 7.

We present a subsection of the trace element records in Fig. 5, where correlation between Mg and other hydrological proxies is highest to illustrate the relationship between hydrological and detrital proxies. In STBB II – 7, maxima in PCP controlled elements Sr, Mg, Ba, and U coincide with minima in the infiltration-controlled elements P and Cu. Thus, wet seasons are identified as maxima in infiltration proxies ($\mathrm{Cu}/\mathrm{Ca}$ and $\mathrm{P}/\mathrm{Ca}$) and minima in the PCP proxies Sr, Mg, Ba, and U. The opposite is true for dry seasons. The relationship between detrital infiltration and PCP proxies is less apparent in STBB I – 1 (Fig. 5), where we have demonstrated that detrital proxies are not controlled by organic matter influx (Fig. S6). However, strong cyclicity is observed in all PCP proxies, as well as P, which shows a similar anticorrelation with PCP proxies to that in STBB II – 7.

Spectral analysis reveals dominant cycles corresponding to distances of 0.3 and 0.5 and 0.2 mm in STBB I – 1 and STBB II – 7, respectively (Fig. S3). Our assertion that these cycles are annual is supported by the observed alignment of the $\mathrm{Sr}/\mathrm{Ca}$ record with annual banding in STBB II – 7 (Fig. S7), identified by changes in the visual greyscale. We use free ImageJ software to extract grey values (the intensity of light carried in a single pixel) along the growth length of a high-resolution composite image STBB II – 7 thin sections (Breitenbach and Marwan, 2023). In this way, annual growth bands were identified. $\mathrm{Sr}/\mathrm{Ca}$ peaks coincide with greyscale peaks (brighter layers). Thus, our trace element periods reflect growth rates of 200 µm yr⁻¹ (STBB II – 7) and 300 µm yr⁻¹ (STBB I – 1), similar to those observed in modern temperate regions (Johnson et al., 2006; Sherwin and Baldini, 2011).

We compare Miocene and modern-day hydrological seasonality using a 20-year instrumental record of rainfall from the nearby Samoylov Island Research Station. We combined common hydrological proxies $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Ba}/\mathrm{Ca}$, $\mathrm{U}/\mathrm{Ca}$, $\mathrm{Mg}/\mathrm{Ca}$, and $\mathrm{P}/\mathrm{Ca}$ to obtain a single representative proxy-stack-average time series and aligned it to the instrumental record using dynamic time warping (Berndt and Clifford, 1994) (see Fig. S8 in the Supplement for details). While this methodology does not eliminate the considerable challenges of comparing palaeoseasonal proxy data to modern seasonality, it yields a pseudo-seasonal time axis that allows us to estimate precipitation seasonality in a comparable fashion for both sets of records. We compared the aligned modern and Miocene hydrological seasons using a spectral seasonality measure. Precipitation seasonality is quantified by the total spectral power |P|_s summed over the spectral band that closely encapsulates the seasonal peak (Fig. S9). Higher power in the seasonal band is interpreted as a more pronounced/stable seasonal signal. Detailed methodology is given in the Supplement.

The combined hydrological proxy records suggest a stronger seasonal precipitation regime during the Tortonian (i.e. a more pronounced seasonal cycle) compared with modern-day conditions (Fig. 6). This is true of both samples, with significantly higher median spectral power in the seasonal band in the Miocene compared with the modern level. We recognize the limitations in the modern record, which, whilst it is the closest available modern precipitation time series, measures only the liquid fraction and not snowfall, which constitutes ca. 30 % total precipitation at the Samoylov Research Station (Boike et al., 2013). However, exclusion of the winter snow component within the modern record will create an apparent stronger seasonal contrast. Given that our Tortonian hydrological proxy record exhibits a stronger seasonality measure than this modern rainfall record, we argue this reinforces our finding of a more seasonal hydrological regime in the Tortonian.

Figure 6 Reconstructed precipitation seasonality (|P|_s summed over the seasonal peak spectral band) in Tortonian samples STBB I – 1 and STBB II – 7, compared with the modern day. An alignment procedure based on dynamic time warping allows for their comparison. Coloured points show seasonality values for different widths of the seasonal band between 0.2 (0.9–1.1) and 0.5 (0.75–1.25). Black points show the median over all values. Both samples show overall enhanced seasonality in the Miocene compared with the modern day.

Numerous studies have suggested the Miocene was a time of enhanced seasonality in the mid latitudes with warm–wet, cold–dry seasonal cycles reconstructed in central Europe (Bruch et al., 2011; Harzhauser et al., 2011). At high latitudes, the picture is less conclusive. Carpological reconstructions from eastern Siberia suggest the Late Miocene marks the onset of modern atmospheric circulation patterns and the establishment of today's wet summer/dry winter regime, albeit with considerably higher total precipitation between 700 and 900 mm yr⁻¹ (Popova et al., 2012). At Temmirdekh-khaj and Omoloy, Miocene carpological reconstructions suggest the mean annual range in precipitation (the difference between wet and dry season amplitude) was 70–80 and 110–130 mm, respectively. This seems large compared with today's mean summer precipitation of 169 mm and annual precipitation of 309 mm but should be viewed in context of a wetter Tortonian regime given higher moisture availability from a summer ice-free Arctic (Stein et al., 2016). Analysis of fossil wood δ¹³C from Cherskiy (67.7° N, 161.6° E) in northeastern Siberia shows a high variability in precipitation, with the wet season alternating between summer and winter (Schubert et al., 2017). In contrast, our Taba Ba'astakh records suggest an enhanced seasonal precipitation regime compared with the modern day. We propose this likely arose from differing moisture availability from the nearby Arctic Ocean, which was largely sea-ice-covered during the Tortonian winter and ice-free during Tortonian summers (Stein et al., 2016). This would have led to enhanced evaporation and moisture transport into northern Siberia in the summer, compared with winter. This effect appears to outweigh any impact of reduced continentality given global average sea levels ca. 10 m higher during the Late Miocene compared with modern levels (Miller et al., 2005). A likely return to similar ice coverage over the coming century (Sigmond et al., 2018) may drive a shift to more enhanced seasonal precipitation regimes similar to those we infer for the Tortonian.

6 Conclusions

We present the first speleothem-derived temperature and precipitation seasonality records from Arctic Siberia during the Tortonian (11.63–7.246 Ma). Our findings offer insight into high-latitude climate during a period when the global mean surface temperature was similar to that projected over the coming decades under high-emissions scenarios.

Apparent hydrologically driven element cycles suggest a more seasonal precipitation regime compared with the modern day, likely driven by fluctuations in Arctic Sea ice extent. Given probable reductions in summer Arctic ice extent over the coming decades, a return to more volatile precipitation regimes seems likely, with increased summer precipitation driving permafrost thaw and impacting local infrastructure.

We use two independent methodologies to reconstruct the most northerly terrestrial Tortonian temperatures between 6.6 and 11.1 °C from high-latitude central Siberia. These findings provide much needed estimates for future Arctic amplification, implying Arctic MAAT increases of 18.9–23.4 °C above modern levels, given global mean temperature rises of ca. 4.5 °C above the present day. This is in good agreement with previous estimates from the region and supports the consensus of a substantially warmer Late Miocene climate compared with the modern day. Our findings confirm the long-established discrepancy between model simulations and proxy reconstructions, producing considerably higher estimates for Tortonian Arctic amplification than models.