Limited early-industrial warming and strong volcanic imprints in the Caucasus: the first temperature reconstruction based on maximum latewood density

The Caucasus occupies a unique climatic region influenced by European, Mediterranean, and Asian circulation systems, yet it remains underrepresented in tree ring-based Northern Hemisphere temperature proxy networks. Here, we present the first summer temperature reconstruction for the Caucasus region based on maximum latewood density (MXD). We used X-ray micro-computed tomography of tree-ring samples from Pinus sylvestris growing at the upper tree line in the Lesser Caucasus and an ensemble nested regression approach to develop a robust 326 year-long June–September temperature reconstruction (1697–2022). The record explains – regionally unprecedented – 72 % of temperature variance during the instrumental period (1901–2022) and captures distinct interannual and multi-decadal variability including pronounced warming since the 1990s and a strong imprint of major volcanic eruptions. Temperatures in the second half of the 19th century were not significantly colder in the Caucasus than in the second half of the 20th century. The reconstruction highlights the exceptional magnitude and pace of 21st century warming in the region, which is without analogue at least in the past three centuries. Comparison with regional and large-scale temperature reconstructions reveals strong agreement within the Caucasus but negative correlations with Central Europe, indicating distinct temperature variability patterns across Europe and western Asia. Future work should focus on the climate dynamics behind this dipole and the extension of temperature-sensitive tree-ring records in the region.

1 Introduction

Global surface temperatures have risen rapidly over the past century, reaching levels unprecedented in the context of the last two millennia (IPCC, 2023). Understanding this recent warming requires a long-term perspective on natural climate variability and the drivers of past temperature change. High-resolution paleoclimate reconstructions extend our view into the Common Era and provide essential context for evaluating the magnitude, spatial patterns, and causes of recent warming (Esper et al., 2018). Such reconstructions have revealed both natural centennial-scale fluctuations, including the Medieval Climate Anomaly and the Little Ice Age, as well as the exceptional nature of the temperature increase since 1900 CE (Anchukaitis et al., 2017; Esper et al., 2024; Neukom et al., 2019).

Most of the proxy records used for large-scale, high-resolution temperature reconstructions are derived from tree rings, which provide a superior archive of terrestrial climate variability for timescales from years to millennia because of their precise annual age control and the relatively extensive spatial coverage (Büntgen et al., 2021; Esper et al., 2018). Among different parameters derived from tree rings, maximum latewood density (MXD) provides the most reliable proxy of summer temperature with usually high correlations and little impact of physiological memory (Briffa et al., 1998; Esper et al., 2015). Over the past few decades, significant efforts have been made to develop a network of MXD records across the Northern Hemisphere to better understand large-scale and long-term temperature changes (Anchukaitis et al., 2017; Schneider et al., 2015).

In southern Europe, MXD-based temperature reconstructions are available for the Pyrenees (Dorado Liñán et al., 2012), the European Alps (Büntgen et al., 2006), the Western Carpathians (Büntgen et al., 2013), Croatia (Römer et al., 2021) and Greece (Esper et al., 2020; Klesse et al., 2015). In southern Asia, MXD chronologies were developed for the Altai Mountains (Kirdyanov et al., 2024), the Tien Shan (Chen et al., 2019), the Tibetan plateau (Huang et al., 2025) and the Himalaya (Khan et al., 2024). However, along the Alpide belt, a significant spatial gap in the MXD network exists between the Euro-Mediterranean and the Central-/East-Asian regions (Anchukaitis et al., 2017)

A few studies have used other dendrochronological proxies to reconstruct past temperature variability in the Caucasus (Bakhtiyorov et al., 2025; Dolgova, 2016; Holobâcă et al., 2016; Köse et al., 2017; Kvaratskhelia and Gavashelishvili, 2025). Holobâcă et al. (2016) reconstructed 179 years of summer temperatures using tree-ring width (TRW) data, but with a relatively weak temperature signal. Dolgova (2016) produced the first summer temperature reconstruction extending back to 1596 using Latewood Blue Intensity (LWBI). The record correlates strongly with observed temperatures and shows unprecedented warming since the 1980s. Compared to Dolgova (2016), more recent reconstructions were shorter and less skillful (Bakhtiyorov et al., 2025; Köse et al., 2017; Kvaratskhelia and Gavashelishvili, 2025). A network of temperature sensitive LWBI chronologies revealed high coherence between the Lesser and the Greater Caucasus at high frequencies, but divergence on multi-decadal and centennial timescales (Dhyani et al., 2025). However, this disagreement is likely not an expression of local temperature variability but rather associated with the uncertainty inherent in the low-frequency signal of the LWBI proxy (Björklund et al., 2017). Consolidating pre-industrial temperature estimates for the Caucasus requires a more robust proxy, like MXD, which can offer greater precision at longer timescales.

One primary constraint for the widespread application of MXD is the limited availability of specialized measurement equipment and the complexity and time-intensive nature of traditional X-ray densitometry methods (Björklund et al., 2019). Recently, X-ray Computed Tomography (XµCT) has emerged as an innovative alternative, offering precise density measurements. Simplified sample preparation and a fully digital workflow after X-ray imaging minimize the labor costs and allow the production of well-replicated chronologies. The three-dimensional approach overcomes challenges associated with varying grain angles in wood samples affecting MXD and enhances the robustness of resulting chronologies (Björklund et al., 2019; De Mil and Van den Bulcke, 2023; Verschuren et al., 2025).

While previous LWBI chronologies from the Caucasus have been valuable and consistent regarding high-frequency temperature variability (Dhyani et al., 2025; Dolgova, 2016), longer-term trends diverged significantly which can partly be attributed to limitations of the proxy. In this context, we aim to answer whether MXD data can provide a more robust basis for the reconstruction of past temperatures. We hypothesise that a high sample replication combined with the new XµCT approach yields more reliable chronologies and a more solid proxy calibration to instrumental data at high and low frequencies resulting in a better confined temperature history at local to regional scales.

To test this hypothesis, we developed a 409 year chronology from Scots pine (Pinus sylvestris) tree rings from the upper treeline in the Lesser Caucasus. Our XµCT-derived MXD record is assessed with established quality metrics and correlated to local instrumental temperature data at different frequency domains. After conversion to temperature anomalies, it is compared to temperature variability during the last 300 years over the wider region and to other tree-ring reconstructions from mid-latitude high-elevation sites across Europe.

2 Material and methods

2.1 Climate of the study region

The Caucasus region coincides with the climatic transition between the humid subtropics of the eastern Black Sea and the more continental, semi-arid interior of Eurasia (Ambroziak, 1977; Stanev, 1990). Moisture is supplied year-round by prevailing westerly circulation, but the rugged topography creates sharp regional gradients. The Greater Caucasus limits the southward penetration of cold continental air, while the Lesser Caucasus redirects warm southerly flows (Ambroziak, 1977). Seasonal climate variability is controlled by the migration of major pressure systems: in summer, the subtropical high pressure over the eastern Mediterranean and the Asian thermal depression promote continental warming and convective activity; in winter, the western extension of the Siberian High dominates, modulating cold-air outbreaks and storm tracks (Shahgedanova et al., 2009). Recent climatic observations from the Caucasus reveal that the signal of summer warming is particularly pronounced, exceeding the rates of temperature increase recorded in winter, spring, and autumn (Elizbarashvili et al., 2013, 2017; Shahgedanova et al., 2005, 2009; Tashilova et al., 2016; Tielidze et al., 2020). As in many mountain regions, the Caucasus lacks long-term meteorological records at high elevations. Moreover, reliable station measurements in the region began operating in the mid-20th century and some were not maintained after the collapse of the Soviet Union (Hu et al., 2014; Keggenhoff et al., 2014; Schöne et al., 2013).

2.2 Tree-ring sampling and preparation

In two field campaigns (2015 and 2023), we collected tree core samples from Pinus sylvestris trees growing at the upper treeline (41.4° N, 42.3° E; 2200–2300 m above sea level (m a.s.l.)) in the western Lesser Caucasus (Fig. 1). The sampling site is characterized by an open-canopy P. sylvestris forest growing on dry, south-facing slopes in a subalpine environment with rocky outcrops (Martin-Benito et al., 2018). The underlying pedology consists predominantly of shallow, rocky soils derived from volcanic and sedimentary parent materials. The sparse understory features rhododendron and juniper shrubs separated by patches of herbaceous vegetation. A total of 73 increment cores were obtained from 55 trees by extracting one or two cores per tree perpendicular to the slope. To avoid the influence of resin on MXD measurements, samples were chemically treated prior to the analysis. A Soxhlet extraction with ethanol for at least 24 h was performed to remove resins for accurate wood density determination (De Mil and Van den Bulcke, 2023). After extraction, the samples were oven-dried at 103 °C for 24 h and then stored before scanning in paper straws and air-tight bags with silica gel to prevent them from re-moisturization.

Figure 1 (a) Map of the study site showing locations of the sampling site, previous studies on LWBI by Dolgova (2016) and Bakhtiyorov et al. (2025), nearby metrological stations and CRU grids used in the study. (b) Pinus sylvestris trees growing at the upper treeline (2200–2300 m a.s.l.) in the Lesser Caucasus.

2.3 MXD measurements

We scanned the cores using the XµCT toolchain from the UGent-Woodlab at Ghent University (De Mil et al., 2016; De Mil and Van den Bulcke, 2023; Van den Bulcke et al., 2014, 2019; Verschuren et al., 2025). The cores were scanned at a voxel size of 15³ µm³ with a helical scan using the TESCAN CoreTOM. We chose the scan setting “High resolution fast” from Verschuren et al. (2025). The Octopus Reconstruction software reconstructed the voxels into 3D volumes (Vlassenbroeck et al., 2007). The reconstructed scans were processed with three different MATLAB stand-alone applications designed for the analysis of 3D wood density data: CoreProcessor, RingIndicator, CoreComparison (Verschuren, 2025b). Images of the individual cores were extracted and aligned to the vertical fiber direction using the CoreProcessor, which also converts grayscale values into absolute density values. Ring boundaries are marked as lines on both the transverse and radial planes of the 3D digital increment core using RingIndicator (Verschuren, 2025b) Cores were visually and statistically cross-dated using the CoreComparison (Verschuren, 2025b) and the COFECHA software (Holmes, 1983) based on obtained TRW measurements. For each sample, a continuous density profile is generated along the entire length of the core and integrating over 30 % of the core radius. Yearly MXD values were then calculated with the density profile for the respective ring using the average of the 20 % densest values within the last 20 % of the radial width. This step, performed with the XCT.Read function (Verschuren, 2025a) in R (R Core Team, 2016), removes the effect of resin ducts and yields a clear separation towards the earlywood of the next year. Additionally, visual and statistical cross-dating based on the obtained MXD measurements has been performed in CDendro (https://cdendro.se/, last access: 16 July 2025) to ensure correct cross-dating of the samples that were difficult to cross-date solely based on TRW measurements.

2.4 Climate data

We analysed monthly records of temperature variability from June–September using the nearest station observations: Ardahan and Artvin (Fig. S1 in the Supplement) (Menne et al., 2012). Annual JJAS means were calculated by averaging June–September values, excluding years with missing or flagged data. However, the station records were relatively short and incomplete, which limited their application to long-term climate analyses. Thus, we complemented the station data with the three nearest gridded temperature data points from the CRU TS 4.09 dataset, which is a pre-processed product derived through interpolation of monthly climate anomalies from extensive weather station observations providing continuous monthly mean temperatures from 1901 to the presence (Harris et al., 2020). We compared individual station records with corresponding grid cells and analysed the correlation between the mean of all stations and the mean of all grids (Fig. 1). Correlation analysis demonstrated strong agreement between stations and gridded data (Fig. S1). Recent warming trends were also comparable over the overlapping period, with CRU JJAS temperatures increasing at nearly the same rate as the station-based means. This consistency in both interannual variability and long-term trends demonstrate the reliability of the CRU TS 4.09 dataset. For all the subsequent analyses we used the averaged time-series extracted from three grid points of the CRU TS 4.09 temperature dataset as the instrumental temperature record.

2.5 Chronology Development

Standardization, which includes detrending and averaging of raw MXD data, was implemented to minimize the effects of non-climatic growth trends on the resulting chronology while preserving the full frequency spectrum of climate variability. We detrended the individual MXD series taking its ratios from fitted cubic smoothing splines. The age-dependent splines had an initial stiffness of 50 years (50 % variance cut off at the wavelength of 50 years) and progressively increasing stiffness with cambial age (Cook and Peters, 1997; Melvin et al., 2007). To further enhance the capacity of the chronology to retain long-term signals, the detrending was performed within the signal-free framework (Melvin and Briffa, 2008) using the ssf() function in dplR (Bunn, 2008). This approach iteratively removes the common signal prior to detrending, allowing the detrending splines to be more accurately estimated without being influenced by strong climate signals (Melvin et al., 2007). Chronology quality was assessed using subsample signal strength (SSS) (Buras, 2017) and average inter-series correlation ($\stackrel{\mathrm{‾}}{R}$) (Wigley et al., 1984).

2.6 Assessment of the climate signal in MXD

Pearson's correlations between the instrumental temperature record and the detrended MXD chronology were calculated for the period 1901–2022, with monthly temperature data for September of the previous year to October of the current year and the June–September (JJAS) seasonal means. A 30 year running correlation analysis was performed to check the stability of the correlations over time. The impact of sample replication was assessed with a bootstrap resampling analysis. For each replication level (ranging from 1–73 MXD tree-ring series), 1000 iterations were performed by randomly sampling and averaging tree-ring series, followed by correlation with the JJAS temperature record.

2.7 Temperature reconstruction

For reconstructing summer temperature, we applied a nested approach that builds successive transfer functions from a set of 73 individual detrended MXD chronologies with increasing temporal coverage and decreasing replication (Cook et al., 2002; Meko, 1997). Within each “nest” the replication of the chronology, i.e., the number of MXD measurements per year, remains constant. Each time a single sample drops out, a new nest chronology is built from the remaining samples. The first nest chronology spans the period from 1901–2022. The algorithm generated 61 different nests by sequentially eliminating the shortest sample at each iteration, thereby producing nested subsets with progressively earlier start dates but decreasing replication. This nesting approach ensures optimal utilization of available MXD data across the entire chronology period and allows for constant re-assessment of reconstruction skill back in time. Each nest chronology was calibrated against the instrumental temperature record (1901–2022) using a linear regression model and an iterative split-period approach with a fixed calibration window of 61 years. The first reconstruction model was calibrated on observations between 1901 and 1961 and validated against temperatures between 1962 and 2022. We then iteratively shifted the 61 year calibration/validation block one year forward at a time while re-calibrating the reconstruction with the remaining years of climate data. For each model, verification statistics including the coefficient of determination (R²), the reduction of error (RE) and the coefficient of efficiency (CE) were computed to assess reconstruction skill (Fritts, 1976). A final reconstruction nest was developed by averaging all the iterative model outputs. The reconstruction nests were spliced together using the best replicated nest at each time step. To assess the robustness of the reconstruction, we quantified two distinct sources of uncertainty: calibration uncertainty and sampling uncertainty (Esper et al., 2007; Kuhl et al., 2024; Schneider et al., 2015). Calibration uncertainty represents the statistical error in the linear regression between the MXD chronologies and summer temperature during the calibration period and was calculated using the residual standard error of the regression coefficients for each nest. Sampling uncertainty was estimated by expressing it as the additional error variance of each nest compared with the best-replicated nest that is, the residual between the standard error of the current nest and the standard error of the nest with the highest replication (Cook et al., 2002; Kuhl et al., 2024). We correlated reconstructed summer temperatures for the Lesser Caucasus with the gridded temperature reconstruction fields of the Northern Hemisphere Tree-Ring Network Development (N-TREND) from Anchukaitis et al. (2017) over their period of overlap (1697–2014) and compared our record to other regional products (Table S1 in the Supplement). These were z-transformed and low-pass filtered in order to emphasize similarity in the low-frequency domain calculated using 30 year moving averages. For the analysis of volcanic signals, temperature reconstructions were re-calibrated with summer months from local gridded instrumental temperature data and a scaling for each of the reconstructions (Esper et al., 2005), because of small deviations in the targeted season. We assessed the relative expression of high- and low-frequency temperature variability in the different local reconstructions by calculating the ratios of standard deviations in high- and low-pass filtered timeseries (σ_HF and σ_LF).

2.8 Superposed Epoch Analysis

We evaluated the influence of large-scale tropical volcanic eruptions on reconstructed summer temperatures using Superposed Epoch Analysis (Rao et al., 2019). Significant tropical eruptions during the 18th and 19th century (1809, 1883, 1835, 1831, 1815) were identified from the eVolv2k database by selecting those with peak estimated Northern Hemisphere sulfate flux ($>\mathrm{10}\mathrm{kg}{\mathrm{km}}^{-\mathrm{2}}$) and non-zero Southern Hemisphere flux (Toohey and Sigl, 2017). In addition, the 1991 Pinatubo eruption was added to the event list. To determine the statistical significance of the observed response, we compared it to the probability distribution generated from 10 000 random draws of six years from the reconstruction (Rao et al., 2019).

3 Results

3.1 MXD chronology and climate

We developed a 326 year (1697–2022) set of 61 MXD chronologies with decreasing replication and increasing length. For the best replicated period, average inter-series correlation is 0.43 and SSS is 0.97 (Fig. 2). The best replicated MXD chronology shows significant positive correlations with temperatures from February–October (Fig. 3a). Distinctively higher correlations are confined to the months from June to September, with the highest correlation for mean JJAS temperature (r=0.84, p<0.001). In the high-frequency domain (i.e., first differences), the highest correlation value was also obtained with the mean JJAS temperature (r=0.78) (Fig. 3c), proofing a very strong relationship between temperature and proxy even if autocorrelation and trends are removed from the data. The moving-window correlation analysis between the MXD chronology and the JJAS mean temperature revealed stable correlations (p<0.01) over the entire instrumental period (1901–2022). The bootstrapped resampling analysis shows the mean correlation increased steadily with replication, from 0.56 (95 % CI: $r=\mathrm{0.25}-\mathrm{0.79}$) at $n=\mathrm{1}-\mathrm{0.84}$ (95 % CI: $r=\mathrm{0.83}-\mathrm{0.83}$) at n=73 (Fig. S3a in the Supplement). Mean correlation after n>10 was 0.78 (CI: $r=\mathrm{0.70}-\mathrm{0.85}).$ After n>30, mean correlation values approach an asymptote around r=0.84 (Fig. S3a).

Figure 2 Nested MXD chronology and chronology statistics. (a) Detrended maximum latewood density (MXD) chronology (green) with 20 year spline smoothing (brown) and corresponding sample depth (grey, secondary axis). The vertical dashed grey lines show end years of individual nests. (b) Mean inter-series correlation ($\stackrel{\mathrm{‾}}{R}$) and Subsample signal strength (SSS) through time. The blue dashed horizontal line in panel (b) denotes the threshold of SSS=0.85.

Figure 3 MXD chronology vs. JJAS CRU TS 4.09 temperature relationships. (a) Correlation coefficients between the MXD chronology and monthly/seasonal temperatures, with significance levels indicated by asterisks (** p<0.01, *** p<0.001). (b) regression between MXD chronology and JJAS temperature for the period 1901–2022 and (c) Z-scores of the MXD chronology (orange) and JJAS temperature data (blue) with 20 years smoothing spline. The correlation of raw MXD and first difference (Δmxd) with JJAS temperature are shown in the figure.

3.2 Reconstruction skills and characteristics

Our MXD-based summer temperature (June–September) reconstruction for the Lesser Caucasus (LC) demonstrates strong predictive skill, explaining a maximum of 72 % of the variance in summer mean temperatures over the instrumental period from 1901–2022 (Fig. 4). During the shifting calibration periods, the model explained a mean of 65 % of the variance in temperature data. In the verification period, the model retained a high level of skill, with 68 % of explained variance, showing that the reconstruction generalizes well beyond the calibration dataset. RE values were positive across all calibration windows representing high predictive skill beyond climatological baselines, while R² values remained consistently high, indicating stability of the temperature signal (Fig. S4 in the Supplement). Across all moving calibration windows, the reconstruction achieved a mean positive RE skill (0.67). The similarity of RE and R² values across nests shows that our reconstruction reliably estimates JJAS temperatures even for the longest nests with relatively low replication. Nests with high replication (≥10 predictors), calibrated in the 20th century, achieve robust predictive performance, with median CE commonly above 0.6 and R² exceeding 0.7. Skill gradually degrades as replication falls: nests with only 3–8 series still maintain mostly positive CE (0.3–0.5; Fig. S4). The best replicated nests show smaller residual errors, with mean root mean squared error (RSE) around 0.4 °C. In contrast, least-replicated nests with the minimum predictor set at the oldest part of the record show larger calibration error with RSE exceeding 0.7–1.0 °C. Monte Carlo simulations demonstrated that the error declines rapidly with increasing replication, with the largest improvements achieved at n>10, while gains diminish once replication is relatively high (Fig. S3b).

Figure 4 Reconstructed and observed JJAS temperature anomalies (°C) for 1697–2022 relative to 1961–1990. The LC reconstruction is shown as a thin dark-green line with a 20 year spline (bold red line); instrumental temperature appears as the thin purple line with a 20 year spline (bold purple line). Calibration uncertainty is conveyed by the inner ±2σ model band (light green) and the sampling uncertainty is shown by outer nest 5 %–95 % band (yellow). Red and blue circles mark the 10 hottest and coldest reconstructed years, respectively. The bottom panel shows median RE, CE, and R² across nests and the dashed horizontal line indicates zero skill.

Reconstructed temperatures varied strongly from year to year, while multi-decadal and centennial-scale variability was relatively low (Fig. S7 in the Supplement). At the annual scale, the ten coldest years spread over the last three centuries with three each in the 18th (1700, 1701, 1781) and 19th (1810, 1816, 1817) century and four cold extremes in the 20th century (1914, 1956, 1959, 1967), while the nine warmest years (2008, 2010, 2015, 2017, 2018, 2019, 2020, 2021, 2022) are exclusively concentrated in the early 21st century. The coldest sustained 30 year period occurred during 1809–1838 with a mean temperature anomaly of −0.66 °C (Fig. 4). This interval includes several consecutive negative extremes, with 1831–1840 representing the coldest decade of the record. Between the 1850s and the 1990s, the reconstruction shows no clear long-term warming trend, with 30 year mean anomalies in the late 20th century remaining comparable to those in the mid-19th and early 20th centuries. At the recent end, the reconstruction shows pronounced warming of ∼0.52 °C per decade (p<0.001) (Fig. 4).

3.3 Regional climate variability

To highlight the spatial representativeness of our LC reconstruction, the correlation between the reconstruction and the observed JJAS mean temperature field was calculated over the eastern Mediterranean domain and southwestern Asia. Correlations between the LC reconstruction and gridded summer temperatures are highly significant over a wide area that includes North Africa, the Arabian Peninsula and the Eurasian steppe region (Fig. S5 in the Supplement). After high-pass filtering, the correlations remain high over the Caucasus region and particularly towards the south (Arabian Peninsula) and east across the Eurasian Steppe. Comparison between our LC reconstruction and other MXD temperature reconstructions from the low- and mid-latitudes (Fig. 5a) shows strong correlation with the Greater Caucasus (r=0.64, p<0.001), a weak correlation with the Eastern Mediterranean (r=0.26, p<0.001) and insignificant correlations with the Pyrenees (r=0.08) and the Alps (r=0.06). After 20 year low-pass filtering, positive correlations ($r=\mathrm{0.69}-\mathrm{0.76}$) are found with all the considered reconstructions. At high frequencies, a declining correlation from east to west is observed (r=0.64 to −0.24) (Fig. S6 in the Supplement). Similar correlation trends are also found with N-TREND data at higher frequencies (Fig. 5b). The low-to-high-frequency variability ratios for the 18th and 19th centuries are low in the LC reconstruction (${\mathit{\sigma }}_{\text{LF}}/{\mathit{\sigma }}_{\text{HF}}\approx \mathrm{0.14}-\mathrm{0.23}$; Fig. S7). In comparison, higher ratios are observed in the Alps (${\mathit{\sigma }}_{\text{LF}}/{\mathit{\sigma }}_{\text{HF}}\approx \mathrm{0.42}-\mathrm{0.54}$; Fig. S7) and the Pyrenees (${\mathit{\sigma }}_{\text{LF}}/{\mathit{\sigma }}_{\text{HF}}\approx \mathrm{0.21}-\mathrm{0.36}$; Fig. S7).

Figure 5 (a) Comparison of the LC reconstruction with other MXD records from the mid latitudes. All the series are smoothed with a 30 year smoothing spline. Colours are consistent across panels for the Lesser Caucasus, the Greater Caucasus, the Eastern Mediterranean, the Pyrenees and the Alps. The top-left inset map marks site locations using the same colours. (b) Spatial correlations of the LC reconstruction vs. N-TREND temperature field (Anchukaitis et al., 2017) and (c) is as (b), but with first differenced data for the period 1697–2014. Blue filled circles indicate the site location of the present study.

3.4 Response to major volcanic eruptions

The SEA showed strong post-volcanic cooling (+1 year), with the strongest anomalies in our present LC reconstruction (−1.11 °C, p<0.01; Fig. 6). The other reconstructions showed similar behaviours in the SEA but with different magnitudes and timing (Fig. 6). Significant post-eruption cooling occurred one year after eruptions in the Lesser and Greater Caucasus. The records from the Eastern Mediterranean and the Pyrenees showed significant cooling only in year 0. The cooling estimates in the years 0 and +1 are similar in magnitude and range from 0.83–1.00 °C. Remarkably, the LC reconstruction thereby exceeds its significance threshold twice, much more than any of the other records, indicating that volcanic cooling is particularly pronounced compared to background variability in this reconstruction. The record from the Alps showed no significant cooling in the years 0 to +5 (Fig. 6).

Figure 6 Superposed epoch analysis (SEA) showing the response of tree-ring reconstructed temperatures to volcanic events. Solid lines represent the mean temperature anomaly for each record relative to five pre-event years. The horizontal dotted lines denote the thresholds for statistical significance determined using random bootstrapping at the 95 % confidence level. These thresholds were calculated by compositing 10 000 random draws of 6-time windows from each temperature reconstruction.

4 Discussion

4.1 Limited early-industrial warming and strong volcanic forcing in the LC reconstruction

We present the first MXD chronology from the upper treeline in the Caucasus region, developed from 73 individual tree core samples of slow-growing living pine trees using the innovative X-ray µCT technology for density measurements. Our study was designed to maximize the climatic signal in the proxy data. The sampling site is characterized by open-canopy pine stands at the upper treeline (>2200 m) growing on dry south-facing slopes naturally minimizing ecological competition (Fig. 1). The high replication, a moderate age trend and a robust inter-series coherence of the MXD data are strong indicators for a reliable chronology. Weak age-trends and minimal differences between raw and detrended MXD chronologies imply a limited impact of detrending on the low-frequency climate signals (Fig. S2). Compared to the previously developed LWBI chronology from the same site (Dhyani et al., 2025), the MXD reconstruction shows a significantly better fit with observed temperatures, particularly in the low-frequency domain (Fig. S8 in the Supplement). This reinforces the capability of MXD to better retain long-term climate signals (Björklund et al., 2024). In addition, the bootstrap resampling experiment further supports the stability of the relationship between MXD and temperature even at lower replications (Fig. S3a) indicating that MXD can provide robust reconstructions even with low to moderate sample replication. Collectively, these attributes reduce the overall uncertainty of the reconstruction, even when conservative assumptions are applied regarding error propagation.

The MXD chronology aligns exceptionally well with the instrumental temperature record over the period of anthropogenic warming, showing no indication of divergence during the rapid warming of the late 20th and early 21st centuries (D'Arrigo et al., 2008). The most striking features of our LC reconstruction are pronounced short-lived volcanic cooling signals, limited low-frequency variability until the 1990s and an unprecedented 21st century warming. Extreme warm years, including 2010, 2018, 2019, and 2022, coincide with large-scale Eurasian heatwaves (Beck et al., 2013, 2024; Jeong et al., 2025; Rasmijn et al., 2018) and record-breaking global warmth (Dunn et al., 2023). The reconstruction captures pronounced cold anomalies in 1810, 1816 and 1817 associated with major ice-core sulphur deposition peaks (1809 and 1815) and with the “Year Without a Summer” following the 1815 Tambora eruption (Stothers, 1984; Toohey and Sigl, 2017), suggesting a pronounced impact of tropical volcanic forcing in the region. Other tropical eruptions, like Krakatau in 1883 and Pinatubo in 1991, also result in significant cooling one summer after the eruption. In contrast, the eruption of Cosigüina in January 1835 (Toohey and Sigl, 2017) caused cooling in summer of the same year because of the early eruption date (Esper et al., 2025; Kuhl et al., 2024). The short-term volcanic cooling found in the LC reconstruction, exceeding 1.0 °C on an absolute scale and 0.6 times the standard deviation on a relative scale, is larger than in many comparable European reconstructions, yet volcanic eruptions do not induce the same level of multi-year to decadal cooling that is typical for regions coupled to North Atlantic dynamics (Wanner et al., 2022). These differences likely show the continental climatic setting of the Caucasus, where reduced ocean–atmosphere coupling limits the persistence of temperature responses to forcing pulses. Additionally, cold extremes in the mid-20th century occurred in 1956, 1959 and 1967 coinciding with Eurasian blocking events (Fink et al., 2006).

Beyond these short-lived events, the LC reconstruction shows remarkably limited long-term temperature variability prior to the late 20th century. The early 19th century (1807–1836) is the coldest multi-decadal interval (30 year mean) in the reconstruction, potentially related to the Dalton Minimum, a period with prolonged solar activity minima (Silverman and Hayakawa, 2021) as well as to repeated volcanic cooling events (Schneider et al., 2017). In contrast with other single site and regional MXD reconstructions, such as those from the Alps (Esper et al., 2025), Southern Finland (Helama et al., 2014), Polar Urals (Briffa et al., 2013), and large-scale reconstructions (Anchukaitis et al., 2017; Büntgen et al., 2021), the LC reconstruction shows limited early-industrial warming. Temperatures during 1950–2000 were only slightly warmer (ΔT∼0.1 °C) than during 1850–1900. Weak early-industrial warming is also evident in other regional or Mediterranean reconstructions (Esper et al., 2020; Kuhl et al., 2024). While limited low-frequency variability can arise from methodological choices in proxy datasets (Büntgen et al., 2021), the strong calibration statistics and absence of divergence in the LC record suggest that this pattern may reflect genuine regional climatic stability.

4.2 Regional climate variability in the Caucasus (and Black Sea) region

To date, few studies have utilized TRW and LWBI data to reconstruct temperature variability from the Caucasus region (Bakhtiyorov et al., 2025; Dolgova, 2016; Holobâcă et al., 2016; Köse et al., 2017; Kvaratskhelia and Gavashelishvili, 2025). The variance explained by summer temperature is much higher in our MXD-based LC reconstruction (72 %) than reported for other chronologies from the region (40 %–65 %) (Bakhtiyorov et al., 2025; Dolgova, 2016), and even exceeds values from many highly cited MXD-based temperature reconstructions (Esper et al., 2016). The spatial correlation fields with observed JJAS temperatures (Fig. S5) show the relevance of the LC reconstruction for revealing the temperature history of a wider region. Our LC reconstruction, together with the Greater Caucasus records from Dolgova (2016) and Holobaca et al. (2016) agree with the finding that recent warming clearly exceeds the amplitude of earlier warm periods. In contrast, the recently published reconstruction from the northern Lesser Caucasus (Bakhtiyorov et al., 2025) shows warming since the 1990s but within the range of previous variability. The weak warming signal in Bakhtiyorov et al. (2025) seems to be supported by the instrumental record from the Kluhorski Pereval station (2037 m a.s.l.) in the Greater Caucasus and the TRW based spring-temperature reconstruction from western Anatolia (Köse et al., 2017). In contrast, the gridded temperature product CRU TS 4.09 and the nearby high-elevation station records (Ardahan; 1830 m a.s.l.) analysed in this study reveal stronger warming rates (Fig. S1). In addition, Bakhtiyorov et al. (2025) use trees from sites situated below 1900 m a.s.l. with a weaker climate/growth relationship questioning the reliability of the reconstruction. A multispecies conifer LWBI network from upper-treeline sites in the Lesser and Greater Caucasus showed strong summer temperature sensitivity but also did not fully capture the recent warming (Dhyani et al., 2025). Dhyani et al. (2025) conclude that the weaker expression of recent warming in available LWBI records from the region likely indicates methodological limitations of LWBI rather than actual temperature differences. However, despite these discrepancies, our LC reconstruction and the LWBI-based temperature record developed by Dolgova (2016) from the Greater Caucasus show highest calibration skill and strong correlation between each other which supports a reliable representation of regional temperature variability in reconstructions with high signal-to-noise ratios.

4.3 Caucasian temperature history and continental scale climate dynamics

Large-scale comparisons place our LC reconstruction within Eurasian climate dynamics (Fig. 5). The climatic setting of the Caucasus is a result of the combined influence of moist westerly circulation, subtropical high-pressure systems and strong continental conditions, which create a transitional regime between European, northern African and Asian climate systems. The spatial correlation fields of our MXD reconstruction represent this pattern with strong correlations extending eastward into the continental drylands of Kazakhstan and southward into the subtropical deserts of the Arabian Peninsula, while correlations decline more rapidly towards the west (Fig. S5). A rapid correlation decline towards the west is further evident in the spatial correlation pattern with N-TREND temperature fields (Anchukaitis et al., 2017) and direct comparisons with MXD based temperature reconstructions along the Alpide Belt (Fig. 5). Our LC record shows strong coherence with temperature reconstructions from the Greater Caucasus (Dolgova, 2016) and a moderate correlation with the Eastern Mediterranean region (Esper et al., 2020), whereas correlations with reconstructions from the Alps (Kuhl et al., 2024) and the Pyrenees (Dorado Liñán et al., 2015) are weak or absent. Particularly, high-frequency correlations decline westward, shifting from weak positive values with the Eastern Mediterranean region to negative correlations with the Alps (Fig. S6). The correlation patterns indicate that temperature variability in the Caucasus partly opposes Central European temperature variability emphasizing the need to capture such regional characteristics. In addition to these spatial differences, the volcanic cooling signals presented show strongest post volcanic cooling response in the LC (1.1 °C) and Greater Caucasus (0.83 °C) whereas the Eastern Mediterranean, Pyrenees and Alps display less pronounced responses (Fig. 5). Similar to regional correlation patterns, the current record shows a stronger and more coherent volcanic response in the continental interior, while the impact becomes progressively weaker towards the maritime western regions likely due to enhanced oceanic buffering and circulation-related moderation of post-eruption cooling (Brönnimann et al., 2019).

Late-20th century warming as well as cooling associated with major volcanic eruptions are evident in our LC reconstruction and in the other regional records, although the timing and magnitude of these characteristics vary among sites. Our LC reconstruction shows only a weak expression of the multi-decadal cooling in the early 19th century, in contrast to its clear prominence in most Alpine records (Esper et al., 2025). The strong agreement with the Greater Caucasus summer temperature reconstruction of Dolgova (2016), which also shows weak multidecadal cooling during the early 19th century, confirms this signal. Lower values for the low-to-high-frequency variability ratios support the interpretation of stronger continentality in the Caucasus region (Fig. S7).

If the exceptional warming observed in the present study during the 21st century continues, as also anticipated by climate model projections (Shahgedanova et al., 2009), it will likely further impact human societies across the Caucasus region. The reconstruction shows that no previous warm interval in the last three centuries approaches the magnitude or persistence of the current episode, emphasizing that modern warming is without analogue in the proxy record. Snow accumulation, which serves as a crucial natural reservoir sustaining water supply during warm months is also projected to decline with continued warming in the Caucasus (Shahgedanova et al., 2009; Stokes et al., 2007; Tielidze et al., 2025) further intensifying water stress as human demand continues to grow (Mankin et al., 2015).

5 Conclusion

This study presents the first MXD-based summer temperature reconstruction from the Caucasus, extending back more than three centuries. By applying X-ray µCT scanning, we overcame limitations of traditional densitometry and demonstrated the strong sensitivity of tree-line Pinus sylvestris MXD to summer temperatures. The skill and the high proportion of variance explained by JJAS temperatures exceed the statistics for most of the single-site MXD records from the Northern Hemisphere. Our reconstruction captures major volcanic cooling peaks as well as the onset of sustained modern warming in the early 21st century, showing that recent summer warming is unprecedented over at least the past three centuries. More broadly, this record closes a critical spatial gap in mid-latitude tree-ring based temperature reconstructions and highlights the potential of MXD proxies to provide robust, spectrally unbiased estimates of past climate. In addition to its regional relevance, the chronology contributes to large-scale initiatives such as N-TREND, helping to better constrain Northern Hemisphere summer temperature variability. Continued development and temporal extension of MXD datasets in the topographically complex and climatically variable Caucasus region will be crucial for refining our understanding of past climate dynamics and for placing ongoing anthropogenic warming into a long-term perspective.

A primary challenge will be the extension of the Caucasus tree-ring reconstructions beyond the current few centuries to many centuries or even millennia to investigate the full period of the LIA as well as the earlier Medieval Climate Anomaly. The relatively short lifespan of trees in the region restricts the temporal reach of reconstructions based on living wood only. Future efforts must therefore integrate archaeological material and wood from old buildings to provide insights into climate variability in the deeper past and to understand human and landscape evolution in the context of climate change for the Caucasus region.