Regional climatology

The Sahiya Cave ((30°36′N, 77°52′E, ~1,190 m above sea level (m.a.s.l.)) is located ~200 km NNE of New Delhi in the Uttarakhand State of India (Fig. 1; Supplementary Fig. 1). Mean annual temperature and annual precipitation in the region are ~22.0 °C and 1,600 mm, respectively, with ~80% of precipitation falling during the monsoon season (between June and September) (Supplementary Fig. 2). The onset of ISM rainfall in the study area is preceded by the development of the cross-equatorial low-level jet (LLJ) that serves as the primary conduit for transporting moisture across the Indian Ocean and the Arabian Sea onto the Indian sub-continent19. The LLJ generates strong cyclonic vorticity in the planetary boundary layer (typically north of 15°N) leading to large convective instabilities over the central portion of India, a region referred to as the Monsoon Trough (MT)20. At the height of the ISM season (typically during July and August), the LLJ develops into a southeasterly flow, transporting moisture from the Bay of Bengal (BoB) to northwest and northern India21. In addition, the negative pressure anomalies in the MT region trigger the development and propagation of synoptic-scale weather systems, collectively referred to as the low-pressure systems (LPS)20, from the BoB (typically north of 18°N during the months of July and August) onto the Indian landmass. The LPS are instrumental in advancing the ISM as far north as the Lesser Himalayas (that is, north of our study area) and as far west as Pakistan during periods of strong monsoon circulation21. Another important source of precipitation variability in the study area stems from chaotic oscillations in monsoonal precipitation, referred to as active-break spells20,21,22. Distinct from the LPS, these oscillatory modes of monsoon exhibit a hierarchy of quasi-periods (3–7 days, 10–20 days and 30–60 days) and are manifestations of north–south fluctuations in the mean position of the MT20,21,22. During the ‘break’ period, the axis of the MT shifts north of our site, resulting in dry air intrusions of westerly winds into northwest and central-north India, which tend to suppress rainfall through decrease of convective instability and depletion of parcel buoyancy23,24. In summary, the sub-seasonal to interannual precipitation variability in our study area, strategically located at the distal end of the BoB branch of the ISM and at the northern periphery of the MT region, is sensitive to changes in the strength and position of the LLJ and attendant changes in the MT, which in turn, influence the strength of the BoB branch of ISM.

Figure 1: Cave location and ensembles of backward air-parcel trajectories originating from the study area. The HYSPLIT25 trajectory ensembles depict air-parcel/moisture transport routes for JJAS, highlighting the maximum contrast between strong and weak monsoon circulation. The green and black lines in each panel are multi-year monthly composites of daily trajectories originating from the study area (red dot) associated with periods of anomalously high and low rainfall, respectively, during June (a), July (b), August (c) and September (d). Anomalously low rainfall (relative to the 20th century mean) during June, July, August and September in our study area is characterized by near-cessation of 18O-depleted moisture flux from the Bay of Bengal, resulting in relatively higher flux of 18O-enriched moisture from the local sources and Arabian Sea, and dry air intrusions of westerly winds that tend to suppress rainfall over central and northern India through increase of convective stability. On other hand, anomalously high rainfall, particularly during July and August, in our study area is associated with significantly enhanced transport of 18O-depleted moisture from the Bay of Bengal (see Methods and Supplementary Fig. 6). Full size image

Oxygen isotope in modern precipitation

We use the observational data from the Global Network of Isotopes in Precipitation program of the International Atomic Energy Agency (GNIP-IAEA, hereafter GNIP) to characterize the δ18O of precipitation (δ18O p ) in our study area. The GNIP data nearest to the cave site (Uttarkashi, 30.72N, 78.44E; 1,140 m.a.s.l.), however, contains only 3 years of monthly δ18O p observations (AD 2004–2006; n=28) (Supplementary Fig. 3). Consequently, we use the δ18O p data from the GNIP station at New Delhi (28.58N, 77.52E; 212 m.a.s.l., and ~200 km south of Sahiya) to characterize the δ18O p variability in our study area (Supplementary Figs 4–5). Our approach to utilize δ18O p data from New Delhi is justified by the gross similarities between the annual cycles of precipitation and δ18O p at Sahiya and New Delhi (see Methods), which suggest that both sites are influenced by similar precipitation dynamics and moisture sources throughout the year. Furthermore, to examine dynamical linkages between δ18O p and changes in the low-level (~850 hPa) monsoon circulation pattern, we use the hybrid single-particle lagrangian integrated trajectory model (HYSPLIT)25 to generate backward (168 h) moisture trajectory ensembles for those periods where GNIP data are available. Our paired analyses of GNIP δ18O p data and the low-level monsoon moisture trajectory patterns (Fig. 1; Supplementary Fig. 6) suggest that the δ18O p variability in the study area is dominated by large-scale changes in low-level monsoon circulation and moisture transport history upstream of our site, and to a lesser extent from changes in local precipitation intensity. We find that periods of strong (weak) monsoon circulation are characterized by an enhanced (reduced) flux of isotopically depleted BoB moisture and reduced (enhanced) flux of isotopically enriched Arabian Sea moisture to our site (see Methods).

Speleothem oxygen isotope interpretation

The contrasting patterns of strong and weak monsoon circulation produce distinct δ18O p , signatures, which in theory, should predominantly drive changes in δ18O of speleothems from the study area. This inference rests in part on the assumption that non-monsoon rainfall and snowmelt are insignificant sources of surface waters that infiltrate into the cave system. Indeed, this assumption is supported by the δ18O values of the cave’s dripwaters (−6.95±0.31‰, n=7), which are indistinguishable from the amount-weighted δ18O of annual precipitation (−6.97‰) and slightly higher than the δ18O of June, July, August and September (JJAS)-integrated precipitation (−8.19‰) in the study region (Supplementary Fig. 3d). Furthermore, the regional precipitation–evaporation budget is positive only during the summer monsoon season (from June to September) (Supplementary Fig. 2b), implying that effective infiltration of surface waters into the local karst occurs primarily during the ISM season (Supplementary Fig. 3). The karstic overburden over the Sahiya Cave (~30 m) would, however, tend to act as a lowpass filter by mixing δ18O p signals from several seasons, thus damping high-frequency δ18O p variations. Nonetheless, the low-frequency (multi-year to decadal) component of the speleothem δ18O record should reflect time-integrated shifts in δ18O of karstic waters stemming from an asymmetric occurrence of 18O-enriched (depleted) precipitation events during weak (strong) monsoon seasons. Consequently, we interpret temporal changes in the speleothem δ18O at our site to reflect space–time-integrated upstream changes in ISM intensity. Indeed, comparisons of the speleothem δ18O time series with instrumental records of precipitation at local, sub-regional and sub-continental spatial scales (Supplementary Fig. 7) indicate that the strength of the inverse δ18O p –rainfall amount relationship is higher at sub-regional and sub-continental scales. This observation suggests that upstream changes in ISM circulation and rainfall variability rather than changes in the local precipitation amount are the primary drivers of δ18O variability at our site. Our observations are consistent with the results from other studies, which suggest that local δ18O p variability in some monsoonal and tropical regions reflect upstream changes in regional circulation and rainfall26,27.

Speleothem oxygen isotope record

The NI δ18O record from the Sahiya Cave is a composite of two δ18O stalagmite records (SAH-A and SAH-B) consisting of 1848 δ18O measurements with an average temporal resolution of 1.5 years per sample (Supplementary Fig. 8). The chronologic frameworks of SAH-A and SAH-B δ18O profiles are provided by 13 and 15 230Th dates, respectively (Supplementary Fig. 9; Supplementary Data 1). The average 2σ age uncertainty of the composite record is ~22.5 years. The age models are developed by generating multiple cubic interpolations28 between the discrete 230Th dates. The SAH δ18O record (mean=−9.0‰; 1σ=0.24) covers the last two thousand years (BC 141 to AD 2006) (Supplementary Data 2). The δ18O profile is characterized by prominent decadal to multi-decadal variability with an amplitude of ~0.5‰ around the long-term mean (Fig. 2b). The observed weakening trend in ISM rainfall since the 1950s manifests as an ~0.5‰ increase in speleothem δ18O, which is a robust feature of our record even with generous assumptions of age model and analytical uncertainties. Comparison between the composite δ18O profile and a 11-year running mean of AISMR exhibits a significant inverse correlation (r=−0.75 with 95% confidence interval (−0.82; −0.64)) that explains ~50% of variance in the instrument record (Fig. 2a; Supplementary Fig. 7). After correcting for standard instrumental error in δ18O measurements (0.06‰), the slope of the speleothem δ18O-AISMR regression yields a conversion of 1.0‰ in speleothem δ18O values to ~86 mm change in AISMR, the latter equivalent to ~10% of the long-term mean (~850 mm) (Fig. 2a). The δ18O record is characterized by several episodic intervals of anomalously (>2σ) positive and negative δ18O values. Particularly noteworthy are two sub-decadal length intervals of strong pluvial and drier conditions, centred at AD ~1,100 and 1,620, respectively (Fig. 2b). The timing and duration of the latter event, which constitutes the largest positive δ18O anomalies in our reconstruction coincides (within age uncertainties) with widespread droughts and famines in north India29 and with the collapse of the Guge Kingdom in western Tibet30. Incidentally, these two intervals of substantially enhanced and reduced ISM rainfall occur within the periods of peak thermal anomalies during the Medieval Climate Anomaly (nominally, AD 1000–1200) and the Little Ice Age (nominally AD ~1600–1850) (Fig. 2), respectively.

Figure 2: Comparison of the NI speleothem δ 18 O record with instrumental and proxy data. (a) The comparison between the NI speleothem time series (raw data) shown as δ18O anomalies (relative to the mean of the time series) (black) with the y axis reversed and the AISMR anomalies (11-year running mean) (shaded) expressed as % departure from the 20th century mean (~850 mm). (b) The NI speleothem time series shown as raw δ18O anomalies (black) and NH temperature anomalies reconstruction36 (red). Bold lines (purple) are estimated trends43 in the time series over the last 1,000 years. White circles mark the break points (timing of trends change). Full size image

Spectral analysis31 of the record reveals significant power in the 60–80- and 15–30-year bands (Fig. 3c). The time-frequency32 and wavelet analyses33 of the speleothem δ18O record (see Methods and Supplementary Fig. 10) suggest that the amplitudes of aforementioned scale-averaged oscillations are non-stationary. For example, the 60–80-year oscillatory band is dominant during the AD 300–500, 1500–1700 and 1850–2000 intervals and nearly absent during other times. Significant concentration of variance in the 60–80-year band is also common in numerous instrumental and proxy global and regional sea surface temperature (SST) records34. The non-stationary aspect of the ISM multi-decadal oscillatory mode likely arises from the influence of specific (and time–space varying) arrays of different modes of SST variability, which themselves are coupled in complex manners35. On multi-centennial and longer timescales, the δ18O profile over the last 1,000 years is characterized by a pronounced trend towards higher values.

Figure 3: Time series analysis of the NI and CI δ 18 O records. (a) The NI speleothem record shown as δ18O anomalies (relative to the mean of the time series). The δ18O anomalies are shown both as raw data (grey) and smoothed (11-year running mean) (black) along with the regressed AISMR anomalies (%). (b) The comparison between the SSA44 detrended NI (black) and CI14,15,16 (red) speleothem records shown as δ18O anomalies (raw data). The long-term non-stationary trends in both time series are removed by subtracting the first reconstructed component indicated by SSA of the raw data. The two dashed horizontal lines delineate a 10% change in monsoon rainfall amounts, highlighting the magnitude of multi-decadal variability inferred from our NI δ18O record. (c) Power spectrum of the composite NI and (d) CI SSA-detrended δ18O time series obtained using REDFIT31 software. A varying number of Welch overlaps were used to optimize bias/variance properties. Spectral band significant above the 90% level are labelled with their period. Full size image

Comparison with a central India speleothem record

The Sahiya Cave δ18O record from NI shares key similarities and notable differences with the speleothem δ18O record from CI14,15,16 (19°N and 82°E) (Fig. 3b; Supplementary Fig. 11). The latter has been shown to reflect the ISM rainfall variability over the core monsoon zone of India (defined here as ~18°–27°N and 69°–88°E). Importantly, both records exhibit strong multi-decadal variability and share common spectral power (Fig. 3). Although the CI δ18O record has slightly larger amplitude and is more variable than the NI δ18O record, such differences may result, at least partially, from a varying degree of smoothing of isotopic signals imparted by processes unique to local karst at each location. The visual comparison between the two records shows similar long-term patterns but no consistent peak-to-peak relationship on decadal-multi-decadal timescales. For example, the interval with most prominent positive anomalies in the CI record (AD ~1370–1400) has no clear counterpart in the NI record, whereas the interval with most positive δ18O anomalies in the NI record (AD ~1600–1630) is clearly evident in the CI record (Supplementary Fig. 11). While such inconsistencies may result from spatial contrast in ISM rainfall, the combined age uncertainties of both records preclude definitive assessment of phase relationships between these records on short timescales.