Articles | Volume 18, issue 7
Research article
13 Jul 2022
Research article |  | 13 Jul 2022

Investigating stable oxygen and carbon isotopic variability in speleothem records over the last millennium using multiple isotope-enabled climate models

Janica C. Bühler, Josefine Axelsson, Franziska A. Lechleitner, Jens Fohlmeister, Allegra N. LeGrande, Madhavan Midhun, Jesper Sjolte, Martin Werner, Kei Yoshimura, and Kira Rehfeld

The incorporation of water isotopologues into the hydrology of general circulation models (GCMs) facilitates the comparison between modeled and measured proxy data in paleoclimate archives. However, the variability and drivers of measured and modeled water isotopologues, as well as the diversity of their representation in different models, are not well constrained. Improving our understanding of this variability in past and present climates will help to better constrain future climate change projections and decrease their range of uncertainty. Speleothems are a precisely datable terrestrial paleoclimate archives and provide well-preserved (semi-)continuous multivariate isotope time series in the lower latitudes and mid-latitudes and are therefore well suited to assess climate and isotope variability on decadal and longer timescales. However, the relationships of speleothem oxygen and carbon isotopes to climate variables are influenced by site-specific parameters, and their comparison to GCMs is not always straightforward.

Here we compare speleothem oxygen and carbon isotopic signatures from the Speleothem Isotopes Synthesis and Analysis database version 2 (SISALv2) to the output of five different water-isotope-enabled GCMs (ECHAM5-wiso, GISS-E2-R, iCESM, iHadCM3, and isoGSM) over the last millennium (850–1850 CE). We systematically evaluate differences and commonalities between the standardized model simulation outputs. The goal is to distinguish climatic drivers of variability for modeled isotopes and compare them to those of measured isotopes.

We find strong regional differences in the oxygen isotope signatures between models that can partly be attributed to differences in modeled surface temperature. At low latitudes, precipitation amount is the dominant driver for stable water isotope variability; however, at cave locations the agreement between modeled temperature variability is higher than for precipitation variability. While modeled isotopic signatures at cave locations exhibited extreme events coinciding with changes in volcanic and solar forcing, such fingerprints are not apparent in the speleothem isotopes. This may be attributed to the lower temporal resolution of speleothem records compared to the events that are to be detected. Using spectral analysis, we can show that all models underestimate decadal and longer variability compared to speleothems (albeit to varying extents).

We found that no model excels in all analyzed comparisons, although some perform better than the others in either mean or variability. Therefore, we advise a multi-model approach whenever comparing proxy data to modeled data. Considering karst and cave internal processes, e.g., through isotope-enabled karst models, may alter the variability in speleothem isotopes and play an important role in determining the most appropriate model. By exploring new ways of analyzing the relationship between the oxygen and carbon isotopes, their variability, and co-variability across timescales, we provide methods that may serve as a baseline for future studies with different models using, e.g., different isotopes, different climate archives, or different time periods.

1 Introduction

Under the current anthropogenic warming trend (Shukla et al.2019), the interest in understanding its impacts on the mean temperature and precipitation and changes in their variability increases. Evaluating the representation of the mean state and the variability of past climate as simulated by climate models is crucial for reliable future projections (Schmidt et al.2012). The abundance of the heavy oxygen isotope 18O to 16O in stable water isotopologues (SWIs), further denoted as δ18O, is a proxy and tracer of natural variability in the water cycle. Oxygen isotope compositions can be measured in many paleoclimate proxy archives such as trees, ice cores, corals, or marine and lake sediments, which collectively extend our knowledge of climatic changes beyond the instrumental record (Bradley1999). Speleothems are secondary cave carbonate deposits, which form in karst systems globally, most commonly in the low latitudes to mid-latitudes, under a wide range of climate conditions (Fairchild and Baker2012). They provide precisely and absolutely dated (semi-)continuous time series of proxy data (Comas-Bru et al.2019) and have long been used as archives for terrestrial climate (Hendy1971). Oxygen isotopes are incorporated in calcite or aragonite matrices in accumulated growth layers. The stable carbon isotopic ratio (δ13C) is a second useful proxy for climate and vegetation conditions (Wong and Breecker2015).

Broad correspondence between speleothem δ18O and surface temperature (e.g., McDermott et al.2001) or local rainfall strength and seasonality (e.g., Medina-Elizalde et al.2016; Kennett et al.2012; Cheng et al.2016) and between speleothem δ13C and vegetation cover can be resolved in regional and global analyses (Comas-Bru et al.2019; Fohlmeister et al.2020; Baker et al.2019; Lechleitner et al.2021). Modification of these signatures by vadose zone fractionation processes (Tremaine et al.2011; Grossman and Ku1986; Romanek et al.1992), karst hydrology, internal cave conditions (Fairchild and Baker2012; Wackerbarth et al.2010; Jean-Baptiste et al.2019; Fohlmeister et al.2020), and differences in geochronological methods between records can complicate paleoclimatic interpretations (Breitenbach et al.2012; Rehfeld and Kurths2014). Age–model standardization (Comas-Bru et al.2020b), multiproxy approaches (Tremaine and Froelich2013; Warken et al.2018), and cave microclimate and drip water chemistry monitoring (Baker et al.2014; Treble et al.2015; Duan et al.2016), however, allow for statistically robust time series comparisons and can substantially improve our ability to disentangle climatic influences from site-specific processes across disparate climate zones (Fohlmeister et al.2017).

Depending on the specific site, some proxies may be easier to interpret than others. As such, carbon isotopes sometimes need to be pre-constrained through the help of other proxies, e.g., δ18O, to determine dominant processes (Fohlmeister et al.2017). However, speleothem oxygen isotopes can carry a less straightforward signal than carbon isotopes, where overlapping processes in specific regions complicate the δ18O interpretation (Scholz et al.2012; Ridley et al.2015), especially during large climate changes such as the deglaciation (Genty et al.2006). Studies considering both isotopes profited from the two proxies' mutual information on fractionation processes and were able to disentangle the encoded climatic signal, e.g., between speleothem δ18O for temperature and speleothem δ13C for vegetation changes (Fohlmeister et al.2017; Baker et al.2017; Novello et al.2019, 2021; Voarintsoa et al.2017).

The climatic interpretation of isotopes in speleothem records, as well as in other paleoclimate archives, is not always straightforward and provides only incomplete information and constraints on the dynamics of the Earth's climate system. Climate models with incorporated SWIs are another tool to better understand the climate system and are particularly powerful when applied in tandem with geochemical proxy records (Comas-Bru et al.2019; PAGESHydro2k-Consortium2017). Incorporating SWIs within the Earth's hydrological cycle in atmospheric general circulation models (AGCMs), atmosphere–ocean general circulation models (AOGCMs), and the most complex Earth system models (ESMs), is usually done by adding an additional water cycle to the hydrology of the model under explicit consideration of isotopic fractionation processes through H2O phase transitions (e.g., Tindall et al.2009; Yoshimura et al.2008; Werner et al.2016; Brady et al.2019; Lewis and Legrande2015). This opens the possibility to study and analyze the complete information of the modeled climate system consistent with model physics in past and present climates. Comparing modeled climate to the archived isotopic signatures provides an “equal ground” comparison opportunity (e.g., Werner2010; Sturm et al.2010; Xi2014; PAGESHydro2k-Consortium2017).

Extensive multi-model comparisons exist for past, present, and future as the Paleoclimate Model Intercomparison Project (PMIP, particularly for Jungclaus et al.2010; Kageyama et al.2018) under the overarching Climate Model Intercomparison Project (CMIP, particularly for Taylor et al.2012; Eyring et al.2016) to better understand the causes of model spreads in future projections and especially in temperature and precipitation. Simulations of the historical period (1850–2014 CE, as in Eyring et al.2016) or the last millennium (850–1850 CE, as in Eyring et al.2016) that will be the focus of this study, as well as idealized experiments under a range of natural and external forcings, are evaluated using different variables. SWIs, however, have not been included in the CMIP5/CMIP6 assessments (Taylor et al.2012; Eyring et al.2016).

The Stable Water Isotope Intercomparison Group (SWING) formed to bridge this gap of systematic intercomparison between isotope-enabled model simulations and compares them to observations over the historical period, while at the same time providing a large dataset to the scientific community (Risi et al.2012). The historical period is, however, too short to analyze and compare multi-decadal to centennial isotopic variability. To evaluate paleoclimate model simulations while also investigating changes in climate, the combination of SWI from proxy records and model simulations is essential to improve our understanding of climate change and its variability (Phipps et al.2013).

The Speleothem Isotope Synthesis and Analyses (SISAL) working group has collected a large number of speleothem records globally and compiled the database SISALv2 (Atsawawaranunt et al.2018; Comas-Bru et al.2019; Comas-Bru et al.2020b). It has been employed for model–data comparisons of the Last Glacial Maximum, the mid-Holocene, the last millennium, and the historical period using different models (iCESM: Midhun et al.2021; iHadCM3: Bühler et al.2021; ECHAM5-wiso: Comas-Bru et al.2019; Parker et al.2021; and GISS-E1-R: Parker et al.2021). Previous model–data comparisons support the use of the database to evaluate modeled δ18O across different time periods (Comas-Bru et al.2019), although speleothems exhibit a lower δ18O variability than simulated δ18O on interannual to decadal timescales globally (Bühler et al.2021). However, a benchmarking study of model performance in simulating δ18O and its variability, including a multi-model comparison or a model–data comparison with SISALv2, has not yet been performed.

Variability in models can either result from internal interactions and processes within the model or from changes in external radiative forcings (e.g., greenhouse gases (GHGs), volcanoes, and solar irradiance as in Fig. 1). Variability in the speleothem isotopic signal can also be a consequence of climate-related variability in oxygen isotopes as reflected in climate modes (e.g., El-Niño–Southern Oscillation (ENSO), Sun et al.2018; Midhun et al.2021; the North Atlantic Oscillation (NAO), Scholz et al.2012; or the Indian summer monsoon, Fleitmann et al.2007; Neff et al.2001) or changes in radiative forcing. Variations in δ18O and δ13C are routinely attributed to changes in solar radiation as a consequence of its influence on climate modes of variability, temperature, or precipitation (Warken et al.2021; Lone et al.2014; Cosford et al.2008; Neff et al.2001). While modeled variability commonly underestimates measured variability in paleoclimate archives, with increasing discrepancies on longer timescales (Laepple and Huybers2014a), cave-internal variability in speleothems may also overlay the archived signal. Especially on subdecadal to decadal timescales, lag time between the surface rainfall and the cave drip water, as well as the usually slow response of the cave microclimate to the surface climate, dampens the signal.

Here we will present a multi-model comparison of five isotope-enabled last millennium simulations: ECHAM5/MPI-OM (Sjolte et al.2018), GISS ModelE2-R (Lewis and Legrande2015; Colose et al.2016a, b), the isotope-enabled version of the Community Earth System Model (Stevenson et al.2019; Brady et al.2019), the isotope-enabled version 3 of the Hadley Model (Bühler et al.2021), and the stable water-isotope-incorporated Scripps Experimental Climate Prediction Center's Global Spectral Model (Yoshimura et al.2008), with climate characteristics and forcings as depicted in Fig. 1 and listed in Table 1. This multi-model comparison complements previous work (Jungclaus et al.2017; Midhun and Ramesh2016; Conroy et al.2013) through its focus on how different models represent SWI and their variability on different timescales over the entire last millennium. We aim to identify common model biases (Kageyama et al.2018) globally and in different regions, as well as distinguish specific climate drivers for modeled isotope variability on decadal and longer timescales. This allows, for the first time, for the joint intercomparison of SWI variability in climate models and proxy archives in a time period dominated by natural forcing.

The last millennium is a suitable time period for model–data comparisons, as it provides an opportunity to study variability on decadal and longer timescales and to decipher internal from externally forced variability (Kageyama et al.2018). Boundary conditions such as orbital forcing, sea level, and ice sheets are close to present-day levels, and external variability is mostly driven by variations in volcanic eruptions (Schurer et al.2014; Neukom et al.2019; Legrande and Anchukaitis2015). It is a key paleoclimate period for the CMIP5 and CMIP6 experiments (Taylor et al.2012; Eyring et al.2016), and speleothem records are abundant in this period (Bühler et al.2021).

With this study, we aim to contribute to the understanding of both model and data by answering the following questions: (a) first, how do different simulations model oxygen isotopes in the hydrological cycle, and how do they compare to archived speleothem data? (b) Second, what processes influence speleothem isotope composition, and what effects of variability can be captured and later analyzed? We first compare their similarities and differences in the isotopic signatures of precipitation both globally (Sect. 4.1) and specifically at the cave site locations of a large number of high-resolution speleothems from the SISALv2 database (Comas-Bru et al.2020b). Through spatial testing of climate variables, we analyze the relationship between the measured stable isotopes of oxygen and carbon and different modeled climate variables (Sect. 4.2). These relationships will help to determine spatial biases between the models and drivers for modeled and recorded isotopes. Temporal coherence between models and data give insight into internal and externally driven climate variability on different timescales. In a second step, we thus investigate if models can reproduce variability as measured in speleothems on annual to centennial timescales (Sect. 4.3). Finally, we test the timescales of externally forced events, e.g., volcanic eruptions or variations in the solar irradiance, that speleothems archives are able to resolve (Sect. 4.4).

2 Data

In this study, we collected and standardized the output from five different isotope-enabled climate model simulations over the last millennium, as well as oxygen and carbon isotopes in speleothems from the SISALv2 database.

2.1 Isotope-enabled general circulation models

A major advantage given by modeling SWI is its ability to both temporally and spatially resolve the variability of isotopes in precipitation by adding H218O and HDO to the part which already simulates and traces the most abundant stable water isotope, H216O. Simulated δ18O will further be denoted as δ18Osim. In the atmospheric advection scheme of the model, which is generally a part of the model's dynamical core, all three water isotopologues behave identically. In case of a phase transition (such as melting, condensing, evaporating, and freezing), additional fractionation effects are applied to the two less abundant water isotopologues. These phase changes typically occur in evaporation from land and ocean surfaces; condensation during formation of clouds, rain, snow; and re-evaporation of precipitation below cloud (for example, Werner2010; Sturm et al.2010).

The models used in this study range from AGCMs forced with sea surface temperatures and sea ice distribution to atmosphere–ocean general circulation models (AOGCMs). Their basic characteristics and boundary conditions are listed in Table 1. They are both used individually in the analysis and by the ensemble mean of all models. Figure 1 shows the climate as represented by the different models and external forcings used in the simulations. Since SWING2, there has not been a consistent protocol for simulations with isotope-enabled models. Hence, the simulations used in this study largely follow the PMIP version 3 last millennium experiment protocol (Schmidt et al.2011, 2012) with its proposed climate forcing reconstructions, with some variations in vegetation and orography. Of the external forcings used, differences in the volcanic forcing may have the largest influence on differences between the simulations (Colose et al.2016a; Schmidt et al.2011), as different responses on larger eruptions may have a long-term impact. Large eruptions can cause local anomalies to the mean state of δ18O of up to ±1.5 ‰ (Colose et al.2016a), hinting at the magnitude of change that can be caused by different forcings. Volcanic eruptions are among the most prominent drivers of natural climate variability during the last millennium (Jungclaus et al.2017). Compared to volcanic forcing, the choice of solar or orbital forcing has a smaller effect over time in the last millennium. Although the simulations do use forcings based on different reconstructions, which then act on various timescales, differences in response may arise not only from forcings but also from their implementation in the models (Jungclaus et al.2017).

Figure 1Climate as represented by the different models (ECHAM5-wiso: light blue; GISS-E2-R: dark blue; iCESM: grey; iHadCM3: green; and isoGSM: orange) and external forcings over the last millennium. (a) Global mean surface temperature anomaly as represented by the models (in model colors) relative to the period of the last millennium (850–1850 CE), the reconstructed temperature anomaly (PAGES2k, red, PAGES2k-Consortium2019), and observed temperatures (HadCRUT4, black, Morice et al.2012). (b) Isotopic composition of precipitation in the different models (in model colors) at the cave site location of Bunker cave (Germany), including the δ18Ospeleo of entity ID 240 (red), all at the temporal resolution of entity ID 240 (Comas-Bru et al.2020b; Fohlmeister et al.2012). For both temperature and δ18O, we show the annual and downsampled values in lighter colors in the background. Bold colors show the values with a 100-year Gaussian kernel bandpass (Rehfeld and Kurths2014). (c) Atmospheric CO2 concentration (SMT: Schmidt et al.2012; MFM: MacFarling Meure et al.2006). (d) Volcanic forcing in units of aerosol optical depth (AOD) (CRO: Crowley et al.2008; Crowley and Unterman2013; GAO: Gao et al.2008), where the AOD for the Gao et al. (2008) reconstruction was estimated by dividing the sulfate loading by 150 Tg (following Atwood et al.2016). (e) Total solar irradiance (TSI) (STH: Steinhilber et al.2009, MSL: Muscheler et al.2016, VR: Vieira et al.2011). The grey bars mark particular periods of high volcanic forcing.

Sjolte et al. (2018)Lewis and Legrande (2015)Brady et al. (2019)Bühler et al. (2021)Yoshimura et al. (2008)Werner et al. (2016)Colose et al. (2016a, b)Stevenson et al. (2019)Tindall et al. (2009)(Gesch et al.1999)20091993Amante and Eakins2009Berger and Loutre (1991)Berger (1978)1988(Schmidt et al.2011)2006Schurer et al.2014Marland et al. (2003)(Schmidt et al.2011)Schmidt et al.2012(Schmidt et al.2012)Paul et al. (1998)Pongratz et al. (2008)Pongratz et al. (2008)Pongratz et al. (2008)Pongratz et al. (2008)Jungclaus et al. (2010)Hurtt et al. (2011)(Cox2001)Hurtt et al. (2011)Crowley et al. (2008)Crowley et al. (2008)Gao et al. (2008)Gao et al. (2008)2013Muscheler et al.Steinhilber et al. (2009)Vieira et al. (2011)Steinhilber et al. (2009)Vieira et al. (2011)20162007Wang et al. (2005)Wang et al. (2005)Schurer et al. (2014)Lean (2009)Schmidt et al. (2011)Schmidt et al. (2012)

Table 1Basic characterization of the last millennium simulation.

Download Print Version | Download XLSX


We use the isotope-enabled version of the fully coupled general circulation models (GCMs) ECHAM5–MPI-OM (hereafter, ECHAM5-wiso) (Werner et al.2016; Jungclaus et al.2006). The model consists of the atmospheric model ECHAM5 (Roeckner et al.2003) and the ocean model MPI-OM with an embedded sea ice model (Marsland et al.2002). The millennium-long simulation by Sjolte et al. (2018) covers the period 800–2000 CE and uses a similar setup to the E1 COSMOS ensemble by Jungclaus et al. (2010) but with a different solar forcing based on the solar modulation record inferred from combined neutron monitor and tree-ring 14C data (Muscheler et al.2016).

Isotope diagnostics have been implemented for the atmosphere, ocean, and land surface component of the model and are computed throughout the entire water cycle in the ECHAM5 (Werner et al.2016) and MPI-OM (Werner et al.2016) models. The land surface model assumes no fractionation in most of the physical processes (Haese et al.2013). Water tracers are fully mixed and advected in the ocean model, and the SWI total mass is conserved (Werner et al.2016).

ECHAM5-wiso has been used extensively within the paleoclimate field, as well as for modeling of the present climate (for example, Werner et al.2016; Langebroek et al.2011; Goursaud et al.2018). The fully coupled GCM version of the model ECHAM5–MPI-OM has good agreement with both present-day isotope observations from the GNIP database, as well as with ice core and speleothem proxies during mid-Holocene (Comas-Bru et al.2019), Last Glacial Maximum (Werner et al.2016; Comas-Bru et al.2019), and for the last interglacial (Parker et al.2021). Both in the GCM and with the atmospheric component (ECHAM5-wiso), a warm bias in the model is found over high-latitudinal regions, especially over Greenland and Antarctica (Werner et al.2011, 2016). This has been attributed to the coarse spatial resolution in the atmospheric component of the model (Werner et al.2016) resulting in an underestimation of isotope depletion in these regions. The last millennium simulation has not been evaluated globally, but climate reconstructions and variability of SWI have been compared in the North Atlantic region, where the amplitude of the variability was underestimated in the model compared to ice cores (Sjolte et al.2018). Previous studies also stress the isotopic response to volcanic eruptions and phases of NAO (Guðlaugsdóttir et al.2018, 2019).

2.1.2 GISS ModelE2-R

The isotope-enabled AOGCM GISS ModelE2-R (hereafter, GISS-E2-R) (Schmidt et al.2006, 2014) is used with the same physics as in the CMIP5 experiments (Miller et al.2014; Schmidt et al.2014). Water tracers and isotopes are incorporated into the atmosphere, land surface, ocean, and sea ice components of the model (Schmidt et al.2005). Several experiments have been set up for the last millennium with GISS-E2-R due to uncertainties in past forcings and their effects, with different combinations of solar, volcanic, and land use (vegetation) forcings but with the same greenhouse gas and orbital change in every case (Colose et al.2016a, b; Lewis and Legrande2015).

GISS-E2-R has been shown to simulate modern observations of SWI well, except over Antarctica, in terms of changes in convection, clouds, and isotope kinetics (Schmidt et al.2005). For the last millennium, GISS-E2-R has also explored the isotopic responses to volcanic eruptions in South America (Colose et al.2016a), to volcanic forcing in relation to the position of the intertropical convergence zone (Colose et al.2016b), and to ENSO (Lewis and Legrande2015). The model has been shown to have a warm sea surface temperature bias and issues in sea ice concentration around Antarctica, related to the transport of the Antarctic Circumpolar Current. In the tropics, a warm bias is also found over land, together with cooler northern midlatitudes (Schmidt et al.2014).

2.1.3 ICESM1

We use the last millennium run of the isotope-enabled AOGCM iCESM1 version 1.2 model (hereafter, iCESM) (Hurrell et al.2013; Brady et al.2019; Midhun et al.2021), a fully forced simulation out of an eight-member ensemble of different external forcings. As the model is open source and publicly available, it is widely used in the scientific community, and simulations for past (Zhu et al.2017; Zhang et al.2012) and present climate exist (Otto-Bliesner et al.2016).

The model consists of the isotope‐enabled Community Atmosphere Model version 5.3 (iCAM5.3, isotope-enabled version based on Neale et al.2010), a land model CLM4 (Oleson et al.2010), a sea ice model, and an ocean component that is based on the isotope-enabled POP2 (Zhang et al.2017). Isotopes in the water cycle are represented as a new parallel hydrological cycle in all hydrological components in the atmosphere, ocean, land, and sea ice in the form of numerical water tracers and can be tracked in space and time.

The isotope-enabled version captures general global isotopic signatures well over ocean areas but shows small discrepancies across the land surface (Brady et al.2019). This effect has been explained by the model showing a slight negative isotopic bias due to overestimated modeled convection in mid-latitude oceans. Consequently, the transport of SWI mass poleward and landward has been deemed insufficient (Nusbaumer et al.2017). Footprints associated with major climatic modes such as ENSO and Pacific decadal oscillation (PDO) are found to also be well represented in isotopic signatures (Midhun et al.2021).

2.1.4 IHadCM3

We use the last millennium run from the fully coupled isotope-enabled AOGCM iHadCM3 (Bühler et al.2021). iHadCM3 has been widely used to simulate present (Dalaiden et al.2020) and future climate (Sime et al.2008; Tindall et al.2009; IPCC2013), as well as for past climates (Tindall et al.2010; Sime et al.2013; Holloway et al.2018). The model consists of several components: the atmosphere model HadAM3 (Pope et al.2000), the ocean model HadOM3 (Gordon et al.2000), a sea ice model (as described in Valdes et al.2017) and a dynamic land surface and vegetation model (Cox2001).

For the isotope-enabled version, SWI were added as two separate water species in the atmospheric model, and as tracers in the ocean model. Fixed isotope fractions are added to a fixed volume grid box of the ocean and experience changes due to evaporation, precipitation, and runoff through a virtual flux of SWI, altering the δ18Osim ratio in the top level of the ocean accordingly (Tindall et al.2009).

Compared to instrumental observations, the model represents sea surface temperature, ice sheet, and ocean heat content well (Gordon et al.2000). The freshwater hydrological cycle in the model shows only a slight overestimation in the local evaporation (Pardaens et al.2003). The model simulates the major isotopic fractionation effects defined by Dansgaard (1964) (e.g., the latitude effect, the amount effect, and the continental effect) appropriately compared to GNIP data (Zhang et al.2012). Additionally, a broad agreement in isotopic output with GNIP data in the general spatial distribution can be observed (Tindall et al.2009).

2.1.5 IsoGSM

IsoGSM is the isotope-enabled version of the AGCM Scripps Experimental Climate Prediction Center's (ECPC) Global Spectral Model (hereafter, isoGSM) (Yoshimura et al.2008). The model is based on the previous medium-range forecast model used at NCEP, making it well documented in its performance as an operational weather forecast model (Kanamitsu et al.2002; Caplan et al.1997).

The IsoGSM is a standalone atmospheric model. Here, it has been forced with sea surface temperatures and sea ice distributions from the CCSM4 last millennium simulation (Landrum et al.2013). Land surface processes are modeled through the NOAH model, but isotopic fractionation is not considered in these processes (Yoshimura et al.2008).

IsoGSM has shown to represent isotope and precipitation observations globally using a spectral nudging technique captured by the NCEP–DOE Reanalysis dataset (Yoshimura et al.2008). Its last millennium simulation has not been evaluated in previous studies, but isoGSM captures large-scale isotope and climate patterns in present times (Risi et al.2012). The model has been shown to also reproduce observed isotopic and precipitation variability well over the regions of Indian summer monsoon (Berkelhammer et al.2012a), western North America (Berkelhammer et al.2012b), and northwestern Scotland (Baker et al.2012). Recently, IsoGSM showed good consistency with speleothem oxygen isotopes from East Asia (Chiang et al.2020) and from South Asia (Kathayat et al.2021). IsoGSM tends to underestimate the depletion of δ18Osim in dry regions such as Antarctica (Yoshimura et al.2008). A decreasing summer temperature increases the precipitation δ18Osim and is caused by the moisture transport scheme of the model associated with areas of extremely dry conditions (Yoshimura et al.2008; Risi et al.2012).

2.2 SISALv2 database

In this study, we use speleothem data from the Speleothem Isotope Synthesis and Analysis version 2 (SISALv2) database (Atsawawaranunt et al.2018; Comas-Bru et al.2019; Comas-Bru et al.2020a, b). The database includes 691 speleothem records from 294 caves across the globe from all continents except Antarctica. Globally, the last millennium has abundant records with sufficient resolution and reasonable dating uncertainties (Bühler et al.2021). We filtered the database for records that cover at least a 600-year period within the last millennium (850–1850 CE) and exhibit at least 2 dates and 36 stable isotope measurements within the time period to guarantee a minimum resolution of 30 years. We obtain 89 records from 75 different sites for δ18Ospeleo, of which 58 records (65 %) from 50 sites also have δ13Cspeleo measurements (Fig. 2).

Figure 2Site locations of the SISALv2 database on a global karst map (brown shading by Williams and Ford2006). We only consider entities that cover a minimum 600-year period within the last millennium and that include at least 2 dates and 36 isotopic measurements in this period. Shown are all sites within the last millennium meeting these selection criteria that have both oxygen and carbon isotope measurements (red) and only oxygen isotope measurements (blue).

3 Methods

To compare climate simulation outputs and speleothem data, both need some preparation. Modeled output comes in regular monthly resolution, while time series of speleothem proxies are irregularly sampled. Different speleothem mineralogies, as well as different isotopic standards between modeled and recorded data, need to be accounted for before statistical similarity measures are applied.

3.1 Speleothem drip water conversion

The database includes calcite, aragonite, and mixed mineralogy speleothem records. Following Comas-Bru et al. (2020b), we only use pure calcite or aragonite speleothems in our analysis. Stable isotope ratios are usually given in δ notation, with stable oxygen isotopes given as δ18O=(18O16Osample18O16Ostandard-1)1000 ‰, with the standard denoting either the Vienna Standard Mean Ocean Water (V-SMOW, Dansgaard1964; Kendall and Caldwell1998) for water in its liquid and gaseous phases, or Vienna Pee Dee Belemnite (V-PDB, Coplen et al.1983) for CaCO3. The stable carbon isotope ratio is denoted as δ13C=(13C12Csample13C12Cstandard-1)1000 ‰ against the V-PDB (Coplen1995).

To be able to compare precipitation δ18Osim values to those measured in calcite or aragonite, the δ18Ospeleo is converted to its drip water equivalent (δ18Odweq) relative to the V-SMOW standard as in Comas-Bru et al. (2019). This conversion is temperature dependent to different extents for both minerals. Tremaine et al. (2011) provide an empirically based fractionation formula for speleothems of calcite mineralogy:

(1) δ 18 O dweq = δ 18 O calcite - ( ( 16.1 1000 T ) - 24.6 ) ,

where T is the temperature in K and δ18O are given in units of ‰.

Aragonite speleothems form under different conditions, e.g., higher Mg content of the dripping solution or slow drip rate (Fairchild and Baker2012), resulting in a different fractionation factor compared to calcite as described by Grossman and Ku (1986):

(2) δ 18 O dweq = δ 18 O aragonite - ( ( 18.34 1000 T ) - 31.954 ) .

For both conversions, the cave temperature at the time of the fractionation is needed. As these are often not available, especially in a paleoclimate setting, we use annual mean modeled surface temperatures as a surrogate for cave temperatures (Fairchild and Baker2012). For caves in very cold conditions, the annual mean surface temperature may underestimate the mean cave temperature by some degrees due to long-lasting snow packs. This underestimation, however, only corresponds to an overestimation of 1 ‰ in δ18Odweq and is within the range of the simulation ensemble. Additionally, cave-dependent time lags between the surface and the cave temperature are not accounted for, as they have a negligible effect on the time-averaged mean isotopic value. The conversion to drip water equivalent is done for each entity and each simulation individually, where we use the simulated annual mean surface temperature, downsampled to the record's resolution.

In a last step, the V-PDB values are converted to V-SMOW for direct comparison to the modeled SWI, using the conversion from Coplen et al. (1983):

(3) δ 18 O SMOW = 1.03092 δ 18 O PDB + 30.92 .

For carbon isotopes, different fractionation paths exist depending on mineralogy. Following Fohlmeister et al. (2020), we convert the aragonite δ13C values to corresponding calcite values using a fractionation offset of 1.9±0.3 ‰ (δ13Cc=δ13Ccalcite=δ13Carag-1.9±0.3). This offset accounts for the different enrichment factors of the two polymorphs as established in laboratory studies (Romanek et al.1992) and confirmed in a speleothem study (Fohlmeister et al.2018).

From the drip water conversion, we obtain a matrix for each of the isotopes with one row per measurement and six columns, where one represents the observations and the other five the modeled estimates.

3.2 Data processing

The simulation output from the five different climate models is available at monthly resolution but with differing time coverage. All simulation runs are cut to cover the time period from 850 to 1850 CE, an interval that is similar to PMIP's interval in the last millennium experiment. We focus on surface temperature, precipitation, precipitation δ18Osim, and evaporation. For the simulations lacking evaporation as a diagnostic (iCESM, isoGSM, and iHadCM3), we convert latent heat to potential evaporation and use these variables within the simulations. Outliers in the simulation are removed by comparing each modeled value in the 3D-output data matrix with its eight neighbors in time and space. If the value deviates from the mean of these eight values by more than five standard deviations, the value is set to NA. On average 0.001 % of values are set to NA through this method.

δ18Ospeleo forms from drip water that reaches the cave, which can be approximated by calculating the difference between precipitation amount and evaporation. When comparing modeled to speleothem isotopes it is therefore more realistic to weight the modeled δ18Osim by precipitation minus evaporation amount (infiltration-adjusted precipitation weighting, iw) to obtain annual δ18Osim values. The weighting automatically focuses on the local seasonal composition of SWI in rainfall, and prevents bias towards seasons with minimal infiltration due to high evaporation. Weighted δ18Osim values are commonly used in speleothem model–data comparison studies (Comas-Bru et al.2019; Baker et al.2019). The simulated weighed δ18O mean is calculated as follows:

(4) δ 18 O iw = 1 n δ 18 O k iw k 1 n iw k ,

where δ18Oiw is the monthly weighted annual composition of isotopes, δ18Ok refers to monthly simulated δ18O, and iwk is the corresponding monthly amount of iw, where n=12. As isotopic fractionation also occurs during evaporation from the soil, models where δ18Osim is also available for soil layers would be more realistic to compare to speleothem data. However, these were only available for a few simulations. Using δ18Osim of precipitation for the calculation of δ18Oiw therefore offered a more equal handling of the data while maintaining the large ensemble and enabled a better comparison of results.

Where simulation data are compared on a grid box level, we block average all simulations to the same spatial resolution as that of the ECHAM5-wiso run, which has the lowest spatial resolution. Data at the speleothem cave sites are extracted via bilinear interpolation as in Bühler et al. (2021). Here, a grid box of the size of the simulation's original resolution with the cave's location at its center is resampled from the original grid boxes it overlaps with. We note that this bilinear interpolation can, however, influence the temporal variability or the values of the response variables (Latombe et al.2018).

The simulated monthly temperature, precipitation, and evaporation at the speleothem cave locations are averaged to annual mean values. All speleothems in our last millennium subset, however, come as irregularly sampled time series with a median resolution of 5.19 years per entity (90 % CI: 4.13, 6.99). Considering only the speleothems with measurements of both isotopes yields a median resolution of 6.08 years (4.07, 7.85). To adjust the model and speleothem data, the simulated variables of all models are block averaged to the irregular temporal resolution of the individual speleothem. We also include speleothems in our analysis where no δ13C measurements, but only δ18O data are available. In direct comparisons between carbon and oxygen isotopes, we only consider those 58 speleothems that provide samples for both isotopes.

The relationship between δ18Ospeleo and simulated climate variables is determined following three different latitude bands to guarantee enough data points within each climate zone, the tropics, the subtropics, and the extratropics (Holden2012). The tropics are commonly defined as the zone between the Tropic of Cancer and the Tropic of Capricorn (23.44 S to 23.44 N), the subtropical region is defined as 23.44–35 N/S, and the extratropical region is defined as 35–90 N/S (with cave sites only extending to 66 N and 42 S).

Spatial testing between speleothem δ18Odweq, δ13Cc, simulated δ18Oiw, and meteorological variables is done by linear regression of the simulated millennium mean, downsampled to the temporal resolution of each record, and entity mean. To account for the spread between simulated variables and calculated δ18Odweq, the linear regression is done via bootstrapping (n=2000). Confidence intervals for all entity mean variables are also calculated via bootstrapping with a significance level of α=0.1. The p values are calculated through a fit linear regression model (fitlm.m, MATLAB2018) using Pearson’s product moment correlation.

Correlation estimates and p values for regular time series, i.e., the annual resolution output of the simulation, are calculated via the Pearson's product moment correlation (via the function cor.test, R Core Team2020). We use a significance level of α=0.1.

Correlation estimates for irregular time series are calculated via Pearson correlation as adapted by Rehfeld and Kurths (2014) and tested for last millennium speleothem records in Bühler et al. (2021). Here, we also choose a significance level of α=0.1. This level is appropriate for both the regular and irregular time series, considering the number of samples N compared to the strictness and expected level of false positives. Whenever calculating correlation estimates where speleothem data is involved, we use the raw δ18Ospeleo or δ13Cspeleo instead of the converted drip water values to decrease any potential biases.

With all simulated variables downsampled to the irregular resolution of each speleothem record, the use of power spectral analysis of the time series can describe the variation of common signals on a frequency spectrum of all time series (Chatfield2003). The power spectral density (PSD) of a time series describes the power distribution versus frequency over a finite interval of time (Chatfield2003).

To compute spectra of irregularly sampled time series, these are first interpolated to their mean resolution by a double interpolation and filtering process (following Laepple and Huybers2014a; Rehfeld et al.2018; Dolman et al.2021). This interpolation is performed to reduce high-frequency artifacts. The robustness of this spectral estimation process was recently confirmed by Hébert et al. (2021).

3.3 Synchronous extreme events in speleothem isotopes

An alternative similarity measure to correlation estimates, especially given the irregularity of the available time series, is checking for synchronous extreme events within two time series. After distinguishing extreme events, strength and direction of synchronous extreme events are only based on their relative timing (Rehfeld and Kurths2014). In this study, we only focus on the relative timing of the distinguished extreme events to each other.

Within the modeled or measured isotope time series, we distinguish the 5 % extreme values as the values above or below the 97.5 % or 2.5 % quantile of the time series distribution, respectively. Two extreme values occurring in a row are treated as two separate extreme events. Extreme values for time series of solar irradiance are determined in the same manner. For the volcanic forcing time series, we distinguish those events above the 95 % quantile of the distribution.

Two events are considered synchronous if they both occur within a time period around the events limited by a local threshold τ. This local threshold is calculated for each possible pair of extreme events and is chosen as half the minimum time between either extreme event and its preceding or succeeding extreme. For our dataset, the median τ is 4.62 years (90 % CI: 4.37, 5.28), which is of the same order of magnitude as the median resolution of all records with both carbon and oxygen measurements of 6.08 years, (4.07, 7.85). We set a hard threshold limit of 50 years, corresponding to the median age uncertainty considering the original chronologies and the SISALv2 ensemble chronologies (Comas-Bru et al.2020b). When comparing synchronous events between isotopes within one particular speleothem, age uncertainties are negligible in the comparison as δ13Cspeleo and δ18Ospeleo values stem from the same measurement of individual sub-samples.

For comparable extreme event signatures between the modeled and measured isotopes to volcanic or solar forcings, each modeled time series is checked for synchronicity against their respective forcing time series. The analysis is repeated for the speleothems for the same number of forcings and averaged.

When looking at the temporal distribution of global synchronous extreme events, they are sorted into bins of 10 years, which is approximately twice the median local τ. If a pair shows several extreme events within one bin, it is only counted once. We determine significance by randomly permuting one of the time series of a pair and repeating the analysis 2000 times. Within one bin, all counts that are larger than the 95 % quantile of this “mean background noise” can be regarded as significant.

4 Results

4.1 Overview of simulated versus measured speleothem oxygen isotopes

We first compare the mean δ18Oiw signal of the five different last millennium simulations to see potential model biases and large differences between the simulations (Figs. 3 and S1 in the Supplement). The global mean δ18Oiw values are fairly similar in area-weighted global mean of −8.48 ‰ (90 % CI: −8.61, −8.36) and −8.41 ‰ (−8.62, −8.2) for isoGSM and GISS-E2-R, respectively. The ECHAM5-wiso run is less depleted with a global δ18Oiw mean of −7.27 ‰ (−7.46, −7.09) and with visibly less depleted mid-latitude oceans than in the other simulations. The iCESM and iHadCM3 data show a stronger depletion of −9.39 ‰ (−9.51, −9.28) and −9.15 ‰ (−9.29, −9.01) respectively, with iCESM showing stronger depletion in the mid-latitudes and iHadCM3 showing stronger depletion towards the Antarctic compared to the other simulations. Although GISS-E2-R shows strong depletion, especially in the Arctic region, the less depleted mid-latitudes dominate the global mean (see Fig. S1).

Figure 3Simulated δ18Oiw climatology (a–e) of the respective simulation: (a) ECHAM5-wiso, (b) GISS-E2-R, (c) iCESM, (d) iHadCM3, and (e) isoGSM in the background. The time-averaged mean δ18Odweq values using the respective simulated temperatures are depicted as circles. Global means (GM) of δ18Oiw are given in the title of each subplot. (f) The range of δ18Oiw between all simulations for each grid box, as well as the range for the difference between simulation and record. Light colors indicate large agreement between the simulations, while darker colors mark areas where the models differ strongly and the spread between the δ18Oiw is larger. Antarctic δ18Oiw ranges are up to 40 ‰, highlighting the different model performance in this region (white area in f).

This general offset between the global mean δ18Oiw is also visible when comparing the spread of mean values on a grid box level (Fig. 3f), where isotopic signatures differ the most strongly in the Antarctic. Restricting the view to low latitudes and mid-latitudes, the largest model data differences are found in the Sahara, the Arabian Peninsula, the Indian subcontinent, and Siberia, where low humidity, high-precipitation amount, or high continentality are the driving local forces of precipitation δ18O. Also shown in Fig. 3f) is the spread between the offsets to the respective δ18Odweq at the cave locations. A close agreement between the models does, however, not necessarily indicate a good model–data match. It only indicates that the differences between the converted speleothem data are similar in the models. The highest agreement of 2.31 ‰ is obtained at Tham Duon Mai cave in Laos (siteID 293 at 20.75 N, 102.65 E, Wang et al.2019), while the strongest disagreement of 8.14 ‰ between simulations is at Huangye cave in China (siteID 17 at 33.58 N, 105.11 E, Tan et al.2011).

Figure 4Kernel density estimates of (a) the general distributions in simulated and speleothem δ18O at cave locations and (b) differences between simulations and the speleothem last millennium mean (Δδ18O=δ18Oiw-δ18Odweq). Dashed lines represent the median.


Analyzing the differences between simulations and speleothems can indicate how well the model data matches the proxy signal. Here we investigate the differences in δ18O between simulated last millennium mean (δ18Oiw) and speleothems (Δδ18O=δ18Oiw-δ18Odweq) on a global scale through general distributions, where disagreements between each model and speleothem data are shown as kernel density estimates (Fig. 4). The full datasets are acknowledged through the mean value, whereas the median values exclude skewed distributions and extremes. The speleothem dataset has a median general distribution of −6.21 ‰ globally. Of the simulations, ECHAM5-wiso has the closest distribution median with −6.82 ‰, followed by isoGSM (−7.72 ‰), iHadCM3 (−8.20 ‰), GISS-E2-R (−8.25 ‰), and iCESM (−9.79 ‰) (Fig. 4a). The distributions between simulations and speleothem δ18O (Fig. 4b) are fairly symmetrical but vary between the simulations, with medians of 0.72 ‰, −0.86 ‰, −2.01 ‰, −0.67 ‰, and 0.28 ‰, for ECHAM5-wiso, GISS-E2-R, iCESM, iHadCM3 and isoGSM, respectively. While ECHAM5-wiso has the closest median globally, isoGSM has the smallest difference between simulation and speleothem δ18O at cave locations, with a mean of −0.17 ‰ (90 % CI: −0.66, 0.33). GISS-E2-R and iCESM have more negative differences, with a mean of −1.02 ‰ (−1.41, −0.62) and −2.04 ‰ (−2.50, −1.60) respectively, followed by iHadCM3 −0.68 ‰ (−1.18, −0.18). ECHAM5-wiso is the only model that has a mean positive difference between simulation and speleothem data, with a mean of 0.63 ‰ (0.20, 1.05), in line with the less depleted global mean seen in Fig. 3.

The largest positive difference (less depleted δ18Oiw than δ18Odweq) is found at Huagapo cave in Peru (siteID 277 at 11.27 S, 75.79 W, Kanner et al.2013) and Minnetonka cave in the USA (siteID 200 at 42.09 N, 111.52 W, Lundeen et al.2013) for at least four model simulations. The highest negative difference (at least three models agree) are found at Hoq cave on Socotra Island, Yemen (siteID 253 at 12.59 N, 54.35 E, Van Rampelbergh et al.2013), Diva cave in Brazil (siteID 38 at 12.38 S, 41.57 W, Novello et al.2012), and Qunf cave in Oman (siteID 159 at 17.16 N, 54.30 E, Fleitmann et al.2007).

The general patterns of the isotope climatology are similar between the models (Fig. 3) but can differ at cave locations (Fig. S2). Larger differences in the modeled isotopic signatures appear particularly in regions where modeled temperature also spreads widely. On average, the differences to the speleothem records (Fig. 4) appear small and are consistent with the general difference in precipitation δ18Oiw.

Figure 5Correlations between SWI and modeled temperature (a) and precipitation (b) time series in each grid box. The background shows the average over all five simulation correlation estimates between annual δ18Oiw and simulated annual temperature per grid box (a) and for precipitation (b). Crosses indicate grid boxes where correlation estimates for four or more models have the same sign as the averaged estimate over all simulations. Symbols indicate the mean correlation of the simulated temperature (precipitation) of all model simulations to the recorded δ18Ospeleo at record resolution. Crossed circles mark sites where more than four models agree in the mean sign of the correlation to δ18Ospeleo. Black circles indicate the location of those speleothems in the last millennium subset that show no significant correlation to any model. In (c) the map shows positive correlation estimates where |ρ(δ18O,T)| is larger than |ρ(δ18O,P)| (and vice versa for negative values).

Differences between δ18Oiw signatures between the models may arise from different simulation drivers for the oxygen isotopes, e.g., temperature and precipitation, and can hint at different processes that govern the isotopic water cycle at a certain region within the simulations. The mean of the correlation to these main climatic drivers to δ18Oiw shows high agreement between the simulations (Fig. 5, individual correlation fields in Fig. S3). For the correlation to temperature (Fig. 5a), two main domains can be distinguished: there is mainly a positive correlation to temperature in the mid-latitudes to high latitudes and on the continents and negative correlation in the low-latitude ocean. Large-scale agreement between the simulations is, however, limited to the higher latitudes and the tropical ocean.

Two domains are also apparent in the correlation to precipitation (Fig. 5b), which are even more clearly separated than for temperature. We find areas of negative correlation to δ18Oiw in the low latitudes to mid-latitudes and areas of positive correlation only in the high latitudes. The agreement between the simulations, indicated by the crosses, is higher in the low- to mid-latitudes. Combining the information of dominant drivers of δ18Oiw (absolute higher correlation per grid box in Fig. 5c), we see clear domains of temperature dominance in the Southern Ocean, the higher northern latitudes and some continental regions, while precipitation is more dominant in the lower latitudes, the Antarctic, and Greenland.

The inter-model comparison shows more agreement in the correlation fields to temperature than to precipitation when focusing only on cave locations: the sign of correlation between δ18Oiw and simulated temperature agrees for three or more simulations at 60 % of locations and for four or more simulations at 26 % of locations. For precipitation, only 11 % of locations agree in sign for three or more simulations and only 1.1 % with agreement in four or more simulations. The more uniform temperature response to external forcing may increase the total number of significant correlation estimates and thus also the number of locations that agree in sign.

The model–data comparison shows more significant temporal correlation estimates between simulated temperature to δ18Ospeleo and also more sign agreement between these correlation estimates than to simulated precipitation. For precipitation, we find no cave, where more than four models agree in the sign to the mean correlation between modeled and recorded δ18O. For temperature, four speleothems from four different cave sites show a significant correlation of absolute strength |c|>0.15 for at least four simulations.

The data suggests that two main drivers for δ18O can be distinguished in specific regions – temperature is dominant in the high latitudes, while precipitation appears to be the main driver in the low latitudes. This is what we expected following the principles established by Dansgaard (1964).

4.2 Spatial testing for climatic and environmental effects on speleothem δ18O and δ13C

The δ18O in precipitation shows global signatures depending on latitude or altitude amongst other variables (Dansgaard1964). We assess this by looking at the relationships between speleothem δ18Odweq and δ13Cc with the site-specific variables of latitude and altitude.

The δ18Odweq decreases as more northern speleothems are considered (hemispherically: Fig. 6a; globally: R2=0.22, p<0.00, Table S1 in the Supplement). The large spread in the mid-latitude Northern Hemisphere is mostly due to the high number of speleothem records available in both Europe and China, implying a high longitudinal spread within the database. Figure 6b shows a global negative relationship of δ18Odweq to altitude, even though records get scarcer as the altitude increases. No clear pattern is visible for δ13Cc, and the dataset has a large mean spread (Fig. 6c). A statistically significant relationship between altitude and δ13Cc cannot be established. However, the spread in δ13Cc appears to increase with increasing altitude (Fig. 6d), although under decreasing data density.

Figure 6Speleothem δ18Odweq (a, b) and δ13Cc (c, d) against latitude (a, c) and altitude (b, d) as provided by the SISALv2 database. Linear regression lines are shown separately for both the Northern Hemisphere and Southern Hemisphere in panels (a) and (c) with their respective the R2 and p values. In panels (b) and (d), R2 and p values correspond to the respective linear regressions across all latitudes. Confidence bounds are 90 %.


We further explore the simulated meteorological variables and investigate spatial relationships between mean δ18Ospeleo and model variable ensemble mean (Fig. 7, Table S1). We find a significant relationship between δ18Oiw and speleothem δ18Odweq across all latitude bands (Fig. 7a–c), with the strongest correlation in the extratropics. Furthermore, we find a global dependency of δ18Odweq to mean simulated temperature and precipitation at the cave sites. For both temperature and precipitation, we find the strongest relationships to δ18Odweq in the subtropical regions. In all three regions, the relationship to temperature always exceeds that of precipitation. In Fig. 7d–f cave site altitude information is applied by color codes. This helps in the interpretation of possible overlapping effects of temperature and altitude on δ18O. A significant relationship to annual evaporation can only be distinguished in the extratropical regions (Fig. 7j–l).

Figure 7Simulated weighted δ18Oiw (a–c), temperature (d–e), precipitation (g–i), and evaporation (j–l) against speleothem δ18Odweq for model ensemble mean in the tropics, subtropics, and extratropics. The tropical region (23.44 S to 23.44 N) is shown in panels (a), (d), (g), and (j); the subtropical region (23.44–35 N/S) is shown in panels (b), (e), (h), and (k); and the extratropical region (35–90 N/S) is shown in panels (c), (f), (i), and (l). Altitude information is applied as shaded colors in panels (d)(l). We use δ18Oiw for all simulations. Note the semi-logarithmic axes for precipitation and evaporation.


We further compare the same meteorological variables to the speleothem δ13Cc data (Fig. 8, Table S1). A significant relationship is only found with temperature in the extratropical region (Fig. 8c), with increasing temperatures corresponding to more depleted δ13Cc. The δ13Cc is also found to be enriched with altitude (R2=0.23, p<0.02, Fig. S4) in the extratropics. This δ13C–altitude relationship is not found in the other latitudinal bands.

Figure 8Climate dependence of carbon isotope variability. Shown are simulated ensemble-mean temperature (a–c), precipitation (d–f), and evaporation (g–i) plotted against speleothem δ13Cc. Altitude information is applied as shaded colors in (d)(i). We used linear regressions in all plots; however, these appear curved in a semi-logarithmic plot as used for precipitation and evaporation.


The spatial testing shows globally strong relationships between δ18Odweq to environmental factors, in particular to altitude, temperature, and precipitation, which is in line with previous studies (for example, Comas-Bru et al.2019; Baker et al.2019). The spatial relationships between speleothem entity mean δ13Cc and meteorological variables from model ensemble mean (Fig. 8) only show clear relationships in the extratropical region and not on a global scale. This indicates more local influences as discussed by Fohlmeister et al. (2020).

4.3 Variability on different timescales

We compare the variance distribution in oxygen and carbon isotopes over all speleothems. This is a useful measure of how climatic and environmental factors influence the proxies to a different extent. Additionally, different simulations have different representations of variability across different timescales. This behavior can be explored by calculating power spectral densities (PSD) of the simulated and recorded isotopes averaged globally.

Figure 9Spectral ratios of isotopes in speleothem and simulations on different timescales as shown by the ratios or mean power spectral densities (PSDs). (a) Spectral ratio between speleothem isotopes (δ18O/δ13C). (b–c) Spectral ratio over all cave locations for δ18Ospeleo and δ18Oiw per simulation (model colors). In (b) we show the spectral ratios of δ18Ospeleo to δ18Oiw downsampled to the individual record resolution, and in (c) we show the ratios to δ18Oiw at annual resolution. The full spectra are shown in faded colors, and a smoothed spectrum is shown in black or the model colors. (d) Variance of detrended δ18Ospeleo (red) and δ13Cspeleo (black) as measured in speleothem records. The dashed line indicates the median of the distribution.


Figure 9a provides the spectral ratio of the two isotopes after detrending the irregular time series. A flat spectral ratio of ∼1 would indicate same levels of variability for both isotopes on all timescales. The spectral ratio here shows higher variability of δ13Cspeleo on all timescales; however, for periods smaller than 10 years, the variability of both isotopes is more similar. This is supported by the total variance of the isotopes over the complete period. The δ13Cspeleo shows a much higher total variance with a median of 0.46 ‰2 (0.38, 0.6) compared to δ18Ospeleo with a median variance of 0.11 ‰2 (0.08, 0.12).

Figure 9b shows the measured average PSD of δ18Ospeleo divided by the simulated average PSD at annual resolution, and Fig. 9c) by the average PSD of δ18Oiw at record resolution. A spectral ratio larger than 1 indicates higher variability at the timescale of the recorded δ18Ospeleo, whereas spectral ratios smaller than 1 indicate higher variability of the simulated δ18Osim. The spectral ratios between δ18Ospeleo and simulations at the cave locations at annual resolution show lower variability in δ18Ospeleo compared to all models on decadal and shorter timescales although to different extents (Fig. 9b). When considering the simulations downsampled to record resolution, the result is similar but there is much lower variability in δ18Ospeleo at decadal and shorter timescales (Fig. 9c). By downsampling, the simulated spectra lose power on the decadal and shorter timescales, which is then reflected in higher spectral ratios. On decadal to centennial timescales, however, δ18Ospeleo shows much higher variability than the modeled δ18Osim and is unaffected by the downsampling.

Variability of δ18Oiw is modeled differently in the simulations, as represented by the different levels of spectral power in the ratios. This difference is supported by the five area-weighted global variances of different magnitude and the simulated variance in annual and downsampled resolution as listed in Table 2. Small deviations between calculated variance and variance as area under the PSD can arise from the interpolation before the calculation of the spectra. The general order of isoGSM having the highest and iHadCM3 the lowest power on shorter frequencies remains throughout, as can also be seen in the table with small deviations to the order. The global variance, however, is only partly represented at the cave locations. While isoGSM shows the highest variance globally and at cave locations, iCESM is of medium variance globally but has the smallest variance at cave locations. For unweighted isotopic composition, the order of simulations changes (results not shown).

The analysis suggests that variability in the simulated δ18Oiw is represented differently in the simulations and that the order is not frequency dependent. The recorded δ18Ospeleo shows more variability on centennial frequencies and less variability on decadal and smaller frequencies than the simulated values (albeit differing extents depending on the simulation).

Table 2Area-weighted mean global δ18Osim variances as varglobal with 90 % intervals of the distribution are given per simulation. Row two and three give the δ18Osim variance at the cave locations in the annual resolution of the simulation and downsampled to the records resolution, respectively.

Download Print Version | Download XLSX

4.4 Analysis of extreme events

To examine if there are factors that influence both δ18Ospeleo and δ13Cspeleo simultaneously, we analyze the similarity of both signals. A total of 86 % of the speleothems show significant correlation between both isotopes (results not shown). A different test for similarity is provided by checking for synchronous extreme events in the time series.

Figure 10(a) The synchronous extreme events between δ18Ospeleo and δ13Cspeleo (red), (b) the synchronous events between the speleothem isotopes (oxygen purple and carbon green) and volcanic eruptions as reconstructed by Crowley and Unterman (2013) or Gao et al. (2008) (depicted in Fig. 1d), (c, d) the synchronous extreme events between simulated δ18O values at the cave locations of all simulations downsampled to record resolution and annual resolution, respectively. Where occurrence of synchronous extreme events is significant with α=0.05, the bars are shown in dark colors, where it is non-significant the bars are shown in transparent colors. The four light grey bars in the background of each plot show areas of high volcanic activity.


Figure 10 shows the temporal distribution of extreme events globally. In Fig. 10a we test for synchronous extreme events between the two isotopes by studying their relative timing to changes in volcanic forcing. Despite the high number of significantly correlated oxygen and carbon isotopes within one record, global patterns are not visible with a maximum of 10 % of speleothems showing synchronous extreme events between δ18Ospeleo and δ13Cspeleo at the same time over the last millennium. Synchronous events are also not higher in time periods of strong volcanic eruptions (indicated by grey bars).

By analyzing synchronous events between volcanic eruptions as reconstructed by Crowley and Unterman (2013) and Gao et al. (2008) (Fig. 1d) and δ18Ospeleo as in Fig. 10b, up to 20 % of speleothems exhibit extreme events at the same time as extreme volcanic eruptions. This share is higher than for the carbon isotopes, where up to 15 % exhibit extreme events synchronous to volcanic eruptions. Both isotopes also show pronounced peaks occurring with different extreme volcano events (indicated by grey background) but also for minor volcanic events.

To check if the speleothems used here can resolve global extreme events of short duration, we compare δ18Oiw at the cave locations of each simulation to the volcanic forcing of the simulations (see Table 1). While up to 50 % of δ18Oiw at cave locations and annual resolution exhibit an extreme event at the same time as an extreme volcanic eruption (Fig. 10d), this number largely decreases to less than half when the resolution of δ18Oiw is decreased to that of the speleothems (Fig. 10c). The number of pseudo-speleothems experiencing synchronous events to volcanic forcing (Fig. 10c) is more similar to that of the speleothems (Fig. 10b). Summarizing, we see no global temporal pattern of synchronous extreme events in both δ18Ospeleo and δ13Cspeleo (Fig. 10a). A lower temporal resolution strongly decreases the ability of the modeled archive to resolve global events like volcanic eruptions.

5 Discussion

5.1 Comparison of oxygen isotope variability in isotope-enabled models and speleothems during the last millennium

We found that the mean δ18Oiw fields show global differences of up to 2.12 ‰ between the models that could most likely be attributed to the global mean temperature differences of up to 1.8 K between the models. Similarly, most of the strong regional differences in δ18Oiw between models could be explained by regional differences in simulated temperature (Fig. S5). This can be expected following the temperature effect (Dansgaard1964) or the more recently observed changes in isotopic signatures under the current anthropogenic warming trend (Rozanski et al.1992). Temperature differences are, however, not the only explanation for differences in isotopic signatures between models, and individual model biases are also visible (Figs. S1 and S2). Specifically, we found less depletion of δ18Oiw in isoGSM over Antarctica as an artifact of its numerical scheme used for moisture transport and linked to extremely dry regions (Yoshimura et al.2008). The warm bias in high latitudes for ECHAM5-wiso results is known to underestimate isotopic depletion (Werner et al.2011, 2016). Overestimation of fractionation processes in iCESM during re-evaporation processes resulted in generally stronger depletion in δ18Osim (Brady et al.2019). The cool bias in northern mid-latitudes in the GISS-E2-R model as found in Schmidt et al. (2014) resulted in more depleted δ18Oiw in this region. Colder temperatures over Antarctica in iHadCM3 explained partly why isotopic signatures are a lot more depleted than in the other simulations. The iHadCM3 data, however, compares well to historical ice core data (Tindall et al.2009), suggesting that the colder Antarctic conditions modeled by iHadCM3 may be more consistent with reality than the multi-model mean. Even though this study mostly focused on a terrestrial mid-to-low-latitude archive, local differences of modeled isotopic signatures in the Antarctic may have an influence on isotopic representation in the general circulation of the models.

At the cave locations, the spread between the simulations yielded 4.51 ‰ (90 % CI: 3.96, 4.79) in median (Fig. 3), while the median disagreement between simulated and speleothem δ18Odweq was around −0.38 ‰ (−0.8, −0.23) (Fig. 4). This means that even though the simulations differed strongly in some regions, using multiple models can be sufficient to average out the disagreement at individual cave locations. The differences with the speleothems are in agreement with those found by Bühler et al. (2021), who compared the SISALv2 database to the iHadCM3 last millennium simulation. Median differences to the speleothem records were small for all models, where differences spatially around the globe (Fig. 3) reflected both the internal variability of models and global differences between the models. For example, ECHAM5-wiso showed the highest values for δ18Oiw in the global mean and also a more positive distribution than the other models. Differences in global mean δ18Oiw may arise from the weighting by the simulated evaporation and precipitation amounts, and while evaporation in ECHAM5-wiso is similar to the other models (Fig. S6), precipitation amounts are lower (Fig. S7). The lower oceanic precipitation amount can most likely be attributed to the coarser vertical resolution of ECHAM5-wiso (Hagemann et al.2006).

Our analysis suggests that a multi-model approach is advisable whenever one is comparing mean modeled values to measured data. Even though global mean δ18Oiw values may be comparable, local and regional temperature estimates, and therefore modeled δ18Oiw values, can vary strongly and deviate between models. Even though isoGSM displayed the lowest disagreement between δ18Oiw and δ18Ospeleo (Fig. 4), additional processes between meteoric water above the cave and drip water may again influence this difference. The multi-model offset comparison justified the use of the multi-model mean at cave locations in the following spatial analysis.

We found significant relationships between all considered modeled climatic variables (temperature, precipitation, and evaporation) and δ18Oiw with mean δ18Odweq (Fig. 7). In general, the temperature and altitude effects on δ18Ospeleo are visible in all latitude bands. Effects such as the precipitation amount relationship to δ18Ospeleo are more distinct in the tropical and subtropical latitude band, which is in line with the effects described by Dansgaard (1964). In the extratropics, the precipitation effect is weak and overlapped by the altitude effect. We restricted our analysis to temperature, precipitation, and evaporation as they were available for all models. Future analyses that include more model diagnostics across the multi-model ensemble may find relationships to different variables that will improve our understanding of model biases and obtain a more coherent picture of which variables drive δ18O in the model world. The relationships between speleothem mean δ13Cc and meteorological variables from model ensemble mean (Fig. 8) were less clear, which again points to a more local influence.

Global studies that have evaluated δ18Ospeleo isotopic signatures using climate models already exist (Bühler et al.2021; Comas-Bru et al.2019; Midhun et al.2021), where regions with shared climatic features showed stronger relationships. For example, Baker et al. (2019) focused on specific climate zones by comparing drip water measurements to precipitation δ18O measurements and was able to identify temperature zones for which mean measured δ18O or δ18Oiw was most similar to δ18Ospeleo. In our study, we found stronger relationships to climate variables in the individual latitude bands than on a global scale as in Bühler et al. (2021). Analyzing regions with high data density and similar climate patterns, we found even stronger temperature relationships (Figs. S8 and S9). Still, local particularities, such as large elevation differences over short distances, could not be resolved properly by the simulations and explained many of the stronger outliers, especially in the tropics.

Comparing a last century subset of the SISALv1 database by Fohlmeister et al. (2020) with our last millennium SISALv2 subset yielded similar results for precipitation and altitude relationships to δ13Cspeleo (Fig. S10). In contrast to our latitudinal approach, they compared δ13Cspeleo to temperature on a global scale. They neglected clusters of speleothems where factors other than temperature play an important role for the carbon isotope composition (e.g., high amounts of precipitation, known cave-specific particularities and processes, or temperatures close to the natural limit of vegetation, i.e., <5C). The remaining records with temperatures between ∼7 to 27 C showed a positive trend between δ13Cspeleo and temperature. This trend is in contrast to our observation based on clustering the records according to latitudinal bands. With this approach, we find no relationship between δ13Cspeleo and temperature for the tropics and subtropics, but a clear inverse relationship between modeled temperatures and δ13Cspeleo is observed for the extratropical records.

Higher cave site elevation coincided significantly with more depleted δ18Odweq, which is in line with the altitude effect (Dansgaard1964). For δ13Cspeleo, local studies exist that predict an increase in δ13Cspeleo with higher altitude (Johnston et al.2013). However, for the global last millennium subset of the SISALv2 database, more entities with carbon measurements in higher altitudes are needed to see a potential global relationship in addition to the relationship we find in the extratropical latitude bands.

5.2 Can models reproduce variability archived in speleothems?

For all simulations, temperature variability was the dominant driver in δ18Oiw at high latitudes and precipitation variability was dominant at low latitudes and in parts of the Antarctic (Figs. 5c and S5). However, local and regional climate dynamics, such as landward moisture transport and ice sheet changes, can mask and alter these relationships, as found for simulated isotopes in GISS-E2-R in a global study by LeGrande and Schmidt (2009). At the cave sites, model-internal regional variability and record age uncertainties substantially decreased correlation estimates. We observed that the sign of the correlation estimated between simulated temperature and δ18Oiw agreed at 60 % of cave locations for ≥3 simulations, and at 26 % for ≥4 simulations. For correlation estimates of precipitation, this was true only at 11 % or 1 % of locations, respectively. When compared to measured δ18Ospeleo, we found more significant temporal correlation estimates to modeled temperature than to modeled precipitation (Fig. 5). This could in part be explained by global temperature responses to, e.g., volcanic forcing being more uniform between model ensemble runs compared to precipitation responses, which depend strongly on regional particularities. Regions with high inter-model climate variable spread (Fig. S11) also coincided with regions of the fewest significant correlation estimates to simulated temperature and the least agreement between the climatic drivers of δ18Osim.

For variability specifically at the cave site locations, δ13Cspeleo appeared to be more variable than δ18Ospeleo on average on all timescales (Fig. 9d) with an increasingly higher variability compared to δ18Ospeleo towards centennial timescales and longer (Fig. 9a). Within 86 % of all speleothems where both isotopes are provided, δ18Ospeleo and δ13Cspeleo showed significant correlations. Jointly explained variance in the isotopic signal can point towards common climatic drivers. However, the amplified variance on long timescales in δ13Cspeleo could either hint at changes in the water cycle but also in land surface processes such as soil formation or vegetation composition and density. Considering more terrestrial archives and trace elements stored within the speleothems may help to better disentangle the climatic and environmental signals. On decadal and shorter timescales, where the meteoric and seasonal vegetation isotopic signal is mostly smoothed by the karst system, higher variability in δ13Cspeleo may result from the stronger isotopic fractionation for carbon compared to oxygen in precipitated calcite (Polag et al.2010; Hansen et al.2019).

Climate models reflected δ18Oiw variability at cave locations to different extents. A clear offset between the models could be found on all timescales. We found no relationship between spatial resolution of the model and the variability of isotopic composition of precipitation. The higher-resolution run of iCESM and the lower-resolution run of iHadCM3 seem to show similar variability globally, and the lower-resolution ECHAM5-wiso shows even higher variability at the cave locations. Due to the strong impact of temperature and precipitation on δ18Osim variability, we expect that this difference in isotopic variability also stems from the difference in the simulated climate. Further assessments using multivariate statistics are needed to firmly attribute the impact of climate on recorded variability of δ18O.

The lower temporal resolution of speleothem records largely explained the model–data mismatch on decadal and shorter timescales. Slow growth rates and limited sampling resolution lead to averaging effects that in turn lead to lower variability on decadal and shorter timescales. Simple karst-filters of a realistic transit time of ∼2.5 years (as in Fig. S12 and used in Bühler et al.2021; Midhun et al.2021; Dee et al.2015) showed that variations in models and speleothems on these shorter timescales are similar if accounted for. Expert knowledge of the local cave hydrology is, however, needed for a more detailed assessment on which model reflects δ18Oiw variability best on decadal and shorter timescales compared to speleothems and may still be restricted by karst and cave internal processes that effectively limit the sampling of climatic signals.

On decadal and longer timescales, models seemed to underestimate δ18Ospeleo variability, although to different extents, with isoGSM showing the highest variability and the smallest model–data mismatch. The model–data mismatch that we observed between speleothems and δ18Oiw starting at decadal timescales is in agreement with previous studies such as Laepple and Huybers (2014b). However, it is worth mentioning that speleothems may also be capable of enhancing climate-driven changes of δ18O and δ13C by cave-specific processes, resulting in higher variability on decadal and longer timescales. Under the assumption that variability on decadal and longer timescales is recorded correctly by speleothems, the iCESM model showed the strongest variability mismatch and isoGSM the smallest.

Cave locations were not necessarily reflective of mean annual δ18O variability globally – at least not in the model simulations. Simulations that showed generally high variability at cave locations at speleothem resolution also tended to be more variable globally and vice versa for low variability. However, simulations that were less variable at cave locations than others can still be more variable globally. In our case, this trend could likely be attributed to the bias of geographic locations of the cave sites as the models mostly show high variance in δ18Oiw in very dry regions and around the regions influenced by the Intertropical Convergence Zone.

5.3 Can external forcings be resolved by speleothems?

We found that 86 % of speleothems show a significant temporal correlation between speleothem oxygen and carbon isotopes, with 47 % even showing strong significant (anti-) correlations of |c|>0.5. High co-variability between both isotopes can either be caused by kinetic fractionation processes (Hendy1971) in the cave environment or may be externally forced. For example, Fohlmeister et al. (2017) studied a stalagmite in an arid region and found strong correlation between the isotopes. They attribute increased correlation to times of strong variations in cave-internal processes triggered by variations of external conditions. This simultaneity agrees with our findings that generally no extreme event in isotopes precedes the other, which can, however, also be attributed to low sampling resolution. More local cave monitoring studies are necessary to potentially exclude kinetic fractionation effects as the dominant driver.

Temporal coherence between model and data can only be expected in the frequency space. However, variability associated with external forcing, such as volcanic eruptions or changes in solar forcing, may be imprinted in both model and data in a congruent time history (PAGESHydro2k-Consortium2017). Changes in temperature or precipitation due to aerosol-forced cooling have been analyzed in a δ13C record as signs of volcanic signatures of speleothems (Ridley et al.2015). Growth rate changes (Baker et al.1995) or the measurement of trace elements such as sulfur (Frisia et al.2008) are other techniques to detect volcanic signals in speleothems, but they generally require up to sub-annual resolution records. In our global analysis of 58 δ13Cspeleo and 89 δ18Ospeleo records, we saw no significant increase in extreme events in the isotope records coinciding with major volcanic eruptions. The individual isotopes yielded more distinct signatures of volcanic eruptions, with up to 20 % (15 %) of speleothems recording synchronous extreme events in δ18Ospeleo (δ13Cspeleo) and a volcanic eruption, similar to δ18Oiw at record resolution. However, both stayed well below the possible simulated detection at cave locations under annual resolution of up to 50 %. The comparison to the synchronous events to the downsampled δ18Oiw showed that the ability to capture events such as volcanic eruptions strongly decreases with record resolution. The attribution of specific peaks in speleothem data to volcanic events needs caution because of age uncertainties and other possible explanations for the changes, e.g., human settlements close by (Baker et al.1995). Increasing the bin size to the average age uncertainty within the last millennium subset of the SISALv2 yielded the same results (Fig. S13).

Solar variations are another external forcing that is often invoked as an influence on the monsoon cycle and has been investigated using speleothem records (Neff et al.2001; Lone et al.2014; Cosford et al.2008). In these studies, it is also standard to use high-resolution speleothems with a lowest resolution of about 5 years for Lone et al. (2014). We repeated the analysis for Fig. 10 to analyze synchronous extreme events of δ18Ospeleo and δ13Cspeleo to the total solar irradiance input (Fig. S14). While the general effects of decreased detectability with decreased resolution are also visible, the overall detection of extreme solar irradiance was much weaker than for volcanic eruptions. This was not only true at cave locations but globally when comparing simulated surface climate variables to solar variations (results not shown). As solar variation on this timescale is mostly cyclical compared to random extreme volcanic eruptions (compare Fig. 1), these results were to be expected from the methods used and are not in contradiction to the literature.

The impact of different climatic backgrounds on the δ18O signal in speleothem records or paleoclimate simulations in time periods such as the Last Glacial Maximum and the Holocene have been studied extensively (e.g., Comas-Bru et al.2019; Tierney et al.2020; Parker et al.2021). Periods of documented warmer and colder periods within the last millennium are for example the Little Ice Age (1550–1850 CE) or the Medieval Climate Anomaly (850–1250 CE) (definitions from Jones et al.2001). We note that neither the global mean surface temperature nor the simulated δ18O or global δ18O as recorded by the speleothems showed significant changes to the mean state on a global scale within the described periods within the last millennium (results not shown). This is in line with Neukom et al. (2019), who found no global coherence of cold or warm periods over the Common Era, even though local changes are observable (e.g., McDermott et al.2001).

Summarizing, the comparison to modeled valued showed that cave locations in this study are in general suitable to detect δ18Oiw variations due to modeled climatic changes in response to changes in volcanic or solar forcing. Even though speleothems are highly resolved archives with small age uncertainties compared to other archives, the median resolution during the last millennium of 6.08 years (4.07, 7.85) was not enough to resolve changes in δ13Cspeleo or δ18Ospeleo due to potential solar or volcanically induced climatic changes. Karst mixing effects which dampen the signal may decrease the ability to detect these changes even further.

5.4 Limitations

A current weakness of this type of analysis is that we only compared δ18Osim of PMIP version 3 generation models to archived isotopes. However, parameter and tuning choices within one model, especially in the cloud and convection scheme, have a strong imprint on δ18Osim signatures (for example Nusbaumer et al.2017 for iCESM, Field et al.2014 for GISS-E2-R). To further systematically explore and constrain modeled δ18O, a multi-model ensemble under different model setups will be needed. The multi-model comparison in this study only included how different models represent the oxygen isotopes, as δ13C is not included in isotope-enabled models yet. For iCESM, a carbon cycle in the ocean exists (Jahn et al.2015), and great effort is put into the incorporation of a carbon cycle to isotope-enabled GCMs by the scientific community. However, this incorporation into GCMs cannot resolve the cave carbon cycle. Proxy system models of intermediate complexity that explicitly model soil carbon isotopes processes, temperature, moisture, and vegetation above the cave can be a more useful approach.

We also only looked at δ18O and δ13C as possible proxies for climatic changes. In contrast to, e.g., ice cores, speleothems do not directly record precipitation δ18O. Instead, both δ18Ospeleo and δ13Cspeleo undergo fractionation processes that are influenced by various cave-internal processes that may not be climate related (Lachniet2009; Fairchild and Baker2012; Hartmann and Baker2017; Dreybrodt and Scholz2011). Considering additional processes such as prior calcite precipitation (PCP) or other geochemical climate-related proxies can help to decipher the climatic signal from karst- or cave-internal processes (Kaufmann2003; Schwarcz et al.1976; Owen et al.2016; Tremaine and Froelich2013; Noronha et al.2014).

When comparing the speleothem isotopes to volcanic data, we note that there now are more recent volcanic reconstructions available that suggest a modification to the timing or magnitude of last millennium eruptions (Sigl et al.2015). Given the temporal resolution of the speleothems, changes in timing of volcanic events would impact the comparisons of model data to only those speleothems with high temporal resolution. As only the timing of the extreme event is compared, a change in magnitude would be irrelevant. As for the comparison with the simulated data, only the response of the model to any given eruption is important (Colose et al.2016b), as we did not compare the timing of extreme events in speleothems to those in the simulation.

Compared to the capability of speleothems to cover complete glacial–interglacial cycles, we instead looked at short timescales when focusing on the last millennium and what drives variability on decadal to centennial timescales. The last millennium is considered a relatively stable time period, and climatic changes might not be strong enough to be fully captured by speleothems. On longer timescales, speleothems may still be a good archive to capture these larger changes (Genty et al.2006).

6 Conclusion

We presented a multi-model comparison over five last millennium isotope-enabled simulations (ECHAM5-wiso, GISS-E2-R, iCESM, iHadCM3, and isoGSM) and compared their representation of isotopic signatures in mean and variability to paleoclimate data from a large speleothem database (a last millennium subset of SISALv2). We found that δ18Oiw differed substantially between models on a regional scale and at speleothem cave sites. This could in part be attributed to differences in simulated temperature, model biases in implementing SWI or topography, and cave- and site-specific controls on speleothem isotopes. Extreme model values that differ greatly from the rest can be compensated for by using the multi-model mean and thus reducing local spatial biases. The isoGSM simulation showed the lowest absolute mean difference to the speleothems at cave locations, while all other simulations show only slightly higher disagreements. Variability on decadal and longer timescales in the speleothems was higher than indicated by models and was also best represented by isoGSM; however, isoGSM still underestimated variability on these timescales. No relationship was found between the spatial resolution of the models and the variability of the isotopic composition of precipitation. In all models, temperature was driving δ18Oiw variability at high latitudes and precipitation at low latitudes. At cave site locations in particular, which are mostly located at low latitudes to mid-latitudes, models agreed more on temperature (rather than precipitation) being the driving factor of SWI variability. However, temperature signatures in climate models are generally more uniform than those of precipitation, as these depend heavily on how models parameterize convection and cloud dynamics (PAGESHydro2k-Consortium2017).

Dividing the global set into latitude bands, we were able to distinguish temperature and altitude relationships for both the oxygen and the carbon isotopes, as well as significant relationships for δ18Odweq to other simulated climate variables. While most records showed significant correlation between the two isotopes, using both isotopes to gain more information than just from one remained difficult. This is especially difficult for variations in solar and volcanic forcing as they are not imprinted in either the single isotope or in the pair on a global scale. Many archive limitations could, however, be attributed to the low resolution of the dataset compared to the processes expected to be resolved.

This joint intercomparison of stable water isotopologue variability in both models and speleothem data is the first dataset in a time period of natural forcing and allows for further future analysis by the scientific community. Our analysis encourages the use of multi-model means whenever possible, as already suggested by other studies (Colose et al.2016a). From the point of model evaluation, the incorporation of different archives with higher resolution (e.g., corals, trees, ice cores as in the Iso2k database; Konecky et al.2020) and with the help of improved proxy system models may provide further insight into why offsets between models can be so large regionally. From a speleothem perspective, within-cave and between-cave variability comparisons using both isotopes will help us to understand the recorded signal better and give higher confidence in the results. Future multi-model comparisons of isotope-enabled models for other time periods are required to further evaluate biases in models, and comparisons to δ18O archives of all kinds are required to gain a deeper understanding of the underlying concept as to what influences variability and co-variability of isotopes in speleothems.

Code and data availability

Code to reproduce figures and analyses is provided at (last access: 3 June 2022; Bühler and Axelsson2022). Model data with output at cave locations (ECHAM5-wiso, GISS-E2-R, iCESM, iHadCM3, isoGSM) at annual and record resolution (as csv files), as well as monthly fields of surface temperature, precipitation, isotopic composition of precipitation, and evaporation or latent heat for all simulations, are freely available on Zenodo at (Bühler et al.2022). The SISAL (Speleothem Isotopes Synthesis and AnaLysis Working Group) database version 2 (SISALv2) is publicly available through the University of Reading repository at (Comas-Bru et al.2020a). We use R for the data analysis (R Core Team2020;, last access: 1 July 2022). The main packages are tidyverse (Wickham et al.2019,, ncdf4 (Pierce2019,, last access: 1 July 2022), ggplot2 (Wickham2016,, last access: 1 July 2022), raster (Hijmans2020,, last access: 1 July 2022), zoo (Zeileis and Grothendieck2005,, and plyr (Wickham2011,; Wickham et al.2021,, last access: 1 July 2022). We use the nest R package (, last access: 1 July 2022, Rehfeld and Kurths2014,; Rehfeld et al.2011, and the PaleoSpec package (, last access: 1 July 2022; Kunz et al.2020,


The supplement related to this article is available online at:

Author contributions

JB, JA, FL, and KR designed this study. JB and JA wrote the paper and created the figures. FL and JF contributed with data interpretations. KR, JS, KY, MM, ALG, and MW contributed with model data that JB and JA collected and standardized. JB prepared the speleothem data. JB, JA, KR, FS, JF, JS, KY, ALG, MM, and MW contributed to the revisions of the manuscript. All authors approved of the final version of the paper.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This study includes data compiled by SISAL (Speleothem Isotopes Synthesis and Analysis), a working group of the Past Global Changes (PAGES) project. PAGES received support from the Swiss Academy of Sciences and the Chinese Academy of Sciences. The project is in part inspired by discussions at the SISAL 4th workshop: Exploiting the SISALv2 database for evaluating climate processes, Xi'an, China, 14–18 October 2019. Nils Weitzel, Elisa Ziegler, Beatrice Ellerhoff, Qiong Zhang, Nikita Kaushal, Natasha Sekhon, Valdir Felipe Novello, Jon Baker, Ny Riavo Voarintsoa, and Yuval Burstyn are acknowledged for helpful advice, comments, and discussion on text and figures. Parts of the computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre (NSC), partially funded by the Swedish Research Council through grant agreement no. 2018-05973.

Financial support

This research has been supported by the Swedish Research Council (Vetenskapsrådet; grant nos. 2013-06476 and 2017-04232); the Deutsche Forschungsgemeinschaft (grant nos. 316076679 and 395588486); the Bundesministerium für Bildung, Wissenschaft, Forschung und Technologie (grant no. 01LP1926C); and the Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung through the PalMod project (SNSF; grant no. P4P4P2_186693).

The article processing charges for this open-access publication were covered by Stockholm University.

Review statement

This paper was edited by Qiuzhen Yin and reviewed by two anonymous referees.


Amante, C. and Eakins, B. W.: ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis, NOAA National Centers for Environmental Information [data set],, 2009. a

Atsawawaranunt, K., Comas-Bru, L., Amirnezhad Mozhdehi, S., Deininger, M., Harrison, S. P., Baker, A., Boyd, M., Kaushal, N., Ahmad, S. M., Ait Brahim, Y., Arienzo, M., Bajo, P., Braun, K., Burstyn, Y., Chawchai, S., Duan, W., Hatvani, I. G., Hu, J., Kern, Z., Labuhn, I., Lachniet, M., Lechleitner, F. A., Lorrey, A., Pérez-Mejías, C., Pickering, R., Scroxton, N., and SISAL Working Group Members: The SISAL database: a global resource to document oxygen and carbon isotope records from speleothems, Earth Syst. Sci. Data, 10, 1687–1713,, 2018. a, b

Atwood, A., Wu, E., Frierson, D., Battisti, D., and Sachs, J.: Quantifying climate forcings and feedbacks over the last millennium in the CMIP5–PMIP3 models, J. Climate, 29, 1161–1178, 2016. a

Baker, A., Smart, P. L., Barnes, W., Edwards, R. L., and Farrant, A.: The Hekla 3 volcanic eruption recorded in a Scottish speleothem?, The Holocene, 5, 336–342, 1995. a, b

Baker, A., Bradley, C., Phipps, S. J., Fischer, M., Fairchild, I. J., Fuller, L., Spötl, C., and Azcurra, C.: Millennial-length forward models and pseudoproxies of stalagmite δ18O: an example from NW Scotland, Clim. Past, 8, 1153–1167,, 2012. a

Baker, A., Hartmann, A., Duan, W., Hankin, S., Comas-Bru, L., Cuthbert, M. O., Treble, P. C., Banner, J., Genty, D., Baldini, L. M., Bartolomé, M., Moreno, A., Pérez-Mejías, C., and Werner, M.: Global analysis reveals climatic controls on the oxygen isotope composition of cave drip water, Nat. Commun., 10, 2984,, 2019. a, b, c, d

Baker, A. J., Mattey, D. P., and Baldini, J. U.: Reconstructing modern stalagmite growth from cave monitoring, local meteorology, and experimental measurements of dripwater films, Earth Planet. Sc. Lett., 392, 239–249,, 2014. a

Baker, J. L., Lachniet, M. S., Chervyatsova, O., Asmerom, Y., and Polyak, V. J.: Holocene warming in western continental Eurasia driven by glacial retreat and greenhouse forcing, Nat. Geosci., 10, 430–435,, 2017. a

Berger, A.: Long-term variations of daily insolation and Quaternary climatic changes, J. Atmos. Sci., 35, 2362–2367, 1978. a

Berger, A. and Loutre, M.: Insolation values for the climate of the last 10 million years, Quaternary Sci. Rev., 10, 297–317,, 1991. a

Berkelhammer, M., Sinha, A., Stott, L., Cheng, H., Pausata, F. S. R., and Yoshimura, K.: An Abrubt Shift in the Indian Monsoon 4000 Years Ago, in: Geoph. Monog. Series, 75–87,, 2012a. a

Berkelhammer, M., Stott, L., Yoshimura, K., Johnson, K., and Sinha, A.: Synoptic and mesoscale controls on the isotopic composition of precipitation in the western United States, Clim. Dynam., 38, 433–454,, 2012b. a

Bradley, R. S.: Paleoclimate: reconstructing climates of the Quaternary, Elsevier, ISBN: 978-0-121-24010-3, 1999.  a

Brady, E., Stevenson, S., Bailey, D., Liu, Z., Noone, D., Nusbaumer, J., Otto-Bliesner, B. L., Tabor, C., Tomas, R., Wong, T., Zhang, J., and Zhu, J.: The Connected Isotopic Water Cycle in the Community Earth System Model Version 1, J. Adv. Model. Earth Sy., 11, 2547–2566,, 2019. a, b, c, d, e, f

Breitenbach, S. F. M., Rehfeld, K., Goswami, B., Baldini, J. U. L., Ridley, H. E., Kennett, D. J., Prufer, K. M., Aquino, V. V., Asmerom, Y., Polyak, V. J., Cheng, H., Kurths, J., and Marwan, N.: COnstructing Proxy Records from Age models (COPRA), Clim. Past, 8, 1765–1779,, 2012. a

Bretagnon, P. and Francou, G.: Planetary theories in rectangular and spherical variables. VSOP 87 solutions, Astron. Astrophys., 202, 309–315, 1988. a

Bühler, J. C. and Axelsson, J.: SISAL1k_MultiModel, GitHub [code],, last access: 3 June 2022. a

Bühler, J. C., Roesch, C., Kirschner, M., Sime, L., Holloway, M. D., and Rehfeld, K.: Comparison of the oxygen isotope signatures in speleothem records and iHadCM3 model simulations for the last millennium, Clim. Past, 17, 985–1004,, 2021. a, b, c, d, e, f, g, h, i, j, k, l, m

Bühler, J. C., Axelsson, J., Rehfeld, K., LeGrande, A. N., Midhun, M., Sjolte, J., Werner, M., and Yoshimura, K.: Monthly climate variables of isotope-enabled climate model simulations over the last millennium (850–1849 CE), Zenodo [data set],, 2022. a

Caplan, P., Derber, J., Gemmill, W., Hong, S. Y., Pan, H. L., and Parrish, D.: Changes to the 1995 NCEP operational medium-range forecast model analysis-forecast system, Weather Forecast., 12, 581–594,<0581:cttnom>;2, 1997. a

Chatfield, C.: The Analysis of Time Series – An Introduction, Chapman & Hall/CRC, 6th edn., ISBN: 978-1-584-88317-3, 2003. a, b

Cheng, H., Edwards, R. L., Sinha, A., Spötl, C., Yi, L., Chen, S., Kelly, M., Kathayat, G., Wang, X., Li, X., Kong, X., Wang, Y., Ning, Y., and Zhang, H.: The Asian monsoon over the past 640,000 years and ice age terminations, Nature, 534, 640–646,, 2016. a

Chiang, J. C., Herman, M. J., Yoshimura, K., and Fung, I. Y.: Enriched East Asian oxygen isotope of precipitation indicates reduced summer seasonality in regional climate and westerlies, P. Natl. Acad. Sci. USA, 117, 14745–14750,, 2020. a

Colose, C. M., LeGrande, A. N., and Vuille, M.: The influence of volcanic eruptions on the climate of tropical South America during the last millennium in an isotope-enabled general circulation model, Clim. Past, 12, 961–979,, 2016a. a, b, c, d, e, f, g

Colose, C. M., LeGrande, A. N., and Vuille, M.: Hemispherically asymmetric volcanic forcing of tropical hydroclimate during the last millennium, Earth Syst. Dynam., 7, 681–696,, 2016b. a, b, c, d, e

Comas-Bru, L., Harrison, S. P., Werner, M., Rehfeld, K., Scroxton, N., Veiga-Pires, C., and SISAL working group members: Evaluating model outputs using integrated global speleothem records of climate change since the last glacial, Clim. Past, 15, 1557–1579,, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Comas-Bru, L., Atsawawaranunt, K., Harrison, S., and SISAL working group members: SISAL (Speleothem Isotopes Synthesis and AnaLysis Working Group) database version 2.0, University of Reading [data set],, 2020a. a, b

Comas-Bru, L., Rehfeld, K., Roesch, C., Amirnezhad-Mozhdehi, S., Harrison, S. P., Atsawawaranunt, K., Ahmad, S. M., Brahim, Y. A., Baker, A., Bosomworth, M., Breitenbach, S. F. M., Burstyn, Y., Columbu, A., Deininger, M., Demény, A., Dixon, B., Fohlmeister, J., Hatvani, I. G., Hu, J., Kaushal, N., Kern, Z., Labuhn, I., Lechleitner, F. A., Lorrey, A., Martrat, B., Novello, V. F., Oster, J., Pérez-Mejías, C., Scholz, D., Scroxton, N., Sinha, N., Ward, B. M., Warken, S., Zhang, H., and SISAL Working Group members: SISALv2: a comprehensive speleothem isotope database with multiple age–depth models, Earth Syst. Sci. Data, 12, 2579–2606,, 2020b. a, b, c, d, e, f, g

Conroy, J. L., Cobb, K. M., and Noone, D.: Comparison of precipitation isotope variability across the tropical Pacific in observations and SWING2 model simulations, J. Geophys. Res.-Atmos., 118, 5867–5892,, 2013. a

Coplen, T. B.: Reporting of stable carbon, hydrogen, and oxygen isotopic abundances, Tech. rep. IAEA-TECDOC-825, ISSN: 1011-4289, International Atomic Energy Agency (IAEA), (last access: 3 July 2022), 1995. a

Coplen, T. B., Kendall, C., and Hopple, J.: Comparison of stable isotope reference samples, Nature, 302, 28–30, 1983. a, b

Cosford, J., Qing, H., Eglington, B., Mattey, D., Yuan, D., Zhang, M., and Cheng, H.: East Asian monsoon variability since the Mid-Holocene recorded in a high-resolution, absolute-dated aragonite speleothem from eastern China, Earth Planet. Sc. Lett., 275, 296–307,, 2008. a, b

Cox, P. M.: Description of the “TRIFFID” dynamic global vegetation model, Met Office, Hadley Centre, Technical Note 24, 17 pp., 2001. a, b

Crowley, T. J. and Unterman, M. B.: Technical details concerning development of a 1200 yr proxy index for global volcanism, Earth Syst. Sci. Data, 5, 187–197,, 2013. a, b, c, d

Crowley, T. J., Zielinski, G., Vinther, B., Udisti, R., Kreutz, K., Cole-Dai, J., and Castellano, E.: Volcanism and the Little Ice Age, PAGES news, 16, 22–23,, 2008. a, b, c

Dalaiden, Q., Goosse, H., Klein, F., Lenaerts, J. T. M., Holloway, M., Sime, L., and Thomas, E. R.: How useful is snow accumulation in reconstructing surface air temperature in Antarctica? A study combining ice core records and climate models, The Cryosphere, 14, 1187–1207,, 2020. a

Dansgaard, W.: Stable isotopes in precipitation, Tellus, 16, 436–468,, 1964. a, b, c, d, e, f, g

Dee, S., Emile-Geay, J., Evans, M. N., Allam, A., Steig, E. J., and Thompson, D.: PRYSM: An open-source framework for PRoxY System Modeling, with applications to oxygen-isotope systems, J. Adv. Model. Earth Sy., 7, 1220–1247,, 2015. a

Dolman, A. M., Kunz, T., Groeneveld, J., and Laepple, T.: A spectral approach to estimating the timescale-dependent uncertainty of paleoclimate records – Part 2: Application and interpretation, Clim. Past, 17, 825–841,, 2021. a

Dreybrodt, W. and Scholz, D.: Climatic dependence of stable carbon and oxygen isotope signals recorded in speleothems: From soil water to speleothem calcite, Geochim. Cosmochim. Ac., 75, 734–752,, 2011. a

Duan, W., Ruan, J., Luo, W., Li, T., Tian, L., Zeng, G., Zhang, D., Bai, Y., Li, J., Tao, T., Zhang, P., Baker, A., and Tan, M.: The transfer of seasonal isotopic variability between precipitation and drip water at eight caves in the monsoon regions of China, Geochim. Cosmochim. Ac., 183, 250–266, 2016. a

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958,, 2016. a, b, c, d, e

Fairchild, I. J. and Baker, A.: Speleothem science: from process to past environments, John Wiles & Sons, ISBN: 978-1-405-19620-8, 2012. a, b, c, d, e

Field, R. D., Kim, D., LeGrande, A. N., Worden, J., Kelley, M., and Schmidt, G. A.: Evaluating climate model performance in the tropics with retrievals of water isotopic composition from Aura TES, Geophys. Res. Lett., 41, 6030–6036,, 2014. a

Fleitmann, D., Burns, S. J., Mangini, A., Mudelsee, M., Kramers, J., Villa, I., Neff, U., Al-Subbary, A. A., Buettner, A., Hippler, D., and Matter, A.: Holocene ITCZ and Indian monsoon dynamics recorded in stalagmites from Oman and Yemen (Socotra), Quaternary Sci. Rev., 26, 170–188,, 2007. a, b

Fohlmeister, J., Schröder-Ritzrau, A., Scholz, D., Spötl, C., Riechelmann, D. F. C., Mudelsee, M., Wackerbarth, A., Gerdes, A., Riechelmann, S., Immenhauser, A., Richter, D. K., and Mangini, A.: Bunker Cave stalagmites: an archive for central European Holocene climate variability, Clim. Past, 8, 1751–1764,, 2012. a

Fohlmeister, J., Plessen, B., Dudashvili, A. S., Tjallingii, R., Wolff, C., Gafurov, A., and Cheng, H.: Winter precipitation changes during the Medieval Climate Anomaly and the Little Ice Age in arid Central Asia, Quaternary Sci. Rev., 178, 24–36,, 2017. a, b, c, d

Fohlmeister, J., Arps, J., Spötl, C., Schröder-Ritzrau, A., Plessen, B., Günter, C., Frank, N., and Trüssel, M.: Carbon and oxygen isotope fractionation in the water-calcite-aragonite system, Geochim. Cosmochim. Ac., 235, 127–139,, 2018. a

Fohlmeister, J., Voarintsoa, N. R. G., Lechleitner, F. A., Boyd, M., Brandtstätter, S., Jacobson, M. J., and Oster, J. L.: Main controls on the stable carbon isotope composition of speleothems, Geochim. Cosmochim. Ac., 279, 67–87, 2020. a, b, c, d, e

Frisia, S., Borsato, A., and Susini, J.: Synchrotron radiation applications to past volcanism archived in speleothems: An overview, J. Volcanol. Geoth. Res., 177, 96–100,, 2008. a

Gao, C., Robock, A., and Ammann, C.: Volcanic forcing of climate over the past 1500 years: An improved ice core-based index for climate models, J. Geophys. Res., 113, D23111,, 2008. a, b, c, d, e, f

Genty, D., Blamart, D., Ghaleb, B., Plagnes, V., Causse, C., Bakalowicz, M., Zouari, K., Chkir, N., Hellstrom, J., Wainer, K., and Bourges, F.: Timing and dynamics of the last deglaciation from European and North African δ13C stalagmite profiles—comparison with Chinese and South Hemisphere stalagmites, Quaternary Sci. Rev., 25, 2118–2142,, 2006. a, b

Gesch, D. B., Verdin, K. L., and Greenlee, S. K.: New land surface digital elevation model covers the Earth, Eos T. Am. Geophys. Un., 80, 69–70,, 1999. a

Gordon, C., Cooper, C., Senior, C. A., Banks, H. T., Gregory, J. M., Johns, T. C., Mitchell, J. F. B., and Wood, R. A.: The simulation of SST, sea ice extents and ocean heat transports in a version of the Hadley Centre coupled model without flux adjustments, Clim. Dynam., 16, 147–168,, 2000. a, b

Goursaud, S., Masson-Delmotte, V., Favier, V., Orsi, A., and Werner, M.: Water stable isotope spatio-temporal variability in Antarctica in 1960–2013: observations and simulations from the ECHAM5-wiso atmospheric general circulation model, Clim. Past, 14, 923–946,, 2018. a

Grossman, E. L. and Ku, T. L.: Oxygen and carbon isotope fractionation in biogenic aragonite: Temperature effects, Chemical Geology: Isotope Geoscience Section, 59, 59–74,, 1986. a, b

Guðlaugsdóttir, H., Steen-Larsen, H. C., Sjolte, J., Masson-Delmotte, V., Werner, M., and Sveinbjörnsdóttir, E.: The influence of volcanic eruptions on weather regimes over the North Atlantic simulated by ECHAM5/MPI-OM ensemble runs from 800 to 2000 CE, Atmos. Res., 213, 211–223,, 2018. a

Guðlaugsdóttir, H., Sjolte, J., Sveinbjörnsdóttir, Á. E., Werner, M., and Steen-Larsen, H. C.: North Atlantic weather regimes in δ18O of winter precipitation: isotopic fingerprint of the response in the atmospheric circulation after volcanic eruptions, Tellus B, 71, 1–19,, 2019. a

Haese, B., Werner, M., and Lohmann, G.: Stable water isotopes in the coupled atmosphere–land surface model ECHAM5-JSBACH, Geosci. Model Dev., 6, 1463–1480,, 2013. a

Hagemann, S., Arpe, K., and Roeckner, E.: Evaluation of the hydrological cycle in the ECHAM5 model, J. Climate, 19, 3810–3827,, 2006. a

Hansen, M., Scholz, D., Schöne, B. R., and Spötl, C.: Simulating speleothem growth in the laboratory: Determination of the stable isotope fractionation (δ13C and δ18O) between H2O, DIC and CaCO3, Chem. Geol., 509, 20–44,, 2019. a

Hartmann, A. and Baker, A.: Modelling karst vadose zone hydrology and its relevance for paleoclimate reconstruction, Earth-Sci. Rev., 172, 178–192,, 2017. a

Hébert, R., Rehfeld, K., and Laepple, T.: Comparing estimation techniques for temporal scaling in palaeoclimate time series, Nonlin. Processes Geophys., 28, 311–328,, 2021. a

Hendy, C.: The isotopic geochemistry of speleothems–I. The calculation of the effects of different modes of formation on the isotopic composition of speleothems and their applicability as palaeoclimatic indicators, Geochim. Cosmochim. Ac., 35, 801–824,, 1971. a, b

Hijmans, R. J.: raster: Geographic Data Analysis and Modeling, r package version 3.0-12, (last access: 2 July 2022), 2020. a

Holden, J. (Ed.): An introduction to physical geography and the environment, Pearson Education, 3rd edn., ISBN: 9780273740698, 2012. a

Holloway, M. D., Sime, L. C., Singarayer, J. S., Tindall, J. C., and Valdes, P. J.: Simulating the 128-ka Antarctic Climate Response to Northern Hemisphere Ice Sheet Melting Using the Isotope-Enabled HadCM3, Geophys. Res. Lett., 45, 11911–921929,, 2018. a

Hurrell, J. W., Holland, M. M., Gent, P. R., Ghan, S., Kay, J. E., Kushner, P. J., Lamarque, J.-F., Large, W. G., Lawrence, D., Lindsay, K., Lipscomb, W. H., Long, M. C., Mahowald, N., Marsh, D. R., Neale, R. B., Rasch, P., Vavrus, S., Vertenstein, M., Bader, D., Collins, W. D., Hack, J. J., Kiehl, J., and Marshall, S.: The community earth system model: a framework for collaborative research, B. Am. Meteorol. Soc., 94, 1339–1360, 2013. a

Hurtt, G. C., Chini, L. P., Frolking, S., Betts, R., Feddema, J., Fischer, G., Fisk, J., Hibbard, K., Houghton, R., Janetos, A., Jones, C. D., Kindermann, G., Kinoshita, T., Klein Goldewijk, K., Riahi, K., Shevliakova, E., Smith, S., Stehfest, E., Thomson, A., Thornton, P., van Vuuren, D. P., and Wang, Y. P.: Harmonization of land-use scenarios for the period 1500–2100: 600 years of global gridded annual land-use transitions, wood harvest, and resulting secondary lands, Climatic Change, 109, 117–161, 2011. a, b

IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press,, 2013. a

Jahn, A., Lindsay, K., Giraud, X., Gruber, N., Otto-Bliesner, B. L., Liu, Z., and Brady, E. C.: Carbon isotopes in the ocean model of the Community Earth System Model (CESM1), Geosci. Model Dev., 8, 2419–2434,, 2015. a

Jean-Baptiste, P., Genty, D., Fourré, E., and Régnier, E.: Tritium dating of dripwater from Villars Cave (SW-France), Appl. Geochem., 107, 152–158,, 2019. a

Johnston, V. E., Borsato, A., Spötl, C., Frisia, S., and Miorandi, R.: Stable isotopes in caves over altitudinal gradients: fractionation behaviour and inferences for speleothem sensitivity to climate change, Clim. Past, 9, 99–118,, 2013. a

Jones, P. D., Osborn, T. J., and Briffa, K. R.: The Evolution of Climate Over the Last Millennium, Science, 292, 662–667,, 2001. a

Jungclaus, J. H., Keenlyside, N., Botzet, M., Haak, H., Luo, J. J., Latif, M., Marotzke, J., Mikolajewicz, U., and Roeckner, E.: Ocean circulation and tropical variability in the coupled model ECHAM5/MPI-OM, J. Climate, 19, 3952–3972,, 2006. a

Jungclaus, J. H., Lorenz, S. J., Timmreck, C., Reick, C. H., Brovkin, V., Six, K., Segschneider, J., Giorgetta, M. A., Crowley, T. J., Pongratz, J., Krivova, N. A., Vieira, L. E., Solanki, S. K., Klocke, D., Botzet, M., Esch, M., Gayler, V., Haak, H., Raddatz, T. J., Roeckner, E., Schnur, R., Widmann, H., Claussen, M., Stevens, B., and Marotzke, J.: Climate and carbon-cycle variability over the last millennium, Clim. Past, 6, 723–737,, 2010. a, b, c

Jungclaus, J. H., Bard, E., Baroni, M., Braconnot, P., Cao, J., Chini, L. P., Egorova, T., Evans, M., González-Rouco, J. F., Goosse, H., Hurtt, G. C., Joos, F., Kaplan, J. O., Khodri, M., Klein Goldewijk, K., Krivova, N., LeGrande, A. N., Lorenz, S. J., Luterbacher, J., Man, W., Maycock, A. C., Meinshausen, M., Moberg, A., Muscheler, R., Nehrbass-Ahles, C., Otto-Bliesner, B. I., Phipps, S. J., Pongratz, J., Rozanov, E., Schmidt, G. A., Schmidt, H., Schmutz, W., Schurer, A., Shapiro, A. I., Sigl, M., Smerdon, J. E., Solanki, S. K., Timmreck, C., Toohey, M., Usoskin, I. G., Wagner, S., Wu, C.-J., Yeo, K. L., Zanchettin, D., Zhang, Q., and Zorita, E.: The PMIP4 contribution to CMIP6 – Part 3: The last millennium, scientific objective, and experimental design for the PMIP4 past1000 simulations, Geosci. Model Dev., 10, 4005–4033,, 2017. a, b, c

Kageyama, M., Braconnot, P., Harrison, S. P., Haywood, A. M., Jungclaus, J. H., Otto-Bliesner, B. L., Peterschmitt, J.-Y., Abe-Ouchi, A., Albani, S., Bartlein, P. J., Brierley, C., Crucifix, M., Dolan, A., Fernandez-Donado, L., Fischer, H., Hopcroft, P. O., Ivanovic, R. F., Lambert, F., Lunt, D. J., Mahowald, N. M., Peltier, W. R., Phipps, S. J., Roche, D. M., Schmidt, G. A., Tarasov, L., Valdes, P. J., Zhang, Q., and Zhou, T.: The PMIP4 contribution to CMIP6 – Part 1: Overview and over-arching analysis plan, Geosci. Model Dev., 11, 1033–1057,, 2018. a, b, c

Kanamitsu, M., Ebisuzaki, W., Woollen, J., Yang, S. K., Hnilo, J. J., Fiorino, M., and Potter, G. L.: NCEP-DOE AMIP-II reanalysis (R-2), B. Am. Meteorol.l Soc., 83, 1631–1644,, 2002. a

Kanner, L. C., Burns, S. J., Cheng, H., Edwards, R. L., and Vuille, M.: High-resolution variability of the South American summer monsoon over the last seven millennia: insights from a speleothem record from the central Peruvian Andes, Quaternary Sci. Rev., 75, 1–10,, 2013. a

Kathayat, G., Sinha, A., Tanoue, M., Yoshimura, K., Li, H., Zhang, H., and Cheng, H.: Interannual oxygen isotope variability in Indian summer monsoon precipitation reflects changes in moisture sources, Communications Earth & Environment, 2, 1–10,, 2021. a

Kaufmann, G.: Stalagmite growth and palaeo-climate: the numerical perspective, Earth Planet. Sc. Lett., 214, 251–266,, 2003. a

Kendall, C. and Caldwell, E. A.: Chapter 2 – Fundamentals of Isotope Geochemistry, in: Fundamentals of Isotope Geochemistry, Elsevier, Amsterdam, 51–86,, 1998. a

Kennett, D. J., Breitenbach, S. F. M., Aquino, V. V., Asmerom, Y., Awe, J., Baldini, J. U., Bartlein, P., Culleton, B. J., Ebert, C., Jazwa, C., Macri, M. J., Marwan, N., Polyak, V., Prufer, K. M., Ridley, H. E., Sodemann, H., Winterhalder, B., and Haug, G. H.: Development and Disintegration of Maya Political Systems in Response to Climate Change, Science, 338, 788–791,, 2012. a

Konecky, B. L., McKay, N. P., Churakova (Sidorova), O. V., Comas-Bru, L., Dassié, E. P., DeLong, K. L., Falster, G. M., Fischer, M. J., Jones, M. D., Jonkers, L., Kaufman, D. S., Leduc, G., Managave, S. R., Martrat, B., Opel, T., Orsi, A. J., Partin, J. W., Sayani, H. R., Thomas, E. K., Thompson, D. M., Tyler, J. J., Abram, N. J., Atwood, A. R., Cartapanis, O., Conroy, J. L., Curran, M. A., Dee, S. G., Deininger, M., Divine, D. V., Kern, Z., Porter, T. J., Stevenson, S. L., von Gunten, L., and Iso2k Project Members: The Iso2k database: a global compilation of paleo-δ18O and δ2H records to aid understanding of Common Era climate, Earth Syst. Sci. Data, 12, 2261–2288,, 2020. a

Kunz, T., Dolman, A. M., and Laepple, T.: A spectral approach to estimating the timescale-dependent uncertainty of paleoclimate records – Part 1: Theoretical concept, Clim. Past, 16, 1469–1492,, 2020 (data available at:, last access: 1 July 2022). a

Lachniet, M. S.: Climatic and environmental controls on speleothem oxygen-isotope values, Quaternary Sci. Rev., 28, 412–432,, 2009. a

Laepple, T. and Huybers, P.: Global and regional variability in marine surface temperatures, Geophys. Res. Lett., 41, 2528–2534,, 2014a. a, b

Laepple, T. and Huybers, P.: Ocean surface temperature variability: Large model–data differences at decadal and longer periods, P. Natl. Acad. Sci. USA, 111, 16682–16687,, 2014b. a

Landrum, L., Otto-Bliesner, B. L., Wahl, E. R., Conley, A., Lawrence, P. J., Rosenbloom, N., and Teng, H.: Last millennium climate and its variability in CCSM4, J. Climate, 26, 1085–1111,, 2013. a

Langebroek, P. M., Werner, M., and Lohmann, G.: Climate information imprinted in oxygen-isotopic composition of precipitation in Europe, Earth Planet. Sc. Lett., 311, 144–154,, 2011. a

Latombe, G., Burke, A., Vrac, M., Levavasseur, G., Dumas, C., Kageyama, M., and Ramstein, G.: Comparison of spatial downscaling methods of general circulation model results to study climate variability during the Last Glacial Maximum, Geosci. Model Dev., 11, 2563–2579,, 2018. a

Lean, J.: Calculations of Solar Irradiance: monthly means from 1882 to 2008, annual means from 1610 to 2008, Freie Universität Berlin [data set], (last access: 7 July 2022), 2009. a

Lechleitner, F. A., Day, C. C., Kost, O., Wilhelm, M., Haghipour, N., Henderson, G. M., and Stoll, H. M.: Stalagmite carbon isotopes suggest deglacial increase in soil respiration in western Europe driven by temperature change, Clim. Past, 17, 1903–1918,, 2021. a

Legrande, A. N. and Anchukaitis, K. J.: Volcanic eruptions and climate, PAGES Magazine, Vol. 23, 2015. a

LeGrande, A. N. and Schmidt, G. A.: Sources of Holocene variability of oxygen isotopes in paleoclimate archives, Clim. Past, 5, 441–455,, 2009. a

Lewis, S. C. and LeGrande, A. N.: Stability of ENSO and its tropical Pacific teleconnections over the Last Millennium, Clim. Past, 11, 1347–1360,, 2015. a, b, c, d, e

Lone, M. A., Ahmad, S. M., Dung, N. C., Shen, C.-C., Raza, W., and Kumar, A.: Speleothem based 1000-year high resolution record of Indian monsoon variability during the last deglaciation, Palaeogeogr. Palaeocl., 395, 1–8,, 2014. a, b, c

Lundeen, Z., Brunelle, A., Burns, S. J., Polyak, V., and Asmerom, Y.: A speleothem record of Holocene paleoclimate from the northern Wasatch Mountains, southeast Idaho, USA, Quatern. Int., 310, 83–95,, 2013. a

MacFarling Meure, C., Etheridge, D., Trudinger, C., Steele, P., Langenfelds, R., Van Ommen, T., Smith, A., and Elkins, J.: Law Dome CO2, CH4 and N2O ice core records extended to 2000 years BP, Geophys. Res. Lett., 33, L14810,, 2006. a, b

Marland, G., Boden, T. A., and Andres, R. J.: Global, Regional, and National Fossil Fuel CO2 Emissions, in: Trends: A Compendium of Data on Global Change, Tech. rep., Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tenn., USA, 2003. a

Marsland, S. J., Haak, H., Jungclaus, J. H., Latif, M., and Röske, F.: The Max-Planck-Institute global ocean/sea ice model with orthogonal curvilinear coordinates, Ocean Model., 5, 91–127,, 2002. a

MATLAB: Version (R2018b), The MathWorks Inc., Natick, Massachusetts, 2018. a

McDermott, F., Mattey, D. P., and Hawkesworth, C.: Centennial-scale holocene climate variability revealed by a high-resolution speleothem δ18O record from SW Ireland, Science, 294, 1328–1331,, 2001. a, b

Medina-Elizalde, M., Burns, S. J., Polanco-Martínez, J. M., Beach, T., Lases-Hernández, F., Shen, C.-C., and Wang, H.-C.: High-resolution speleothem record of precipitation from the Yucatan Peninsula spanning the Maya Preclassic Period, Global Planet. Change, 138, 93–102,, 2016. a

Midhun, M. and Ramesh, R.: Validation of δ18O as a proxy for past monsoon rain by multi-GCM simulations, Clim. Dynam., 46, 1371–1385,, 2016. a

Midhun, M., Stevenson, S., and Cole, J.: Oxygen isotopic signatures of major climate modes and implications for detectability in speleothems, Geophys. Res. Lett., 48, e2020GL089515,, 2021. a, b, c, d, e, f

Miller, R. L., Schmidt, G. A., Nazarenko, L. S., Tausnev, N., Bauer, S. E., Delgenio, A. D., Kelley, M., Lo, K. K., Ruedy, R., Shindell, D. T., Aleinov, I., Bauer, M., Bleck, R., Canuto, V., Chen, Y., Cheng, Y., Clune, T. L., Faluvegi, G., Hansen, J. E., Healy, R. J., Kiang, N. Y., Koch, D., Lacis, A. A., Legrande, A. N., Lerner, J., Menon, S., Oinas, V., Pérez García-Pando, C., Perlwitz, J. P., Puma, M. J., Rind, D., Romanou, A., Russell, G. L., Sato, M., Sun, S., Tsigaridis, K., Unger, N., Voulgarakis, A., Yao, M. S., and Zhang, J.: CMIP5 historical simulations (1850–2012) with GISS ModelE2, J. Adv. Model. Earth Sy., 6, 441–477,, 2014. a

Morice, C. P., Kennedy, J. J., Rayner, N. A., and Jones, P. D.: Quantifying uncertainties in global and regional temperature change using an ensemble of observational estimates: The HadCRUT4 data set, J. Geophys. Res., 117, D08101,, 2012. a

Muscheler, R., Joos, F., Beer, J., Müller, S. A., Vonmoos, M., and Snowball, I.: Solar activity during the last 1000 yr inferred from radionuclide records, Quaternary Sci. Rev., 26, 82–97,, 2007. a

Muscheler, R., Adolphi, F., Herbst, K., and Nilsson, A.: The Revised Sunspot Record in Comparison to Cosmogenic Radionuclide-Based Solar Activity Reconstructions, Sol. Phys., 291, 3025–3043,, 2016. a, b, c, d

National Geophysical Data Center: 5-minute Gridded Global Relief Data (ETOPO5), National Geophysical Data Center [data set], NOAA,, 1993. a

Neale, R. B., Chen, C. C., Gettelman, A., Lauritzen, P. H., Park, S., Williamson, D. L., Conley, A. J., Garcia, R., Kinnison, D., Lamarque, J. F., and Marsh, D.: Description of the NCAR community atmosphere model (CAM 5.0), NCAR Tech. Note NCAR/TN-486+ TR, 1, 1–12, 2010. a

Neff, U., Burns, S. J., Mangini, A., Mudelsee, M., Fleitmann, D., and Matter, A.: Strong coherence between solar variability and the monsoon in Oman between 9 and 6 kyr ago, Nature, 411, 290–293,, 2001. a, b, c

Neukom, R., Steiger, N., Gómez-Navarro, J. J., Wang, J., and Werner, J. P.: No evidence for globally coherent warm and cold periods over the preindustrial Common Era, Nature, 571, 550–554,, 2019. a, b

NOAA National Geophysical Data Center: ETOPO1 1 Arc-Minute Global Relief Model, NOAA National Centers for Environmental Information [data set],, 2009. a

Noronha, A. L., Johnson, K. R., Hu, C., Ruan, J., Southon, J. R., and Ferguson, J. E.: Assessing influences on speleothem dead carbon variability over the Holocene: Implications for speleothem-based radiocarbon calibration, Earth Planet. Sc. Lett., 394, 20–29,, 2014. a

Novello, V., Cruz, F., Karmann, I., Burns, S., Stríkis, N., Vuille, M., Cheng, H., Lawrence Edwards, R., Santos, R., Frigo, E., and Barreto, E.: Multidecadal climate variability in Brazil's Nordeste during the last 3000 years based on speleothem isotope records, Geophys. Res. Lett., 39, L23706,, 2012. a

Novello, V. F., Cruz, F. W., McGlue, M. M., Wong, C. I., Ward, B. M., Vuille, M., Santos, R. A., Jaqueto, P., Pessenda, L. C., Atorre, T., Ribeiro, L. M., Karmann, I., Barreto, E. S., Cheng, H., Edwards, R. L., Paula, M. S., and Scholz, D.: Vegetation and environmental changes in tropical South America from the last glacial to the Holocene documented by multiple cave sediment proxies, Earth Planet. Sc. Lett., 524, 115717,, 2019. a

Novello, V. F., William da Cruz, F., Vuille, M., Pereira Silveira Campos, J. L., Stríkis, N. M., Apaéstegui, J., Moquet, J. S., Azevedo, V., Ampuero, A., Utida, G., Wang, X., Paula-Santos, G. M., Jaqueto, P., Ruiz Pessenda, L. C., Breecker, D. O., and Karmann, I.: Investigating δ13C values in stalagmites from tropical South America for the last two millennia, Quaternary Sci. Rev., 255, 106822,, 2021. a

Nusbaumer, J., Wong, T. E., Bardeen, C., and Noone, D.: Evaluating hydrological processes in the Community Atmosphere Model Version 5 (CAM5) using stable isotope ratios of water, J. Adv. Model. Earth Sy., 9, 949–977,, 2017. a, b

Oleson, K. W., Bonan, G. B., Feddema, J., Vertenstein, M., and Kluzek, E.: Technical Description of an Urban Parameterization for the Community Land Model (CLMU), No. NCAR/TN-480+STR, University Corporation for Atmospheric Research,, 2010. a

Otto-Bliesner, B. L., Brady, E. C., Fasullo, J., Jahn, A., Landrum, L., Stevenson, S., Rosenbloom, N., Mai, A., and Strand, G.: Climate variability and change since 850  CE: An ensemble approach with the Community Earth System Model, B. Am. Meteorol. Soc., 97, 735–754, 2016. a

Owen, R., Day, C., Hu, C.-Y., Liu, Y.-H., Pointing, M., Blättler, C., and Henderson, G.: Calcium isotopes in caves as a proxy for aridity: Modern calibration and application to the 8.2 kyr event, Earth Planet. Sc. Lett., 443, 129–138,, 2016. a

PAGES2k-Consortium: Consistent multidecadal variability in global temperature reconstructions and simulations over the Common Era, Nat. Geosci., 12, 643–649,, 2019. a

PAGES Hydro2k Consortium: Comparing proxy and model estimates of hydroclimate variability and change over the Common Era, Clim. Past, 13, 1851–1900,, 2017. a, b, c, d

Pardaens, A. K., Banks, H. T., Gregory, J. M., and Rowntree, P. R.: Freshwater transports in HadCM3, Clim. Dynam., 21, 177–195,, 2003. a

Parker, S. E., Harrison, S. P., Comas-Bru, L., Kaushal, N., LeGrande, A. N., and Werner, M.: A data–model approach to interpreting speleothem oxygen isotope records from monsoon regions, Clim. Past, 17, 1119–1138,, 2021. a, b, c, d

Paul, J., Fortuin, F., and Kelder, H.: An ozone climatology based on ozonesonde and satellite measurements, J. Geophys. Res., 103, 31709–31734,, 1998. a

Phipps, S. J., Mcgregor, H. V., Gergis, J., Gallant, A. J., Neukom, R., Stevenson, S., Ackerley, D., Brown, J. R., Fischer, M. J., and Van Ommen, T. D.: Paleoclimate data-model comparison and the role of climate forcings over the past 1500 years, J. Climate, 26, 6915–6936,, 2013. a

Pierce, D.: ncdf4: Interface to Unidata netCDF (Version 4 or Earlier) Format Data Files, r package version 1.17, R Foundation, (last access: 2 July 2022), 2019. a

Polag, D., Scholz, D., Mühlinghaus, C., Spötl, C., Schröder-Ritzrau, A., Segl, M., and Mangini, A.: Stable isotope fractionation in speleothems: Laboratory experiments, Chem. Geol., 279, 31–39,, 2010. a

Pongratz, J., Reick, C., Raddatz, T., and Claussen, M.: A reconstruction of global agricultural areas and land cover for the last millennium, Global Biogeochem. Cy., 22, GB3018,, 2008. a, b, c, d

Pope, V. D., Gallani, M. L., Rowntree, P. R., and Stratton, R. A.: The impact of new physical parametrizations in the Hadley Centre climate model: HadAM3, Clim. Dynam., 16, 123–146,, 2000. a

R Core Team: R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, (last access: 2 July 2022), 2020. a, b

Rehfeld, K. and Kurths, J.: Similarity estimators for irregular and age-uncertain time series, Clim. Past, 10, 107–122,, 2014 (code available at:, last access: 1 July 2022). a, b, c, d, e

Rehfeld, K., Marwan, N., Heitzig, J., and Kurths, J.: Comparison of correlation analysis techniques for irregularly sampled time series, Nonlin. Processes Geophys., 18, 389–404,, 2011. a

Rehfeld, K., Münch, T., Ho, S. L., and Laepple, T.: Golbal paterns of declining temperature variability from the Last Glacial Maximum to the Holocene, Nature, 554, 356–359,, 2018. a

Ridley, H. E., Asmerom, Y., Baldini, J. U., Breitenbach, S. F., Aquino, V. V., Prufer, K. M., Culleton, B. J., Polyak, V., Lechleitner, F. A., Kennett, D. J., Zhang, M., Marwan, N., Macpherson, C. G., Baldini, L. M., Xiao, T., Peterkin, J. L., Awe, J., and Haug, G. H.: Aerosol forcing of the position of the intertropical convergence zone since AD 1550, Nat. Geosci., 8, 195–200, 2015. a, b

Risi, C., Noone, D., Worden, J., Frankenberg, C., Stiller, G., Kiefer, M., Funke, B., Walker, K., Bernath, P., Schneider, M., Wunch, D., Sherlock, V., Deutscher, N., Griffith, D., Wennberg, P. O., Strong, K., Smale, D., Mahieu, E., Barthlott, S., Hase, F., García, O., Notholt, J., Warneke, T., Toon, G., Sayres, D., Bony, S., Lee, J., Brown, D., Uemura, R., and Sturm, C.: Process-evaluation of tropospheric humidity simulated by general circulation models using water vapor isotopologues: 1. Comparison between models and observations, J. Geophys. Res., 117, D05303,, 2012. a, b, c

Roeckner, E., Bäuml, G., Bonaventura, L., Brokopf, R., Esch, M., Giorgetta, M., Hagemann, S., Kirchner, I., Kornblueh, L., Manzini, E., Rhodin, A., Schlese, U., Schulzweida, U., and Tompkins, A.: The Atmospheric General Circulation Model ECHAM5. Part I: Model description, Tech. rep., Max Planck Institute for Meteorology, Hamburg, Germany, Report No. 349, ISSN: 0937-1060, 2003. a

Romanek, C. S., Grossman, E. L., and Morse, J. W.: Carbon isotopic fractionation in synthetic aragonite and calcite: Effects of temperature and precipitation rate, Geochim. Cosmochim. Ac., 56, 419–430,, 1992. a, b

Rozanski, K., Araguás-Araguás, L., and Gonfiantini, R.: Relation between long-term trends of oxygen-18 isotope composition of precipitation and climate, Science, 258, 981–985,, 1992. a

Schmidt, G. A., Hoffmann, G., Shindell, D. T., and Hu, Y.: Modeling atmospheric stable water isotopes and the potential for constraining cloud processes and stratosphere-troposphere water exchange, J. Geophys. Res., 110, D21314,, 2005. a, b

Schmidt, G. A., Ruedy, R., Hansen, J. E., Aleinov, I., Bell, N., Bauer, M., Bauer, S., Cairns, B., Canuto, V., Cheng, Y., Del Genio, A., Faluvegi, G., Friend, A. D., Hall, T. M., Hu, Y., Kelley, M., Kiang, N. Y., Koch, D., Lacis, A. A., Lerner, J., Lo, K. K., Miller, R. L., Nazarenko, L., Oinas, V., Perlwitz, J., Perlwitz, J., Rind, D., Romanou, A., Russell, G. L., Sato, M., Shindell, D. T., Stone, P. H., Sun, S., Tausnev, N., Thresher, D., and Yao, M. S.: Present-day atmospheric simulations using GISS ModelE: Comparison to in situ, satellite, and reanalysis data, J. Climate, 19, 153–192,, 2006. a

Schmidt, G. A., Jungclaus, J. H., Ammann, C. M., Bard, E., Braconnot, P., Crowley, T. J., Delaygue, G., Joos, F., Krivova, N. A., Muscheler, R., Otto-Bliesner, B. L., Pongratz, J., Shindell, D. T., Solanki, S. K., Steinhilber, F., and Vieira, L. E. A.: Climate forcing reconstructions for use in PMIP simulations of the last millennium (v1.0), Geosci. Model Dev., 4, 33–45,, 2011. a, b, c, d, e

Schmidt, G. A., Jungclaus, J. H., Ammann, C. M., Bard, E., Braconnot, P., Crowley, T. J., Delaygue, G., Joos, F., Krivova, N. A., Muscheler, R., Otto-Bliesner, B. L., Pongratz, J., Shindell, D. T., Solanki, S. K., Steinhilber, F., and Vieira, L. E. A.: Climate forcing reconstructions for use in PMIP simulations of the Last Millennium (v1.1), Geosci. Model Dev., 5, 185–191,, 2012. a, b, c, d, e, f

Schmidt, G. A., Kelley, M., Nazarenko, L., Ruedy, R., Russell, G. L., Aleinov, I., Bauer, M., Bauer, S. E., Bhat, M. K., Bleck, R., Canuto, V., Chen, Y. H., Cheng, Y., Clune, T. L., Del Genio, A., De Fainchtein, R., Faluvegi, G., Hansen, J. E., Healy, R. J., Kiang, N. Y., Koch, D., Lacis, A. A., Legrande, A. N., Lerner, J., Lo, K. K., Matthews, E. E., Menon, S., Miller, R. L., Oinas, V., Oloso, A. O., Perlwitz, J. P., Puma, M. J., Putman, W. M., Rind, D., Romanou, A., Sato, M., Shindell, D. T., Sun, S., Syed, R. A., Tausnev, N., Tsigaridis, K., Unger, N., Voulgarakis, A., Yao, M. S., and Zhang, J.: Configuration and assessment of the GISS ModelE2 contributions to the CMIP5 archive, J. Adv. Model. Earth Sy., 6, 141–184,, 2014. a, b, c, d

Scholz, D., Frisia, S., Borsato, A., Spötl, C., Fohlmeister, J., Mudelsee, M., Miorandi, R., and Mangini, A.: Holocene climate variability in north-eastern Italy: potential influence of the NAO and solar activity recorded by speleothem data, Clim. Past, 8, 1367–1383,, 2012. a, b

Schurer, A. P., Tett, S. F. B., and Hegerl, G. C.: Small influence of solar variability on climate over the past millennium, Nat. Geosci., 7, 104–108,, 2014. a, b, c

Schwarcz, H. P., Harmon, R. S., Thompson, P., and Ford, D. C.: Stable isotope studies of fluid inclusions in speleothems and their paleoclimatic significance, Geochim. Cosmochim. Ac., 40, 657–665,, 1976. a

Shukla, P., Skea, J., Calvo Buendia, E., et al.: IPCC, 2019: Climate Change and Land: an IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems, in press, 2019. a

Sigl, M., Winstrup, M., McConnell, J. R., Welten, K. C., Plunkett, G., Ludlow, F., Büntgen, U., Caffee, M., Chellman, N., Dahl-Jensen, D., Fischer, H., Kipfstuhl, S., Kostick, C., Maselli, O. J., Mekhaldi, F., Mulvaney, R., Muscheler, R., Pasteris, D. R., Pilcher, J. R., Salzer, M., Schüpbach, S., Steffensen, J. P., Vinther, B. M., and Woodruff, T. E.: Timing and climate forcing of volcanic eruptions for the past 2,500 years, Nature, 523, 543–549, 2015. a

Sime, L. C., Tindall, J. C., Wolff, E. W., Connolley, W. M., and Valdes, P. J.: Antarctic isotopic thermometer during a CO2 forced warming event, J. Geophys. Res., 113, D24119,, 2008. a

Sime, L. C., Kohfeld, K. E., Le Quéré, C., Wolff, E. W., de Boer, A. M., Graham, R. M., and Bopp, L.: Southern Hemisphere westerly wind changes during the Last Glacial Maximum: Model-data comparison, Quaternary Sci. Rev., 64, 104–120,, 2013. a

Sjolte, J., Sturm, C., Adolphi, F., Vinther, B. M., Werner, M., Lohmann, G., and Muscheler, R.: Solar and volcanic forcing of North Atlantic climate inferred from a process-based reconstruction, Clim. Past, 14, 1179–1194,, 2018. a, b, c, d

Steinhilber, F., Beer, J., and Fröhlich, C.: Total solar irradiance during the Holocene, Geophys. Res. Lett., 36, L19704,, 2009. a, b, c

Stevenson, S., Otto-Bliesner, B. L., Brady, E. C., Nusbaumer, J., Tabor, C., Tomas, R., Noone, D. C., and Liu, Z.: Volcanic Eruption Signatures in the Isotope-Enabled Last Millennium Ensemble, Paleoceanography and Paleoclimatology, 34, 1534–1552,, 2019. a, b

Sturm, C., Zhang, Q., and Noone, D.: An introduction to stable water isotopes in climate models: benefits of forward proxy modelling for paleoclimatology, Clim. Past, 6, 115–129,, 2010. a, b

Sun, Z., Yang, Y., Zhao, J., Tian, N., and Feng, X.: Potential ENSO effects on the oxygen isotope composition of modern speleothems: Observations from Jiguan Cave, central China, J. Hydrol., 566, 164–174,, 2018. a

Tan, L., Cai, Y., An, Z., Edwards, R. L., Cheng, H., Shen, C.-C., and Zhang, H.: Centennial- to decadal-scale monsoon precipitation variability in the semi-humid region, northern China during the last 1860 years: Records from stalagmites in Huangye Cave, The Holocene, 21, 287–296,, 2011. a

Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An overview of CMIP5 and the experiment design, B. Am. Meteorol. Soc., 93, 485–498, 2012. a, b, c

Tierney, J. E., Zhu, J., King, J., Malevich, S. B., Hakim, G. J., and Poulsen, C. J.: Glacial cooling and climate sensitivity revisited, Nature, 584, 569–573, 2020. a

Tindall, J., Flecker, R., Valdes, P., Schmidt, D. N., Markwick, P., and Harris, J.: Modelling the oxygen isotope distribution of ancient seawater using a coupled ocean-atmosphere GCM: Implications for reconstructing early Eocene climate, Earth Planet. Sc. Lett., 292, 265–273,, 2010. a

Tindall, J. C., Valdes, P. J., and Sime, L. C.: Stable water isotopes in HadCM3: Isotopic signature of El Niño-Southern Oscillation and the tropical amount effect, J. Geophys. Res., 114, D04111,, 2009. a, b, c, d, e, f

Treble, P. C., Fairchild, I. J., Griffiths, A., Baker, A., Meredith, K. T., Wood, A., and McGuire, E.: Impacts of cave air ventilation and in-cave prior calcite precipitation on Golgotha Cave dripwater chemistry, southwest Australia, Quaternary Sci. Rev., 127, 61–72,, 2015. a

Tremaine, D. M. and Froelich, P. N.: Speleothem trace element signatures: A hydrologic geochemical study of modern cave dripwaters and farmed calcite, Geochim. Cosmochim. Ac., 121, 522–545,, 2013. a, b

Tremaine, D. M., Froelich, P. N., and Wang, Y.: Speleothem calcite farmed in situ: Modern calibration of δ18O and δ13C paleoclimate proxies in a continuously-monitored natural cave system, Geochim. Cosmochim. Ac., 75, 4929–4950,, 2011. a, b

Valdes, P. J., Armstrong, E., Badger, M. P. S., Bradshaw, C. D., Bragg, F., Crucifix, M., Davies-Barnard, T., Day, J. J., Farnsworth, A., Gordon, C., Hopcroft, P. O., Kennedy, A. T., Lord, N. S., Lunt, D. J., Marzocchi, A., Parry, L. M., Pope, V., Roberts, W. H. G., Stone, E. J., Tourte, G. J. L., and Williams, J. H. T.: The BRIDGE HadCM3 family of climate models: HadCM3@Bristol v1.0, Geosci. Model Dev., 10, 3715–3743,, 2017. a

Van Rampelbergh, M., Fleitmann, D., Verheyden, S., Cheng, H., Edwards, L., De Geest, P., De Vleeschouwer, D., Burns, S. J., Matter, A., Claeys, P., and Keppens, E.: Mid- to late Holocene Indian Ocean Monsoon variability recorded in four speleothems from Socotra Island, Yemen, Quaternary Sci. Rev., 65, 129–142,, 2013. a

Vieira, L. E. A., Solanki, S. K., Krivova, N. A., and Usoskin, I.: Evolution of the solar irradiance during the Holocene, Astron. Astrophys., 531, A6,, 2011. a, b, c

Voarintsoa, N. R. G., Wang, L., Railsback, L. B., Brook, G. A., Liang, F., Cheng, H., and Edwards, R. L.: Multiple proxy analyses of a U/Th-dated stalagmite to reconstruct paleoenvironmental changes in northwestern Madagascar between 370 CE and 1300 CE, Palaeogeogr. Palaeocl., 469, 138–155,, 2017. a

Wackerbarth, A., Scholz, D., Fohlmeister, J., and Mangini, A.: Modelling the δ18O value of cave drip water and speleothem calcite, Earth Planet. Sc. Lett., 299, 387–397,, 2010. a

Wang, J. K., Johnson, K. R., Borsato, A., Amaya, D. J., Griffiths, M. L., Henderson, G. M., Frisia, S., and Mason, A.: Hydroclimatic variability in Southeast Asia over the past two millennia, Earth Planet. Sc. Lett., 525, 115737,, 2019. a

Wang, Y., Lean, J. L., and Sheeley Jr, N. R.: Modeling the Sun's Magnetic Field and Irradiance since 1713, Astrophys. J., 625, 522–538,, 2005. a, b

Warken, S. F., Fohlmeister, J., Schröder-Ritzrau, A., Constantin, S., Spötl, C., Gerdes, A., Esper, J., Frank, N., Arps, J., Terente, M., Riechelmann, D. F., Mangini, A., and Scholz, D.: Reconstruction of late Holocene autumn/winter precipitation variability in SW Romania from a high-resolution speleothem trace element record, Earth Planet. Sc. Lett., 499, 122–133,, 2018. a

Warken, S. F., Schorndorf, N., Stinnesbeck, W., Hennhoefer, D., Stinnesbeck, S. R., Förstel, J., Steidle, S. D., Avilés Olguin, J., and Frank, N.: Solar forcing of early Holocene droughts on the Yucatán peninsula, Scientific Reports, 11, 13885,, 2021. a

Werner, M.: Modelling stable water isotopes: Status and perspectives, EPJ Web Conf., 9, 73–82,, 2010. a, b

Werner, M., Langebroek, P. M., Carlsen, T., Herold, M., and Lohmann, G.: Stable water isotopes in the ECHAM5 general circulation model: Toward high-resolution isotope modeling on a global scale, J. Geophys. Res., 116, D15109,, 2011. a, b

Werner, M., Haese, B., Xu, X., Zhang, X., Butzin, M., and Lohmann, G.: Glacial–interglacial changes in H218O, HDO and deuterium excess – results from the fully coupled ECHAM5/MPI-OM Earth system model, Geosci. Model Dev., 9, 647–670,, 2016. a, b, c, d, e, f, g, h, i, j, k

Wickham, H.: The Split-Apply-Combine Strategy for Data Analysis, J. Stat. Softw., 40, 1–29,, 2011. a

Wickham, H.: ggplot2: Elegant Graphics for Data Analysis, Springer-Verlag New York, (last access: 2 July 2022), 2016. a

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., Takahashi, K., Vaughan, D., Wilke, C., Woo, K., and Yutani, H.: Welcome to the tidyverse, Journal of Open Source Software, 4, 1686,, 2019. a

Wickham, H., François, R., Henry, L., and Müller, K.: dplyr: A Grammar of Data Manipulation, r package version 1.0.7, (last access: 2 July 2022), 2021. a

Williams, P. W. and Ford, D. C.: Global distribution of carbonate rocks, Zeitschrift für Geomorphologie Supplementband, 147, p. 1, 2006. a

Wong, C. I. and Breecker, D. O.: Advancements in the use of speleothems as climate archives, Quaternary Sci. Rev., 127, 1–18,, 2015. a

Xi, X.: A Review of Water Isotopes in Atmospheric General Circulation Models: Recent Advances and Future Prospects, Int. J. Atmos. Sci., 2014, 250920,, 2014. a

Yoshimura, K., Kanamitsu, M., Noone, D., and Oki, T.: Historical isotope simulation using Reanalysis atmospheric data, J. Geophys. Res., 113, D19108,, 2008. a, b, c, d, e, f, g, h, i

Zeileis, A. and Grothendieck, G.: zoo: S3 Infrastructure for Regular and Irregular Time Series, J. Stat. Softw., 14, 1–27,, 2005. a

Zhang, J., Liu, Z., Brady, E. C., Oppo, D. W., Clark, P. U., Jahn, A., Marcott, S. A., and Lindsay, K.: Asynchronous warming and δ18O evolution of deep Atlantic water masses during the last deglaciation, P. Natl. Acad. Sci. USA, 114, 11075–11080,, 2017. a

Zhang, X., Sun, Z., Guan, H., Zhang, X., Wu, H., and Huang, Y.: GCM simulations of stable isotopes in the water cycle in comparison with GNIP observations over east asia, Acta Meteorol. Sin., 26, 420–437,, 2012. a, b

Zhu, J., Liu, Z., Brady, E., Otto-Bliesner, B., Zhang, J., Noone, D., Tomas, R., Nusbaumer, J., Wong, T., Jahn, A., and Tabor, C.: Reduced ENSO variability at the LGM revealed by an isotope-enabled Earth system model, Geophys. Res. Lett., 44, 6984–6992,, 2017. a

Short summary
We collected and standardized the output of five isotope-enabled simulations for the last millennium and assess differences and similarities to records from a global speleothem database. Modeled isotope variations mostly arise from temperature differences. While lower-resolution speleothems do not capture extreme changes to the extent of models, they show higher variability on multi-decadal timescales. As no model excels in all comparisons, we advise a multi-model approach where possible.