Interannual oxygen isotope variability in Indian summer monsoon precipitation reflects changes in moisture sources

The primary influences on the spatio-temporal variability of oxygen isotope compositions in precipitation over the Indian summer monsoon domain are inadequately constrained by the limited observational record. Consequently, the climatic significance of isotopic signatures of precipitation preserved in proxy archives from the region remains unclear. Here we present simulations with an isotope-enabled climate model (IsoGSM2) with the moisture-tagging capability to investigate the role of relative contributions of moisture from oceanic and terrestrial sources to the interannual variability in oxygen isotope composition in summer monsoon rainfall. During weak monsoon years, the moisture contribution from the Arabian Sea dominates precipitation over the Indian subcontinent while the remote oceanic and terrestrial sources have a greater influence during strong monsoon years. We suggest that changes in monsoon circulation, moisture source, and precipitation intensity are interrelated and that speleothem oxygen isotope records from the region can potentially help reconstruct interannual to decadal monsoon rainfall variability. A combination of monsoon circulation, moisture source dynamics, and precipitation intensity influences speleothem oxygen isotope records in the Indian monsoon region, according to simulations using an isotope-enabled model with moisture-tagging capability.

T he oxygen isotopic signatures of precipitation (δ 18 O p ) preserved in natural archives such as ice cores, tree ring cellulose, and speleothems are widely used for reconstructing past climatic conditions. Nonetheless, the relative importance of various hydrological cycle processes in driving the spatial-temporal δ 18 O p variability in the Asian monsoonal regions remains inadequately constrained. A case in point is the contentious debate surrounding the climatic significance of the speleothem δ 18 O records from the East Asian summer monsoon (EASM) region [1][2][3][4][5][6][7] . The earlier (c. the 2000s) interpretations of these records focused on the role of summer monsoon "intensity" or the ratio of summer to winter precipitation in driving δ 18 O changes on orbital timescales 8,9 . However, as more observational δ 18 O p data have since become available, together with the refinements in numerical modeling of water isotopes in general circulation model (GCM) [10][11][12] , the EASM δ 18 O p and speleothem δ 18 O variability are now suggested to arise from "upstream" processes such as seasonal changes in circulation 13,14 , moisture sources 15,16 , transport history 17 , and not necessarily from changes in the amount of precipitation at the local or regional scale (see review in ref. 18 ).
In contrast to their EASM counterparts, the speleothem δ 18 O records from the Indian summer monsoon (ISM) regions ( Fig. 1; Supplementary Fig. 1 and Supplementary Note 1) have largely escaped scrutiny but key elements of these records require a critical reappraisal 19 . For example, some speleothem δ 18 O records from India exhibit a rather large δ 18 O change (~6-8‰) across the millennial-scale abrupt climate events and glacial terminations 20,21 , which is two to three times larger than the coeval EASM speleothem records 22 . If interpreted solely in the framework of changes in upstream circulation and/or moisture sources, large millennial-scale δ 18 O shifts imply a massive reorganization of the ISM dynamics. Similarly, the late Holocene speleothem δ 18 O records from India exhibit protracted (subdecadal to multidecadal) intervals of extreme positive δ 18 O anomalies (~2-3 standard deviation excursions from the mean), which are interpreted to reflect periods of drastically reduced monsoon rainfall. Such reconstructions are, however, at odds with the ISM variability observed during the instrumental period, where there have never been more than five droughts in any continuous 10-year interval in India [23][24][25][26][27] .
Previous attempts to cast the observed spatio-temporal δ 18 O p patterns over the Indian subcontinent in the context of changes in circulation and moisture source-transport-sink relationships have largely relied on Lagrangian models 26,[28][29][30][31] . While this approach has provided some mechanistic insights into the δ 18 O p variability at synoptic scales 26,31 , its broader applications to seamlessly investigate long-term variations in moisture origin are limited by the brevity of observational δ 18 O p data. The Eulerianbased approach implemented within the isotope-enabled GCMs provides an alternative method to elucidate the role of moisture source dynamics on the δ 18 O p . This method involves "tagging" water vapor from predefined regions and tracing their movement during the hydrological cycle [32][33][34][35][36][37] . In this study, we have used an isotope-incorporated Global Spectral Model version 2 (IsoGSM2) 38 with water vapor-tagging capability (see "Methods") to illustrate the relative importance of large-scale moisture transport from various oceanic and terrestrial sources in driving the interannual δ 18 O p variability over the Indian subcontinent. Intercomparisons with other isotope-enabled models indicate that IsoGSM2 reproduces the annual cycles in precipitation and δ 18 O p over the ISM 39 and EASM regions 7 , arguably, with the highest fidelity. Importantly, the model simulates the intra-seasonal to interannual variability in δ 18 O p over the Indian subcontinent with high fidelity, which is directly comparable to the observational data (Supplementary Figs. 2-5; Supplementary Table 1, "Methods"), indicating that IsoGSM2 has high skill in reproducing moisture transport processes and hydrological cycle over the ISM domain.

Results
Regional setting. The onset of summer monsoon rainfall over the Indian subcontinent is marked by the development and intensification of cross-equatorial low-level jet (LLJ). The LLJ with a core at the height of~1.5 km (corresponding to~850 hPa level), is the dominant pathway of moisture transport across the equatorial Indian Ocean and the Arabian Sea onto the Indian landmass 1 (Supplementary Fig. 1). The LLJ-induces strong cyclonic vorticity in the planetary boundary layer over central India that leads to moisture convergence and the formation of the monsoon trough (MT across central India to the Bay of Bengal (BoB) (e.g., ref. 41 ). After exiting India, the LLJ reconstitutes over the northern BoB, turning into a southeasterly flow, transporting moisture to eastern and northeastern India, and occasionally as far northwest as Pakistan. Besides these two moisture transport pathways, the negative pressure anomalies over the MT draw synoptic-scale low-pressure systems (LPS) from the BoB 42 , which further transport oceanic vapor across central to northwest India. The LLJ mediated moisture transport has been shown to vary on both intraseasonal to interannual timescales [43][44][45] . Typically, during weaker monsoon, moisture transport along the LLJ climatological axis (~15°N) weakens and the LLJ also morphs into two branches, with one branch passing south of India and the other passes through north India, with its core at~25°N (Supplementary Fig. 1) 40 .
Tagging design. In this study, we tagged moisture from the following source regions: (1) The Bay of Bengal (BoB), (2) the Arabian Sea (including the Red Sea and the Gulf of Oman (hereafter, ARAB), (3) northern Indian Ocean (NIO), (4) southern Indian Ocean (SIO), (5) the western part of the Pacific Ocean (PCO), and (6) land surfaces (LND) (Fig. 1a). The atmospheric column precipitable water derived from tagged regions was used to identify the source region of the water vapor. Comparisons with climatological values of vertically integrated moisture flux divergence using the Fifth Generation European Centre for Medium-Range Weather Forecasts Reanalysis (ERA-5) suggests that our tagged moisture source regions encompass all potential sources of oceanic moisture based on the E-P threshold of 400 mm/season (Supplementary Fig. 6) 46,47 . The climatological contributions of moisture from each tagged region were evaluated over a common "sink" region over central India (~15°N to 28°N and 70-94°E), commonly referred to as the core monsoon zone (CMZ) 48 . The moisture contributions to three additional sink regions (4°× 4°area centered over each of the proxy site locations) in north, northeast, and central India were also assessed ( Supplementary Fig. 1c). The combined contributions of moisture originating from all six tagged sources to the total precipitable water (TPW) over the CMZ (Fig. 1c) and other sinks range between 96 and 100% during the pre-monsoon (April-May), monsoon (JJAS), and post-monsoon (October) seasons, indicating that the identified moisture source regions adequately represent the dominant sources of precipitation at each of our sink regions ( Supplementary Fig. 7). Extensive validation of model-simulated isotope data including reproducing the observed interannual, seasonal, and intra-seasonal variations in δ 18 O p 49 suggests that IsoGSM2 has high skill in reproducing moisture transport processes and hydrological cycle over the ISM domain.
Climatology of moisture contributions. The annual cycles in simulated moisture contributions to the TPW (longitudinally averaged from 70°E to 94°E) from each tagged region between −5°S to 40°N are illustrated by a Hovmöller diagram (Fig. 1b). During the pre-monsoon season (April-May), moisture contributions to the CMZ are dominated by the ARAB (~50%) and LND source regions (i.e., recycled vapor originating from land surfaces (~20%) (Fig. 1b, Supplementary Fig. 7). Between late May and early June, contributions from the NIO and SIO abruptly increases over the CMZ marking the intensification of the LLJ, with the SIO exhibiting the largest relative increase with respect to the pre-monsoon level ( Fig. 1b and Supplementary  Fig. 7). The JJAS climatological moisture contribution to the CMZ and other sink regions are dominated by the ARAB (~40%), LND (~26%), and SIO (2-0%) ( Fig. 1b and Supplementary   Fig. 7), peaking in the early (June), late (September), and mid (July-August) part of the monsoon season, respectively (Supplementary Fig. 8a). Between late August and early October, as the ISM withdraws, the relative moisture contributions from the ARAB, NIO, and SIO decrease whereas the contribution from the LND increases sharply (~45%) over the CMZ (Fig. 1b and Supplementary Fig. 7). Linear correlations between the de-seasoned anomalies of IsoGSM2 simulated monthly isotopic composition of precipitable water (δ 18 O pw ) and moisture contributions from each tagged source region over the CMZ and other sink regions are summarized in Fig. 1 and Supplementary Fig. 7 ARAB (r~0.6) and LND (r~−0.3) over the CMZ as well for the other sink regions ( Supplementary Fig. 7). In summary, the following key observations from our tagging simulation are noteworthy: (1) contributions from the ARAB, SIO, and LND (i.e., local recycling) tagged regions are important moisture sources of ISM precipitation in the early, mid, and late monsoon season, respectively. These results are in agreement with those estimated within the framework of Lagrangian dispersion models, although the boundaries of the source regions and their relative importance vary among studies 46,50-52 . (2) From a climatological perspective, moisture contribution from the BoB tagged region is not an important source of JJAS precipitation over much of continental India (except over northeast India). This finding, however, does not contradict the notion that a significant portion of the moisture is transported from the monsoon's southeasterly circulation and the LPS originating from the BoB because the transported moisture may have origins from other sources 53 . (3) JJAS moisture contributions from ARAB and LND significantly influence the δ 18 O pw and δ 18 O p over the Indian subcontinent.
Moisture contributions during weak and strong monsoon years. We used composite analysis to examine relative changes in moisture contributions from each tagged source over the Indian subcontinent during the 'weak and strong' monsoon years. The weak (1982, 1987, 2002, 2004, and 2009) and strong (1983, 1988, and 1994) years are identified as departures of one standard deviation from the All-India Rainfall (AIR) mean (~850 mm) between 1979 and 2010. The AIR is an area-weighted mean summer monsoon rainfall index based on a homogeneous rainfall data set of 306 rain gauges in India 54 and is widely considered as a reliable index of summer monsoon activity over the Indian region (e.g., ref. 24 ).
The differences between the weak and strong composites of JJAS precipitation and vertically integrated moisture flux (the sign of those differences is indicative of a weaker monsoon) highlight the following key features (Fig. 2a). The weak monsoon years are characterized by (i) negative precipitation anomalies over much of continental India except over northeast, northern BoB, and the coastal Myanmar regions; (ii) reduced vertically integrated zonal moisture flux along the climatological axis of the LLJ over the Arabian Sea, compensated by anomalously enhanced northwestsoutheast trending zonal moisture flux across north India (with its core at~25°N) and near the equator; (iii) negative meridional moisture flux anomalies, particularly east of 75°E over continental India. During the weak monsoon years, the moisture originating from the ARAB tagged region converges over northeast India and BoB due to the anomalous moisture flux and low-level circulation 55 . On other hand, moisture contribution from the SIO is substantially reduced over continental India and all along the LLJ. There is also a marked reduction in recycled moisture over north India but a large increase over the BoB and southeast Asia during the weak monsoon years (Fig. 2b). In addition, the moisture contributions from the BoB, NIO, and PCO regions are reduced to a varying degree over continental India. It is important to note that Fig. 2b highlights the spatial patterns of moisture contributions from each tagged region on a standalone basis. However, to assess how the δ 18 O pw and δ 18 O p vary during the weak and strong monsoon years, it is more instructive to examine the relative changes in moisture contribution ratio. When viewed in this context, the spatial patterns of moisture contribution ratios from the ARAB and LND during the weak monsoon years are uniformly higher over continental India (except over northwest India in case of LND), and lower for all other moisture sources to the TPW (Fig. 2c). These spatial patterns can be primarily attributed to a shifting balance of ARAB and LND moisture contribution between the weak and strong monsoon years. The higher values of ARAB over the Indian subcontinent during the weak years are mirrored by substantially enriched δ 18 O pw and δ 18 O p (Supplementary Fig. 8b) across the ISM domain, which can be attributed to the higher prevalence of isotopically enriched moisture from the proximal ARAB source region.
Interannual variability of moisture contributions, precipitation, and isotopes. We further examined how variations in moisture contributions from ARAB and LND, which show the largest change between weak and strong monsoon years (Fig. 2c), relate to interannual changes in the ISM precipitation and its δ 18 O pw and δ 18 O p . This is approached by examining the spatial pattern of correlation coefficients between the ARAB and LND with precipitation amount and δ 18 O p at each grid cell as well as through comparisons of their time series values averaged over a large rectangular-shaped region over continental India (Fig. 3). Here, we used IsoGSM2 simulated precipitation for this analysis, however, the results of this study are not sensitive to the choice of precipitation dataset (Supplementary Fig. 9). Spatial correlation analysis indicates an east-west dipole like-structure, where the JJAS moisture contributions from ARAB are negatively and positively correlated with precipitation amount west and east of 85°E, respectively (Fig. 3a). This spatial pattern is qualitatively similar to the spatial pattern of precipitation anomalies during the weak minus strong monsoon years (Fig. 2a). On the other hand, the spatial pattern of correlations between ARAB (%) and δ 18 O p at each grid cell are uniformly positive east of 64°E longitude ( Supplementary Fig. 10). These observations suggest that during weak monsoon years, a northwest-southeast trending band of anomalously enhanced zonal moisture flux exports 18 O enriched moisture from the ARAB source region across the Indian subcontinent and into Southeast Asia and Southern China (Fig. 2a and Supplementary Fig. 1). Over continental India (west of 85°E) time-series comparisons show that precipitation amount is negatively (r~−0.5) correlated with the ARAB contributions, δ 18 O p and δ 18 O pw (Fig. 3c) whereas the ARAB contributions are positively correlated (r~0.7) with δ 18 O p and δ 18 O pw (Supplementary Fig. 10a). These results are supported by a multivariate empirical orthogonal function (MV-EOF) 56 analysis using the seasonal mean of JJAS ARAB and precipitation fields (Supplementary Figure 10c-d). The MV-EOF analysis extends conventional EOF by extracting coupled patterns of variability among multiple fields.
The first mode of MV-EOF using the detrended and seasonal mean (normalized by its standard deviation) data illustrates a clear inverse relationship between the ARAB (%) and precipitation over the Indian subcontinent, explaining~26% of the shared covariance.
The relationships between moisture contributions from LND with precipitation amount, δ 18 O p and δ 18 O pw for the mid-to-late part of the monsoon season (AS = August to September) suggests the LND contributions during AS are positively (Fig. 3b) and negatively ( Supplementary Fig. 10b) correlated with precipitation amount and δ 18 O p and δ 18 O pw , respectively over the Indian subcontinent (west 85°E longitude). These observations suggest that higher precipitation amounts during the strong monsoon years lead to greater transpired moisture flux and vice versa. It is interesting to note that IsoGSM2 simulates strong negative correlations (~−0.7) between moisture contributions from LND with the δ 18 O p and δ 18 O pw over continental India, which act to strengthen the inverse relationship between the precipitation amount and δ 18 O p driven by moisture contributions from the ARAB. The negative relationship between LND and δ 18 O p (and δ 18 O pw ) likely arises due to the waning and waxing influence of moisture contributions from the ARAB and LND during the mid and late part of the monsoon season.
In sum, our analysis suggests that the interannual variability in δ 18 O p and δ 18 O pw over the ISM domain during the monsoon season is influenced by variations in the moisture contributions from the ARAB and terrestrial recycling. During weak monsoon years, the anomalous increase in zonal moisture flux via a northward shifted LLJ over the Arabian Sea and continental India serves as a conduit for transporting 18 O enriched moisture from the ARAB source region across a northwest-southeast trending region that extends from the Arabian Peninsula to the southeast Asia ( Supplementary Fig. 11). These model-based inferences are also supported by the positive correlations of JJAS zonal moisture flux anomalies with the observational monthly δ 18 O p data from the GNIP sites located within this band ( Supplementary Fig. 11). Together, the model-instrumental observations suggest that proxy records of δ 18 O p from sites situated within this corridor of anomalous moisture flux during weak monsoon years are   . Circles in all panels mark the locations of the cave records discussed in the text. The strong (1983, 1988, and 1994) and weak (1982, 1987, 2002, 2004, and 2009) monsoon seasons are defined as positive (negative) departures of one standard deviation from the mean of All India Rainfall 54 . The weak monsoon years are characterized by (i) negative precipitation anomalies over much of continental India except over northeast, northern BoB, and the coastal Myanmar regions; (ii) reduced vertically integrated zonal moisture flux along the climatological axis of the LLJ over the Arabian Sea, compensated by anomalously enhanced northwest-southeast trending zonal moisture flux across north India (with its core at~25°N) and near the equator; (iii) negative meridional moisture flux anomalies, particularly east of 75°E over the Indian subcontinent.  (% TPW) and precipitation amount at each grid cell. The correlation was performed after removing the climatological mean and trend from data (b) same as (a) but for LND (% TPW) vs. August to September (AS) precipitation amount. The black rectangle in each panel marks the region used for extracting time-series comparisons in panels (c) and (d). Stippling indicates regions of significant correlations at 95% significance level obtained after accounting for serial correlations and followed by the application of false discovery rate (FDR) 67 procedure with a 5% threshold (see "Methods"). FDR is the expected proportion of rejected hypotheses when the null hypothesis is actually true for those tests. particularly well-suited for reconstructing the ISM variability ( Fig. 1 and Supplementary Fig. 11).

Discussion
The results from our study have implications for interpreting the existing speleothem oxygen isotope records from the region ( Fig. 1 and Supplementary Note 1). The IsoGSM2 38 simulated JJAS δ 18 O p data extracted from the model's grid cells enclosing the speleothem sites in north and central India ( Supplementary  Fig. 1) exhibit a moderate inverse relationship with regional precipitation amount over the land and the oceanic regions upstream of the 85°E longitude (Fig. 4). Similarly, the time-series comparisons of the published speleothem δ 18 O profiles from Jhumar and Sahiya caves 26,57 with smoothed instrumental precipitation data (~1901-2006) also exhibit inverse relationships (Fig. 4), which are broadly consistent with the model simulated inverse δ 18 O p -precipitation relationship. However, the regional precipitation amount at both cave sites can only explain, at the best, up to half of the isotopic variance in δ 18  around Sahiya cave (yellow square) and model-simulated precipitation amount at all other grid points. The correlation was performed after removing the climatological mean and trend from the data. Stippling indicates regions of significant correlations at a 95% significance level obtained after accounting for serial correlations in data at each grid point followed by the application of FDR 67 procedure with a 5% threshold (see "Methods"). FDR is the expected proportion of rejected hypotheses when the null hypothesis is actually true for those tests. The black rectangle (17°N-32°N and 70°E-82°E source dynamics and circulation as well as from other processes such as organized convection over the ocean moisture source regions, rainout during the moisture transport, and local convection 13,14,17,49,58 . It is also important to note that the inverse δ 18 O p -precipitation relationships observed in the simulated and proxy data do not necessarily stem from the classic "amount effect" as previously suggested [23][24][25] but rather from the relative changes in moisture contributions and attendant changes in circulation. These observations imply that the speleothem δ 18 O records from north and central India can be broadly interpreted in a framework of overall changes in the upstream circulation, moisture source dynamics, and regional precipitation intensity. Reliable interpretation of speleothem data, however, must also take into account the various cave and karst-specific processes, which can modify the δ 18 O p signal to a varying degree (e.g., ref. 26 ).
While our analysis with IsoGSM2 tagging simulation data primarily elucidates the potential mechanisms of δ 18 O p interannual variability, it can provide potential insights to understanding the large amplitude and low-frequency changes observed in various ISM speleothem records on multidecadal-to-millennial-and-orbital timescales. For example, a large increase in δ 18 O (~6-8‰) during the millennial-scale cold climate events is a striking feature in some speleothem records from India 20,21 , which is not easily explained within the framework of modern climatological ISM dynamics. Observational data show that a large fraction of ISM intraseasonal variability is tied to quasiperiodic episodes of active and break spells-two distinct oscillatory modes of monsoon that have radically different synopticscale circulation and precipitation patterns 40,48 . A characteristic feature of the break phase is large-scale negative precipitation anomalies over much of continental India that are accompanied by anomalously positive zonal moisture transport with its core centered at~25°N at the expense of its climatological position at 15°N (Supplementary Fig. 1d)-features that are reminiscent of the weak minus strong years as discussed above. The difference of the break minus active composites of the amount weighted simulated δ 18 O p from IsoGSM2 38 show positive δ 18 O p anomalies of up to 6-8‰ over a wide swath of region extending from northwest India to southeast Asia that can be attributed to higher moisture contributions from the ARAB source region (Supplementary Fig. 12). Admittedly speculative, a scenario that calls for preferential locking of the ISM circulation into a "break mode" during the millennial cold events could potentially explain a large portion of the observed increase in δ 18 O values in the speleothem records 20,21 in conjunction with changes in the δ 18 O p resulting from other processes such as variations in the global ice volume, sea level, and regional temperatures. The same mechanism but working on shorter timescales (i.e., via the more frequent asymmetric occurrence of break spells in each monsoon season) 24 can also potentially explain the protracted episodes of anomalously enriched δ 18 O values reported from some speleothem records from India 23-27 . On orbital timescales, recent simulations with an Earth System model with water isotope tracers 36 indicate reduced ISM rainfall, heavier δ 18 O p , and increased moisture contributions from proximal sources (e.g., the Arabian Sea and the Red Sea) to the Indian subcontinent during precession maxima (low Northern Hemisphere summer insolation) and vice versa 36 thus, supporting the idea that δ 18 O p variations over the Indian subcontinent reflect large-scale changes in circulation, moisture sources, and rainfall amounts.

Methods
We used Isotope-incorporated Global Spectral Model version 2 (IsoGSM2) 38 , which is a version of the Scripps Experimental Climate Prediction Center's global spectral model. The atmosphere in IsoGSM2 has a horizontal resolution of~200 km (T62) with 28 vertical levels. The model is forced with prescribed sea surface temperatures and sea ice conditions from the optimal interpolation daily dataset provided by the National Centers for Environmental Prediction (NCEP) 59 . The default land surface scheme is the Noah model 60 . The IsoGSM2 applies a global downscaling technique, which nudges temperature and wind speed (at scales larger than 1000 km) toward the National Center for Environmental Prediction Reanalysis 2 (NCEP R2) 38,61 . The nudging (performed at every 6 h for all sigma levels) provides dynamical constraints to the simulated atmospheric circulation and is directly comparable with the observations. In IsoGSM, most isotopic fractionations are assumed to occur at thermodynamic equilibrium, except for open water evaporation, condensation in supersaturation conditions (vapor to ice), raindrop reevaporation, and air-rain isotopic exchange, where kinetic fractionation occurs. IsoGSM2 assumes no fractionation when water evapotranspires over land. All precipitation is fully mixed into a simple single bucket-type model for treatment of soil water isotopes, which provides storage of all water species, and consequently, the isotope ratio of evapotranspiration is assumed to be the same as stored values.
Simulated precipitation from IsoGSM2 broadly replicates the observations over South Asia although it underestimates precipitation amount over northwest India, Indo-Gangetic plains, and the Western Ghats and overestimates over the BoB and coastal Myanmar (Supplementary Fig. 2). IsoGSM2 simulated δ 18 O p has been previously validated with observational data from South Asia at event based 25 , daily 31,49 , and monthly timescales 26 . Here, we further validate the IsoGSM2 simulated data with some of the longest daily (Supplementary Table 1) and monthly GNIP δ 18 O p data from several additional sites located in south and southeast Asia. The model can reproduce the annual cycle and interannual variability in δ 18 O p with high fidelity for GNIP stations at Dhaka, New Delhi, Bangkok, and Kunming . In the tracer mode, IsoGSM2 can track a given moisture mass back to its last time at sea level by setting the surface evaporative fractionation factor to 1 for the target region (and 0 for all other regions) and by turning off all the other isotopic fractionation in the atmosphere. The evaporated water is marked by its origin with a tag. The tagged vapor once added to the atmosphere is allowed to undergo the same atmospheric hydrological cycle as normal water vapor until it leaves the atmosphere as precipitation 49,62,63 . The tagged simulation output (1979-2010) was generated at daily intervals and interpolated into 17 vertical pressure levels and monthly averaged.
We used the least-squares linear regression to determine linear Pearson correlation coefficients between independent and dependent variables. Confidence intervals of Pearson correlation coefficients were determined using a pairwise moving-block bootstrap method with an average block length proportional to the estimated data autocorrelation using a PearsonT3 software 64 . This method preserves the serial dependence of time series and provides the corrected 95% confidence intervals that are valid in the presence of autocorrelation 64,65 . The coverage accuracy is increased by applying calibration to standard error-based bootstrap Student's t confidence interval. The effective degree of freedom (EDOF) was estimated by the following equation where N denotes the length of time series and r1 and r2 are lag-one autocorrelations of each timeseries 66 . Unlike pairwise comparisons (where the significance level α-the probability of falsely rejecting the null hypothesis of zero correlation provides the basis of hypothesis testing), the statistical significance of field correlations (example, Figs. 3 and 4 in the main text and other supplementary figures) requires conducting multiple tests at different locations simultaneously. Increasing numbers of simultaneous tests (in accordance with increasing resolution of grid points) may lead to a higher number of correlations that may be incorrectly deemed as "significant"-the test multiplicity problem. We applied the false discovery rate (FDR) procedure 67 to control the proportion (q = 5%) of erroneously rejected null hypotheses using a MATLAB code 68 , where q guarantees that 5% or fewer of the locations where the null hypothesis is rejected are false detections. FDR generally provides a more powerful method for correcting for multiple comparisons (e.g., refs. 69,70 ) than procedures like Bonferroni correction that provide strong control of the familywise error rate (i.e., the probability that one or more null hypotheses are mistakenly rejected).

Data availability
The simulation data from isotope-incorporated Global Spectral Model version 2 (IsoGSM2) are available at http://isotope.iis.u-tokyo.ac.jp/~kei/tmp/isogsm2/. Moisture tagging simulation data are available at the Zenodo repository at https://doi.org/10.5061/ dryad.4b8gthtc7. Data are also available upon request from the corresponding authors. The All-India Rainfall (AIR) data from the Indian Institute of Tropical Meteorology (IITM) meteorological is available at https://mol.tropmet.res.in/. The precipitation isotope data by Global Network of Isotopes in Precipitation (GNIP) stations are available at https://nucleus.iaea.org/wiser. The Global Precipitation Climatology Centre (GPCC) version 7 precipitation data is available at https://psl.noaa.gov/data/gridded/data.gpcc. html. The Global Precipitation Climatology Project (GPCP) version 2.3 data are available at https://psl.noaa.gov/data/gridded/data.gpcp.html. The Climatic Research Unit Timeseries (CRU-TS) version 4.0 data are available at https://crudata.uea.ac.uk/cru/data/ hrg/. The moisture flux divergence, vertically integrated moisture flux, and precipitation