South American Monsoon variability over the last millennium in paleoclimate records and isotope-enabled climate models

. The South American Summer Monsoon System (SAMSSASM) is the main driver of regional hydroclimate variability across tropical and subtropical South America. It is best recorded on paleoclimatic timescales by stable oxygen isotope proxies, which are more spatially representative of regional hydroclimate than proxies for local precipitation alone. Network studies of proxies that can isolate regional influences lend particular insight into various environmental characteristics that modulate hydroclimate, such as atmospheric circulation variability and changes in the regional energy budget as well as understanding the climate system sensitivity to external forcings. ToWe extract the coherent modes of variability of the SAMSSASM over the Last Millennium (LM), we use using a Monte Carlo Empirical Orthogonal Function (MCEOF) decomposition of 14 d 18 O proxy records and compare them with modes extracted from a similar decomposition using decomposed from isotope-enabled climate models data. The two leading modes reflect the isotopic expressionvariability associated with 1) thermodynamic changes driving of the upper-tropospheric

climate models, complicating the attribution of changes on these timescales to specific forcings and pointing toward areas of important model development.

Introduction
The annual hydrologic cycle in northern and central South America is dominated by distinct wet and dry seasons following the development and demise of the South American Summer Monsoon System (SAMSSASM) (Vera et al., 2006).
The mature phase of the monsoon occurs during the austral summer in December, January, and February (DJF) when the region receives the majority of its annual precipitation (Zhou and Lau, 1998;Garreaud et al., 2009).Monsoon rainfall is critical to many environmental services, the integrity of indigenous traditions and practices, drinking water, agriculture, and hydroelectric power generation, all of which are vulnerable to changes in the monsoon system.Although historical hydroclimate observations do exist from this region over the 20 th century, such records are not long enough to characterize multidecadal to centennial changes, which are critical for interpreting and constraining projections of future change.
Observational records have, however, aided in defining some of the characteristics of present-day SAMSSASM dynamics.The onset of the SAMSSASM is initiated by the increase in austral summer insolation over the southern hemisphere driving the seasonal southward progression of the Intertropical Convergence Zone (ITCZ) and strengthening the onshore flow of moisture from the tropical Atlantic toward the South American continent (Zhou and Lau, 1998).As the monsoon matures, near-surface easterly winds are channeled southeastward along the Andes mountains into the South American Low-Level Jet (SALLJ).The SALLJ transports moisture into the subtropics and feeds the South Atlantic Convergence Zone (SACZ), an axis of NW-SE convective activity in southeastern Brazil (Kodama, 1992, 1993, Carvalho et al., 2004).
Convection over the Amazon, Andes Mountains, and to a lesser extent, inover the SACZ region generates diabatic heating anomalies via the release of latent heat.This triggers an easterly Rossby wave response known as the Bolivian High -Nordeste Low (BH-NL) system (Lenters and Cook, 1997;Chen et al., 1999) and leads to subsidence to the east of the monsoon heating region over northeastern Brazil and the tropical South Atlantic, intensifying the South Atlantic Subtropical Anticyclone (SASA) (Rodwell and Hoskins, 2001).The strength of this monsoonal circulation varies on multiple timescales and is modulated interannually by the El Niño-Southern Oscillation (ENSO) phenomenon (Vuille and Werner, 2005;Garreaud et al., 2009;Sulca et al., 2018;Cai et al., 2020) and on decadal to multidecadal timescales by modes of variability in both the Atlantic and Pacific, via perturbations of the local Walker and Hadley circulation (He et al., 2021).
Because these dynamics are linked to hydroclimatic change, the recent development of paleo-hydroclimatic records from the SAMSSASM region has extended our characterization of monsoon variability in the past with increasing spatial and temporal resolution.Stable oxygen isotope records form the basis for much of our understanding about past variability of the SAMSSASM.The existing sub-decadally resolved oxygen isotope records spanning the last millennium (LM, 850-1850 CE) are known to be representative of regional variability related to the summer monsoon (Vuille et al., 2003;Vuille and Werner, 2005).High-resolution oxygen isotope paleorecords are widely available from archives in South America: ice cores in Andean glaciers (e.g.Thompson et al., 1985;Vimeux et al., 2009), lake sediment cores from closed basin lakes (e.g., Bird et al. 2011), speleothems from karstic regions (e.g., Novello et al., 2018;Campos et al. 2019), and tree rings in high-altitude regions (e.g., Baker et al., 2015;Rodriguez-Caton et al., 2022).
In tropical South America, isotopic fractionation through Rayleigh distillation occurs during rainout along moisture trajectories where it is governed by convective processes rather than temperature changes.Local anti-correlation of precipitation amount with isotopic values due to the amount effect has also been documented, especially in close proximity to the Atlantic coast (Vuille et al., 2003;Moquet et al., 2016).More recent studies document the importance of an interplay of fractionation processes, finding that oxygen isotopes integrate various controls on tropical hydroclimate, such as circulation changes, energetic processes, microphysical processes, and local environmental conditions (Risi et al., 2008;Konecky et al., 2019).
The current interpretation of isotopic signatures in South American paleorecords therefore allows insight into both local climate and regional atmospheric circulation.Along the moisture trajectory, convective rainout follows a Rayleigh Distillation model (Dansgaard, 1964;Vuille et al., 2012) with the remaining water vapor becoming isotopically depleted downstream of the source.Due to evapotranspiration within the Amazon basin, a reduced depletion gradient is seen across the Amazon (Salati et al., 1979;Ampuero et al., 2020).Secondary effects also modify isotopic variability during the monsoon season, such as mixing of air masses with different moisture sources (e.g. from the tropical North or South Atlantic) or significant rainfall contributions from other circulation systems outside of the monsoon season (Cruz et al., 2005;Apaéstegui et al., 2018).Ultimately, the degree of rainout upstream associated with the SAMS is the dominant common denominator influencing all records within the monsoon domain, explaining why records from different archives located in different environments separated by thousands of miles share a common signal (Vuille et al., 2012;Campos et al., 2019;Ampuero et al., 2020;Novello et al., 2021).
The current interpretation of isotopic signatures in tropical South American paleorecords allows insight into both local climate and regional atmospheric circulation.Isotopic fractionation occurs to a first order via Rayleigh distillation during rainout along moisture trajectories where it is governed by hydroclimate rather than temperature changes.Deep convective activity drives condensation and precipitation following a Rayleigh Distillation model (Dansgaard, 1964;Grootes et al., 1989;Hoffman et al., 2003;Vimeux et al., 2005;Vuille et al., 2012;Ampuero et al., 2020) across northern South America, the Amazon basin, and in the eastern flank of the Andes, leaving remaining water vapor isotopically depleted downstream of the tropical North Atlantic.Further downstream in southeastern Brazil, the role of Rayleigh distillation in driving isotopic variability decreases.Local anti-correlation of precipitation amount with isotopic values due to the amount effect has also been documented, especially in close proximity to the Atlantic coast (Vuille et al., 2003;Moquet et al., 2016).
A number of secondary processes which are integrated into paleorecords must be considered as well to capture the full interplay of fractionation processes, including circulation changes, energetic and microphysical processes, and local environmental conditions (Risi et al., 2008;Konecky et al., 2019).Local precipitation recycling, enhanced by evapotranspiration within the Amazon basin, is evident from the reduced depletion gradient seen across the Amazon (Salati et al., 1979;Ampuero et al., 2020).Circulation changes also modify isotopic variability during the monsoon season, such as mixing of air masses with different moisture sources (e.g. from the tropical North or South Atlantic) and significant rainfall contributions from other circulation systems outside of the monsoon season can influence paleorecord variability (Cruz et al., 2005;Apaéstegui et al., 2018).Satellite observations on intraseasonal and interannual timescales show that depletion of heavy isotopes in precipitation can sometimes exceed expected values modeled by Rayleigh distillation alone due to the impact of convective processes such as dry air entrainment and high humidity levels reducing rain re-evaporation (Samuels-Crow et al., 2014;Aron et al., 2021).However, the dominant common denominator influencing all records within the monsoon domain is the degree of rainout upstream associated with deep SASM convection, explaining why records from different archives located in different environments separated by thousands of miles share a common signal (Vuille et al., 2012;Campos et al., 2019;Ampuero et al., 2020;Novello et al., 2021).
The wealth of information archived by paleoclimate records is clear, but the balance between site-specific and regional information is unique to each record.Today, the growing catalog of existing records is ripe for a more systematic analysis that can extract coherent regional changes.To isolate the influence of the SAMSSASM requires the application of a network analysis such as Monte Carlo Empirical Orthogonal Function (MCEOF) decomposition (Anchukaitis and Tierney, 2013).A coherency analysis, demonstrated by previous network studies both in South America and in other regions, is uniquely capable of highlighting the regional system variability, while minimizing the influence of local climate at individual proxy sites (Anchukaitis and Tierney, 2013;Tyler et al., 2015;Deininger et al., 2017;Falster et al., 2018;Campos et al., 2019;Novello et al., 2021).Campos et al. (2019) describe characteristics of the monsoon system with an emphasis on the temporal variability as captured by leading modes of variability from a multi-proxy network.Their analysis characterizes the first two modes of variability as well as critical time periods detected by the principal component analysis and compares these to EOFs of present-day precipitation and other hydroclimate records from the region.Atwood et al. (2021) also perform an EOF analysis across the global tropics using hydroclimate proxies for the LM, including over tropical South America.While both these studies provide valuable new insight into South American hydroclimate during the LM, they are both based on a network of sites that lacks data in the core monsoon region of the southwestern Amazon as well as at the most distal monsoon sites in the southern tropical Andes.
Additionally, while the dynamical characteristics of the SAMSSASM climatology are well known, a mechanistic understanding of the factors underlying the regional isotopic variability has yet to be established.The emergence of new tropical isotopic records enables the reconstruction of a more regionally coherent picture of past hydroclimatic change than the relatively heterogeneous modes described by proxies recording local precipitation only.Combining this network approach with an analysis of isotope-enabled climate model simulations allows mechanisms of past change to be identified and for proxy records to be dynamically contextualized within past variability.
We present an analysis of the SAMSSASM using a proxy network derived solely from oxygen isotope records covering the LM, all located within the SAMSSASM domain.Included in this analysis are three records that have not been used in any prior proxy network analysis (Boto, PIM4, MV; see Table 1) and help improve the spatial coverage of monsoon impacts in hitherto undersampled regions.We show the first two common modes of isotopic variability and contextualize them within our current understanding of present-day monsoon dynamics.In parallel, an identical analysis was performed on simulations from two isotope-enabled climate models and compared with the results from the proxy analysis.Finally, we interpret the monsoon behavior during the Medieval Climate Anomaly (MCA, 850-1160 CE) and Little Ice Age (LIA, 1489-1850 CE) (time periods as defined in Campos et al., (2019)), which emerge as periods of multi-centennial deviations from the mean state of the LM.

MCEOF analysis approach
The variability of the SAMSSASM was investigated by using an Empirical Orthogonal Function (EOF) decomposition of ing oxygen isotope data from both a network of paleoclimate records networks and isotope-enabled climate model simulations.into the leading modes of variability using a MCEOF analysis originally developed for the analysis of West African lacustrine sediments (Anchukaitis and Tierney, 2013).Since its development, this technique has been applied to a variety of proxy networks, including in South America by Campos et al. (2019) and Novello et al. (2021).
The EOF decomposition allows for discrete characteristicsaspects of the monsoon system variability to be isolated and displayed asinto the leading modes.A Monte Carlo EOF (MCEOF) wasused forapplied to the network of paleorecords network, with the applying a Monte Carlo resampling ensuring that while the application of a Monte Carlo method conserves uncertainties associated with the age models of the individual paleoclimate records used as inputwere conserved.This method was originally developed for the analysis of West African lacustrine sediments (Anchukaitis and Tierney, 2013).
Since its development, this technique has been applied to a variety of proxy networks, including in South America by Campos et al. (2019) and Novello et al. (2021).For all EOF analyses, the first two modes were considered significant using North's rule of thumb (North et al., 1982) and the fact that they aligned with a dynamically and physically plausible interpretation.
Each age model ensemble member was linearly interpolated between depth/age points and fit by the depth/time series from each proxy sample to establish an ensemble of proxy time series (Deininger et al., 2017;Campos et al., 2019;Novello et al., 2021).In cases where individual samples did not span the entire time period of study, multiple shorter samples were splicedmerged togetherto form composite records.RecordsOxygen isotope time series were resampled to annual resolution and truncated to 850 -1850 CE, resulting in a 1,000-member ensemble of fourteen complete records spanning the LM.Each of the 1,000 ensemble sets of oxygen isotope records was subsequently decomposed into leading modes of spatial and temporal variability, and the median of the ensemble modes was analyzed in this study (Anchukaitis and Tierney, 2013;Deininger et al., 2017).

Record processing
The spatial distribution of the proxy archives is shown in Fig. 1 and more detailed information is provided in Table 1.Much of the data used here has been published, with the exception of the MV1 and MV30 samples, which wereas collected in the Mata Virgem Cave in Brazil (see Azevedo et al., 2019 for site information).Published data were accessed through the National Centers for Environmental Information (NCEI) Paleo Data server (https://www.ncdc.noaa.gov/paleosearch/)or through the University of São Paulo Geosciences speleothem database.Institute of Geosciences.Eleven of tThe records used for this analysis have sub-decadal to annual resolution, as determined by the number of isotopic samples available in between U/Th dates that established record age chronologies.using U/Th radioisotope dating and were found to be sub-decadally to annually resolved.This dating method provides an uncertainty of approximately 1% at 2 sigma confidence levels, yielding errors on the order of a few years (less than 10 years in the vast majority of samples used in this study).The remaining three records in the network were published with annual resolution (Quelccaya, Pumacocha, ALHO6), determined from visible stratigraphy and confirmed with radiometric dating where possible (Bird et al., 2011;Thompson et al., 2013;Novello et al., 2016).These records were linearly interpolated to exact calendar years where necessary, with the finalized record duplicated for each ensemble member.For those samples that do not cover the entire LM, composite records were synthesized by merging multiple samples from the same cave systemconsistent with original data in the original publications (indicated in Table 1) and samples with hiatuses were interpolated to ensure a continuous record for analysis.Samples actively growing upon collection (P00-H1, CR1, SBE3, PIM4) were assigned an age tie point without error to the year of collection at zero depth to anchor the chronology (Anchukaitis and Tierney, 2013).Some age points of the TM0 record were published with multiple samples indicated for a given depth; the employed final dates were calculated as the mean of the samples analyzed.Age models were constructed with strict superposition enforcement.Those age ties with processing errors or dates leading to age reversals were excluded from analysis to ensure superposition of the age chronologies (JAR1, JAR4, MV30).In the Monte Carlo resampling, simulated age tie values that violated superposition were rejected and recalculated.Isotopic measurements excluded from previously published analyses were likewise excluded here (SBE3, PIM4), as were individual data points exceeding 3 standard deviations (MV30).Multiple isotopic measurements from a single depth were averaged to establish one value per depth (PAL3).

Gap interpolation
Samples that contained hiatuses (SBE3, PIM4) or multiple discontinuous samples (MV1, MV30) had gaps filled after the pre-and post-gap segments underwent ensemble generation.Gap filling proceeded in two steps: (1) a linear interpolation was used to fill gaps in the time series; and (2) the isotopic value for each year of the gap was resampled from a normal distribution centered on the interpolated point from a normal distribution scaled to the average standard deviation of the thirty years before and after the gap.This procedure bridges a gap in the data in order to ensure isotopic time series have sharedcommon time steps when applying the EOF analysismeet the continuous data requirement of the MCEOF procedure while also maintaining the intrinsic characteristics of variability present in a given record.For bridging the gap between the MV1 and MV30 samples, MV30 was bias-corrected to have the same mean as MV1; the two records were standardized, and the gap was filled as described above.As a result of the Monte Carlo age model resampling of the MV1 and MV30 samples, a small number of ensemble members returned samples that overlapped or doid not have a gap in data, in which case the samples were merged or taken to be continuous, respectively.

Sample merging
Six of the isotopic records used in this analysis were composited from multiple raw samples (MV, JAR, PAR, SBE3+SMT5, PAL, Boto).With the exception of the unpublished MV record, this compositing is consistent with the original publications of these records.As in Novello et al. (2021), composite records were formed by merging individual ensemble members using three steps: (1) ensemble members in the period of overlap were fit with a cubic spline function and synchronized to annual resolution; (2) each time series was standardized by its respective mean and standard deviation during the LM; and (3) samples were averaged together to create one time series that was reconstructed by inverting the normalization of the record with the longest coverage during the LM.Finally, records were truncated to 850 -1850 CEestrictly to cover the LM and linearly interpolatedresampled to haveachieve annual data resolution.

Isotope-enabled model analysis and pseudoproxy experiments
Output from tTwo previously published isotope-enabled coupled climate models experiments were used to investigate mechanisms associated with climate variability over the LM: (1) the water isotope-enabled NCAR Community Earth System Model (iCESM) (Brady et al., 2019); and (2) the Goddard Institute for Space Studies coupled oceanatmosphere model (GISS-E2-R) (Schmidt et al, 2007;Schmidt et al., 2014;Colose et al., 2016a,b).In these isotope-enabled simulations, stable oxygen and hydrogen isotopes have been added as tracers of the hydrological cycle, enabling a more direct comparison of model and proxy data and a better understanding of water cycle dynamics.The iCESM was run with 1.9° x 2° horizontal resolution and 30 vertical levels.Analyzed simulations were from the isotope-enabled Last Millennium Ensemble project, including full and single forcing (volcanic, solar, and orbital) experiments following the CESM Last Millennium Ensemble (iLME, Otto-Bliesner et al., 2016;Brady et al., 2019).The GISS-E2-R model was run at 2° x 2°h orizontal resolution with 40 vertical levels.We analyzed the simulation from the E4rhLMgTKck experiment (Blake et al., 2018), which contains a full suite of transient external forcings.Both models include forcings selected from the PMIP3 Last Millennium forcing set v1.0 (Schmidt et al., 2011), although they rely on somewhat different forcing reconstructions.
The iCESM is known to have a negative isotope bias (Brady et al., 2019), overestimating the strength of the relationship between precipitation and d 18 O.This is particularly true at higher rainfall rates, though spatial patterns are consistent with observations and values are similar to the proxy data (Brady et al., 2019;Stevenson et al., 2019).Isotopes in the GISS-ER2-R are highly correlated with observations, but show a weak dependence to changes in temperature and precipitation (Schmidt et al., 2007); the GISS-E2-R reproduces key isotopic features of the South American monsoon (Colose et al., 2016a).
We extracted the leading EOFs from two model-derived datasets that are complimentary in their comparison with proxy records of the monsoon domain: one based on spatially continuous data across the full monsoon domain and one using a pseudoproxy experiment (PPE) approachframework (Smerdon, 2012(Smerdon, , 2016)).In both cases, the EOFs wereas based on precipitation-weighted d 18 O during DJF, which corresponds to the mature monsoon season when isotopic proxies receive the largest fraction of their total annual precipitation (Fig. 1).
For the spatially continuous EOF approach, data from both the full-forcing ensemble members of the iCESM and GISS-E2-R model as well as the single forcing simulations of the iLME using the volcanic, solar, and orbital forcings were used.The land use/land cover and greenhouse gas forced ensemble members do not exhibit significant change in the forcing data sets during the LM and were thus excluded from the analysis.
For the pseudoproxy approach, the iLME full-forcing experiment was used and a spatial sub-sampling was applied by extracting an ensemble of pseudoproxies from the latitude-longitude grid cells of the climate model simulations corresponding to real-world proxy locations.The use of PPEs as an independent way to evaluate the signal-to-noise statistical profiles of proxy data was based on examples used in past literature (Smerdon et al., 2016;Yun et al., 2021).
Randomly generated Gaussian white noise was then added to tThese subsampledsub-sampled time series were then added to randomly generated Gaussian white noise, effectively scaling the model-simulated isotopic time series to more realistically replicatereflect the internal variability of the proxy systems.In our specific experiment, the white noise was used to represent unresolved sub-grid scale processes (Berner et al., 2017) which constitute the noise inherent to the climate system in both the real world and model simulations.While this representation cannot fully account for the cascading influence of sub-grid scale processes resulting from external forcings, it does allow the EOF analysis to decompose spatial modes from model data based on a broader representation of multi-frequency (i.e.multi-scale) information as captured by the proxy record.We interpret the PPE as a filter for analyzing the variability of the climate signal in both the models and the network of paleorecords.Following Smerdon (2012), this white noise perturbation was drawn from a Gaussian distribution and added to model-derived pseudoproxy isotopic time series using signal-to-noise (SNR) ratios of 0.5 and 0.25.Previous PPE studies focused on temperature have found SNR = 0.5 to be a realistic estimate (Mann et al., 2007;Smerdon, 2012;Wang et al., 2014).Given the larger variability of precipitation compared to temperature, SNR = 0.5 and SNR = 0.25 were considered here to be bounding estimates of SNR for hydrological proxies, though further work to understand noise in hydroclimate proxies is important for a more sophisticated hydrological PPE.

Present-day composite analysis
Our physical interpretations of the leading MCEOFs are based on modern observations for the DJF season.The upper tropospheric (200 hPa) wind field data are based on ERA5 Reanalysis for present-day (1979-2019) (Hersbach et al., 2020).Precipitation station data were extracted from GHCN (Menne et al., 2012) and supplemented with data obtained from the national weather services of Bolivia, Peru, and Brazil.A time series of seasonal (DJF) precipitation at the site of each proxy record was created based on the median of all stations available within a 2° radius.The data period was 1975-2014 for all locations, except for those near the PIM4 site, where four years of data are missing (1998,2002,2003,2004).Outgoing longwave radiation (OLR), a commonly employed proxy for tropical convection, was analyzed using the NOAA interpolated OLR data set (Liebmann and Smith, 1996) for the same period as the precipitation data.OLR data for the SACZ region was decomposed using an EOF decomposition and the resulting PC1 correlated in space with the gridded OLR product and precipitation data.
A composite analysis was performed for the MCA (850-1160 CE) and LIA (1489-1850 CE) periods based on both paleoclimatic records and climate model data.The MCA and LIA periods were defined using the Campos et al. (2019) change point analysis, which highlighted major departures in South American hydroclimate.Anomalies were calculated for these periods relative to the LM.Atmospheric circulation and hydroclimate diagnostics of mid-tropospheric vertical velocity (ω500, 500 hPa), precipitation, and precipitation-weightedd 18 O of precipitation were analyzed by averaging these variables for the mature monsoon season (DJF) for the three full-forcing ensemble members from the iLME experiment.Proxy composites were calculated during the above time periods from the average of annually interpolated composite records generated for the MCEOF analysis.Statistical significance was calculated using a two-tailed students t-test and evaluated at p ≤ 0.1.

MCEOF results from proxy network
The results of the proxy-based MCEOF analysis are shown in  The leading mode of isotopic variability explains approximately 24% of the total variance.This mode is characterized byis interpreted as representative of the in-phase isotopic variability in across the core monsoon region and is shown in PC1 to vary on centennial timescales.The majority of the proxy records in the region of strong monsoonal convection (the center and western region of the monsoon domain) vary in phase, while the northeastern/eastern records are negatively correlated with PC1, implying that they vary out of phase (Fig. 2a).The PC time series indicates that records over the southwestern SAMSSASM domain exhibited anomalously positive isotopic values during the MCA and anomalously negative values during the LIA, while records to the northeast show opposite isotopic anomalies (Figs.2a-b).
The second mode explains approximately 12% of the total variance and exhibits the strongest cluster of loadings in central-eastern Brazil, a region influenced by the SACZ (Fig. 2c).Proxy records that project strongly onto this mode influenced by SACZ activity include ALHO6, SBE3+SMT5, TM0 and DV2, while records over the southern central Andes (Quelccaya and especially Boto) vary out of phase.There is a prolonged negative excursion of PC2 during the core of the LIA (approximately 1550 -1750 CE) (Fig. 2d).The Boto record, a regional anomaly, loads strongly onto this mode, indicating a strong local response inversely correlated with spatiotemporal dynamics that are primarily dominant in southeastern Brazil.isotope dipole (Cruz et al., 2009;Cheng et al., 2013).This characteristic aspect of the SAMSSASM is notably absent in the GISS EOF1, where positive loadings extend across the entire domain analyzed.The largest loadings present in EOF2 are located in central-eastern Brazil, the region corresponding to SACZ activity, again consistent with the proxy-based EOF2 (Fig. 2c),.Modeled although the loadings over the northern tropical Andes differ slightly from the proxy-based EOF2.The second model-derived mode explains 14% and 17% of the data in GISS and iCESM models, respectively (Figs. 3h-i), a somewhat larger explained variance than in the proxy-based MCEOF2 (~12%).Nevertheless, the representation of SACZ variability present in all models indicates models represent the dynamics of this region well.The emergence of isotopic variability over the SACZ-region in this second mode is more apparent in the iCESM results than the GISS model-derived EOF2.The GISS model shows a decreased explained variance of this second mode relative to the iCESM that is more on par with the partitioning of the proxy-based MCEOF2.

EOF1
The EOF decomposition of the iLME single-forcing experiments (Figs.3c-e, j-l) yields only minor departures from patterns of the full-forcing experiment and replicates the proxy-based MCEOF patterns equally well.This underscores the inherent stability of the spatial characteristics of these modes of the SAMSSASM, regardless of the external forcing applied to the experiment.The presence of almost identical leading modes across all single-forcing experiments indicates that these modes are intrinsic to the system itself, rather than a consequence of external forcings and likely result from time-invariant boundary conditions such as the shape and latitudinal extent of the continent, the Andean orography, and the intersection of the landmass with the equatorial trade winds.

Pseudoproxy experiment analysis
The EOF decomposition of the pseudoproxy network derived from the addition of white noise to the modeled isotopic time series results in a reduction of explained variance relative to the proxy-based modes, likely indicating the experimental SNRs applied to this region were too small relative to the real climate system .Despite the decrease in explained variance, the spatial patterns of EOF1 and 2 from the pseudoproxy network replicate those described above from both proxies and models.The main monsoon and SACZ footprints in EOF1 and EOF2 are readily apparent using a SNR = 0.5.Increasing the amount of white noise added to the pseudoproxy records to SNR = 0.25 results in a decrease of the total variance explained (EOF1: ~16 to ~10%; EOF2: ~12% to ~9%) and the overall spatial patterns of the modes becoming more muted.Despite this, the first mode remains representative of SASM activity (Fig. 3g) and the second mode retains the regional loadings in the region of central-eastern Brazil associated with the SACZ (Fig. 3n).This SNR sensitivity test confirms the robustness of the modeled results and the model's ability to reproduce the main modes of South American monsoon variability in a pseudoproxy framework.The larger explained variance of the proxy-based modes relative to those from the pseudoproxy network of SNR = 0.5 suggests that the noise characteristics of the SAMS proxy network may be very region-and proxy-specific, and that isotopic proxy networks in some cases may exhibit less noise than currently established within the literature for other multiproxy networks (Mann et al., 2007;Smerdon, 2012).More detailed estimates of environmental SNRs would open the possibility for closer calibration between proxy network analysis and estimation of regional climate variability via pseudoproxy experiments.

MCEOF1
The twofirst modes extracted from the proxy network via the MCEOF decomposition corresponds to distinct patterns that areis inherent to the SAMSSASM.MCEOF1 represents the isotopic responsevariability corresponding to changes in the Bolivian High-Nordeste Low (BH-NL) system: an upper tropospheric Rossby wave induced by diabatic heating from convection over the core monsoon region in the western Amazon (Lenters and Cook, 1997;Chen et al., 1999) (Fig. 4).This stationary wave develops downstream of deep tropical convection, establishing an upper-level anticyclone over Bolivia, known as the 'Bolivian High' and is balanced by enhanced subsidence over northeastern Brazil where a deep trough develops, known as the 'Nordeste Low.'This mechanism results in a precipitation dipole between the tropical Andes/western Amazon basin and northeastern Brazil, which has been documented on multiple timescales, including orbital (Cruz et al., 2009;Cheng et al., 2013), centennial (Apaéstegui et al., 2018), multidecadal (Campos et al., 2019), interannual (Vuille and Werner, 2005) and intraseasonal (Sulca et al., 2016).Variability of the BH-NL system, excited by the diabatic forcing from latent heat release during convection, is therepresents a dynamic response to changes in the waterhydrologic cycle, which also imprints directly on the variability of the oxygen isotopice composition recorded in paleorecords.Thus, variability of oxygen isotopes as captured by the first mode represents is essentially a proxy for the monsoondynamic response to changes in internal or external forcingthe water cycle.Long-termMulti-centennial departures exhibited in PC1 support the current understanding of the MCA as a period characterized by relatively weak large-scale monsoon circulation, while the SASM was intensified during the LIA (Novello et al., 2012;Rojas et al., 2016).The reduced departures from the mean of PC1 during the MCA when compared to the LIA support evidence from individual proxy records that the monsoon intensification during the LIA represented a stronger regional departure from the monsoon mean state than the relatively subdued weakening during the MCA (Vuille et al., 2012).Here, monsoon intensification (weakening) refers to the intensification (weakening) of the upper-tropospheric BH-NL system (Chen et al., 1999;Sulca et al. 2016;Campos et al., 2019) and its surface expression of enhanced monsoon precipitation over the western Amazon/tropical Andes and enhanced aridity over northeastern Brazil, while a reduction in intensity implies the inverse.Changes in Rayleighisotopic fractionation of water vapor from condensation and rainout during strong/weak convectionve periods resultsleads to in this precipitation dipole being prominently displayed in stable oxygen isotopes of precipitation (Vuille and Werner, 2005;Sturm et al., 2007).

MCEOF2
The second mode captures the variability of the SACZ during the LM as seen through oxygen isotopes in our proxy network, consistent with results presented by Campos et al. (2019) (Figs. 4c-d).The present-day footprint of SACZ variability is represented by correlations of both the OLR field and precipitation derived from station data with PC1 of the OLR EOF decomposition over the SACZ region.Loadings of opposite sign are seen in the core SACZ region and the subsidence region to the southwest in the proxy-derived mode as well as the present-day analysis.The only difference emerges over the Peruvian Andes where the MCEOF2 loadings are weak, yet the precipitation data projects negatively onto the present-day SACZ mode.
The interpretation of MCEOF2 as isolating isotopic variability related to SACZ shifts is consistent with previous studies (Novello et al., 2018;Campos et al., 2019).While records located close to the node of the dipole (SBE3+SMT5, TM0) have shown little variation over the LM (Wortham et al., 2017;Novello et al., 2018), records located to the southwest (CR1, ALHO6, Boto) are anticorrelated to those records in the northeast (MV, DV2), i.e. central-eastern Brazilian records are becoming more depleted (enriched) as southwestern records become more enriched (depleted) in heavy isotopes (Novello et al., 2018).The similarity between the correlation of the leading EOF of OLR variability over the SACZ domain with local station precipitation and the isotope-derived MCEOF pattern strongly supports existing evidence of a precipitation amount effect driving isotopic variability in this region (Moquet et al., 2016;Novello et al., 2018).High moisture flux into this region during enhanced SACZ activity leads to the local decrease of d 18 O anomalies recorded in proxy archives.
Convective activity in the SACZ domain is fueled by moisture transported from the SALLJ and anticyclonic low-level flow originating over the tropical South Atlantic, converging with cold fronts entrained in the sub-tropical westerly jet.A more southerly SACZ location reflects a more north-south oriented SALLJ, invigorating convection along the foothills of the eastern Bolivian Andes and could explain why the Boto record projects so strongly onto this mode.At the same time, the southern location of the SACZ also results in stronger interaction with extratropical Rossby waves emanating north from the westerly jet (Kodama, 1992).
SACZ activity within the dipole structure suggested by MCEOF2 , and underscored by the precipitation dependence on OLR, is influenced bya function of the SAMSSASM strength.Neilsen et al. (2019) document that the northernmost position of the SACZ is coincident with the onset and demise of the SAMSSASM, while a more southern position is consistent with the peak phase of the SAMSSASM.Though these results represent an intraseasonal perspective, the dynamics are consistent with an enhanced monsoon corresponding to the peak phase of monsoon activity corresponding to enhanced monsoon activity.Hence, aA strengthened monsoon is accompanied by a shift of the SACZ rainfall to the southwest of its mean position and a decrease of rainfall along the northeastern SACZ margin, accompanied by the consequent depletion change inof isotopic values co-located with shifts in convectionrainfall (Novello et al., 2018).The regional loading pattern in MCEOF2 is indicative of aThe results here show a southwestern shift of the SACZ during periods of enhanced monsoon activityin the LIA rather than a widening of the convective zone as has been previously suggested (Campos et al., 2019).This interpretation agrees with the dominant view in the literature that during the LIA, intensification of the SAMSSASM was accompanied by reduced (enhanced) precipitation along the northeastern (southwestern) edge of the SACZ (Novello et al., 2012(Novello et al., , 2018;;Deininger et al., 2017).Coincident with an intensification of the SAMSSASM core driving the SACZ to higher latitudes, this displacement has been shown on seasonal timescales to be linked to reduced SACZ convective intensity (Barros et al., 2000).This dynamic mechanism may reduce the magnitude of amount-effect driven δ¹⁸O distillation measured in proxies found in southeastern Brazil.
Though MCEOF2 is orthogonal to MCEOF1, the dominant shift in the mean state of PC2 and a similar shift in PC1 during the LIA period reflect a departure toward decreasing values of δ 18 O, indicating a period of stronger Rayleigh distillation and fractionation in the SACZ region driven by a stronger SAMS.This enhancement also corresponds to a period when the ITCZ was displaced south of its mean location, and bringing additionalresulting in enhanced moisture influx and convergence overinto the monsoon region (Lechleitner et al., 2017).This was therefore a period when the SAMSSASM was enhanced overall and both the ITCZ and the SACZ were displaced to the south of their mean locations.While the monsoon core loadings are positively correlated with PC1, central-eastern Brazil proxy loadings (MV, SBE3+SMT5, TM0, DV2) in the northern center of the SACZ dipole are negatively correlated with PC2, and thus during the LIA become relatively enriched in heavy isotopes compared to their mean state.The change in isotopic composition captured by proxies along the northeastern edge of the SACZ indicates a reduction in the local rainfall.The opposite is exhibited by the proxy loadings to the southwest, which are positively correlated with PC2.
In addition to being physically linked to Amazonian convection and moisture transport by the SALLJ, SACZ displacements have also been shown to be modulated by a variety of other external changes on different timescales.Changes in the SST of the South Atlantic Subtropical Dipole have influenced rainfall in northeastern Brazil and southeastern South America on millennial timescales in the past (Wainer et al., 2014) and also on intraseasonal timescales in the present-day climate (Chaves and Nobre, 2004).Other drivers of SACZ variability are linked to changes in anti-cyclonic flow associated with latitudinal shifts of the SASA (Zilli et al., 2019), increases in baroclinicity driven by subtropical incursions of transient synoptic systems from higher latitudes enhancing frontal activity in the region (Kodama, 1992;Novello et al., 2018), and propagating southern hemisphere Rossby wave trains (Gelbrecht et al., 2018).

Model-EOF mode 1
The leading EOFs derived from the isotope-enabled GCMs show remarkable consistency with those extracted from the proxy network, lending confidence to the use of the current generation of climate models in analyzing the dynamics giving rise to these spatial patterns of isotopic variability.There are, however, some discrepancies between the leading proxy-and model-based modes.The model-derived spatial patterns show a larger region of positive correlations associated with the first PC (Figs. 3a-b) relative to the proxy-derived loading pattern.In regions of complex terrain, the model grid resolution of ~2° resolution may limit detailed comparison of proxy-derived and model-derived loadings, though the emergent regional coherence is not anticipated to be sensitive to these regional features.The largest discrepancies are found, especially for records located along the margin of the SAMSSASM domain (PAR, Boto).Those records are more sensitive to large-scale circulation changes and related non-monsoonal influences outside the mature monsoon season (DJF).The PAR record, at the northern margin of the SAMSSASM, currently receives <40% of its annual rainfall in DJF and is influenced by convective precipitation in the monsoon onset and demise phases prior to and after the DJF season.Furthermore, this site may also be influenced by past displacements in the convective margin that characterizes the boundary between Amazonian convection and NE Brazil subsidence in EOF 1 (Wang et al., 2017).The loadings of the Boto record, located on the eastern slopes of the Bolivian Andes along the exit region of the SALLJ, also show some discrepancies with the model simulations.
Even though this location today receives precipitation almost exclusively during the peak monsoon season, it is located at a distal site of the SAMSSASM.and can show considerable spatial heterogeneities in relation to records from the core of the Andes region.This feature has been observed in other reconstructions (e.g.Apaéstegui et al., 2018;Jara et al., 2020;Kock et al., 2020), which argue for the influence of the South Atlantic Ocean delivering moisture to the continent and influencing precipitation patterns over this region in South America.A change in regional circulation modifying moisture transport could be driven by increased Atlantic sea surface temperatures that weaken the SACZ and promote the northward intrusion of extratropical cold fronts through the la Plata basin (Campos et al., 2019).Small changes in the meridional extent of largeand meso-scale circulation systems, such as extratropical cold air incursions propagating north along the eastern Andes, have the potential to drastically change moisture source contributions and the isotopic signal (Apaéstegui et al., 2018).Such changes, especially in a region of complex topography and distal to the monsoon core, are difficult to resolve in a GCMthe coherency analysis and may explain some of the differences observed between the proxy and model EOF loadings in this region.

Model-EOF mode 2
The second model-derived mode represents SACZ variability, consistent with the findings of the proxy-based MCEOF2.The model results present a large loading center, while the proxy-based MCEOF2 displays some indication of a NE-SW dipole.The NE-SW shift of SACZ activity seen in MCEOF2 is driven by periods of weaker/stronger SAMSSASM activity.Though the SACZ is an intrinsic component of the monsoon system, its behavior will vary in responseresponds to both external forcings and regional internal variability that modify the background climate.The lack of this dipole activity in the model-based EOF2 may point toward an inability of the models to correctly simulate the transient behavior of the SACZ, although iCESM appears to correctly simulate the southward displacement of the SACZ during the LIA (see section 4.3).

Model Principal Component characteristics
The temporal characteristics represented in the proxy PC time series characteristics  are not well replicated in the model-based EOF analysis (not shown).The models, particularly iCESM, exhibit principal component time series that are dominated by high-frequency variability around a stationary mean state, but lack power at lower frequencies documented in the proxy data.The proxy network-derived PCs indicate robust changes in the mean state, persisting for decades or centuries, but the models poorly reproduce the observed amplitude or persistence of the monsoon response to internal variability or external forcing.This suggests missing or misrepresented feedbacks within the climate models, which are key to this sustained transient monsoon response.While diagnosing this performance is beyond the scope of the analysis presented herein, it is consistent with similar conclusions derived in previous analyses of the simulated SAMSSASM response to external forcing during the LM (e.g., Rojas et al., 2016).

Pseudoproxy Experiments
The fidelity of the two PPEs in reproducing the main modes of South American monsoon variability based on fourteen heterogeneously distributed time series strengthensenhances our confidence that the monsoon modes are a robust feature of South American summertime climate, not merely apparent from the spatially continuous climate model analysis.
Furthermore, this also confirms the robustness of the proxy MCEOF results as a regional climate signal, not merelyrather than being representative an emergent feature of the noise inherent to the climate system.Though the PPE accurately represent the SASM modes of variability, the reduction of explained variance in the PPE relative to the proxy-based modes likely indicates that the experimental SNRs applied to this region were too small relative to the real -world climate system.
The larger explained variance of the proxy-based modes relative to those from the pseudoproxy network of SNR = 0.5 suggests that the noise characteristics of the SASM proxy network may be very region-and proxy-specific, and that isotopic proxy networks in some cases may exhibit less noise than currently established within the literature for other multiproxy networks (Mann et al., 2007;Smerdon, 2012).More detailed estimates of environmental SNRs would open the possibility for closer calibration between proxy network analysis and estimation of regional climate variability via pseudoproxy experiments.

MCA/LIA changes
To further illustrate the application of isotope-enabled climate models for diagnosing past climate dynamics, we analyze spatial anomalies of three diagnostic variables during the MCA and LIA periods relative to the full LM.The centennial-scale departures from the mean state during these periods are observed in both individual records across the monsoon basin and extracted from the network by means of the MCEOF analysis as well as documented in the CESM LM ensemble experiments (Otto-Bliesner et al., 2016).Here we employ model diagnostics from the iLME full-forcing simulation ensemble average (n=3) to examine dynamical changes associated with hydroclimatic perturbations observed during these time periods.During the MCA, upward motion in the mid-troposphere (500 hPa) is reduced over the tropical Andes, Peruvian Amazon, central Brazil, and in the monsoon trough region (South Atlantic Convergence Zone) (Fig. 5a).
Anomalies during the LIA are roughly opposite, with enhanced vertical motion over the Amazon basin and central Brazil and increased rising motion along the southwestern margin of the South Atlantic Convergence Zone (Fig. 5b).These results support the hypothesis of a weak (strong) SASM during the MCA (LIA) and an corresponding adjustment shifts in the downstream isotopic records.and a shift in tThe SACZ axis was displaceddipolea to the southwesterly SACZ shift during the LIA (Novello et al., 2018) rather than experiencing a stationary enhancement of convective activity (Novello et al., 2018).LIA give more insight into the isotopic variability of the SASM during these periods than the model simulations.The MCA was a period of overall positive δ 18 O anomalies as seen in both the model and paleorecords, especially over the western Amazon, Andes, and downstream monsoon records.The Boto record stands out as a regional anomaly characterized by more negative values during the MCA (Fig. 5c).As discussed above, the isotopic signal at this distal site provides a challenge for the model to capture correctly, highlighting the value of increased observational studies to better understand regional variability in isotopic composition.The LIA period on the other hand shows a clear east-west dipole of d 18 O anomalies in both the proxies and model (Fig. 5d), reminiscent of EOF1 and consistent with the large departure in PC1 during this time.
Strong isotopic enrichment of heavy isotopes is observed over northeastern Brazil and associated with enhanced subsidence and reduced precipitation.Negative isotopic anomalies over the Andes are associated with enhanced upstream convection and high rainfall rates, while depleted isotopic anomalies in the downstream monsoon region are associated with enhanced local convection (Fig. 5b).Notably, significant precipitation anomalies during the MCA and LIA (Figs. 5e-f) are primarily located in the northern (upstream) monsoon region and are coincident with changes in vertical motion but not local δ 18 O changes.Instead, the isotopic signals are displaced downstream of the largest changes seen in vertical motion.This supports our interpretation of changes in oxygen isotopes in tropical precipitation being dominantly controlled by processes associated with deep tropical convective activity upstream of the paleorecord archives.While this study does not decompose convective processes at play, satellite observations have identified air mass trajectory history, continental precipitation recycling, and convective intensity as key processes linked to deep convection that influence δ¹⁸O of precipitation in the austral summer (Samuels-Crow et al., 2014;Vimeux et al., 2021).

Conclusions
We employed a MCEOF decomposition to extract the main modes of hydroclimate variability that are shared across a network of oxygen isotope proxy records and separate them from local environmental factors that are unique to individual sites.The same decomposition was performed for isotope-enabled climate models, which were found to that share the same modes as extracted from the oxygen isotope network.The two main modes of variability are interpreted in the context of SAMSSASM dynamics over the LM, and are characterized by broad regional coherence across the proxy network.The physical underpinnings of the modes are explained by analogs that represent well-understood dynamic components of the present-day South American climate.This newly developed understanding of proxies as integrating regional dynamics promotes caution when interpreting individual records as representing changes in precipitation amount, which may be an mode as a single center of activity rather than a dipole as shown in the paleorecord, it does correctly simulate a more southerly position of the SACZ during the LIA as compared to the MCA, consistent with proxy evidence (Novello et al., 2018).
Notwithstanding the spatial stability of these modes, the temporal evolution of proxy-derived PCs is characterized by clear multidecadal to centennial-scale departures from the mean state.The monsoon is shown to be responsive to external forcings and internal variability, though partitioning the attribution of low-frequency variability between particular drivers has not been parsed out.While the low-frequency regional response to internal or external forcing is evident in the paleoclimate records, it is notably absent in model simulations.More detailed diagnosis of how feedback processes are represented at a regional scale in climate models will therefore be important future work for improving model simulations of monsoonal and regional processes in models.This may, in part, result from sub-grid scale processes that are not well represented in climate models.Indeed, Rojas et al. (2016) have documented that while PMIP ensemble members correctly simulate circulation changes in the LM, corresponding precipitation changes are not recorded well, suggesting that the implementation of some microphysics and convective parameterization schemes may require further study.Nonetheless, the lack of low-frequency variations in the model-derived PCs suggests that a more detailed diagnosis of how feedback processes are represented at a regional scale in climate models will be important future work for improving model simulations of monsoonal and regional processes in models.Additionally, the generation of new high-resolution paleorecords, particularly from the eastern slopes of the Andes along the trajectory of the SALLJ near Boto, would help to further document multi-centennial changes in South America and in the SACZ region over the LM.
Our results are in general agreement with earlier work by Campos et al. (2019), although here we show that the comparison of isotopic analyses with modern continental modes of precipitation as in Campos et al. (2019) is complicated by the fact that a multitude of factors besides precipitation govern isotopic variability across South America.The simplified interpretation of negative (positive) d 18 O anomalies directly corresponding to enhanced (reduced) precipitation across the continent does not take into account the richness of theintegrated isotopic processes and atmospheric dynamics to which isotopic analysisproxies can provide a wealth of insight.Further study with both modeling and observational frameworks would help to identify processes associated with deep convection and sub-grid scale convective activity driving isotopic variability captured in paleoclimate records.
It is also noteworthy that the main modes of variability do not appear to be linked to ENSO or Pacific SST forcing in general.ENSO leads to an in-phase precipitation and isotopic response between the tropical Andes and NE Brazil (drying and isotopic enrichment during the ENSO warm phase), while the response over SE Brazil is of the opposite phase (wetting and negative d 18 O anomalies).This spatial response to ENSO is inconsistent with both MCEOF 1 and MCEOF 2. While the teleconnection mechanisms linking ENSO with tropical South American hydroclimate may not have been stationary throughout the LM, a more likely explanation involves the fact that the ENSO power spectrum resides mostly in the interannual band, which is not sufficiently resolved in our analysis.In observations and in annually resolved and precisely dated archives, the ENSO signal can easily be identified in the isotopic composition of precipitation (Vuille and Werner, 2005;Hurley et al., 2019).These studies show that the influence of ENSO is indirect, via modulation of the monsoon mean state, rather than directly affecting the precipitation or temperature at the proxy sites.Here we focused instead on multidecadal to centennial-scale changes in the monsoon, such as occurred between the MCA and LIA.Hence oOur results do not imply that ENSO was irrelevant for tropical South American climate variability during the LM, but suggest instead that mean state changes in ENSO were not the dominant driver of centennial-scale monsoon changes in the LM.
By reconstructing the variability of these modes through the LM, we have established a baseline envelope of variability against which future change can be compared.The shared presence and spatial characteristics of these modes of isotopic variability in paleorecords and isotope-enabled climate models (both full-forcing and single-forcing experiments) affirms that their existence is intrinsic to the SAMSSASM itself and results from stable boundary conditions such as the shape and latitudinal extent of the continent, Andean orography in particular, rather than as a consequence of external forcings.
Notwithstanding the spatial stability of these modes, their temporal evolution is characterized by clear multidecadal to centennial-scale departures from the mean state.The monsoon is shown to be responsive to external forcings and internal variability, though partitioning the attribution of low-frequency variability between particular drivers as not been parsed out.
While the low-frequency regional response to internal or external forcing is evident in the paleoclimate records, it is notably absent in model simulations.More detailed diagnosis of how feedback processes are represented at a regional scale in climate models will therefore be important future work for improving model simulations of monsoonal and regional processes in models.

Figure 1 :
Figure 1: Oxygen isotope proxy network spanning the South American Monsoon domain (details presented in Table 1).Color shading represents the percent of annual precipitation falling in December to February (DJF, austral summer) based on GPCP data (1979-2019).Vectors represent the low-level (850 hPa) DJF wind field based on ERA5 data (1979-2019).Major dynamical features are highlighted.

Fig. 2 .
The magnitude of the MCEOF Lpoint loadings are weighted byin the spatial patterns represent the correlation coefficient between the median of the Monte Carlo ensemble of the individual proxy records at those locations and the corresponding principal component (PC) time series for that mode.

Figure 2 :
Figure 2: MCEOF decomposition of the oxygen isotope proxy network.a) MCEOF1.Spatial pattern corresponds to ensemble median point loading correlations with the PC time seriesThe magnitude of the point loadings in the spatial patterns represent the correlation coefficient between the median of the Monte Carlo ensemble of proxy records at those locations and the corresponding principal component (PC) time series for that mode; b) PC1.Blue Sshading indicates the total envelope of ensemble variability; c) same as in a) but for MCEOF2; d) same as in b) but for PC2.Pink/blue shading in b) and d) highlight the MCA/LIA periods, respectively.
in the full-forcing simulations with GISS and iCESM shows the same monsoon footprint as present in the proxy-based MCEOF1 albeit with the region of positive correlations associated with the first PC (Figs. 3a-b) encompassing a broader area,.The explained variance of the proxy-based MCEOF1 is matched by the GISS model and exceeded by the iCESM EOF1 (33% in iCESM versus 24% in the proxy data).The first mode across all experiments (full-forcing and single forcings) represents the region characterized by convective activity of the monsoon core, reproducing well the proxy-based MCEOF1.In iCESM, the highest loadings are located over the central Amazon, while in the GISS model, they are centered farther west, over the eastern flank of the Andes.An antiphased region over northeastern Brazil is present in the iCESM experiments, consistent with the proxy-based MCEOF1 and indicative of the BH-NL-induced east-west precipitation and

Figure 3 :
Figure 3: Comparison of model-derived modes (shading) overlaid with proxy-derived modes (circles) leading modes model-derived and of isotopic variability.Color gradients correspondshading indicates to the correlation of the mode EOF loadings pattern with the corresponding PC time series.Models plotted are a) EOF1 from GISS full forcing experiment (shading) overlaid with proxybased MCEOF1 (colored circles), , b) as in a) but for iCESM full forcing experiment, c) as in a) but for iCESM volcanic forcing experiment, d) as in a) but for iCESM solar irradiance experiment, e) as in a) but for iCESM orbital forcing experiment, f) iCESM PPE EOF1 derived from grid cells corresponding to proxy locations only, with white noise added (SNRsignal-to-noise ratio = 0.5), g) as in f) but for SNR signal-to-noise ratio = 0.25.h-n) as in a-e) but for EOF2.

Figure 4 :
Figure 4: MCEOF patterns and their corresponding dynamic interpretations.a) MCEOF1 loading patterns; b) upper-level (200 hPa) DJF streamlines from ERA5 reanalysis (1979-2019); c) MCEOF2 loading patterns; d) spatial correlation of DJF OLR field with the PC1 of OLR variability derived from the SACZ region (black box).Colored circles indicate the correlation of local precipitation derived from station data with PC1 of the OLR decomposition over the SACZ region.

Figure 5 :
Figure 5: iCESM full-forcing multi-model mean anomalies during MCA (850-1160 CE; left column) and LIA (1489-1850 CE; right column) relative to the full LM (850-1850 CE) mean in austral summer (DJF) for a-b) mid-tropospheric vertical velocity (500 hPa, Pa/s, note that negative anomalies signify enhanced upward motion / reduced subsidence), c-d) precipitation-weighted δ 18 O (permil; colored circles represent values from proxy MC 1000-member ensemble mean; note that color scales differ between proxy and model estimates) and e-f) precipitation (mm/day).Stippling indicates regions where anomalies are significant at p ≤ 0.1.Isotopic changes during the MCA and LIA are simulated surprisingly well in iCESM as far as the observed spatial pattern across SA are concerned (Figs.5c-d), although the amplitude of the anomalies is lower in iCESM compared to the paleorecord (circles overlaid on the model-simulated contours).The anomalies seen in the paleorecord during the MCA and oversimplification for regions of the monsoon where the model-based MCA/LIA analysis suggests strong influence of upstream effectsconvective processes rather than local controls.The first mode of variability corresponds to an uppertropospheric Rossby wave response to condensational heating over the Amazon basin during the mature phase of the monsoon, exciting or relaxing the BH-NL system.The enhancement (weakening) of the South American Monsoon during the LIA (MCA) is captured by the leading PC1 from the first MCEOF of the proxy records, resulting in modulation of the E-W precipitation dipole across the continent.While the spatial characteristics of this mode are quite accurately captured by climate models, they do not simulate the multi-centennial shifts in the mean state seen in the isotopic composition of the proxy records.The second mode of variability corresponds to changes in intensity and location of the monsoon convective axis (SACZ).The NE-SW shift in the SACZ captured by the MCEOF2 reaffirms the prevailing theory that the changes in monsoon strength are reflected by displacements in the mean SACZ position and intensity.While the model represents this