# Speleothem records of monsoon interannual-interdecadal variability through the Holocene

Sarah E. Parker, Sandy P. Harrison, Pascalé Braconnot

► To cite this version:

Sarah E. Parker, Sandy P. Harrison, Pascalé Braconnot. Speleothem records of monsoon interannual-interdecadal variability through the Holocene. Environmental Research Communications, 2021, 3 (12), 10.1088/2515-7620/AC3EAA. hal-03524583

HAL Id: hal-03524583

https://hal.science/hal-03524583

Submitted on 14 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LETTER • OPEN ACCESS

Speleothem records of monsoon interannual-interdecadal variability through the Holocene

To cite this article: Sarah E Parker et al 2021 *Environ. Res. Commun.* 3 121002

View the [article online](#) for updates and enhancements.

You may also like

- *Generation of obliquely incident ions using phase-shifted RF voltages applied on rod electrodes*  
  Akio Ui, Yosuke Sato, Toshiyuki Sasaki et al.

- *Ion energy and angular distributions in planar Ar/O₂ inductively coupled plasmas: hybrid simulation and experimental validation*  
  Yu-Ru Zhang, Zhuan-Zhuan Zhao, Chan Xue et al.

- *Bit selection using field drilling data and mathematical investigation*  
  M S Momeni, S Ridha, S J Hosseini et al.
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

13 August 2021

REVISED

24 November 2021

ACCEPTED FOR PUBLICATION

30 November 2021

PUBLISHED

13 December 2021

Original content from this work may be used under the terms of the *Creative Commons Attribution 4.0 licence*.

Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<figure><img src="image_3.png" /></figure>

LETTER

# Speleothem records of monsoon interannual-interdecadal variability through the Holocene

Sarah E Parker¹ Ⓓ, Sandy P Harrison¹ Ⓓ and Pascale Braconnot²

¹ Department of Geography and Environmental Science, University of Reading, Reading, RG6 6AB, United Kingdom

² Laboratoire des Sciences du Climat et de l'Environnement-IPSL, Unité Mixte CEA-CNRS-UVSQ, Université Paris-Saclay, Orme des Merisiers, Gif-sur-Yvette, France

E-mail: s.parker@pgr.reading.ac.uk

**Keywords:** monsoon, palaeoclimate, speleothems, climate variability, interannual variability, Holocene

Supplementary material for this article is available [online](#)

## Abstract

Modern observations show considerable interannual to interdecadal variability in monsoon precipitation. However, there are few reconstructions of variability at this timescale through the Holocene, and there is therefore less understanding of how changes in external forcing might have affected monsoon variability in the past. Here, we reconstruct the evolution of the amplitude of interannual to interdecadal variability (IADV) in the East Asian, Indian and South American monsoon regions through the Holocene using a global network of high-resolution speleothem oxygen isotope ($\delta^{18}$O) records. We reconstruct changes in IADV for individual speleothem records using the standard deviation of $\delta^{18}$O values in sliding time windows after correcting for the influence of confounding factors such as variable sampling resolution, growth rates and mean climate. We then create composites of IADV changes for each monsoon region. We show that there is an overall increase in $\delta^{18}$O IADV in the Indian monsoon region through the Holocene, with an abrupt change to present-day variability at ~2 ka. In the East Asian monsoon, there is an overall decrease in $\delta^{18}$O IADV through the Holocene, with an abrupt shift also seen at ~2 ka. The South American monsoon is characterised by large multi-centennial shifts in $\delta^{18}$O IADV through the early and mid-Holocene, although there is no overall change in variability across the Holocene. Our regional IADV reconstructions are broadly reproduced by transient climate-model simulations of the last 6 000 years. These analyses indicate that there is no straightforward link between IADV and changes in mean precipitation, or between IADV and orbital forcing, at a regional scale.

## 1. Introduction

More than two thirds of the global population live in regions which are dependent on monsoon rainfall; interannual variability in precipitation has a significant impact on the livelihoods of these people (Wang et al. 2021). The variability of precipitation in the tropics on interannual to interdecadal timescales (interannual to interdecadal variability: IADV) is influenced by sea surface temperature (SST) changes: the El Niño Southern Oscillation (ENSO) has a major influence on monsoon precipitation IADV (McPhaden et al. 2006, Messié and Chavez 2011, Wang et al. 2012). SST variability in the extratropical Pacific (Pacific Decadal Oscillation, PDO, Krishnan and Sugi 2003, Yoon and Yeh 2010), Indian Ocean (Indian Ocean Dipole, IOD, Ashok et al. 2001, Saji and Yamagata 2003, Ummenhofer et al. 2011;) and Atlantic Ocean (Atlantic Multidecadal Oscillation, AMO, Grimm and Zilli 2009) also influence the IADV of monsoon precipitation. These modes interact with one another to drive complex and considerable interannual fluctuations in precipitation. They are often associated with extreme events over monsoon regions (Kane 1999, Kirono et al., 1999) with major impacts on water resources, crop production (Phillips et al. 2012, Iizumi et al. 2014, Ray et al. 2015) and long-term impacts on health (Bouma and Kaay, 1996, Gagnon et al. 2001, Hashizume et al. 2012).
<figure><img src="image_4.png" /></figure>

Although there has been considerable focus on changes in monsoon IADV in the recent period (Wang et al 2012, Yim et al 2014), observational records are too short to sample IADV sufficiently and it is therefore difficult to examine how this variability responds to external forcing. Future changes in monsoon precipitation IADV in response to increasing CO2 concentrations are difficult to predict because of the large uncertainties related to internal variability (Wang et al 2021). Climate models disagree, for example, on the how the variability of both ENSO (Stevenson 2012, Chen et al 2017, Brown et al 2020) and the IOD (Hui and Zheng 2018, Zheng et al 2013, Ng et al 2018) will change in the future. The palaeorecord provides an alternative way of examining how changes in external forcing have influenced IADV in monsoon regions. Specifically, the Holocene (11.7 thousand years ago, ka, to present) provides an opportunity to examine how monsoon IADV responds to changes in the seasonal and latitudinal distribution of insolation resulting from changes in the Earth's orbit. In the early Holocene, northern hemisphere summer insolation was greater than today whilst southern hemisphere summer insolation was reduced compared to present. Through the Holocene, summer insolation gradually declined (increased) to present-day levels in the northern (southern) hemisphere. These insolation changes altered the mean climate (Wang et al 2014, Zhang et al 2021) and likely had a marked impact on IADV.

There are only a few reconstructions of Holocene changes in precipitation IADV from tropical continents. Lake records of hydrological variability from the circum-Pacific region show changing precipitation variability through the Holocene, interpreted as ENSO variability. Records show lower precipitation variability in the mid-Holocene than today, with higher variability during the late Holocene (Riedinger et al 2002, Conroy et al 2008, Zhang et al 2014, Thompson et al 2017). The time when IADV increases from the MH minimum varies between records. For example, in the Galapagos Islands, Bainbridge Crater Lake shows increasing ENSO variance from 2 ka onwards (Riedinger et al 2002) and El Junco Lake shows an increase at ~4 ka (Conroy et al 2008, Zhang et al 2014). The red-light reflectivity record (Moy et al 2002) from lake Pallcacocha, Ecuador, indicates the increase in variability begins at ~7 ka, whilst a grey-scale record from the same site (Rodbell 1999) shows this change beginning at ~5 ka. Reconstructions of early Holocene IADV also vary: an algal lipid concentration and hydrogen isotope record from El Junco Lake shows fluctuations between high and low variability (Zhang et al 2014), whilst sediment records from El Junco and Pallcacocha lakes show low variability through the early Holocene (Rodbell 1999, Moy et al 2002, Conroy et al 2008).

The differences between existing reconstructions highlights the limitations of lake sediment records. For example, the interpretation of the Pallcacocha records has been questioned because the highly non-linear response of sediment accumulation to precipitation could bias the reconstructions (Emile-Geay and Tingley 2016). The lack of a clear wet-dry pattern to ENSO events in the area (Schneider et al 2018, Kiefer and Karamperidou 2019) and the influence of catchment differences on lake sediment records (Schneider et al 2018) also reduces the reliability of these records as indicators of IADV and ENSO. Furthermore, the spatial complexity of IADV necessitates regional-scale reconstructions to understand differences in IADV changes.

The Gallapagos and Pallcacocha records were a focus for reconstructions of IADV because, compared to most other terrestrial records, they have quasi-annual resolution. Speleothems provide an alternative source of high-resolution data in monsoon regions and can therefore be used to investigate precipitation IADV over land. Speleothems record changes in the strength of the monsoon, via a combination of changes in regional precipitation and circulation changes that cause changes in moisture source and transport pathway (Sinha et al 2015, Cheng et al 2019, Parker et al 2021). However, there are several factors that influence the oxygen isotopic composition of speleothems and might complicate the interpretation of short-term climate variability, including changing growth rates and the mixing of meteoric water in the epikarst (Lachniet 2009) although compositing of speleothem records can be used to minimise the impact of such local factors when deriving regional climate signals (Comas-Bru et al 2019, Parker et al 2021).

In this paper, we use speleothem oxygen isotope (δ¹⁸O) records from the second version of the Speleothem Isotopes Synthesis and AnaLysis (SISAL) database (Atsawawaranunt et al 2018, Comas-Bru et al 2020a, 2020b) to reconstruct changes in IADV in monsoon regions through the Holocene. We first investigate the potential influence of confounding factors on δ¹⁸O variability using multiple regression. We then reconstruct changes in the amplitude of IADV through the Holocene using the standard deviation (s.d.) of δ¹⁸O values. We construct regional-scale composite changes in IADV through time for individual monsoon regions and compare these to trends in IADV shown by four transient climate model simulations. Previous data-model comparisons of Holocene IADV have focused on ENSO variability (e.g. Carré et al 2021) but have not examined the wider implications for monsoon precipitation IADV. Although we do not associate the changes in IADV with specific climate modes, we compare regional monsoon IADV changes to changes in summer insolation and long-term trends in monsoon precipitation to investigate the relationship with external climate forcing.
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Spatial distribution of speleothem records used in this study, shown with the ensemble mean summer precipitation (mm/d) across the IPSLCM5A-LR, IPSL-Sr, AWI-ESM2 and MPI-ESM1.2 simulations, for 1800–1850. Summer is defined as May to September in the northern hemisphere and November to March in the southern hemisphere. Boxes show the limits used to extract regional precipitation from the Earth System models: SAM = South American Monsoon (latitude: −10° to 0°; longitude: −80° to −64° and latitude: −24° to −10°; longitude: −68° to −40°), ISM = Indian Summer Monsoon (latitude: 10° to 35°; longitude: 75° to 95°) and EAM = East Asian Monsoon (latitude: 20° to 43°; longitude: 100° to 127°).</figcaption></figure>

## 2. Methods

### 2.1. Speleothem oxygen isotope data

We use speleothem $\delta^{18}\text{O}$ data from the SISAL database (Atsawawaranunt et al 2018, Comas-Bru et al 2020a, 2020b), selected using the following criteria:

* The record is located in a monsoon region, between 40 °N and 35 °S;

* The mineralogy is known and does not vary because oxygen isotope fractionation depends on speleothem mineralogy, introducing a non-climatic influence on $\delta^{18}\text{O}$ s.d.;

* The record covers at least 3 000 years of the Holocene;

* The record has a mean sampling resolution of at least 20 years;

This resulted in the selection of 144 records from 79 sites (figure 1). We interpret speleothem $\delta^{18}\text{O}$ records ($\delta^{18}\text{O}_{\text{spel}}$) as changes in the strength of the summer monsoon, via regional precipitation and moisture source pathway and source changes (including oceanic versus land-derived moisture), which likely reinforce one another to give more negative $\delta^{18}\text{O}$ when the summer monsoon is stronger (Vuille et al 2003, Yang et al 2016, Cai et al 2017, Kathayat et al 2021).

### 2.2. Investigating the impact of potential confounding factors

The use of $\delta^{18}\text{O}_{\text{spel}}$ variability as an indicator of IADV in regional precipitation relies on showing that other factors have a minimal impact on this variability. Smoothing of the $\delta^{18}\text{O}$ signals occurs when groundwater of different ages mixes in the epikarst above the cave (Lachniet 2009). Differences in precipitation amount through time could influence the time taken for water to reach the cave, with epikarst transmission times being faster in wetter climates. Apparent $\delta^{18}\text{O}_{\text{spel}}$ variability could also be affected by changes in sampling frequency across the record and changes in speleothem growth rate. We used multiple linear regression (MLR) to investigate if confounding factors such as mean climate, speleothem growth rate or sampling frequency show a significant relationship with $\delta^{18}\text{O}_{\text{spel}}$ variability. The s.d. of $\delta^{18}\text{O}_{\text{spel}}$ values, calculated for a sliding 100-year window (with 50% overlap) across individual records, was used to represent $\delta^{18}\text{O}$ interannual to interdecadal variability (IADV$_{spel}$). We used mean growth rate, number of samples (n-samples) and mean $\delta^{18}\text{O}_{\text{spel}}$ (as an index of mean climate) for each window as predictors. Mean growth rate was calculated as the gradient of the relationship between distance along a speleothem and age. Initial analyses showed that growth rate and n-samples are correlated ($R^2$ = 0.33), whereas mean $\delta^{18}\text{O}$ is not strongly correlated with either growth rate ($R^2$ = −0.04) or n-samples ($R^2$ = 0.07). We therefore constructed two separate MLR models, one using growth rate and mean $\delta^{18}\text{O}$ and the other using n-samples and mean $\delta^{18}\text{O}$ as predictors. Northern and southern hemisphere (NH, SH) records were analysed separately since the trajectory of the mean climate differs between the NH and SH as a result of differences in insolation forcing. The use of multiple geographically separated records in the MLR model captures the influence of confounding factors on composite regional IADV reconstructions. We examined whether variables were normally distributed using the Shapiro-Wilk test and log transformed non-normal variables. Examination of the residuals from the fitted values showed no pattern and confirms that an
<figure><img src="image_7.png" /></figure>

MLR is appropriate. We used P-values to assess the significance of each predictor and R² to determine the influence on δ¹⁸O<sub>spel</sub> variability. Partial residual plots were used to show the relationship between each predictor variable and IADV<sub>spel</sub> when other variables were held constant.

## 2.3. Regional-scale IADV amplitude evolution

We reconstructed IADV amplitude through the Holocene for individual monsoon regions by calculating IADV<sub>spel</sub> for a 100-year sliding window with 50% overlap (as described in section 2.2). We corrected for the influence of confounding factors by subtracting the MLR slope and intercept values from the IADV<sub>spel</sub> values. Regional composites were constructed by fitting a running median and interquartile range through the IADV<sub>spel</sub> values. These composites represent the evolution of the overall amplitude of IADV represented by numerous sites in each region.

We analysed IADV<sub>spel</sub> changes through the Holocene using breakpoint analysis to distinguish intervals with coherent behaviour (Bai and Perron 2003). Breakpoint analysis was carried out using the strucchange package in R (Zeileis et al 2002, Zeileis et al 2003). We identified breakpoints in mean and slope to identify periods with significantly higher or lower IADV (changes in mean) and whether IADV was stable or showed a trend (changes in slope) during these intervals. We determined whether the IADV<sub>spel</sub> of each segment was significantly different using t-tests.

## 2.4. Relationship between IADV and monsoon long-term evolution

The relationship between IADV and mean climate was investigated by comparing the long-term changes in δ¹⁸O<sub>spel</sub> with the IADV<sub>spel</sub> for each region through the Holocene. Long-term changes were calculated by converting δ¹⁸O<sub>spel</sub> to z-scores, with mean and variance standardised relative to a base period of 3 to 5 ka, then calculating the mean δ¹⁸O<sub>spel</sub> z-score for each 100-year bin (across each site) to standardise sampling resolution. Mean z-scores were calculated for 100-year windows for each region and a 500-year half-window locally weighted regression (Cleveland and Devlin 1988) was constructed to emphasise multi-millennial scale variability.

## 2.5. Simulated monsoon IADV evolution

We use transient simulations from four climate models: version 2 of the AWI (Alfred Wegener Institute) Earth System model (Sidorenko et al 2019), version 1.2 of the MPI (Max Planck Institute) Earth System model (Dallmeyer et al 2020) and two versions of the IPSL (Institut Pierre Simon Laplace) Earth system model. The first IPSL simulation uses the IPSLCM5A-LR version, with prescribed vegetation, whilst the second simulation (here termed IPSL-Sr) uses a modified version of IPSLCM5A with a dynamical vegetation module (Dufresne et al 2013, Braconnot et al 2019b). Simulated IADV from IPSLCM5A-LR and IPSL-Sr have been discussed in Braconnot et al (2019a) and Crétat et al (2020). The IPSL and AWI simulations were run from 6 ka to 1950 CE; the MPI simulation from 7.95 ka to 1850 CE. All simulations were spun up using boundary conditions appropriate for the start year. Appropriate orbital parameters (Berger 1978) and greenhouse gases were prescribed annually through the transient runs (See Supplementary Information). The use of an ensemble of four different model simulations provides a way of identifying signals that are consistent across the models (Carré et al 2021). Model performance has been evaluated for these simulations previously, considering climate mean state and variability (Dufresne et al 2013, Duvel et al 2013, Hourdin et al 2013, Mauritsen et al 2019, Sidorenko et al 2019, Braconnot et al 2019b). Despite common model biases (e.g. Tian and Dong, 2020, Good et al, 2021), these comparisons indicate that they capture the tropical climate variability.

We calculated area-averaged summer precipitation over the regional monsoons for each simulation using the region limits in figure 1, where summer was defined as May to September for the NH and November to March for the SH (Wang and Ding 2008). We calculated the s.d. of area-averaged precipitation for 100-year sliding windows with 50% overlap to determine the changes in IADV (termed IADV<sub>precip</sub>). We analysed the overall evolution of IADV<sub>precip</sub> from the mid-Holocene onwards using linear regression. We compared changes in IADV<sub>precip</sub> to IADV<sub>spel</sub> by imposing the breakpoints detected from the speleothem records for each region and comparing whether coeval changes in slope and mean were similar. We calculated the relationship between IADV<sub>precip</sub> and long-term changes by correlating IADV against mean summer regional precipitation for 100-year windows.

# 3. Results

## 3.1. Influence of confounding factors on δ¹⁸O<sub>spel</sub> variability

The MLR model using growth rate and mean δ¹⁸O and the alternative model, using n-samples and mean δ¹⁸O, produced similar results. The models show a significant (P < 0.001) relationship between IADV<sub>spel</sub> and both
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 2.** Partial residual plots from multiple linear regression analysis, showing the relationship between the standard deviation of speleothem $\delta^{18}$O ($\delta^{18}$Ospel s.d.) and two predictor variables: growth rate and mean $\delta^{18}$O (as an indicator of mean climate). Speleothem data is separated into northern hemisphere (top panel) and southern hemisphere (bottom panel) sites. All relationships are significant with P values < 0.001.</figcaption></figure>

<table><caption>**Table 1.** Results of the multiple linear regression analysis, separated into northern hemisphere (NH) and southern hemisphere (SH) sites. Significant relationships (P <0.01) are shown in bold.</caption><thead><tr><th></th><th>NH</th><th>SH</th></tr></thead><tbody><tr><td>Growth rate</td><td><strong>0.08</strong></td><td><strong>0.23</strong></td></tr><tr><td>Mean δ<sup>18</sup>O</td><td><strong>-0.03</strong></td><td><strong>0.08</strong></td></tr></tbody></table>

growth rate (or n-samples) and mean $\delta^{18}$O<sub>spel</sub> (table 1, table S1 available online at stacks.iop.org/ERC/3/121002/mmedia). Growth rate has a significant positive relationship with IADV<sub>spel</sub> (figures 2, S1). There is a significant negative relationship in the NH and a significant positive relationship in the SH between IADV<sub>spel</sub> and mean $\delta^{18}$O<sub>spel</sub> (figure 2). The models have low R<sup>2</sup> (NH: 0.06; SH: 0.27) suggesting that confounding factors only have a limited influence on IADV<sub>spel</sub>.

### 3.2. Regional $\delta^{18}$O<sub>spel</sub> variability evolution

Three regions have sufficient data to reconstruct composites of IADV<sub>spel</sub>: the Indian Monsoon (ISM), East Asian Monsoon (EAM) and South American Monsoon (SAM) (figure 1). The Central American and Indonesian-Australian monsoon regions are excluded because they contain too few records of a sufficient length. The regional evolution of corrected and uncorrected IADV<sub>spel</sub> is very similar (figure 3), confirming that changes in growth rate (or n-samples) and mean $\delta^{18}$O<sub>spel</sub> are not the main cause of the reconstructed changes in $\delta^{18}$O<sub>spel</sub> variability during the Holocene. Any uncertainty arising from combining records of different lengths is small, based on bootstrap resampling (figure S2). Regional composites show a significant amount of centennial-scale variability, nevertheless broad trends can be seen, with differing IADV<sub>spel</sub> evolution between the three regions.

The EAM IADV<sub>spel</sub> record is decomposed into four segments (figure 4(a), table S2): The early Holocene interval (12 to 8.9 ka) has significantly higher IADV<sub>spel</sub> than any subsequent period. The following interval between 8.9 and 7.1 ka has significantly lower IADV<sub>spel</sub> than the early Holocene interval. The mid-Holocene interval (7.1 to 2 ka) shows a significant increase in IADV<sub>spel</sub>. The last 2 000 years has significantly lower IADV<sub>spel</sub> than any preceding interval. The ISM IADV<sub>spel</sub> record is decomposed into three segments (figure 4(b), table S2): the early Holocene (12 to 9 ka) has significantly lower IADV<sub>spel</sub> values than the other two periods. The mid-Holocene period (9 to 1.7 ka) has significantly lower IADV<sub>spel</sub> than the period after 1.7 ka. Breakpoint analysis of
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 3.** Evolution of regional speleothem $\delta^{18}\text{O}$ standard deviation ($\delta^{18}\text{O}_{\text{spel}}$ s.d.) values through the Holocene. Trends are shown with uncorrected values and values with a correction applied, using the relationships constrained in the multiple linear regression model. Regional monsoons are EAM = East Asian Monsoon, ISM = Indian Summer Monsoon, SAM = South American Monsoon. $\delta^{18}\text{O}_{\text{spel}}$ s.d. values are extracted using a 100-year windows with 50% overlap. Regional composites are expressed as a running median with confidence intervals represented by 5% and 95% quantiles.</figcaption></figure>

the slopes (figure S3, table S3) show stable $\text{IADV}_{\text{spel}}$ values, except between 2.05 and 0 ka when they show a significant increase over time. The SAM $\text{IADV}_{\text{spel}}$ record is decomposed into four segments (figure 4(c), table 3). The highest values of $\text{IADV}_{\text{spel}}$ occur in the mid-Holocene (6 to 4.2 ka) and the early Holocene (12 to 10.15 ka). The interval between these two peaks is characterised by $\text{IADV}_{\text{spel}}$ values that are statistically indistinguishable from the past 4 000 years. Although the recent interval shows stable variability (figure S3, table S3), the earlier part of the SAM record shows increases in variability with rapid transitions to lower variability at 10.15, 7.9 and 4.75 ka.

### 3.3. Comparison with simulated monsoon variability

There is no overall change in $\text{IADV}_{\text{spel}}$ over the last 6 000 years in the EAM (table 2), reproduced by the AWI and MPI simulations (table 2). The IPSL simulations show a small but significant increase in $\text{IADV}_{\text{precip}}$ through time. The significant increase in $\text{IADV}_{\text{spel}}$ from 6 onwards in the ISM is reproduced by all four simulations (figure 5, table 2). There is no significant change in $\text{IADV}_{\text{precip}}$ through time in the SAM in any simulation (figure 5, table 2), whereas $\text{IADV}_{\text{spel}}$ increases significantly through time. However, this is driven by the peak at 4.2 to 6 ka; a linear regression from 5 ka onwards shows no significant change.

There are significant changes in simulated $\text{IADV}_{\text{precip}}$ between the segments identified from the speleothems (table 4). In the ISM, all four simulations show significantly higher $\text{IADV}_{\text{precip}}$ in the 1.7 to 0 ka period than the preceding period (table 4), consistent with the $\text{IADV}_{\text{spel}}$ record. All simulations show a significant increase in $\text{IADV}_{\text{precip}}$ values from 1.7 ka onwards, but no significant trend prior to this (Table S4). This is consistent with the marked increase in $\text{IADV}_{\text{spel}}$ during the most recent 1 700 years. In the EAM, the AWI and MPI simulations show no significant difference in either slope or magnitude between the two intervals. (tables 4, S4), whilst the IPSL simulations show higher $\text{IADV}_{\text{precip}}$ values for the most recent 2 000 years compared to the prior period. There is no significant difference between the intervals before and after 4.2 ka in SAM $\text{IADV}_{\text{precip}}$ and neither interval shows significant changes in IADV through time (tables 4, S4). This is not consistent with the speleothem record, which shows higher variability before 4.2 ka than after.

### 3.4. Relationship between monsoon variability and long-term evolution

Speleothems and simulations show no consistent relationship between IADV and long-term changes in monsoon precipitation. There is declining summer monsoon strength and regional precipitation from the mid-Holocene onwards in the EAM and ISM, with a slower decline from ~ 2 ka, following NH summer insolation (figure 6). However, IADV evolution differs between these two regions. In the EAM, the decreasing monsoon strength is mirrored by a decrease in $\text{IADV}_{\text{spel}}$ (figure 7(a), table 5). This relationship is reproduced by the AWI
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 4.** Breakpoint models ($IADV_{spel}$ ~ 1) of regional $IADV_{spel}$ evolution. Breakpoints are shown with their 2.5 and 97.5% confidence intervals. Step changes in $IADV_{spel}$ are shown. Regional monsoons are EAM = East Asian Monsoon, ISM = Indian Summer Monsoon, SAM = South American Monsoon.</figcaption></figure>

<table><caption>**Table 2.** Linear regression relationships of regional $\delta^{18}O_{spel}$ s.d. versus time and simulated precipitation s.d. trends versus time (mm/month yr$^{-1}$) from the mid-Holocene (6 000 years BP) to 0 years BP. Negative relationships represent lower variability (s.d.) at the mid-Holocene than at present. Significant relationships (P < 0.0001) are shown in bold. EAM = East Asian monsoon, ISM = Indian Summer monsoon, SAM = South American monsoon.</caption><thead><tr><th></th><th>spel 6–0ka</th><th>IPSLCM5A-LR</th><th>IPSL-Sr</th><th>AWI-ESM2</th><th>MPI-ESM1.2</th></tr></thead><tbody><tr><th>EAM</th><td>0.015</td><td><strong>-1.17E-05</strong></td><td><strong>-1.01E-05</strong></td><td>9.96E-07</td><td>-2.34E-07</td></tr><tr><th>ISM</th><td><strong>-0.061</strong></td><td><strong>-1.22E-05</strong></td><td><strong>-2.51E-05</strong></td><td><strong>-1.96E-05</strong></td><td><strong>-1.26E-05</strong></td></tr><tr><th>SAM</th><td>0.020</td><td>1.49E-07</td><td>-4.85E-06</td><td><strong>-6.40E-06</strong></td><td>-1.89E-06</td></tr></tbody></table>

and MPI simulations, although only MPI shows a significant relationship between mean precipitation and $IADV$. In contrast, the declining summer monsoon strength from the mid-Holocene onwards is accompanied by increasing $IADV$ amplitude in the ISM. This relationship is reproduced by all four simulations. SAM observations and simulations show a strengthening monsoon from the mid-Holocene onwards (figure 6(e)), mirroring SH summer insolation (figure 6(d)). There is no significant relationship between the long-term changes and $IADV_{spel}$ in the SAM. The lack of a relationship is reproduced by three of the four simulations (table 5).

## 4. Discussion and conclusions

We have shown that variable growth rate, sampling frequency and mean climate have a minimal impact on $\delta^{18}O_{spel}$ variability overall. Although cave monitoring studies have emphasised the influence of site-specific
<figure><img src="image_14.png" /><figcaption>**Table 3.** P-values for pairwise t-test, with the null hypothesis that the mean difference between two segments is zero. Segments are separated for each regional monsoon by breakpoint analysis, where the timings of abrupt changes ($IADV_{spel}$ ~ 1) were determined. Regional monsoons are: EAM, East Asian monsoon, ISM, Indian summer monsoon, SAM, South American monsoon. Bold values indicate P-values <= 0.001, where the null hypothesis can be rejected.</figcaption></figure>

<table><thead><tr><th>Segments (years BP)</th><th>Mean s.d.</th><th>P-values</th><th></th><th></th><th></th></tr></thead><tbody><tr><td>EAM</td><td></td><td>0–2000</td><td>2000–5950</td><td>5950–8800</td><td>8800–12000</td></tr><tr><td>0–2000</td><td>−0.049</td><td>/</td><td></td><td></td><td></td></tr><tr><td>2000–7100</td><td>0.015</td><td><strong>7.88E-07</strong></td><td>/</td><td></td><td></td></tr><tr><td>7100–8900</td><td>−0.045</td><td>7.78E-01</td><td><strong>6.13E-06</strong></td><td>/</td><td></td></tr><tr><td>8900–12000</td><td>0.069</td><td><strong>5.20E-16</strong></td><td><strong>2.28E-06</strong></td><td><strong>3.44E-14</strong></td><td>/</td></tr><tr><td>ISM</td><td></td><td>0–2150</td><td>2150–9000</td><td>9000–12000</td><td></td></tr><tr><td>0–1700</td><td>0.198</td><td>/</td><td></td><td></td><td></td></tr><tr><td>1700–9000</td><td>0.022</td><td><strong>1.10E-17</strong></td><td>/</td><td></td><td></td></tr><tr><td>9000–12000</td><td>−0.024</td><td><strong>2.21E-20</strong></td><td><strong>4.62E-03</strong></td><td>/</td><td></td></tr><tr><td>SAM</td><td></td><td>0–4200</td><td>4200–7900</td><td>7900–10150</td><td>10150–12000</td></tr><tr><td>0–4200</td><td>0.051</td><td>/</td><td></td><td></td><td></td></tr><tr><td>4200–6000</td><td>0.1942</td><td><strong>1.40E-06</strong></td><td>/</td><td></td><td></td></tr><tr><td>6000–10150</td><td>0.046</td><td>8.20E-01</td><td><strong>8.83E-07</strong></td><td>/</td><td></td></tr><tr><td>10150–12000</td><td>0.234</td><td><strong>1.04E-09</strong></td><td>1.49E-01</td><td><strong>8.65E-11</strong></td><td>/</td></tr></tbody></table>

factors on $\delta^{18}O_{spel}$ (Luo et al 2014, Wu et al 2014, Pu et al 2016), these factors are unimportant when records from different speleothems and cave sites are combined. This is expected because, unlike climate signals, site-specific factors that might affect growth rate are unlikely to be consistent across sites. It might be useful to take account of the residence time of water in the epikarst directly in reconstructing $\delta^{18}O_{spel}$ variability, but this is hard to quantify under the changing climate conditions of the geological past and may be unnecessary given that climate conditions that influence residence time may be indirectly incorporated via mean $\delta^{18}O_{spel}$. Despite the simplicity of our MLR model, the relationships between $IADV_{spel}$ and the predictor variables are physically plausible and corrected $IADV_{spel}$ evolution is similar between speleothems and sites.

Regional $IADV_{spel}$ changes cannot be compared with the Galapagos Islands and Ecuador lake records which are interpreted as ENSO frequency rather than IADV amplitude. However, the hydrogen isotope record from El Junco (Zhang et al 2014), which is interpreted as a record of El Niño amplitude, shows large multi-centennial scale fluctuations in the early to mid-Holocene like those in SAM $IADV_{spel}$. Indeed, numerous records show fluctuations of a similar nature, including tropical Pacific reconstructions (Carré et al, 2021) and the Lake Pallcacocha records (Moy et al 2002). The low amplitude variability recorded at El Junco during the mid-Holocene and the abrupt shift to present-day variability at ~3.5 ka is also similar to the ISM $IADV_{spel}$ record. ISM and SAM $IADV_{spel}$ shows similarities with regional reconstructions of IADV amplitude from the tropical Pacific (Carré et al 2021). The western Pacific shows lower IADV through the early to mid-Holocene, then an abrupt shift to modern levels at ~2 ka, similar to ISM $IADV_{spel}$, whilst central Pacific records show a peak in IADV at ~6 ka, similar to the peak in SAM $IADV_{spel}$ at this time. Thus, the speleothem composites provide a coherent picture of regional changes in IADV that are broadly consistent with other Holocene IADV reconstructions.

Some features of the $IADV_{spel}$ records are reproduced by the models. Simulated IADV evolution reproduces the higher IADV of the last 1 700 years compared to the mid-Holocene seen in ISM $IADV_{spel}$ and the lack of an overall long-term Holocene changes in $IADV_{spel}$ in SAM. The AWI and MPI simulations show differing trends between the ISM and EAM, as seen in the $IADV_{spel}$ records. There are, however, differences between speleothem and simulated IADV changes. The IPSL simulations show increasing IADV from the mid-Holocene to present for the EAM, whereas $IADV_{spel}$ decreases. However, these simulations show smaller changes in IADV in the EAM than the ISM and thus there is a difference in ISM and EAM simulated IADV even though this differentiation is less marked than the opposite changes seen in regional $IADV_{spel}$. The models do not show the large multi-centennial scale fluctuations seen in SAM $IADV_{spel}$. Comparisons of simulated and reconstructed SST IADV in the tropical Pacific also show that these models exhibit less multi-centennial variability than coral and bivalve records (Carré et al 2021). These discrepancies between simulated and reconstructed IADV, and indeed between different models, strongly suggests that the controls on monsoon precipitation are only indirectly tied to changes in external forcing.

The long-term evolution of regional precipitation follows summer insolation (figure 6), but changes in regional monsoon IADV do not show a direct relationship with insolation. Long-term precipitation changes in the EAM and ISM are similar, but IADV evolution differs and even the abrupt change in IADV at ~ 2 ka has a different sign in the two regions. Analyses of marine records from the Pacific (Carré et al 2021) suggests that
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 5.** Evolution of regional speleothem $\delta^{18}$O standard deviation ($\delta^{18}$O<sub>spe</sub> s.d.) values through the Holocene, compared to regional precipitation standard deviation simulated by four transient Earth System Models. Regional monsoons are: (a) EAM, East Asian Monsoon, (b) ISM, Indian Summer Monsoon, (c) SAM, South American Monsoon. Precipitation and $\delta^{18}$O<sub>spe</sub> s.d. values are extracted for 100-year windows with 50% overlap. Precipitation s.d. values are calculated as anomalies to the last 1 000 years. Trends are fitted with a smoothed loess fit (500-year window).</figcaption></figure>

<table><caption>**Table 4.** Mean regional precipitation s.d. values for Holocene segments, identified in the breakpoint analysis of the speleothem trends. Bold values show segments (within an individual region and model) that are significantly different from the other (P < 0.0001), determined using pairwise t-test. Regional monsoons are: EAM, East Asian monsoon, ISM, Indian summer monsoon and SAM, South American monsoon.</caption><thead><tr><th>Region</th><th>Segments</th><th>IPSLCM5A-LR</th><th>IPSL-Sr</th><th>AWI-ESM2</th><th>MPI-ESM1.2</th></tr></thead><tbody><tr><td rowspan="2">EAM</td><td>0–2000</td><td><strong>0.310</strong></td><td><strong>0.357</strong></td><td>0.274</td><td>0.271</td></tr><tr><td>2000–6000</td><td><strong>0.274</strong></td><td><strong>0.327</strong></td><td>0.271</td><td>0.277</td></tr><tr><td rowspan="2">ISM</td><td>0–1700</td><td><strong>0.355</strong></td><td><strong>0.465</strong></td><td><strong>0.479</strong></td><td><strong>0.413</strong></td></tr><tr><td>1700–6000</td><td><strong>0.313</strong></td><td><strong>0.389</strong></td><td><strong>0.418</strong></td><td><strong>0.370</strong></td></tr><tr><td rowspan="2">SAM</td><td>0–4200</td><td>0.442</td><td>0.448</td><td>0.277</td><td>0.305</td></tr><tr><td>4200–6000</td><td>0.441</td><td>0.434</td><td>0.255</td><td>0.293</td></tr></tbody></table>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 6.** Comparison of summer insolation with long-term evolution of monsoon mean state through the Holocene. Regional speleothem $\delta^{18}$O trends are represented by smoothed (500-year half-window loess fit) z-scores with respect to a base period of 3 000 to 5 000 years and 100 year mean z-score values (dark red line). Simulated monsoon evolution is represented by region-averaged summer precipitation (May to September for Northern Hemisphere monsoons, November to March for Southern hemisphere monsoons), as anomalies to 0–1 000 years BP. Insolation values are calculated as summer means (May to September for Northern hemisphere monsoons, November to March for Southern hemisphere monsoons). Regional monsoons are b) East Asian monsoon (EAM), c) Indian summer monsoon (ISM) and e) Southern American monsoon (SAM). Note that the y-axes for $\delta^{18}$O$_{spel}$ z-scores have been reversed to match the negative relationship of $\delta^{18}$O$_{spel}$ with monsoon strength and summer insolation.</figcaption></figure>

there is a link between orbital forcing and ENSO variability, but ENSO is not the only cause of changes in rainfall IADV in monsoon regions.

Differences in the relationship between simulated long-term evolution and precipitation IADV between the Indian and west African monsoons have been explained as reflecting the important role of ENSO over India and Atlantic modes over west Africa for precipitation variability (Braconnot et al 2019a). Differences in regional IADV$_{spel}$ evolution presumably also reflects the importance of different modes for IADV in each region. In the present-day, IADV over the EAM is influenced by ENSO signals that may also be modulated by the Pacific Decadal Oscillation (Feng et al 2014, Zhang et al 2018), whilst the ISM may also be influenced by changes in the Indian ocean (Ashok et al 2004, Krishnaswamy et al 2015). SAM IADV is influenced by changes in the tropical Atlantic and Indian Ocean, as well as ENSO (Grimm and Zilli 2009, Yoon and Zeng 2010, Jimenez et al 2021). Furthermore, each monsoon region has a distinct relationship with ENSO, with the nature and strength of the precipitation responses depending strongly on the location of SST anomalies in the tropical Pacific (Kumar et al 2006, Tedeschi et al 2015, Zhang et al 2016). All these regions show a changing relationship with ENSO in recent decades (Samanta et al 2020, Seetha et al 2019, Zhang et al 2018, Wang et al 2020) and this was likely the case during the Holocene. In the ISM region, for example, a weakening of the monsoon precipitation-ENSO relationship in recent decades has been attributed to a warming Eurasian continent (Kumar et al 1999), changes in tropical Atlantic SSTs (Kucharski et al 2007) and influence from the IOD (Ashok et al 2001). Analysis of ISM IADV evolution in the IPSL simulations (Crétat et al 2020) emphasised the interaction between ENSO, IOD and ISM precipitation variability. In the East Asian monsoon, the observed changing relationship between monsoon precipitation and ENSO over the last millennium is attributed to interaction with the PDO (Zhang et al 2018). Reconstructions and model simulations show changing PDO-like variability through the Holocene, including a positive phase during the early Holocene (Chen et al 2021). It is therefore likely that the relationship between
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 7.** Relationships between the long-term evolution of the regional monsoons (x-axes) and their variability (y-axes). We use regional speleothem $\delta^{18}\text{O}$ z-score data and mean precipitation for 100 year sliding windows to represent long-term changes in monsoons for speleothem data and modelled precipitation respectively. Variability is represented by the standard deviation of $\delta^{18}\text{O}_{\text{spel}}$ and precipitation for 100 year sliding windows. Regions are a) EAM, East Asian monsoon, b) ISM, Indian Summer monsoon and c) SAM, South American monsoon. For each region, we show the results from four climate model simulations. Note that the axes for $\delta^{18}\text{O}_{\text{spel}}$ z-scores have been reversed to match the negative relationship between $\delta^{18}\text{O}_{\text{spel}}$ and monsoon strength.</figcaption></figure>

<table><caption>**Table 5.** Linear regression relationships of s.d. values ($\delta^{18}\text{O}_{\text{spel}}$ s.d., regional precipitation s.d.) versus long-term evolution ($\delta^{18}\text{O}_{\text{spel}}$ z-score, mean 100-year regional precipitation). Significant relationships (P < 0.0001) are shown in bold. Regions are EAM = East Asian monsoon, ISM = Indian Summer monsoon, SAM = South American monsoon. Note that more negative $\delta^{18}\text{O}$ z-scores and higher regional precipitation values represent a stronger monsoon, therefore agreement between speleothem data and models will be opposite direction trends.</caption><thead><tr><th>Region</th><th>Spel</th><th>IPSLCM5A-LR</th><th>IPSL-Sr</th><th>AWI-ESM2</th><th>MPI-ESM1.2</th></tr></thead><tbody><tr><td>EAM</td><td><strong>-0.0245</strong></td><td><strong>-0.102</strong></td><td><strong>-0.055</strong></td><td>0.016</td><td><strong>0.027</strong></td></tr><tr><td>ISM</td><td><strong>0.0556</strong></td><td><strong>-0.083</strong></td><td><strong>-0.178</strong></td><td><strong>-0.105</strong></td><td><strong>-0.051</strong></td></tr><tr><td>SAM</td><td>0.105</td><td>0.004</td><td>0.03</td><td><strong>0.059</strong></td><td>0.017</td></tr></tbody></table>

EAM precipitation variability and Pacific SSTs was significantly different during this period. Since regional monsoon precipitation variability is influenced by numerous modes that can interact with one another, it is likely that the lack of a direct and linear relationship between monsoon IADV and insolation through the Holocene reflects the complexity of these influences on regional precipitation and the potential for changes in the importance of different climate modes through time. More emphasis should therefore be placed on modes of variability other than ENSO to improve our understanding of the evolution of monsoon climate during the Holocene, using high resolution terrestrial records and model simulations.

## Acknowledgments

SP and SPH acknowledge support from the ERC-funded project GC2.0 (Global Change 2.0: Unlocking the past for a clearer future, grant number 694481). SPH and PB acknowledge support from JPI-Belmont Forum project
<figure><img src="image_21.png" /></figure>

Palaeo-climate Constraints on Monsoon Evolution and Dynamics (PACMEDY). We acknowledge the modelling groups of Institut Pierre Simon Laplace, Alfred Wegener Institute and Max Planck Institute for producing and sharing the climate simulations. We thank Laia Comas-Bru for helpful discussion on the speleothem analysis.

## Data availability statement

The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.5281/zenodo.5084426 (Parker 2021).

## Data and code availability

The code used to generate the results and figures in this study are available at (doi.org/10.5281/zenodo.5084426; Parker 2021). SISAL version 2 is available at http://doi.org/10.17864/1947.256. We use R for analyses (R Core Team 2019). The main packages used are tidyverse (Wickham et al 2021), ggplot2 (Wickham 2016), RMySQL (Ooms et al 2020) and ncdf4 (Pierce 2019).

## ORCID iDs

Sarah E Parker https://orcid.org/0000-0002-2467-2624

Sandy P Harrison https://orcid.org/0000-0001-5687-1903

## References

Ashok K, Guan Z, Saji N H and Yamagata T 2004 Individual and combined influences of ENSO and the Indian Ocean dipole on the Indian summer monsoon *J. Clim.* **17** 3141–55

Ashok K, Guan Z and Yamagata T 2001 Impact of the Indian Ocean dipole on the relationship between the Indian monsoon rainfall and ENSO *Geophys. Res. Lett.* **28** 4499–502

Atsawawanunt K et al 2018 The SISAL database: A global resource to document oxygen and carbon isotope records from speleothems *Earth System Science Data*, **10** 1687–713

Bai J and Perron P 2003 Computation and analysis of multiple structural change models *Journal of Applied Econometrics* **18** 1–22

Berger A 1978 Long-Term Variations of Daily Insolation and Quaternary Climatic Changes *Journal of the Atmospheric Sciences* **35** 2362–7

Bouma M J and Kaay H J 1996 The El Niño Southern Oscillation and the historic malaria epidemics on the Indian subcontinent and Sri Lanka: an early warning system for future epidemics? *Tropical Medicine & International Health* **1** 86–96

Braconnot P, Crétat J, Marti O, Balkanski Y, Caubel A and Cozic A 2019a Impact of multiscale variability on last 6000 years Indian and West African monsoon rain *Geophys. Res. Lett.* **46** 14021–9

Braconnot P, Zhu D, Marti O and Servonnat J 2019b Strengths and challenges for transient Mid- to Late Holocene simulations with dynamical vegetation *Climate of the Past* **15** 997–1024

Brown J R et al 2020 Comparison of past and future simulations of ENSO in CMIP5/PMIP3 and CMIP6/PMIP4 models *Climate of the Past*, **16** 1777–805

Cai Z, Tian L and Bowen G J 2017 ENSO variability reflected in precipitation oxygen isotopes across the Asian Summer Monsoon region *Earth Planet. Sci. Lett.* **475** 25–33

Carré M et al 2021 High-resolution marine data and transient simulations support orbital forcing of ENSO amplitude since the mid-Holocene *Quat. Sci. Rev.* **268** 107125

Chen C, Zhao W and Zhang X 2021 Pacific decadal oscillation-like variability at a millennial timescale during the holocene *Global Planet. Change* **199** 103448

Chen L, Li T, Yu Y and Behera S K 2017 A possible explanation for the divergent projection of ENSO amplitude change under global warming *Clim. Dyn.* **49** 3799–811

Cheng H, Zhang H, Zhao J, Li H, Ning Y and Kathayat G 2019 Chinese stalagmite paleoclimate researches: A review and perspective *Science China Earth Sciences* **62** 1489–513

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

Comas-Bru L, Atsawawanunt K, Harrison S P and SISAL working group members (2020a) SISAL (Speleothem Isotopes Synthesis and AnaLysis Working group) Database Version 2.0. University of Reading (https://doi.org/10.17864/1947.256)

Comas-Bru L, Harrison S P, Werner M, Rehfeld K, Scroxton N and Veiga-Pires C 2019 Evaluating model outputs using integrated global speleothem *Climate of the Past Discussions*. **15** 1557–79

Comas-Bru L et al 2020b SISALv2: a comprehensive speleothem isotope database with multiple age–depth models *Earth System Science Data* **12** 2579–606

Conroy J L, Overpeck J T, Cole J E, Shanahan T M and Steinitz-Kannan M 2008 Holocene changes in eastern tropical Pacific climate inferred from a Galápagos lake sediment record *Quat. Sci. Rev.* **27** 1166–80

Crétat J, Braconnot P, Terray P, Marti O and Falasca F 2020 Mid-Holocene to present-day evolution of the Indian monsoon in transient global simulations *Clim. Dyn.* **55** 2761–84

Dallmeyer A, Claussen M, Lorenz S J and Shanahan T 2020 The end of the African humid period as seen by a transient comprehensive Earth system model simulation of the last 8000 years *Climate of the Past* **16** 117–40

Dufresne J-L et al 2013 Climate change projections using the IPSL-CM5 Earth System Model: from CMIP3 to CMIP5 *Climate Dynamics*, **40** 2123–65
<figure><img src="image_22.png" /></figure>

Duvel J P, Bellenger H, Bellon G and Remaud M 2013 An event-by-event assessment of tropical intraseasonal perturbations for general circulation models *Clim. Dyn.* **40** 857–73

Emile-Geay J and Tingley M 2016 Inferring climate variability from nonlinear proxies: Application to palaeo-ENSO studies *Climate of the Past* **12** 31–50

Feng J, Wang L and Chen W 2014 How does the East Asian summer monsoon behave in the decaying phase of El Niño during different PDO phases? *J. Clim.* **27** 2682–98

Gagnon A, Bush A and Smoyer-Tomic K 2001 Dengue epidemics and the El Niño Southern Oscillation *Climate Research* **19** 35–43

Good P, Chadwick R, Holloway C E, Kennedy J, Lowe J A, Roehrig R and Rushley S S 2021 High sensitivity of tropical precipitation to local sea surface temperature *Nature* **589** 408–14

Grimm A M and Zilli M T 2009 Interannual variability and seasonal evolution of summer monsoon rainfall in South America *J. Clim.* **22** 2257–75

Hashizume M, Chaves L F and Minakawa N 2012 Indian Ocean dipole drives malaria resurgence in East African highlands *Sci. Rep.* **2** 1–6

Hourdin F *et al* 2013 Impact of the LMDZ atmospheric grid configuration on the climate and sensitivity of the IPSL-CM5A coupled model. *Climate Dynamics* **40** 2167–92

Hui C and Zheng X 2018 Uncertainty in Indian Ocean Dipole response to global warming: the role of internal variability *Climate Dynamics* **51** 3597–611

Iizumi T, Luo J-J, Challinor A J, Sakurai G, Yokozawa M, Sakuma H, Brown M E and Yamagata T 2014 Impacts of El Niño Southern Oscillation on the global yields of major crops *Nat. Commun.* **5** 3712

Jimenez J C, Marengo J A, Alves L M, Sulca J C, Takahashi K, Ferrett S and Collins M 2021 The role of ENSO flavours and TNA on recent droughts over Amazon forests and the Northeast Brazil region *Int. J. Climatol.* **41** 3761–80

Kane R P 1999 Some characteristics and precipitation effects of the El Niño of 1997–1998 *J. Atmos. Sol. Terr. Phys.* **61** 1325–46

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

Kiefer J and Karamperidou C 2019 High-resolution modeling of ENSO-induced precipitation in the tropical Andes: Implications for proxy interpretation *Paleoceanography and Paleoclimatology* **34** 217–36

Kirono D G C, Tapper N J and McBride J L 1999 Documenting Indonesian rainfall in the 1997/1998 El Niño event *Physical Geography* **20** 422–35

Krishnan R and Sugi M 2003 Pacific decadal oscillation and variability of the Indian summer monsoon rainfall *Clim. Dyn.* **21** 233–42

Krishnaswamy J, Vaidyanathan S, Rajagopalan B, Bonell M, Sankaran R, Bhalla R S and Badiger S 2015 Non-stationary and non-linear influence of ENSO and Indian Ocean Dipole on the variability of Indian monsoon rainfall and extreme rain events *Clim. Dyn.* **45** 175–84

Kucharski F, Bracco A, Yoo J H and Molteni F 2007 Low-frequency variability of the Indian monsoon–ENSO relationship and the tropical Atlantic: The ‘weakening’ of the 1980s and 1990s *J. Clim.* **20** 4255–66

Kumar K K, Rajagopalan B, Hoerling M, Bates G and Cane M 2006 Unraveling the mystery of Indian monsoon failure during El Niño *Science* **314** 115–9

Kumar K K, Rajagopalan B and Cane M A 1999 On the weakening relationship between the Indian monsoon and ENSO *Science* **284** 2156–9

Lachniet M S 2009 Climatic and environmental controls on speleothem oxygen-isotope values *Quat. Sci. Rev.* **28** 412–32

Luo W, Wang S, Zeng G, Zhu X and Liu W 2014 Daily response of drip water isotopes to precipitation in Liangfeng Cave, Guizhou Province, SW China *Quat. Int.* **349** 153–8

Mauritsen T *et al* 2019 Developments in the MPI-M Earth System Model version 1.2 (MPI-ESM1.2) and Its Response to Increasing CO2 *Journal of Advances in Modeling Earth Systems* **11** 998–1038

McPhaden M J, Zebiak S E and Glantz M H 2006 ENSO as an integrating concept in earth science *Science* **314** 1740–5

Messié M and Chavez F 2011 Global modes of sea surface temperature variability in relation to regional climate indices *J. Clim.* **24** 4314–31

Moy C M, Seltzer G O, Rodbell D T and Anderson D M 2002 Variability of El Niño/Southern Oscillation activity at millennial timescales during the Holocene epoch *Nature* **420** 162–5

Ng B, Cai W, Cowan T and Bi D 2018 Influence of internal climate variability on Indian Ocean Dipole properties *Sci. Rep.* **8** 1–8

Ooms J, James D, DebRoy S, Wickham H and Horner J 2020 RMySQL: Database Interface and ‘MySQL’ Driver for R. R package version 0.10.22. https://CRAN.R-project.org/package=RMySQL

Parker S E 2021 Zenodo Code for ‘Speleothem records of monsoon interannual-interdecadal variability through the Holocene’ [Accessed 28 November 2021] (https://doi.org/10.5281/zenodo.5084426)

Parker S E, Harrison S P, Comas-Bru L, Kaushal N, LeGrande A N and Werner M 2021 A data–model approach to interpreting speleothem oxygen isotope records from monsoon regions *Climate of the Past* **17** 1119–38

Phillips T, Nerem R S, Fox-Kemper B, Famiglietti J S and Rajagopalan B 2012 The influence of ENSO on global terrestrial water storage using GRACE *Geophys. Res. Lett.* **39**

Pierce D 2019 ncdf4: Interface to Unidata netCDF (Version 4 or Earlier) Format Data Files. R package version 1.18. https://CRAN.R-project.org/package=ncdf4

Pu J, Wang A, Shen L, Yin J, Yuan D and Zhao H 2016 Factors controlling the growth rate, carbon and oxygen isotope variation in modern calcite precipitation in a subtropical cave, Southwest China *J. Asian Earth Sci.* **119** 167–78

Ray D K, Gerber J S, MacDonald G K and West P C 2015 Climate variation explains a third of global crop yield variability *Nat. Commun.* **6** 5989

Riedinger M A, Steinitz-Kannan M, Last W M and Brenner M 2002 A ~6100 14C yr record of El Niño activity from the Galápagos Islands *J. Paleolimnol.* **27** 1–7

Rodbell D T 1999 An 15,000-Year record of El Niño-driven alluviation in southwestern Ecuador *Science* **283** 516–20

Saji N and Yamagata T 2003 Possible impacts of Indian Ocean dipole mode events on global climate *Climate Research* **25** 151–69

Samanta D, Rajagopalan B, Karnauskas K B, Zhang L and Goodkin N F 2020 La Niña’s Diminishing Fingerprint on the Central Indian Summer Monsoon *Geophysical Research Letters* **47** e2019GL086237

Schneider T, Hampel H, Mosquera P V, Tylmann W and Grosjean M 2018 Paleo-ENSO revisited: Ecuadorian Lake Pallcacocha does not reveal a conclusive El Niño signal *Global Planet. Change* **168** 54–66

Seetha C J, Varikoden H, Babu C A and Kuttippurath J 2019 Significant changes in the ENSO-monsoon relationship and associated circulation features on multidecadal timescale *Climate Dynamics* **54** 1491–506

Sidorenko D *et al* 2019 Evaluation of FESOM2.0 Coupled to ECHAM6.3: Preindustrial and HighResMIP Simulations *Journal of Advances in Modeling Earth Systems* **11** 3794–815
<figure><img src="image_23.png" /></figure>

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

Stevenson S L 2012 Significant changes to ENSO strength and impacts in the twenty-first century: Results from CMIP5 *Geophys. Res. Lett.* **39** 17

R Core Team 2019 *R: A Language and Environment for Statistical Computing*. (Vienna, Austria)

Tian B and Dong X 2020 The double-ITCZ Bias in CMIP3, CMIP5, and CMIP6 models based on annual mean precipitation *Geophys. Res. Lett.* **47** e2020GL087232

Tedeschi R G, Grimm A M and Cavalcanti I F 2015 Influence of Central and East ENSO on extreme events of precipitation in South America during austral spring and summer *Int. J. Climatol.* **35** 2045–64

Thompson D M et al 2017 Tropical Pacific climate variability over the last 6000 years as recorded in Bainbridge Crater Lake, Galápagos *Paleoceanography* **32** 903–22

Ummenhofer C C, Gupta A S, Li Y, Taschetto A S and England M H 2011 Multi-decadal modulation of the El Nino–Indian monsoon relationship by Indian Ocean variability *Environ. Res. Lett.* **6** 034006

Vuille M, Bradley R S, Werner M, Healy R and Keimig F 2003 Modeling $\delta^{18}$O in precipitation over the tropical Americas: 1. Interannual variability and climatic controls *J. Geophys. Res. D: Atmos.* **108** D6

Wang B and Ding Q 2008 Global monsoon: Dominant mode of annual variation in the tropics *Dyn. Atmos. Oceans* **44** 165–83

Wang B, Liu J, Kim H J, Webster P J and Yim S Y 2012 Recent change of the global monsoon precipitation (1979-2008) *Clim. Dyn.* **39** 1123–35

Wang B, Luo X and Liu J 2020 How robust is the Asian Precipitation-ENSO relationship during the industrial warming period (1901-2017)? *J. Clim.* **33** 2779–92

Wang B et al 2021 Monsoons climate change assessment *Bull. Am. Meteorol. Soc.* **102** E1–19

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

Wickham H 2016 *ggplot2: Elegant Graphics for Data Analysis* (New York, NY: Springer) 978-3-319-24277-4

Wickham H, François R, Henry L and Müller K 2021 dplyr: A Grammar of Data Manipulation. R package version 1.0.5 https://CRAN.R-project.org/package=dplyr

Wu X, Zhu X, Pan M and Zhang M 2014 Seasonal variability of oxygen and hydrogen stable isotopes in precipitation and cave drip water at Guilin, southwest China *Environmental Earth Sciences* **72** 3183–91

Yang H, Johnson K R, Griffiths M L and Yoshimura K 2016 Interannual controls on oxygen isotope variability in Asian monsoon precipitation and implications for paleoclimate reconstructions *J. Geophys. Res.* **14** 8410–8428

Yim S Y, Wang B, Liu J and Wu Z 2014 A comparison of regional monsoon variability using monsoon indices *Clim. Dyn.* **43** 1423–37

Yoon J-H and Zeng N 2010 An Atlantic influence on Amazon rainfall *Clim. Dyn.* **34** 249–64

Yoon J and Yeh S-W 2010 Influence of the Pacific decadal oscillation on the relationship between El Niño and the northeast Asian summer monsoon *J. Clim.* **23** 4525–37

Zeileis A, Kleiber C, Kraemer W and Hornik K 2003 Testing and dating of structural changes in Ppractice *Comput. Stat. Data Anal.* **44** 109–23

Zeileis A, Leisch F, Hornik K and Christian K 2002 strucchange: An R package for testing for structural change in linear regression models *Journal of Statistical Software* **7** 1–38

Zhang H et al 2021 A data-model comparison pinpoints Holocene spatiotemporal pattern of East Asian summer monsoon *Quat. Sci. Rev.* **261** 106911

Zhang X, Wu M, Liu Y, Hao Z and Zheng J 2018 The relationship between the East Asian Summer Monsoon and El Niño-Southern Oscillation revealed by reconstructions and a control simulation for millennium *Quat. Int.* **493** 106–13

Zhang Q, Wang Y, Singh V P, Gu X, Kong D and Xiao M 2016 Impacts of ENSO and ENSO Modoki + A regimes on seasonal precipitation variations and possible underlying causes in the Huai River basin, China *J. Hydrol.* **533** 308–19

Zhang Z, Leduc G and Sachs J P 2014 El Niño evolution during the Holocene revealed by a biomarker rain gauge in the Galápagos Islands *Earth Planet. Sci. Lett.* **404** 420–34

Zheng X-T, Xie S-P, Du Y, Liu L, Huang G and Liu Q 2013 Indian Ocean dipole response to global warming in the CMIP5 multimodel ensemble *J. Clim.* **26** 6067–80