Research ArticleCLIMATOLOGY

Remote and local drivers of Pleistocene South Asian summer monsoon precipitation: A test for future predictions

See allHide authors and affiliations

Science Advances  04 Jun 2021:
Vol. 7, no. 23, eabg3848
DOI: 10.1126/sciadv.abg3848


South Asian precipitation amount and extreme variability are predicted to increase due to thermodynamic effects of increased 21st-century greenhouse gases, accompanied by an increased supply of moisture from the southern hemisphere Indian Ocean. We reconstructed South Asian summer monsoon precipitation and runoff into the Bay of Bengal to assess the extent to which these factors also operated in the Pleistocene, a time of large-scale natural changes in carbon dioxide and ice volume. South Asian precipitation and runoff are strongly coherent with, and lag, atmospheric carbon dioxide changes at Earth’s orbital eccentricity, obliquity, and precession bands and are closely tied to cross-equatorial wind strength at the precession band. We find that the projected monsoon response to ongoing, rapid high-latitude ice melt and rising carbon dioxide levels is fully consistent with dynamics of the past 0.9 million years.


The South Asian (including Indian) monsoon region is home to a fifth of the world’s population who rely heavily on seasonal precipitation for residential, agricultural, and industrial purposes (1, 2). Hence, a great deal of effort is spent assessing the future of monsoonal precipitation (35). Anthropogenic aerosols have weakened the dynamical monsoon circulation and associated precipitation in past decades (6). Nevertheless, models indicate that the thermodynamic impact of increased greenhouse gases (GHGs), as well as consequent increased atmospheric moisture content, will result in increased summer monsoon precipitation overall and extreme precipitation events (36). Hence, the moisture source regions and transport paths have become a focus of current research (7, 8). Mei et al. (8) analyzed 18 coupled model intercomparison project 5 (CMIP5) simulations to assess the dominant controls on increased 21st century continental precipitation. The four current sources of Indian summer monsoon precipitation (land: 5° to 30°N, 70° to 90°E) include the southern (S) Indian Ocean, local recycling, the Arabian Sea, and the Bay of Bengal (BoB). Of the four, only the remote (S Indian Ocean) source was shown to increase throughout the 21st century. Analysis of 15 CMIP6 models support the CMIP5 results; increased GHG radiative forcing leads to increased southern hemisphere (SH) sea level pressure and enhanced cross-equatorial flow of moisture and energy into the South Asian system (4) via the southeast trades, which, upon crossing the equator, recurve into the southwest monsoon winds and Somali Jet.

We seek to test the extent to which changes in GHG and the SH moisture source are also evident in Pleistocene reconstructions of the South Asian summer monsoon. To this end, we analyze several new 1.5–million year (Ma)–long records, including biomarker and Mg/Ca-based sea surface temperature (SST), planktonic oxygen isotopes (δ18O), benthic δ18O, seawater δ18O (δ18Osw), leaf wax carbon isotopes (δ13Cwax), and scanning x-ray fluorescence (XRF) Rb and Rb/Ca from Site U1446 on the Indian continental margin, northwest (NW) BoB (Figs. 1A and 2). These variables provide critical constraints; they offer a multiproxy means of assessing precipitation and runoff from the Mahanadi and the Ganges-Brahmaputra-Meghna (GBM), the largest drainage basin in South Asia, using δ18Osw, Rb, and Rb/Ca, as well as India vegetation structure, using δ13Cwax. Moreover, these data allow us to assess potential differences relative to runoff reconstructed from the Andaman Sea, which monitors the Irrawaddy and Salween catchments of Southeast Asia. Our NW BoB multiproxy records are reconstructed from the same site, on the same age model, greatly facilitating insights into monsoon response and phase (timing) over the Pleistocene. This allows us to more confidently assess critical differences in the phase of the precipitation isotopic response relative to precipitation amount. We accomplish these goals by quantifying the phase of these records relative to one another, external (insolation) and internal [global ice volume (IV) and GHG] forcing at Earth’s orbital eccentricity, obliquity, and precession bands. We interpret these new results in the context of existing precipitation isotope reconstructions from speleothem δ18O (δ18Oprecip) and leaf wax hydrogen isotopes (δ2Hprecip) as well as model simulations characterizing the monsoon response to external insolation forcing. We find that summer monsoon precipitation and runoff lag the precipitation isotopic response and are sensitive to wind-driven SH moisture export as well as global IV and GHG forcing, demonstrating that reconstructed Pleistocene dynamics are fully consistent with predicted future dynamics.

Fig. 1 Site locations and climatology.

(A) Location of proxy sites showing major rivers and geographic designations used in the text and figures. Site location colors denote records reflecting terrestrial precipitation isotopes (orange), wind strength (blue), and seawater isotopes (green). The same color scheme is used in Figs. 3, 4, and 6. (B) Seasonal precipitation climatology for purple box (80° to 99.5°E, 18° to 29°N). (C) BoB salinity (blue time series and box; 84.5° to 89.5°E, 16.5° to 20.5°N) compared to precipitation over land (purple time series and box) and precipitation over Mahanadi basin (D) (green time series; 82.5° to 85.5°E and 19.5° to 22.5°N). Precipitation climatology (96) and salinity climatology (97). Black outline is the Indian core monsoon zone.

Fig. 2 U1446 proxy records and spectra.

Nine hundred thousand years marks shifts in the mean, variance, or spectral nature of various proxies, occurring near the midpoint of the mid-Pleistocene transition (MPT) between ~1150 and 700 ka ago. Spectra are for the 0- to 900-ka ago interval only. Red spectrum is the same in all cases, denoting ETP, an averaged record of normalized eccentricity, obliquity (tilt; T), and precession (June 21 perihelion) spanning the past 900 ka ago. Benthic and planktic δ18O [‰ Vienna Pee Dee belemnite (VPDB)], δ18Osw (‰ standard mean ocean water), and δ13Cwax (‰ VPDB). δ18Osw is corrected for both temperature and sea level change (see Materials and Methods).

Background and motivation

The classical view of monsoons as driven by land-ocean thermal gradients may be sufficient to explain seasonal wind reversals but is insufficient to explain monsoonal precipitation (9). Deep convective precipitation requires surface air masses with sufficient buoyancy (moisture) to rise through the troposphere, during which buoyancy is further enhanced through latent heating associated with condensation. In this case, subcloud moist entropy are primary drivers, and tropical, deep convective South Asian precipitation is associated with extreme regional migration of the intertropical convergence zone beyond its typical location (within 10° of the equator), into the BoB and the South Asian continent [(9, 10) and references therein]. When this occurs, large-scale Hadley circulation is transformed into ascending flow in the northern hemisphere (NH) and descending flow in the SH with strong northerly, cross-equatorial transport of moisture and energy in the surface wind field.

Rigorous debate exists regarding the extent to which cross-equatorial transport was important in the past. On the basis of speleothem δ18Oprecip, Zhang et al. (11) and Kathayat et al. (12) contend that the South Asian summer monsoon predominantly and directly responds to changes in external forcing (NH summer insolation) with SH processes playing no significant role at orbital time scales. This is in contrast to reconstructions derived from other archives that indicate significant sensitivity to internal forcings and SH climate dynamics (1316). These divergent interpretations lead to substantial uncertainty about the extent to which the predicted drivers of future changes in South Asian monsoon precipitation operated in the past as well. Specifically, (i) the extent to which orbital-scale variability in South Asian monsoon precipitation is driven primarily by external insolation forcing or by internal feedbacks within the climate system including global IV and GHG concentrations, and (ii) the extent to which the South Asian monsoon precipitation is, or is not, influenced by cross-equatorial moisture and energy transport from the SH. Site U1446 was drilled in the core convective region of the Indian summer monsoon to specifically address these questions; here, we present new records and interpret them in the context of extant, spatially distributed South Asian records. Resolution of the divergent interpretations stems from our finding of clear proxy groupings; those recording changes in the isotopic composition of precipitation, primarily reflecting changes in moisture source area and transport path dynamics, and those primarily recording more regional precipitation amount and runoff from river drainage basins.

Site U1446 climatology and oceanography

The BoB, surrounded on three sides by land, is the world’s least saline ocean basin, receiving ~2950 km3 of runoff annually (17). The entire system is strongly seasonal (Fig. 1B) with 80 to 90% of precipitation resulting from summer monsoon (June to September) precipitation (18). Evaporation plays a minor role with precipitation plus runoff minus evaporation [(P + R) – E] being positive throughout the year (19). Runoff is dominated by GBM discharge peaking at ~1 × 105 m3/s between August and September (20), approximately a month after peak precipitation over the surrounding landmass in July (21). Hence, precipitation and runoff are strongly coupled to salinity (Fig. 1, C and D).

Ensemble model data indicate that the S Indian Ocean contributes the most summer season moisture to India, representing 46%, followed by local recycling (30.4%), the Arabian Sea (19.6%), and the BoB (4%) (8). Reanalysis data also indicate that the S Indian Ocean contributes the most summer season moisture to India, representing 55%, followed by local recycling (23%), the Arabian Sea (15%), and the BoB (7%) (8). Mei et al. (8) point out that these results are not in conflict with those describing the considerable influence of BoB synoptic low-pressure systems and monsoon depressions as contributing 40 to 50% of summer monsoon precipitation (22), differentiating the ultimate SH source of the moisture from the systems delivering it to continental India. More recent investigations indicate an even greater role of synoptic lows and depressions, accounting for 60 to 80% of total precipitation in the core monsoon zone (23) and up to 80% of total Indian precipitation (22, 24). BoB precipitation begins between late April and mid-May, reaching ~10°N between mid-May and mid-June, with strong precipitation over the northern (N) BoB and India from June to September [(10) and references therein].

Site U1446 (19°5.01′ N, 85° 44.08′ E) is located in the Mahanadi offshore basin at 1430 m below sea level on the Indian margin (Fig. 1) (25). The site is on a NW-southeast trending ridge, shielding it from turbidite deposition. The 204-m composite section, constructed from three offset holes, spans the past 1.46 Ma with an average sedimentation rate of ~14 cm per thousand year (ka). Sediments are composed of dark gray to gray hemipelagic clay with nannofossils and foraminifers. The lithogenic component is derived from the Mahanadi River (26) basin with over 90% of the annual sediment delivered during the summer monsoon (June to September) (27).

The site is strategically located 70-km offshore of the Mahanadi River Delta beneath the N limb of the seasonally reversing East India Coastal Current (EICC). The EICC flows northward before and during the summer monsoon, reversing direction after the summer monsoon and providing for transport of fresh water out of the BoB. Interannual salinity variability within the N BoB and the ~100-km-wide EICC is directly linked to runoff, having salinities up to 10 practical salinity units fresher than those offshore, with little impact from direct precipitation (21, 28). Within the N BoB, (north of 17°N), Ganges runoff alone yields close to 3 m/year of sea level equivalent (29). Regional ocean modeling indicates that the spatial and temporal variability in seasonal sea surface salinity, sea level anomalies, the organization and strength of the equatorward flowing EICC, and aspects of coastal Kelvin wave dynamics are not well replicated without including river runoff, highlighting the critical impact of runoff on BoB oceanography (20).

Proxy records

The U1446 time scale is derived by six radiocarbon ages on the upper 8 m of the section and by correlation of benthic δ18O to the LR04 global benthic stack (Fig. 2) (30), yielding high coherence (≥0.91) and near-zero phase (≤1 ka) across all orbital bands (see Materials and Methods). On this age model, we present six new records including three that independently monitor aspects of summer monsoon precipitation and runoff from the largest drainage basin in S Asia.

NW BoB δ18Osw is reconstructed by removing the local SST and global δ18Osw signals from the planktonic Globigerinoides ruber δ18O (see Materials and Methods). BoB δ18Osw is linearly related to salinity (19, 31, 32), reflecting both direct precipitation and runoff from surrounding catchment basins (Fig. 1). The strong influence of runoff is evident from the steep zonal and meridional salinity gradients observed presently (21, 28) and in paleoreconstructions (33). For the N BoB, where most runoff enters, the impact of direct precipitation on interannual salinity variability is minimal, and runoff dominates (21, 28). The impact of potential changes in isotopic composition of runoff is small. Runoff volume at the mouth of the Ganges (−6‰) (34) would have to change by only 17% to alter EICC δ18O by 1‰, as is observed at the precession time scale. In contrast, the δ18O of runoff [currently 19% of EICC volume (35)] would have to change by 5.3‰ (−6 to −11.3‰) to alter the EICC δ18O by 1‰. This would require that runoff at the mouth approaches the isotopic composition (−11 to −16‰) of the headwaters above 2000-m elevation that are currently composed of ~44% glacial meltwater (36). Hence, the U1446 δ18Osw record is interpreted as salinity, dominated by continental precipitation and runoff from the extensive South Asian river catchments (GBM and Mahanadi) with smaller contributions from direct precipitation or changes in the isotopic composition of runoff. This interpretation is supported in the following sections by comparison with trace-element XRF and δ13Cwax reconstructions that are independent of changing isotopic composition of precipitation and runoff.

XRF records (see Materials and Methods) are interpreted following Phillips et al. (37) from nearby site NGHP-01-19. Changes in Ca are dominated by changes in surface ocean biogenic CaCO3 production (37); Ca, CaCO3 percentage, and CaCO3 mass accumulation rate (MAR) increase during glacial relative to interglacial intervals. Phillips et al. (37) interpret increased CaCO3 MAR to reflect reduced surface ocean stratification; decreased fresh-water input during weak glacial-age monsoons leads to reduced stratification, increased nutrient access, and increased primary production. Lithogenic MAR and sedimentation rates increase during interglacials (high sea level), indicating that continental slope sedimentation is influenced more by monsoon precipitation and runoff than by shelf exposure and erosion during glacial sea level low stands. Rb is an element common in felsic rocks and varies as a function of fine-grained lithogenic input, increasing with increased runoff. As with lithogenic MAR, Rb increases during interglacial intervals, in opposition to Ca (fig. S1, a and c). Hence, the Rb/Ca ratio (Fig. 2) represents the contribution ratio of fine grained lithogenic input to marine biogenic carbonate and is interpreted as a proxy for runoff and runoff-induced stratification. Both Rb and Rb/Ca are independent of factors that potentially affect water isotope proxies such as changing moisture source areas and changing isotopic composition of precipitation and runoff. Rb/Ca and δ18Osw are coherent at the eccentricity and precession bands, whereas Rb and δ18Osw are coherent at the eccentricity, obliquity, and precession bands (table S1), further supporting the interpretation as runoff and runoff-induced stratification proxies.

δ13Cwax reflects the relative contribution of C3 (−37.1 ± 2 ‰) and C4 (−19.5 ± 1.8‰) plant endmembers (38). In general, the proportion of C3 and C4 vegetation is affected by air temperature, water availability, and CO2 concentration (39); C4 plants have a lower CO2 compensation point for carbon fixation relative to C3 plants (40). In this region, at glacial-interglacial time scales (eccentricity and obliquity), this proxy is strongly sensitive to SST and atmospheric CO2 as demonstrated by high coherence (table S1) and near-zero phase (see Discussion). At the precession band, water availability (including seasonality) plays a greater role, as demonstrated by coherence and phase relationships with δ18Osw, Rb, and Rb/Ca (see Discussion). Table S2 provides a summary of these and other proxies discussed in subsequent sections.

Interpretive paradigm: Spectral structure and phase

The study of Pleistocene orbital-scale climate change is unique. The signal-to-noise ratio is large, and the external forcing (insolation) is known (41). In addition, the two critical internal drivers (GHG concentration and global IV) are well constrained in the paleo record on the basis of ice cores and benthic δ18O (30, 42, 43). Hence, it is possible to derive the relative sensitivity of monsoon proxy records to these potential drivers by examining the spectral structure, coherence, and phase of time series reconstructions (44).

From a modeling perspective, Erb et al. (45) showed that first-order aspects of monsoons are reasonably well captured by linear reconstructions using idealized single-forcing simulations that include obliquity, precession, CO2, and ice sheets. Results indicate that precipitation over the Asian landmass increases with increased summer insolation, increased GHG concentrations, and decreased global IV; hereafter, we refer to the internal forcing related to the coupled IV minima and GHG maxima as IVGHG. Modeling of the South Asian monsoon (46) on Pleistocene time scales indicates that SST is equally sensitive to IVGHG and insolation; continental temperatures are found to be 37% more sensitive to IVGHG than to insolation and continental precipitation are found to be 67% more sensitive to insolation than to IVGHG. Although these simulations included large portions of the Tibetan Plateau and, hence, may not be directly applicable to the more southerly regions discussed in this work, they nevertheless indicate that the South Asian system is strongly affected by both internal (IVGHG) and external (insolation) drivers. We test these sensitivities by reconstructing precipitation and runoff into the BoB, quantitatively assessing spectral coherence and phase relationships.

At orbital time scales, interpretation of proxy records derives from analysis of spectral structure, as well as coherence and phase relationships. Spectral structure involves the relative concentrations of variance within the three primary orbital bands including eccentricity (124- and 95-ka periods), obliquity (also referred to as tilt; 41-ka period), and precession (23- and 19-ka periods) (41); note that the 124- and 95-ka eccentricity bands are most often referred to as the ~100-ka band, as will be the case here. Coherence and phase are measures of the linear correlation and timing between time series at each frequency band. For example, climate proxies that have similar relative amounts of variance in the eccentricity, obliquity, and precession bands as is found in global IV (and are in phase with or slightly lag IV) are interpreted as being strongly influenced by changes in high-latitude IV. Global IV, in turn, is known to be strongly coupled with changes in GHG; the records are coherent at all three orbital bands, with IV minima occurring slightly after maxima in GHG (table S1). In contrast, records that are dominated by variance in the precession band and are in phase with or slightly lag precession maxima (Pmax) (warm SH summers) or minima (warm NH summers) are interpreted as being forced dominantly by external insolation, based on the fact that precession-band variance dominates the summer insolation spectrum at all latitudes up to ~65° (fig. S2).

Coherence and phase relationships among monsoon records and potential internal and external drivers are quantitatively established by conventional time series analytical techniques (see Materials and Methods). The results of these analyses are summarized on phase wheels (Fig. 3A), providing a useful means of assessing the underlying physics of the climate system; the timing of a proxy response relative to known external (i.e., insolation) and internal (e.g., IVGHG) forcing provides a measure of sensitivity to those climate drivers. For example, if a climate proxy responds dominantly and directly to NH summer insolation, then its vector would plot at or near 0° on the precession and obliquity phase wheels (Fig. 3, C and D). Phase is easily converted to time; for example, an obliquity-band phase lag of 70° is equivalent to a 7.9-ka lag (70°/360° × 41 ka = 7.9 ka). Significant leads or lags relative to external forcing maxima indicate the influence of other (e.g., internal) factors.

Fig. 3 Phase wheel summaries of cross-spectral coherence and phase relationships at orbital frequency bands.

Zero phase is set at June 21 perihelion for precession, denoting the timing of minimum precession (Pmin), maximum obliquity, and maximum eccentricity. (A) Positive phase (counterclockwise) represents leads, and negative phase (clockwise) represents lags. Red text and dots represent the timing of external forcing factors (hemispheric and seasonal insolation maxima and minima; at 0° and 180°) and internal forcing factors (global IV minima and GHG maxima; CO2). Vectors represent the timing of monsoon proxy responses within the (B) 100-ka eccentricity cycle, (C) the 23-ka precession cycle, and (D) the 41-ka obliquity cycle (averaged over the number of cycles present in a given record). All phase vectors plotted are coherent (>80 CI) with eccentricity, obliquity, and precession, with one exception; NW BoB δ18Osw coherence does not exceed 80 CI at eccentricity (0.78) but is strongly coherent with IV (0.90), so this phase is plotted relative to ice mininum. Phase errors are typically ±10° to 20° (table S3). Color-coded text at 0° and 180° (as in Figs. 1, 4, and 6) indicates the timing of model simulated responses to insolation-only forcing, green for maximum precipitation, blue for maximum winds, and orange for precipitation isotope minima. Proxy response vectors are similarly color coded (e.g., green for maximum precipitation/runoff). The gray vector and dashed line indicate the timing of SST maxima and G. ruber δ18O minima, respectively. SITIG, summer intertropical insolation gradient.


All NW BoB records exhibit the well-known shift in dominant periodicity from 41-ka cycles to 100-ka cycles across the mid-Pleistocene transition (MPT) at ~900 ka ago (47), with the 100-ka variance dominant since then (Fig. 2). The G. ruber planktonic δ18O record has a structure similar to the global benthic record but with 3.2 times the variance after ~900 ka ago and 5.9 times the variance before ~900 ka ago, signifying a great deal of BoB surface water variability throughout (fig. S1a). The planktonic δ18O mean shifts by +0.25‰ at ~900 ka ago, mostly attributed to enrichment of the heavier isotopes during glacial intervals; interglacial values remain at ~−3.5‰ over the entire record. SST has a glacial-interglacial range of ~3°C, similar to that documented for the Andaman Sea, eastern (E) BoB (13). SST has no substantial change in the mean ~900 ka ago. δ18Osw has a range of ~2‰, exhibiting a shift (+0.19‰) at ~900 ka ago, slightly less than the +0.25‰ shift in planktonic δ18O. Neither Rb (fig. S1c) nor Rb/Ca (Fig. 2) exhibit substantial mean shifts at ~900 ka ago although the amplitude of Rb/Ca increases. Further discussion, including cross spectral coherence and phase analyses, are restricted to the past 900 ka ago such that variance before the MPT is not included. This is because the records available for comparison to these data are all shorter and do not include pre-MPT variance.

The spectra of all precipitation and runoff proxies are dominated by 100-ka eccentricity band variance with lesser amounts of 41- and 23-ka variance, similar to that of CO2 and global IV (benthic δ18O) but very different from NH summer insolation, which is dominated by precession throughout the low and mid-latitudes, with increasing amounts of obliquity in the high latitudes (fig. S2). Coherence and phase results are reported in tables S1 and S3 and summarized in Fig. 3. The three runoff and precipitation proxies [δ18Osw, XRF (Rb and Rb/Ca averaged), and δ13Cwax] form clusters at all three orbital bands, the cluster averages (±1σ, n = 3) lag eccentricity maxima by 21° ± 28° (6 ± 8 ka), tilt maxima (Tmax) by 97° ± 20° (11 ± 2 ka), and precession minima (Pmin) by 132° ± 6° (8 ± 0.4 ka). At all orbital bands, the SST and G. ruber δ18O records are in phase with or very slightly lag IVGHG, whereas, after removing the SST and global IV signal from the G. ruber δ 18O record, the resultant δ18Osw record significantly lags IVGHG but is in phase with δ13Cwax, Rb, and Rb/Ca at the precession band and in phase with Rb and Rb/Ca at the obliquity band.


Several orbital-scale multiproxy reconstructions exist for South Asia, monitoring aspects of wind strength, precipitation isotopic composition, and precipitation/runoff (table S2). These results are added to the phase wheels (Fig. 4) and integrated into the following discussion.

Fig. 4 As in Fig. 3 with the addition of proxy vectors from sites other than U1446 (Fig. 1).

Proxy responses interpreted as responding to winds (blue), precipitation/runoff (green), and precipitation isotopic composition (orange) follow color coding of insolation-only simulation results (text at 0° and 180°) and that in Figs. 1, 3, and 6. Black vectors denote phase of the first derivative of the IV record, reflecting sea level (SL) change. The modifier ‘-TC’ in panel D indicates measurement at the thermocline, as opposed to measurement within the mixed-layer.

Spectral structure

All NW BoB (U1446) precipitation and runoff records (δ18Osw, Rb, Rb/Ca, and δ13Cwax) are dominated by variance at the eccentricity (~100 ka) band (Fig. 2 and fig. S2) with lesser amounts of variance at the tilt (41 ka) band, as are the Pleistocene global IV and GHG records. Furthermore, all are coherent with global IV [>0.90 CI (confidence interval)] and with CO2 (> 0.80 CI) at both the eccentricity and obliquity bands (table S1). Dominance of the 100-ka variance indicates that summer monsoon precipitation is strongly influenced by large-scale global boundary conditions associated with coupled IVGHG feedbacks, consistent with model simulations (45). Dominance of eccentricity and obliquity in the spectral structure of the E BoB (Andaman Sea) δ18Osw precipitation/runoff record (Fig. 1) also demonstrates sensitivity to these drivers (13), documenting a robust, reproducible response spanning the Indian core monsoon zone through the E BoB.

These findings are inconsistent with interpretation of the Bittoo and Xiaobailong (XBL) South Asian speleothem δ18O reconstructions (Fig. 1). Speleothem δ18O reflects the isotopic composition of local precipitation and changes as a function of source area, evaporation, rainout, and convection type along the transport path, as well as local precipitation intensity and the interplay of meteoric waters with the soil and karst environment [(11, 48) and references therein]. Bittoo and XBL δ18O are interpreted as representative of both the South and East Asian (EA) summer monsoons and as responding directly to NH summer insolation (12). In contrast, a longer, more continuous NW BoB reconstruction of South Asian precipitation isotopic composition derived from leaf wax δ2H (δ2Hprecip) documents strong sensitivity to internal drivers including CO2 and IV (49). Below, we examine phase relationships and model simulations, providing evidence that South Asian precipitation and runoff lag changes in the isotopic composition of precipitation (speleothem δ18O and leaf wax δ2H) and are sensitive to both IVGHG boundary conditions and SH moisture export associated with large-scale wind patterns.

Eccentricity phase: Termination dynamics

The direct impact of eccentricity on insolation is negligible (41); hence, phase is interpreted in the context of linkages to coupled IVGHG dynamics, which are responsible for the bulk of variance across Pleistocene climate. U1446 proxies plot near the top of quadrant 1 on the eccentricity phase wheel (Fig. 4B). The proxies linked to vegetation (δ13Cwax), terrestrial precipitation isotopic composition (δ2Hprecip), fine-grained lithogenic runoff (Rb), and stratification (Rb/Ca) have zero or small phase lags with respect to IVGHG, whereas the NW and E BoB δ18Osw (precipitation/runoff) proxies have larger lags. The difference is anticipated given that the IV signal has been quantitatively removed to derive the δ18Osw records. Overall, the precipitation isotope proxies (orange vectors) and the rainfall/runoff proxies (green vectors) are not well differentiated, indicating that precipitation isotopic minima and precipitation/runoff maxima both occur at or slightly after the coupled IVGHG forcing.

The paired nature of the U1446 data allows us to assess the dynamics of South Asian δ18Osw relative to IV terminations (benthic δ18O) within the same core (Fig. 5A) and, therefore, facilitate comparison with those demarcated on the basis of EA speleothem δ18Oprecip from the Yangtze River Valley (YRV) [see figure 1 in (50)]. The weakest EA monsoons (speleothem δ18Oprecip) occur, on average, 3 ka after glacial maxima, and the midpoint of the ensuing rapid transition toward interglacial values occurs rapidly thereafter, on average, ~1 ka after the midpoint of IV terminations (Fig. 5C). Cheng et al. (50) demonstrated that each IV termination is associated with a precession cycle. Since the EA YRV record is almost exclusively dominated by precession-band variance, the rapid transition to interglacial values is part of the precession-band variance, although influenced by IV as well at terminations. The dynamics have been interpreted in the context of changing and varied moisture sources, trajectories, and transport over wide continental shelf areas (51); moisture sources include the Arabian Sea, BoB, South China Sea, and North Pacific with trajectories that cross the wide exposed (glacial) and submerged (interglacial) continental shelf regions of Indonesian and the East China Sea.

Fig. 5 Termination dynamics.

(A) NW BoB precipitation/runoff (δ18Osw; red) relative to IV (benthic δ18O; blue). Terminations (T) 1 to 10 are denoted by vertical dashed blue lines. The midpoints of subsequent precipitation/runoff transitions to interglacial values are denoted by vertical dashed red lines. Red and blue dots indicate IV maxima and precipitation/runoff minima. (B) Same as in (A) but focusing on the ±15 ka surrounding IV terminations, for which the midpoints have been set to zero, with negative values indicating younger time. Thin dashed lines are individual records that have been averaged to yield thick lines. Arrow indicates the average lag between midpoints of IV terminations and midpoints of the subsequent proxy transitions. (C) Same as in (B) but showing EA YRV speleothem δ18O for T1 to T7 (50).

In contrast, the weakest South Asian precipitation/runoff (δ18Osw) occurs, on average, 8 ka after glacial maxima, and the midpoints of the ensuing rapid transitions toward interglacial values occur, on average, 6 ka after IV terminations (Fig. 5B). Moisture sources and trajectories for the South Asian system are less variable and more independent of shelf emergence and submergence (the Indian shelf is narrow). Hence, the delayed increase in South Asian precipitation/runoff relative to terminations is regional in nature, and the two subsystems are uncoupled in this context. We show in the next section that South Asian precipitation/runoff has a phase lag of 132° (8 ka) relative to Pmin. Hence, the average 6-ka lag in intensification following IV terminations is largely accounted for by processes responsible for the precession band phase lag as opposed to dynamics related to glacial terminations.

Precession phase

δ18Osw is constructed by removing the calcification temperature and global IV (benthic δ18O) signals from the G. ruber δ18O record (see Materials and Methods). The impact on phase is the same across all orbital bands; the “raw” G. ruber δ18O phase (dashed lines in Fig. 4, B to D) only slightly lags IVGHG and the terrestrial precipitation isotope records (δ18Oprecip and δ2Hprecip). Removal of the local SST and global δ18O signals (recorded in the same samples) causes the phase lag to increase significantly, falling in phase with δ13Cwax, Rb, and Rb/Ca at the precession band and in phase with Rb and Rb/Ca at the obliquity band. The finding that δ18Osw reconstructions from the NW and E BoB (13) are in phase with independent monsoon proxies not affected by the isotopic composition of precipitation and runoff supports their interpretation as salinity proxies driven by runoff. If meltwater were a dominant driver, then the phase would fall near the first derivative of the benthic δ18O signal (maximum ice melt, at −7°); instead, both δ18Osw records fall within the cluster centered on −132°. If the δ13Cwax response were dominated by GHG, then it would be approximately in phase with or slightly lag GHG maxima (~−70°), whereas if it were dominated by temperature, then it would likely have the opposite phase (~110°). Again, δ13Cwax falls within the cluster centered on −132°, consistent with the δ18Osw and Rb/Ca proxies.

Precession-band phases fall into two distinct clusters, delineating continental precipitation isotope proxies (orange vectors; δ2Hprecip and δ18Oprecip) that plot in quadrant 1 from all other South Asian summer monsoon proxies (blue and green vectors) that plot in quadrant 2 (Fig. 4C). The three continental precipitation isotope records cluster near IVGHG (−71° ± 9°), indicating strong sensitivity to these internal drivers (49). This phasing is consistent with previous interpretations of these records and with simulations indicating more distal sources, longer transport paths, and lighter precipitation δ18O during Pmin compared to Pmax (52), predominantly reflecting large-scale circulation versus more local to regional rainfall amount. Similar source area and transport dynamics have been put forth to interpret East Asia YRV speleothem δ18O variability (53, 54). In contrast, NW and E BoB proxies that reflect lithogenic runoff and surface-water stratification (Rb and Rb/Ca), seawater isotopic composition (δ18Osw), and Indian vegetation structure (δ13Cwax) all cluster in quadrant 2 (−132° ± 6°), between IVGHG and Pmax (Fig. 4C), indicating different sets of drivers relative to the precipitation isotope proxies. These precipitation/runoff proxies are also in phase with two stacked wind reconstructions generated from four records in the N Arabian Sea [(15, 55) and references therein] and five records in the western (W) Arabian Sea [(14) and references therein] (Fig. 1). These wind reconstructions reflect the strong influence of the low-level southwesterly (Somali) Jet that drives Ekman-induced upwelling, strong biological productivity, denitrification in the intermediate waters beneath, and eolian transport of lithogenic materials from the surrounding desert to these Arabian Sea sites. Although there is negligible precipitation in these regions, modern dynamics suggest a close association between the onset of summer monsoon rains over India and the abrupt strengthening of the Somali Jet over the Arabian Sea (56). We therefore interpret the clustering of these multiple independent proxies (at an average of −132°) as indicating that the large-scale winds in the N Arabian Sea, W Arabian Sea, and S BoB (16) are linked to SH extraction and cross-equatorial transport of moisture into the South Asian system driving increased precipitation and runoff, akin to modern precipitation dynamics (7, 8, 2224, 57). The inferred linkage of winds and precipitation in the proxy records is consistent with the Jalihal et al. (58) finding linking changes in BoB wind strength to latent heat flux and precipitation in the EC (European Community)–Earth simulations.

The multiproxy phase response at −132°, between the timing of IVGHG (~−73°) and December 21 perihelion (−180°), is interpreted within the context of model results (58) as the timing of maximum total column energy (Qdiv) over the BoB, if the impact of IV and GHGs were included. Total column energy over the BoB is found to be dominated by latent heat flux, which, in turn, is driven by the changes in the wind-field. In this context, the multiproxy phase response at −132° is interpreted as the combined response to maximum surface sensible heating (IVGHG) and the maximum import of moisture (latent heat) from the SH Indian Ocean at December 21 perihelion (−180°), as follows. The seasonal cycle at December 21 perihelion (Pmax) is characterized by warm SH summers and cold SH winters. The warm SH summer results in increased energy storage across the SH Indian Ocean, preconditioning the system for enhanced latent heat export at the onset of the NH summer monsoon in May. Decreased insolation during the SH winter (June to August) results in cold, dry winds over warm ocean waters, maximizing the ocean-to-atmosphere transfer of latent heat available for cross-equatorial transport and release over South Asia during the summer monsoon (14, 15). This hypothesized mechanism is supported by model sensitivity experiments (59) and modern observations (7) and is similar in nature to those found in CMIP5 and CMIP6 models that project future monsoon intensification. The −132° phase precludes NH summer insolation as the direct driver of South Asian summer monsoon precipitation. This interpretation is in contrast to that based on Bittoo and XBL South Asian monsoon speleothem δ18O, as responding directly to NH summer insolation (NHSI) with no significant phase lag (11, 12, 51) and no contribution of moisture from the SH (12).

Proxy spatial coverage is now sufficient to make site-specific proxy-model comparisons. Bosmans et al. (60) synthesized results of insolation-only model simulations (precession and obliquity) run on the EC-Earth, GFDL, CESM, and HadCM3 models (obliquity experiment only). Our results are interpreted in the context of these idealized model experiments, keeping in mind two important points. First, the ability of climate models to capture subgrid-scale dynamics deteriorates at coarser resolutions due to deficiencies in simulating land-surface feedbacks, water vapor and cloud feedback processes, and factors related to topography [(61) and references therein]. While various downscaling approaches and parameterizations have been proposed to compensate for these issues, they have been met with limited success; Indian precipitation biases remain (62), even within the core monsoon zone using the latest (CMIP6) multimodel analyses (63). Simulating synoptic-scale convective responses to large-scale forcing remains a fundamental challenge in climate models (64). Second, these idealized precession and obliquity experiments do not incorporate IV or GHG as boundary conditions, both of which significantly lag NH summer-insolation forcing at the obliquity and precession bands. Hence, it follows that simulations forced only with orbital-driven insolation changes will be deficient in eliciting longer-term lags in the climate system; idealized model responses (maximum winds, precipitation, etc.) are necessarily constrained to be in phase with orbital maxima or minima. Nevertheless, these idealized experiments are very useful in contextualizing the response to the external forcing alone, helping to evaluate both the model results and the observed proxy phases.

All three models (EC-Earth, GFDL, and CESM) show increased summer monsoon precipitation over N India at Pmin, particularly in the Himalayan forefront, an orographic effect common to all models (Fig. 6). GFDL has slightly increased precipitation at Pmax in the Indian core monsoon zone, while CESM shows slightly increased precipitation at Pmax in S peninsular India. Nevertheless, the overall insolation-only model result for continental areas under the influence of the South Asian summer monsoon is increased precipitation at Pmin (0° in Fig. 4C). In contrast, all three models indicate maximum precipitation over the BoB at Pmax, corresponding to 180° on the precession phase wheel. Hence, idealized simulations indicate a 180° phase difference between the precession-band timing of maximum summer monsoon precipitation over the ocean versus the nearby continents. Jalihal et al. (58, 65) discuss the mechanisms underlying these simulation results. They find that precipitation over land is more sensitive to radiation forcing (net top of atmosphere long wave and short wave), whereas surface latent heat flux (dominated by winds) and vertical stability are more important over the BoB, again, keeping in mind that IV and GHG forcing are not included in these idealized simulations.

Fig. 6 Differences in average precipitation (June to August, mm per day) for precession (Pmin-Pmax) and obliquity (Tmax-Tmin) for each model.

(A, B) European Community (EC). (C, D) Geophysical Fluid Dynamics Laboratory (GFDL). (E, F) Community Earth Systems Model (CESM). (G) Hadley Centre Coupled Model (HadCM3). Contours indicate values for Pmax and Tmin. Thick contour indicates 4000-m elevation for the Tibetan Plateau. Positive values (blue shading) indicate increased precipitation at Pmin and Tmax. Negative values (red shading) indicate increased precipitation at Pmax and Tmin. Figure after Bosmans et al. (60) Fig. 3 with proxy site locations superimposed and Indian core monsoon zone (63) outlined in (C) and (D). Site location colors denote records reflecting terrestrial precipitation isotopes (orange), wind strength (blue), and seawater isotopes (green), consistent with Figs. 1, 3, and 4.

NW BoB proxy indicators for runoff and stratification (Rb and Rb/Ca) and vegetation structure (δ13Cwax) have similar phasing as the δ18Osw salinity proxies. Hence, we conclude that maximum continental precipitation does not occur near Pmin (0°) but, rather, at ~−132°, closer to the idealized model result for the timing of maximum precipitation over the BoB at Pmax. This may point toward model deficiencies in accounting for the contribution of BoB synoptic low-pressure systems and monsoon depressions in driving continental precipitation, possibly due to insufficient model resolution, and deficiencies in tropospheric temperatures or SST biases [(62, 66) and references therein]. Acosta and Huber (66), for example, find that only high-resolution models (½° or smaller spatial resolution) capture the Indo-Gangetic low-level jet responsible for bringing moisture into India from the BoB. We suggest that high-resolution models run in transient mode, and incorporating known IVGHG boundary conditions would simulate less divergent timing of precipitation over land and ocean, converging on a phase near -132° at the precession band. We note that GFDL is most consistent with our proxy phase results, simulating a weak increase in precipitation over the India core monsoon zone and Indochina Peninsula at Pmax.

The proxy result of strengthened Arabian Sea and BoB summer monsoon winds at −132° (closer to Pmax than Pmin) is difficult to assess relative to idealized model results. All three models indicate strengthened 850 hPa winds at Pmin in the N Arabian Sea and Arabian Peninsula but in different locations and with transitions from stronger at Pmin to stronger at Pmax taking place in the vicinity of the proxy core sites (fig. S3 a, c, e). Hence, proxy-model wind comparisons are not definitively convergent or divergent in the N Arabian Sea region. However, all models do show strengthened winds at Pmax over the S Arabian Sea and BoB, consistent with the S BoB wind reconstruction (Fig. 4C) (16).

Obliquity phase

The NW BoB proxies that reflect lithogenic runoff and surface-water stratification (Rb and Rb/Ca), seawater isotopic composition (δ18Osw), and Indian vegetation structure (δ13Cwax) cluster at −97° ± 20° (Fig. 4D). The NW BoB δ13Cwax record (in phase with the δ18Osw and Rb/Ca at the precession band) plots early, in quadrant 1; the high coherence with, and close proximity to, IVGHG indicates sensitivity to IVGHG at the eccentricity and obliquity (tilt) bands. The E BoB δ18Osw record from the Andaman Sea has a phase only slightly lagging obliquity minima (Tmin), significantly different from the NW BoB cluster. This may be accounted for by the impact of obliquity-band sea level changes in the Andaman, as hypothesized by Gebregiorgis et al. (13); it is in phase with the maximum rate of sea level fall and physical isolation of the Andaman Sea, which would reduce the export of runoff to the larger BoB and reduce Andaman salinity. In any case, the runoff/stratification proxy cluster (−97° ± 20°) significantly lags the South Asian (SA) precipitation isotope records, which plot at −52° ± 3°. The −97° phase indicates that strengthened Indian precipitation/runoff is not a direct response to NH summer insolation at the obliquity band nor is it directly driven by the strengthened summer intertropical insolation gradient (67). Our interpretation, in the context of simulations and other proxy records, follows.

Simulated precipitation responses to changes in obliquity are weaker and more variable than at the precession band. Obliquity maxima (Tmax) are characterized by increased precipitation over N India in all four models, again, an orographic effect held in common (Fig. 6). GFDL shows a weak tendency for increased precipitation in the Indian core monsoon zone and Indochina Peninsula at Tmin, while CESM and EC-Earth show increased precipitation in central and S peninsular India at Tmin, a pattern somewhat different from HadCM3 and GFDL.

The −97° proxy phase indicated by NW BoB δ18Osw, Rb/Ca, and Rb is more consistent with the timing of model precipitation over the BoB where all except HadCM3 show increased BoB precipitation at Tmin (Fig. 6). Increased BoB precipitation at Tmin is internally consistent, with all models showing BoB-wide Tmin responses characterized by increased evaporation [figure C.8 in (60)], increased SST [figure C.14 in (60)], and decreased sea level pressure [figures 6a, C.4a, C.5a, and C.6a in (60)].

Similar to our precession-band interpretation, we interpret the −97° Indian margin phase as direct BoB precipitation and continental runoff having similar obliquity-band timing. In the context of model results, we interpret the −97° phase as consistent with stronger simulated BoB precipitation at Tmin but taking place earlier in the tilt cycle due to the influence of IVGHG forcing. We hypothesize that maximum summer monsoon precipitation at −97° is sensitive to both (i) internal IVGHG radiative forcing at −65° (warm surface ocean and increased atmospheric water content) and (ii) external (obliquity) forcing toward increased BoB precipitation at Tmin, possibly associated with tropical summer insolation maxima at Tmin. Again, these results suggest that the increased BoB precipitation is propagated onto the continents via synoptic lows and depressions, as in the modern. As with precession, the GFDL model response is most consistent with the proxy phase response in terms of increased precipitation over land closer to Tmin than Tmax.

Unlike the precession-band results, BoB precipitation proxies are not in phase with wind-strength proxies in the N and W Arabian Sea (both of which are found in quadrant 1). Nevertheless, the Arabian Sea wind strength proxies are consistent with model results; all four models show increased N Arabian Sea wind strength at Tmax (fig. S3). We interpret this as indicating a decoupling of N Arabian Sea winds and moisture transport into India at the obliquity band, as has been observed during anomalous monsoon years on interannual time scales (68). NW BoB δ18Osw and Rb/Ca have the same phase as the S BoB wind proxy suggesting a link to BoB wind forcing (Fig. 4D). While weak, the model simulations also show increased S BoB winds at Tmin, consistent with the proxy phase.


Our results indicate that Pleistocene precipitation/runoff amount is decoupled from changes in the isotopic composition of that precipitation at orbital time scales, a result that is consistent with isotope-enabled model simulations. Seawater δ18O is uncoupled from terrestrial water isotope proxies (speleothem δ18Oprecip and leaf wax δ2Hprecip) at the obliquity and precession bands. At these bands, the consistent response of δ18Osw and non-isotopic indicators of terrestrial runoff and vegetation leads to the conclusion that these proxies primarily reflect changes in terrestrial precipitation and runoff, as opposed to changes in the isotopic composition of precipitation, driven primarily by source area and transport path variability now thought to dominate the timing of Asian terrestrial water-isotope proxies. Precipitation maxima occur well after NH summer insolation maxima at both orbital bands (−132° relative to Pmin and −97° relative to Tmax), precluding external insolation as the primary driver of changes in South Asian summer monsoon precipitation and runoff. Proxy-inferred precipitation maxima occur between the timing of coupled IVGHG forcing and the timing of model-inferred precipitation maxima over the BoB (at Tmin and Pmax); BoB precipitation maxima are, in turn, driven by increased total column energy, latent heat flux, and large-scale wind fields. Consistent timing with an array of proxies for wind strength in the Arabian Sea and BoB supports the strong link between precipitation and large-scale winds in cross-equatorial moisture transport. The dominance of 100-ka spectral power in precipitation/runoff proxies further supports sensitivity of South Asian summer monsoon precipitation to decreased IV and increased CO2. These summer monsoon dynamics, evident in our Pleistocene multiproxy reconstructions, are fully consistent with CMIP5 and CMIP6 simulations linking future increased monsoonal precipitation and moisture source area changes to increased GHG concentrations.


Stable planktonic isotopes

U1446 was sampled at 30-cm resolution from 0 to 204 composite meters below sea floor (mbsf). Samples were freeze-dried, wet-sieved at 63 μm, and dried in an oven at 50°C. The >63-μm fraction was sieved into four size fractions: 63 to 150, 150 to 250, 250 to 355, and >355 μm. Approximately 8 ± 3 G. ruber individuals (white, sensu stricto) were picked from the 250- to 355-μm fraction for stable isotope analysis and, where sufficient numbers were present, 35 ± 10 individuals for Mg/Ca analysis. Samples were analyzed on an MAT252 isotope ratio mass spectrometer (IRMS) coupled to a Kiel III carbonate device. Samples were reacted by individual acid addition (99% H3PO4 at 70°C). A total of 811 G. ruber samples and a total of 107 standards [Brown Yule Marble (BYM) and Carrara] were analyzed. Repeated analyses of BYM (n = 43, 1σ) yield −2.25 ± 0.03 for δ13C and −6.46 ± 0.08 for δ18O. Repeated analyses of Carrara Marble (n = 64, 1σ) yield 2.05 ± 0.02 for δ13C and −1.88 ± 0.04 for δ18O. Replicate analysis of foraminifera samples (n = 83, 1σ) yields a median half range of 0.07 for δ18O and 0.11 for δ13C. All results were calibrated to National Institute of Standards and Technology (Gaithersburg, MD) carbonate isotope standard NBS19 and are reported as ‰ Vienna Pee Dee belemnite (VPDB).

Stable benthic isotopes

Individual Uvigerina spp. and Cibicidoides wuellerstorfi (depending on availability) were picked from the 250 to 355 μm and >355 μm (when necessary) fractions. Samples were analyzed on an MAT252 IRMS coupled to a Kiel III carbonate device, reacted by individual acid addition (99% H3PO4 at 70°C). A total of 1323 samples and 180 standards (BYM and Carrara) were analyzed. Repeated analyses of BYM (n = 72, 1σ) yields −2.24 ± 0.03 for δ13C and −6.40 ± 0.08 for δ18O. Repeated analyses of Carrara Marble (n = 108, 1σ) yield 2.06 ± 0.02 for δ13C and −1.85 ± 0.05 for δ18O. Replicate analysis of foraminifera samples (n = 42, 1σ) yields a median half range of 0.03 for δ18O and 0.05 for δ13C. Isotopic offsets between C. wuellerstorfi and Uvigerina spp. were accounted for by subtracting 0.61 and adding 0.77 to Uvigerina spp. values for δ18O and δ13C, respectively, based on median offsets of paired samples (n = 452). All results were calibrated to National Institute of Standards and Technology (Gaithersburg, MD) carbonate isotope standard NBS19 and are reported as ‰ VPDB.


Age control is established by 53 tie points, including four accelerator mass spectrometry (AMS) radiocarbon dates between 1.7 and 6.1 m, the base of the Toba ash at 15.6 mbsf (74.5 ka ago), and 48 tie points correlating structure in the benthic δ18O record to the LR04 (30) global benthic stack using AnalySeries (Fig. 2) (69). For eccentricity, obliquity, and precession, coherence is 0.98, 0.99, and 0.91, respectively, with phase offsets of less than 1 ka (less than our average 2-ka sample resolution). For AMS, specimens of upper-ocean planktonic foraminiferal species primarily consisting of G. ruber (white variety), Trilobatus sacculifer, and Globigerinoides conglobatus species were picked for radiocarbon analysis, which was performed at the Woods Hole Oceanographic Institution National Ocean Sciences Accelerator Mass Spectrometry (WHOI NOSAMS) facility. Each radiocarbon date was individually calibrated into calendar years using the R code BACON (70) using the standard marine reservoir correction and the Marine13 curve (71).


Lipids were extracted (in triplicate) from approximately 3 g of freeze-dried sediment using the Dionex Accelerated Solvent Extractor ASE-200 with dichloromethane-methanol (6:4). The extract was separated into neutral and acid fractions, followed by aminopropyl silica gel chromatography (72). The neutral fraction was separated into four fractions by column SiO2 chromatography: F1, 3 ml of hexane; F2, 3 ml of hexane-toluene (3:1); F3, 4 ml of toluene; F4, 3 ml of toluene-CH3OH (3:1). The F4 was analyzed using an Agilent 1260 high-performance liquid chromatography system coupled to 6130 quadrupole mass spectrometer following the method of Schouten et al. (73). The analysis was duplicated for all samples.

The tetraether index of 86 carbon atoms (TEX86) was calculated from the concentrations of glycerol dialkyl glycerol tetraethers (GDGT)-1, GDGT-2, GDGT-3, and a regioisomer of crenarchaeol using the following expressions (74, 75)TEX86=([GDGT2]+[GDGT3]+[crenarchaeol regioisomer])/([GDGT1]+[GDGT2]+[GDGT3]+[crenarchaeol regioisomer])TEX86H=log (TEX86)

Temperatures were calculated according to the following equations based on a global core top calibration (75)T=68.4 TEX86H+38.6 (when T > 15°C)where T = temperature (in degrees Celsius); analytical accuracy was 0.45°C. Because the contribution of Archaea other than pelagic Thaumarchaeota modifies the TEX86 temperature signal (76), we excluded the data showing an anomalous methane index (>0.23 and <0.17).

Leaf wax δ13C

The acid fraction of extracted lipids was methylated with methanol–acetyl chloride (95:5) and purified with SiO2 column chromatography. The δ13C of the methylated n-fatty acid was analyzed using an Agilent 6890 series gas chromatograph combined with a Finnigan MAT Delta Plus mass spectrometer. The δ13C values of fatty acids (FA) were obtained from the measured values of fatty acid–methyl esters by correcting for methyl carbon (−34.1‰). The reproducibility of the measurement based on repeated analyses was better than ±0.1‰. The δ13CFA was calculated by averaging the δ13C of C26, C28, C30, and C32 n-fatty acids. Values are reported in ‰ relative to VPDB.

Compound-specific radiocarbon analysis of surface sediments from the Santa Monica Basin indicates that the Δ14C values of long-chain n-fatty acids are the same as those of phytoplankton biomarkers (77). By contrast, the Δ14C values of long-chain n-alkanes are more depleted than those of phytoplankton biomarkers. This indicates that n-fatty acids are better for monitoring paleovegetation changes than n-alkanes due to their smaller continental retention time.

In the BoB, French et al. (78) estimated the continental reservoir age based on the depth profile of the radiocarbon concentration in Bengal Fan sediments. A two-endmember model indicated that 79 to 83% of the long-chain n-fatty acids were stored in continental reservoirs for an average of 1000 to 1200 calendar years (slow-cycling component), and the remainder was stored for an average of 15 years (fast-cycling component). However, the millennial-scale retention time of the slow-cycling component does not markedly delay the timing of multimillennial scale δ13CFA events because it represents a broad distribution of n-fatty acids with different ages that tend to smear each other out (78).


A total of 539 Mg/Ca analyses from Site U1446 have been performed at three different laboratories. Here, we use the subset of the samples that have paired TEX86 SST, allowing a direct comparison of δ18Osw generated from the two SST proxies. Of the 243 samples presented here, 25 samples from 0 to 7.5 m (~0 to 20 ka ago) were analyzed at Rutgers University (Rosenthal), 53 samples from 16.07 to 33.57 m [70 to 140 ka ago] were analyzed at the Open University, UK (Kerr and Anand), and 165 samples from three time intervals (18 to 70, 143 to 204, 215 to 390, and 415 to 673 ka ago) were analyzed at United States Geological Survey (USGS) St. Petersburg (Richey). The 165 University States Geological Survey (USGS) samples were from splits of the Brown University samples and required no interpolation to match the depths (age) of the TEX86H data. The 25 samples from Rutgers are from a high-resolution study to be published separately and required nominal interpolation to TEX86H depths; offset from nearest sample is 0.08 ± 0.22 ka (1σ, n = 25). The 53 Open University samples also required minimal interpolation; offset from nearest sample is 0.58 ± 1.1 ka (1σ, n = 53).

Typical samples contained 6 to 50 G. ruber sensu stricto individuals (250- to 355-μm fraction). Samples analyzed at Rutgers and the Open University went through rigorous clay removal and a full-trace element cleaning (reductive + oxidative) method (79, 80). Samples analyzed at the USGS laboratory were subjected only to oxidative cleaning without the added reductive cleaning step. Samples were dissolved in diluted HNO3 to achieve constant [Ca] concentration to minimize matrix effects. Samples were analyzed using an Agilent Technologies Triple-Quad inductively coupled plasma–mass spectrometer (ICP-MS) 8800 and Element XR Sector Field ICP-MS at the Open University and Rutgers, respectively. Element to Ca ratios were calculated using a calibration derived from matrix-matched synthetic standards with varying trace element concentrations (81). At USGS, samples were analyzed on a PerkinElmer 7300 DV inductively coupled optical emission spectrometer. Mg/Ca was corrected for instrumental drift using the internal gravimetric standard method devised by Schrag (82). Reproducibility of Mg/Ca measurements, as determined by repeated analysis of synthetic CaCO3 standards, was better than 1.8% (1 σ) in all three laboratories. Contaminant ratios (Al/Ca and Fe/Ca) were monitored, and samples with Al/Ca (>300 μmol/mol) and Fe/Ca (>10 μmol/mol) were discarded if their Mg/Ca ratios were significantly different from the samples next to them.

Following previous studies (79, 83, 84), Mg/Ca data obtained with the full cleaning method were corrected by +6%, to be consistent with samples undergone only oxidative cleaning including the calibration studies. No regional Mg/Ca calibration exists for the BoB. We estimated SST from the Mg/Ca data using the calibration of Tierney et al. (85) and the Paleo-Seawater Uncertainty Solver (PSU Solver) code of Thirumalai et al. (86) that accounts for both changes in temperature and salinity. This calibration uses all available culture data in a multivariate calibration that accounts for both salinity and temperature.

X-ray fluorescence

XRF bulk geochemistry was measured using Itrax XRF microscanners at WHOI and University of Massachusetts-Amherst. Elemental counts were recorded for 10 s at either 5- or 10-mm intervals depending on time availability using molybdenum x-ray tube sources. To account for different core conditions and x-ray tube ages at the time of measurement, raw elemental counts (per second) for each scanner were unity-normalized. Here, we report Rb and Ca to construct the Rb/Ca ratio as a proxy for terrestrial runoff and upper ocean stratification. Rb is associated with fine-grained terrigenous sediment discharge increasing with runoff [e.g., (87)], whereas Ca represents changes in biogenic CaCO3 that decreases with runoff as freshwater-induced stratification in the BoB leads to reduced primary production (37).

Seawater δ18O

δ18Osw is derived using PSU Solver, removing both the global sea level and local temperature components (86). SST is reconstructed using TEX86H and Mg/Ca. On the basis of paired analyses, we subtracted 2°C from the TEX86H SST values to match those of the Mg/Ca SST values, the result being that both core-top values match modern mean annual SST at U1446 (29.7°C) (88). The 2°C offset is consistent with the unique barrier-layer structure of the BoB, characterized by an inversion whereby temperatures at or above the thermocline are up to 3°C warmer than surface waters (8991). We use the seawater δ18O temperature relationship (low light) of Bemis et al. (92). The δ18O-sea surface salinity (SSS) relationship (δ18Osw = 0.36 × SSS – 12.02) is derived from surface measurements within the EICC (31, 32). While there are a number of these regressions for the BoB (32), they are all linear, serving only to shift the resultant δ18Osw mean; hence, they have no impact on the timing of maxima and minima and therefore no impact on the orbital-band phase, the focus of this work. Quantitative reconstruction of salinity is not attempted, given the unconstrained nature of the δ18Osw-salinity relationship over time, although modeling suggests that this relationship is relatively constant in the Indian Ocean over both spatial and temporal scales (93). We use, indirectly, the global sea level curve of Waelbroeck et al. (94). Instead of using the Waelbroeck curve directly, the U1446 benthic δ18O record was smoothed (three-point Gaussian) and scaled to match the Waelbroeck curve such that the age model and sample resolution at U1446 remains intact; this is important as it prevents introduction of spurious structure into the seawater δ18O record due to differences in age models and is necessary because the Waelbroeck curve does not cover the entire 1.46-Ma length of our record. Scaling is accomplished by subtracting the core top benthic δ18O value (2.33‰) from the heaviest Last Glacial Maximum (LGM) value (3.64‰) to arrive at an offset value (1.3‰). This offset value is divided by the Waelbroeck LGM to present sea level difference (−123 m) to yield the conversion factor (−0.01068‰/m). For each U1446 δ18O benthic sample, the core top value is subtracted, and the result is divided by the conversion factor to yield the sea level curve input to PSU Solver. Propagated uncertainty in the δ18Osw result (fig. S1d) is assessed using Monte Carlo (n = 10,000) and the following 2σ error estimates on the underlying parameters (δ18Oplanktonic ± 0.12‰, Mg/Ca ± 0.2 mmol/mol, age ± 4 ka, sea level ± 5 m, and TEX86H SST ± 2°C). We test the impact of using TEX86H SST versus Mg/Ca SST by comparing reconstructed δ18Osw for samples where both exist (fig. S1e). The similarity of the resulting curves supports the use of the TEX86H SST record in reconstructing δ18Osw. We further tested the temperature sensitivity by comparing δ18Osw that was reconstructed using the TEX86H SST and the Shakun et al. (95) global SST stack (fig. S1f). The Shakun stack was sampled at U1446 resolution and adjusted to match the core top SST. The two δ18Osw reconstructions are coherent above the 95% confidence interval and in phase. Last, we note that the U1446 precession band phase of δ18Osw that was reconstructed using TEX86H SST is the same as the one found in the Andaman Sea using paired G. ruber Mg/Ca SST estimates.

Time series analysis

Cross-spectral analyses were performed with the Blackman-Tukey approach using AnalySeries software (69). All spectra used a Bartlett window and a 50% lag after interpolation to a constant time step of 2 ka, very near that of the true data resolution. Individual linear spectra were calculated using the AnalySeries periodogram function and a Bartlett window. The periodogram produces an unsmoothed spectrum (equivalent to the Blackman-Tukey with a 100% lag). The coherence of spectral peaks relative to orbital parameters and interproxy comparisons, are given in tables S1 and S3.


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: This research used samples provided by the International Ocean Discovery Program Expedition 353. Funding: S.C.C. and S.M.M. were supported by U.S. NSF OCE1634774. M.Y. was funded by JSPS grants JPMXS05R2900001 and 19H05595 and JAMSTEC Exp. 353 postcruise study. K.N.-K. and P.A. were supported by UK-IODP, Open University, and NERC (NE/L002493/1), K.T. was supported by the Technology and Research Initiative Fund, Arizona Board of Regents. Author contributions: S.C.C. conceptualized and designed the study. S.C.C. generated foraminifera δ18O data and age-depth model. M.Y. generated δ13CFA and TEX86H SST data. K.T. produced AMS dates and modern climatological analyses. L.G. generated XRF data. J.N.R., K.N.-K., Y.R., P.A., and K.T. generated Mg/Ca SST data. S.C.C. wrote the manuscript with reviews, edits, and significant intellectual contributions from all authors. The manuscript was improved by two anonymous reviewers and review by J. Rodysill of the U.S. Geological Survey. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. government. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Data necessary to understand and evaluate the conclusions of this paper are archived at Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article