Trends and oscillations in the Indian summer monsoon rainfall over the last two millennia

Abstract

Observations show that summer rainfall over large parts of South Asia has declined over the past five to six decades. It remains unclear, however, whether this trend is due to natural variability or increased anthropogenic aerosol loading over South Asia. Here we use stable oxygen isotopes in speleothems from northern India to reconstruct variations in Indian monsoon rainfall over the last two millennia. We find that within the long-term context of our record, the current drying trend is not outside the envelope of monsoon's oscillatory variability, albeit at the lower edge of this variance. Furthermore, the magnitude of multi-decadal oscillatory variability in monsoon rainfall inferred from our proxy record is comparable to model estimates of anthropogenic-forced trends of mean monsoon rainfall in the 21st century under various emission scenarios. Our results suggest that anthropogenic-forced changes in monsoon rainfall will remain difficult to detect against a backdrop of large natural variability.

ARTICLE

Received 14 Aug 2014 | Accepted 14 Jan 2015 | Published 17 Feb 2015

Ashish Sinha1, Gayatri Kathayat2, Hai Cheng2,3, Sebastian F.M. Breitenbach4,5, Max Berkelhammer6, Manfred Mudelsee7, Jayant Biswas8 & R.L. Edwards3

Observations show that summer rainfall over large parts of South Asia has declined over the past ve to six decades. It remains unclear, however, whether this trend is due to natural variability or increased anthropogenic aerosol loading over South Asia. Here we use stable oxygen isotopes in speleothems from northern India to reconstruct variations in Indian monsoon rainfall over the last two millennia. We nd that within the long-term context of our record, the current drying trend is not outside the envelope of monsoons oscillatory variability, albeit at the lower edge of this variance. Furthermore, the magnitude of multi-decadal oscillatory variability in monsoon rainfall inferred from our proxy record is comparable to model estimates of anthropogenic-forced trends of mean monsoon rainfall in the 21st century under various emission scenarios. Our results suggest that anthropogenic-forced changes in monsoon rainfall will remain difcult to detect against a backdrop of large natural variability.

DOI: 10.1038/ncomms7309

Trends and oscillations in the Indian summer monsoon rainfall over the last two millennia

1 Department of Earth Sciences, California State University Dominguez Hills, Carson, California 90747, USA. 2 Institute of Global Environmental Change, Xian Jiaotong University, Xian 710049, China. 3 Department of Earth Sciences, University of Minnesota, Minneapolis, Minnesota 55455, USA. 4 Department of Earth Sciences, Eidgenssische Technische Hochschule, CH-8092 Zrich, Switzerland. 5 Institute for Geology, Mineralogy & Geophysics, Ruhr-Universitat Bochum, Universitatsstrasse 150, D-44801 Bochum, Germany. 6 Department of Earth and Environmental Sciences, University of Illinois, Chicago, Illinois 60607, USA. 7 Alfred Wegener Institute for Polar and Marine Research, 27570 Bremerhaven, Germany. 8 National Cave Research and Protection Organization, Raipur, CG 492001, India. Correspondence and requests for materials should be addressed to A.S. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 6:6309 | DOI: 10.1038/ncomms7309 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7309

Among the most serious impacts of climate change to South Asia are those that relate to water availability and drought due to changes in Indian summer monsoon (ISM) rainfall

amount and spatio-temporal pattern. Although, the spatial structure of the mean seasonal rainfall over South Asia exhibits complex patterns1, the area-weighted all-India summer monsoon rainfall (AISMR) time series2 serve as an adequate index of ISM intensity over continental India on interannual to decadal timescales3,4. Between AD 1870 and 1990, the AISMR appears to exhibit quasi-oscillatory multi-decadal (B6070 years)

variability, with amplitudes that are generally within 5.0% of historical normals (B850 mm)46. However, in context of the last 5060 years, the AISMR has been undergoing a statistically signicant decreasing trend, which is also observed in various gridded precipitation data sets79. This decline stands in contrast to the near-unanimous results from climate models, which indicate increasing rainfall in a warming climate due to increasing greenhouse gas concentrations10,11. In some model simulations, the observed summertime drying trend has been attributed to increased anthropogenic aerosol loading over South Asia8,9,12, but considerable uncertainty remains in modelling the role of aerosols on monsoon rainfall13. A fundamental question therefore emerges whether the weakening of ISM rainfall, particularly over northern India (NI), is anomalous within the long-term context of natural variability. Deciphering whether the recent decade-long ISM drying trend is caused by human-induced aerosol emissions or a manifestation of ISM natural variability is critical to understanding both the near and long-term behaviour of the ISM.

Although the AISMR is among the longest instrumental records of precipitation on Earth, it is still too short (B145 years)

to allow a robust characterization of multi-decadal-centennial ISM variability from which, one can examine the uniqueness of the observed drying trend. Existing proxy reconstructions from South Asia suggest that the ISM exhibits complex spatio-temporal variations, the ner details of which are highlighted by speleothem oxygen isotope (d18O) records from central India (CI) and northeastern India1418. Building upon our previous work, here we present a new high-resolution and radiometric-dated speleothem d18O-based reconstruction of ISM rainfall over the last two millennia. The new record is developed from two stalagmites (SAH-A and SAH-B) collected from the Sahiya Cave in NI. The climatologically unique location of our site allows us to cast temporal changes in speleothem d18O in terms of decadal-averaged variations in ISM rainfall, and thus, the low-frequency aspect of our record can be construed as a paleo-counterpart of the instrumental AISMR time series. Our reconstruction reveals a wide range of ISM behaviour, consisting of multi-centennial trends forced by Northern Hemisphere (NH) temperatures, distinct multi-decadal oscillatory variability and sporadic subdecadal intervals of anomalously high and low rainfall. Within the long-term context of our record, the observed decline in ISM rainfall since the middle of the 20th century is not outside the envelope of natural variability, albeit at the lower edge of this variance. In addition, the magnitude of ISM multi-decadal oscillatory variability inferred from our record is comparable to model estimates of anthropogenic-forced trends of mean monsoon rainfall in the 21st century under various representative concentration pathway (RCP) scenarios11. Our results underscore the challenge of detecting anthropogenic-forced changes in ISM against a background of large internal variability.

ResultsRegional climatology. The Sahiya Cave ((30360N, 77520E, B1,190 m above sea level (m.a.s.l.)) is located B200 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 B22.0 C and 1,600 mm, respectively, with B80% 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 15N) 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 ow, 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 18N 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 spells2022. Distinct from the LPS, these oscillatory modes of monsoon exhibit a hierarchy of quasi-periods (37 days, 1020 days and 3060 days) and are manifestations of northsouth uctuations in the mean position of the MT2022. 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, inuence the strength of the BoB branch of ISM.

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 d18O of precipitation (d18Op) 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 d18Op observations (AD 20042006; n 28) (Supplementary Fig. 3). Consequently,

we use the d18Op data from the GNIP station at New Delhi(28.58N, 77.52E; 212 m.a.s.l., and B200 km south of Sahiya) to characterize the d18Op variability in our study area (Supplementary Figs 45). Our approach to utilize d18Op data from New Delhi is justied by the gross similarities between the annual cycles of precipitation and d18Op at Sahiya and New Delhi (see Methods), which suggest that both sites are inuenced by similar precipitation dynamics and moisture sources throughout the year. Furthermore, to examine dynamical linkages between d18Op and changes in the low-level (B850 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 d18Op data and the low-level monsoon moisture trajectory patterns (Fig. 1; Supplementary Fig. 6) suggest that the d18Op variability in the

2 NATURE COMMUNICATIONS | 6:6309 | DOI: 10.1038/ncomms7309 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7309 ARTICLE

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 nd that periods of strong (weak) monsoon circulation are characterized by an enhanced (reduced) ux of isotopically depleted BoB moisture and reduced (enhanced) ux 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 d18Op, signatures, which in theory, should predominantly drive changes in d18O of speleothems from the study area. This inference rests in part on the assumption that non-monsoon rainfall and snowmelt are insignicant sources of surface waters that inltrate into the cave system. Indeed, this assumption is supported by the d18O values of the caves dripwaters ( 6.950.31%, n 7), which are indistinguishable from the

amount-weighted d18O of annual precipitation ( 6.97%) and

slightly higher than the d18O 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 inltration of surface waters into the local karst occurs primarily during the ISM season (Supplementary Fig. 3). The karstic overburden over the Sahiya Cave (B30 m)

would, however, tend to act as a lowpass lter by mixing d18Op signals from several seasons, thus damping high-frequency d18Op variations. Nonetheless, the low-frequency (multi-year to

decadal) component of the speleothem d18O record should reect time-integrated shifts in d18O 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 d18O at our site to reect spacetime-integrated upstream changes in ISM intensity. Indeed, comparisons of the speleothem d18O 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 d18Oprainfall 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 d18O variability at our site. Our observations are consistent with the results from other studies, which suggest that local d18Op variability in some monsoonal and tropical regions reect upstream changes in regional circulation and rainfall26,27.

Speleothem oxygen isotope record. The NI d18O record from the Sahiya Cave is a composite of two d18O stalagmite records (SAH-A and SAH-B) consisting of 1848 d18O measurements with an average temporal resolution of 1.5 years per sample (Supplementary Fig. 8). The chronologic frameworks of SAH-A and SAH-B d18O proles are provided by 13 and 15 230Th dates, respectively (Supplementary Fig. 9; Supplementary Data 1). The average 2s age uncertainty of the composite record is B22.5 years. The age models are developed by generating multiple cubic interpolations28 between the discrete 230Th dates. The SAH d18O record (mean 9.0%; 1s 0.24) covers the last two thousand

40E

60E 80E 100E 120E 40E 60E 80E 100E 120E

40E 60E 80E 100E 120E 40E 60E 80E 100E 120E

June

July

10N

August

30N

September

30N

10N

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 ux from the Bay of Bengal, resulting in relatively higher ux 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 signicantly enhanced transport of 18O-depleted moisture from the Bay of Bengal (see Methods and Supplementary Fig. 6).

NATURE COMMUNICATIONS | 6:6309 | DOI: 10.1038/ncomms7309 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7309

years (BC 141 to AD 2006) (Supplementary Data 2). The d18O prole is characterized by prominent decadal to multi-decadal variability with an amplitude of B0.5% around the long-term mean (Fig. 2b). The observed weakening trend in ISM rainfall since the 1950s manifests as an B0.5% increase in speleothem d18O, which is a robust feature of our record even with generous assumptions of age model and analytical uncertainties. Comparison between the composite d18O prole and a 11-year running mean of AISMR exhibits a signicant inverse correlation (r 0.75 with 95% condence interval ( 0.82; 0.64)) that

explains B50% of variance in the instrument record (Fig. 2a; Supplementary Fig. 7). After correcting for standard instrumental error in d18O measurements (0.06%), the slope of the speleothem d18O-AISMR regression yields a conversion of 1.0% in speleothem d18O values to B86 mm change in AISMR, the latter equivalent to B10% of the long-term mean (B850 mm) (Fig. 2a). The d18O record is characterized by several episodic intervals of anomalously (42s) positive and negative d18O values. Particularly noteworthy are two sub-decadal length intervals of strong pluvial and drier conditions, centred at AD B1,100 and 1,620, respectively (Fig. 2b). The timing and duration of the latter event, which constitutes the largest positive d18O 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 10001200) and the Little Ice Age (nominally AD B16001850) (Fig. 2), respectively.

Spectral analysis31 of the record reveals signicant power in the 6080- and 1530-year bands (Fig. 3c). The time-frequency32 and wavelet analyses33 of the speleothem d18O record (see Methods and Supplementary Fig. 10) suggest that the amplitudes of aforementioned scale-averaged oscillations are non-stationary. For example, the 6080-year oscillatory band is dominant during the AD 300500, 15001700 and 18502000 intervals and nearly absent during other times. Signicant concentration of variance in the 6080-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 inuence of specic (and timespace varying) arrays of different modes of SST variability, which themselves are coupled in complex manners35. On multi-centennial and longer timescales, the d18O prole over the last 1,000 years is characterized by a pronounced trend towards higher values.

Comparison with a central India speleothem record. The Sahiya Cave d18O record from NI shares key similarities and notable differences with the speleothem d18O record from CI1416 (19N and 82E) (Fig. 3b; Supplementary Fig. 11). The latter has been shown to reect the ISM rainfall variability over the core monsoon zone of India (dened here as B1827N and 6988E). Importantly, both records exhibit strong multi-decadal variability and share common spectral power (Fig. 3).

AISMR anomaly (%)

Year (AD)

1900

1920

1940

1960

1980

2000

0.4

4

NI18O anomaly ()

0.2

0.4

0.2

2

0

0

2

4

10.0

9.5

0.00

NH temperature anomaly (C)

NI18O anomaly ()

9.0

0.25

8.5

8.0

0.50

7.5

200 200 400 600 1000 1200 1400 1600 1800 2000

800

Year (AD/BC)

Figure 2 | Comparison of the NI speleothem d18O record with instrumental and proxy data. (a) The comparison between the NI speleothem time series (raw data) shown as d18O 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 (B850 mm). (b) The NI speleothem time series shownas raw d18O 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).

4 NATURE COMMUNICATIONS | 6:6309 | DOI: 10.1038/ncomms7309 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7309 ARTICLE

Year AD/BC

200 200 400 600 800 1000 1200 1400 1600 1800 2000

18O anomaly () 18 O anomaly ()

0.5

+5.0

5.0

0.0

0.0

AISMR anomaly (%)

0.5

0.5

+5.0

5.0

0.0

0.0

Change in JJAS rainfall (%)

0.5

200 200 400 600 800 1000 1200 1400 1600 1800 2000

Year AD/BC

10

5

NI 18O time series

CI 18O time series

8

4

95% 90%

Spectral amplitute

Spectral amplitute

3

~70 year

6

95% 90%

~65 year

~33 year

~20 year~15 year

2

4

~20 year

1

2

0

0

0.0 0.02 0.04 0.06 0.08 0.10

0.00 0.02 0.04 0.06 0.08 0.10

Frequency (year1)

Frequency (year1)

Figure 3 | Time series analysis of the NI and CI d18O records. (a) The NI speleothem record shown as d18O anomalies (relative to the mean of the time series). The d18O 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 CI1416 (red) speleothem records shown as d18O anomalies (raw data).

The long-term non-stationary trends in both time series are removed by subtracting the rst 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 d18O record. (c) Power spectrum of the composite NI and (d) CI SSA-detrended d18O time series obtained using REDFIT31 software. A varying number of Welch overlaps were used to optimize bias/variance properties. Spectral band signicant above the 90% level are labelled with their period.

Although the CI d18O record has slightly larger amplitude and is more variable than the NI d18O 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 B13701400) has no clear counterpart in the NI record, whereas the interval with most positive d18O anomalies in the NI record (AD B16001630) 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 denitive assessment of phase relationships between these records on short timescales.

DiscussionThe low-frequency (multi-centennial) aspect of our ISM reconstruction from NI is broadly similar to NH temperature

reconstructions36 over the last two millennia (Fig. 2b). The close correspondence between these reconstructions on such a timescale supports the notion that NH temperature is an important driver of ISM variability via changes in meridional surface/tropospheric thermal gradients12,37 and the inter-hemispheric thermal contrast6,37. However, neither proxy nor instrumental data comparisons between ISM rainfall and NH temperatures reveal consistent relationships on decadal to multi-decadal timescales. This suggests that the thermodynamic effect on ISM rainfall maybe more important on longer timescales via changes in boundary conditions at a larger spatial scale (global or hemispherical)37. At shorter timescales, the ISM is more strongly inuenced by changes in local and/or remote SST elds and other aspects of external boundary conditions46. Nonetheless, given that the magnitude of the NH temperature rise since the 1950s is unprecedented in the context of the last 2,000 years, the decreasing trend in the ISM rainfall seems anomalous in light of changes expected from a simple thermodynamical-scaling argument38. The latter postulates that atmospheric water vapour

NATURE COMMUNICATIONS | 6:6309 | DOI: 10.1038/ncomms7309 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7309

should increase in line with the Clausius Clapeyron relationship at B7%C 1 rise in atmospheric temperature due to increasing concentrations of greenhouse gases. The lack of an upward trend in ISM rainfall since the 1950s supports results from coupled climate models that include aerosol forcing8,9,12. These model experiments indicate that the increased anthropogenic aerosol loading over south Asia may have substantially counterbalanced a precipitation increase over continental India that would otherwise have occurred in response to singular increase in global temperatures.

Parsing the relative contribution of aerosol forcing to the observed decline in ISM rainfall requires establishing an envelope of ISM baseline natural variability on decadal to multi-decadal timescales. We use the observed inverse speleothem d18O

AISMR relationship in deriving, admittedly rough, quantitative estimates of decadal-scale variations in ISM. Using the standard linear regression tted by means of the ordinary least-square method, we use the slope of the speleothem d18OAISMR relationship (B1.0%/86 mm) to estimate decadal-scale ISM rainfall variability over the last 2,000 years (Fig. 3). Although more sophisticated transfer functions employing a forward modelling approach may ultimately be more robust, any model for which the d18O of monsoon precipitation (and consequently, the d18O of speleothem) is inversely linked to area-weighted changes in ISM rainfall, should lead us, at least qualitatively, to the conclusion we reach as follows. Our analysis shows that the magnitude of ISM rainfall decrease since the mid-20th century is not unique in the long-term context of our record, but nonetheless right near the edge of the envelope of ISM long-term variability (Fig. 3a). This observation implies that the observed drying can be due to ISM natural variability alone. However, it will be imprudent to ignore the possible inuence of aerosol forcing on ISM rainfall, especially considering that previous excursions to comparably drier intervals inferred from our record typically occurred in the context of cooler NH temperatures and in absence of signicant human-induced aerosols emissions. It is also instructive to compare the maximum expression (peak-to-trough) of ISM multi-decadal variability in our record (B10%), for example, from AD 1720 to B1800 (Fig. 3b) with anthropogenic-forced trends of mean monsoon rainfall projected from climate models for the 21st century. For example, a recent evaluation of ISM rainfall for the period 18502100 from 20 climate models within the Coupled Model Intercomparison Project 5 showed11 a consistent increase in mean monsoon rainfall for the period AD 20702100 (relative to the historic period AD 18701900) under different RCP scenarios, with the highest increase in the strongest warming scenario RCP 8.5. However, only half of these models can realistically simulate the ISM11. The projected increase in mean ISM rainfall for the period AD 20702100 (relative to AD 18701900) from all but one of the realistic models is within 10% for all warming scenarios, which is similar to the maximum magnitude of ISM natural multi-decadal variability inferred from our record. These observations highlight the fundamental problem of detecting forced trends in ISM rainfall against the background of relatively large internal variability.

Understanding the current and future changes in ISM rainfall spatio-temporal pattern and amount has enormous societal implications for more than 1.5 billion people in the region. The paleoclimate data presented here provide a baseline of ISM natural variability on multi-decadal to centennial timescales. It is on precisely these timescales that the anthropogenic-forced changes are likely to be superimposed on monsoons natural climate variability. Detection and attribution of climate change on ISM rainfall will require a coordinated effort through improvements in both climate models and proxy reconstructions of ISM

rainfall. The paper presented here highlights the promise of this approach and provides impetus for further renement of ISM rainfall reconstructions.

Methods

Cave characteristics and sample collection. In April 2007, we collected a pair of columnar-shaped stalagmites (SAH-A and SAH-B) from a small and poorly ventilated chamber, B100 m from the caves entrance. The stalagmite SAH-A was found wet (and presumably actively growing) at the time of collection. In situ measurements indicate relative humidity in excess of 95%, numerous active feeder drips and a narrow range of ambient temperatures (Tcave air 21.80.6 C, n 7)

throughout the cave, with the latter being similar to the regional mean annual air temperature (B22.0 C) (Supplementary Fig. 2). Unlike higher altitude regions to its north, the study area (elevation B1,0001,200 m) does not experience snowfall. Consequently, spring snowmelt is not a source of meteoric water into the local karst above the Sahiya Cave. The regional precipitationevaporation is positive only during the summer monsoon season (Supplementary Fig. 2b), suggesting that effective inltration of meteoric waters into the local karst is only signicant during the ISM season.

Sample preparation. The stalagmites were cut with a thin diamond blade along the growth axes. A digitally controlled triaxial micromill was used to obtain sub-samples for d18O and d13C analyses by continuously milling in 100-mm increments along the stalagmites growth axes, yielding sub-samples weighing typically between 100 and 130 mg for isotopic analyses. Sub-samples for 230Th dating (weighing typically between 150 and 300 mg) were obtained by drilling with a 0.9-mm carbide dental burr.

Stable isotope analysis. Sub-samples for d18O and d13Canalyses were measured using a Multi-Prep Dual Inlet System attached to an Iso-prime Stable Isotope Ratio Mass Spectrometer at the University of Southern California. The long-term precision of the Ultissima standard for d18O and d13C measurements over the course of this study was B0.06% and 0.03%, respectively.

Isotopic equilibrium. Two lines of evidence support that calcite deposition in the Sahiya cave occurs at or near isotopic equilibrium. First, strong replication between the SAH-A and SAH-B d18O records during the period of overlap (r 0.60 with

95% condence interval (0.17, 0.84)) (Supplementary Fig. 8) suggest that both stalagmites have recorded a primary climate signal that has not been signicantly altered by disequilibrium processes associated with conditions unique to each stalagmite, such as different ow pathways through the overlying bedrock, groundwater mixing and residence time, pCO2 and so on. Second, the theoretical range of the caves ambient temperature (20.823.9 C) under isotopic equilibrium conditions39 calculated from the measured range in d18O of dripwaters( 6.950.31%, n 7) and the d18O of modern calcite ( 8.53%) is well within

the measured range of ambient temperature in the Sahiya cave (21.80.6 C, n 7).

230Th dating and the age model. 230Th dating was performed at the Minnesota Isotope Laboratory on a Multi-Collector Inductively Coupled Plasma Mass Spectrometer (Thermo-Finnigan Neptune). We use standard chemical procedures to separate uranium and thorium for 230Th dating, as reported in Edwards et al.40 The standardization and half-lives are reported in Cheng et al.41 All U and Th isotopes were measured on a MasCom multiplier behind the retarding potential quadrupole in the peak-jumping mode. We followed similar procedures of characterizing the multiplier as described in Shen et al.42 Initial 230Th values were corrected and all ages are reported with a 2s error and are in correct stratigraphic order (Supplementary Data 1). The chronologic frameworks of stalagmite SAH-A and SAH-B d18O proles are constrained by 13 and 15 230Th dates, respectively.

We use Monte Carlo simulations (N 2,000) to account for the dating

uncertainties (2.5% and 97.5% quantile condence limits) within the time series28 (Supplementary Fig. 9; Supplementary Data 1).

Time series analyses. Condence intervals of correlation coefcients were determined using the pairwise moving-block bootstrap method to preserve the serial dependence of time series and interval calibration to increase the accuracy43. To identify periodic components in the spectrum of unevenly spaced SAH-A, SAH-B and composite SAH d18O time series, we applied the LombScargle Fourier transform in combination with Welchs Overlapped Sequence Averaging procedure using the REDFIT software31. REDFIT produces rst-order autoregressive (AR1) time series with sampling times and characteristic timescales matching those of the input climate data. A varying number of overlapping (50%) detrended segments were used to optimize bias/variance properties of spectrum estimates. The univariate spectra were bias-corrected using 1,000 Monte Carlo simulations. To assess the statistical signicance of a spectral peak, REDFIT estimates the upper condence interval of the AR1 noise for various signicance levels based on a w2-distribution. Temporal changes in the average amplitude of signal components

6 NATURE COMMUNICATIONS | 6:6309 | DOI: 10.1038/ncomms7309 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7309 ARTICLE

were estimated using a modied harmonic-ltering algorithm32, which ts sinusoidal waves to an unevenly spaced series by means of least squares (analogous to the LombScargle Fourier transform). We used a sliding window of length four times the signal period and ascribe to the windows midpoint the average signal amplitude; the window is then shifted in steps of 0.1 year. This leads to a truncation of the time interval for the resulting estimate of the average amplitude of 4 80/2

years 160 years, from both the lower and upper time interval bound. In addition,

we also applied the conventional continuous transform wavelet analysis to annually interpolated d18O data being mindful of the caveat that here interpolation might introduce artefacts, potentially biasing the derived spectra. The singular spectrum analysis-Multitaper Method (SSA-MTM) Toolkit44 was used to remove non-stationary linear trends from the d18O time series by examining the variability of the largest eigenvalue from an autocovariance matrix. SSA decomposes time series into several modes of variability (monotonic and oscillatory) and the remaining low-amplitude noise. Time series trends and associated break points (timing of trends change) are quantied using a parametric, nonlinear regression technique43. This model is a continuous function consisting of two linear parts that are joined at the break point. The break model is tted to data using a weighted least-squares method with a brute-force search for the break point. Statistical uncertainties of the trend and break point are evaluated using 2,000 block bootstrap simulations, which preserve the distribution and serial dependence of the data over the length of a block.

Statistical analysis of observational d18Op data. We extensively utilized observational monthly d18Op data from two sitesUttarkashi and New Delhi (data sets are accessible from http://www.naweb.iaea.org

Web End =www.naweb.iaea.org ) to characterize the d18Op variability and understand the dynamical controls of the d18Opprecipitation (d18OpP) relationship in our study area. The d18Op data from Uttarkashi (B50 km from the Sahiya Cave) contains only 3 years of monthly d18Op observations (AD 20042006; n 28)45, whereas the GNIP station at New Delhi (B200 km south of the Sahiya

Cave), contains 341 monthly d18Op observations (available with gaps from AD 1961 to 2008). Despite the large difference in number of observations, the annual cycles of the d18Op and rainfall amount from both sites are remarkably similar (Supplementary Fig. 2ab). This close correspondence is not surprising considering that both sites are under the inuence of similar precipitation dynamics due to their proximity. In addition, we nd strong correlation between the observational d18Op data from New Delhi and simulated d18Op data (extracted from the grid point nearest to the Sahiya Cave) from Experimental Climate Prediction Centers Global Spectral Model (ECPC-GSM) tted with isotope tracers (IsoGSM)46 (Supplementary Fig. 3). In light of these observations, we argue that the GNIP d18Op database from New Delhi data can be reliably used to understand the nature of the d18OpP relationship in our study area.

The New Delhi GNIP database is characterized by a large range in both monthly d18Op and precipitation amount (Supplementary Fig. 4). The distribution of JJAS-averaged d18Op and precipitation amount illustrates a negative d18OpP relationship (r 0.37, with 95% condence interval ( 0.62; 0.04), n 35).

The inverse d18OpP relationship for the individual months JJAS (Supplementary Fig. 5) is most pronounced at the onset (June) and end (September) of the ISM season. In contrast, the inverse d18OpP relationship during the height of the monsoon season (JA) is relatively weaker despite July and August being typically the wettest months. One plausible explanation for this observation is that lower humidity, higher evaporation and unsaturated soil conditions in June and September lead to faster re-evaporation of precipitation below cloud base and from the ground, thus explaining the observed inverse d18OpP relationship through the classic amount effect. However, a critical examination of the GNIP data set, together with composite trajectory analyses of low-level monsoon winds during strong and weak monsoon seasons, suggest that sub-seasonal changes in monsoon circulation and associated changes in the moisture source and rainout history of air parcels, are more dominant contributors to the observed inverse d18OpP relationship in our study area. This can be explained as follows. The monthly observations of d18Op and precipitation amount in the GNIP data exhibit a large range with some values representing extreme ends of both, d18Op and rainfall amount ranges. We dene the extreme values as those that fall in the lower and the upper quartiles of observed range of d18Op and precipitation amount (Supplementary Fig. 5). On the basis of this scheme, we classify these extreme d18OpP pairs (and their relative abundance in the GNIP database) into four groups: (a) negative-wet (NW, relatively common); (b) negative-dry (rare);(c) positive-wet (rare) and (d) positive-dry (PD, relatively common). We nd that the presence or absence of the inverse d18OpP relationship is strongly sensitive to the abundance of the NW and PD pairs in the data set, and if these data pairs are removed from the regression analysis, the inverse d18OpP relationship during

June weakens and vanishes for the remaining summer monsoon season. These observations suggest that the primary source of the inverse d18OpP relationship in the GNIP database arises from occurrence of anomalously strong (18O depleted) and weak (18O enriched) precipitation events.

HYSPLIT trajectory analysis. We utilized the HYSPLIT25 model provided by the National Oceanographic and Atmospheric Administration Air Resources Laboratory together with the National Center for Environmental Prediction and National Center for Atmospheric Research Reanalysis version 1 data47 to generate

quasi-isentropic backward (168 h) trajectory ensembles of strong and weak monsoon seasons (Fig. 1; Supplementary Fig. 6). We are cognizant that HYSPLIT-derived trajectories can have large computational errors that are inherited from uncertainties in the reanalysis products. Therefore, we increase our condence in trajectory analysis by computing hundreds of trajectories over an extended period to identify major moisture source areas and broad circulation features that produce precipitation in our study area. For each month of the monsoon season (JJAS), we produce the wet and dry trajectory ensembles (Supplementary Fig. 6). The wet and dry months are dened here as those that fall in the upper and lower quartile of the observed range in precipitation amount in the GNIP database. For each wet and dry composite, we computed B600 isentropic back-trajectories(168 h)(sampled four times daily at UTC 00, 06, 12 and 18) initialized at 1,500 metres (above ground level corresponding to B850 hPa). Results from a smaller sub-set of randomly selected trajectories initialized at 1.0 or 1.5 km above ground level (corresponding to B900 and 850 hPa levels, respectively) and computed for either 120 or 168 h duration, were nearly identical.

References

1. Hoyos, C. D. & Webster, P. J. The role of intraseasonal variability in the nature of Asian monsoon precipitation. J. Climate 20, 44024424 (2007).

2. Parthasarathy, B., Munot, A. A. & Kothawale, D. R. All-India monthly and seasonal rainfall series18711993. Theor. Appl. Climatol. 49, 217224 (1994).

3. Kripalani, R. H. & Kulkarni, A. Monsoon rainfall variations and teleconnections over South and East Asia. Int. J. Climatol. 21.5, 603616 (2001).

4. Krishnamurthy, V. & Goswami, B. N. Indian monsoon-ENSO relationship on interdecadal time scales. J. Climate 13, 579595 (2000).

5. Kriplani, R. H., Kulkarni, A. A., Sabade, S. S. & Khandekar, M. L. Indian monsoon variability in a global warming scenario. Nat. Hazards 29, 189206 (2003).

6. Goswami, B. N., Madhusoodanan, M. S., Neema, C. P. & Sengupta, D.

A physical mechanism for North Atlantic SST inuence on the Indian summer monsoon. Geophys. Res. Lett. 33, L02706 (2006).

7. Turner, A. G. & Annamalai, H. Climate change, and the South Asian summer monsoon. Nat. Clim. Change 2, 587595 (2012).

8. Ramanathan, V. et al. Atmospheric brown clouds: Impacts on South Asian climate and hydrological cycle. Proc. Natl Acad. Sci USA 102, 53265333 (2005).

9. Bollasina, M. A., Ming, Y. & Ramaswamy, V. Anthropogenic aerosols and the weakening of the South Asian summer monsoon. Science 334, 502505 (2011).

10. Kripalani, R. H. et al. South Asian summer monsoon precipitation variability: coupled climate model simulations and projections under IPCC AR4. Theor. Appl. Climatol. 90, 133159 (2007).

11. Menon, V., Levermann, A., Schewe, J., Lehmann, J. & Frieler, K. Consistent increase in Indian monsoon rainfall and its variability across CMIP-5 models. Earth Syst. Dynam. 4, 287300 (2013).

12. Polsen, D., Bollasina, M., Hegerl, G. C. & Wilcox, L. J. Decreased monsoon precipitation in the Northern Hemisphere due to anthropogenic aerosols. Geophy. Res. Lett. 41, 60236029 (2014).

13. Gautam, R., Hsu, N. C., Lau, K. M. & Kafatos, M. Aerosol and rainfall variability over the Indian monsoon region: distributions, trends and coupling. Ann. Geophys. 27, 36913703 (2009).

14. Sinha, A. et al. A 900-year (600 to 1500 AD) record of Indian summer monsoon precipitation from the core monsoon zone of India. Geophys. Res. Lett. 34, L16707 (2007).

15. Berkelhammer, M. et al. Persistent multidecadal power of the Indian Summer Monsoon. Earth Planet. Sci. Lett. 290, 166172 (2010).

16. Sinha, A. et al. The leading mode of Indian Summer Monsoon precipitation variability during the last millennium. Geophys. Res. Lett. 38, L15703 (2011).

17. Berkelhammer, M. et al. in Climates, Landscapes, and Civilizations (eds Giosan L., Fuller D. Q., Nicoll K., Flad R. K. & Clift P. D.) 7588 (American Geophysical Union, 2012).

18. Berkelhammer, M. et al. On the low-frequency component of the ENSOIndian monsoon relationship: a paired proxy perspective. Clim. Past. 10, 733744 (2014).

19. Joseph, P. V. & Sijikumar, S. Intraseasonal variability of the low-level jet stream of the Asian summer monsoon. J. Climate 17, 14491458 (2004).

20. Sikka, D. R. Some aspects of the large scale uctuations of summer monsoon rainfall over India in relation to uctuations in the planetary and regional scale circulation parameters. Proc. Indian Acad. Sci. (Earth Planet. Sci.) 89, 179195 (1980).

21. Goswami, B. N., Ajayamohan, R. S., Xavier, P. K. & Sengupta, D. Clustering of synoptic activity by Indian summer monsoon intraseasonal oscillations. Geophys. Res. Lett. 30, 1431 (2003).

22. Webster, P. J. et al. Monsoons: processes, predictability, and the prospects for prediction. J. Geophys. Res. 103, 1445114510 (1998).

23. Krishnan, R., Vinay, K., Sugi, M. & Yoshimura, J. Internal feedbacks from monsoonmidlatitude interactions during droughts in the Indian summer monsoon. J. Atmos. Sci. 66, 553578 (2009).

NATURE COMMUNICATIONS | 6:6309 | DOI: 10.1038/ncomms7309 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7309

24. Krishnan, R. et al. Will the South Asian monsoon overturning circulation stabilize any further. Clim. Dyn. 40, 187211 (2012).

25. Draxler, R. R. & Rolph, G. D. HYSPLIT (Hybrid Single-Particle Lagrangian Integrated Trajectory) Model access via NOAA ARLREADY Website http://www.arl.noaa.gov/ready/hysplit4.html

Web End =http:// http://www.arl.noaa.gov/ready/hysplit4.html

Web End =www.arl.noaa.gov/ready/hysplit4.html NOAA Air Resources Laboratory (2003).

26. Kurita, N., Ichiyanagi, K., Matsumoto, J., Yamanaka, M. D. & Ohata, T. The relationship between the isotopic content of precipitation and the precipitation amount in tropical regions. J. Geochem. Explor. 102, 113122 (2009).

27. Moerman, J. W. et al. Diurnal to interannual rainfall d18O variations in northern Borneo driven by regional hydrology. Earth Planet. Sci. Lett. 369-370, 108119 (2013).

28. Breitenbach, S. F. M. et al. Constructing proxy-record age models (COPRA). Clim. Past. 8, 23692408 (2012).

29. Maharatna, A. The Demography of Famines: an Indian Historical Perspective xviii 317 (Oxford Univ. Press, 1996).

30. Cheng, H., Edwards, R. L. & Haug, G. H. Comment on On linking climate to Chinese dynastic change: Spatial and temporal variations of monsoonal rain. Chinese Sci. Bull. 55, 37343737 (2010).

31. Schulz, M. & Mudelsee, M. REDFIT: estimating red-noise spectra directly from unevenly spaced paleoclimatic time series. Comp. Geosci. 28, 421426 (2002).

32. Feraz-Mello, S. Estimation of periods from unequally spaced observations. Astron. J. 86, 619624 (1981).

33. Torrence, C. & Compo., G. P. A practical guide to wavelet analysis. Bull. Am. Meteor. Soc. 79, 6178 (1998).

34. Delworth, T. L. & Mann, M. E. Observed and simulated multidecadal variability in the Northern Hemisphere. Clim. Dyn. 16, 661676 (2000).

35. Wyatt, M. G., Kravtsov, S. & Tsonis, A. A. Atlantic multidecadal oscillation and Northern Hemispheres climate variability. Clim. Dyn. 38, 929949 (2012).

36. Mann, M. E. & Jones, P. D. Global surface temperatures over the past two millennia. Geophys. Res. Lett. 30, 1820 (2003).

37. Wang et al. Northern Hemisphere summer monsoon intensied by mega-El Nino/southern oscillation and Atlantic Multidecadal oscillation. Proc. Natl Acad. Sci. USA 110, 53475352.

38. Held, I. & Soden, B. Robust responses of the hydrological cycle to global warming. J. Climate 19, 56865699 (2006).

39. Kim, S. T. & ONeil, J. R. Equilibrium and nonequilibrium oxygen isotope effects in synthetic carbonates. Geochim. Cosmochim. Acta 61, 34613475 (1997).

40. Edwards, R. L., Chen, J. H. & Wasserburg, G. J. 238U, 234U, 230Th, 232Th systematics and the precise measurement of time over the past 500,000 years. Earth Planet. Sci. Lett. 81, 175192 (1987).

41. Cheng, H. et al. Ice age terminations. Science 326, 24825214 (2009).

42. Shen, C.-C. et al. High-precision and high-resolution carbonate 230Th dating by MC-ICP-MS with SEM protocols. Geochim. Cosmochim. Acta 99, 7186 (2012).

43. Mudelsee, M. Climate Time Series Analysis: Classical Statistical and Bootstrap Methods Vol. 42, 474 (Springer, 2010).

44. Ghil, M. et al. Advanced spectral methods for climatic time series. Rev. Geophys. 40, 1003 (2002).

45. Kumar, B. et al. Isotopic characteristics of Indian precipitation. Water Resour. Res. 46, W12548 (2010).

46. Yoshimura, K., Kanamitsu, M., Noone, D. & Oki, T. Historical isotope simulation using reanalysis atmospheric data. J. Geophys. Res. 113, D19108 (2008).

47. Kalnay, E. et al. The NCEP/NCAR 40-year reanalysis project. Bull. Am. Meteor. Soc. 77, 437471 (1996).

Acknowledgements

We gratefully acknowledge the Chinese National Science Foundation grants to H.C. (NBRP 2013CB955902 and NSFC 41230524) and the National Science Foundation grants to A.S. (ATM: 0823554), R.L.E. and H.C. (1211299, 1103403 and 1337693) for funding this research. S.F.M.B. acknowledges nancial support from the Swiss National Foundation (Sinergia grant CRSI22-132646/1). We thank Lowell Stott, Miguel Rincon and Tannia Ochoa for the stable isotopic measurements.

Author contributions

A.S. H.C. and G.K. conceived the project and A.S. wrote the rst draft of the paper. H.C. R.L.E. and G.K. contributed to 230Th dating. S.F.M.B. and M.M. performed age modelling and statistical analyses. All authors discussed the results, edited and commented on the manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions

Web End =reprintsandpermissions/

How to cite this article: Sinha, A. et al. Trends and oscillations in the Indian summer monsoon rainfall over the last two millennia. Nat. Commun. 6:6309 doi: 10.1038/ncomms7309 (2015).

8 NATURE COMMUNICATIONS | 6:6309 | DOI: 10.1038/ncomms7309 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.