The Influence of Volcanic Eruptions on the Climate of Tropical South America during the Last Millennium in an Isotope-enabled General Circulation Model

Currently, little is known on how volcanic eruptions impact large-scale climate phenomena such as South American paleo-Intertropical Convergence Zone (ITCZ) position and summer monsoon behavior. In this paper, an analysis of observations and model simulations is employed to assess the influence of large volcanic eruptions on the climate of tropical South America. This problem is first considered for historically recent volcanic episodes for which more observations are available but where fewer events exist and the confounding effects of El Niño–Southern Oscillation (ENSO) lead to inconclusive interpretation of the impact of volcanic eruptions at the continental scale. Therefore, we also examine a greater number of reconstructed volcanic events for the period 850 CE to present that are incorporated into the NASA GISS ModelE2-R simulation of the last millennium. An advantage of this model is its ability to explicitly track water isotopologues throughout the hydrologic cycle and simulating the isotopic imprint following a large eruption. This effectively removes a degree of uncertainty associated with error-prone conversion of isotopic signals into climate variables, and allows for a direct comparison between GISS simulations and paleoclimate proxy records. Our analysis reveals that both precipitation and oxygen isotope variability respond with a distinct seasonal and spatial structure across tropical South America following an eruption. During austral winter, the heavy oxygen isotope in precipitation is enriched, likely due to reduced moisture convergence in the ITCZ domain and reduced rainfall over northern South America. During austral summer, however, more negative values of the precipitation isotopic composition are simulated over Amazonia, despite reductions in rainfall, suggesting that the isotopic response is not a simple function of the " amount effect ". During the South American monsoon season, the amplitude of the temperature response to volcanic forcing is larger than the rather weak and spatially less coherent precipitation signal, complicating the isotopic response to changes in the hydrologic cycle.


Volcanic forcing on climate
Plinian (large, explosive) volcanic eruptions are a dominant driver of naturally forced climate variability during the last millennium (LM, taken here to be 850 CE to present; e.g., Stothers and Rampino, 1983;Hansen et al., 1992;Crowley et al., 2000;Robock et al., 2000Robock et al., , 2003;;Goosse et al., 2005;Yoshimori et al., 2005;Emile-Geay et al., 2008;Cole-Dai, 2010;Timmreck, 2012;Iles et al., 2013;Schurer et al., 2014).In addition to their importance for 20th century climate, they are the largest magnitude external forcing during last 1000 years of the preindustrial period, the most recent key interval identified by the Paleoclimate Modelling Intercomparison Project Phase III (PMIP3).As such, these eruptions serve as a natural test bed to assess the skill of climate models in simulating how climate responds to external perturbations.
Although the most significant climate impacts of eruptions are realized over just a few years following the eruption, they provide the source of the largest amplitude perturbations Published by Copernicus Publications on behalf of the European Geosciences Union.Aerosol optical depth (AOD) used to force the NASA GISS ModelE2-R over the last millennium and (bottom) zoomed in on the period 1950-1999 (Sato et al., 1993;Crowley and Unterman, 2013) as discussed in the text.AOD is the vertically integrated (15-35 km) and latitudinal average from 30 • S to 30 • N. Note difference in vertical scale between graphs.Orange dashed line marks the AOD threshold for defining a LM eruption in the present study.Eruption events defined in the text must sustain the threshold AOD for at least 1 year, so not all events above the orange dashed line are used in the composites.
to Earth's energy budget during the LM.For example, the eruption of Mt.Pinatubo in June 1991, although transitory, exerted a radiative forcing comparable to an instantaneous halving of atmospheric CO 2 (Hansen et al., 1992;Minnis et al., 1993; see also Driscoll et al., 2012, for models in the Coupled Model Intercomparison Project Phase 5 (CMIP5)); several paleo-eruptions during the LM likely had an even larger global impact (Fig. 1).
The principle climate impact from volcanic eruptions results from the liberation of subsurface sulfur-containing gases such as sulfur dioxide, which are injected into the stratosphere and react with water to form sulfate aerosols (e.g., Harshvardhan and Cess, 1976;Coakley and Grams, 1976;Pollack et al., 1976Pollack et al., , 1981;;Lacis et al., 1992).The most pronounced impact of large tropical eruptions includes a radiatively cooled troposphere and heated stratosphere (e.g., Lacis et al., 1992;Robock and Mao, 1995;Stenchikov et al., 1998).Sulfate aerosols from the Mt.Pinatubo eruption grew from a background effective radius of ∼ 0.2 µm up to ∼ 0.8 µm, strongly scattering incoming solar radiation.For sulfate aerosols in this size range, this shortwave scattering is 5-10 times larger than the increase in infrared opacity from the aerosols, and results in a warming stratosphere and cooling of Earth's surface (Turco et al., 1982;Lacis et al., 1992).
Studies on the impacts of volcanic eruptions have generally focused on global or Northern Hemisphere metrics (e.g., Lucht et al., 2002;Gillett et al., 2004;Shindell et al., 2004;Oman et al., 2005Oman et al., , 2006;;Anchukaitis et al., 2010;Peng et al., 2010;Evan et al., 2012;Zhang et al., 2013;Man et al., 2014;Stoffel et al., 2015), for instance in examining responses to the East Asian monsoon or the Arctic Oscillation (e.g., Ortega et al., 2015).Comparatively little attention has been given to the Southern Hemisphere, or to South America specifically (although see Joseph andZeng, 2011, andWilmes et al., 2012).Some previous work has focused on the Southern Annular Mode in the ERA-40 and NCEP/NCAR reanalysis, in addition to a previous version of NASA Goddard Institute for Space Studies (GISS) Model-E (Robock et al., 2007) and in a subset of CMIP3 models (Karpechko et al., 2010) or in CMIP5 (Gillett and Fyfe, 2013).
How volcanic forcing is expressed over South America remains an important target question for several reasons.First, recognition of the South American monsoon system (SAMS) as an actual monsoon system is less than two decades old (Zhou and Lau, 1998), and thus study of SAMS dynamics is still relatively young (Sect.1.3) and very little work has been done specifically focused on volcanic eruptions.For instance, should we expect to see a reduction in austral summer rainfall?Secondly, the largest volcanic eruptions during the late 20th century (e.g., Mt. Agung, 1963, Indonesia;El Chichón, 1982, Mexico;Mt. Pinatubo, 1991, island of Luzon in the Philippines -hereafter, these three events are referred to as L20 eruptions) occur quasi-simultaneously with an anomalous El Niño-Southern Oscillation (ENSO) state, and in general represent a small sample size in a noisy system.This limits the prospect of robust hypothesis testing and guidance for what impacts ought to be expected following large eruptions at the continental scale.Finally, South America offers promise for a comparatively dense network of highresolution proxy locations relative to other tropical regions (see below), offering the potential to detect whether South American hydroclimate signals to large eruptions are borne out paleoclimatically.
In this study, we explore the post-volcanic response of South American climate operating through the vehicle of unique model simulations (spanning the LM) using the recently developed GISS ModelE2-R (LeGrande et al., 2016;Schmidt et al., 2014a), which allows for the sampling of a greater number of events than is possible over the instrumental period.Emphasis is placed on temperature and precipitation, but a novel part of this study extends to the response of water isotopologues (e.g., H 18 2 O) (colloquially referred to hereafter as "isotopes" and expressed as δ 18 O in units per mil (‰) vs. Vienna Standard Mean Ocean Water).The isotopic composition of precipitation (δ 18 O p ) is a key variable that is directly derived from proxy data used in tropical paleoclimate reconstructions.
The aim of this paper is to create a potentially falsifiable prediction for the isotopic imprint that a volcanic eruption should tend to produce across the South American continent.The ability to explicitly model the isotopic response allows for a less ambiguous comparison of simulations and paleoclimate records and for hypothesis testing.It is unclear whether or not the current proxy archives are suitable to test such a prediction with high confidence, given dating uncertainties (in both proxies and in the actual timing of eruptions), or the level of noise in proxy data and the real world.Additionally, the prevailing high-resolution archives in South America only feature a few tropical records (Vimeux et al., 2009;Neukom and Gergis, 2012;Vuille et al., 2012).Nonetheless, the growing number of high-resolution records offers hope that testing the modeled response to high-frequency volcanic signals will be an avenue for future research.This can also better inform debate centered on the inverse problem in interpreting isotopic signals (i.e., what do observed changes in proxy data imply about past climate changes?), which remains contentious (Sect. 1.4).
The structure of this article is as follows.In the remaining part of Sect. 1, we summarize previous literature on the impact of large volcanic eruptions on paleoclimate, in addition to a discussion of South American climate.Section 2 presents data and methodology, including how volcanic forcing is implemented in ModelE2-R.Section 3 discusses our results and we end with conclusions in Sect. 4.
Two volcanic forcing data sets (Gao et al., 2008;Crowley and Unterman, 2013) relying on ice core reconstructions of volcanism are used as input in the LM ModelE2-R simulations (and are the CMIP5/PMIP3 LM standard), as discussed in Sect. 2.

Tropical South American Climate
South America is home to nearly 390 million people.The continent spans a vast meridional extent (from ∼ 10 • N to 55 • S), contains the world's largest rainforest (the Amazon), in addition to one of the driest locations on Earth (the Atacama Desert).The continent has diverse orography, spanning the high Andes along the Pacific to Laguna del Carbón in Argentina, the lowest point in the Southern Hemisphere.Because of this, South America hosts a rich diversity of climate zones and biodiversity, all of which may respond in unique ways to external forcing.
The most prominent climatic feature of tropical and subtropical South America is the South American monsoon system (Zhou and Lau, 1998;Marengo et al., 2001Marengo et al., , 2012;;Vera et al., 2006;Garreaud et al., 2009).Much of South America is in a monsoon regime, with tropical/subtropical rainfall over the continent exhibiting a pronounced seasonal cycle.Unlike other monsoon systems such as that in Asia, low-level easterly winds prevail during the entire year in tropical South America, although the wind anomalies do change direction when the annual mean wind field is removed from winter and summer composites (Zhou and Lau, 1998).
During austral winter, the maximum in continental precipitation is largely restricted to north of the Equator, in a bandlike pattern associated with the oceanic Intertropical Convergence Zone (ITCZ).During austral summer, convection is displaced from northwestern South America, and a band of heavy precipitation covers much of the continent, from the southern Amazon Basin to central Brazil and northern Argentina.A distinctive feature of the SAMS is the South Atlantic Convergence Zone (SACZ), a band of cloudiness and precipitation sourced primarily from the tropical Atlantic that extends diagonally (southeastward) from the Amazon towards southeastern Brazil (Fig. 2).
The SAMS onset occurs around the end of October and the demise between the end of March and April (e.g., Nogués-Paegle et al., 2002;Vera et al., 2006;Silva and Carvalho, 2007).The dominant mode of intraseasonal precipitation variability over South America during summer exhibits a dipole pattern (Nogués-Paegle and Mo, 1997), seesawing between the SACZ region and southeastern South America, the latter including the densely populated La Plata Basin with local economies strongly dependent on agricultural activities.
The SAMS is strongly modulated by ENSO behavior on interannual timescales (Vuille and Werner, 2005;Garreaud et al., 2009).In general, SAMS-affected regions of tropical South America tend to experience drier than normal conditions during El Niño, while conditions in subtropical latitudes are anomalously humid, including the southeastern part of the continent.Surface air temperatures tend to be anomalously warm in tropical and subtropical South America during El Niño events.These relationships depend somewhat on the time of year, and during La Niña events, the pattern is essentially reversed.

Recent South American monsoon reconstructions from isotopic proxies
SAMS variability spanning most of the Holocene has been diagnosed from speleothem records in the Peruvian Andes (Kanner et al., 2013) and a review focused on the last 1000- 2000 years was given in Bird et al. (2011) and Vuille et al. (2012).In all cases, a critical piece of information that is required to properly diagnose paleo-SAMS variability is the ability to translate oxygen isotope variability from natural recorders into a physical climate signal of interest.
Early work on isotopes in ice core records from the tropical Andes detected a Little Ice Age (LIA) signal in the oxygen isotope composition of the ice, with results initially interpreted to reflect variations in local temperature due to their resemblance to ice core records from Greenland (e.g., Thompson et al., 1995Thompson et al., , 1998) ) and due to their isotopic enrichment over the past 150 years, in parallel with rising global mean temperatures (Thompson et al., 2006).A temperature-dependence to oxygen isotope variability has been long known and is particularly important in mid-to high latitudes (Dansgaard, 1964) and is most directly related to the ratio of initial and final water vapor content of a parcel that is transported horizontally, rather than the temperaturedependence of fractionation itself (Hoffman and Heimann, 1997).
This interpretation in the tropics has been challenged through a number of observational and modeling efforts (Hardy et al., 2003;Vuille and Werner 2005;Vimeux et al., 2005Vimeux et al., , 2009;;Kanner et al., 2012) which suggest that the isotopic signal is more closely related to the degree of rainout upstream in regions of intense convection (in the case of South America, over the Amazon Basin).Additionally, since sea surface temperatures (SST) in the Pacific have a large influence on SAMS intensity on interannual timescales in the present, oxygen isotope variability over much of tropical South America is linked to the state of the equatorial Pacific (Bradley et al., 2003;Vuille et al., 2003a, b).
In regimes that are highly convective in nature as in tropical South America, empirical evidence shows that the amount of precipitation (the so-called "amount effect", Dansgaard, 1964) rather than the condensation temperature correlates most strongly with δ 18 O p variability, at least on seasonal to interannual timescales.In reality, however, the rainout most relevant for the oxygen isotope signal may be at a significant distance from the site where the proxy is derived, potentially complicating the use of local calibrations to climatology as a guide for δ 18 O p interpretations (Schmidt et al., 2007).Isotopic concentrations are explainable as being a function of isotopic concentration of evaporative fluxes, rainout along the moisture transport path, and mixing.
The influence of precipitation amount on δ 18 O p , in addition to changes in the partitioning of precipitation sources, has also been identified on decadal to orbital timescales through speleothem records and lake sediments (Cruz et al., 2005;Van Breukelen et al., 2008;Bird et al., 2011;Kanner et al., 2012).These studies have also highlighted the role of latitudinal displacements of the ITCZ, which is ultimately the main moisture conduit for precipitation over the South American continent.Furthermore, many records collected throughout South America now provide evidence for enriched δ 18 O p values during the Medieval Climate Anomaly, which is indicative of weakened SAMS convection and rainout, followed by depleted δ 18 O p values, suggesting heavier rainfall during the LIA in tropical South America (Bird et al., 2011;Apaéstegui et al., 2014) with an opposite response in Northeast Brazil (Novello et al., 2012).This, in turn, has been interpreted in terms of North Atlantic SST anomalies (Vuille et al., 2012;Ledru et al., 2013) and the position of the Atlantic ITCZ.
Nonetheless, oxygen isotopes respond in unique ways depending on the climate forcing of interest.Indeed, a unique, quantitative local relationship between an isotope record and any particular climate variable of interest is unlikely to hold for all timescales and prospective forcing agents (Schmidt et al., 2007) thus motivating the use of forward modeling to work in conjunction with proxy-based field data.For the re-mainder of this paper, we focus specifically on the volcanic forcing response.

Data
The primary tool used in this study is the water-isotopeenabled GISS ModelE2-R.ModelE2-R is a fully coupled atmosphere-ocean general circulation model (GCM) (LeGrande et al., 2016;Schmidt et al., 2014a) that explicitly tracks stable water isotopes.The version used here is the same as the non-interactive atmospheric composition (NINT) physics version used in the CMIP5 experiments (Miller et al., 2014).The current model features 2 • latitude × 2.5 • longitude horizontal resolution and 40 vertical levels in the atmosphere up to 0.1 hPa, and is coupled to the Russell Ocean that conserves heat, water mass, and salt (Russell et al., 1995) at 1 • × 1.25 • resolution with 32 vertical levels.ModelE2-R includes stratospheric dynamics and prescribed ozone and aerosol species.
Due to uncertainties in past radiative forcing, a suite of LM simulations using ModelE2-R have been run with different combinations of plausible solar, volcanic, and anthropogenic land-use histories (Schmidt et al., 2011(Schmidt et al., , 2012) ) but with identical greenhouse gas and orbital evolution.These simulations span the period 850-2005 CE.There are two reconstructions of past volcanic activity (Gao et al., 2008;Crowley and Unterman, 2013) that are used in six combinations of the ModelE2-R LM simulations (see the "past1000" experimental design at http://data.giss.nasa.gov/modelE/ar5/,and below).
For the LM, three forcing combinations are available in the GISS ModelE2-R simulations that use the Crowley reconstruction for volcanic perturbations.These include Pongratz et al. (2008;land)/Krivova et al. (2007;solar), Kaplan et al. (2010;land)/Krivova et al. (2007;solar), and Pongratz et al. (2008;land)/Steinhilber et al. (2009;solar) (see Schmidt et al., 2011Schmidt et al., , 2012)).We focus only on results from the Crowley reconstruction prior to 1850 CE due to a mis-scaling of the Gao forcing in the model that roughly doubled the appropriate radiative forcing.For the historical period (1850present), the volcanic forcing history is based on Sato et al. (1993) and is equivalent among the different (six) simulation members.
Crowley and Unterman (2013) discuss the details behind the LM aerosol optical depth (AOD) reconstruction that defines the volcanic forcing time series in ModelE2-R (Fig. 1).This estimate is derived from sulfate peaks in ice cores, which are relatively well dated and referenced to the historical record during the satellite era.Crowley and Unterman (2013) provide an AOD history over four latitude bands (from 0 to 30 and 30 to 90 • in both hemispheres).ModelE2-R uses a cubic spline to interpolate this forcing data set over 24 latitude bands.The choice of volcanic eruptions used for the LM analysis (Sect.2.2 below) is based on the AOD data set from this 24-latitude grid.
We note that there are more recent volcanic reconstructions available (e.g., Sigl et al., 2015) suggesting modifications to the timing or magnitude of LM eruptions, as well as developments of data sets focusing on sulfur injection and microphysics-based evolution of the aerosol forcing (e.g., Arfeuille et al., 2014).In this contribution, we are agnostic concerning the veracity of the forcing data sets that were standard for CMIP5/PMIP3 but stress that timing of eruptions is irrelevant in our modeling context and that the model results should be interpreted as a self-consistent response to the imposed AOD and particle size.
Water isotope tracers are incorporated into the model's atmosphere, land surface, sea ice, and ocean.These isotopes are advected and tracked through every stage of the hydrologic cycle.At each phase change (including precipitation, evaporation, ice formation or melting) an appropriate fractionation factor is applied (Schmidt et al., 2005) and all freshwater fluxes are tagged isotopically.Stable isotope results from the lineage of GISS models have a long history of being tested against observations and proxy records (e.g., Schmidt et al., 2007;LeGrande andSchmidt, 2008, 2009;Lewis et al., 2010Lewis et al., , 2013Lewis et al., , 2014;;Field et al., 2014).
In addition to the model, we briefly explore the observed instrumental record to assess responses to eruptions occurring during the 20th century.To do this, we take advantage of the NASA GISS Surface Temperature analysis (GISTEMP) land-ocean index (Hansen et al., 1999), and Global Precipitation Climatology Centre (GPCC) v6, a monthly precipitation data set over land (Schneider et al., 2011).For Figs. 2 and 3m where ocean climatological data are shown, we use the Global Precipitation Climatology Project (GPCP) version 2.2 (Adler et al., 2003), a combined land station and satellite product available since 1979.These data sets are called upon to gauge the tropical climate response following the three L20 eruptions.We use the 2.5 • resolution GPCC data set, as that is comparable to the GISS model and what is justified by the station coverage in this part of the world.The GPCC product offers considerably better global and South American coverage than other precipitation data sets, although observational density for rainfall is still considerably more problematic over South America than for many other regions of the globe.There is a sharp drop-off in the number of rain gauge stations used earlier in the 20th century over much of the South American continent.Figure S1 shows the station density at the time of each L20 eruption, as well as the total number of land stations over South America with time.
Finally, in Sect.3.1 we present data from the Global Network of Isotopes in Precipitation (GNIP) accessible from the International Atomic Energy Agency (IAEA) for δ 18 O p as a test of the model's ability to track the seasonal hydrologic cycle in the form of its isotopic response over South America before discussing the LM results.Unfortunately, there is considerable spatial and temporal heterogeneity in the GNIP data over South America.In fact, only a few stations have data overlap with one or two eruptions and with a sufficient number of δ 18 O p data points to establish reasonable seasonal or annual statistics.Additionally, the post-volcanic (L20) anomalous isotope field over South America strongly resembles the ENSO expression on the isotope field (Vuille et al., 2003a) and with large spread between events (not shown).This suggests that internal variability (ENSO) dominates the forced (volcanic) response in this very small historical sample size, thereby leaving little hope that the prevailing network of observations is suitable for hypothesis testing and model validation in our context.

Superposed epoch and composite analysis
We present the spatial pattern of observed and simulated response for temperature and precipitation over land for two L20 eruptions (El Chichón and Mt.Pinatubo).Results are shown for annual means in 1983 and 1992.We choose only two for brevity, as our argument that assessing the signal in any specific region is difficult in a small sample of eruptions is unaffected.Because of the dominant influence of unforced variability on tropical South American climate (Garreaud et al., 2009) overriding the volcanic signal during the L20 eruptions, we instead present a superposed epoch anomaly composite of the tropical-mean temperature anomaly, zonally averaged from 30 • S to 30 • N. Results are shown for years −3 to +5, with zero defining the time of the eruption.This composite is formed for all three L20 eruptions.In all cases, the 5 years prior to the eruption were subtracted from the superposed composite.Other sensible choices for the non-eruption reference period do not significantly change the results.
For the full LM spatial composites, we use only eruptions where vertically integrated (15 to 35 km) stratospheric AOD averaged from 30 • N to 30 • S exceeds 0.1 for at least 12 consecutive months in the simulation (top panel in Fig. 1).For the LM composites, we focus only on seasonal (DJF and JJA) composites, and a given season will enter the composite if at least 2 out of 3 months meet the AOD threshold; this criterion yields 15 eruptions since 850 CE.The selection of events used in the LM composite is very weakly sensitive to this choice of latitude band.Mt.Pinatubo is the only L20 eruption in this composite, and is actually one of the smallest eruptions in this selection based on the maximum AOD encountered near the time of the eruption (see Table 1 for dates of each event).We believe sampling a larger number of events with greater forcing is a better way to understand the volcanic response in this model, rather than increasing the ensemble size for the L20 events.
For the LM "non-eruption" fields, we use 15 years prior to the eruption as a reference period to calculate the anomaly for each event, unless another event occurs during that time (overlap occurs only once for eruptions in 1809 and 1815), in which case the pre-1809 climatology is used twice.The exception is for Mt.Pinatubo, which again uses the previous 5 years to calculate the anomaly.When constructing seasonal averages of δ 18 O p , the oxygen isotope value for each month is weighted by the precipitation amount during that month, at each grid cell.
Since each post-eruption difference field is computed using the immediate response minus a local 15-year climatology, time is not relevant in this analysis and so we use all three members with the Crowley forcing (representing over 3000 years of simulation time) to generate a composite that features 45 volcanic "events" (15 eruptions in each of the three members).In the historical (post-1850) extension of these runs, the coding error that resulted in a misimplementation of the Gao forcing is not an issue, and so we use six ensemble members each (three volcanic events in six ensemble members) for the L20 results.
The ensemble-mean composite results displayed for the LM eruptions include contributions from three members that differ not just in the internal variability but also in their solar and land-use forcing.Similarly, the L20 results are from model runs that also include other transient historical forcings occurring at the time of the eruption, including greenhouse gas increases throughout the duration of the event (although these forcings are the same among all ensemble members).However, in all cases we focus only on the immediate years after the eruption.Since the primary signal of interests is expected to be large compared to the impact of more slowly varying and smaller-amplitude forcings, the ensemble spread for a given eruption can be interpreted as a sampling of the model internal variability coincident with the event.We have tested our composite results using the same dates as our volcanic events in simulations with other varying forcings with no volcanoes (there are no volcano-only runs with this model version for the LM), and the results are indistinguishable from noise (not shown).The LM composite results are discussed in Sect.3.2.
Finally, it is now well appreciated that any climate response under investigation will be shackled to the spatial structure of the forcing imposed on a model.For example, preferential heating/cooling of one hemisphere will induce different tropical precipitation responses than a well-mixed gas that behaves CO 2 -like (Kang et al., 2008(Kang et al., , 2009;;Frierson and Hwang, 2012;Haywood et al., 2013).Figures S2  and S3 in the Supplement show the latitudinal AOD distribution structure for all eruptions used in the generation of the LM composites within ModelE2-R.The mean of all events is rather symmetric between hemispheres (though somewhat skewed toward the Southern Hemisphere tropics).Thus, the resulting climate responses outlined in this paper ought to be viewed as a response consistent with a forcing that is relatively symmetric about the Equator.

Influence of ENSO on the late 20th century (L20) eruptions
For the L20 volcanic events, El Niño events are occurring quasi-simultaneously with the eruption.This introduces a pervasive issue when attempting to isolate the volcanic signal (e.g., Robock, 2003;Trenberth and Dai, 2007;Joseph and Zeng, 2011) and is particularly important over South America (e.g., Garreaud et al., 2009).
In order to remove the effects of ENSO from the superposed epoch and spatial composite analyses described above in the GISTEMP and GPCC data, we first perform a multiple regression with the variable of interest over the period 1951-2005 using a linear time trend and the Niño 3 www.clim-past.net/12/961/2016/Clim.Past, 12, 961-979, 2016 index as predictors (5 • N-5 • S, 150 • W-90 • W; data from http://www.cpc.ncep.noaa.gov/data/indices/using ERSSTv4 data) over the same period, excluding 2 years of data after each L20 eruption.At each grid cell, the Niño 3 index is lagged from 0 to 6 months and the correlation coefficient with the maximum absolute value (since a positive index can induce a negative anomaly in the variable of interest) is found.This is similar to the approach used in Joseph and Zeng (2011), allowing the maximum ENSO influence to be removed at each grid point at the expense of contemporaneous relationships.The lagged Niño index is then regressed against the time series of each variable and the residual from this regression is retained.This approach assumes a linear relationship between ENSO and the climate response over South America, an assumption that appears justified on interannual to decadal timescales (Garreaud et al., 2009).
For each of the six ensemble members used in the model L20 composite, a similar procedure is performed in which the Niño 3 index (consistent with the realization of the Niño 3 domain SSTs in that model simulation) is calculated and regressed out in the same manner.For the full LM computations, the number of larger-amplitude events in the threeensemble member composite should help average out the influence of Pacific SST variability, and no ENSO removal procedure is applied.

L20
Figure 3 illustrates that ModelE2-R reproduces the seasonal cycle of climatological rainfall (comparing Fig. 3a with b) and oxygen isotope distribution (comparing Fig. 3c with  d) with some fidelity over South America.This includes a meridional migration of the ITCZ toward the summer hemisphere and an intensification of the South American monsoon during DJF.Where data permit (Fig. 3c), there is good agreement between model and observations, both displaying oxygen isotope DJF enrichment relative to JJA in the tropics north of the Equator and the higher latitudes south of 30 • S, and depletion in the continental interior south of the Equator associated with the monsoon wet season.ModelE2-R (Fig. 3b) tends to produce too much precipitation over northeastern Brazil, although the gross features of the seasonal migration in rainfall are well captured.This ability to accurately simulate the seasonality of δ 18 O p over the tropical Americas has also been noted in two atmospheric GCMs with no coupled ocean (NASA-GISS II and ECHAM-4; see Vuille et al., 2003a).The model also provides a skillful representation of δ 18 O p in response to ENSO (not shown), including increased δ18O p values over tropical South America in response to El Niño (Vuille and Werner, 2005).
Figure 4 shows the ENSO-removed superposed epoch analysis for tropical temperature associated with the recent three L20 eruptions.There is good agreement between the observed and modeled temperature response, both in amplitude and recovery timescale.The tropical-mean cooling is on the order of several tenths of a degree, and larger after Mt.Pinatubo (not shown individually).
The spatial structure of the post-El Chichón and Pinatubo events in land observations and the individual model realizations are shown in Figs. 5 and 6, respectively.Observations exhibit cooling over much of the globe, especially after Mt.Pinatubo, that is largely reproduced by the model.However, there is considerable spread among the individual ensemble members and between the two events, indicating a large role for internal variability in dictating the observed spatial pattern following these events.This is also true over South America.
In GISTEMP, the high latitudes of South America cool more than the tropical region of the continent after Mt.Pinatubo.There is still a residual signal from ENSO in tropical South America following both L20 eruptions that is not reproduced by the model.This is not unexpected, since ENSO events comparable to the magnitude of the historic  realizations do not occur coincident with the volcanic forcing in the individual ensemble members.The magnitude of signal is sensitive to the Niño index used in the regression method described above.Without ENSO removal, tropical South America warms following the two eruptions (not shown).The influence of ENSO appears minimal over the higher-latitude sectors of the continent.
The precipitation pattern following the L20 eruptions exhibits substantial variability in space and across eruptions, with a general drying pattern over land in tropical latitudes.South America experiences less precipitation near the Equator after Mt.Pinatubo (see also Trenberth and Dai, 2007), a pattern reproduced in some of the ensemble realizations.
It should be noted that model-observation comparison is hindered not just by internal variability but also by the spec- ified historical volcanic forcing in the model.In fact, the Stratospheric Aerosol and Gas Experiment (or SAGE) II satellite sensor was saturated by the aerosol cloud after Mt.Pinatubo; subsequent work (Santer et al., 2014;Schmidt et al., 2014c) suggests that the forcing following Pinatubo is too large in the CMIP5 generation of models.
Because of the considerable variability seen in observations (following historical eruptions) and also across ensemble members, it is evident that a larger signal-to-noise ratio than is available from the L20 eruptions alone is required to help isolate any volcanic signal.ModelE2-R is the laboratory from which we proceed to sample a larger number of events, some of which contain larger amplitude than the L20 eruptions.

Temperature and precipitation
Figure 7 shows the LM post-volcanic temperature composite for all 45 events.During both seasons, cooling is statistically significant over virtually the entire continent (stippling indicates significance at the 90 % level, t test).The temperature response is strongest in the interior of the continent, particularly during the austral winter.The enhanced high-latitude cooling exhibited in the observations after Mt.Pinatubo does not emerge in the model composite.
The precipitation anomalies for the LM composite are shown in Fig. 8.As expected, there is a distinct seasonal structure in the response, with the largest anomaly concentrated in a narrow region north of the Equator during austral winter, coincident with the location of climatological rainfall maxima in the region.During JJA, precipitation increases in the North Atlantic region following volcanic eruptions, while very strong and statistically significant precipitation reductions occur just north of the Equator (including over northern Brazil, Ecuador, Venezuela, Colombia, and Guyana) and encompassing the northern Amazon Basin.This signal is consistent with a weakening of the moisture flux owing to the decrease in saturation vapor pressure due to cooling that is demanded by the Clausius-Clapeyron relationship (Held and Soden, 2006).During this season, the precipitation response is significant virtually everywhere in northern South America.Figure S5 further illustrates that the JJA precipitation response is remarkably robust to all eruptions that enter into the composite.
Figure 9b illustrates the relationship between areaaveraged precipitation from 20 to 0 • , 77.5 to 45 • W (for DJF) and 0 to 10 • N, 77.5 to 52.5 • W (for JJA) and the maximum AOD encountered for each eruption.These two regions were selected to reflect the seasonal migration of rainfall (Figs. 2, 3).Fifteen eruptions are displayed with the threemember ensemble spread given for each.Precipitation only increases north of the Equator during austral winter in a few model realizations.Moreover, the magnitude of the precipitation response during JJA scales with the size of the eruption, particularly for very large eruptions (e.g., comparing five eruptions with AOD > 0.3 vs. those with smaller perturbations, although the spread amongst the ensemble members is large).The spatial composite for each individual eruption (each averaged over the three ensemble members) is shown in Fig. S5.The precipitation response during austral summer is more difficult to interpret (Fig. 8a).During this season, the zonally oriented Atlantic ITCZ migrates southward and the SACZ becomes more intense as it is connected with the area of convection over the central and southeastern part of the continent.It is noteworthy that the land cools substantially more than the surrounding ocean (Fig. 7), which one could expect to weaken the monsoon-sourced precipitation during DJF.While precipitation is indeed reduced over the tropical continent, the response is weaker than in JJA and less spatially coherent, with many areas failing to meet statistical significance.An analysis of the individual responses reveals that the signal is more eruption-dependent during DJF than during JJA (see Fig. S4), with a few events actually exhibiting modest increases in precipitation.Nonetheless, there is a clear tendency for reduced DJF precipitation within the SAMS region, although there is little to no dependence of the mean rainfall anomaly on the magnitude of the AOD perturbation, at least above the 0.1 threshold used in this study (Fig. 9b), unlike for equatorial South America during JJA.Conversely, the temperature response (Fig. 9a) depends on the size of the eruption in both seasons, as is expected given its dependence on the size of the radiative forcing.

Tropical hydroclimate response
Since the South American climate is intimately linked to large-scale tropical dynamics, the global precipitation composite is shown in Fig. S6 to better inform the model response.The most robust signal is characterized by a reduction in tropically averaged precipitation and the tendency for wet regions to become drier, and dry regions to become wetter (see also Iles et al., 2013;Iles and Hegerl, 2014), in contrast to the anticipated hydrologic response in a future, higher-CO 2 world (Held and Soden, 2006).
This pattern is a thermodynamic effect linked to reduced moisture convergence within the convergence zones and to reduced moisture divergence in the descending zones of the Hadley cell, which reduces the contrast in values of precipitation minus evaporation (P − E) between moisture convergence and divergence regions (Chou et al., 2009).The complete hydrologic response of the P − E field (not shown) has the same spatial structure as the P field, since evaporation is decreasing nearly everywhere in the tropics.Because both P and E are decreasing on the equatorward flank of the ITCZ, the P − E signal is rather weak in the deep tropics, while P − E increases more rapidly than P in the subtropics.
The tendency for modest precipitation anomalies over the continent during DJF appears to be part of a pattern that spans a broad swath of longitudes across the entire deep tropics in association with the seasonal cycle.Nonetheless, the response during DJF is weaker over land.

Oxygen isotope anomalies
In order to relate the responses discussed in the previous sections back to a potentially observable paleoclimate metric, we show the composite δ 18 O p field for the DJF and JJA seasons in South America (Fig. 10).It should be cautioned that much of the isotopic variability that can be observed in proxies within the continental interior or high-elevation  glacier sites will likely be seasonally biased toward the wetseason months (Hardy et al., 2003).
During the JJA season, there is a strong enrichment of the δ 18 O p pattern that is zonally extended over equatorial South America.In addition, there is a corresponding δ 18 O p depletion in the adjacent North Atlantic sector.This response is inextricably coincident with the strong change in precipitation in the ITCZ domain that was assessed in Fig. 8, and is broadly consistent with a "rainfall amount" control on the isotopic imprint (Dansgaard, 1964).South of approximately 15 • S, the sign of the anomaly reverses to a depletion of the heavy isotope.
During the austral summer, volcanic eruptions lead to a clear negative excursion in δ 18 O p over virtually the entire SAMS region, including the Amazon Basin, tropical Andes, and eastern Brazil.The statistical significance of the resulting isotopic anomaly extends throughout most of the landmass within the tropics and in the North Atlantic.There are small but non-significant exceptions (positive δ 18 O p excursions) such as in eastern Brazil.The negative excursions also include regions outside of the SAMS belt in the subtropics and mid-to high latitudes of South America.
The austral summer δ 18 O p depletion is the opposite sign from what one would expect if the reduced precipitation were driving the isotopic response.Thus, it may well be that the strong temperature response to volcanic eruptions dominates the continent-wide oxygen isotope depletion during the DJF season and in the extratropics during JJA over the relatively weak precipitation response.Precipitation, on the other hand, appears to be the primary control knob of δ 18 O p during JJA within the ITCZ region.
The correlation between δ 18 O p and temperature or precipitation, based on a regression using all 45 volcanic events, is reported in Fig. 9, using the same domains for DJF and JJA described in Sect.3.2.1.In the case of volcanic forcing it appears that the amplitude of the temperature response to volcanic eruptions over tropical South America is larger than the rather weak and spatially incoherent precipitation signal, although both the temperature and precipitation coefficients must be considered to characterize the isotopic variability during this season (Fig. S7).This may explain why the DJF isotopic signal related to volcanic eruptions seems to respond to atmospheric cooling, even in the tropics, where isotopic variability is usually more closely associated with changes in the hydrologic cycle.During JJA, the isotopic enrichment is much more closely associated with precipitation reduction north of the Equator, whereas the JJA δ 18 O p -temperature relationship is weak and non-significant.
Taken together, these results suggest that the primary controls on oxygen isotope variability may vary by forcing agent, rather than being determined inherently by the latitude of interest (e.g., "precipitation-driven" in the tropics and "temperature-driven" in the extratropics).This conclusion is compelled by the fact that the precipitation production and distribution in proxy records are the result of an interaction between multiple scales of motion in the atmosphere, the temperature of air in which the condensate was embedded, and exchange processes operating from source to sink of the parcel deposited at a site.Thus, a consistent description of how to interpret oxygen isotopes into a useful climate signal cannot be given without considering all of these processes and the target process of interest.
To further complement the spatial analysis, a composite Hovmöller diagram is utilized (Fig. 11) in order to illustrate the time evolution of the temperature, precipitation, and oxygen isotope response.For this plot, the start of each eruption is defined as the closest January to the first month in which AOD reaches 0.1 in order to illustrate the seasonal evolution (rather than compositing by "month from each eruption" as in Fig. 3).Therefore, for all 45 events in the composite, the local AOD may reach this threshold within 5 months (before or after) of the January baseline point (eruptions in June are rounded up to the following January).The Hovmöller composites are plotted for 10 years (beginning January 3 years prior to the eruption).The closest January point to the start of each eruption occurs in the 37th month of the Hovmöller (solid black line in Fig. 11a, b, d).Results are zonally averaged from 77.5 to 45 • W.
Figure 11a demonstrates a substantial temperature anomaly that peaks south of 10 • S (compare also to Fig. 7).The cooling lasts for several years following the eruption, and decays gradually until most of the signal is lost (∼ 4 years after the eruption in the South American sector), but remains 0.1-0.2• C colder than the pre-eruption climatology.The zonally averaged peak reductions in South American precipitation anomalies occur over the tropical latitudes and last for a comparable period of time as the maximum temperature response.The precipitation anomaly itself migrates synchronously with the seasonal cycle (red line in Fig. 11c maps out the latitude of maximum climatological precipitation averaged over all 15-year climatologies of each 45member event, as a function of time of year).Figure 11b indicates that the largest precipitation response is confined to the equatorial and northern regions during JJA, with weak protrusion into higher tropical latitudes only 1-2 years after the eruption.The JJA isotopic enrichment in northern South America lasts for two seasons in our composite, while there is sustained isotopic depletion during DJF in the SAMS region for about 3 years (Fig. 11d).
Figure 12 provides additional statistical insight into the magnitude of the excursions described in this section.Here, we sampled 100 random 45-event composites in a control simulation with no external forcing (each "event", two seasons in length, is defined as an anomaly expressed relative to a pre-eruption climatology as done previously).The anomalies were averaged over the same areas as in Fig. 9, with different domains for DJF and JJA.Notably, for both seasons and for all three variables examined, the single 45-event postvolcanic composite (purple square) lies outside the distribution of all sampled 45-event composites constructed with no external forcing.Nonetheless, the distribution for a smaller sample of events (black circles denote the data for each (15) eruption, each averaged over the three ensemble members) shows considerable spread.
The δ 18 O p anomalies discussed above result from changes in the isotopic content of precipitation, which may be due to changes in precipitation amount or to other changes in the isotopic composition of the water vapor that condensed to form the precipitate.The changes are not determined by changes in the seasonality of the precipitation.To illustrate this (Fig. S8), we decomposed the δ 18 O p field (see Liu and Battisti, 2015) by weighting the monthly oxygen isotope field by the pre-eruption precipitation values.The results are indistinguishable from the total δ 18 O p field, suggesting that any changes in monsoon seasonality are negligible in contributing to the isotopic signal, unlike the orbital case considered in Liu and Battisti (2015).

Conclusions
In this study, we have analyzed the response of temperature, precipitation, and δ 18 O p over South America to volcanic forcing associated with large tropical eruptions during the last millennium.It is now well known that volcanic eruptions lead to large-scale cooling throughout the tropics, and this result extends to most of the South American continent as well, except in regions that may be simultaneously affected by opposing ENSO behavior.In general, the precipitation response has been more enigmatic, though our results are in broad agreement with numerous other studies showing that there is a substantial decline in tropical-mean precipitation.
However, the immediate post-volcanic impact over South America has a complex seasonal and spatial structure.During the austral winter, the precipitation response over the continent is slaved to the response of the large-scale circulation, including a weakening of rainfall intensity within the ITCZ that is migrating northward.In the extratropics, the continent cools and exhibits slight precipitation declines nearly everywhere.Our results suggest the seasonal monsoon precipitation (during DJF) in ModelE2-R exhibits a fairly weak response that is scattered across the continent.It appears that volcanic forcing preconditions the tropical rainfall over the continent to decline during the wet season but that this response is likely to be eruption-dependent and may be overwhelmed by internal variability.
A unique aspect of this study was to probe the δ 18 O p response to volcanic eruptions.During JJA, the precipitation isotopic composition is less negative in northern South America as convective activity produces substantially less precipitation.No such relation was found during the monsoon season, even within the tropics, where the large cooling appears to lead to more depleted δ 18 O p , despite a weakened hydrologic cycle and reduced monsoon precipitation.In the extratropics, it appears that the temperature decline is driving isotopes toward more depleted values.
Unfortunately, validation of our model results is hindered by the paucity of observational stable isotope data and by the coincidence of volcanic eruptions with ENSO events over the 20th century.Nonetheless, our results may provide some guidance in the search of volcanic signals in highresolution isotopic or other temperature-and precipitationsensitive proxy data from South America.Given the importance of volcanic forcing for climate variability over the past millennium, and in particular the LIA period, which has been identified as a period of significant climatic perturbation in isotopic proxies from South America, a better understanding of the climatic response to volcanic forcing over this region is urgently needed.
Figure1.Aerosol optical depth (AOD) used to force the NASA GISS ModelE2-R over the last millennium and (bottom) zoomed in on the period 1950-1999(Sato et al., 1993;Crowley and Unterman, 2013) as discussed in the text.AOD is the vertically integrated (15-35 km) and latitudinal average from 30 • S to 30 • N. Note difference in vertical scale between graphs.Orange dashed line marks the AOD threshold for defining a LM eruption in the present study.Eruption events defined in the text must sustain the threshold AOD for at least 1 year, so not all events above the orange dashed line are used in the composites.

CFigure 2 .
Figure 2. Top: observed Climatological Precipitation for DJF (shading, in mm day −1 ).SAMS box is drawn over the domain from 20 to 0 • , 77.5 to 45 • W and used for Figs. 9 and 12. Data from the GPCP product, long-term climatological rainfall derived from years 1981-2010.Bottom: as above, except for JJA.Box from 0 to 10 • N, 77.5 to 52.5 • W used in averaging for Figs. 9 and 12.

Figure 3 .
Figure 3. Seasonal cycle (DJF minus JJA) of precipitation in the (a) GPCP precipitation product, from data in Fig. 2, (b) in ModelE2-R, (c) δ 18 O p in GNIP data, and (d) δ 18 O p in ModelE2-R.GNIP data only shown for stations with at least 90 reported δ 18 O p values at a given station from 1960 to present, in addition to at least 10 data values for each month: December, January, February, June, July, and August.Stations with seasonal differences of less than ±1.0 ‰ are also omitted in panel (c).

Figure 4 .
Figure 4. Composite tropical (30 • S to 30 • N) temperature response following the L20 volcanic eruptions.Fill colors denote observed monthly anomalies using GISTEMP, with 24-month running average shown as a solid black line.ModelE2-R ensemble mean is shown as solid orange line and dashed grey lines indicate the six individual ensemble members.Anomalies are referenced to 5 years prior to eruptions (years −5 to 0).Dashed purple lines encompass the 5-95 % interval for monthly tropical-mean temperature anomalies (relative to the previous 5-year mean) in the GISTEMP product from 1950 to present.The calculation of this range omits data 2 years after the L20 eruptions.The range is not symmetric about zero due to the tropical warming trend during this interval.All data are based on the ENSO-removal technique discussed in the text.

Figure 5 .
Figure 5. Annual-mean temperature change ( • C, ocean masked) for each L20 eruption (labeled on plot) in GISTEMP (top row) and each ModelE2-R ensemble member, as discussed in the text.All plots use ENSO-removal procedure described in the text.

Figure 7 .
Figure 7.Last millennium post-volcanic temperature composite ( • C) averaged over all 45 events during (a) DJF and (b) JJA from GISS ModelE2-R using procedure described in the text.Stippling highlights areas with anomalies significant at p<0.1.

Figure 8 .
Figure 8.Last millennium post-volcanic precipitation composite (mm day −1 ) with all eruption events during (a) DJF and (b) JJA from GISS ModelE2-R using procedure described in the text.Stippling highlights areas with anomalies significant at p<0.1.

Figure 9 .
Figure 9. (a) Average temperature anomaly during DJF within the SAMS region (red; 20 to 0 • S, 77.5 to 45 • W) and equatorial South America during JJA (blue; 0 to 10 • N, 77.5 to 52.5 • W) plotted against the peak AOD for all 15 eruptions (each point averaged over three ensemble members with the three-member spread shown as horizontal bars).(b) As in (a), but for precipitation.Dashed horizontal lines indicate the 5-95 % range for a two-season average of each variable relative to the previous 15 years (averaged over the same domain) in the entire control simulation with no external forcing.In both panels, the correlation coefficient and p value are reported for (a) temperature and (b) precipitation vs. δ 18 O p in each season and over the same domain.The regression uses all 45 volcanic events.

Figure 10 .Figure 12 .
Figure 10.Last millennium post-volcanic oxygen isotope in precipitation ( δ 18 O p ) composite (per mil) with all eruption events during (a) DJF and (b) JJA from GISS ModelE2-R using procedure described in the text.Stippling highlights areas with anomalies significant at p<0.1.

Table 1 .
Time of eruptions and global aerosol optical depth (AOD) from Crowley and Unterman (2013).List of eruptions used in study.List of LM and L20 eruptions.Start of eruption dates based on when they can be identified in the Crowley/Sato time series averaged over the latitude band from 30 • S to 30 • N. May be slightly different than actual eruption date.b Maximum AOD over the 30 • S to 30 • N latitude band encountered in monthly time series during the duration of each event.
a c December in year prior to listed date.d Mt.Agung and El Chichón included in L20 but not LM composites.