Articles | Volume 17, issue 3
Clim. Past, 17, 1119–1138, 2021
Clim. Past, 17, 1119–1138, 2021

Research article 04 Jun 2021

Research article | 04 Jun 2021

A data–model approach to interpreting speleothem oxygen isotope records from monsoon regions

A data–model approach to interpreting speleothem oxygen isotope records from monsoon regions
Sarah E. Parker1, Sandy P. Harrison1, Laia Comas-Bru1, Nikita Kaushal2, Allegra N. LeGrande3, and Martin Werner4 Sarah E. Parker et al.
  • 1School of Archaeology, Geography and Environmental Science, University of Reading, Reading, UK
  • 2Asian School of the Environment, Nanyang Technological University, Singapore, Singapore
  • 3NASA Goddard Institute for Space Studies and Center for Climate Systems Research, Columbia University, New York, USA
  • 4Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany

Correspondence: Sarah E. Parker (


Reconstruction of past changes in monsoon climate from speleothem oxygen isotope (δ18O) records is complex because δ18O signals can be influenced by multiple factors including changes in precipitation, precipitation recycling over land, temperature at the moisture source, and changes in the moisture source region and transport pathway. Here, we analyse >150 speleothem records of the Speleothem Isotopes Synthesis and AnaLysis (SISAL) database to produce composite regional trends in δ18O in monsoon regions; compositing minimises the influence of site-specific karst and cave processes that can influence individual site records. We compare speleothem δ18O observations with isotope-enabled climate model simulations to investigate the specific climatic factors causing these regional trends. We focus on differences in δ18O signals between the mid-Holocene, the peak of the Last Interglacial (Marine Isotope Stage 5e) and the Last Glacial Maximum as well as on δ18O evolution through the Holocene. Differences in speleothem δ18O between the mid-Holocene and the Last Interglacial in the East Asian and Indian monsoons are small, despite the larger summer insolation values during the Last Interglacial. Last Glacial Maximum δ18O values are significantly less negative than interglacial values. Comparison with simulated glacial–interglacial δ18O shows that changes are principally driven by global shifts in temperature and regional precipitation. Holocene speleothem δ18O records show distinct and coherent regional trends. Trends are similar to summer insolation in India, China and southwestern South America, but they are different in the Indonesian–Australian region. Redundancy analysis shows that 37 % of Holocene variability can be accounted for by latitude and longitude, supporting the differentiation of records into individual monsoon regions. Regression analysis of simulated precipitation δ18O and climate variables show significant relationships between global Holocene monsoon δ18O trends and changes in precipitation, atmospheric circulation and (to a lesser extent) source area temperature, whereas precipitation recycling is non-significant. However, there are differences in regional-scale mechanisms: there are clear relationships between changes in precipitation and δ18O for India, southwestern South America and the Indonesian–Australian regions but not for the East Asian monsoon. Changes in atmospheric circulation contribute to δ18O trends in the East Asian, Indian and Indonesian–Australian monsoons, and a weak source area temperature effect is observed over southern and central America and Asia. Precipitation recycling is influential in southwestern South America and southern Africa. Overall, our analyses show that it is possible to differentiate the impacts of specific climatic mechanisms influencing precipitation δ18O and use this analysis to interpret changes in speleothem δ18O.

1 Introduction

The oxygen isotopic (δ18O; the 18O/16O ratio relative to a standard, in per mille, ‰) composition of speleothems is widely used to infer past regional climates (Bar-Matthews et al., 1997; McDermott, 2004; Wang et al., 2008). Speleothem oxygen isotope (δ18Ospel) signals are inherited from δ18O in precipitation (δ18Oprecip) above the cave, which in turn is determined by the initial δ18O of water vapour as it evaporates at the oceanic moisture source region, the degree of rainout and evaporation from source to cave site, and air temperature changes encountered throughout the moisture transport pathway (Fairchild and Baker, 2012; Lachniet, 2009). Understanding the effects and contribution of each of these climate processes to δ18Oprecip (and therefore δ18Ospel) is essential to inferring palaeoclimate from speleothem δ18O records.

Initial δ18O is determined by oceanic δ18O at the evaporative moisture source region (Craig and Gordon, 1965), which varies spatially (LeGrande and Schmidt, 2006) and through time (e.g. Waelbroeck et al., 2002). During evaporation from the moisture source, 16O is preferentially incorporated into the vapour, whilst subsequent fractionation during atmospheric transport occurs by Rayleigh distillation; as air masses cool and moisture condenses, heavier 18O is enriched in the liquid phase and removed by precipitation. With progressive rainout along a moisture pathway, precipitation becomes gradually more depleted (Dansgaard, 1964). Within this framework, δ18Oprecip is controlled by two variables: temperature and the amount of precipitation along a moisture pathway. The temperature effect stems from the cooling required for progressive rainout during Rayleigh distillation (Dansgaard, 1964; Rozanski et al., 1993). The temperature-δ18O impact is dominant at mid to high latitudes, whilst observations suggest that changes in upstream and local precipitation dominate changes in the δ18Oprecip signal at tropical latitudes. The negative relationship between local precipitation and δ18Oprecip, often referred to as the “amount effect” (Bailey et al., 2018; Dansgaard, 1964), results from the re-evaporation and diffusive exchange between precipitation and water vapour during deep convective precipitation (Risi et al., 2008). However, Rayleigh distillation is complicated by changes in atmospheric circulation and moisture recycling. Changes in the area from which the moisture is sourced will modify δ18Oprecip because the initial δ18O values differ between sources (Cole et al., 1999; Friedman et al., 2002), whilst changes in the moisture transport pathway and/or distance between the source and cave site can result in differing degrees of fractionation associated with condensation and evaporation (Aggarwal et al., 2012; Bailey et al., 2018). The isotopic composition of atmospheric water vapour may also be modified by precipitation recycling over land, as evapotranspiration returns moisture from precipitation back to the atmosphere, thereby minimising the δ18Oprecip/distance gradient along an advection path that occurs with Rayleigh distillation (Gat, 1996; Salati et al., 1979).

Speleothem δ18O records from monsoon regions show multi-millennial variability that has been interpreted as documenting the waxing and waning of the monsoons in response to changes in summer insolation, often interpreted predominantly as a change in the absolute amount of precipitation (Cheng et al., 2013; Fleitmann et al., 2003) or a change in the ratio of more negative δ18O summer precipitation to less negative δ18O winter precipitation (Dong et al., 2010; Wang et al., 2001). However, the multiplicity of processes that influence δ18O before its incorporation in the speleothem makes it difficult to attribute the climatic causes of changes in individual speleothem records unambiguously. In the East Asian monsoon, for example, speleothem δ18O records have been interpreted as a summer monsoon signal, manifested as a change in the amount of water vapour removed along the moisture trajectory (Yuan et al., 2004), and/or as a change in the contribution of summer precipitation to annual totals (Cheng et al., 2006, 2009, 2016; Wang et al., 2001) based on the relationship between modern δ18Oprecip and climate. Other interpretations of Chinese monsoon δ18Ospel have included rainfall source changes (Tan, 2009, 2011, 2014) or local rainfall changes (Tan et al., 2015). Maher (2008) interpreted δ18Ospel as reflecting changes in moisture source area, based on differences between δ18Ospel and loess/palaeosol records of rainfall and the strong correlation between East Asian and Indian monsoon speleothems. Maher and Thompson (2012) used a mass balance approach to show that the changes in precipitation (either local or upstream) or rainfall seasonality required to reproduce δ18Ospel trends would be unreasonably large. Thus, they argued that changes in moisture source were required to explain shifts in δ18O both on glacial–interglacial timescales and during interglacials. Overall, there are several plausible climate mechanisms that could contribute to δ18Ospel on multi-millennial timescales. East Asian monsoon speleothem records are often interpreted as a combination of several of these processes (Cheng et al., 2016; Dykoski et al., 2005) which generally represent monsoon intensity (Cheng et al., 2019). There are also multiple interpretations of the causes of δ18Ospel variability in other monsoon regions. In the Indian monsoon region, speleothem δ18O records are interpreted primarily as an amount effect signal (Berkelhammer et al., 2010; Fleitmann et al., 2004), supported by modern δ18Oprecip and climate observations (e.g. Battacharya et al., 2003). However, other studies have suggested that δ18Oprecip changes in this region are driven primarily by large-scale changes in monsoon circulation; hence, Indian monsoon δ18Ospel should be interpreted as a moisture source and/or trajectory signal (Breitenbach et al., 2010; Sinha et al., 2015). In the Indonesian–Australian monsoon region, δ18Ospel variability has been interpreted as a precipitation amount signal (Carolin et al., 2016; Krause et al., 2019) or a precipitation seasonality signal (Ayliffe et al., 2013; Griffiths et al., 2009), based on modern δ18Oprecip and climate observations (Cobb et al., 2007; Moerman et al., 2013), and/or as a moisture source and/or trajectory signal (Griffiths et al., 2009; Wurtzel et al., 2018). South American speleothem records have been interpreted as records of monsoon intensity, due to changes in the amount of precipitation over the region (Cruz et al., 2005; Wang et al., 2006; Cheng et al., 2013), changes in the degree of upstream precipitation and evapotranspiration (Cheng et al., 2013), or changes in the ratio of precipitation sourced from the low-level jet versus the Atlantic (Cruz et al., 2006; Wang et al., 2006).

These interpretations generally rely on modern δ18Oprecip and climate observations, which may not have remained constant through time. The sources of δ18O variability can also be explored using isotope-enabled climate models (e.g. Hu et al., 2019), which incorporate known isotope effects and, therefore, provide plausible explanations for δ18Ospel trends. Modelling studies suggest that changes in East Asian monsoon δ18Oprecip, during Heinrich events (Lewis et al., 2010; Pausata et al., 2010) and on an orbital timescale (Battisti et al., 2014; LeGrande and Schmidt, 2009), do not reflect local rainfall variability but instead reflect changes in the δ18O of vapour delivered to the region. Variability in the δ18O of vapour delivered to East Asia on orbital timescales has been diagnosed as being due to changes in precipitation upstream of the region (Battisti et al., 2014), changes in moisture source location (Hu et al., 2019; Tabor et al., 2018) or changes in the strength of monsoon winds (LeGrande and Schmidt, 2009; Liu et al., 2014). δ18Oprecip variability in the East Asian monsoon during Heinrich events has also been attributed to non-local isotope fractionation (Lewis et al., 2010; Pausata et al., 2011). Modelling results suggest that changes in precipitation amount are the predominant source of δ18O variability in the Indian monsoon during the Holocene (LeGrande and Schmidt, 2009) and in the glacial (Lewis et al., 2010) as well as in the South American and Indonesian–Australian regions during Heinrich events (Lewis et al., 2010) and the Last Interglacial (Sjolte and Hoffman, 2014).

In this study, we combine speleothem δ18O records of the Speleothem Isotopes Synthesis and AnaLysis (SISAL) database with isotope-enabled palaeoclimate simulations from two climate models to investigate the plausible mechanisms driving changes in δ18O in monsoon regions through the Holocene (last 11 700 years) and between the mid-Holocene, the peak of the Last Interglacial and the Last Glacial Maximum. We compare δ18Ospel signals across geographically separated cave sites to extract a regional signal, thereby minimising the influence of karst and in-cave processes, such as the mixing of groundwaters from different precipitation events or changes in cave ventilation, that can be important for the δ18Ospel of individual records. We use principal coordinate analysis (PCoA) to identify regions with geographically coherent δ18Ospel records and then examine how these regions behave on glacial–interglacial timescales and through the Holocene. We use isotope-enabled model simulations to investigate the potential causes of δ18Ospel variability in regions where the models reproduce the large-scale δ18O changes shown by observations. We exploit the fact that models produce internally physically consistent changes to explore potential and plausible causes of the trends observed in speleothem records across specific monsoon regions, using multiple regression analysis.

2 Methods

2.1 Speleothem oxygen isotope data

Speleothem δ18O records were obtained from the SISAL (Speleothem Isotopes Synthesis and AnaLysis) database (Atsawawaranunt et al., 2018; Comas-Bru et al., 2020a, b). Records were selected based on the following criteria:

  • they are located in monsoon regions, between 35 S and 40 N;

  • the mineralogy is known and does not vary (i.e. between calcite and aragonite) through time, because oxygen isotope fractionation during speleothem precipitation is different for calcite and aragonite;

  • for the analysis of mid-Holocene (MH), Last Glacial Maximum (LGM) and Marine Isotope Stage 5e during the Last Interglacial (LIG) δ18O signals, the records contain samples within at least one of these time periods, defined as 6000±500 years BP for the MH, 21 000±1000 years BP for the LGM and 125 000±1000 years BP for the LIG, where BP (before present) is 1950 CE;

  • for the PCoA, the records have a temporal coverage of at least 4000 years in the Holocene;

  • for Holocene trend analyses, speleothems have a record of the period from 7000 to 3000 years BP;

  • they are the most recent update of the record from a site available in version 2 of the SISAL database.

This resulted in the selection of 125 records from 44 sites for the PCoA analysis, 64 records from 38 sites for the analysis of MH, LGM and LIG signals, and 79 records from 40 sites for the Holocene trend analysis (Fig. 1). Although the SISALv2 database contains multiple age models for some sites, we use the published age models given by the original authors for all records.

Figure 1Spatial distribution of speleothem records used is this study. Colours indicate the sites used in principal coordinates analysis (PCoA) and redundancy analysis (RDA) to separate monsoon regions, and sites not used in PCoA and RDA but used in subsequent analyses. The individual regional monsoons are shown by boxes: CAM denotes the Central American monsoon (latitude from 10 to 33, longitude from −115 to −58), SW-SAM denotes the southwestern South American monsoon (latitude from −10 to 0, longitude from −80 to −64; latitude from −30 to −10, longitude −68 to −40), NE-SAM denotes the northeastern South American monsoon (latitude from −10 to 0, longitude: −60 to −30), SAfM denotes the southern African monsoon (latitude from −30 to −17, longitude from 10 to 40), ISM denotes the Indian summer monsoon (latitude from 11 to 32, longitude from 50 to 95), EAM denotes the East Asian monsoon (latitude from 20 to 39, longitude from 100 to 125) and IAM denotes the Indonesian–Australian monsoon (latitude from −24 to 5, longitude from 95 to 135). Source region limits used in the multiple linear regression analysis are also shown. The background carbonate lithology is from the WOrld Karst Aquifer Mapping (WOKAM) project (Goldschneider et al., 2020).

2.2 Climate model simulations

There are relatively few palaeoclimate simulations made with models that incorporate oxygen isotope tracers, and the available simulations do not necessarily focus on the same periods or use the same modelling protocols. Here, we use simulations of opportunity from two isotope-enabled climate models: ECHAM5 (version 5 of the European Centre for medium-range weather forecasting model in HAMburg) and GISS ModelE-R (Goddard Institute for Space Studies Model version E-R). The ECHAM5 simulations provide an opportunity to examine large-scale changes between glacial and interglacial states, using simulations of the MH, LGM and LIG. The GISS ModelE-R ocean–atmosphere coupled general circulation model was used to investigate the evolution of δ18O evolution during the Holocene, using eight time slice (9, 6, 5, 4, 3, 2, 1 and 0 ka) experiments. Although simulations of the MH 6 ka time slice are available with both models, there are differences in the protocol used for the two experiments which preclude direct comparison of the simulations.

The ECHAM5-wiso MH experiment (Wackerbarth et al., 2012; Werner, 2019) was forced by orbital parameters (based on Berger and Loutre, 1991) and greenhouse gas (GHG) concentrations (CO2=280 ppm, CH4=650 ppb, N2O=270 ppb) appropriate to 6 ka. Changes in sea surface temperature (SST) and sea ice were derived from a transient Holocene simulation (Varma et al., 2012). The control simulation for the MH experiment was an ECHAM-wiso simulation of the period from 1956 to 1999 (Langebroek et al., 2011), using observed SSTs and sea ice cover. This control experiment was forced by SSTs and sea ice only, with atmospheric circulation free to evolve. The ECHAM5-wiso LGM experiment (Werner, 2019; Werner et al., 2018) was forced by orbital parameters (Berger and Loutre, 1991), GHG concentrations (CO2=185 ppm, CH4=350 ppb, N2O=200 ppb), land–sea distribution and ice sheet height and extent appropriate to 21 ka; SST and sea ice cover were prescribed from the GLAMAP dataset (Schäfer-Neth and Paul, 2003). Sea surface water and sea ice δ18O were uniformly enriched by 1 ‰ at the start of the experiment. The control simulation for the LGM experiment used present-day conditions, including orbital parameters and GHG concentrations set to modern values, and SSTs and sea ice cover from the last 20 years (1979–1999). Both the MH and LGM simulations were run at T106 horizontal grid resolution, approximately 1.1×1.1 at the Equator. Comparison of the MH and LGM simulations with speleothem data globally (Figs. S1 and S2 in the Supplement; Comas-Bru et al., 2019) show that the ECHAM model reproduces the broadscale spatial gradients and the sign of isotopic changes at the majority of cave sites (MH: 72 %; LGM: 76 %). However, the changes compared with present are generally more muted in the simulations than shown by the speleothem records.

The LIG experiment (Gierz et al., 2017b, a) was run using the ECHAM5/MPI-OM Earth system model, with stable water isotope diagnostics included in the ECHAM5 atmosphere model (Werner et al., 2011), the dynamic vegetation model JSBACH (Haese et al., 2012) and the MPI-OM ocean–sea-ice module (Xu et al., 2012). This simulation was run at a T31L19 horizontal grid resolution, approximately 3.75×3.75. The LIG simulation was forced by orbital parameters derived from Berger and Loutre (1991) and GHG concentrations (CO2=276 ppm, CH4=640 ppb, N2O=263 ppb) appropriate to 125 ka, but it was assumed that the ice sheet configuration and land–sea geography remain unchanged from modern; therefore no change was made to the isotopic composition of sea water. The LIG simulation is compared to a pre-industrial (PI) control with appropriate insolation, GHG and ice sheet forcing for 1850 CE. The sign of simulated isotopic changes in the LIG is in good agreement with ice core records from Antarctica and Greenland and speleothem records from Europe, the Middle East and China (Gierz et al., 2017b), although, as with the MH and LGM, the observed changes tend to be larger than the simulated changes (Fig. S3).

There are GISS ModelE-R (LeGrande and Schmidt, 2009) simulations for eight time slices during the Holocene (9, 6, 5, 4, 3, 2, 1 and 0 ka). The 0 ka experiment is considered as the pre-industrial control (ca. 1880 CE). Orbital parameters were based on Berger and Loutre (1991), and GHG concentrations were adjusted based on ice core reconstructions (Brook et al., 2000; Indermühle et al., 1999; Sowers, 2003) for each time slice. A remnant Laurentide ice sheet was included in the 9 ka simulation, following Licciardi et al. (1998), and the corresponding adjustment was made to mean ocean salinity and ocean water δ18O to account for this (Carlson et al., 2008). The ice sheet in all the other experiments was specified to be the same as modern; therefore, no adjustment was necessary. The simulations were run using the M20 version of GISS ModelE-R, which has a horizontal resolution of 4×5. Each experiment was run for 500 years, and we use the last 100 simulated years for the analyses. Comparison of the simulated trends in δ18O show good agreement with Greenland ice core records, marine records from the tropical Pacific and Chinese speleothem records (LeGrande and Schmidt, 2009). However, as is the case with the ECHAM simulations, the model tends to produce changes that are less extreme than those shown by the observations (Figs. S4, S5, S6).

2.3 Principal coordinate analysis and redundancy analysis

We used PCoA to identify regionally coherent patterns in the speleothem δ18O records for the Holocene. PCoA is a multivariate ordination technique that uses a distance matrix to represent inter-object (dis)similarity in reduced space (Gower, 1966; Legendre and Legendre, 1998). Speleothem records from individual sites are often discontinuous; missing data are problematic for many ordination techniques. PCoA is more robust to missing data than other methods (Kärkkäinen and Saarela, 2015; Rohlf, 1972). We used a correlation matrix of speleothem records as the (dis)similarity measure. The temporal resolution of speleothem records was first standardised by calculating a running average mean with non-overlapping 500-year windows. This procedure produces a single composite record when there are several records for a given site. PCoA results were displayed as a biplot, where sites ordinated close to one another (i.e. with similar PCoA scores) show similar Holocene trends, and sites ordinated far apart have dissimilar trends. We used the “broken stick” model (Bennett, 1996) to identify which PCoA axes were significant. We used redundancy analysis (RDA: Legendre and Legendre, 1998; Rao, 1964) with latitude and longitude as predictor variables to identify if PCoA (dis)similarities were related to geographical location, and principal components analysis (PCA) to identify the main patterns of variation. As these explanatory variables are not dimensionally homogeneous, they were centred on their means and standardised to allow direct comparison of the gradients. PCoA and RDA analyses were carried out using the “vegan” package in R (Oksanen et al., 2019).

2.4 Glacial–interglacial changes in δ18O

We examined shifts in δ18Ospel observations and in annual precipitation-weighted mean δ18Oprecip from ECHAM-wiso in regions influenced by the monsoon, between the MH, LGM and LIG. Values are given as anomalies with respect to the present day for speleothems or the control simulation experiment for model outputs. Comas-Bru et al. (2019) have shown that differences in speleothem δ18O data between the 20th century and the pre-industrial period (i.e. 1850±15 CE) are within the temporal and measurement uncertainties of the data; thus, the use of different reference periods (i.e. PI for the ECHAM LIG experiment, 20th century for ECHAM MH, LGM experiments) should have little effect on our analyses. We used mean site δ18Ospel values for each period for the regions identified in the PCoA analysis. Where there are multiple speleothem δ18O records for a site in a time period, they were averaged to calculate mean δ18Ospel. Three sites above 3500 m were excluded from the calculation of the means because high-elevation sites have more negative δ18O values than their low-elevation counterparts, and their inclusion would distort the regional estimates.

There are relatively few speleothems covering both the present day and the period of interest (i.e. MH, LGM or LIG), precluding the calculation of δ18Ospel anomalies from the speleothem data. Therefore, we calculated anomalies with respect to the modern period (1960–2017 CE) using the Online Isotopes in Precipitation Calculator (OIPC: Bowen, 2018; Bowen and Revenaugh, 2003), a global gridded dataset of interpolated mean annual precipitation-weighted δ18Oprecip data, as reference. This dataset combines data from 348 stations from the Global Network of Isotopes in Precipitation (IAEA/WMO, 2018), covering part or all of the 1960–2014 period, and other records available at the Waterisotopes Database (Waterisotopes Database, 2017).

OIPC δ18Oprecip was converted to its speleothem equivalent assuming that (i) precipitation-weighted mean annual δ18Oprecip is equivalent to mean annual drip-water δ18O (Yonge et al., 1985) and that (ii) precipitation of calcite is consistent with the empirical speleothem-based kinetic fractionation factor of Tremaine et al. (2011) and precipitation of aragonite follows the fractionation factor from Grossman and Ku (1986), as formulated by Lachniet (2015):

(1) δ 18 O calcite _ SMOW = w δ 18 O precip _ SMOW + 16.1 1000 T - 24.6 ( T in K ) ,

(2) δ 18 O aragonite _ SMOW = w δ 18 O precip _ SMOW + 18.34 1000 T - 31.954 ( T in K ) ,

where δ18Ocalcite_SMOW and δ18Oaragonite_SMOW are the respective speleothem isotopic composition for calcite and aragonite speleothems with reference to the VSMOW (Vienna Standard Mean Ocean Water) standard (in per mille), wδ18Oprecip is the OIPC precipitation-weighted annual mean isotopic composition of precipitation with respect to the VSMOW standard and T is the mean annual cave temperature (in kelvin). We used the long-term (1960–2016) mean annual surface air temperature from the CRU-TS4.01 database (Harris and Jones, 2017; Harris et al., 2020) at each site as a surrogate for mean annual cave air temperature. The resolution of the gridded data means that wδ18Oprecip_SMOW and T may be the same for nearby sites.

We use the VSMOW to VPDB (Vienna Pee Dee Belemnite) conversion from Coplen et al. (1983), which is independent of speleothem mineralogy:

(3) δ 18 O PDB = 0.97001 δ 18 O SMOW - 29.29 ,

where δ18OPDB is relative to the VPDB standard, and δ18OSMOW is relative to VSMOW standard.

Average uncertainties in the speleothem age–depth models are ∼50 years during the Holocene. This interval is smaller than the time windows used in this analysis; therefore, the age uncertainty is expected to have a negligible impact on the results. We investigated the influence of age uncertainties on the LGM and LIG δ18Ospel anomalies by examining the impact of using different window widths (±500, ±700, ±1000, ±2000 years) on the regional mean δ18Ospel anomalies.

We used anomalies of wδ18Oprecip, mean annual surface air temperature (MAT) and mean annual precipitation (MAP) from the ECHAM5-wiso simulations to investigate the changes in δ18Ospel between the MH, LGM and LIG as well as their association with changes in climate. Values were calculated from land grid cells (>50 % land) ±3 around each speleothem site. This distance was chosen with reference to the coarsest-resolution simulation (LIG, ca. 3.75×3.75). Gridded values of MAT and MAP were weighted by the proportion of each grid cell that lies within ±3 of the site, and linear distance-weighted means were calculated for each site and time slice. We only considered regions with at least one speleothem record for each of the three time periods, although these were not required to be the same sites, and where the observed shifts in δ18Ospel were in the same direction and of a similar magnitude to the simulated wδ18Oprecip.

2.5 Holocene and Last Interglacial regional trends

Regional speleothem δ18O changes through the Holocene were examined by creating composite time series for each region identified in the PCoA analysis with at least four Holocene records (>5000 years long). Regional composites were constructed using a four-step procedure, modified from Marlon et al. (2008):

  1. The δ18O data for individual speleothems were transformed to z scores, so all records have a standardised mean and variance:

    (4)z score i = δ 18 O i - δ 18 O ( base period ) s δ 18 O ( base period ) ,

    where δ18O is the mean, and sδ18O is the standard deviation of δ18O for a common base period. A base period of 7000 to 3000 years BP was chosen to maximise the number of records included in each composite.

  2. The standardised data for a site were resampled by applying a 100-year non-overlapping running mean with the first bin centred at 50 years BP, in order to create a single site time series while ensuring that highly resolved records do not dominate the regional composite.

  3. Each regional composite was constructed using locally weighted regression (Cleveland and Devlin, 1988) with a window width of 3000 years and fixed target points in time.

  4. Confidence intervals (5th and 95th percentiles) for each composite were generated by bootstrap resampling by site over 1000 iterations.

There are too few sites to construct regional composites for the peak of the LIG (Marine Isotope Stage 5e); thus, the trends in δ18Ospel were examined using records from individual sites covering the period from 130 to 116 ka BP.

We calculated Holocene regional composites from annual precipitation-weighted mean δ18Oprecip anomalies simulated by the GISS model. Simulated δ18Oprecip trends were calculated using linear distance-weighted mean δ18Oprecip values from land grid cells (>50 % land) within ±4 around each site. This distance was determined by the grid resolution of the model. Regional composites were then produced using bootstrap resampling in the same way as for the speleothem data. The simulated anomalies are relative to the control run rather than the specified base period used for the speleothem-based composites, so absolute values of simulated and observed Holocene trends are expected to differ. Preliminary analyses showed that neither the mean values nor trends in δ18Oprecip were substantially different if the sampled area was reduced to match the sampling used for the ECHAM-based box plot analysis, or was increased to encompass the larger regions shown in Fig. 1 and used in the multiple regression analysis.

2.6 Multiple regression analysis

We investigate the underlying relationships between regional δ18Oprecip (and by extension δ18Ospel and monsoon climate through the Holocene using multiple linear regression (MLR). We use annual precipitation-weighted mean δ18Oprecip anomalies and climate variables from GISS ModelE-R. Climate variables were chosen to represent the four potential large-scale drivers of regional changes in the speleothem δ18O records. Specifically, we use changes in mean precipitation and precipitation recycling over the monsoon regions, and changes in mean surface air temperature and surface wind direction over the moisture source regions. Whereas the influence of changes in precipitation, recycling and temperature are relatively direct measures, the change in surface wind direction over the moisture source region is used as an index of potential changes in the moisture source region and transport pathway. The boundaries of each monsoon region (Fig. 1) were defined to include all the speleothem sites used to construct the Holocene δ18Ospel composites. Moisture source area limits (Fig. 1) were defined based on moisture tracking studies (Bin et al., 2013; Breitenbach et al., 2010; D'Abreton and Tyson, 1996; Drumond et al., 2008, 2010; Durán-Quesada et al., 2010; Kennett et al., 2012; Nivet et al., 2018; Wurtzel et al., 2018) and GISS simulated summer surface winds. All climate variables were extracted for the summer months, defined as May to September (MJJAS) for Northern Hemisphere regions and November to March (NDJFM) for Southern Hemisphere regions (Wang and Ding, 2008) on the basis that these regions are dominated by summer season precipitation (Fig. S7). Only grid cells with >50 % land were used to extract variables over monsoon regions, and only grid cells with <50 % land were used to extract variables over moisture source regions. The inputs to the MLR for each time interval were calculated as anomalies from the control run.

Precipitation recycling was calculated as the ratio of locally sourced precipitation versus total precipitation. Although the GISS ModelE-R mid-Holocene experiment explicitly estimates recycling using vapour source distribution tracers (Lewis et al., 2014), this was not done for all the Holocene time slice simulations. Therefore, we calculate a precipitation recycling index (RI), following Brubaker et al. (1993):

(5) RI = P R P = E 2 Q H + E ,

where locally sourced (recycled) precipitation (PR) is estimated using total evaporation over a region (E), and total precipitation (P) is estimated as the sum of total evaporation and net incoming moisture flux integrated across the boundaries of the region (QH). Therefore, RI expresses the change in the contribution of local, recycled precipitation independently of any overall change in precipitation amount.

We incorporate mean meteorological variables and δ18Oprecip for all Holocene time slices (1 to 9 ka) and all monsoon regions into the MLR model. Thus, the relationships constrained by the overall (global) MLR model represent the combined response across all monsoon regions. We use the pseudo-R2 to determine the goodness of fit for the global MLR model and t values (the regression coefficient divided by its standard error) to determine the strength of each relationship. Partial residual plots were used to show the relationship between each predictor variable and δ18Oprecip when the effects of the other variables are held constant.

All statistical analyses were performed in R (R Core Team, 2019), and plots were generated using ggplot (Wickham, 2016).

3 Results

3.1 Principal coordinate analysis and redundancy analysis

PCoA shows the (dis)similarity of Holocene δ18Ospel evolution across individual records and, thus, allows an objective regionalisation of these records. The first two PCoA axes are significant, according to the broken stick test, and account for 65 % and 20 % of δ18Ospel variability respectively (Table 1). The PCoA scores differentiate records geographically (Fig. 2a): Southern Hemisphere monsoon regions such as the southwestern South American monsoon (SW-SAM) and South African monsoon (SAfM) are characterised by low PCoA1 scores, whereas Northern Hemisphere monsoons such as the Indian summer monsoon (ISM) and the East Asian monsoon (EAM) are characterised by higher PCoA1 scores. This indicates that regions can be differentiated based on their temporal evolution as captured by the first PCoA axis. Most Southern Hemisphere regions also have lower PCoA2 scores, although this is not consistent over time. Speleothem records from the Central American monsoon (CAM) and the Indonesian–Australian monsoon (IAM) have PCoA scores that are an intermediate between the Northern Hemisphere and Southern Hemisphere regions. PCoA clearly separates the South American records into a northeastern region (NE-SAM), with scores similar to other Northern Hemisphere monsoon regions, and a southwestern region (SW-SAM), with scores similar to other Southern Hemisphere regions. The RDA supports a geographical control on the (dis)similarity of speleothem δ18O records over the Holocene (Fig. 2b). RDA1 explains 37 % of the variability and is significantly correlated with both latitude and longitude (Table 2).

Table 1Results of the principal coordinates analysis (PCoA). Significant axes, as determined by the broken stick method (Bennett, 1996), are shown in bold.

Download Print Version | Download XLSX

Figure 2Results of the principal coordinate analysis (PCoA) and redundancy analysis (RDA). (a) PCoA biplot showing the loadings of each site on the first two axes, which represent 85 % of the total variance. Shapes indicate the Holocene coverage of each site, where sites with a coverage ≥8000 years represent most or all of the Holocene (Hol). Sites with a temporal coverage of <8000 years are coded to show whether they represent the early Holocene to mid-Holocene (EH to MH; record midpoint >8000 years BP), the mid-Holocene (MH; record midpoint between 8000 and 5000 years BP) or the mid-Holocene to late Holocene (LH to MH; midpoint <5000 years BP). (b) RDA triplot, where the response variables are the PCoa1 and PCoA2 axes explained by latitude and longitude. The direction of the PCoA axes have been fixed so that they align with the explanatory variables.


Table 2Results of the redundancy analysis (RDA). Variables that are significantly correlated (P <0.01) with the RDA axes are shown in bold.

Download Print Version | Download XLSX

3.2 Regional interglacial–glacial differences

To investigate the causes of shifts in δ18O between the MH, LGM and LIG, we compare simulated and observed regional δ18O signals during these periods with shifts in climate variables (precipitation and temperature). Only the ISM, EAM and IAM regions have sufficient speleothem data (i.e. at least one record from every time period) to allow comparisons across the MH, LGM and LIG (Fig. 3) and have similar shifts in observed δ18Ospel and simulated δ18Oprecip. The regional mean δ18Ospel anomalies calculated for different time windows (±500, ±700, ±1000, ±2000 years) vary by less than 0.35 ‰ for the LGM (ISM: <0.16 ‰; EAM: <0.35 ‰; IAM: <0.22 ‰) and 0.48 ‰ for the LIG (ISM: <0.16 ‰; EAM: <0.48 ‰; IAM: <0.11 ‰), indicating that age uncertainties have a minimal impact on these mean values. The most positive δ18Ospel anomalies in all three regions occur at the LGM, with more negative anomalies for the MH and LIG. The simulated δ18Oprecip anomalies show a similar pattern, with more positive anomalies during the LGM than during the MH or the LIG. The amplitude of this pattern is also similar between δ18Oprecip and δ18Ospel, when the observations are converted to their drip water equivalent (Fig. S8). The differences in regional δ18Ospel anomalies between MH and LIG differ across the three regions. In both the ISM and the EAM, differences in δ18Ospel values between the MH and LIG are small (Fig. 3a, b), although ISM LIG δ18Ospel values are slightly more negative than MH values. In the IAM, MH values are less negative than the LIG (Fig. 3c). However, there are only a limited number of speleothem records from the ISM and IAM during the LIG, so the apparent differences between the two intervals in these regions may not be meaningful. Glacial–interglacial shifts are also seen in simulated temperature and precipitation, with warmer and wetter conditions during interglacials and cooler and drier conditions during the LGM in all three regions. Differences in simulated precipitation between the MH and the LIG could help explain the differences between δ18Ospel in the ISM and IAM, as the LIG is wetter than the MH in the ISM and drier than the MH in the IAM. However, the LIG is also drier than the MH in the EAM, a feature that appears inconsistent with the lack of differentiation between the δ18O signals in this region.

Figure 3Speleothem δ18O anomalies compared to anomalies of δ18Oprecip, precipitation and temperature from the ECHAM simulations for the (a) East Asian (EAM), (b) Indian (ISM) and (c) Indonesian–Australian (IAM) monsoons. The boxes show the median value (line) and the interquartile range, and the whiskers show the minimum and maximum values, with outliers represented by grey dots. Note that the isotope axes are reversed, so that the most negative anomalies are at the top of the plot, to be consistent with the assumed relationship with the direction of change in precipitation and temperature.


3.3 Regional-scale interglacial δ18O evolution

There are four monsoon regions with sufficient data to examine regional Holocene δ18O trends: EAM, ISM, SW-SAM and IAM (Fig. 4). The IAM region has the fewest records (n=7), whereas the EAM has the largest number of records (n=14). The regional composites are expressed as z scores, i.e. changes with respect to the mean and variance of δ18O for the base period (3000–7000 years BP). The confidence intervals on the regional composites are small for all regions, except SW-SAM in the early Holocene. The EAM and ISM regions (Fig. 4a–e) show the most positive δ18Ospel z scores around 12 ka followed by a rapid decrease towards their most negative values at ∼9.5 and ∼9 ka respectively. The δ18Ospel z scores in the EAM are relatively constant from 9.5 to ∼7 ka, whereas this plateau is present but less marked in the ISM. There is a gradual trend towards more positive δ18Ospel z scores towards the present in both regions thereafter. The SW-SAM records (Fig. 4i) have their most positive δ18Ospel z scores in the early Holocene with a gradual trend to more negative scores towards the present. By contrast, the IAM z scores (Fig. 4g) are most positive at 12 ka, gradually decrease until ca. 5 ka and are relatively flat thereafter.

Figure 4Evolution of regional speleothem δ18O signals through the Holocene compared to δ18Oprecip simulated by the GISS model. Panels (a–e) show Northern Hemisphere monsoons (EAM denotes the East Asian monsoon; ISM denotes the Indian summer monsoon) and summer (May through September) insolation at 30 N (Berger, 1978). Panels (f–j) show Southern Hemisphere monsoons (SW-SAM denotes the southwestern South American monsoon; IAM denotes the Indonesian–Australian monsoon) and summer (November through March) insolation for 20 S (Berger, 1978). The speleothem δ18O changes are expressed as z scores, with a smoothed loess fit (3000-year window), and confidence intervals obtained by bootstrapping by site. δ18Oprecip values are expressed as anomalies from the pre-industrial control simulation. Note that the isotope axes are reversed, so that the most negative anomalies are at the top of the plot, to be consistent with the assumed relationship with the changes in insolation.


There are insufficient data to create composite curves for the LIG, but individual records from the four regions (Fig. 5) show similar features to the Holocene trends. Records from the ISM and EAM (Fig. 5 left), for example, are characterised by an initial sharp decrease in δ18Ospel values of about 4 ‰ between 130 and 129 ka, and most of the records (Dykoski et al., 2005; Kathayat et al., 2016; Wang et al., 2008) then show little variability for several thousand years. Despite the fact that the Tianmen record (Cai et al., 2010, 2012) shows considerable variability between 123 and 127 ka, there is nevertheless a similar plateau in the average observed value before the rapid change to less negative values after 127 ka. Similar to the Holocene, the SW-SAM record (Cheng et al., 2013) shows increasingly negative δ18Ospel values through the LIG. The trend shown for Whiterock Cave (Carolin et al., 2016) also displays similar features to the IAM Holocene composite, with a gradual trend towards more negative values initially and a relatively complacent curve towards the end of the interglacial (Fig. 5 right).

Figure 5Comparison of changes in summer insolation and δ18Ospel through the peak of the Last Interglacial (Marine Isotope Stage 5e) from the (b, c) East Asian monsoon (EAM), (d, e) Indian summer monsoon (ISM), (g) southwestern South American monsoon (SW-SAM) and (h) Indonesian–Australian monsoon (IAM) regions. The U–Th dates and uncertainties are shown for each record. The summer insolation curves (Berger, 1978) are for May through September at 30 N in the Northern Hemisphere (a) and for November through March for 20 S in the Southern Hemisphere (f). Note that the isotope axes are reversed, so that the most negative anomalies are at the top of the plot, to be consistent with the assumed relationship with the changes in insolation. The LIG (Marine Isotope Stage 5e) time slice used in the analysis in Sect. 2.4 is shown by the dark-grey bar.


3.4 Multiple regression analysis of Holocene δ18Oprecip

The MLR analyses of simulated δ18Oprecip trends identify the impact of an individual climate variable on δ18Oprecip in the absence of changes in other variables. The global MLR model includes the Holocene (1–9 ka) δ18Oprecip trends combined across all monsoon regions (CAM, ISM, EAM, SW-SAM, NE-SAM, SAfM, IAM). This global monsoon MLR model has a pseudo-R2 of 0.80 and shows statistically significant relationships between the anomalies in δ18Oprecip and anomalies in regional precipitation, temperature and surface wind direction (Table 3). The partial residual plots (Fig. 6) show there is a strong negative relationship with regional precipitation (t value of −8.75) and a strong positive relationship (t value of 8.03) with surface wind direction over the moisture source region, an index of changes in either source area or moisture pathway. This indicates that increases in regional precipitation alone will lead to a decrease in δ18O, whereas changes in source area and/or moisture pathway, in the absence of changes in other variables, will lead to a significant change in δ18O. The relationship with temperature over the moisture source region is weaker, but positive (t value of 2.05), i.e. an increase in temperature over the moisture source region will lead to an increase in δ18O if there are no changes in other climate variables. Precipitation recycling is not significant in this global analysis. The exact choice of source region has a negligible impact on the model – for example, expanding the ISM source region to include the Bay of Bengal does not change the outcome of this analysis (Fig. S9, Table S1).

Figure 6Partial residual plots from the multiple linear regression analysis, showing the relationship between anomalies in simulated δ18Oprecip and the four predictor variables, after considering the fitted partial effects of all the other predictors. The simulated δ18Oprecip values are anomalies relative to the pre-industrial control simulation and are annual values weighted by precipitation amount. The predictor variables are precipitation in the delineated monsoon region (mm d−1), temperature in the source region (C), surface wind direction over the source region () as an index of potential changes in source region and the ratio of precipitation recycling to total precipitation over the monsoon region (RI, unitless). The predictor variables are summer mean values, representing the summer monsoon, where summer is defined as May to September for the Northern Hemisphere monsoons and November to March for the Southern Hemisphere monsoons.


Table 3Results of the multiple linear regression analysis. Significant relationships (P <0.01) are shown in bold.

Download Print Version | Download XLSX

There are too few data points to make regressions for individual monsoon regions, but the distribution of data points for each region in the partial residual plots (Fig. 6) is indicative of the degree of conformity to the global MLR model (representing the combined response across all monsoon regions). Data points from the ISM, SW-SAM, IAM and SAfM are well aligned with the overall relationship with regional precipitation (Fig. 6a), indicating that precipitation is an important control on changes in δ18Oprecip in these regions. The NE-SAM, EAM and CAM values deviate somewhat from the overall relationship and, although there are relatively few points, this suggests that changes in precipitation are a less important influence on δ18Oprecip changes in these regions. The impact of temperature changes (Fig. 6b) in the ISM, EAM and SW-SAM is broadly consistent with the overall relationship. The slope of the relationship with temperature is negative for the IAM and NE-SAM, and as this is physically implausible, it suggests that some factor not currently included in the MLR is influencing these records. However, the inconsistencies between the regional signals helps to explain why the global relationship between anomalies in temperature and δ18Oprecip is weak (Fig. 6b) and probably reflects the fact that tropical temperature changes during the Holocene are small. Data points from the EAM, ISM and IAM are well aligned with the overall relationship between changes in δ18Oprecip and changes in wind direction (Fig. 6c), indicating that changes in source area or moisture pathway are an important control on changes in δ18Oprecip in these regions. However, values for CAM, SW-SAM, NE-SAM and SAfM deviate strongly from the overall relationship. Recycling does not appear to be an important contributor to changes in δ18Oprecip except in SW-SAM and SAfM (Fig. 6d).

4 Discussion

We have shown that it is possible to derive an objective regionalisation of speleothem records based on the PCoA of the oxygen-isotope trends through the Holocene (Fig. 2). This approach separates out regions with a distinctive Northern Hemisphere signal (e.g. ISM, EAM, NE-SAM) from regions with a distinctive Southern Hemisphere signal (e.g. SW-SAM, SAfM), reflecting the fact that the evolution of regional monsoons in each hemisphere follows, to some extent, insolation forcing. It also identifies regions that have an intermediate pattern (e.g. IAM). The robustness of the regionalisation is borne out by the fact that Holocene composite trends from each region have tight confidence intervals (Fig. 4), showing that the signals of individual records across a region show broad similarities. The monsoon regions identified by PCoA are consistent with previous studies (Wang et al., 2014). The tracking of Northern Hemisphere insolation is a recognised feature of monsoon systems in India and China (see reviews by Kaushal et al., 2018; Zhang et al., 2019). The separation of speleothem records from NE-SAM from those in SW-SAM is consistent with the precipitation dipole that exists between northeastern Brazil (Nordeste) and the continental interior (Berbery and Barros, 2002; Boers et al., 2014). The anti-phasing of speleothem records from the two regions during the Holocene has been recognised in previous studies (Cruz et al., 2009; Deininger et al., 2019). The intermediate nature of the records from the maritime continent is consistent with the fact that the Indonesian–Australian (IAM) summer monsoon is influenced by cross-equatorial air flow and, hence, can be influenced by Northern Hemisphere conditions (Trenberth et al., 2000). Palaeoenvironmental records from this region show mixed signals for the Holocene: some have been interpreted as showing enhanced (Beaufort et al., 2010; Mohtadi et al., 2011; Quigley et al., 2010; Wyrwoll and Miller, 2001) and others reduced precipitation (Kuhnt et al., 2015; Steinke et al., 2014) during the early and mid-Holocene. Modelling studies have shown that this region is highly sensitive to SST changes in the Indian Ocean and South China Sea, which in turn reflect changes in the Northern Hemisphere winter monsoons. Although most climate models produce a reduction in precipitation across the IAM during the mid-Holocene in response to orbital forcing, this is less than might be expected in the absence of ocean feedbacks associated with changes in the Indian Ocean (Zhao and Harrison, 2012).

The separation of northern and southern monsoon regions is consistent with the idea that changes in monsoon rainfall are primarily driven by changes in insolation (Ding and Chan, 2005; Kutzbach et al., 2008). Indeed, regional δ18Ospel composites from the EAM, ISM and SW-SAM show a clear relationship with the long-term trends in local summer insolation (Fig. 4). Similar patterns are seen in individual speleothem records from each region confirming that the composite trends are representative. However, the composite trends are not an exact mirror of the insolation signal over the Holocene. For example, the ISM and EAM composites show a more rapid rise during the early Holocene than implied by the insolation forcing. The maximum wet phase in these two regions lasts for ca. 3000 years, again contrasting with the gradual decline in insolation forcing after its peak at ca. 11 ka. Both the rapid increase and the persistence of wet conditions for several thousand years is also observed in other palaeohydrological records across southern and central China, including pollen (Zhao et al., 2009; Li et al., 2018) and peat records (Hong et al., 2003; Zhou et al., 2004). These features are also characteristic of lake records from India (Misra et al., 2019). The lagged response to increasing insolation is thought to be due to the presence of Northern Hemisphere ice sheets in the early Holocene (Zhang et al., 2018). The persistence of wetter conditions through the early and mid-Holocene is thought to reflect the importance of land surface and ocean feedbacks in sustaining regional monsoons (Dallmeyer et al., 2010; Kutzbach et al., 1996; Marzin and Braconnot, 2009; Rachmayani et al., 2015; Zhao and Harrison, 2012). The evolution of regional monsoons during the LIG shows patterns similar to those observed during the Holocene, including the lagged response to insolation and the persistence of wet conditions after peak insolation. This is again consistent with the idea that internal feedbacks play a role in modulating the monsoon response to insolation forcing. We have also shown that there is little difference in the isotopic values between the MH and the LIG in the ISM and EAM regions, which is also observed in individual speleothem records (Kathayat et al., 2016; Wang et al., 2008). The LIG (125 ka) period was characterised by higher summer insolation, higher CO2 concentrations (Otto-Bliesner et al., 2017) and lower ice volumes (Dutton and Lambeck, 2012) than the MH, suggesting that the LIG ISM and EAM monsoons should be stronger than the MH monsoons. The lack of a clear differentiation in the isotope signals between the LIG and MH suggests that other factors play a role in modulating the monsoon response to these forcings and may reflect the importance of global constraints on the externally forced expansion of the tropical circulation (Biasutti et al., 2018).

Global relationships between δ18Oprecip and climate variables (precipitation amount, temperature and surface wind direction; Fig. 6) are consistent with existing studies: a strong relationship with precipitation and a weaker temperature effect have been widely observed at tropical and subtropical latitudes in modern observations (Dansgaard, 1964; Rozanski et al., 1993). The significant global relationship between δ18Oprecip and surface winds supports the idea that changes in moisture source and pathway are also important for explaining δ18O variability over the Holocene. The multiple regression analysis also provides insights into the relative importance of different influences at a regional scale. In the ISM, the results support existing speleothem studies which suggest that changes in precipitation amount (Cai et al., 2015; Fleitmann et al., 2004) and to a lesser extent moisture pathway (Breitenbach et al., 2010) drive δ18Ospel variability. The δ18O variability in the IAM region through the Holocene also appears to be strongly driven by changes in precipitation and moisture pathway, consistent with the interpretation of Wurtzel et al. (2018). Changes in regional precipitation (where the cave sites are located) do not seem to explain the observed changes in δ18Ospel in the EAM during the Holocene, where Holocene δ18Oprecip evolution is largely driven by changes in atmospheric circulation (indexed by changes in surface winds). This is consistent with existing studies that emphasise changes in moisture source and/or pathway rather than local precipitation changes (Maher, 2016; Maher and Thompson, 2012; Tan, 2014; Yang et al., 2014). Speleothem δ18O records in the SW-SAM clearly reflect regional-scale changes in precipitation, consistent with interpretations of individual records (Cruz et al., 2009; Kanner et al., 2013). However, this is a region where changes in precipitation recycling also appear to be important. Based on regional water budget estimates, recycling presently contributes ca. 25 %–35 % of the precipitation over the Amazon (Brubaker et al., 1993; Eltahir and Bras, 1994); these figures increase up to ca. 40 %–60 % based on moisture tagging studies (Risi et al., 2013; Yoshimura et al., 2004).

The LGM is characterised by lower Northern Hemisphere summer insolation, globally cooler temperatures, expanded global ice volumes and lower GHG concentrations than either the MH or the LIG. The MH and LIG (Marine Isotope Stage 5e) periods represent peaks in the present and last interglacial periods, whereas the LGM represents maximum ice extent during the Last Glacial Period. Hence, comparison of these time periods provides a snapshot view of glacial–interglacial variability. The δ18Ospel anomalies are more positive during the LGM than during the MH or LIG, suggesting drier conditions in the ISM, EAM and IAM, supported by simulated changes in δ18Oprecip and precipitation (Fig. 3). Cooler SSTs of approximately 2 C (relative to the MH and LIG) in the ISM and EAM and of approximately 3 C in IAM source areas, along with a ca. 5 % decrease in relative humidity (Yue et al., 2011), would result in a water vapour δ18O signal at the source that is ca. 1 ‰ more depleted than seawater. This depletion results from the temperature dependence of equilibrium fractionation during evaporation and kinetic isotope effects related to humidity (Clark and Fritz, 1997). This fractionation counteracts any impact from enriched seawater δ18O values during the LGM (ca. +1 ‰ relative to the MH or LIG; Waelbroeck et al., 2002). Cooler air temperatures will also result in a depletion of δ18Ospel during the LGM of ca. 0.4 ‰ and 0.6 ‰ for the ISM/EAM and IAM respectively, as a result of water–calcite (or water–aragonite) fractionation (Grossman and Ku, 1986; Tremaine et al., 2011). This has the effect of slightly reducing the regional LGM δ18Ospel signals, although the change is small and within the uncertainty of the regional signals. Enriched δ18Oprecip and δ18Ospel values during the LGM must therefore be caused by a significant decrease in atmospheric moisture and precipitation that resulted from the cooler conditions.

We have used version 2 of the SISAL database (Atsawawaranunt et al., 2018; Comas-Bru et al., 2020a) in our analyses. Despite the fact that SISALv2 includes more than 70 % of known speleothem isotope records, there are still too few records from some regions (e.g. Africa, the Caribbean) to make meaningful analyses. The records for older time periods are also sparse. For example, there are only 14 records from monsoon regions covering the LIG in SISALv2. Nevertheless, our analyses show that there are robust and explicable patterns for most monsoon regions during the Holocene and sufficient records to make meaningful analyses of the LGM and LIG. Whilst there is a need for the generation of new speleothem records from key regions such as northern Africa, further expansion of the SISAL database will certainly provide additional opportunities to analyse the evolution of the monsoons through time.

The impact of age uncertainties, included in SISALv2, are not taken into account in our analyses. Age uncertainties during the Holocene are smaller than the interval used for binning records and the width of the time windows used and, thus, should not have a significant effect on our conclusions. The mean age uncertainty at the LGM and LIG is ca. 430 and 1140 years respectively. However, varying the window length for the selection of LGM and LIG samples from ±500 to ±2000 years, thereby encompassing this uncertainty, has a negligible effect (<0.5 ‰) on the average δ18O values. Thus, the interglacial–glacial contrast in regional δ18Ospel is also robust to age uncertainties.

Isotope-enabled climate models are used in this study to explore observed regional-scale trends in δ18Ospel. There is a limited number of isotope-enabled models, and there are no simulations of the same time period using the same experimental protocol. Although there are simulations of the MH from both ECHAM5-wiso and GISS, for example, these models have different grid resolutions and used different boundary conditions. This could help to explain why the two models yield different estimates of the change in regional δ18Oprecip (of 0.5 ‰) at the MH. However, both models show trends in δ18Oprecip that reproduce the observed changes in regional δ18Ospel (Figs. 3, 4), and this provides a basis for using these models to explore the causes of these trends on different timescales. The failure to reproduce the LGM δ18Ospel signal in SW-SAM in the ECHAM5-wiso model, which precluded a consideration of interglacial–glacial shifts in this region, is a common feature of other isotope-enabled simulations (Caley et al., 2014; Risi et al., 2010). Identifying the underlying relationships between δ18Oprecip and monsoon climate variables using multiple linear regression allows us to identify plausible mechanistic controls on δ18O variability in the monsoon regions. Correlations between δ18O and specific climate variables do not explicitly indicate causality. However, the relationships identified in the MLR model are consistent with the theoretical understanding of oxygen isotope systematics, and the findings of this paper are consistent with existing studies, suggesting that these relationships provide a plausible explanation for observed changes.

This study illustrates a novel data–model approach to investigate the relationship between δ18Ospel and monsoon climate under past conditions: we compare composite regional records and then use multiple linear regression of isotope-enabled palaeoclimate simulations to determine the change in individual climate variables associated with these trends. This obviates the need to use modern δ18Oprecip–climate relationships to explain changes under conditions considerably different from today or to rely on coherency between different palaeohydrological archives which may respond to different climate variables. This model interrogation approach could be employed to address questions about the regional drivers of speleothem records outside the monsoon regions.

5 Conclusions

Geographically distributed speleothem δ18O records and isotope-enabled climate models can be used together to understand the underlying relationships between δ18Ospel and monsoon climate in the past and, therefore, elucidate possible drivers of δ18O variability. Speleothem records, objectively grouped into monsoon regions by record correlation and multivariate ordination techniques, show regional trends that are consistent with changes in summer insolation but modulated by land surface and ocean feedbacks. LGM δ18Ospel signals are best explained by a large decrease in precipitation, as a consequence of lower atmospheric moisture content driven by global cooling. The evolution of δ18Ospel through the Holocene across the global monsoon domain is closely correlated with changes in precipitation, atmospheric circulation and temperature. At the regional scale, our analyses support the increasing number of studies suggesting that East Asian monsoon speleothem δ18O evolution through the Holocene relates to changes in atmospheric circulation (i.e. changes in moisture pathway and/or source). Changes in regional precipitation are the predominant driver of Holocene δ18Ospel evolution in the Indian, southwestern South American and Indonesian–Australian monsoons, although changes in atmospheric circulation also contribute in the Indian and Indonesian–Australian monsoon regions and changes in precipitation recycling appear to be important in southwestern South America.

Code and data availability

The SISAL (Speleothem Isotopes Synthesis and AnaLysis) database version 2 is available through the University of Reading Research Data Archive at (Comas-Bru et al., 2020a). The ECHAM5-wiso MH and LGM simulations are available at (Werner, 2019). The ECHAM LIG simulation is available at (Gierz et al., 2017a). The OIPC mean annual δ18Oprecip data are available at (last access: 10 May 2021, Waterisotopes Database, 2017). CRU-TS4.01 mean annual temperature data are available at (Harris and Jones, 2017). The GISS simulations and code used to generate the figures in this paper are available at (Parker, 2020).


The supplement related to this article is available online at:

Author contributions

The study was designed by SEP, SPH, LCB and NK. MW and ANLG provided climate model outputs. SEP ran the analyses. The first draft of the paper was written by SEP, SPH and LCB, and all authors contributed to the final article.

Competing interests

The authors declare that they have no conflict of interest.


Ideas in this paper were developed at a meeting of the SISAL (Speleothem Isotopes Synthesis and AnaLysis) working group of the Past Global Changes (PAGES) programme. We thank PAGES for their support of this meeting and our colleagues in SISAL for the useful discussions. We thank Gabriele Messori for help with the code for calculating precipitation recycling. We also thank Steven Clemens and two anonymous reviewers for their helpful comments on this paper.

Sarah E. Parker, Sandy P. Harrison and Laia Comas-Bru acknowledge funding support from the ERC-funded project GC2.0 (Global Change 2.0: Unlocking the past for a clearer future). Sandy P. Harrison also acknowledges support from the JPI-Belmont Forum PaCMEDy (Palaeoclimate Constraints on Monsoon Evolution and Dynamics) project funded through NERC.

Financial support

This research has been supported by the European Research Council (GC2.0; grant no. 694481 to Sarah E. Parker, Sandy P. Harrison and Laia Comas-Bru) and the Natural Environment Research Council (grant no. NE/P006752/1 to Sandy P. Harrison).

Review statement

This paper was edited by Zhengtang Guo and reviewed by Steven Clemens and two anonymous referees.


Aggarwal, P. K., Alduchov, O. A., Froehlich, K. O., Araguas-Araguas, L. J., Sturchio, N. C., and Kurita, N.: Stable isotopes in global precipitation: A unified interpretation based on atmospheric moisture residence time, Geophys. Res. Lett., 39, L11705,, 2012. 

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. 

Ayliffe, L. K., Gagan, M. K., Zhao, J. X., Drysdale, R. N., Hellstrom, J. C., Hantoro, W. S., Griffiths, M. L., Scott-Gagan, H., St Pierre, E., Cowley, J. A., and Suwargadi, B. W.: Rapid interhemispheric climate links via the Australasian monsoon during the last deglaciation, Nat. Commun., 4, 1–6,, 2013. 

Bailey, A., Posmentier, E., and Feng, X.: Patterns of evaporation and precipitation drive global isotopic changes in atmospheric moisture, Geophys. Res. Lett., 45, 7093–7101,, 2018. 

Bar-Matthews, M., Ayalon, A., and Kaufman, A.: Late Quaternary paleoclimate in the eastern Mediterranean region from stable isotope analysis of speleothems at Soreq Cave, Israel, Quaternary Res., 47, 155–168,, 1997. 

Battacharya, S. K., Froehlich, K., Aggarwal, P. K., and Kulkarni, K. M.: Isotopic variation in Indian Monsoon precipitation: Records from Bombay and New Delhi, Geophys. Res. Lett., 30, 2285,, 2003. 

Battisti, D. S., Ding, Q., and Roe, G. H.: Coherent pan-Asian climatic and isotopic response to orbital forcing of tropical insolation, J. Geophys. Res.-Atmos., 119, 11997–12020,, 2014. 

Beaufort, L., van der Kaars, S., Bassinot, F. C., and Moron, V.: Past dynamics of the Australian monsoon: precession, phase and links to the global monsoon concept, Clim. Past, 6, 695–706,, 2010. 

Bennett, K. D.: Determination of the number of zones in a biostratigraphical sequence, New Phytol., 132, 155–170,, 1996. 

Berbery, E. H. and Barros, V. R.: The hydrologic cycle of the La Plata basin in South America, J. Hydrometeorol., 3, 630–645,<0630:THCOTL>2.0.CO;2, 2002. 

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

Berger, A. L.: Long-term variations of daily insolation and Quaternary climatic changes, J. Atmos. Sci., 35, 2362–2367,<2362:ltvodi>2.0.CO;2, 1978. 

Berkelhammer, M., Sinha, A., Mudelsee, M., Cheng, H., Edwards, R. L., and Cannariato, K.: Persistent multidecadal power of the Indian Summer Monsoon, Earth Planet. Sc. Lett., 290, 166–172,, 2010. 

Biasutti, M., Voigt, A., Boos, W., Braconnot, P., Hargreaves, J. C., Harrison, S. P., Kang, S. M., Mapes, B. E., Scheff, J., Schumaker, C., Sobel, A. H., and Xie, S. P.: Global energetics and local physics as drivers of past, present and future monsoons, Nat. Geosci., 11, 392–400,, 2018. 

Bin, C., Xiang-De, X., and Tianliang, Z.: Main moisture sources affecting lower Yangtze River Basin in boreal summers during 2004–2009, Int. J. Climatol., 33, 1035–1046,, 2013. 

Boers, N., Rheinwalt, A., Bookhagen, B., Barbosa, H. M. J., Marwan, N., Marengo, J., and Kurths, J.: The South American rainfall dipole: A complex network analysis of extreme events, Geophys. Res. Lett., 41, 7397–7405,, 2014. 

Bowen, G.: Gridded maps of the isotopic composition of meteoric waters,, The University of Utah, available at:, last access: 10 May 2021. 

Bowen, G. J. and Revenaugh, J.: Interpolating the isotopic composition of modern meteoric precipitation, Water Resour. Res., 39, 1299,, 2003. 

Breitenbach, S. F. M., Adkins, J. F., Meyer, H., Marwan, N., Kumar, K. K., and Haug, G. H.: Strong influence of water vapor source dynamics on stable isotopes in precipitation observed in Southern Meghalaya, NE India, Earth Planet. Sc. Lett., 292, 212–220,, 2010. 

Brook, E. J., Harder, S., Severinghaus, J., Steig, E. J., and Sucher, C. M.: On the origin and timing of rapid changes in atmospheric methane during the Last Glacial Period, Global Biogeochem. Cy., 14, 559–572,, 2000. 

Brubaker, K. L., Entekhabi, D., and Eagleson, P. S.: Estimation of continental precipitation recycling, J. Climate, 6, 1077–1089,<1077:EOCPR>2.0.CO;2, 1993. 

Cai, Y., Cheng, H., An, Z., Lawrence Edwards, R., Wang, X., Tan, L., and Wang, J.: Large variations of oxygen isotopes in precipitation over south-central Tibet during Marine Isotope Stage 5, Geology, 38, 243–246,, 2010. 

Cai, Y., Zhang, H., Cheng, H., An, Z., Lawrence Edwards, R., Wang, X., Tan, L., Liang, F., Wang, J., and Kelly, M.: The Holocene Indian monsoon variability over the southern Tibetan Plateau and its teleconnections, Earth Planet. Sc. Lett., 335, 135–144,, 2012. 

Cai, Y., Fung, I. Y., Edwards, R. L., An, Z., Cheng, H., Lee, J. E., Tan, L., Shen, C. C., Wang, X., Day, J. A., Zhou, W., Kelly, M. J., and Chiang, J. C. H.: Variability of stalagmite-inferred Indian monsoon precipitation over the past 252,000 y, P. Natl. Acad. Sci. USA, 112, 2954–2959,, 2015. 

Caley, T., Roche, D. M., Waelbroeck, C., and Michel, E.: Oxygen stable isotopes during the Last Glacial Maximum climate: perspectives from data–model (iLOVECLIM) comparison, Clim. Past, 10, 1939–1955,, 2014. 

Carlson, A. E., LeGrande, A. N., Oppo, D. W., Came, R. E., Schmidt, G. A., Anslow, F. S., Licciardi, J. M., and Obbink, E. A.: Rapid early Holocene deglaciation of the Laurentide ice sheet, Nat. Geosci., 1, 620–624,, 2008. 

Carolin, S. A., Cobb, K. M., Lynch-Stieglitz, J., Moerman, J. W., Partin, J. W., Lejau, S., Malang, J., Clark, B., Tuen, A. A., and Adkins, J. F.: Northern Borneo stalagmite records reveal West Pacific hydroclimate across MIS 5 and 6, Earth Planet. Sc. Lett., 439, 182–193,, 2016. 

Cheng, H., Edwards, R. L., Wang, Y., Kong, X., Ming, Y., Kelly, M. J., Wang, X., Gallup, C. D., and Liu, W.: A penultimate glacial monsoon record from Hulu Cave and two-phase glacial terminations, Geology, 34, 217–220,, 2006. 

Cheng, H., Edwards, R. L., Broecker, W. S., Denton, G. H, Kong, X., Wang, Y., Zhang, R., and Wang, X.: Ice Age Terminations, Science, 326, 248–252,, 2009. 

Cheng, H., Sinha, A., Cruz, F. W., Wang, X., Edwards, R. L., d'Horta, F. M., Ribas, C. C., Vuille, M., Stott, L. D., and Auler, A. S.: Climate change patterns in Amazonia and biodiversity, Nat. Commun., 4, 1411,, 2013. 

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

Cheng, H., Zhang, H., Zhao, J., Li, H., Ning, Y., and Kathayat, G.: Chinese stalagmite paleoclimate researches: A review and perspective, Sci. China Earth Sci., 62, 1489–1513,, 2019. 

Clark, I. and Fritz, P.: Environmental Isotopes in Hydrology, Lewis Publishers, New York, ISBN: 1-56670-249-6, 1997. 

Cleveland, W. S. and Devlin, S. J.: Locally weighted regression: An approach to regression analysis by local fitting, J. Am. Stat. Assoc., 83, 596–610,, 1988. 

Cobb, K. M., Adkins, J. F., Partin, J. W., and Clark, B.: Regional-scale climate influences on temporal variations of rainwater and cave dripwater oxygen isotopes in northern Borneo, Earth Planet. Sc. Lett., 263, 207–220,, 2007. 

Cole, J. E., Rind, D., Webb, R. S., Jouzel, J., and Healy, R.: Climatic controls on interannual variability of precipitation δ18O: Simulated influence of temperature, precipitation amount, and vapor source region, J. Geophys. Res., 104, 14223–14235,, 1999. 

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. 

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. 

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. 

Coplen, T. B., Kendall, C., and Hopple, J.: Comparison of stable isotope reference samples, Nature, 302, 236–238,, 1983. 

Craig, H. and Gordon, L. I.: Deuterium and oxygen 18 variations in the ocean and the marine atmosphere, in: Stable isotopes in Oceanographic Studies and Paleotemperatures, edited by: Tongiorgi, E., Consiglio Nazionake delle Ricerche Laboratoria di Geologia Nucleare, Pisa, 9–130, 1965. 

Cruz, F. W., Burns, S. J., Karmann, I., Sharp, W. D., Vuille, M., Cardoso, A. O., Ferrari, J. A., Silva Dias, P. L., and Viana, O.: Insolation-driven changes in atmospheric circulation over the past 116,000 years in subtropical Brazil, Nature, 434, 63–66,, 2005. 

Cruz, F. W., Burns, S. J., Karmann, I., Sharp, W. D., and Vuille, M.: Reconstruction of regional atmospheric circulation features during the late Pleistocene in subtropical Brazil from oxygen isotope composition of speleothems, Earth Planet. Sc. Lett., 248, 495–507,, 2006. 

Cruz, F. W., Vuille, M., Burns, S. J., Wang, X., Cheng, H., Werner, M., Lawrence Edwards, R., Karmann, I., Auler, A. S., and Nguyen, H.: Orbitally driven east-west antiphasing of South American precipitation, Nat. Geosci., 2, 210–214,, 2009. 

D'Abreton, P. C. and Tyson, P. D.: Three-dimensional kinematic trajectory modelling of water vapour transport over Southern Africa, Water SA, 22, 297–306, 1996. 

Dallmeyer, A., Claussen, M., and Otto, J.: Contribution of oceanic and vegetation feedbacks to Holocene climate change in monsoonal Asia, Clim. Past, 6, 195–218,, 2010. 

Dansgaard, W.: Stable isotopes in precipitation, Tellus, 16, 436–468,, 1964. 

Deininger, M., Ward, B. M., Novello, V. F., and Cruz, F. W.: Late Quaternary variations in the South American monsoon system as inferred by speleothems–New perspectives using the SISAL database, Quaternary, 2, 6,, 2019. 

Ding, Y. and Chan, J. C. L.: The East Asian summer monsoon: an overview, Meteorol. Atmos. Phys., 89, 117–142,, 2005. 

Dong, J., Wang, Y., Cheng, H., Hardt, B., Lawrence Edwards, R., Kong, X., Wu, J., Chen, S., Liu, D., Jiang, X., and Zhao, K.: A high-resolution stalagmite record of the Holocene East Asian monsoon from Mt Shennongjia, central China, Holocene, 20, 257–264,, 2010. 

Drumond, A., Nieto, R., Gimeno, L., and Ambrizzi, T.: A Lagrangian identification of major sources of moisture over Central Brazil and la Plata Basin, J. Geophys. Res., 113, D14128,, 2008. 

Drumond, A., Nieto, R., Trigo, R., Ambrizzi, T., Souza, E., and Gimeno, L.: A Lagrangian identification of the main sources of moisture affecting northeastern Brazil during its pre-rainy and rainy seasons, PLoS One, 5, e11205,, 2010. 

Durán-Quesada, A. M., Gimeno, L., Amador, J. A., and Nieto, R.: Moisture sources for Central America: Identification of moisture sources using a Lagrangian analysis technique, J. Geophys. Res., 115, D05103,, 2010. 

Dutton, A. and Lambeck, K.: Ice volume and sea level during the last interglacial, Science, 337, 216–219, 2012. 

Dykoski, C. A., Edwards, R. L., Cheng, H., Yuan, D., Cai, Y., Zhang, M., Lin, Y., Qing, J., An, Z., and Revenaugh, J.: A high-resolution, absolute-dated Holocene and deglacial Asian monsoon record from Dongge Cave, China, Earth Planet. Sc. Lett., 233, 71–86,, 2005. 

Eltahir, E. A. B. and Bras, R. L.: Precipitation recycling in the Amazon basin, Q. J. Roy. Meteor. Soc., 120, 861–880,, 1994. 

Fairchild, I. J. and Baker, A.: Speleothem Science: From Process to Past Environments, 1st edn., Wiley-Blackwell, Chichester, 2012. 

Fleitmann, D., Burns, S. J., Mudelsee, M., Neff, U., Kramers, J., Mangini, A., and Matter, A.: Holocene Forcing of the Indian Monsoon Recorded in a Stalagmite from Southern Oman, Science, 300, 1737–1739,, 2003. 

Fleitmann, D., Burns, S. J., Neff, U., Mudelsee, M., Mangini, A., and Matter, A.: Palaeoclimatic interpretation of high-resolution oxygen isotope profiles derived from annually laminated speleothems from Southern Oman, Quaternary Sci. Rev., 23, 935–945,, 2004. 

Friedman, I., Harris, J. M., Smith, G. I., and Johnson, C. A.: Stable isotope composition of waters in the Great Basin, United States 1. Air-mass trajectories, J. Geophys. Res., 107, 4400,, 2002. 

Gat, J. R.: Oxygen and hydrogen isotopes in the hydrological cycle, Annu. Rev. Earth Pl. Sc., 24, 225–262,, 1996. 

Gierz, P., Werner, M., and Lohmann, G.: Simulating climate and stable water isotopes during the last interglacial using a coupled climate-isotope model, link to model results, PANGAEA [data set],, 2017a. 

Gierz, P., Werner, M., and Lohmann, G.: Simulating climate and stable water isotopes during the Last Interglacial using a coupled climate-isotope model, J. Adv. Model. Earth Sy., 9, 2027–2045,, 2017b. 

Goldschneider, N., Chen, Z., Auler, A. S., Bakalowicz, M., Broda, S., Drew, D., Hartmann, J., Jiang, G., Moosdorf, N., Stevanovic, Z., and Veni, G.: Global distribution of carbonate rocks and karst water resources, Hydrogeol. J., 28, 1661–1677,, 2020. 

Gower, J. C.: Some distance properties of latent root and vector methods used in multivariate analysis, Biometrika, 53, 325–338,, 1966. 

Griffiths, M. L., Drysdale, R. N., Gagan, M. K., Zhao, J. X., Ayliffe, L. K., Hellstrom, J. C., Hantoro, W. S., Frisia, S., Feng, Y. X., Cartwright, I., Pierre, E. S., Fischer, M. J., and Suwargadi, B. W.: Increasing Australian-Indonesian monsoon rainfall linked to early Holocene sea-level rise, Nat. Geosci., 2, 636–639,, 2009. 

Grossman, E. L. and Ku, T. L.: Oxygen and carbon isotope fractionation in biogenic aragonite: Temperature effects, Chem. Geol. Isot. Geosci. Sect., 59, 59–74,, 1986. 

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,, 2012. 

Harris, I. C. and Jones, P. D.: CRU TS4.01: Climatic Research Unit (CRU) Time-Series (TS) version 4.01 of high-resolution gridded data of month-by-month variation in climate (Jan. 1901-Dec. 2016), Centre for Environmental Data Analysis [data set],, 2017. 

Harris, I., Osborn, T. J., Jones, P. and Lister, D.: Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset, Sci. Data, 7, 109,, 2020. 

Hong, Y. T., Hong, B., Lin, Q. H., Zhu, Y. X., Shibata, Y., Hirota, M., Uchida, M., Leng, X. T., Jiang, H. B., Xu, H., Wang, H., and Yi, L.: Correlation between Indian Ocean summer monsoon and North Atlantic climate during the Holocene, Earth Planet. Sc. Lett., 211, 371–380,, 2003. 

Hu, J., Emile-Geay, J., Tabor, C., Nusbaumer, J., and Partin, J.: Deciphering oxygen isotope records from Chinese speleothems with an isotope-enabed climate model, Paleoceanography and Paleoclimatology, 43, 2098–2112,, 2019. 

IAEA/WMO: Global Network of Isotopes in Precipitation, The GNIP Database, available at: (last access: 20 May 2021), 2018. 

Indermühle, A., Stocker, T. F., Joos, F., Fischer, H., Smith, H. J., Wahlen, M., Deck, B., Mastroianni, D., Tschumi, J., Blunier, T., Meyer, R., and Stauffer, B.: Holocene carbon-cycle dynamics based on CO2 trapped in ice at Taylor Dome, Antarctica, Nature, 398, 121–126,, 1999. 

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. 

Kärkkäinen, T. and Saarela, M.: Robust Principal Component Analysis of Data with Missing Values, in: Machine Learning and Data Mining in Pattern Recognition. MLDM 2015, edited by: Perner, P., Springer, Cham, Lecture Notes in Computer Science, vol. 9166,, 2015. 

Kathayat, G., Cheng, H., Sinha, A., Spötl, C., Edwards, R. L., Zhang, H., Li, X., Yi, L., Ning, Y., Cai, Y., Lui, W. L., and Breitenbach, S. F. M.: Indian monsoon variability on millennial-orbital timescales, Sci. Rep., 6, 24374,, 2016. 

Kaushal, N., Breitenbach, S. F. M., Lechleitner, F. A., Sinha, A., Tewari, V. C., Ahmad, S. M., Berkelhammer, M., Band, S., Yadava, M., Ramesh, R., and Henderson, G. M.: The Indian Summer Monsoon from a Speleothem δ18O Perspective–A Review, Quaternary, 1, 29,, 2018. 

Kennett, D. J., Breitenbach, S. F. M., Aquino, V. V., Asmerom, Y., Awe, J., Baldini, J. U. L., 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. 

Krause, C. E., Gagan, M. K., Dunbar, G. B., Hantoro, W. S., Hellstrom, J. C., Cheng, H., Edwards, R. L., Suwargadi, B. W., Abram, N. J., and Rifai, H.: Spatio-temporal evolution of Australasian monsoon hydroclimate over the last 40,000 years, Earth Planet. Sc. Lett., 513, 103–112,, 2019. 

Kuhnt, W., Holbourn, A., Xu, J., Opdyke, B., De Deckker, P., Röhl, U., and Mudelsee, M.: Southern Hemisphere control on Australian monsoon variability during the late deglaciation and Holocene, Nat. Commun., 6, 1–7,, 2015. 

Kutzbach, J., Bonan, G., Foley, J., and Harrison, S. P.: Vegetation and soil feedbacks on the response of the African monsoon to orbital forcing in the early to middle Holocene, Nature, 384, 623–626,, 1996. 

Kutzbach, J. E., Liu, X., Liu, Z., and Chen, G.: Simulation of the evolutionary response of global summer monsoons to orbital forcing over the past 280,000 years, Clim. Dynam., 30, 567–579,, 2008. 

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

Lachniet, M. S.: Are aragonite stalagmites reliable paleoclimate proxies? Tests for oxygen isotope time-series replication and equilibrium, Bull. Geol. Soc. Am., 127, 1521–1533,, 2015. 

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. 

Legendre, P. and Legendre, L.: Numerical Ecology, 2nd edn., Elsevier B.V., Amsterdam, 1998. 

LeGrande, A. N. and Schmidt, G. A.: Global gridded data set of the oxygen isotopic composition in seawater, Geophys. Res. Lett., 33, L12604,, 2006. 

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

Lewis, S. C., LeGrande, A. N., Kelley, M., and Schmidt, G. A.: Water vapour source impacts on oxygen isotope variability in tropical precipitation during Heinrich events, Clim. Past, 6, 325–343,, 2010. 

Lewis, S. C., LeGrande, A. N., Schmidt, G. A., and Kelley, M.: Comparison of forced ENSO-like hydrological expressions in simulations of the preindustrial and mid-Holocene, J. Geophys, Res.-Atmos., 119, 7064–7082,, 2014. 

Li, J., Dodson, J., Yan, H., Wang, W., Innes, J. B., Zong, Y., Zhang, X., Xu, Q., Ni, J., and Lu, F.: Quantitative Holocene climatic reconstructions for the lower Yangtze region of China, Clim. Dynam., 50, 1101–1113,, 2018. 

Licciardi, J. M., Clark, P. U., Jenson, J. W., and Macayeal, D. R.: Deglaciation of a soft-bedded Laurentide Ice Sheet, Quaternary Sci. Rev., 17, 427–448,, 1998. 

Liu, Z., Wen, X., Brady, E. C., Otto-Bliesner, B., Yu, G., Lu, H., Cheng, H., Wang, Y., Zheng, W., Ding, Y., Edwards, R. L., Cheng, J., Liu, W., and Yang, H.: Chinese cave records and the East Asia Summer Monsoon, Quaternary Sci. Rev., 83, 115–128,, 2014. 

Maher, B. A.: Holocene variability of the East Asian summer monsoon from Chinese cave records: A re-assessment, Holocene, 18, 861–866,, 2008. 

Maher, B. A.: Palaeoclimatic records of the loess/palaeosol sequences of the Chinese Loess Plateau, Quaternary Sci. Rev., 154, 23–84,, 2016. 

Maher, B. A. and Thompson, R.: Oxygen isotopes from Chinese caves: Records not of monsoon rainfall but of circulation regime, J. Quat. Sci., 27, 615–624,, 2012. 

Marlon, J. R., Bartlein, P. J., Carcaillet, C., Gavin, D. G., Harrison, S. P., Higuera, P. E., Joos, F., Power, M. J., and Prentice, I. C.: Climate and human influences on global biomass burning over the past two millennia, Nat. Geosci., 1, 697–702,, 2008. 

Marzin, C. and Braconnot, P.: Variations of Indian and African monsoons induced by insolation changes at 6 and 9.5 kyr BP, Clim. Dynam., 33, 215–231,, 2009. 

McDermott, F.: Palaeo-climate reconstruction from stable isotope variations in speleothems: A review, Quaternary Sci. Rev., 23, 901–918,, 2004. 

Misra, P., Tandon, S. K., and Sinha, R.: Holocene climate records from lake sediments in India: Assessment of coherence across climate zones, Earth-Sci. Rev., 190, 370–397,, 2019. 

Moerman, J. W., Cobb, K. M., Adkins, J. F., Sodemann, H., Clark, B., and Tuen, A. A.: Diurnal to interannual rainfall δ18O variations in northern Borneo driven by regional hydrology. Earth Planet. Sc. Lett., 369, 108–119,, 2013. 

Mohtadi, M., Oppo, D. W., Steinke, S., Stuut, J. B. W., De Pol-Holz, R., Hebbeln, D., and Lückge, A.: Glacial to Holocene swings of the Australian-Indonesian monsoon, Nat. Geosci., 4, 540–544,, 2011. 

Nivet, F., Bergonzini, L., Mathé, P. E., Noret, A., Monvoisin, G., Majule, A., and Williamson, D.: Influence of the balance of the intertropical front on seasonal variations of the isotopic composition in rainfall at Kisiba Masoko (Rungwe Volcanic Province, SW, Tanzania), Isotopes Environ. Health Stud., 54, 352–369,, 2018. 

Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., O'Hara, R. B., Simpson, G. L., Solymos, P., Stevens, M. H. H., Szoecs, E., and Wagner, H.: vegan: Community Ecology Package, R package version 2.5-2, Cran R, 2019. 

Otto-Bliesner, B. L., Braconnot, P., Harrison, S. P., Lunt, D. J., Abe-Ouchi, A., Albani, S., Bartlein, P. J., Capron, E., Carlson, A. E., Dutton, A., Fischer, H., Goelzer, H., Govin, A., Haywood, A., Joos, F., LeGrande, A. N., Lipscomb, W. H., Lohmann, G., Mahowald, N., Nehrbass-Ahles, C., Pausata, F. S. R., Peterschmitt, J.-Y., Phipps, S. J., Renssen, H., and Zhang, Q.: The PMIP4 contribution to CMIP6 – Part 2: Two interglacials, scientific objective and experimental design for Holocene and Last Interglacial simulations, Geosci. Model Dev., 10, 3979–4003,, 2017. 

Parker, S. E.: Speleothem_orbital_monsoon_var, Zenodo [code],, 2020. 

Pausata, F. S. R., Battisti, D. S., Nisancioglu, K. H., and Bitz, C. M.: Chinese stalagmite δ18O controlled by changes in the Indian monsoon during a simulated Heinrich event, Nat. Geosci., 4, 474–480,, 2011. 

Quigley, M. C., Horton, T., Hellstrom, J. C., Cupper, M. L., and Sandiford, M.: Holocene climate change in arid Australia from speleothem and alluvial records, Holocene, 20, 1093–1104,, 2010. 

R Core Team: R: A Language and Environment for Statistical Computing, available at: (last access: 31 January 2021), 2019. 

Rachmayani, R., Prange, M., and Schulz, M.: North African vegetation–precipitation feedback in early and mid-Holocene climate simulations with CCSM3-DGVM, Clim. Past, 11, 175–185,, 2015. 

Rao, C. R.: The use and interpretation of principal component analysis in applied research, Sankhya Indian J. Stat. Ser. A, 329–358, 1964. 

Risi, C., Bony, S., and Vimeux, F.: Influence of convective processes on the isotopic composition (δ18O and δD) of precipitation and water vapor in the tropics: 2. Physical interpretation of the amount effect, J. Geophys. Res., 113, D19306,, 2008. 

Risi, C., Bony, S., Vimeux, F., and Jouzel, J.: Water-stable isotopes in the LMDZ4 general circulation model: Model evaluation for present-day and past climates and applications to climatic interpretations of tropical isotopic records, J. Geophys. Res., 115, D12118,, 2010. 

Risi, C., Noone, D., Frankenberg, C., and Worden, J.: Role of continental recycling in intraseasonal variations of continental moisture as deduced from model simulations and water vapor isotopic measurements, Water Resour. Res., 49, 4136–4156,, 2013. 

Rohlf, F. J.: An Empirical Comparison of Three Ordination Techniques in Numerical Taxonomy, Syst. Biol., 21, 271–280,, 1972. 

Rozanski, K., Araguás-Araguás, L., and Gonfiantini, R.: Isotopic patterns in modern global precipitation, Geophys Monogr., 78, 1–36,, 1993. 

Salati, E., Dall'Olio, A., Matsui, E., and Gat, J. R.: Recycling of water in the Amazon Basin: An isotopic study, Water Resour. Res., 15, 1250–1258,, 1979. 

Schäfer-Neth, C. and Paul, A.: The Atlantic Ocean at the Last Glacial Maximum: 1. Objective mapping of the GLAMAP sea-surface conditions, in: The South Atlantic in the Late Quaternary, edited by: Wefer, G., Mulitza, S., and Ratmeyer, V., Springer, Berlin, Heidelberg, 531–548, 2003. 

Sinha, A., Kathayat, G., Cheng, H., Breitenbach, S.F., Berkelhammer, M., Mudelsee, M., Biswas, J., and Edwards, R.L.: Trends and oscillations in the Indian summer monsoon rainfall over the last two millennia, Nat. Commun., 6, 1–8,, 2015. 

Sjolte, J. and Hoffmann, G.: Modelling stable water isotopes in monsoon precipitation during the previous interglacial, Quaternary Sci. Rev., 85, 119–135,, 2014. 

Sowers, T.: Ice core records of atmospheric N2O covering the last 106,000 years, Science, 301, 945–948,, 2003. 

Steinke, S., Mohtadi, M., Prange, M., Varma, V., Pittauerova, D., and Fischer, H. W.: Mid- to Late-Holocene Australian-Indonesian summer monsoon variability, Quaternary Sci. Rev., 93, 142–154,, 2014. 

Tabor, C. R., Otto-Bliesner, B. L., Brady, E. C., Nusbaumer, J., Zhu, J., Erb, M. P., Wong, T. E., Liu, Z., and Noone, D.: Interpreting Precession-Driven δ18O Variability in the South Asian Monsoon Region, J. Geophys. Res.-Atmos., 123, 5927–5946,, 2018. 

Tan, M.: Circulation effect: climatic significance of the short term variability of the oxygen isotopes in stalagmites from monsoonal China, Quat. Sci., 29, 851–862, 2009. 

Tan, M.: Tread-wind driven inverse coupling between stalagmite δ18O from monsoon region of China and large scale temperature-circulation effect on decadal to precessional timescales, Quat. Sci., 31, 1086–1097, 2011. 

Tan, M.: Circulation effect: Response of precipitation δ18O to the ENSO cycle in monsoon regions of China, Clim. Dynam., 42, 1067–1077,, 2014. 

Tan, L. C., Cai, Y., Cheng, H., Edwards, R. L., Shen, C., Gao, Y., and An, Z.: Climate significance of speleothem δ18O from central China on decadal timescale, J. Asian Earth Sci., 106, 150–155,, 2015. 

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. 

Trenberth, K. E., Stepaniak, D. P., and Caron, J. M.: The global monsoon as seen through the divergent atmospheric circulation, J. Climate, 13, 3969–3993,<3969:TGMAST> 2.0.CO;2, 2000. 

Varma, V., Prange, M., Merkel, U., Kleinen, T., Lohmann, G., Pfeiffer, M., Renssen, H., Wagner, A., Wagner, S., and Schulz, M.: Holocene evolution of the Southern Hemisphere westerly winds in transient simulations with global climate models, Clim. Past, 8, 391–402,, 2012. 

Wackerbarth, A., Langebroek, P. M., Werner, M., Lohmann, G., Riechelmann, S., Borsato, A., and Mangini, A.: Simulated oxygen isotopes in cave drip water and speleothem calcite in European caves, Clim. Past, 8, 1781–1799,, 2012. 

Waelbroeck, C., Labeyrie, L., Michel, E., Duplessy, J. C., McManus, J. F., Lambeck, K., Balbon, E., and Labracherie, M.: Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records, Quaternary Sci. Rev., 21, 295–305,, 2002. 

Wang, B. and Ding, Q.: Global monsoon: Dominant mode of annual variation in the tropics, Dynam. Atmos. Oceans, 44, 165–183,, 2008. 

Wang, P. X., Wang, B., Cheng, H., Fasullo, J., Guo, Z. T., Kiefer, T., and Liu, Z. Y.: The global monsoon across timescales: coherent variability of regional monsoons, Clim. Past, 10, 2007–2052,, 2014. 

Wang, X., Auler, A. S., Edwards, R. L., Cheng, H., Ito, E., and Solheid, M.: Interhemispheric anti-phasing of rainfall during the last glacial period, Quaternary Sci. Rev., 25, 3391–3403,, 2006. 

Wang, Y., Cheng, H., Edwards, R. L., Kong, X., Shao, X., Chen, S., Wu, J., Jiang, X., Wang, X., and An, Z.: Millennial- and orbital-scale changes in the East Asian monsoon over the past 224,000 years, Nature, 451, 1090–1093,, 2008. 

Wang, Y. J., Cheng, H., Edwards, R. L., An, Z. S., Wu, J. Y., Shen, C. C., and Dorale, J. A.: A high-resolution absolute-dated late Pleistocene monsoon record from Hulu Cave, China, Science, 294, 2345–2348,, 2001. 

Waterisotopes Database: OIPC mean annual δ18Oprecip data,, University of Utah [data set], available at:, last access: 10 May 2021. 

Werner, M.: ECHAM5-wiso simulation data – present-day, mid-Holocene, and Last Glacial Maximum, PANGAEA [data set],, 2019. 

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. 

Werner, M., Jouzel, J., Masson-Delmotte, V., and Lohmann, G.: Reconciling glacial Antarctic water stable isotopes with ice sheet topography and the isotopic paleothermometer, Nat. Commun., 9, 1–10,, 2018. 

Wickham, H.: ggplot2: Elegant Graphics for Data Analysis, Springer-Verlag, New York, available at: (last access: 31 January 2021), 2016. 

Wurtzel, J. B., Abram, N. J., Lewis, S. C., Bajo, P., Hellstrom, J. C., Troitzsch, U., and Heslop, D.: Tropical Indo-Pacific hydroclimate response to North Atlantic forcing during the last deglaciation as recorded by a speleothem from Sumatra, Indonesia, Earth Planet. Sci. Lett., 492, 264–278,, 2018. 

Wyrwoll, K. H. and Miller, G. H.: Initiation of the Australian summer monsoon 14,000 years ago, Quatern. Int., 83, 119–128,, 2001. 

Xu, X., Werner, M., Butzin, M., and Lohmann, G.: Water isotope variations in the global ocean model MPI-OM, Geosci. Model Dev., 5, 809–818,, 2012. 

Yang, X., Liu, J., Liang, F., Yuan, D., Yang, Y., Lu, Y., and Chen, F.: Holocene stalagmite δ18O records in the East Asian monsoon region and their correlation with those in the Indian monsoon region, Holocene, 24, 1657–1664,, 2014. 

Yonge, C. J., Ford, D. C., Gray, J., and Schwarcz, H. P.: Stable isotope studies of cave seepage water, Chem. Geol. Isot. Geosci. Sect., 58, 97–105,, 1985. 

Yoshimura, K., Oki, T., Ohte, N., and Kanae, S.: Colored moisture analysis estimates of variations in 1998 Asian Monsoon water sources, J. Meteorol. Soc. Jpn., 82, 1315–1329,, 2004. 

Yuan, D., Cheng, H., Edwards, R. L., Dykoski, C. A., Kelly, M. J., Zhang, M., Qing, J., Lin, Y., Wang, Y., Wu, J., Dorale, J. A., An, Z., and Cai, Y.: Timing, duration, and transitions of the Last Interglacial Asian monsoon, Science, 304, 575–578,, 2004. 

Yue, X., Wang, H., Liao, H., and Jiang, D.: Simulation of the Direct Radiative Effect of Mineral Dust Aerosol on the Climate at the Last Glacial Maximum, J. Climate, 24, 843–858,, 2011. 

Zhang, H., Ait Brahim, Y., Li, H., Zhao, J., Kathayat, G., Tian, Y., Baker, J., Wang, J., Zhang, F., Ning, Y., Edwards, R. L., and Cheng, H.: The Asian Summer monsoon: Teleconnections and forcing Mechanisms–A review from Chinese speleothem δ18O records, Quaternary, 2, 26,, 2019. 

Zhang, X., Jin, L., Chen, J., Lu, H., and Chen, F.: Lagged response of summer precipitation to insolation forcing on the northeastern Tibetan Plateau during the Holocene, Clim. Dynam., 50, 3117–3129,, 2018. 

Zhao, Y. and Harrison, S. P.: Mid-Holocene monsoons: A multi-model analysis of the inter-hemispheric differences in the responses to orbital forcing and ocean feedbacks, Clim. Dynam., 39, 1457–1487,, 2012. 

Zhao, Y., Yu, Z., Chen, F., Zhang, J., and Yang, B.: Vegetation response to Holocene climate change in monsoon-influenced region of China, Earth-Sci. Rev., 97, 242–256,, 2009. 

Zhou, W., Yu, X., Jull, A. J. T., Burr, G., Xiao, J. Y., Lu, X., and Xian, F.: High-resolution evidence from southern China of an early Holocene optimum and a mid-Holocene dry event during the past 18,000 years, Quaternary Res., 62, 39–48,, 2004. 

Short summary
Regional trends in the oxygen isotope (δ18O) composition of stalagmites reflect several climate processes. We compare stalagmite δ18O records from monsoon regions and model simulations to identify the causes of δ18O variability over the last 12 000 years, and between glacial and interglacial states. Precipitation changes explain the glacial–interglacial δ18O changes in all monsoon regions; Holocene trends are due to a combination of precipitation, atmospheric circulation and temperature changes.