Research ArticleCLIMATOLOGY

Hydroclimate footprint of pan-Asian monsoon water isotope during the last deglaciation

See allHide authors and affiliations

Science Advances  22 Jan 2021:
Vol. 7, no. 4, eabe2611
DOI: 10.1126/sciadv.abe2611


Oxygen isotope speleothem records exhibit coherent variability over the pan-Asian summer monsoon (AM) region. The hydroclimatic representation of these oxygen isotope records for the AM, however, has remained poorly understood. Here, combining an isotope-enabled Earth system model in transient experiments with proxy records, we show that the widespread AM δ18Oc signal during the last deglaciation (20 to 11 thousand years ago) is accompanied by a continental-scale, coherent hydroclimate footprint, with spatially opposite signs in rainfall. This footprint is generated as a dynamically coherent response of the AM system primarily to meltwater forcing and secondarily to insolation forcing and is further reinforced by atmospheric teleconnection. Hence, widespread δ18Op depletion in the AM region is accompanied by a northward migration of the westerly jet and enhanced southwesterly monsoon wind, as well as increased rainfall from South Asia (India) to northern China but decreased rainfall in southeast China.


With nearly half of the world’s population affected by the pan-Asian summer monsoon (AM) (Fig. 1 and fig. S1), it is important to improve our understanding of the AM response to climate change. One important advance in our knowledge of past AM changes is the development of a large number of well-dated, high-resolution oxygen isotope speleothem (δ18Oc, for δ18O in the cave) records that, combined, span the past 640,000 years (1, 2). One notable feature of these δ18Oc records is their coherent temporal variability across the entire AM region (Fig. 2, C to E, and fig. S2). This AM δ18Oc signal during the last deglaciation [21 to 11 thousand years (ka)] can be seen as characterized by enrichment in Heinrich Stadial 1 (HS1; ~18 to 14.5 ka), depletion in the Bølling-Allerød interstadial (BA; 14.5 to 12.9 ka), and further enrichment in the Younger Dryas stadial (YD; 12.9 to 11.7 ka).

Fig. 1 iTRACE hydroclimate-climatology at the LGM.

δ18O (A) Model summer [June, July, August (JJA)] precipitation (Prcp, shading), 850-hPa wind (vector), and surface MSE (contours, in K, regions higher than 315 K plotted); (B) δ18Op (shading) and zonal wind velocity averaged between 500 and 200 hPa (contour, solid for positive and dash for negative) over the AM region at LGM. In (A), red boxes mark the regions of northern China (N. China), southern China (S. China), and South Asia (S. Asia) in Fig. 2 (C to E), respectively, while purple crosses and red stars are proxy sites in Fig. 2 (C to E). This figure also serves to illustrate the gross features of the AM climatology for the present day and LGM, which are similar.

Fig. 2 Hydroclimate-δ18O evolution in iTRACE and observation.

(A) Forcing: June insolation at 60°N (red), atmospheric CO2 concentration (green) (17), and meltwater fluxes in the Northern (blue) and Southern (orange) Hemisphere; (B) Greenland (GISP2): Observed ice-core δ18Op (black) (38) and model annual δ18Op (purple), surface air temperature (red), and precipitation (blue); (C) N. China: Calcite δ18Oc in model (purple) and observation (black) (39) at Kulishu Cave, annual precipitation in northern China in model (37°N to 43°N, 108°E to 120°E; blue) and in Gonghai Lake in reconstruction (green) (20), as well as lake level from Dali Lake (gray dots) (21); (D) S. China: δ18Oc in model (purple) and observation (black) (1) at Hulu Cave; model annual precipitation in southern China (22°N to 35°N, 108°E to 120°E; blue) and the leading PC of Mg/Ca, Sr/Ca, and Ba/Ca records observed in Haozhou, a proxy for southern China monsoonal precipitation (green) (22); (E) S. Asia: δ18Oc in model (purple) and observation (black) (40) in Mawmluh Cave, model annual precipitation in South Asia (10°N to 25°N, 70°E to 95°E; blue), model Arabian Sea upwelling (23°N to 24.3°N, 62.8°E to 65°E; gray), and observed Arabian Sea sediment total reflectance as a proxy for SAM wind strength (green) (41); (F) Model AMOC intensity (blue) and observed 231Pa/230Th in sediment core GGC5 as a proxy for AMOC intensity (black) (42). All model results have been smoothed with a 50-year running mean. Note that in each panel, the δ18Oc scales are the same in the model and observation but with an offset (Materials and Methods). In (E), the Arabian Sea upwelling (20- to 50-m average) is constructed by removing upwelling responses in ICE run in which the local bathymetry changes at 14 and 12 ka produce artificial jumps. In (C) to (E), monsoon regions are marked as red boxes and the observational sites as stars and crosses in Fig. 1A. ka BP, thousand years before the present; masl, height above sea level; NHSI, Northern Hemisphere Solar Insolation; PDB, Pee Dee belemnite; ppm, parts per million; SMOW, standard mean ocean water; TS, surface temperature.

This pioneering work in developing AM δ18Oc records has provided important insights into changes in the AM associated with the effects of millennial-scale North Atlantic climate changes and orbital-scale insolation changes (1, 3, 4). In particular, this earlier work proposed two mechanisms that link the observed δ18Oc variability in the East Asian Monsoon (EAM) to these forcings: changing proportions of low δ18O summer-monsoon rainfall in annual totals (referred to as monsoon intensity) (1, 4) and changing fractions of rainout and precipitation δ18O from air masses advected from tropical Indo-Pacific source regions to EAM cave sites (3). The proposed remote-depletion effect on δ18Oc variability over the EAM region has been confirmed by numerical experiments (58). In contrast, the spatial/temporal features of monsoon rainfall that may also contribute to the δ18Oc signal and their underlying physical mechanisms remain unclear (712), with some studies attributing the δ18Op signal to remote depletion with no local rainfall change (6, 8). In fact, if the δ18Oc signal in the EAM region is primarily from the remote-depletion effect, the next fundamental question is: what are the local rainfall changes over the EAM, and more generally, the AM region?

Here we address the question of whether there is a continental-scale hydroclimate response pattern, or footprint, associated with this continental-scale AM δ18Oc signal and if so, what is the mechanism that maintains such a continental-scale coherence of the hydroclimate variability, especially in rainfall. Based on transient experiments in an isotope-enabled Earth system model and their comparison with proxy records, we propose that the widespread AM δ18Oc signal during the last deglaciation is accompanied by a continental-scale, coherent hydroclimate footprint with spatially opposite signs in rainfall. This coherent pattern is caused by atmospheric teleconnections as well as global-scale external climate forcing.


We conducted the first set of transient simulations of the evolution of global climate and water isotopes (δ18O, δD) during the last deglaciation in a state-of-the-science isotope-enabled [Community Earth System Model (iCESM); Supplementary Text] (13), largely following the modeling strategy of TRACE-21ka (Supplementary Text) (14). The iCESM captures the major observed features of the present climatology in climate and water isotopes in the AM region, including monsoon winds, precipitation, and δ18O of precipitation (δ18Op) (Supplementary Text) (Fig. 1 and fig. S1). Our deglacial simulations are forced by four forcing factors that are applied additively (Supplementary Text). Starting from the Last Glacial Maximum (LGM; 20 ka), the first experiment is integrated with changing ice sheets and ocean bathymetry (ICE) (15). We next added insolation forcing (ICE + ORB) (16) followed by greenhouse gases (ICE + ORB + GHG) (17) and, lastly, meltwater fluxes (ICE + ORB + GHG + MWF) (Fig. 2A) or the isotope-enabled transient climate experiment (iTRACE). Insofar as the forced responses are roughly linear, these transient sensitivity experiments allow us to isolate each forcing effect, approximately, as follows: the ice sheet and bathymetric effect (ICE), the orbital effect (ICE + ORB − ICE), the GHG effect (ICE + ORB + GHG − ICE + ORB), and the meltwater effect (iTRACE − ICE + ORB + GHG). Our transient simulations provide the most comprehensive model-data comparison study to date on AM δ18O variability and the associated climate changes during the last deglaciation, allowing us to directly assess the magnitude and timing of the full deglacial variability observed in the high-fidelity δ18Oc records from the AM region.

Model-data comparison of δ18O-hydroclimate evolution

We first describe our model simulation over Greenland because of the clear linkage between ice-core δ18Op (Fig. 2B) and AM δ18Oc (Fig. 2, C to E) (1, 2). Model Greenland δ18Op evolution closely resembles the GISP2 (Greenland Ice Sheet Project 2) δ18Op record (Fig. 2B). Largely following the annual temperature in positive correlation (Fig. 2B) as expected in the temperature effect, model δ18Op first decreases into the HS1, then increases by ~4 per mil (‰) for a warming of ~11°C during BA, and, at the end of YD, giving a temporal δ18Op-temperature slope of ~0.36‰/°C, consistent with independent calibrations (18). Our sensitivity experiments show that this δ18Op-temperature coevolution is caused mainly by increasing GHG and meltwater forcing (fig. S3A) with an additional impact from sea ice. In particular, increased meltwater fluxes during HS1 and YD cause the Atlantic Meridional Overturning Circulation (AMOC) to collapse (Fig. 2F and fig. S3B), with attendant cooling over the North Atlantic, suppressing the warming trend from the GHG forcing (19).

Our modeled AM δ18Oc variability also captures the coherent changes in reconstructed δ18Oc variability across the entire AM region during the last deglaciation (Fig. 2, C and D, fig. S2, and Supplementary Text) and the roughly out-of-phase relationship with Greenland δ18Op (Fig. 2B). In particular, our model successfully simulates the magnitude of the observed δ18Oc variability across the AM region (Fig. 2, C to E, and fig. S2), in contrast to many previous simulations that tend to underestimate the magnitude of the δ18Oc response (58). A decomposition of the AM δ18Oc response to those from different seasons, and furthermore from the δ18Op value and the monthly precipitation weight, shows a dominant contribution by summer monsoon in the δ18Op value (fig. S4). We will thus focus herein on the summer monsoon and the associated δ18Op. To compare the observed AM δ18Oc records directly, we calculated model δ18Oc from annual temperature and δ18Op (Supplementary Text).

Our model AM δ18Oc variability evolves roughly out of phase with Greenland δ18Op (Fig. 2B). Over most sites, the δ18Oc variability also resembles closely, but out of phase with, the local temperature change (fig. S2, A to G, I, and K). Furthermore, across the AM region, the AM δ18Oc variability is contributed primarily by the variability of precipitation δ18Op, instead of cave temperature correction (fig. S2). Therefore, the AM δ18Oc variability is opposite to what would be expected from the temperature effect, leaving no obvious reason that the AM δ18Oc variability is caused by the temperature effect locally. Instead, temperature varies coherently in phase across the entire Northern Hemisphere, as seen in the similar temperature evolution between across AM cave sites and Greenland (fig. S2 versus Fig. 2B). Thus, this AM δ18Oc response could be caused by the hydrological variability along with the hemisphere-wide temperature change associated with, mainly, AMOC changes, as proposed in previous studies (58). Regionally, the AM δ18Op signal is caused predominantly by the amount effect in the South Asian monsoon (SAM; mainly India) region and the subsequent en route depletion in the EAM region; furthermore, our analysis shows that the Indian Ocean source overwhelms the Pacific and acts as the dominant moisture source for deglacial variability (Fig. 3 and figs. S5 to S7). The detailed mechanism of the AM δ18Op will be described in more details in Materials and Methods.

Fig. 3 Tracking δ18Op and its response at HS1.

Climatological annual δ18Op map at LGM (20 to 19 ka) originating from the source region (A) Indian Ocean, (B) Pacific Ocean, (C) East Asia, (D) South Asia, and (E) other regions in the LGM tagging experiment. (F to J) Same as (A) to (E) but for the difference between HS1 (15.8 to 15.5 ka) and LGM. Gray contours mark the tagging source regions of (A and B) oceans and (C and D) land. Also shown are the contributions from different sources to the annual climatological δ18Op in (K) N. China (37°N to 43°N, 108°E to 120°E), (L) S. China (22°N to 35°N, 108°E to 120°E), and (M) S. Asia (10°N to 25°N, 70°E to 95°E) at LGM (blue) and HS1 (red), with the three rainfall regions marked in Fig. 1A.

Hereafter, we focus on the continental-scale hydroclimate variability that accompanies this AM δ18Op variability. As will be addressed more fully later, we note that this hydroclimate variability is not necessarily the cause of local δ18Op variability everywhere in the AM region. In contrast to the δ18Oc signal that is broadly the same sign across the AM region, model rainfall response and, in turn, the rainfall-δ18Oc relation change sign within the AM region. Notably, model rainfall over South Asia (Fig. 2E and fig. S2, H to K) and nearby region (fig. S2, L and M) correlates negatively with δ18Oc as in the amount effect. Model rainfall over northern China is dominated by a wetting trend, punctuated by drying episodes during HS1 and YD (Fig. 2C and fig. S2, A and B), bearing some resemblance to South Asia (Supplementary Text), also with a somewhat negative rainfall-δ18Oc correlation. However, model rainfall in southern (more precisely, southeast) China evolves almost out of phase with that in South Asia and, to some extent, in northern China, despite having similar δ18Oc variability, leading to a positive rainfall-δ18Oc correlation (Fig. 2D).

Since the δ18Oc response is widespread of the same sign, the different model rainfall-δ18Oc relations across the AM region are caused by different rainfall variability, which appears to be consistent with other hydroclimate proxies independent of δ18Oc. Model rainfall over South Asia is consistent with the Arabian Sea upwelling (Fig. 2E), which is driven by monsoon winds via coastal upwelling. Over the northern China area of East Asia, the wetting trend of model rainfall largely agrees with lake status records (7), while the millennial variability during BA and YD agrees with high-resolution hydrological reconstructions (Fig. 2C) (20, 21). Last, our model rainfall in southern China is largely consistent with a recent multiproxy reconstruction of moisture variability (Fig. 2D) (22).

The difference in regional evolution behavior between δ18Oc and rainfall can also be seen in their responses to different forcing in our sensitivity experiments. The AM δ18Oc variability is caused primarily by meltwater forcing via the AMOC and secondarily by orbital forcing, of the same sign across the AM region (Fig. 4). The rainfall variability over South Asia and southern China is also driven primarily by meltwater and secondarily by insolation but of the opposite signs. In contrast, rainfall variability in northern China is driven primarily by summer insolation and secondarily by meltwater (Fig. 4 and fig. S3, C to H).

Fig. 4 Hydroclimate and δ18Oc variability over Asia monsoon region in four deglacial sensitivity experiments.

(A) N. China: Cave δ18Oc (top) and annual precipitation (Prcp; bottom). (B) Same as (A) but for S. China. (C) Same as (A) but for S. Asia. Model results have been smoothed with a 50-year running mean for ICE (blue), ICE + ORB (red), ICE + ORB + GHG (orange), and ICE + ORB + GHG + MWF (iTRACE) (black). Regions are the same as in Fig. 2 (C to E).


The hydroclimate footprint

Although the rainfall-δ18Oc relation across the AM region appears complex, this relation can be understood in terms of a continental-scale hydroclimate response, which has a coherent and distinctive footprint that coevolves with the widespread δ18Oc variability (Supplementary Text). We identify this footprint as the leading maximum covariance analysis (MCA) mode (Supplementary Text) between summer δ18Op and rainfall fields (explaining 93% of the covariance) (Fig. 5). The time coefficient of this MCA mode is characterized by an increasing trend punctuated by negative excursions in HS1 and YD (Fig. 5D), similar to the δ18Oc variability (Fig. 2, C to E). This footprint appears robust because it can also be extracted as the leading combined empirical orthogonal function mode of the combined δ18Op and summer rainfall field or additionally with 200-hPa zonal wind and 850-hPa meridional wind (explaining ~65% of total variance).

Fig. 5 Asian summer monsoon hydroclimate footprint and the accompanying δ18Op, as illustrated by the first MCA mode between summer (JJA) δ18Op and precipitation over the AM region in iTRACE.

(A) Summer δ18Op (shading) and 200-hPa wind (vector) regressed on the normalized time coefficient of δ18Op [black in (D)]. (B) As in (A) but for summer precipitation (shading) and 850-hPa wind (vector). (C) As in (A) but for zonal wind u (shading), meridional wind v, and vertical velocity w (vector) zonally averaged over East Asia (108°E to 120°E). (D) Time series of the normalized time coefficients of the first MCA mode (black for δ18Op, green for precipitation, both using the same black axis), 850-hPa meridional wind v strength over East Asia (22°N to 43°N, 108°E to 120°E; yellow), 200-hPa westerly jet latitude centroid (red) (108°E to 120°E), and Indian Ocean (0° to 20°N, 40°E to 100°E) summer precipitation (Prcp, blue). Regression fields at 99% (P < 0.01) confidence level in two-tailed Student’s t test are plotted. In (A) and (B), wind vectors with the strength smaller than 0.2 ms−1 are masked. In (C), contours are for climatological u (positive, solid red; negative, dash blue) at LGM, and w is amplified by a factor of 500 for visualization. In (D), the time coefficients of δ18Op and precipitation are almost identical, further highlighting the close coupling of the two fields; the jet latitude centroid is calculated as the latitude weighted by the strength of the westerly wind (u > 15 ms−1). Note that the response pattern of (A) to (C) corresponds to the periods of positive coefficient, such as BA. For the period of negative time coefficient, such as HS1, the response pattern is of the opposite sign.

Further analysis of this MCA mode shows a widespread δ18Op signal (Fig. 5A) accompanied by a continental-scale hydroclimate response in rainfall, monsoon winds (Fig. 5B), and the westerly jet (Fig. 5A). Overall, a negative δ18Op loading over the AM and Indian Ocean regions accompanies a positive rainfall loading (Fig. 5, A and B) such that the rainfall-δ18O correlation appears negative as expected in the amount effect. Further examination, however, shows that the South Asia rain belt extends eastward only into northern China, leaving a “drought hole” over southern China (Fig. 5B). The opposite rainfall changes between northern and southern China thus lead to a rainfall-δ18O correlation that changes from negative in northern China (and South Asia) to positive in southern China, as identified in the time series (Fig. 2, C and D, and fig. S2, A to G).

The rainfall increase in South Asia/northern China can be understood, dynamically, as being associated with intensified southwesterly monsoon wind and resulting moisture convergence (Fig. 5B), while the decreased rainfall in southern China is associated with moisture divergence (7). This monsoon wind and the rain belt over East Asia are also coupled with the westerly jet in the upper troposphere (23, 24). The EAM is located just south of the entrance region of the western Pacific westerly jet (Fig. 1B), where the ascent associated with the secondary circulation of ageostrophic wind provides a favorable environment for monsoon rainfall (23, 25). The westerly jet also steers transient eddies along the Meiyu front, generating synoptic conditions conducive to convective instability and ascent for EAM (26). Latent heating further induces ascent and then low-level southerly wind through the Sverdrup relation by stretching the air column (27). The heated atmospheric column also increases the temperature gradient to the north, shifting the westerly jet northward via thermal wind (9). The close coupling among δ18Op, rainfall, and winds in East Asia can also be seen in the almost identical evolution of the MCA coefficients, the latitude of the local westerly jet, and the intensity of southerly monsoon wind (Fig. 5D). The coupled response of monsoon rainfall and winds thus forms a distinctive AM hydroclimate footprint accompanying the δ18Op variability, such that δ18Op depletion is associated with increased rainfall from South Asia to northern China but decreased rainfall in southern China, along with a northward migration of the westerly jet (9, 22) and enhanced southwesterly monsoon wind (Fig. 5, A to C) (7).

Mechanism for the coherence of the hydroclimate footprint

Next, from climate dynamics perspective, we further identify the mechanisms that explain why the AM hydroclimate footprint exhibits such a continental-scale coherent pattern, while, at the same time, exhibiting spatially opposite signs in rainfall (Fig. 5) using sensitivity experiments. First, this continental-scale footprint is synchronized across regions by the local response to global-scale forcing, notably the meltwater and insolation, at millennial and orbital scales (fig. S8). This is in contrast to interannual rainfall variability that has a small decorrelation distance of ~500 km because it is induced predominantly by atmospheric internal variability (28). At the orbital time scale, increased summer heating increases surface moist static energy (MSE) north of the anomalous rain belt from South Asia to northern China (fig. S8A3), extending the region of maximum surface MSE northward (Fig. 1A) and thus the AM monsoon rain belt northward (29, 30). This increases monsoon rainfall from South Asia to northern China along with a southwesterly monsoon wind and a decreasing rainfall in southern China (fig. S8, A3 and B3) (7). Heating of the atmospheric column in the AM region also shifts the westerly jet northward via the thermal wind, further enhancing the southerly monsoon wind and shifting the rain belt northward over East Asia (fig. S8C3) (9).

At the millennial time scale, the weakening AMOC by increased meltwater during HS1 and the YD results in opposite responses to those induced by increased insolation. Specifically, the reduced AMOC cools the atmosphere and decreases surface MSE north of the rain belt from South Asia to northern China, shifting the westerly jet southward over East Asia (fig. S8, A1, B1, and C1; but with a negative sign because HS1 and YD are negative excursions in the correlation time series in Fig. 5D). This, in turn, decreases southerly monsoon winds and rainfall in the rain belt extending from South Asia to northern China, while increasing rainfall in southern China (again, negative sign of fig. S8, A1 to C2) (31). This EAM response can also be understood from the MSE budget of the atmospheric column as the result of a southward shift of the westerly jet and its horizontal advection of MSE (fig. S9 and Supplementary Text).

The synchronization between climate changes in SAM and EAM is further enhanced by atmospheric teleconnection. In summer, anomalous latent heat released by deep convection over the Indian Ocean/SAM region generates disturbances in the upper tropopause that propagate into East Asia via the so-called Silk Road teleconnection (32), enhancing rainfall over northern China (33), as in present-day interannual and interdecadal climate variability (34). This teleconnection response can be also shown in our sensitivity experiments. As an example, a diabatic heating is imposed over the northern Indian Ocean/SAM region (fig. S10A). The heated atmosphere column generates an anomalous anticyclone in the upper troposphere, which then propagates downstream in the westerly waveguide, forming an anticyclone over northeastern Asia and shifting the local westerly jet northward (fig. S10, A and C). The shifted jet allows for a northward penetration of the southerly monsoon wind via the secondary circulation and transient eddies, increasing rainfall over northern China, while decreasing rainfall in southern China (fig. S10B). This complex pattern of heating-induced rainfall response is associated with a widespread δ18Op depletion over the AM region (fig. S10D), largely consistent with the coupled δ18Op-rainfall pattern of the footprint (Fig. 5, A and B).

We present the first transient simulation of the coevolution of water isotopes and hydroclimate over the globe during the last deglaciation in a state-of-the-science isotope-enabled Earth system model. Our model is able to quantitatively reproduce the observed variability directly in the high-fidelity δ18O records across the Asian monsoon region from speleothems and over Greenland from ice cores. Different from previous work, which tended to focus on the formation mechanism of the δ18O signal, our focus here is the hydroclimate footprint that accompanies the δ18O signal. Our model-data study indicates that the widespread AM δ18Op variability is accompanied by a continental-scale hydroclimate footprint that is generated coherently by large-scale climate dynamics. Meanwhile, δ18Op is changed by rainfall variability in the Indian Ocean and SAM region and is extended downstream into the EAM region. This leads to local rainfall coevolving with δ18Op over the AM region, directly via isotopical processes (the amount effect) in the Indian Ocean and SAM region but indirectly via large-scale climate dynamics downstream in the EAM region.

The mechanism for the hydroclimate footprint and its relation with the δ18Op signal is illustrated schematically in Fig. 6. During increased summer insolation or enhanced AMOC, such as in the case of the BA relative to HS1, the warming over the AM region tends to intensify deep convection and rainfall over the Indian Ocean and SAM region. With the moisture source contributed primarily from the Indian Ocean, this rainfall response generates a widespread δ18Op depletion locally in the SAM region and then remotely in the EAM region. Meanwhile, over the EAM region, the westerly jet is forced to shift northward associated with increased moisture transport in monsoon winds, causing moisture convergence and increased rainfall in northern China but moisture divergence and decreased rainfall in southeast China. The overall climate response therefore forms a continental-scale rainfall response that is characterized by a rain belt extending from South Asia to northern China and a drought hole in southeast China. This rainfall response is further reinforced by the atmospheric teleconnection that propagates in the westerly waveguide from the SAM region into the EAM region. Hence, the AM δ18O signal appears negatively correlated with the local rainfall over the SAM region and northern China but positively correlated with the local rainfall over southeast China. Accompanying this hydroclimate footprint is a hemisphere-wide temperature response to the climate forcing, either AMOC or orbital forcing, which leads to a similar temperature response between the AM region and Greenland. This temperature response leads to an apparent negative correlation between cave δ18O and local temperature as noted previously in fig. S2.

Fig. 6 Schematic figure for the mechanism of the AM δ18Op and hydroclimate footprint.

The coupled δ18Op-hydroclimate response is forced by a warming response over the AM region to either enhanced summer insolation or strengthened AMOC (orange arrows). This response causes widespread δ18Op depletion (gray shading on 850 hPa) by enhanced deep convection (purple cloud) and rainfall (green shading on 850 hPa) locally over the Indian Ocean and SAM region but by downstream depletion along the route of moisture transport (blue arrows on 850 hPa) from the major source of India and secondarily western North Pacific over the EAM region. The warming also shifts the westerly jet northward and enhances the northward moisture transport over the EAM, increasing rainfall in northern China but decreasing rainfall in southern China (green and dark red shading on 850 hPa). The enhanced convective heating over the Indian Ocean and SAM further induces an atmospheric wave train (red and cyan circles on 200 hPa) propagating into East Asia in the westerly jet waveguide in the upper troposphere (gray arrow on 200 hPa), further synchronizing the rainfall between the SAM region and EAM region. (The signs of the responses here are consistent largely with the case of BA when AMOC recovers and insolation is high and should be reversed for the meltwater impact in HS1 and YD.)

The distinct hydroclimate footprint of the AM δ18Oc provides a dynamically consistent framework to understand the climate implication of the AM δ18Oc records. It enables these δ18Oc records to provide a first-order picture of the overall hydroclimate feature associated with the δ18Oc variability. In the meantime, the strong heterogeneity of the hydroclimate footprint highlights the distinct nature of each regional monsoon subsystem that may not be directly related to the local δ18Oc change, with implication to local monsoon climate changes not only in the past but also in the future.


Mechanism for water isotope response

The overall homogeneous response of δ18Op is caused mainly by its transport effect originating from the Indian Ocean, as suggested by previous studies (6). This is confirmed by our tagging experiments, which enable us to quantitatively assess the impacts on AM δ18Op variability from the Pacific versus Indian Ocean sources, precipitation weight versus δ18Op value, and en route depletion versus evaporation/condensation.

Tagging experiments. Two water tagging experiments are carried out at the LGM (20 ka) and late HS1 (15.5 ka), respectively, in the atmospheric component model iCAM5.3. The LGM and HS1 tagging experiments are forced by the sea ice distribution, sea surface temperature, and sea surface δ18O and δD extracted from the iTRACE experiment as well as the continental ice sheets, orbital parameters, and GHG concentrations at 20 and 15.5 ka, respectively. Each tagging experiment is integrated for 40 years with the last 20 years used for analysis. The tagging experiments reproduce the HS1-LGM pan-Asian δ18O response in the coupled iTRACE simulation (Fig. 3), as expected.

In the tagging experiments, the life cycle of H216O (vapor) and H218O are tracked from each tagging region (source region) where they evaporate and then follow the hydrological processes in the model to the region where they rain out (sink region). We have divided the global source regions into 25 subregions, with 13 covering the ocean and the rest covering land. Here, for the purpose of studying the AM, we group the 25 subregions into five regions: the Indian Ocean, the Pacific Ocean, East Asia (land, 105°E to 140°E, 10°N to 45°N), South Asia (land, 60°E to 105°E, 20°S to 25°N) (regions marked in Fig. 3, A to D, or fig. S5, A to D), and the rest of the globe including both ocean and land that is denoted as the “others” region, numbered as regions i = 1, 2, 3, 4, and 5, respectively. Thus, the precipitation P (=H216O) and δ18Op at a grid point are the sum from all the five source regions (35)P=Σi=15Pi,δ18Op=Σi=15δ18Oi(PiP)where Pi and δ18Oi represent the precipitation and δ18Op from the tagging region “i.

Decomposition of δ18Op response. The response of δ18Op during HS1 or the difference between HS1 and the LGM, Δδ18Op = δ18Op, HS1 − δ18Op, LGM, can be decomposed in two steps. In the first step, Δδ18Op is decomposed into two parts: the change in the value of the isotope ratio Δδ18Op and the change associated with the precipitation weight Δ(PiP) (35)Δδ18Op=Σi=15(PiP)Δδ18Opi+Σi=15δ18OpiΔ(PiP)(1)

The part of (PiP)Δδ18Opi for each region and their sum are shown in the fourth row, while the other part δ18OpiΔ(PiP) is shown in the bottom row of fig. S6. Here, we point out an important constraint on this decomposition that, to our knowledge, has not been noticed in previous publications. Since the sum of the total precipitation weights is always 1, the sum of the changes of all the precipitation weights is always 0, i.e., Σi=15Δ(PiP)=0. This simply states that, if there are increases of precipitation weights in some regions, the increases must be compensated by the decreases from some other regions. This is indeed the case for the compensation effects between the Pacific and Indian Ocean sources, as will be discussed later.

In the second step, the δ18Op from a source region i can be decomposed following its life cycle along the trajectory from the source region to the sinking region as a sum of three processes, approximately (35)δ18Opi=(δ18Ov,source,i)+(δ18Ov,sink,iδ18Ov,source,i)+(δ18Opiδ18Ov,sink,i)

Here, on the right-hand side, the first term denotes the vapor δ18O in the source region, which is produced by the net effect of local precipitation and evaporation at the source region. The second term denotes the en route depletion of vapor δ18O due to rainout along its trajectory from the source region to the sink region. The third term denotes the local condensation enrichment due to the transition from vapor to rainfall in the sink region, along with any post-condensation processes such as rain evaporation. Therefore, the change of the source region δ18Opi value (Δδ18Opi) can also be decomposed into three processes: the changes in the source region, the en route depletion, and condensation enrichmentΔδ18Opi=Δ(δ18Ov,source,i)+Δ(δ18Ov,sink,iδ18Ov,source,i)+Δ(δ18Opiδ18Ov,sink,i)(2)

The three processes are displayed in fig. S6 in rows 1, 2, and 3, respectively, for the five source regions and their sum.

Sources for the AM δ18Op response. The tagging experiments show that the climatological δ18Op over the AM region originates mainly from the Indian Ocean and Pacific Ocean (Fig. 3, A to E and K to M), largely consistent with the moisture tracking (fig. S5, A to E and K to M) and moisture budget there. In climatological δ18Op (as the average of LGM and HS1), the combined contribution of the Indian Ocean and Pacific Ocean to SAM, southern China, and northern China are, respectively, 64, 64, and 40%. From the LGM to HS1, δ18Op is enriched across the AM with the major source contribution similar to the climatology (Fig. 3, K to M): the Indian Ocean for the entire AM region (Fig. 3F) and, additionally, the Pacific Ocean for East Asia (Fig. 3G), with land recycling further enriching northern China (Fig. 3H). The combined Indian Ocean and Pacific Ocean contribution to the change of δ18Op to SAM, southern China, and northern China is, respectively, 97, 54, and 25%. The enrichment over South Asia (Fig. 3, F and M) is caused by the rainfall reduction from the dominant moisture source of the nearby Indian Ocean (fig. S5, F and M), consistent with the amount effect; the enrichment over East Asia (Fig. 3, K and L), however, is accompanied by increased local rainfall from all moisture sources (fig. S5, K and L), opposite to the expectation from the amount effect (also see discussion on time series in Fig. 2, C to E).

In northern China, it is expected that the enrichment is also contributed significantly by land recycling over the EAM region (Fig. 3H). However, the recycling moisture originates mainly from the oceanic sources, and therefore, the contribution on the δ18Op evolution is also similar to those from the oceanic sources. The seasonal cycle of δ18Op peaks in summer and therefore appears to follow the temperature effect in the model and observation (36). The deglacial evolution of δ18Op also appears, at the first sight, to follow closely that of annual temperature, but it is in the opposite sign to the temperature effect (fig. S2, A and B) and therefore is not caused by the temperature effect. This shows that the δ18Op-temperature relation can differ significantly for different time scales (36).

The widespread δ18Op enrichment is caused predominantly by the enrichment effect from the Indian Ocean δ18Op value, through the en route depletion, consistent with (6). This conclusion is derived below in two steps. First, we use decomposition Eq. 1 to show that the Pacific source effect is compensated by the Indian Ocean source due to the change of precipitation weight, an effect that seems to have not been discussed in previous works. Second, we use decomposition Eq. 2 to show that the remaining effect from the Indian Ocean is dominated by the en route depletion.

Pacific effect compensated by Indian Ocean in precipitation weight. The dominant effect from the Indian Ocean source is actually not obvious at the first sight. The net enrichment effect from the Pacific source is comparable with that from the Indian Ocean source, especially for EAM (Fig. 3, F, G, and K to M). However, a decomposition of the contribution into the changes due to precipitation weight and δ18Op value following Eq. 1 (rows 4 and 5, respectively, in fig .S6) shows that the enrichment from the Pacific (Fig. 3G) is caused mainly by the contribution of reduced precipitation weight δ18Op1Δ(P1P)>0 (due to Δ(P1P)<0 and δ18Op1 < 0; fig. S6B5), which is nevertheless compensated by the depletion effect from the Indian Ocean due to increased precipitation weight δ18Op2Δ(P2P)<0 (due to Δ(P2P)>0and δ18Op2 < 0; fig. S6A5), leaving little net effect associated with the weight change: δ18Op1Δ(P1P)+δ18Op2Δ(P2P)0 (fig. S6F5). The reduced precipitation weight from the Pacific occurs despite an increased precipitation amount from the Pacific source (ΔP1 > 0) because this increase from the Pacific moisture source is smaller than that from total sources and, in turn, the Indian Ocean moisture source (ΔP2 > ΔP > ΔP1 > 0; fig. S5, K and L and F and G), leading to a reduced weight from the Pacific Δ(P1P)<0and increased weight from the Indian Ocean Δ(P2P)>0.

En route depletion from the Indian Ocean. The δ18Op value enrichment can be further decomposed following Eq. 2 into three processes as (i) changes in the source, (ii) en route depletion, and (iii) condensation enrichment at the sink region (fig. S6, rows 1 to 3). This decomposition shows that the dominant contribution to δ18Op enrichment is the reduced en route depletion from the Indian Ocean (fig. S6A2), caused by decreased rainfall along the moisture pathway mainly over the Indian Ocean and SAM regions (Fig. 3F and fig. S5F). This dominant role of Indian Ocean source rainfall depletion is consistent with the strong resemblance between the evolution of Indian Ocean precipitation and the footprint coefficients, jet migration, and EAM winds (Fig. 5D).

The en route depletion from the Indian Ocean can be further assessed quantitatively in the combined hydroclimate-δ18Ov evolution in iTRACE. Figure S7 shows the evolution of regional δ18Ov and specific humidity q roughly along the mean moisture pathway from the Indian Ocean into East Asia, from 80°E to 110°E along 10°N (fig. S7A) and then from 20°N to 40°N along 115°E (fig. S7B). This evolution shows interesting features in two aspects: the spatial pattern (“spatial slope”) at a given time and the temporal evolution (“temporal slope”) of this pattern. First, the spatial pattern first follows the “super-Rayleigh depletion” in the SAM and Indian Ocean regions (follow longitude along 10°N) (fig. S7A) and then follows the Rayleigh depletion in East Asia (with latitude along 115°E) (fig. S7B). The super-Rayleigh depletion is steeper than the Rayleigh depletion because of deep convection and rain evaporation effects, reflecting the dominant deep convection process in the source region in the tropics, including the nearby SAM region (37). The Rayleigh process in East Asia further supports the en route depletion mechanism there, as identified in the tagging experiments (fig. S6). These two depletion processes also explain why the amount effect is dominant in South Asia while remote, instead of local, depletion by the amount effect becomes dominant for East Asia. Second, this spatial pattern tends to evolve temporally following the source region of the Indian Ocean, experiencing an overall trend of warming (Fig. 2 and fig. S2) and moisture increase with time (from left to right in fig. S7), with a δ18Ov enrichment into HS1, depletion into the BA, and then modest enrichment and depletion again through the YD, consistent with the time series in Fig. 2 (C to E) and fig. S2. The temporal slope, for both enrichment and depletion, also exhibits a super-Rayleigh depletion that is steeper than just Rayleigh, suggesting that the temporal slope of the entire AM region tends to follow the source region, again consistent with the dominant source from the Indian Ocean.

Last, the AM δ18Op response seems not to be caused by a major change in the precipitation seasonal cycle. This follows because the seasonal cycle of the coupled δ18Op precipitation remains qualitatively similar at the LGM, HS1, and P1(pre-industry) in the AM region (fig. S1, B to D). In addition, the δ18Op responses associated with the change of precipitation seasonal weight are not the major contribution to total δ18Op in both South and East Asia (fig. S4).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Our simulations were performed at the NCAR supercomputing facility. Funding: This work is supported by Chinese NSFC 41630527; U.S. NSF 1810682, 1810681, 1401803, and 1656907; Chinese NSFC 41888101; MOST2017YFA0603801; Nanjing Univ. Inf. Sci. & Tech. and Nanjing Normal University. The CESM project is supported primarily by the NSF. This material is based on work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the NSF under Cooperative Agreement no. 1852977. Computing and data storage resources, including the Cheyenne supercomputer (doi: 10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. D.N. was supported by NSF award 1502806. Author contributions: Z.L., B.L.O.-B., and P.U.C. proposed the iTRACE research and designed the experiments; Z.L. and C.H. conceived this study; Z.L. wrote the paper with contributions from C.H. and P.U.C.; C.H., E.C.B., R.T., and C.Z. performed the experiments; C.H. performed the analysis with contributions from M.Y. and Y.B.; E.C.B., R.T., J.Zhu, A.J., S.G., J.Zha., J.N., and D.N. contributed to the model development; P.U.C., H.C., and Y.W. contributed to proxy interpretations. All authors discussed the results and commented on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data supporting our findings are archived at The code for present-day study is available in XCESM ( and xMCA repository (, MCA analysis). Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article