A data-model approach to interpreting speleothem oxygen isotope records from monsoon regions on orbital timescales

Reconstruction of past changes in monsoon climate from speleothem oxygen isotope (δO) records is complex 10 because δO 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 from version 2 of the Speleothem Isotopes Synthesis and Analysis (SISAL) database to produce composite regional trends in δO in monsoon regions; compositing minimises the influence of site-specific karst and cave processes that can influence individual site records. We compare speleothem δO observations with isotope-enabled climate 15 model simulations to investigate the specific climatic factors causing these regional trends. We focus on differences in δO signals between interglacial (mid-Holocene and Last Interglacial) and glacial (Last Glacial Maximum) states, and on δO evolution through the Holocene. Differences in speleothem δO between the mid-Holocene and Last Interglacial in the East Asian and Indian monsoons are small, despite the larger summer insolation values during the Last Interglacial. Last Glacial Maximum δO values are significantly less negative than interglacial values. Comparison with simulated glacial-interglacial 20 δO shows that changes are principally driven by global shifts in temperature and regional precipitation. Holocene speleothem δO records show distinct and coherent regional trends. Trends are similar to summer insolation in India, China and southwestern South America, but 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 δO and climate variables show that global Holocene 25 monsoon δO trends are driven by changes in precipitation, atmospheric circulation and (to a lesser extent) source area temperature, whilst precipitation recycling is non-significant. However, there are differences in regional scale mechanisms; there are clear relationships between changes in precipitation and in δO for India, southwestern South America and the Indonesian-Australian regions, but not for the East Asian monsoon. Changes in atmospheric circulation contributes to δO trends in the East Asian, Indian and Indonesian-Australian monsoons, and a weak source area temperature effect is observed 30 over southern and central America and Asia. Precipitation recycling is influential in southwestern South America and southern https://doi.org/10.5194/cp-2020-78 Preprint. Discussion started: 13 July 2020 c © Author(s) 2020. CC BY 4.0 License.

precipitation to less negative δ 18 O winter precipitation (Cruz et al., 2005;Dong et al., 2010;Wang et al., 2001). However, the multiplicity of processes that influence δ 18 O before incorporation in the speleothem make it difficult to attribute the climatic 65 causes of changes in individual speleothem records unambiguously. In the East Asian monsoon, for example, speleothem δ 18 O records have been interpreted as regional rainfall signals (Cheng et al., , 2009Wang et al., 2001;Yuan et al., 2004) based on the relationship between modern δ 18 Oprecip and climate. However, Maher (2008) interpreted δ 18 Ospel as reflecting changes in moisture source area, based on differences between δ 18 Ospel 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 70 that the change in precipitation required to reproduce the observed changes in δ 18 Ospel was unreasonably large even if the effects of seasonality and temperature were considered. They therefore argued that changes in moisture source were required to explain shifts in δ 18 O both on glacial/interglacial time scales and during interglacials. There are also conflicting interpretations of the causes of δ 18 Ospel variability in other monsoon regions. In the Indonesian-Australian monsoon region, for example, δ 18 Ospel variability has been interpreted as both a precipitation amount signal (Carolin et al., 2016;Krause et al., 2019) 75 and as a result of changing moisture source regions (Griffiths et al., 2009;Wurtzel et al., 2018). South American speleothem records have been interpreted as records of monsoon intensity, either due to a change in the amount of precipitation, based on relationships derived from instrumental data (Hardy et al., 2003) or due to changes in the ratio of precipitation sourced from the low-level jet versus the Atlantic (Cruz et al., 2005;Wang et al., 2006).
The sources of δ 18 Ospel variability have been explored using isotope-enabled climate models. Modelling studies (e.g. LeGrande 80 and Schmidt, 2009;Lewis et al., 2010;Pausata et al., 2011) suggest that past changes in δ 18 Oprecip do not reflect local rainfall variability in the East Asian monsoon but instead reflect changes in δ 18 O of vapour delivered to a region, consistent with interpretation by Maher (2008). Variability in the δ 18 O of vapour on orbital timescales has been attributed 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;. δ 18 Oprecip variability in the 85 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 generally support the qualitative interpretations of speleothem records from other regions, suggesting that changes in precipitation amount are the predominant source of δ 18 O variability in the Indian monsoon during the Holocene (LeGrande and Schmidt, 2009) and in the glacial (Lewis et al., 2010), and in the South American and Indonesian/Australian regions during Heinrich events (Lewis et al., 2010) and the Last Interglacial (Sjolte and Hoffman, 2014). 90 In this study, we combine speleothem δ 18 O records from version 2 of the Speleothem Synthesis and Analysis (SISAL) database with isotope-enabled palaeoclimate simulations from two climate models to investigate the causes of changes in δ 18 O in monsoon regions through the Holocene (last 11,700 years) and between interglacial (mid-Holocene and Last Interglacial) and glacial (Last Glacial Maximum) states. We compare δ 18 Ospel signals across geographically separated cave sites to extract a regional signal, thus minimising the influence of karst and in-cave processes, such as the mixing of groundwaters from different 95 precipitation events or changes in cave ventilation, that can be important for the δ 18 Ospel of individual records. We use Principal Coordinate Analysis (PCoA) to identify regions with geographically coherent δ 18 Ospel records, and then examine how these https://doi.org/10.5194/cp-2020-78 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
regions behave on glacial-interglacial time scales and through the Holocene. We use isotope-enabled model simulations, where the processes that influence δ 18 Oprecip are explicitly simulated, to investigate the main drivers of regional δ 18 Ospel variability.
We also use multiple regression analysis to identify the likely causes of the observed δ 18 Oprecip changes and trends in specific 100 monsoon regions.

Speleothem oxygen isotope data
Speleothem δ 18 O records were obtained from the SISAL (Speleothem Isotopes Synthesis and Analysis) database (Atsawawaranunt et al., 2018;Comas-Bru et al., 2020a, 2020b 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.

Climate model simulations 120
We use two isotope-enabled climate models: ECHAM5 and GISS E-R. We used published simulations from the ECHAM5 climate model, run either in atmosphere-only  or in fully coupled mode (Gierz et al., 2017b;Werner et al., 2016)  Intercomparison Project) protocols, using reasonable and recommended boundary conditions (trace gases, orbital configuration etc).
The ECHAM5-wiso MH experiment (Wackerbarth et al., 2012;Werner, 2019) was forced by orbital parameters (based on 130 Berger and Loutre, 1991) and greenhouse gas (GHG) concentrations (CO2 = 280 ppm, CH4 = 650 ppb, N2O = 270 ppb) https://doi.org/10.5194/cp-2020-78 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License. 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 1956-1999 , 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) 135 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 δ 18 O 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 . 140 Both the MH and LGM simulations were run at T106 horizontal grid resolution, approximately 1.1° by 1.1°. The LIG experiment (Gierz et al., 2017b(Gierz et al., , 2017a was run using the ECHAM5/MPI-OM Earth System Model, with stable water isotope diagnostics included in the ECHAM5 atmosphere model , 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° by 3.75°. The LIG simulation was forced by orbital parameters derived from the Berger and 145 Loutre (1991) solution and GHG concentrations (CO2 = 276 ppm, CH4 = 640 ppb, N2O = 263 ppb) appropriate to 125 ka, but it was assumed that ice sheet configuration and land-sea geography is unchanged from modern and 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.
There are GISS ModelE-R (LeGrande and Schmidt, 2009) simulations for eight time slices during the Holocene (9 ka, 6 ka, 5 150 ka, 4 ka, 3 ka, 2 ka, 1 ka and 0 ka). The 0 ka experiment is considered as the pre-industrial control (ca 1880). Orbital parameters were based on the Berger and Loutre (1991) solution 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 δ 18 O to account for this (Carlson et al., 2008). The ice sheet in all the other experiments was specified to be 155 the same as modern, and therefore no adjustment was necessary. The simulations were run using the M20 version of GISS ModelE-R, which has a horizontal resolution of 4° by 5°. Each experiment was run for 500 years and we use the last 100 simulated years for the analyses.

Principle Coordinate Analysis and Redundancy Analysis
We used Principle Coordinate Analysis (PCoA) to identify regionally coherent patterns in the speleothem δ 18 O records for the 160 Holocene. PCoA is a multivariate ordination technique that uses a distance/dissimilarity 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 is problematic for many ordination techniques. PCoA is more robust to missing data that other methods (Kärkkäinen and Saarela, 2015;Rohlf, 1972). We used a correlation matrix of speleothem records as the https://doi.org/10.5194/cp-2020-78 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
(dis)similarity measure. The temporal resolution of speleothem records was first standardised by calculating a running average 165 mean with non-overlapping 500-year windows. This procedure produces a single composite record when there are several records for a given site. 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 170 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).

Glacial-interglacial changes in δ 18 O
We examined glacial-interglacial shifts in 18 Ospel and δ 18 Oprecip in regions influence by the monsoon, focusing on regional differences between MH, LGM and LIG with respect to the present-day for speleothems or the control simulation experiment 175 for model outputs. Comas-Bru et al. (2019) have shown that differences in speleothem δ 18 O data between the 20 th century and the pre-industrial period (i.e. 1850±15 CE) are within the temporal and measurement uncertainties of the data, and thus the use of different reference periods (i.e. PI for the ECHAM LIG experiment, 20 th century for ECHAM MH, LGM experiments) should have little effect on our analyses. We used mean site δ 18 Ospel values for each period for the regions identified in the PCoA analysis. Where there are multiple speleothem δ 18 O records for a site in a time period, they were averaged to calculate 180 mean δ 18 Ospel. Three sites above 3500m were excluded from the calculation of the means because high elevation sites have more negative δ 18 O 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 δ 18 Ospel anomalies from the speleothem data. We therefore calculated anomalies with respect to modern (1960-2017 CE) using as reference a global gridded dataset of interpolated mean annual precipitation-weighted 185 δ 18 Oprecip data (Bowen, 2018;Bowen and Revenaugh, 2003). This dataset combines data from 348 stations from the Global OIPC δ 18 Oprecip was converted to its speleothem equivalent assuming that: (i) precipitation-weighted mean annual δ 18 Oprecip is equivalent to mean annual drip-water δ 18 O (Yonge et al., 1985) and (ii) precipitation of calcite is consistent with the empirical 190 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): where δ 18 Ocalcite_SMOW and δ 18 Oaragonite_SMOW are the speleothem isotopic composition for calcite and aragonite speleothems with reference to the V-SMOW standard (in permil); wδ 18 Oprecip is the OIPC precipitation-weighted annual mean isotopic composition of precipitation with respect to the V-SMOW standard and T is the mean annual cave temperature (in degrees Kelvin). We used the long-term (1960-2016) mean annual surface air temperature from the CRU-TS4.01 database (Harris et al., 2014) at each site as a surrogate for mean annual cave air temperature. The resolution of the gridded data means that 200 wδ 18 Oprecip_SMOW and T may be the same for nearby sites.
We use the V-SMOW to V-PDB conversion from Coplen et al. (1983), which is independent of speleothem mineralogy: where δ 18 OPDB is relative to the V-PDB standard and δ 18 OSMOW is relative to V-SMOW standard.
Average uncertainties in the speleothem age models are ~50 years during the Holocene. This interval is smaller than the time 205 windows used in this analysis, and age uncertainty is therefore expected to have a negligible impact on the results. We investigated the influence of age uncertainties on the LGM and LIG δ 18 Ospel anomalies by examining the impact of using different window widths (± 500, ± 700, ± 1000, ± 2000 years) on the regional mean δ 18 Ospel anomalies.
We used anomalies of wδ 18 Oprecip, mean annual surface air temperature (MAT) and mean annual precipitation (MAP) from the ECHAM5-wiso simulations to investigate the drivers of glacial-interglacial δ 18 Ospel variability. Values were calculated from 210 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 x 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 δ 18 Ospel were in the same direction and of a similar magnitude to the simulated wδ 18 Oprecip. 215

Holocene and Last Interglacial regional trends
Speleothem δ 18 O 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 4-step procedure, modified from Marlon et al. (2008): (i) the δ 18 O data for individual speleothems were standardised to zscores by subtracting the mean δ 18 Ospel for a defined base period (7,000 to 3,000 years BP) then dividing these anomalies by 220 the base period standard deviation. The base period was chosen to maximise the number of records included in each composite; (ii) the standardised data for a site were re-sampled 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; (iii) each regional composite was constructed using locally weighted regression (Cleveland and Devlin, 1988) with a window width of 3,000 years and fixed target points in time; and (iv) confidence intervals (5 th and 95 th percentiles) 225 for each composite were generated by bootstrap resampling by site over 1,000 iterations. There are too few sites to construct regional composites for the LIG and thus the trends in δ 18 Ospel were examined using records from individual sites covering the period 130-116 ka BP.
We calculated Holocene regional composites from simulated δ 18 Oprecip from the GISS model. Simulated δ 18 Oprecip trends were calculated using linear distance-weighted mean δ 18 Oprecip values from land grid cells (>50% land) within ±4° around each site. 230 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 δ 18 Oprecip were substantially different if the sampled area was reduced to match the sampling used for the ECHAM-based box plot analysis, or was increased 235 to encompass the larger regions shown in Fig. 1 and used in the multiple regression analysis.

Multiple regression analysis
We investigate the drivers of regional δ 18 Oprecip, and by extension δ 18 Ospel, through the Holocene using multiple linear regression (MLR) of simulated δ 18 Oprecip and climate variables chosen to represent the four potential large-scale drivers of regional changes in the speleothem δ 18 O records. Specifically, we use changes in mean precipitation and precipitation recycling 240 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 δ 18 Ospel composites. Moisture source area limits ( Fig. 1) were defined based on moisture tracking 245 studies (Bin et al., 2013;Breitenbach et al., 2010;D'Abreton and Tyson, 1996;Drumond et al., 2008Drumond et al., , 2010Durá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). 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 250 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 E-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 255 Brubaker et al. (1993): 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 https://doi.org/10.5194/cp-2020-78 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
(QH). RI therefore expresses the change in the contribution of local, recycled precipitation independently of any overall change 260 in precipitation amount.
We use pseudo-R 2 to determine the goodness-of-fit for the overall 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 δ 18 Oprecip when the effects of the other variables are removed.
All statistical analyses were performed in R (R Core Team, 2019) and plots were generated using ggplot (Wickham, 2016). 265

Principle Coordinate Analysis and Redundancy Analysis
The first two PCoA axes are significant, according to the broken stick test, and account for 65% and 20% of δ 18 Ospel variability respectively ( Table 1) The RDA supports a geographical control on the (dis)similarity of speleothem δ 18 O records over the Holocene (Fig. 2b). RDA1 explains 37% of the variability and is significantly correlated with both latitude and longitude (Table 2).

Regional interglacial-glacial differences
Only the ISM, EAM and IAM regions have sufficient data (i.e. at least one record from every time period) to allow comparisons 280 across the MH, LGM and LIG (Fig. 3) and have similar shifts in observed δ 18 Ospel and simulated δ 18 Oprecip. The regional mean δ 18 Ospel anomalies calculated for different time windows (± 500, ± 700, ± 1000, ± 2000 years) vary by less than 0. In both the ISM and the EAM, differences in δ 18 Ospel values between the MH and LIG are small (Fig. 3a,b), although ISM LIG δ 18 Ospel values are slightly more negative than MH values. In the IAM, MH values are less negative than the LIG (Fig. 3c). 290 https://doi.org/10.5194/cp-2020-78 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
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. The glacial-interglacial changes in δ 18 Oprecip are consistent with the simulated temperature and precipitation changes, 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 δ 18 Ospel in the ISM and IAM, since the LIG is wetter than the MH in 295 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 δ 18 O signals in this region.

Regional-scale interglacial δ 18 O evolution
There are four regions with sufficient data to examine Holocene trends: EAM, ISM, SW-SAM and IAM (Fig. 4). The IAM region has the fewest records (n=7) whilst the EAM has the largest number (n=14). The regional composites are z-scores, i.e. 300 anomalies with respect to the base period (3000-7000 yr 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. 4 left) show the most positive δ 18 Ospel z-scores around 12 ka followed by a rapid decrease towards their most negative values at ~9.5 ka and ~9 ka, respectively. The δ 18 Ospel 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 δ 18 Ospel z-scores towards the present in both regions thereafter. The SW-305 SAM records (Fig. 4i) have their most positive δ 18 Ospel 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 12ka, gradually decreasing until ca 5 ka and are relatively flat thereafter.
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 310 initial sharp decrease in δ 18 Ospel values of about 4 ‰ between 130-129 ka and then most of the records (Dykoski et al., 2005;Kathayat et al., 2016;Wang et al., 2008) show little variability for several thousand years. Despite the fact that the Tianmen record (Cai et al., 2010(Cai et al., , 2012 shows considerable variability between 123-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  shows increasingly negative δ 18 Ospel values through the LIG. The trend shown for Whiterock cave 315 (Carolin et al., 2016) also shows 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).

Multiple regression analysis of Holocene δ 18 Oprecip
The global multiple linear regression model that includes all of the Holocene (1 to 9ka) δ 18 Oprecip regional trends has a pseudo-R 2 of 0.80 and shows statistically significant relationships between the anomalies in δ 18 Oprecip and anomalies in precipitation 320 amount, temperature and surface wind direction (Table 3). There is a strong negative relationship with precipitation amount (t value = -8.75) and a strong positive relationship (t value = 8.03) with surface wind direction over the moisture source region, https://doi.org/10.5194/cp-2020-78 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
an index of changes in either source area or moisture pathway. The relationship with temperature over the moisture source region is weaker, but positive (t value = 2.05). Precipitation recycling is not significant in this global analysis.
There are too few data points to make regressions for individual monsoon regions, but the distribution of data points for each 325 region in the partial residual plots (Fig. 6) is indicative of the degree of conformity to the global model. Data points from the ISM, SW-SAM, IAM and SAfM are well aligned with the global relationship with precipitation amount (Fig. 6a), indicating that precipitation amount is an important control on changes in δ 18 Oprecip in these regions. The NE-SAM, EAM and CAM values deviate somewhat from the global relationship and, although there are relatively few points, this suggests that changes in precipitation are a less important influence on δ 18 Oprecip changes in these regions. The impact of temperature changes (Fig.  330 6b) in the ISM, EAM and SW-SAM is broadly consistent with the global model. The slope of the relationship with temperature is negative for the IAM and NE-SAM, and since 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 δ 18 Oprecip 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 335 with the global relationship between changes in δ 18 Oprecip and changes in wind direction (Fig. 6c), indicating that changes in source area or moisture pathway are an important control on changes in δ 18 Oprecip in these regions. However, values for CAM, SW-SAM, NE-SAM and SAfM deviate strongly from the global relationship. Recycling does not appear to be an important contributor to changes in δ 18 Oprecip except in SW-SAM and SAfM (Fig. 6d).

Discussion 340
We have shown that it is possible to derive an objective regionalisation of speleothem records based on PCoA of the oxygenisotope 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 345 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 350 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). https://doi.org/10.5194/cp-2020-78 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
Palaeoenvironmental records from this region show mixed signals for the Holocene: some have been interpreted as showing 355 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 southern 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 360 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 δ 18 Ospel 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 365 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 3,000 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 370 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, 375 2012). We have 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). Given that the increase in summer insolation is much larger during the LIG than the MH, this finding is again consistent with the idea that other factors play a role in modulating the monsoon response to insolation forcing.
Global relationships between δ 18 Oprecip and its drivers (precipitation amount, temperature and surface wind direction; Fig. 6) 380 are consistent with existing studies: a strong relationship with precipitation and a weaker temperature effect has been widely observed at tropical and sub-tropical latitudes in modern observations (Dansgaard, 1964;Rozanski et al., 1993). The significant global relationship between δ 18 Oprecip and surface winds supports the idea that changes in moisture source and pathway are also important for explaining δ 18 O 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 that 385 suggest changes in precipitation amount (Cai et al., 2015;Fleitmann et al., 2004) and to a lesser extent moisture pathway (Breitenbach et al., 2010) drive δ 18 Ospel variability. The δ 18 O 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 precipitation amount do not seem to explain the observed changes in δ 18 Ospel in the EAM during the Holocene, where Holocene δ 18 Oprecip evolution is largely driven by changes in atmospheric circulation (indexed by changes in 390 surface winds). This is consistent with existing studies that emphasise changes in moisture source and/or pathway rather than a precipitation amount (Maher, 2016;Maher and Thompson, 2012;Tan, 2014;Yang et al., 2014). Speleothem δ 18 O 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 appears to be important, perhaps unsurprisingly given that recycling presently contributes ca 25-35% of the precipitation over the Amazon 395 (Brubaker et al., 1993;Eltahir and Bras, 1994).
The LGM is characterised by a similar orbital configuration to today, however global ice volume was at a maximum and GHG concentrations were lower than present. The δ 18 Ospel anomalies are more positive during the LGM than the MH or LIG, suggesting drier conditions in the ISM, EAM and IAM, supported by simulated changes in δ 18 Oprecip 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 400 source areas, together with a ca 5% decrease in relative humidity (Yue et al., 2011) would result in a water vapour δ 18 O signal at the source 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 δ 18 O values during the LGM (ca. +1 ‰ relative to the MH or LIG; Waelbroeck et al., 2002). Enriched δ 18 Oprecip and δ 18 Ospel values during the LGM must therefore be caused by a significant decrease in 405 atmospheric moisture and precipitation.
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. There are only 14 records from monsoon regions covering the LIG in SISALv2, for example. Nevertheless, our analyses show that there 410 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 415 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 δ 18 O values. Thus, the interglacialglacial contrast in regional δ 18 Ospel is also robust to age uncertainties. 420 Isotope-enabled climate models are used in this study to explain observed regional-scale trends in δ 18 Ospel. 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 δ 18 Oprecip (of 0.5 ‰) at the MH. However, both models show trends in δ 18 Oprecip that reproduce the 425 observed changes in regional δ 18 Ospel (Figs 3 and 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 δ 18 Ospel 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 isotopeenabled simulations (Caley et al., 2014;Risi et al., 2010).
This study illustrates a novel data-model approach to investigate the drivers of δ 18 Ospel under past conditions, by comparing 430 composite regional records and then using multiple linear regression of isotope-enabled palaeoclimate simulations to determine the role of individual drivers of these trends. This obviates the need to use modern δ 18 Oprecip-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. 435

Conclusions
Geographically distributed speleothem δ 18 O records and isotope-enabled climate models can be used together to understand the drivers of δ 18 Ospel on orbital timescales. 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 Glacial δ 18 Ospel signals are best explained by a large decrease in 440 precipitation, as a consequence of lower atmospheric moisture content driven by global cooling. The evolution of δ 18 Ospel through the Holocene across the global monsoon domain is driven by 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 δ 18 O evolution through the Holocene relates to changes in atmospheric circulation (i.e. changes in moisture pathway and/or source). Changes in precipitation amount are the predominant driver of Holocene δ 18 Ospel evolution in the 445 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 in southwestern South America.

Code and data availability 450
The SISAL (Speleothem Isotopes Synthesis and AnaLysis) database version 2 is available through the University of Reading

Author contributions
The study was designed by SP, SPH, LCB and NK. MW and ALG provided climate model outputs. SP ran the analyses. The first draft of the manuscript was written by SP, SPH and LCB and all authors contributed to the final draft. 460

Competing interests
The authors declare no competing interests.      (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 840 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.