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 ﬂood and drought cycles. Therefore, explor-ing what drives the millennial-scale variability is of signiﬁcant importance for future prediction of extreme cli-mates. 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 2 during intermediate glaciations lead to increased variability in North Atlantic stratiﬁcation and Atlantic meridional overturning circulation, propagat-ing abrupt climate changes into East Asia via the westerlies. Within the 100 kyr cycle, precession variability further en-hances the response, showing that stronger insolation and increased atmospheric GHG cause increases in the MMV of EASM rainfall. These ﬁndings indicate increased extreme precipitation events under future warming scenarios, consistent with model results.


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., 1996Sun et al., 2012Sun et al., , 2021a, 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. Un-derstanding 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 δ 18 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 δ 18 O record (δ 18 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 δ 18 O sp and postulated that it is, instead, modulated by internal processes related to the cryosphere. This work also raised the possibility that δ 18 O sp is decoupled from regional Asian monsoon rainfall over millennial timescales . The latest research suggests that the MMV through the Pleistocene is influenced by both glacial boundary conditions and orbital configurations . 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 . The Linxia profile is well-suited for reconstructing rapid monsoon changes because it is located in the monsoon frontal zone and sensitive to highand low-latitude climate variability . To address the above questions, we have generated a highresolution summer monsoon proxy (Ca/Ti) from Linxia on the western CLP (Fig. 1). The Ca/Ti ratio is a precipitationsensitive proxy linked to summer monsoon rainfall . Low values of Ca/Ti indicate stronger Ca leaching associated with intensified summer rainfall. The new precipitation proxy (Ca/Ti) and δ 18 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 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 δ 18 O sp ; δ 18 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).

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 .

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 δ 18 O stack (Porter and An, 1995;Hao et al., 2012;Sun et al., 2021a). A similar approach utilizes correlation to speleothem δ 18 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 δ 18 O records back to 650 ka . Detailed information of loess-paleosol boundaries is shown in Fig. 2. The first set of control points tie the loesspaleosol boundaries S 6 to S 0 to the timing of the glacial terminations and inceptions in speleothem δ 18 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 δ 18 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.
The composite speleothem δ 18 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 δ 18 O. The age differences of most glacial terminations are around 2-4 kyr (Fig. 2b, Sun et al., 2021a). The speleothem δ 18 O synchronized age model is compared with the benthic δ 18 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 δ 18 O and to speleothem δ 18 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 millennialscale 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).

Spectrum and wavelet coherence analyses
In order to estimate the MMV, loess Ca/Ti and δ 18 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., 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 (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 (Vander-Plas, 2018).
The normalized and combined orbital parameters eccentricity, tilt and negative precession (ETP); GHG; insolation; and benthic δ 18 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 2 ] and [CH 4 ]) and the pre-industrial greenhouse gas level ([CO 2 ] 0 = 280 ppm, [CH 4 ] 0 = 700 ppb) (Ramaswamy et al., 1991). While CH 4 contributes < 5 %, we calculated the RF GHG using both CO 2 and CH 4 . The equation used to determine RF GHG is as follows (Lo et al., 2017).

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 , with relatively low values during the interglacials. The highresolution δ 18 O 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 δ 18 O (Lisiecki and Raymo, 2005) and speleothem δ 18 O (Cheng et al., 2016), respectively (Fig. 2b).
Both Ca/Ti and δ 18 O sp records show clear millennialscale 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  (Cheng et al., 2016) and benthic δ 18 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 . Light blue bands denote the interglacial times. The short green rectangles represent the age differences between the two age models. 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.
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 δ 18 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 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). δ 18 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).

Orbital-scale modulation factors for MMV of the EASM
Geological records and modeling results indicate that highlatitude 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 δ 18 O and global benthic δ 18 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 δ 18 O stack is relatively weak and variable at the 41 kyr band (δ 18 O sp ; 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 δ 18 O 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 δ 18 O 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 δ 18 O (Fig. S3b and d). Thus, GHG RF does play a role in modulating Ca/Ti MMV but not that of δ 18 O sp at the precession band, indicating a difference in the millennial-scale response of these two proxies. We investi- gate 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 δ 18 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 δ 18 O is modulated by local win-ter insolation, opposite to the hypothesis calling on Northern Hemisphere summer insolation (Cheng et al., 2016).

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., 2012Sun et al., , 2021a. 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 2 concentrations can trigger rapid climate transitions (Zhang et al., 2014. 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 2 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 2 during Heinrich stadial cold events can trigger abrupt transitions to warm conditions. Decreasing CO 2 during warm events leads to abrupt cooling transitions . Therefore, CO 2 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 2 , 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 δ 18 O sp 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 δ 18 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 δ 18 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 δ 18 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 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.

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 δ 18 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.
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 http://grinsted.github.io/wavelet-coherence/ (Grinsted et al., 2004).
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.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.