Greenhouse gases modulate the strength of millennial-scale subtropical rainfall, consistent with future predictions

Millennial-scale East Asian monsoon variability is closely associated with natural hazards through long-term variability in flood and drought cycles. Therefore, exploring what drives the millennial-scale variability is of significant importance for future prediction of extreme climates. Here we present a new East Asian summer monsoon (EASM) rainfall reconstruction from the northwest Chinese Loess Plateau (CLP) spanning the past 650 kyr. The magnitude of millennial-scale variability (MMV) in EASM rainfall is linked to ice volume and greenhouse gas (GHG) at the 100 kyr eccentricity band and to GHG and summer insolation at the precession band. At the glacial–interglacial timescale, gradual changes in CO₂ during intermediate glaciations lead to increased variability in North Atlantic stratification and Atlantic meridional overturning circulation, propagating abrupt climate changes into East Asia via the westerlies. Within the 100 kyr cycle, precession variability further enhances the response, showing that stronger insolation and increased atmospheric GHG cause increases in the MMV of EASM rainfall. These findings indicate increased extreme precipitation events under future warming scenarios, consistent with model results.

1 Introduction

Chinese loess is a unique terrestrial archive that can document East Asian monsoon (EAM) variability well at tectonic to millennial timescales (Porter and An, 1995; An, 2000). High-resolution loess records have revealed persistent millennial-scale (1–10 kyr periodicity) EAM fluctuations spanning the last several glacial cycles (Guo et al., 1996, 2021a; Sun et al., 2012, 2021a, b), which are dynamically linked with high-latitude abrupt changes in the north Atlantic including Heinrich (H) (Heinrich, 1988) and Dansgaard–Oeschger (DO) events (Dansgaard et al., 1993). This millennial-scale monsoon variability is superimposed on glacial–interglacial variations (Ding et al., 1999; Clemens et al., 2018). Abrupt summer monsoon changes are closely linked to natural hazards such as flood and drought events (Huang et al., 2007), since the summer monsoon plays a leading role in transporting water vapor from low-latitude oceans to middle- and high-latitude continents of the Northern Hemisphere (Liu et al., 2013; An et al., 2015). Abrupt rainfall events associated with short-term summer monsoon variations strongly influence agriculture, food production, water supply and social economic development (Huang et al., 2007; Cook et al., 2010; Li et al., 2017). However, how these flood and drought events are affected by both natural and anthropogenic factors remains poorly constrained. Understanding the mechanisms that modulate the magnitude of millennial-scale variability (MMV) is of critical importance for the scientific community as well as policy makers. Here we use the term “modulate” in the context of lower-frequency components of the climate system influencing or determining the amplitude of a higher-frequency component.

A number of well-dated, high-resolution speleothem δ¹⁸O records have been developed in recent years (Wang et al., 2008; Cheng et al., 2016), providing the opportunity to examine the underlying relationship(s) between the MMV and potential longer-term (orbital-scale) modulators. Cheng et al. (2016) hypothesized, on the basis of an East Asian composite speleothem δ¹⁸O record (δ¹⁸O_sp), that periods of maximum Northern Hemisphere summer insolation correspond to weaker millennial-scale variability. Subsequently, however, Thirumalai et al. (2020) showed that precession does not modulate the MMV of δ¹⁸O_sp and postulated that it is, instead, modulated by internal processes related to the cryosphere. This work also raised the possibility that δ¹⁸O_sp is decoupled from regional Asian monsoon rainfall over millennial timescales (Zhang et al., 2018). The latest research suggests that the MMV through the Pleistocene is influenced by both glacial boundary conditions and orbital configurations (Sun et al., 2021b). As such, two important outstanding questions remain: is there a reliable proxy for East Asian summer monsoon (EASM) rainfall at the millennial timescale, and what factors modulate the MMV thereof?

Due to weak pedogenesis and high sedimentation rates, millennial-scale oscillations are well preserved in the northwestern Chinese Loess Plateau (CLP) over the past glacial cycles (Sun et al., 2021a). The Linxia profile is well-suited for reconstructing rapid monsoon changes because it is located in the monsoon frontal zone and sensitive to high- and low-latitude climate variability (Guo et al., 2021a). To address the above questions, we have generated a high-resolution summer monsoon proxy ($\mathrm{Ca}/\mathrm{Ti}$) from Linxia on the western CLP (Fig. 1). The $\mathrm{Ca}/\mathrm{Ti}$ ratio is a precipitation-sensitive proxy linked to summer monsoon rainfall (Guo et al., 2021a). Low values of $\mathrm{Ca}/\mathrm{Ti}$ indicate stronger Ca leaching associated with intensified summer rainfall. The new precipitation proxy ($\mathrm{Ca}/\mathrm{Ti}$) and δ¹⁸O_sp are evaluated to elucidate the modulating drivers of these two proxy records. As discussed in the Results section, we find that the MMV of $\mathrm{Ca}/\mathrm{Ti}$ is mainly modulated by ice volume and greenhouse gases (GHG) at the eccentricity band. Both GHG and summer insolation modulate the MMV of $\mathrm{Ca}/\mathrm{Ti}$ at the precession band but not that of δ¹⁸O_sp; δ¹⁸O_sp MMV is modulated by winter insolation at the eccentricity and obliquity bands. The interpretations of these results are presented in the Discussion section.

Figure 1 The location of the Linxia loess profile and Sanbao–Hulu cave records. The Linxia profile, located on the edge of convergence zone of the alpine Qinghai–Tibet Plateau, with the northwest being arid and the southeast being monsoon area, is very sensitive to the migration of desert regions and monsoonal rainfall. Sanbao–Hulu cave is located in the monsoon-influenced Yangtze River valley, sensitive to monsoon-induced precipitation changes. The black dashed line represents the scope of Chinese Loess Plateau. The base map is drawn using GMT software, and the elevation data are from http://www.ngdc.noaa.gov/mgg/global/global.html (last access: 14 June 2022).

2 Materials and methods

2.1 Sampling and measurements

The Linxia (LX; 35.15^∘ N, 103.63^∘ E, 2200 m a.s.l.) loess record is located in the western edge of the CLP (Fig. 1). At present, mean annual temperature and precipitation in this region are about 8.1 ^∘C and 484 mm, respectively, with ∼80 % of the annual precipitation falling during the summer season (May to September). The 203.8 m long core A (LXA) consists of 185 m of eolian loess–paleosol sequences, underlain by 17 m of fluvial loess and 1.8 m of sandy gravel layers. The 72 m long core B (LXB) and a 7 m pit were excavated in 2017. Powder samples were collected at 2 cm intervals for analyzing mean grain size (MGS) with the resolution ranging from 10 to 200 yr cm⁻¹. Also, 2 cm resolution samples were dried at 40 ^∘C overnight and ground to 200 mesh size (about <75 µm) with an agate mortar and pestle, and then they were pressed into a plastic sheet (4 cm × 4 cm × 0.3 cm), creating a flat and homogeneous slide. The plastic slides were then placed on a wood pallet for XRF scanning to obtain elemental intensities (Guo et al., 2021a).

2.2 Age model and evaluations of age uncertainties

The chronological framework for Chinese loess–paleosol sequences is widely constructed by matching grain size to the benthic LR04 δ¹⁸O stack (Porter and An, 1995; Hao et al., 2012; Sun et al., 2021a). A similar approach utilizes correlation to speleothem δ¹⁸O in recent years (Beck et al., 2018; Sun et al., 2021a; Zhang et al., 2022). The age model of the Linxia profile is generated by synchronizing Chinese loess and speleothem δ¹⁸O records back to 650 ka (Sun et al., 2021a). Detailed information of loess–paleosol boundaries is shown in Fig. 2. The first set of control points tie the loess–paleosol boundaries S₆ to S₀ to the timing of the glacial terminations and inceptions in speleothem δ¹⁸O (Cheng et al., 2016). The second and third sets of age control points tie the timing of precessional transition boundaries and abrupt cooling events in the MGS record to those in speleothem δ¹⁸O (Fig. 2b), based on the assumption that the East Asian summer and winter monsoons co-vary at orbital timescales, and millennial-scale abrupt events are synchronous in the Northern Hemisphere (Hemming et al., 2004; Sun et al., 2012; Clemens et al., 2018). The tie points are shown in Fig. 2b.

Figure 2 (a) Strata and down-core variations in mean grain size (MGS), magnetic susceptibility (MS), $\mathrm{Ca}/\mathrm{Ti}$ and sedimentation rate against depth (brown in benthic δ¹⁸O age model and dark brown in speleothem δ¹⁸O age model). Brown red, orange and yellow rectangles represent paleosol layers, weakly pedogenic paleosols and loess layers, respectively. The timing of dashed lines and glacial–interglacial transitions is control points of benthic δ¹⁸O chronology. (b) Variations in MGS, $\mathrm{Ca}/\mathrm{Ti}$ over last 650 kyr and age model of the Linxia loess section. Comparison of MGS and $\mathrm{Ca}/\mathrm{Ti}$ in the Linxia section with Sanbao–Hulu (Cheng et al., 2016) and benthic δ¹⁸O stack (Lisiecki and Raymo, 2005). The dark brown squares, blue triangles and red dots represent the first (glacial–interglacial transition), second (precession cycles) and third (millennial-scale events) class age control points at the corresponding position of the cave record, respectively (Sun et al., 2021a). Light blue bands denote the interglacial times. The short green rectangles represent the age differences between the two age models.

The composite speleothem δ¹⁸O record is well resolved by absolute U–Th dating and applied as the target for regional synchronization. The errors are less than 2 kyr for the last 450 kyr and increase to 4–8 kyr before 450 ka (Cheng et al., 2016). The MGS record yields good correlation between loess–paleosol boundaries and glacial–interglacial transitions of LR04 δ¹⁸O. The age differences of most glacial terminations are around 2–4 kyr (Fig. 2b, Sun et al., 2021a). The speleothem δ¹⁸O synchronized age model is compared with the benthic δ¹⁸O age model to evaluate the influence of age uncertainties on the wavelet coherence analysis. The small differences in the two age models (correlation to marine δ¹⁸O and to speleothem δ¹⁸O) make little difference in the MMV and associated wavelet coherence (WTC). This is because tie points are separated by 20–30 kyr, and the small differences in the tie points among the two age models have insignificant influence on the amplitude (MMV) of millennial-scale peaks. Only minor differences in the MMV WTC phase are observed at the obliquity (450–550 ka) and precession bands (100–200 ka; compare Fig. S1 in the Supplement with Fig. 4).

2.3 Spectrum and wavelet coherence analyses

In order to estimate the MMV, loess $\mathrm{Ca}/\mathrm{Ti}$ and δ¹⁸O_sp are linearly interpolated at 0.1 kyr intervals. The WTC results remain unchanged unless the cutoff threshold is reduced to 6 kyr or increased to 12 kyr; then original time series are filtered using a Butterworth filter at a cutoff threshold of 10 kyr (e.g., $\mathrm{Ca}/\mathrm{Ti}$-hi-10 kyr). The moving standard deviation of millennial-scale variability is calculated to ascertain the orbitally related modulation and its association with internal and external forcing using a 2 kyr sliding window (calculation method follows Thirumalai et al., 2020). The spectra of all proxies were calculated using the Lomb–Scargle periodogram (https://exoplanetarchive.ipac.caltech.edu/cgi-bin/Pgram/nph-pgram, last access: 10 April 2022), which has the advantage of analyzing discontinuous time series and removal of spurious spectral characteristics (VanderPlas, 2018).

The normalized and combined orbital parameters eccentricity, tilt and negative precession (ETP); GHG; insolation; and benthic δ¹⁸O were evaluated by WTC to extract maximal phase and amplitude correlations with astronomical, ice volume and greenhouse gas forcing over the past 650 kyr. WTC between time series was performed in a Monte Carlo framework (n=1000) following Grinsted et al. (2004). The WTC evaluates the time evolution of coherence between two time series across a range of frequency bands (independent of the power). The black arrows in the figures represent the phase relationship between the two time sequences with rightward, upward and downward arrows indicating in phase, leading phase and lagging phase, respectively. The color scale indicates the amplitude correlations between the two datasets.

In this paper, the parameter ΔRF_GHG is a measure of GHG radiative forcing. ΔRF_GHG is reconstructed by referencing the content of EPICA ice core greenhouse gases to the modern value. ΔRF_GHG is defined as the difference between a certain past GHG level ([CO₂] and [CH₄]) and the pre-industrial greenhouse gas level ([CO₂]₀=280 ppm, [CH₄]₀=700 ppb) (Ramaswamy et al., 1991). While CH₄ contributes <5 %, we calculated the ΔRF_GHG using both CO₂ and CH₄. The equation used to determine ΔRF_GHG is as follows (Lo et al., 2017).

$$\begin{array}{rl}\mathrm{\Delta }{\mathrm{RF}}_{\mathrm{GHG}}& =\mathrm{\Delta }{\mathrm{RF}}_{{\mathrm{CO}}_{\mathrm{2}}}+\mathrm{\Delta }{\mathrm{RF}}_{{\mathrm{CH}}_{\mathrm{4}}}\\ & =\mathrm{4.841}\ln \left(\left[{\mathrm{CO}}_{\mathrm{2}}\right]/\left[{\mathrm{CO}}_{\mathrm{2}}\right]{}_{\mathrm{0}}\right)\\ & +\mathrm{0.0906}\left(\sqrt{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}-\sqrt{\left[{\mathrm{CO}}_{\mathrm{2}}\right]{}_{\mathrm{0}}}\right)\\ & +\mathrm{0.036}\ln \left(\sqrt{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}-\sqrt{\left[{\mathrm{CH}}_{\mathrm{4}}\right]{}_{\mathrm{0}}}\right)\end{array}$$

3 Results

The MGS reflects grain-size sorting and is very sensitive to winter monsoon variations (Porter and An, 1995; Sun et al., 2006), with coarser particle size during the glacials. The $\mathrm{Ca}/\mathrm{Ti}$ ratio reflects precipitation-induced leaching intensity linked to summer monsoon rainfall (Guo et al., 2021a), with relatively low values during the interglacials. The high-resolution δ¹⁸O of the Sanbao–Hulu speleothem is an indicator of EASM changes at orbital to millennial timescales (Cheng et al., 2016). The MGS and $\mathrm{Ca}/\mathrm{Ti}$ exhibit distinct glacial–interglacial and precessional variations over the last 650 kyr as seen in LR04 δ¹⁸O (Lisiecki and Raymo, 2005) and speleothem δ¹⁸O (Cheng et al., 2016), respectively (Fig. 2b).

Both $\mathrm{Ca}/\mathrm{Ti}$ and δ¹⁸O_sp records show clear millennial-scale fluctuations overlaying orbital-scale variations. The high-frequency millennial signals were persistent over the last 650 kyr, but the amplitude varies from proxy to proxy (Figs. 3a and S2a). Spectral analysis of the raw records and MMV for loess and speleothem records displays variable associations with eccentricity (∼100 kyr), obliquity (∼41 kyr) and precession (∼23 and ∼19 kyr) scales over the past 650 kyr. Loess $\mathrm{Ca}/\mathrm{Ti}$ variance is mainly concentrated in obliquity, with lesser variance in the eccentricity and precession bands (Fig. 3b), indicating prominent ice volume (eccentricity and obliquity) and isolation (precession) forcing. The speleothem δ¹⁸O shows predominant precession-scale variance (Fig. S2b), suggesting strong links to insolation forcing (Cheng et al., 2016). These results indicate ice volume and insolation play dominant roles in driving changes in loess $\mathrm{Ca}/\mathrm{Ti}$ and speleothem δ¹⁸O, respectively (Cheng et al., 2016; Sun et al., 2021a).

Figure 3 Raw datasets, millennial-scale components (10 kyr high-pass-filtering signals) and MMV of the Linxia loess $\mathrm{Ca}/\mathrm{Ti}$ record over the past 650 kyr with their corresponding spectra. The orbital bands are marked with red dashed lines (eccentricity: 100 kyr, obliquity: 41 kyr, precession: 23 and 19 kyr). Clear eccentricity, obliquity and precession variance as well as persistent millennial-scale components are observed for loess $\mathrm{Ca}/\mathrm{Ti}$ and MMV.

Millennial-scale fluctuations co-exist with long-term orbital- and ice-volume variability. We seek to assess the potential linkages among them and, in particular, the extent to which MMV is modulated by these longer-term orbital and internal climate parameters. The spectrum of $\mathrm{Ca}/\mathrm{Ti}$ MMV shows dominant eccentricity with less strong precession and weak obliquity variance (Fig. 3d). The spectrum of δ¹⁸O_sp MMV has a small peak near 100 kyr and an offset 41 kyr peak with little to no variance at the 23 kyr period (Fig. S2d). Thus, while both proxies are similarly modulated at the 100 kyr period (i.e., larger MMV during glacial intervals relative to interglacial times), the MMV modulation is variable for the two proxies at other orbital bands.

As with the spectral differences in the raw records, the MMV spectra also imply different MMV modulating drivers, potentially associated with insolation, ice volume and/or GHG (Thirumalai et al., 2020). How do internal and external drivers interact with each other and modulate the MMV of these records at the orbital timescale? We performed wavelet coherence and phase analyses of both MMV records relative to ETP, ice volume, ΔRF_GHG, summer insolation and winter insolation to identify which variables might modulate the MMV of these EASM records. The MMV in $\mathrm{Ca}/\mathrm{Ti}$ is strongly coherent with ice volume and GHG at the 100 kyr eccentricity band and with GHG and summer insolation at the 23 kyr precession band (Fig. 4c, d and g). δ¹⁸O_sp MMV is most strongly coherent with GHG and ice volume at the 100 kyr band and with winter insolation at the eccentricity and obliquity bands (Fig. S3c, d and g).

Figure 4 Comparison of (a) 35^∘ N summer insolation forcing, (b) GHG radiative forcing (black dashed line denotes the precession band-pass filtering results of ΔRF_GHG) and (e) ice volume and (f) ETP for MMV of Linxia loess $\mathrm{Ca}/\mathrm{Ti}$. Wavelet coherence between (c) 35^∘ N summer insolation, (d) GHG radiative forcing, (g) ice volume, (h) ETP and MMV of loess $\mathrm{Ca}/\mathrm{Ti}$ over the past 650 kyr. The orbital bands are marked with red dashed lines. The orange color indicates strong correlation for the two time series. The black lines plot coefficients of determination greater than 0.76. The black arrows represent the phase relationship, with rightward, upward and downward arrows indicating in phase, leading phase and lagging phase, respectively. Strong eccentricity, weak obliquity and precession band ice volume modulation are observed for MMV of loess $\mathrm{Ca}/\mathrm{Ti}$. Strong eccentricity and precession band GHG modulation as well as weak summer insolation forcing are detected for MMV of loess $\mathrm{Ca}/\mathrm{Ti}$.

4 Discussion

4.1 Orbital-scale modulation factors for MMV of the EASM

Geological records and modeling results indicate that high-latitude ice volume and ice sheet topography play important roles in triggering abrupt climate changes (Broecker et al., 1994; Clark et al., 2001). In particular, abrupt climate changes are highly sensitive to ice volume variations; ice sheets are widely hypothesized to motivate and amplify these high-frequency signals within a constrained benthic oxygen isotope “ice volume threshold” between 3.5 ‰ and 4.5 ‰ (McManus et al., 1999; Bailey et al., 2010; Naffs et al., 2013; Zhang et al., 2014). Wavelet coherence between the MMV of loess $\mathrm{Ca}/\mathrm{Ti}$, speleothem δ¹⁸O and global benthic δ¹⁸O stack shows excellent coherence and near-zero phase with ice volume at the 100 kyr band (Figs. 4e, g and S3e, g). This in-phase variation demonstrates that EASM MMV primarily follows the glacial–interglacial rhythm of ice volume variations, enlarged during glacial times and dampened during interglacial times. However, coherence of the MMV for these two proxies with the benthic δ¹⁸O stack is relatively weak and variable at the 41 kyr band (δ¹⁸O_sp; Fig. S3e and g) and 23 kyr band ($\mathrm{Ca}/\mathrm{Ti}$; Fig. 4e and g). These relationships demonstrate that ice volume directly modulates the MMV of the EASM, predominantly at the 100 kyr band, with high ice volume corresponding to larger MMV.

GHG concentration is another potential driver of abrupt climate changes (Alvarez-Solas et al., 2011; Zhang et al., 2017). Wavelet results between the MMV of loess $\mathrm{Ca}/\mathrm{Ti}$, speleothem δ¹⁸O and the record of GHG RF show excellent coherence and $\sim \mathrm{180}{}^{\circ }$ phase at the 100 kyr eccentricity band (Fig. 4b, d and S3b, d), indicating strong MMV at times of low GHG. Given the coupled nature of global ice volume and atmospheric GHG, it is suggested that during the late Pleistocene glacial–interglacial cycles, these two factors can amplify the MMV of the EASM as recorded by $\mathrm{Ca}/\mathrm{Ti}$ and speleothem δ¹⁸O during times of high ice volume and low GHG concentration. However, this is not the case for the precession band. MMV of loess $\mathrm{Ca}/\mathrm{Ti}$ displays discrete intervals of high coherence and near-zero phase with GHG RF at the precession band (Fig. 4b and d), which is not the case for speleothem δ¹⁸O (Fig. S3b and d). Thus, GHG RF does play a role in modulating $\mathrm{Ca}/\mathrm{Ti}$ MMV but not that of δ¹⁸O_sp at the precession band, indicating a difference in the millennial-scale response of these two proxies. We investigate this further by assessing the response to local insolation forcing.

The MMV of $\mathrm{Ca}/\mathrm{Ti}$ shows discontinuous relatively weak coherence with 35^∘ N summer insolation at the precession band with even weaker coherence at the 41 kyr band (Fig. 4a and c) and relatively strong modulation of the summer insolation compared to that of GHG at the precession band (Fig. 4b and d). In contrast, the MMV of δ¹⁸O_sp displays high coherence and zero phase with 35^∘ N winter insolation at the 100 kyr period; relatively weaker coherence, with a lagging phase, at the 41 kyr band; and negligible coherence at the 23 kyr band (Fig. S3a and c). These results indicate that the MMV of speleothem δ¹⁸O is modulated by local winter insolation, opposite to the hypothesis calling on Northern Hemisphere summer insolation (Cheng et al., 2016).

4.2 Mechanisms and implication for modulation of EASM MMV

At the glacial–interglacial timescale, the MMV is amplified under the glacial boundary conditions. This indicates dynamic linkages with high-latitude North Atlantic Heinrich and DO events (Cheng et al., 2016; Sun et al., 2012, 2021a, b). Heinrich and DO variability are linked to Northern Hemisphere ice sheet (NHIS) perturbations via its influence on freshwater flux into the North Atlantic Ocean and consequent Atlantic meridional overturning circulation (AMOC) changes (McManus et al., 1999; Hemming, 2004; Naffs et al., 2013). At times of intermediate ice sheet volume, minor changes in NHIS height and atmospheric CO₂ concentrations can trigger rapid climate transitions (Zhang et al., 2014, 2017). Altering the height of NHIS leads to changes in the gyre circulation and sea-ice coverage by shifting the northern westerlies (Zhang et al., 2014). The maximum westerly wind stress shifts northwards associated with a gradual increase in the Northern Hemisphere ice volume. The northward westerly, in turn, encourages the EASM rain belt to move northward (He et al., 2021) and results in increases in the MMV of EASM rainfall (especially northern China).

In addition, CO₂ acts as an internal feedback agent to AMOC changes (Barker et al., 2016). Under intermediate glacial conditions, when the AMOC reaches a regime of bi-stability, rising CO₂ during Heinrich stadial cold events can trigger abrupt transitions to warm conditions. Decreasing CO₂ during warm events leads to abrupt cooling transitions (Zhang et al., 2017). Therefore, CO₂ generally provides a negative feedback on MMV of EASM rainfall. During interglacial times, decreasing ice volume, accompanied by reduced sea ice and more frequent freshwater perturbation, is correlated with lower-frequency and smaller-amplitude variability in abrupt climate events. The co-evolving GHG concentrations would further alter the sea surface temperature by greenhouse forcing, subsequently modulating the MMV.

Within the 100 kyr cycle, precession-band variability (four to five cycles), characterized by increased insolation and atmospheric GHG, further heightens the positive response, leading to larger MMV of subtropical rainfall. Recent transient sensitivity experiments suggest that millennial-scale rainfall variability is driven primarily by meltwater and secondarily by insolation (He et al., 2021). During interglacial times under the combined influence of insolation and CO₂, model simulation shows that when insolation reaches the lower “threshold” value (between 358.2 and 352.1 W m⁻²), it triggers a strong abrupt weakening of the AMOC and results in abrupt cooling transitions over the last 800 kyr (Yin et al., 2021). Increased insolation could warm sea surface temperature and accelerate freshwater input from high-latitude ice sheets as well as alter GHG concentration in the atmosphere (Lewkowicz and Way, 2019), which could, in turn, modulate MMV changes in the low-latitude monsoon regions.

If both millennial-scale $\mathrm{Ca}/\mathrm{Ti}$ and δ¹⁸O_sp represent subtropical rainfall amount, the modulation factors should be consistent. However, MMV modulators at eccentricity, obliquity and precession bands differ for loess $\mathrm{Ca}/\mathrm{Ti}$ and δ¹⁸O_sp, indicating they monitor different aspects of millennial-scale monsoon circulations. Modern observations and Lagrangian trajectories of air parcels in China during the summer monsoon indicate that moisture-induced precipitation does not derive from the strongest water vapor pathways (Sun et al., 2011; Jiang et al., 2017); local water vapor recycling contributes significantly to regional precipitation in east China (over 30 %) and north China (exceeding 55 %, Shi et al., 2020). Hence, we speculate that δ¹⁸O_sp MMV monitors changes in the isotopic composition of rainfall, varying with changes in westerly transport paths associated with North Atlantic cooling events, consistent with the MMV of δ¹⁸O_sp being closely linked to winter insolation at 100 and 41 kyr periods and the absence of MMV modulation at the precession band. We further hypothesize that $\mathrm{Ca}/\mathrm{Ti}$ mainly represents the MMV in local rainfall amount, consistent with the MMV of tropical rainfall being more dynamically related to GHG and summer insolation at the precession band.

In recent decades, atmospheric GHG concentration has been accelerating due to anthropogenic contribution of fossil fuels, suggesting that extreme EASM-induced precipitation will increase as well. This inference is consistent with model simulations indicating that the number of extreme daily precipitation events and mean precipitation overall will increase significantly in response to higher GHG concentration (Dairaku and Emori, 2006). The anthropogenic GHG-evoked warming is projected to increase the lower-tropospheric water vapor content and enhance the thermal contrast between land and ocean (Kitoh et al., 1997). This will give rise to a northward shift of lower tropospheric monsoon circulation and an increase in rainfall during the East Asian summer monsoon (Vecchi and Soden, 2007). Our results indicate that factors modulating EASM precipitation MMV in the past are consistent with those predicted to influence future changes in monsoonal precipitation, lending further confidence in those projections.

5 Conclusions

Our high-resolution loess $\mathrm{Ca}/\mathrm{Ti}$ record displays millennial monsoon oscillations that persist over the last 650 kyr. The MMV of loess $\mathrm{Ca}/\mathrm{Ti}$ and speleothem δ¹⁸O is modulated by different orbital factors, implying that these two proxies document different climatic response of millennial-scale monsoon circulation. The inferred mechanisms of how these internal and external factors modulate the MMV call on dynamic linkages to variability in AMOC at both eccentricity and precession bands. In recent decades, atmospheric GHG concentration has been dramatically increasing due to anthropogenic contribution of fossil fuels (Bousquet et al., 2006; Davis et al., 2010), resulting in accelerated melting of ice sheets in bi-polar regions (Swingedouw et al., 2008; Golledge et al., 2019). Their combined effects lead to more frequent occurrences of extreme rainfall (Dairaku and Emori, 2006; Hoegh-Guldberg et al., 2018). Our results indicate that the MMV EASM rainfall is modulated by ice volume, GHG and insolation factors, consistent with those predicted to influence future changes in monsoonal precipitation.