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

The South American Summer Monsoon (SASM) 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. We extract the coherent modes of variability of the SASM over the last millennium (LM) using a Monte Carlo empirical orthogonal function (MCEOF) decomposition of 14 δ¹⁸O proxy records and compare them with modes decomposed from isotope-enabled climate model data. The two leading modes reflect the isotopic variability associated with (1) thermodynamic changes driving the upper-tropospheric monsoon circulation (Bolivian High–Nordeste Low waveguide) and (2) the latitudinal displacement of the South Atlantic Convergence Zone (SACZ). The spatial characteristics of these modes appear to be robust features of the LM hydroclimate over South America and are reproduced both in the proxy data and in isotope-enabled climate models, regardless of the nature of the model-imposed external forcing. The proxy data document that the SASM was characterized by considerable temporal variability throughout the LM, with significant departures from the mean state during both the Medieval Climate Anomaly (MCA) and the Little Ice Age (LIA). Model analyses during these periods suggest that the local isotopic composition of precipitation is primarily a reflection of upstream rainout processes associated with monsoon convection. Model and proxy data both point to an intensification of the monsoon during the LIA over the central and western parts of tropical South America and indicate a displacement of the South Atlantic Convergence Zone (SACZ) to the southwest. These centennial-scale changes in monsoon intensity over the LM are underestimated in climate models, complicating the attribution of changes on these timescales to specific forcings and pointing toward areas of important model development.

1 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 (SASM; 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 20th 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 SASM dynamics. The onset of the SASM 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, over 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 (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 SASM 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 SASM. 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).

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, across 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 SASM 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 lack 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 SASM 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 SASM using a proxy network derived solely from oxygen isotope records covering the LM, all located within the SASM 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.

Table 1 Oxygen isotope proxy network.

2 Data and methods

2.1 MCEOF analysis approach

The variability of the SASM was investigated by using an empirical orthogonal function (EOF) decomposition of oxygen isotope data from a network of paleoclimate records and isotope-enabled climate model simulations. The EOF decomposition allows for discrete aspects of monsoon variability to be isolated and displayed as leading modes. A Monte Carlo EOF (MCEOF) was applied to the network of paleorecords, with the Monte Carlo resampling ensuring that uncertainties associated with the age models of the individual paleoclimate records were 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 (North et al., 1982) and the fact that they aligned with a dynamically and physically plausible interpretation.

The 14 records used in this analysis were derived from raw samples from varied archives of speleothems (n=12), an ice core (n=1), and authigenic calcite in a lake sediment core (n=1). These records were identified as PAR (composed of PAR01 + PAR03), PAL (composed of PAL3 + PAL4), Pumacocha, P00-H1, MV (composed of MV1 + MV3), PIM4, DV2, SBE3 + SMT5, Quelccaya, ALHO6, TM0, Boto (composed of Boto1 + Boto3 + Boto7 + Boto10), JAR (composed of JAR1 + JAR4), and CR1. For each sample, 1000-member ensembles of age models were developed based on a Monte Carlo resampling from a Gaussian distribution of the 1-sigma measurement uncertainty of each age tie point. 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 merged to form composite records. Oxygen isotope time series were resampled to annual resolution and truncated to 850–1850 CE, resulting in a 1000-member ensemble of 14 complete records spanning the LM. Each of the 1000 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).

2.2 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 have been published, with the exception of the MV1 and MV30 samples, which were 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/paleo-search/, last access: 30 January 2022) or through the University of São Paulo Institute of Geosciences. The 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 establish record age chronologies. 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). 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 consistent with data in the original publications (indicated in Table 1), and samples with hiatuses were interpolated to ensure a continuous record for analysis.

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.

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).

2.3 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 30 years before and after the gap. This procedure bridges a gap in the data in order to ensure isotopic time series have common time steps when applying the EOF analysis while also maintaining the intrinsic characteristics of variability present in a given record. 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 did not have a gap in data, in which case the samples were merged or taken to be continuous, respectively.

2.4 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 CE to cover the LM and linearly interpolated to achieve annual resolution.

2.5 Isotope-enabled model analysis and pseudoproxy experiments

Output from two previously published isotope-enabled fully coupled global climate model (GCM) 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 ocean–atmosphere model (GISS-E2-R; Schmidt et al., 2007, 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 $\mathrm{1.9}{}^{\circ }\times \mathrm{2}{}^{\circ }$ 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 $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ horizontal 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).

The iCESM is known to have a negative isotope bias (Brady et al., 2019), overestimating the strength of the relationship between precipitation and δ¹⁸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 on 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 data sets that are complementary 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) framework (Smerdon, 2012, 2016). In both cases, the EOFs were based on precipitation-weighted δ¹⁸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 and the single-forcing simulations of the iLME using the volcanic, solar, and orbital forcings were used. The ensemble members forced by land use–land cover and greenhouse gas changes 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 ratio 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 these sub-sampled time series, effectively scaling the model-simulated isotopic time series to more realistically reflect 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 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 ratios (SNRs) 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.

2.6 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 the present day (1979–2019; Hersbach et al., 2020). Precipitation station data were extracted from Global Historical Climatology Network (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 4 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 were 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 δ¹⁸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 Student t test and evaluated at p≤0.1.

3 Results

3.1 MCEOF results from proxy network

The results of the proxy-based MCEOF analysis are shown in Fig. 2. The magnitude of the MCEOF 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.

Figure 2 MCEOF decomposition of the oxygen isotope proxy network. (a) MCEOF1. The magnitude of the 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 shading indicates the total envelope of ensemble variability. (c) Same as in (a) but for MCEOF2. (d) Same as in (b) but for PC2. Pink and blue shading in (b) and (d) highlight the MCA and LIA periods, respectively.

The leading mode of isotopic variability explains approximately 24 % of the total variance. This mode is characterized by in-phase isotopic variability 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 regions of the monsoon domain) vary in phase, while the northeastern and 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 SASM 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 (Fig. 2a and 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 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.

3.2 Model results

3.2.1 Continuous domain results

EOF1 in the full-forcing simulations with GISS and iCESM shows the same monsoon footprint as present in the proxy-based MCEOF1. 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- and single-forcing experiments) represents the region characterized by convective activity of the monsoon core, reproducing the proxy-based MCEOF1 well. 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 isotope dipole (Cruz et al., 2009; Cheng et al., 2013). This characteristic aspect of the SASM 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 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 (Fig. 3h and 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.

The EOF decomposition of the iLME single-forcing experiments (Fig. 3c–e and 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 SASM, 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 land mass with the equatorial trade winds.

3.2.2 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. Despite the decrease in explained variance, the spatial patterns of EOF1 and EOF2 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 in 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 loadings in the region of central-eastern Brazil associated with the SACZ (Fig. 3n).

Figure 3 Comparison of model-derived modes (shading) overlaid with proxy-derived modes (circles) of isotopic variability. Color shading indicates the correlation of the EOF loadings with the corresponding PC time series. Models plotted are (a) EOF1 from the GISS full-forcing experiment, (b) as in (a) but for the iCESM full-forcing experiment, (c) as in (a) but for the iCESM volcanic forcing experiment, (d) as in (a) but for the iCESM solar irradiance experiment, (e) as in (a) but for the iCESM orbital forcing experiment, (f) iCESM PPE EOF1with signal-to-noise ratio =0.5, and (g) as in (f) but for signal-to-noise ratio =0.25. (h)–(n) as in (a)–(e) but for EOF2.

4 Discussion

4.1 Proxy-based modes of monsoon variability

4.1.1 MCEOF1

The first mode extracted from the proxy network via the MCEOF decomposition corresponds to a distinct pattern that is inherent to the SASM. MCEOF1 represents the isotopic variability 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).

Figure 4 MCEOF patterns and their corresponding dynamic interpretations. (a) MCEOF1 loadings; (b) upper-level (200 hPa) DJF streamlines from ERA5 reanalysis (1979–2019); (c) MCEOF2 loadings; (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.

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, represents a dynamic response to changes in the hydrologic cycle, which also imprints directly on the oxygen isotopic composition recorded in paleorecords. Thus, variability of oxygen isotopes as captured by the first mode represents a proxy for the monsoon response to changes in internal or external forcing. Multi-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 refers to the intensification 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 isotopic fractionation of water vapor from strong or weak convective periods lead to this precipitation dipole being prominently displayed in stable oxygen isotopes of precipitation (Vuille and Werner, 2005; Sturm et al., 2007).

4.1.2 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) (Fig. 4c and 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 project 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 in δ¹⁸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 orientation of the SALLJ, invigorating convection along the foothills of the eastern Bolivian Andes and possibly explaining 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 is influenced by SASM strength. Neilsen et al. (2019) document that the northernmost position of the SACZ is coincident with the onset and demise of the SASM, while a more southern position is consistent with the peak phase of the SASM. Though these results represent an intraseasonal perspective, the dynamics are consistent with the peak phase of monsoon activity corresponding to enhanced monsoon activity. A strengthened monsoon is accompanied by a shift of the SACZ rainfall to the southwest of its mean position and a decrease in rainfall along the northeastern SACZ margin, accompanied by the consequent change in isotopic values co-located with shifts in rainfall (Novello et al., 2018). The results here show a southwestern shift of the SACZ during periods of enhanced monsoon activity 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 SASM was accompanied by reduced (enhanced) precipitation along the northeastern (southwestern) edge of the SACZ (Novello et al., 2012, 2018; Deininger et al., 2017). Coincident with an intensification of the SASM 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 δ¹⁸O, indicating a period of stronger Rayleigh distillation and fractionation in the SACZ region driven by a stronger SASM. This enhancement also corresponds to a period when the ITCZ was displaced south of its mean location, resulting in enhanced moisture influx and convergence over the monsoon region (Lechleitner et al., 2017; Steinman et al., 2022). 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.

4.2 Model-derived MCEOF discussion

4.2.1 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 (Fig. 3a and 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 for records located along the margin of the SASM 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 SASM, 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 EOF1 (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 SASM 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). Such changes, especially in a region of complex topography and distal to the monsoon core, are difficult to resolve in the coherency analysis and may explain some of the differences observed between the proxy and model EOF loadings in this region.

4.2.2 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 or stronger SASM activity. Though the SACZ is an intrinsic component of the monsoon system, its behavior responds 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 Sect. 4.3).

4.2.3 Model principal component characteristics

The proxy PC time series characteristics (Fig. 2b–d) 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 they 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 SASM response to external forcing during the LM (e.g., Rojas et al., 2016).

4.2.4 Pseudoproxy experiments

The fidelity of the two PPEs in reproducing the main modes of South American monsoon variability based on 14 heterogeneously distributed time series enhances 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, rather than being representative of noise inherent to the climate system. Though the PPEs accurately represent the SASM modes, 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.

4.3 MCA and 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, and 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 adjustment in the downstream isotopic records. The SACZ axis was displaced to the southwest during the LIA rather than experiencing a stationary enhancement of convective activity (Novello et al., 2018).

Figure 5 iCESM full-forcing multi-model mean anomalies during the MCA (850–1160 CE; a, c, e) and LIA (1489–1850 CE; b, d, f) 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 and reduced subsidence), (c, d) precipitation-weighted δ¹⁸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 d⁻¹). Stippling indicates regions where anomalies are significant at p≤0.1.

Isotopic changes during the MCA and LIA are simulated well in iCESM as far as the observed spatial pattern across SA is concerned (Fig. 5c and 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 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 δ¹⁸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 δ¹⁸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 (Fig. 5e and f) are primarily located in the northern (upstream) monsoon region and are coincident with changes in vertical motion but not local δ¹⁸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).

5 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 share the same modes as extracted from the oxygen isotope network. The two main modes of variability are interpreted in the context of SASM 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 oversimplification for regions of the monsoon where the model-based MCA–LIA analysis suggests strong influence of upstream convective processes rather than local controls. The first mode of variability corresponds to an upper-tropospheric 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. 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 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 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. 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) δ¹⁸O anomalies directly corresponding to enhanced (reduced) precipitation across the continent does not take into account the richness of integrated isotopic processes and atmospheric dynamics to which isotopic proxies 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 sea surface temperature 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 δ¹⁸O anomalies). This spatial response to ENSO is inconsistent with both MCEOF1 and MCEOF2. 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 multi-decadal to centennial-scale changes in the monsoon, such as occurred between the MCA and LIA. Our 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) affirm that their existence is intrinsic to the SASM 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.