Insolation and Greenhouse Gas Forcing of the South American Monsoon System Across Three Glacial‐Interglacial Cycles

Precipitation extremes with devastating socioeconomic consequences within the South American Monsoon System (SAMS) are expected to become more frequent in the near future. The complexity in SAMS behavior, however, poses severe challenges for reliable future projections. Thus, robust paleomonsoon records are needed to constrain the high spatiotemporal variability in the response of SAMS rainfall to different climatic drivers. This study uses Ti/Ca ratios from X‐ray fluorescence scanning of a sediment core retrieved off eastern Brazilian to trace precipitation changes over the past 322 Kyr. The results indicate that despite the spatiotemporal complexity of the SAMS, insolation forcing is the primary pacemaker of variations in the monsoonal system. Additional modulation by atmospheric pCO2 suggests that SAMS intensity over eastern Brazil will be suppressed by rising CO2 emissions in the future. Lastly, our record reveals an unprecedented strong and persistent wet period during Marine Isotope Stage 6 driven by anomalously strong trade winds.


Introduction
Over a broad region of South America, low-level atmospheric circulation and convective precipitation during austral summer are governed by the South American Monsoon System (SAMS) (Liebmann & Mechoso, 2011;Zhou & Lau, 1998). SAMS precipitation affects the basins of several major rivers in Brazil (e.g., Amazon, São Francisco, and Paraná), some of which are densely populated and intensively exploited. As such, monsoon-related hydrological extremes may have disastrous consequences for human livelihood, agriculture, energy production, and infrastructure (Marengo et al., 2012). A recent review by Pascale et al. (2019) suggests that while the core of the SAMS will likely migrate southward under continuing greenhouse gas (GHG) emissions, there is still uncertainty regarding projected precipitation changes for much of the SAMS domain. Therefore, it is an urgent task to better understand how SAMS rainfall responds to different climatic drivers, including GHG forcing.
The onset of peak SAMS activity occurs during austral summer when strong heating over central South America enhances the prevailing northeast (NE) trade winds which transport moist air from the tropical North Atlantic Ocean toward the continent (Silva & Kousky, 2012;Zhou & Lau, 1998). This moist air mass is deflected by the Andes and then heads southeastward toward the South Atlantic Convergence Zone (SACZ) (Gan et al., 2004;Zhou & Lau, 1998), a convective belt that extends from the western Amazon basin toward southeastern (SE) South America and into the western South Atlantic Ocean (Carvalho et al., 2002(Carvalho et al., , 2004 (Figure 1). Monsoon development begins in austral spring, reaches its mature phase during austral summer, and recedes in austral fall (Marengo et al., 2001;Vera et al., 2006). Across much of tropical South America, over 40% of the total annual precipitation occurs during the mature SAMS phase (Figure 1).
Available long-term records of SAMS precipitation reveal the dominant influence of precession-paced austral summer insolation changes on monsoonal rainfall intensity (Cheng et al., 2013;Cruz et al., 2005;Govin et al., 2014). However, distinct out-of-phase precipitation changes between western Amazonia and SE Brazil versus eastern Amazonia and NE Brazil over the past 250 Kyr indicate high spatial variability in rainfall intensity within the core SAMS domain (Cheng et al., 2013;Govin et al., 2014).
This study aims to disentangle the complexity in the SAMS by investigating the climatic drivers of long-term monsoonal rainfall variability at the NE extent of the SACZ, where rainfall intensity is likely sensitive to spatial shifts in the monsoon system. We present the longest record of the SAMS yet published by using Ti/Ca ratios from a marine sediment core retrieved off the coast of tropical eastern Brazil, coupled with planktic foraminiferal stable isotopes to reconstruct changes in monsoonal precipitation intensity across three glacial-interglacial cycles. The new record allows us to assess the potential roles of austral summer insolation, atmospheric pCO 2 , and latitudinal SST gradients in modulating SAMS rainfall intensity.

Regional Setting
Site M125-55-7 lies within the southward flowing Brazil Current, near the mouth of the Doce River, which is the dominant source of terrigenous sediments to the core location ( Figure 1). Climate in the Doce basin is humid (sub)tropical with temperate or hot summers and dry winters, except in the coastal lowlands where a tropical monsoon climate prevails (Alvares et al., 2013). Sediment outflow from the Doce River is dependent on rainfall amount and is highest during the austral summer rainy season (Oliveira & da Silva Quaresma, 2017). Along the coast, precipitation proceeds throughout austral winter due to the landward advection of moist air by the prevailing SE trade winds (Vera et al., 2002).

Core M125-55-7
Piston Core M125-55-7 was raised from 1,960.8 m depth off the coast of eastern Brazil (20°21.807′S, 38°37.387′W) in March/April 2016 during RV METEOR cruise M125 (Bahr et al., 2016) (Figure 1). The 1,175 cm long piston core contains bioclast-bearing clay to silty clay, alternates between dark and light gray in color, and lacks visible hiatuses (Bahr et al., 2016). We use the age model from Hou et al. (2020), which is based on 14 C dating and tuning of the benthic δ 18 O record to the LR04 benthic δ 18 O stack (associated with a 4 Kyr age uncertainty) (Lisiecki & Raymo, 2005). Core M125-55-7 spans the past 322 Kyr and has a mean sedimentation rate of 4 cm/Kyr.

XRF Scanning
The archive halves of the split cores were scanned using the fourth generation AVAATECH X-ray fluorescence (XRF) core scanner at the Institute of Earth Sciences, Heidelberg University. Elemental intensities were measured every 5 mm for 20 s at 10 kV using a slit size of 5 mm downcore and 12 mm crosscore. The surface of the sediment was smoothed and covered with a 4 μm thick Ultralene® foil to avoid contamination of the detector window and desiccation of the sediment. Processing of the raw XRF data was performed using the bAxilBatch software (Version 1.4, July 2016; www.brightspec.be). The mean temporal resolution of the XRF data is 125 years.

Stable Isotope Analysis
Ten cubic centimeters of sediment was sampled at~2.5 cm resolution. The samples were freeze dried, wet sieved over a 63 μm mesh, and oven dried overnight. Thirty to 40 Globigerinoides ruber (pink) specimens were selected from each sample and gently crushed to open shell chambers. The aliquot of cracked G. ruber (pink) tests designated for stable isotope analysis was rinsed three times with ultrapure methanol and briefly ultrasonicated between rinses. Stable isotopes were measured with a Thermo Fisher Scientific MAT 253plus mass spectrometer coupled to a Kiel IV carbonate preparation line at the Institute of Earth Sciences, Heidelberg University (Germany). An in-house carbonate standard (Solnhofen limestone) calibrated against the reference material IAEA-603 (calcite; δ 18 O VPDB = −2.37 ± 0.04‰) was used to normalize measured raw δ 18 O values to the VPDB isotope reference scale. The analytical precision for δ 18 O is better than 0.054‰ (1σ) based on replicate measurements (n = 23) of the in-house standard.

River and Seawater δ 18 O Analyses
δ 18 O analysis of seawater and river samples was performed at the GeoZentrum, Nordbayern, University of Erlangen-Nürnberg, Germany. Samples were analyzed by an automated equilibration unit (Thermo Fisher Scientific Gasbench 2) in continuous flow mode coupled to a Thermo Fisher Scientific Delta plus XP isotope ratio mass. All samples were measured in triplicates, and the reported value is the mean. Values are reported in the standard delta notation in per mil (‰) versus Vienna Standard Mean Ocean (VSMOW) following Coplen (2011). The data were corrected for instrumental drift and normalized to the VSMOW/Standard Light Antarctic Precipitation scale by assigning a value of 0‰ and −55.5‰ (δ 18 O) to VSMOW2 and Standard Light Antarctic Precipitation 2, respectively (Brand et al., 2014). External reproducibility based on repeated analyses of a control sample was better than 0.1‰ (±1σ).

Spectral Analyses
Multitaper method spectral analysis (Thomson, 1982) and Evolutive Harmonic Analysis of the evenly interpolated (130 year resolution) and linearly detrended M125-55-7 ln (Ti/Ca) record were performed using the R package "astrochron" (Meyers, 2014). Using a taner window (minimum frequency = 0.056, maximum frequency = 0.04, and roll-off rate = 10 5 ), we performed band-pass filtering of precession frequencies from the ln (Ti/Ca) record and notch filtering of the precession frequencies from the ln (Ti/Ca) record to produce a precession-free ln (Ti/Ca) record. Using the R package "biwavelet" (Gouhier et al., 2019), we analyzed the wavelet transform coherence between the evenly interpolated (600 year resolution) M125-55-7 ln (Ti/Ca) and δ 18 O IVF-SW time series over 1,000 Monte Carlo randomizations.

Interpretation of ln (Ti/Ca) Ratios
Marine sediment cores from near-shore settings are ideally suited for reconstructing past variations in terrigenous sediment input related to precipitation changes in the hinterland (Govin et al., 2012(Govin et al., , 2014. The Ti/Ca ratio of marine sediments can be used to trace the relative proportion of terrigenous siliciclastic components versus marine biogenic carbonate (Arz et al., 1998;Peterson et al., 2000;Rothwell & Croudace, 2015). Previous studies have successfully employed the Ti/Ca ratio as a proxy for fluvial terrigenous runoff, for example, to reconstruct continental precipitation changes offshore Northeastern Brazil (Arz et al., 1998;Jaeschke et al., 2007;Mulitza et al., 2017). Using Ti to represent terrigenous input by the Doce River is further corroborated by a principal component analysis of the elemental composition of the sediment performed using the program PAST 3.20 (Hammer et al., 2001; see supporting information Text S1 and Figures S1 and S2). Generally, low ln (Ti/Ca) values co-occur with warm Marine Isotope Stage (MIS) substages, while high ln (Ti/Ca) values coincide with cold MIS substages (Figure 2). The variability in ln (Ti/Ca) ratios prominently displays pacing that is in phase with orbital precession, albeit with notably damped amplitude during MIS 6 (see section 4.4).

Geophysical Research Letters
Given the near-coast locality of Site M125-55-7, terrestrial input variability may be sensitive to changes in sea level (Govin et al., 2014;Rühlemann et al., 1996). We expect ln (Ti/Ca) values to be lower during sea level high stands, when there is less erodible shelf area and the sediment source is distant from the continental slope due to retrograding sedimentary systems (and vice versa for sea level lowstands). However, the seismic stacking pattern at Site M125-55-7 and the comparison of our ln (Ti/Ca) record with benthic δ 18 O values from the same core (Hou et al., 2020) do not indicate a strong influence of sea level fluctuations on ln (Ti/ Ca) values (see Text S2 and Figures S3 and S4). ln (Ti/Ca) ratios may also be affected by changes in marine productivity, which alters the relative proportion of biogenic carbonates in marine sediments (Govin et al., 2012(Govin et al., , 2014. However, given the dominance of oligotrophic Tropical Water in the photic zone near Site M125-55-7 (de Castro et al., 2006) and the lack of correlation between ln (Ti/Ca) values and the productivity-dependent G. ruber (pink) δ 13 C record from the same core ( Figure S3), we propose that the variability in ln (Ti/Ca) ratios is not strongly affected by changes in the biogenic carbonate component. Carbonate dissolution is also unlikely to influence our ln (Ti/Ca) ratios since Site M125-55-7 (1,960.8 m) is located well above the modern (~4,000 m) and glacial (~3,500 m) lysocline in the western South Atlantic (Volbers & Henrich, 2004).
To further evaluate the robustness of our ln (Ti/Ca) record as a terrigenous runoff proxy, we compared it with the δ 18 O IVF-SW record based on G. ruber (pink) specimens from the same core ( Figure 2). The records

10.1029/2020GL087948
Geophysical Research Letters reveal a statistically robust correlation ( Figure S5a) indicating that when sediment input from the Doce River is high, offshore sea surface δ 18 O IVF-SW (a proxy for surface salinity [SSS]) is low. Today, the influence of Doce River fluvial discharge on the δ 18 O water composition of oceanic surface water is confined to the near-shore region (see Text S3, Figure S6, and Table S1) and therefore plays a minimal role in SSS changes near Site M125-55-7. Instead, we infer that the reconstructed variability in SSS at our site is more dependent on changes in the rate of surface evaporation (Bahr et al., 2013), which, among other factors, is strongly influenced by cloudiness. Therefore, an increase in convective cloudiness over eastern Brazil and the nearby oceanic region may reduce offshore evaporation rates and SSS while simultaneously intensifying continental precipitation and runoff. Oceanic evaporation and SSS are also influenced by changes in surface wind stress (Yu, 2007), where an increase in wind stress leads to an increase in the rate of evaporation and higher SSS. Notably, we observe a diversion from the otherwise out-of-phase relationship in the precession band (~21 Kyr) between ln (Ti/Ca) and δ 18 O ivf-sw during~125-175 ka (based on the wavelet transform coherence; Figure 2b), when the two records are almost perfectly in-phase. This anomalous relationship during MIS 6 may be a consequence of exceptionally strong SE trade winds during MIS 6 (Hou et al., 2020) (see ; section 4.4).

Insolation Forcing
In agreement with previous proxy-and modeling-based monsoonal studies (e.g., Cheng et al., 2012;Kutzbach, 1981), we observe a tight coupling between ln (Ti/Ca) values and precession-paced austral summer insolation cycles (Laskar et al., 2004) and find dominant peaks in power at 23 Kyr based on spectral analysis of the ln (Ti/Ca) record (Figures 3a, 3b, and S7). Our results also reveal that the sensitivity of precipitation over the Doce Basin to insolation forcing is higher when the amplitude of precession is large (e.g., MIS 5) and lower when precession is weak (e.g., MIS 2-4). Thus, we infer that both the amplitude and frequency of fluctuations in the ln (Ti/Ca) record are modulated by precession-paced insolation changes over the past 322 Kyr. During maxima in austral summer insolation, differences in land-sea heat capacity produce stronger thermal gradients, which lead to intense convection and precipitation over the continent (Silva & Kousky, 2012). During MIS 6e-c, however, the amplitude of the ln (Ti/Ca) fluctuations are muted compared to those of the insolation cycles, indicating that continental precipitation and runoff changes occurring in eastern Brazil during this period are likely modulated by other drivers (see section 4.4).
Comparison of our Southern Hemisphere ln (Ti/Ca) record to the Northern Hemisphere East Asian Summer Monsoon (Cheng et al., 2009;Dykoski et al., 2005;Kelly et al., 2006;Wang et al., 2001Wang et al., , 2008Yuan et al., 2004) reveals a long-term precession-paced antiphasing, consistent with insolation as the primary modulator of the intensity of (sub)tropical monsoons (Cheng et al., 2012;Kutzbach et al., 2008) (Figure 3). However, the periodic mismatch (e.g., during MIS 3) between our ln (Ti/Ca) record and the speleothem δ 18 O record from Botuverá Cave (subtropical SE Brazil) (Cruz et al., 2005) reveals precipitation differences within the SAMS domain, despite the shared insolation influence (Figure 3b). Thus, while our data suggest that insolation has been the primary driver of SAMS variability on orbital time scales for the past 322 Kyr, deviations in the shape and amplitude of our ln (Ti/Ca) record from orbital precession elude to modulation by other, superimposed drivers of precipitation variability in eastern Brazil. As discussed in the following, we infer that these overprinting factors might be directly and indirectly coupled to GHG forcing.

GHG Forcing
To investigate the influence of past changes in GHG concentration on monsoonal rainfall over the Doce Basin, we compared our precession-free ln (Ti/Ca) record (Figure 4d) to atmospheric pCO 2 changes (EPICA Dome C: Monnin et al., 2001;Voskok: Petit et al., 1999) (Figure 4a). Although we find that an increase in GHG concentrations is associated with a statistically significant decrease in precipitation, the large differences in amplitude between the ln (Ti/Ca) and atmospheric pCO 2 records point to nonlinear feedbacks ( Figure S5b). Since modeling studies stipulate that a doubling of CO 2 concentration lowers latitudinal temperature gradients (Rind, 1998(Rind, , 2000, we assess the ability of GHG forcing to indirectly modulate SAMS intensity through changes to the interhemispheric and intrahemispheric latitudinal temperature gradients. To test this hypothesis, we computed the correlations between the precession-free ln (Ti/Ca) record from M125-55-7, GHG concentrations, and the interhemispheric and intrahemispheric SST gradients . The interhemispheric temperature gradient is based on the difference between the IODP Site U1313 (midlatitude North Atlantic) alkenone-based SST record (Naafs et al., 2012) and the Site M125-55-7 (western tropical South Atlantic) G. ruber (pink) Mg/Ca-based SST record (Hou et al., 2020) ( Figures S8a and S8b). For computing the intrahemispheric temperature gradient, we utilized Site M125-55-7 SSTs and Site PS2489 (subantarctic zone) faunal SSTs (Becquey & Gersonde, 2003) ( Figures S8b and S8c).
In line with modeling results (Rind, 1998(Rind, , 2000, we find that the interhemispheric and intrahemispheric SST gradients increase (decrease) when atmospheric pCO 2 is lower (higher) (Figures S5c and S5d). Moreover, the increase in these latitudinal temperature gradients is associated with enhanced precipitation over the Doce Basin ( Figures S5e and S5f). This agrees with model simulations which indicate that stronger O records from Hulu, Dongge, and Sanbao caves in eastern, southcentral, and central China, respectively (green line) (Cheng et al., 2009;Dykoski et al., 2005;Kelly et al., 2006;Wang et al., 2001Wang et al., , 2008Yuan et al., 2004). Gray bars denote peaks in austral summer insolation. The orange bar denotes an interval of Marine Isotope Stage 6 during which variability in M125-55-7 ln (Ti/Ca) ratios deviate from insolation changes. Numbers at the top indicate Marine Isotope Stages 1-9.
latitudinal temperature gradients intensify the Hadley circulation (Held & Hou, 1980;Seo et al., 2014) and increase moisture availability in the (sub)tropics (Rind, 1998). A steeper intrahemispheric temperature gradient leads to intensified wintertime SE trade winds (Liu & Yang, 2003), thereby enhancing landward moisture advection and resulting in an increase in austral winter precipitation along of the coast of eastern Brazil (e.g., during the interval of weak SAMS response to insolation forcing; Figure 4c). The interhemispheric SST gradient has been shown to modulate the latitudinal position of the SACZ (Stríkis et al., 2015), which affects austral summer precipitation in the Doce Basin, since it is located at the NE extent of the SACZ and therefore sensitive to displacements in its location. This may explain the periodic differences in SAMS rainfall between the Doce Basin and the subtropical Botuverá Cave (e.g., over the last 50 Kyr), where a northward displacement of the SACZ promotes enhanced rainfall at the NE SACZ boundary and subtropical drying (and vice versa for a southward shift of the SACZ) (Cruz et al., 2005(Cruz et al., , 2006 (Figures 3b and 3c). Thus, we propose that GHG may indirectly modulate SAMS precipitation through initial modifications to both the interhemispheric and intrahemispheric temperature gradients. Our findings

10.1029/2020GL087948
Geophysical Research Letters therefore suggest that monsoonal precipitation in tropical eastern Brazil is likely to decrease in the future, owing largely to the indirect effects of rising GHG concentrations.

MIS 6
As previously discussed, during MIS 6, the long-term anticorrelation between M125-55-7 ln (Ti/Ca) ratios and δ 18 O ivf-sw values reverses, and the tight coupling between the ln (Ti/Ca) record and austral summer insolation weakens (Figures 2 and 3a). We propose that exceptionally high Southern Hemisphere intrahemispheric SST gradients during this period led to anomalously strong wintertime SE trade winds, which caused prolonged SST warming in the western tropical South Atlantic due to the accumulation of wind-advected warm-water masses (Hou et al., 2020). The intensified SE trades also increased both wintertime moisture advection from the western South Atlantic to the eastern Brazilian coast and wind-driven evaporation in the western tropical South Atlantic. This led to more positive δ 18 O IVF-SW values and higher ln (Ti/Ca) ratios due to enhanced wintertime coastal rainfall (Figure 4c). Thus, austral winter precipitation may have had a stronger influence on ln (Ti/Ca) ratios than SAMS precipitation during this period, explaining the anomalous ln (Ti/Ca) ratio response to austral summer insolation forcing during MIS 6.

Conclusions
Our 322 Kyr ln (Ti/Ca) record traces continuous changes in rainfall intensity across tropical eastern Brazil and finds insolation forcing to be the dominant pacemaker of long-term SAMS variability, despite its high spatial complexity. Additionally, we find that the intensity of SAMS rainfall responds negatively to GHG forcing, through indirect mechanisms. We propose that GHG forcing directly influences the interhemispheric and intrahemispheric latitudinal temperature gradients, which in turn modify the strength of tropical atmospheric circulation and precipitation. During MIS 6, we discern an anomalous muted SAMS response to insolation forcing and a shift from the typical out-of-phase relationship between ln (Ti/Ca) and δ 18 O IVF-SW values. We attribute this to exceptionally strong wintertime SE trade winds which caused austral winter precipitation to dominate the ln (Ti/Ca) record. Our findings suggest that SAMS rainfall intensity over tropical eastern Brazil will likely decrease due to the anthropogenic-related rise in GHG concentrations. This will likely have important implications for government decisions regarding the timely adoption of mitigation strategies in the near future.