Cave setting

The two stalagmites used in this study (LR06-B1 and LR06-B3) were collected from Liang Luar, an ∼1.7-km-long cave situated on the east Indonesian island of Flores (8° 32’N, 120° 26’E; 550 m above sea level; Fig. 1 and Supplementary Figs 1 and 2). Both specimens were located in a large chamber ∼600 m from the cave entrance and were actively growing at the time of collection in August 2006. Liang Luar resides close to the southerly limit of the austral summer ITCZ and most of the annual rainfall (∼70%) occurs during the AISM season (December to March)11.

Figure 1: ENSO influence on tropical hydroclimate. Field correlation map of annual rainfall and east equatorial Pacific sea surface temperatures. Correlation coefficients (r) were calculated between GPCP v2.2 precipitation (at each 1° × 1° grid point) and HadISST1 NINO3.4 SSTs (1979–2010). The coloured dots show the location of our study site and other paleoclimate records mentioned in the text. Yellow indicates relatively dry conditions during the Medieval Climate Anomaly5,6,9,44,46,48 (950–1250 CE), whereas cyan/purple indicates relatively wet conditions7,8/intense precipitation events10,47. Full size image

Stable isotope and trace element records

The hydroclimate record derived from stalagmites LR06-B1 and LR06-B3 is fixed in time by 35 uranium-thorium (U-Th) ages (Supplementary Figs 4–6) and spans the period 3 to 2005 CE (Common Era, Supplementary Table 1). All ages lie in stratigraphic order (within 2σ uncertainty) and have typical uncertainties of 1–3%. The age models for both records were first calculated using a Bayesian–Monte Carlo approach12 and then refined using the intra-site correlation age modelling programme (iscam)13 written in MATLAB (see the Methods for details). The age models give an average growth rate of ∼125–135 μm per year for both stalagmites.

The δ18O and δ13C records are based on 1,157 measurements performed on calcite powders micromilled along the central growth axes of LR06-B1 and LR06-B3. The Mg/Ca and Sr/Ca ratios in stalagmite LR06-B1, previously presented in Griffiths et al.14, were analysed on H 3 PO 4 residues following calcite dissolution for stable isotope analysis. The strong degree of replication between the two stable isotope profiles (Fig. 2, Supplementary Fig. 7 and Supplementary Table 2), along with ‘Hendy’ screening tests15 conducted on these specimens11, rules out significant disruption of the records by kinetic isotope fractionation effects.

Figure 2: Multi-proxy reconstruction of east Indonesian hydroclimate. Comparison of stalagmites LR06-B3 and LR06-B1 (a) δ18O, (b) δ13C, (c) Mg/Ca (cyan) and Sr/Ca (red), (d) LL PC1 (black line) and [234U/238U] i (purple dots: LR06-B3; blue dots: LR06-B1; dashed line: cubic spline fitted to the data) and (e) LL PC2 . We interpret higher (lower) values in LL PC1 and [234U/238U] i in d to reflect drier (wetter) conditions. Mg/Ca and Sr/Ca records were previously reported in Griffiths et al.14. The PC factor loadings for each time series of PC1 and PC2 are shown in parentheses. PC1 accounts for 50% of the variance and PC2 accounts for 20%. Colour-coded symbols along the bottom indicate U-Th ages (purple: LR06-B3; blue: LR06-B1) and their associated 2σ uncertainties. The grey shading in d,e indicates the 95% bootstrap confidence interval. Vertical bars indicate the approximate timing of the Medieval Climate Anomaly (MCA; yellow), Little Ice Age (LIA; blue) and Current Warm Period (CWP; pink) in Flores, defined by decreased AISM rainfall during the MCA and CWP, and increased AISM rainfall during the LIA. Full size image

The δ18O values for LR06-B1 and LR06-B3 range from −5.6 to −6.9‰ and average −6.2‰ between 1 and 2005 CE, whereas the δ13C values range from −7.8 to −13.4‰ and average −10.8‰. Broadly speaking, the δ18O and δ13C records display multidecadal- to centennial-scale oscillations with generally lower values from ∼500 to 900 CE and ∼1300 to 1600 CE, and generally higher values from ∼1000 to 1300 CE The Mg/Ca and Sr/Ca profiles show patterns similar to the stable isotopes over the last millennium, with generally higher values from ∼1000 to 1300 CE followed by a trend to lower values until ∼1600 CE Between ∼1600 CE and ∼2004 values have generally increased, signifying an overall reduction in rainfall over this period with the trend being most significant during the twentieth century.

Interpretation of the stalagmite δ18O record

The seasonal pattern of monsoon rainfall at Liang Luar is reflected in the δ18O of precipitation (δ18O p ), with summer rainfall (DJFM) ∼6–7‰ lighter than its winter monsoon (JJA) counterpart11. Hence, a change in the fraction of the year dominated by the summer monsoon (that is, change in seasonality), or a shift in monsoon rainfall amount in summer, should have a significant impact on the average annual δ18O value. The strong relationship between Indonesian rainfall amount and δ18O p has been identified in modern isotope studies11 and CGCM simulations of modern16,17 and past18 hydroclimates in Indonesia. Model simulations of precipitation δ18O p using the National Center for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) nudged IsoGSM17 for the grid point closest to Liang Luar reveal a significant correlation with sea-surface temperatures (SSTs) over the NINO3.4 region and, to a lesser extent, in the western Indian Ocean; this pattern is similar to the observed pattern of rainfall amount at Liang Luar (Supplementary Fig. 3). Moreover, comparison of Liang Luar IsoGSM δ18O p with observed rainfall amount, outgoing longwave radiation and zonal 850 mb winds reveals that δ18O p is strongly influenced by western Pacific rainfall variability, which is governed by the strength of the AISM, the position of the ITCZ and ENSO (Supplementary Fig. 3).

Based on observed and modelled precipitation δ18O systematics, we interpret variations in the Flores speleothem δ18O records to primarily reflect regional rainfall amount. However, on relatively short timescales (centuries), other factors may affect stalagmite δ18O, such as changes in cave temperature (due to shifts in regional SSTs), the δ18O of surface seawater5, the proportion of summer- and winter-sourced moisture11 and multidecadal changes in the relationship between δ18O p and ocean-atmosphere climate modes, such as ENSO19 and the Indian Ocean Zonal Mode20. In addition, mixing of meteoric water in the soil-karst overburden (∼50–100 m) above Liang Luar reduces the amplitude of the highest frequencies21 (for example, seasonal-annual variations), which likely varies with time. Therefore, some of the observed variability may be an artefact of filtering through the karst network. However, given the amount of cave recharge through the epikarst each monsoon season, it is likely that the δ18O of the cave drip-water represents decadal (if not multi-year) shifts in precipitation.

A composite multi-proxy stalagmite record

To constrain the Liang Luar paleomonsoon record, we compare the stalagmite δ18O with the contemporaneous carbon isotope (δ13C), Mg/Ca and Sr/Ca ratios, which primarily reflect Liang Luar karst hydrology14,22. Recent geochemical modelling of the soil-karst system above Liang Luar demonstrated that the δ13C in the late Holocene section of stalagmite LR06-B1 was controlled by the degree of prior calcite precipitation (PCP)22 in the epikarst ‘upstream’ of the stalagmite. Under this process, periods of reduced recharge result in partial dewatering of the karst fracture system overlying the cave, which causes a greater proportion of the CO 2 degassing to occur before the drips reach the growing stalagmite23. The waters therefore become supersaturated in CaCO 3 before their emergence in the cave, causing upstream (prior) calcite precipitation in fractures and on stalactite tips. In the case of δ13C, carbonate precipitation results in a higher stalagmite δ13C owing to the greater removal of 12CO 2 during the degassing process24,25. The opposite occurs during wetter intervals25 (Fig. 2). Additional climate-related factors may have influenced the stalagmite δ13C, such as changes in the ratio of C3:C4 plants, soil respiration rates, contributions from atmospheric CO 2 (refs 26, 27) and the relative contribution of bedrock carbon due to changes in open- versus closed-system dissolution22. However, these factors are all susceptible to hydrological changes and fluctuate in the same sense as PCP (that is, higher δ13C values indicate drier conditions). Therefore, regardless of their relative contributions to the carbon-isotope mass balance, shifts in stalagmite δ13C are inextricably linked to hydroclimate variations at this site.

Additional support for a PCP signal is provided by the pattern of Sr/Ca and Mg/Ca changes over the last 2,000 years (Supplementary Fig. 8). Covariation between these two ratios is diagnostic of a hydrologic control23,28: increases in both indicate drier conditions because the preferential loss of Ca2+ from solution (relative to both Mg2+ and Sr2+) in the percolation waters during PCP raises the Mg/Ca and Sr/Ca of the drip waters. This is transmitted to the stalagmite as higher Mg/Ca and Sr/Ca. During wetter intervals, PCP is either reduced or absent, leading to lower Mg/Ca and Sr/Ca in the stalagmite23. The significant correlation between δ13C and both trace element ratios provides even more compelling evidence of hydroclimate-driven changes in the degree of PCP over the last 2,000 years (Supplementary Fig. 8). It not only reinforces the similar findings recorded in Liang Luar for the entire Holocene14, but is in good agreement with previous studies of short-term (seasonal to annual)26 and long-term (century to orbital)29 hydroclimate change.

To effectively combine the climate proxies for stalagmites LR06-B1 and LR06-B3 into one composite hydroclimate record, we conducted an unrotated principal component (PC) analysis on the six individual stable isotope and trace element time series. The leading PC (from here referred to as ‘LL PC1 ’) represents the dominant hydrologic signal embedded in the records30 (Fig. 2). The stalagmite δ18O does not significantly load on LL PC1 indicating that the processes important in altering precipitation δ18O (vapour transport, oceanic source and regional temperatures)31 may not strongly affect the amount of moisture at the cave site. These additional effects on δ18O p are best captured in the second PC (Fig. 2). Therefore, at century timescales, δ18O p probably reflects the over-riding influence of large-scale ocean-atmospheric processes that do not strongly affect the total amount of precipitation at the site32. These findings highlight the importance of utilizing multiple speleothem geochemical tracers to reconstruct tropical hydroclimate during the relatively stable background climate state of the last two millennia.

Additional evidence to support LL PC1 as a hydrologic proxy is provided by initial uranium isotope activity ratios ([234U/238U] i ) in the Liang Luar stalagmites (Supplementary Table 1 and Fig. 2d). Holocene [234U/238U] i ratios of LR06-B1 and LR06-B3 reflect the extent of dissolution related to the rate of groundwater movement in the epikarst of Liang Luar14,30. The [234U/238U] i values increase during drier periods when relatively slow-moving percolation waters (in approximately chemical equilibrium with the host rock) preferentially remove 234U (produced by alpha recoil or preferential leaching and loosely bound in limestone bedrock)33. In contrast, during wetter periods, the more rapid traverse of the percolation water through the karst network can result in greater dissolution rates causing both U isotopes to be leached approximately equally from the host rock. This can lead to relatively low [234U/238U] i values in a stalagmite. Comparison of the [234U/238U] i values with LL PC1 reveals similar multicentury trends over the past 2,000 years.

A final confirmation that LL PC1 reflects a regional hydroclimate signal over the western Pacific is provided by the significant correlation with other paleohydrologic records from the region, including a seawater δ18O record5 (r=0.5, 50-year smoothed, significant at the 99% level) and a leaf wax δD record6 (r=0.28, 50-year smoothed, significant at the 95% level) from marine sediment cores drilled in the Makassar Strait immediately west of Sulawesi (Fig. 1).

Meridional shifts in the Australasian ITCZ

The LL PC1 record exhibits pronounced century-scale excursions, with positive values (drier climate) at ∼0–400 CE, ∼1000–1400 CE and ∼1900–2004 CE, and negative excursions (wetter climate) at ∼400–1000 CE and ∼1400–1900 CE (Fig. 3). The dry period at around the turn of the last millennium coincides with the approximate timing of the Medieval Climate Anomaly (MCA; ∼950–1250 CE)34. This ∼400-year reduction in AISM rainfall begins to recover in strength at around 1300, and by 1500 is above average during the Little Ice Age (LIA). Maximum rainfall in Flores at ∼1600 CE, among the wettest periods of the past 2,000 years, is synchronous (within dating uncertainty) with peak cooling in the Northern Hemisphere34 and maximum ice discharge in the North Atlantic35 (Fig. 3a).

Figure 3: Tropical and high northern latitude teleconnections. (a) North Atlantic drift ice index35. (b) Stalagmite δ18O record from Wanxiang Cave, China7. (c) Mean grain size (MGS) of lake core DY6 from Dongdao Island in the South China Sea (SCS)9. (d) Northwest Australian flood events as recorded in a stalagmite from KNI-51 cave42. (e) Total solar irradiance reconstructed from cosmogenic nuclides43. The black line and grey shading in all panels shows the LL PC1 record and corresponding 95% bootstrap confidence interval. The y axis for the LL PC1 record is oriented to show a weaker AISM with warmer North Atlantic temperatures. Positions of the U-Th dates (with 2σ uncertainties) for both stalagmites are shown along the bottom. Vertical bars indicate the approximate timing of the MCA (yellow), LIA (blue) and CWP (pink) in Flores, defined by decreased AISM rainfall during the MCA and CWP, and increased AISM rainfall during the LIA. Full size image

Comparison of the Flores LL PC1 record with subtropical paleomonsoon records from China7 (Fig. 3b) reveals synchronous (antiphased) multicentury fluctuations in monsoon rainfall across the hemispheres. The cross-equatorial antiphase relationship is consistent with both CGCMs36 and glacial paleoclimate records37,38, showing that cooling of the high northern latitudes increases the interhemispheric thermal gradient, displacing the Australasian ITCZ southward. However, zonal changes in the PWC are also known to strongly influence western Pacific hydroclimate39, and are thought to be a significant driver of low-latitude climate change on century scales40. Indeed, an in-phase relationship between LL PC1 and paleohydrologic records from the South China Sea9 (Fig. 3c) challenges the paradigm of regional north–south ITCZ migrations as the primary control of lower-frequency monsoon shifts in the western Pacific over the past millennium. Critically, a recent synthesis of paleohydrologic records for the Australasian monsoon region41, including a new northwest Australian flood record42 (Fig. 3d), demonstrated that, rather than moving southward during the LIA, the latitudinal range of monsoon-ITCZ migration probably contracted equatorward, perhaps in response to lower solar irradiance43 (Fig. 3e).

Zonal teleconnections via the PWC

Evidence for a potential link between the PWC and multi-century changes in western Pacific rainfall is provided by the inverse relationship between the AISM5,6 and rainfall in the central8-eastern10,44,45,46,47,48 equatorial Pacific (Fig. 4) over the past millennium. Specifically, the western Pacific and Peruvian Andes (Fig. 4a–e) were generally wetter during the LIA, whereas the central and eastern equatorial Pacific (EEP) experienced drier conditions (Fig. 4f) or reduced heavy precipitation events (Fig. 4g,h), similar to La Niña events today when the PWC strengthens (Fig. 1). Conversely, drier conditions in Indonesia and Peru (along with Panama48) during the MCA were matched by wetter conditions in the central and EEP, signifying a more ‘El Niño-like’ mean state. Wavelet-transform analysis of the Flores LL PC1 shows significant multidecadal to centennial periodicities in the AISM record (Supplementary Fig. 8) that have also been identified in reconstructions of ENSO amplitude from North American tree rings49 and NINO3.4 SSTs40 (Supplementary Fig. 9).

Figure 4: Hydroclimate records for the tropical western and eastern Pacific. (a) Flores LL PC1 record. (b) Marine foraminifera δ18O sw (ref. 5) and (c) terrestrial δD leaf-wax (ref. 6) records recovered from marine sediment cores located in the Makassar Strait on the Sulawesi margin. (d) δ18O of lake sediment calcite in Laguna Pumacocha in the central Peruvian Andes (proxy for the strength of the South American summer monsoon)44. (e) Speleothem δ18O record from Cascayunga cave in northeast Peru46. (f) δD leaf-wax record from Washington Island in the central equatorial Pacific8. (g) Red-colour intensity from Laguna Pallcacocha, southern Ecuador47. (h) Percent sand in El Junco lake, Galápagos Islands10. For clarity, all records have been converted to standard (z) scores with blue indicating wetter conditions (a–f) or heavier precipitation events (g–h) and vice versa for red. Vertical bars indicate the approximate timing of the MCA (yellow), LIA (blue) and CWP (pink) in Flores. Full size image

The multi-proxy record of eastern Indonesian hydroclimate presented here (LL PC1 ) lends support to a LIA characterized by a stronger PWC, whereas the LL PC1 shows drier conditions during the MCA, indicating a weaker PWC. These findings are compatible with the ITCZ expansion (LIA)/contraction (MCA) hypothesis of Yan et al.41 The temporal changes in zonal precipitation across the equatorial Pacific5,6,8,9,10,40 are also consistent with some terrestrial reconstructions of ENSO40,45,50 activity (Fig. 5), which show strengthening of the PWC during the LIA, and weakening during the MCA. These findings run counter to inferences from some SST-based ENSO reconstructions from the western5,51, central52 and EEP53, which suggest that the tropical Pacific was characterized by ‘La Niña-like’ conditions during the MCA compared with the LIA. However, there is still no consensus on the marine paleo-ENSO signals; a recent reconstruction of EEP SSTs54 indicates that the MCA was marked by relatively warm SSTs in the EEP, contradicting the findings of earlier studies.

Figure 5: Western Pacific rainfall and ocean-atmosphere climate modes from paleoclimate records and climate model simulations. Panels a–d show paleoclimate records34,40,50,59,60 and e–h show CCSM4 LM+hist CGCM simulations. Comparison of Flores precipitation (a,e), SOI (b,f), NINO3.4 SSTs (c,g) and land temperatures (red: NH; blue: global; d,h). The records have been smoothed with a 50-year (thick coloured lines) and 7-year (thin grey lines) loess filter. All time series were converted to standard (z) scores to facilitate comparison of proxy reconstructions with model simulations. Vertical bars indicate the approximate timing of the MCA (yellow), LIA (blue) and CWP (pink) in Flores. Full size image

Theoretical and physical models of the PWC are also not entirely consistent within the context of Indo-Pacific climate change during the past millennium. Modern dynamical studies suggest that under warmer conditions, the PWC should weaken55,56, and at least one paleoclimate study supports the inference of a stronger PWC during the LIA41. On the other hand, the ‘ocean dynamical thermostat’ mechanism57 supports a potential weakening of the PWC during the LIA. Modern observations do not yet distinguish the PWC’s response to global warming. Yan et al.41 recently noted that there are numerous interacting components that may contribute to the discrepancies among paleoclimate records: (i) misinterpretation of SST or hydrological signals with respect to meridional versus zonal changes in ocean-atmosphere circulation; (ii) decoupling of ocean and atmospheric processes in response to different forcings over century time scales and (iii) the interacting influence of multidecadal climate modes (for example, Atlantic Multidecadal Oscillation) with ENSO, which are difficult to identify in the proxy records that originate from outside these system’ centres of action.

Palaeoclimate proxy-model comparisons

To explore the potential role of multidecadal- to century-scale changes in the PWC on global climate variability, we examined the relationship between Flores rainfall, NINO3.4 SSTs and the Southern Oscillation Index (SOI) with similar indices in CGCM simulations (Fig. 5). To do so, we employed the Paleoclimate Modelling Intercomparison Project Phase 3 (PMIP3) last millennium (LM; 850–1849 CE) and historical (1850–2005 CE) CGCM runs from the US NCAR CCSM4 (ref. 58) model. The CCSM4 model was employed because its simulation of the modern hydroclimate dynamics of east Indonesia compares well with observations. Importantly, model preindustrial control runs (500-year-long) accurately simulate the influence of ENSO on rainfall variability around Flores (Supplementary Fig. 10 and Supplementary Table 3). Although the association of SST with rainfall is somewhat stronger and more widespread in models compared with observations, the regional pattern of SST is consistent. The dominant influence of ENSO on east Indonesian rainfall is also apparent for the LM and historical (LM+hist; 850–2005 CE) simulations (Supplementary Table 3), although the correlation is not as strong as for the preindustrial control runs.

If variations in global temperatures over the LM were amplified by decadal to multidecadal variations in the tropical Pacific, as is the case during the instrumental period2,3,4, then model skill in the tropical Indo-Pacific region is clearly important for projecting future climate change. Therefore, we compared Flores precipitation, the SOI and NINO3.4 SSTs captured in the proxy records40,50 (Fig. 5a–c) with the same indices in the LM+hist CCSM4 simulations (Fig. 5e–g). In the LM+hist CGCM experiments, the variance in Pacific SST and Indonesian rainfall is concentrated at interannual/decadal time scales, whereas the paleoclimate record for Indonesia exhibits strongest variance at multidecadal to centennial periodicities (Supplementary Fig. 8). However, the coherent patterns of lower-frequency variability evident in the network of western Pacific hydroclimate records suggest that the data-model mismatch likely stems from model deficiencies with respect to century-scale variations in the PWC (Fig. 5b,c,f,g), particularly during the MCA and LIA. Consequently, the models may be underestimating the magnitude and duration of reconstructed global59 and hemispheric34,60 Medieval warming, and the persistent LIA cooling between ∼1400 and ∼1850 CE (Fig. 5d,h).

Our results highlight significant discrepancies between the proxy records and model simulations for the past millennium. Critically, these discrepancies coincide with century-scale anomalies in the strength of the PWC. We cannot rule out the possibility that some of the low-frequency Pacific variability was a forced response to variable solar intensity36 and changing teleconnections to higher latitudes61 that are not simulated by the models, or that non-climatic processes have influenced the proxies62. However, of particular importance is that the paleodata-model mismatch supports the possibility that unforced, low-frequency internal climate variability (that is difficult for models to simulate) was responsible for at least some of the global temperature change of the past millennium62,63.