Articles | Volume 18, issue 7
Clim. Past, 18, 1675–1684, 2022
Clim. Past, 18, 1675–1684, 2022
Research article
20 Jul 2022
Research article | 20 Jul 2022

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

Greenhouse gases modulate the strength of millennial-scale subtropical rainfall, consistent with future predictions
Fei Guo1,2,3, Steven Clemens3, Yuming Liu2,4, Ting Wang2,4, Huimin Fan2, Xingxing Liu2, and Youbin Sun2,5,6 Fei Guo et al.
  • 1Institute of Marine Science and Technology, Shandong University, Qingdao 266237, China
  • 2State Key Laboratory of Loess and Quaternary Geology, Institute of Earth Environment, Chinese Academy of Sciences, Xi'an 710061, China
  • 3Department of Earth, Environmental, and Planetary Sciences, Brown University, Providence, RI 02912-1846, USA
  • 4University of Chinese Academy of Sciences, Beijing 100049, China
  • 5Institute of Global Environmental Change, Xi'an Jiaotong University, Xi'an 710061, China
  • 6Open Studio for Oceanic-Continental Climate and Environment Changes, Pilot National Laboratory for Marine Science and Technology (Qingdao), Qingdao 266200, China

Correspondence: Fei Guo ( and Steven Clemens (


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 CO2 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 δ18O 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 δ18O record (δ18Osp), 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 δ18Osp and postulated that it is, instead, modulated by internal processes related to the cryosphere. This work also raised the possibility that δ18Osp 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 (Ca/Ti) from Linxia on the western CLP (Fig. 1). The Ca/Ti ratio is a precipitation-sensitive proxy linked to summer monsoon rainfall (Guo et al., 2021a). Low values of Ca/Ti indicate stronger Ca leaching associated with intensified summer rainfall. The new precipitation proxy (Ca/Ti) and δ18Osp are evaluated to elucidate the modulating drivers of these two proxy records. As discussed in the Results section, we find that the MMV of Ca/Ti is mainly modulated by ice volume and greenhouse gases (GHG) at the eccentricity band. Both GHG and summer insolation modulate the MMV of Ca/Ti at the precession band but not that of δ18Osp; δ18Osp MMV is modulated by winter insolation at the eccentricity and obliquity bands. The interpretations of these results are presented in the Discussion section.

Figure 1The 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 (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−1. 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 δ18O stack (Porter and An, 1995; Hao et al., 2012; Sun et al., 2021a). A similar approach utilizes correlation to speleothem δ18O 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 δ18O 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 S6 to S0 to the timing of the glacial terminations and inceptions in speleothem δ18O (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 δ18O (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), Ca/Ti and sedimentation rate against depth (brown in benthic δ18O age model and dark brown in speleothem δ18O 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 δ18O chronology. (b) Variations in MGS, Ca/Ti over last 650 kyr and age model of the Linxia loess section. Comparison of MGS and Ca/Ti in the Linxia section with Sanbao–Hulu (Cheng et al., 2016) and benthic δ18O 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 δ18O 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 δ18O. The age differences of most glacial terminations are around 2–4 kyr (Fig. 2b, Sun et al., 2021a). The speleothem δ18O synchronized age model is compared with the benthic δ18O 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 δ18O and to speleothem δ18O) 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 Ca/Ti and δ18Osp 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., Ca/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 (, 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 δ18O 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 ΔRFGHG is a measure of GHG radiative forcing. ΔRFGHG is reconstructed by referencing the content of EPICA ice core greenhouse gases to the modern value. ΔRFGHG is defined as the difference between a certain past GHG level ([CO2] and [CH4]) and the pre-industrial greenhouse gas level ([CO2]0=280 ppm, [CH4]0=700 ppb) (Ramaswamy et al., 1991). While CH4 contributes <5 %, we calculated the ΔRFGHG using both CO2 and CH4. The equation used to determine ΔRFGHG is as follows (Lo et al., 2017).

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 Ca/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 δ18O of the Sanbao–Hulu speleothem is an indicator of EASM changes at orbital to millennial timescales (Cheng et al., 2016). The MGS and Ca/Ti exhibit distinct glacial–interglacial and precessional variations over the last 650 kyr as seen in LR04 δ18O (Lisiecki and Raymo, 2005) and speleothem δ18O (Cheng et al., 2016), respectively (Fig. 2b).

Both Ca/Ti and δ18Osp 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 Ca/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 δ18O 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 Ca/Ti and speleothem δ18O, respectively (Cheng et al., 2016; Sun et al., 2021a).

Figure 3Raw datasets, millennial-scale components (10 kyr high-pass-filtering signals) and MMV of the Linxia loess Ca/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 Ca/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 Ca/Ti MMV shows dominant eccentricity with less strong precession and weak obliquity variance (Fig. 3d). The spectrum of δ18Osp 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, ΔRFGHG, summer insolation and winter insolation to identify which variables might modulate the MMV of these EASM records. The MMV in Ca/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). δ18Osp 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 4Comparison of (a) 35 N summer insolation forcing, (b) GHG radiative forcing (black dashed line denotes the precession band-pass filtering results of ΔRFGHG) and (e) ice volume and (f) ETP for MMV of Linxia loess Ca/Ti. Wavelet coherence between (c) 35 N summer insolation, (d) GHG radiative forcing, (g) ice volume, (h) ETP and MMV of loess Ca/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 Ca/Ti. Strong eccentricity and precession band GHG modulation as well as weak summer insolation forcing are detected for MMV of loess Ca/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 Ca/Ti, speleothem δ18O and global benthic δ18O 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 δ18O stack is relatively weak and variable at the 41 kyr band (δ18Osp; Fig. S3e and g) and 23 kyr band (Ca/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 Ca/Ti, speleothem δ18O and the record of GHG RF show excellent coherence and 180 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 Ca/Ti and speleothem δ18O during times of high ice volume and low GHG concentration. However, this is not the case for the precession band. MMV of loess Ca/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 δ18O (Fig. S3b and d). Thus, GHG RF does play a role in modulating Ca/Ti MMV but not that of δ18Osp 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 Ca/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 δ18Osp 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 δ18O 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 CO2 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, CO2 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 CO2 during Heinrich stadial cold events can trigger abrupt transitions to warm conditions. Decreasing CO2 during warm events leads to abrupt cooling transitions (Zhang et al., 2017). Therefore, CO2 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 CO2, model simulation shows that when insolation reaches the lower “threshold” value (between 358.2 and 352.1 W m−2), 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 Ca/Ti and δ18Osp represent subtropical rainfall amount, the modulation factors should be consistent. However, MMV modulators at eccentricity, obliquity and precession bands differ for loess Ca/Ti and δ18Osp, 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 δ18Osp 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 δ18Osp 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 Ca/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 Ca/Ti record displays millennial monsoon oscillations that persist over the last 650 kyr. The MMV of loess Ca/Ti and speleothem δ18O 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.

Code availability

The wavelet coherence analyses between MMV records and ETP, GHG, ice volume and insolation are calculated using the MATLAB code. The MATLAB codes for our analyses are available at (Grinsted et al., 2004).

Data availability

Loess mean grain size and magnetic susceptibility data of this paper are available in the East Asian Paleoenvironmental Science Database (, last access: 18 July 2022, data DOI:, Sun et al., 2021a). The loess Ca / Ti is available in ESSOAr (, Guo, 2021b).


The supplement related to this article is available online at:

Author contributions

FG performed data analysis and interpretation and prepared the manuscript. SC and YS designed this study and contributed to the interpretation and writing of the manuscript. HF assisted in the data generation and data processing. TW, YL and XL contributed to the fieldwork and paper discussion.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank Xiaojing Du for offering ideas on potential model tests for this paper.

Financial support

This research has been supported by the Institute of Earth Environment, Chinese Academy of Sciences (grant no. XDB40000000) and the National Natural Science Foundation of China (grant nos. 41525008 and 41977173).

Review statement

This paper was edited by Amaelle Landais and reviewed by two anonymous referees.


Alvarez-Solas, J., Charbit, S., Ramstein, G., Paillard, D., Dumas, C., Ritz, C., and Roche, D. M.: Millennial-scale oscillations in the Southern Ocean in response to atmospheric CO2 increase, Global Planet. Change, 76, 128–136,, 2011. 

An, Z.: The history and variability of the East Asian paleomonsoon climate, Quat. Sci. Rev., 19, 171–187,, 2000. 

An, Z., Wu, G., Li, L., Li, J., Sun, Y., Liu, Y., Zhou, W., Cai, Y., Duan, A., Li, L., Mao, J., Cheng, H., Shi, Z., Tan, L., Yan, H., Ao, H., Chang, H., and Feng, J.: Global monsoon dynamics and climate change, Annu. Rev. Earth Planet., 43, 29–77,, 2015. 

Bailey, I., Bolton, C. T., DeConto, R. M., Pollard, D., Schiebel, R., and Wilson, P. A.: A low threshold for North Atlantic ice rafting from “low-slung slippery” late Pliocene ice sheets, Paleoceanography, 25, PA1212,, 2010. 

Barker, S. and Knorr, G.: A paleo-perspective on the AMOC as a tipping element, PAGES News, 24, 14–15, (last access: 12 July 2022), 2016. 

Beck, J. W., Zhou, W., Li, C., Wu, Z., White, L., Xian, F., Kong, X., and An, Z.: A 550,000-year record of East Asian monsoon rainfall from 10Be in loess, Science, 360, 877–881, (last access: 12 July 2022), 2018. 

Bousquet, P., Ciais, P., Miller, J. B., Dlugokencky, E. J., Hauglustaine, D. A., Prigent, C., and White, J.: Contribution of anthropogenic and natural sources to atmospheric methane variability, Nature, 443, 439–443,, 2006. 

Broecker, W. S.: Massive iceberg discharges as triggers for global climate change, Nature, 372, 421–424,, 1994. 

Cheng, H., Edwards, L., Sinha, A., Spötl, C., Yi, L., Chen, S., Kelly, M., Kathayat, G., Wang, X. F., Li, X. L., Kong, X. G., Wang, Y. J., Ning, Y. F., and Zhang, H. W.: The Asian monsoon over the past 640,000 years and ice age terminations, Nature, 534, 640–646,, 2016. 

Clark, P. U., Marshall, S. J., Clarke, G. K., Hostetler, S. W., Licciardi, J. M., and Teller, J. T.: Freshwater forcing of abrupt climate change during the last glaciation, Science, 293, 283–287,, 2001. 

Clemens, S. C., Holbourn, A., Kubota, Y., Lee, K. E., Liu, Z., Chen, G., Nelson, A., and Fox-Kemper, B.: Precession-band variance missing from East Asian monsoon runoff, Nat. Commun., 9, 1–12,, 2018. 

Cook, E. R., Anchukaitis, K. J., Buckley, B. M., D'Arrigo, R. D., Jacoby, G. C., and Wright, W. E.: Asian monsoon failure and megadrought during the last millennium, Science, 328, 486–489,, 2010. 

Dairaku, K. and Emori, S.: Dynamic and thermodynamic influences on intensified daily rainfall during the Asian summer monsoon under doubled atmospheric CO2 conditions, Geophys. Res. Lett., 33, L01704,, 2006. 

Dansgaard, W., Johnsen, S. J., Clausen, H. B., Dahl-Jensen, D., Gundestrup, N. S., Hammer, C. U., Hvidberg, C. S., Steffensen, J. P., Sveinbjörnsdottir, A. E., Jouzel, J., and Bond, G.: Evidence for general instability of past climate from a 250-kyr ice-core record, Nature, 364, 218–220,, 1993. 

Davis, S. J., Caldeira, K., and Matthews, H. D.: Future CO2 emissions and climate change from existing energy infrastructure, Science, 329, 1330–1333,, 2010. 

Ding, Z. L., Sun, J., Rutter, N. W., Rokosh, D., and Liu, T.: Changes in the sand content of loess deposits along a north to south transect of the Chinese Loess Plateau and the implications for desert variations, Quat. Res., 52, 56–62,, 1999. 

Golledge, N. R., Keller, E. D., Gomez, N., Naughten, K. A., Bernales, J., Trusel, L. D., and Edwards, T. L.: Global environmental consequences of twenty-first-century ice-sheet melt, Nature, 566, 65–72,, 2019. 

Grinsted, A., Moore, J. C., and Jevrejeva, S.: Application of the cross wavelet transform and wavelet coherence to geophysical time series, Nonlin. Processes Geophys., 11, 561–566,, 2004 (data available at:, last access: 18 July 2022). 

Guo, F., Clemens, S. C., Wang, T., Wang, Y., Liu, Y., Wu, F., Jin, Z., and Sun, Y.: Monsoon variations inferred from high-resolution geochemical records of the Linxia loess/paleosol sequence, western Chinese Loess Plateau, Catena, 198, 105019,, 2021a. 

Guo, F., Clemens, S., Liu, Y., Wang, T., Fan, H., Liu, X., and Sun, Y.: Greenhouse gases modulate the strength of millennial-scale subtropical rainfall, consistent with future predictions, ESSOAr [preprint],, 2021b. 

Guo, Z., Liu, T., Guiot, J., Wu, N., Lü, H., Han, J., and Gu, Z.: High frequency pulses of East Asian monsoon climate in the last two glaciations: link with the North Atlantic, Clim. Dynam., 12, 701–709,, 1996. 

Hao, Q., Wang, L., Oldfield, F., Peng, S., Qin, L., Song, Y., Xu, B., Qiao, Y., Bloemendal, J., and Guo, Z.: Delayed build-up of Arctic ice sheets during 400,000-year minima in insolation variability, Nature, 490, 393–396,, 2012. 

He, C., Liu, Z., Otto-Bliesner, B. L., Brady, E. C., Zhu, C., Tomas, R., Clark, P. U., Zhu, J., Jahn, A., Gu, S., and Zhang, J., Nusbaumer, J., Noone, D., Cheng, H., Wang, Y., Yan, M., and Bao, Y.: Hydroclimate footprint of pan-Asian monsoon water isotope during the last deglaciation, Sci. Adv., 7, eabe2611,, 2021. 

Heinrich, H.: Origin and consequences of cyclic ice rafting in the northeast Atlantic Ocean during the past 130,000 years, Quat. Res., 29, 142–152,, 1988. 

Hemming, S. R.: Heinrich events: Massive late Pleistocene detritus layers of the North Atlantic and their global climate imprint, Rev. Geophys., 42, RG1005,, 2004. 

Hoegh-Guldberg, O., Jacob, D., Bindi, M., Brown, S., Camilloni, I., Diedhiou, A., Djalante, R., Ebi, K., Engelbrecht, F., Guiot, J., Hijioka, Y., Mehrotra, S., Payne, A., Seneviratne, S. I., Thomas, A., Warren, R., Zhou, G., Halim, S. A., Achlatis, M., Alexander, L. V., Allen, M., Berry, P., Boyer, C., Byers, E., Brilli, L., Buckeridge, M., Cheung, W., Craig, M., Ellis, N., Evans, J., Fischer, H., Fraedrich, K., Fuss, S., Ganase, A., Gattuso, J. P., Greve, P., Bolaños, T. G., Hanasak, N, Hasegawa, T., Hayes, K., Hirsch, A., Jones, C., Jung, T., Kanninen, M., Krinner, G., Lawrence , D., Lenton, T., Ley, D., Liverman, D., Mahowald, N., McInnes , K., Meissner, K. J., Millar, R., Mintenbeck, K., Mitchell, D., Mix, A. C., Notz, D., Nurse, L., Okem, A., Olsson, L., Oppenheimer, M., Paz, S., Petersen, J., Petzold, J., Preuschmann, S., Rahman, M. F., Rogelj, J., Scheuffele, H., Schleussner , C-F. , Scott, D., Séférian, R., Sillmann, J., Singh, C., Slade, R., Stephenson, K., Stephenson, T., Sylla, M. B., Tebboth, M., Tschakert, P., Vautard, R., Wartenburger, R. , Wehner, M, Weyer, N. M., Whyte, F., Yohe, G., Zhang, X., and Zougmoré, R. B.: Impacts of 1.5 C global warming on natural and human systems, Global warming of 1.5 C, IPCC Special Report, IPCC Secretariat, 175–311, (last access: 12 July 2022), 2018. 

Huang, R., Chen, J., and Huang, G.: Characteristics and variations of the East Asian monsoon system and its impacts on climate disasters in China, Adv. Atmos. Sci., 24, 993–1023,, 2007. 

Jiang, Z., Jiang, S., Shi, Y., Liu, Z., Li, W., and Li, L.: Impact of moisture source variation on decadal-scale changes of precipitation in North China from 1951 to 2010, J. Geophys. Res.-Atmos., 122, 600–613,, 2017. 

Kitoh, A., Yukimoto, S., Noda, A., and Motoi, T.: Simulated changes in the Asian summer monsoon at times of increased atmospheric CO2, J. Meteorol. Soc. Jpn. Ser. II, 75, 1019–1031,, 1997. 

Lewkowicz, A. G. and Way, R. G.: Extremes of summer climate trigger thousands of thermokarst landslides in a High Arctic environment, Nat. Commun., 10, 1–11,, 2019. 

Li, F. X., Zhang, S. Y., Chen, D., He, L., and Gu, L. L.: Inter-decadal variability of the east Asian summer monsoon and its impact on hydrologic variables in the Haihe River Basin, China, J. Resour. Ecol. 8, 174–375,, 2017. 

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records, Paleoceanography, 20, PA1003,, 2005. 

Liu, J., Wang, B., Cane, M. A., Yim, S. Y., and Lee, J. Y.: Divergent global precipitation changes induced by natural versus anthropogenic forcing, Nature, 493, 656–659,, 2013. 

Lo, L., Chang, S. P., Wei, K. Y., Lee, S. Y., Ou, T. H., Chen, Y. C., Chuang, C. K., Mii, H. S., Burr, G. S., Chen, M. T., Tung, Y. H., Tsai, M. C., Hodell, D. A., and Shen, C. C.: Nonlinear climatic sensitivity to greenhouse gases over past 4 glacial/interglacial cycles, Sci. Rep., 7, 1–7,, 2017. 

McManus, J. F., Oppo, D. W., and Cullen, J. L.: A 0.5-million-year record of millennial-scale climate variability in the North Atlantic, Science, 283, 971–975,, 1999. 

Naafs, B. D. A., Hefter, J., and Stein, R.: Millennial-scale ice rafting events and Hudson Strait Heinrich(-like) Events during the late Pliocene and Pleistocene: a review, Quat. Sci. Rev., 80, 1–28,, 2013. 

Porter, S. C. and An, Z. S.: Correlation between climate events in the North Atlantic and China during the last glaciation, Nature, 375, 305–308,, 1995. 

Ramanswamy, V., Shine, K., Leovy, C., Wang, W. C., Rodhe, H., Wuebbles, D. J., Ding, M., Lelieveld, J., Edmonds, J. A., and Mccormick, M. P. : Radiative forcing of climate, Scientific Assessment of Ozone Depletion, NASA, Washington, 349–416, (last access: 12 July 2022), 1991. 

Shi, Y., Jiang, Z., Liu, Z., and Li, L.: A Lagrangian analysis of water vapor sources and pathways for precipitation in East China in different stages of the East Asian summer monsoon, J. Clim., 33, 977–992,, 2020. 

Sun, B., Zhu, Y., and Wang, H.: The recent interdecadal and interannual variation of water vapor transport over eastern China, Adv. Atmos. Sci., 28, 1039–1048,, 2011. 

Sun, Y., Clemens, S. C., An, Z., and Yu, Z.: Astronomical timescale and palaeoclimatic implication of stacked 3.6-Myr monsoon records from the Chinese Loess Plateau, Quat. Sci. Rev., 25, 33–48,, 2006. 

Sun, Y., Clemens, S. C., Morrill, C., Lin, X., Wang, X., and An, Z.: Influence of Atlantic meridional overturning circulation on the East Asian winter monsoon, Nat. Geosci., 5, 46–49,, 2012. 

Sun, Y., Clemens, S., Guo, F., Liu, X., Wang, Y., Yan, Y., and Liang, L.: High-sedimentation-rate loess records: A new window into understanding orbital-and millennial-scale monsoon variability, Earth-Sci. Rev., 220, 103731,, 2021a (data available at: 

Sun, Y., McManus, J. F., Clemens, S. C., Zhang, X., Vogel, H., Hodell, D. A., Guo, F., Wang, T., Liu, X., and An, Z.: Persistent orbital influence on millennial climate variability through the Pleistocene, Nat. Geosci., 14, 812–818,, 2021b.  

Swingedouw, D., Fichefet, T., Huybrechts, P., Goosse, H., Driesschaert, E., and Loutre, M. F.: Antarctic ice-sheet melting provides negative feedbacks on future climate warming, Geophys. Res. Lett., 35, L17705,, 2008. 

Thirumalai, K., Clemens, S. C., and Partin, J. W.: Methane, Monsoons, and Modulation of Millennial-Scale Climate, Geophys. Res. Lett., 47, e2020GL087613,, 2020. 

VanderPlas, J. T.: Understanding the lomb-scargle periodogram, Astrophys. J. Suppl. S., 236, 16,, 2018. 

Vecchi, G. A. and Soden, B. J.: Global warming and the weakening of the tropical circulation, J. Clim., 20, 4316–4340,, 2007. 

Wang, Y., Cheng, H., Edwards, R. L., Kong, X., Shao, X., Chen, S., and An, Z.: Millennial-and orbital-scale changes in the East Asian monsoon over the past 224,000 years, Nature, 451, 1090–1093,, 2008. 

Yin, Q. Z., Wu, Z. P., Berger, A., Goosse, H., and Hodell, D.: Insolation triggered abrupt weakening of Atlantic circulation at the end of interglacials, Science, 373, 1035–1040,, 2021. 

Zhang, H., Griffiths, M. L., Chiang, J.C.H., Kong, W., Wu, S., Atwood, A., Huang, J., Cheng, H., Ning, Y., and Xie, S.: East Asian hydroclimate modulated by the position of the westerlies during Termination I, Science, 362, 580–583,, 2018. 

Zhang, X., Lohmann, G., Knorr, G., and Purcell, C.: Abrupt glacial climate shifts controlled by ice sheet changes, Nature, 512, 290–294,, 2014. 

Zhang, X., Knorr, G., Lohmann, G., and Barker, S.: Abrupt North Atlantic circulation changes in response to gradual CO2 forcing in a glacial climate state, Nat. Geosci., 10, 518–523,, 2017. 

Zhang, Z., Li, G., Cai, Y., Cheng, X., Sun, Y., Zhao, J., Shu, P., Ma, L., and An, Z.: Millennial-Scale Monsoon Variability Modulated by Low-Latitude Insolation During the Last Glaciation. Geophys. Res. Lett., 49, e2021GL096773,, 2022. 

Short summary
Our high-resolution loess Ca/Ti record displays millennial monsoon oscillations that persist over the last 650 kyr. Wavelet results indicate the ice volume and GHG co-modulation at the 100 kyr band and GHG and local insolation forcing at the precession band for the magnitude of millennial monsoon variability of loess Ca/Ti. The inferred mechanism calls on dynamic linkages to variability in AMOC. At the precession band, combined effects of GHG and insolation lead to increased extreme rainfall.