Hemispherically Asymmetric Volcanic Forcing of Tropical Hydroclimate during the Last Millennium

Volcanic aerosols exert the most important natural radiative forcing of the last millennium. State-of-the-art paleoclimate simulations of this interval are typically forced with diverse spatial patterns of volcanic forcing, leading to different responses in tropical hydroclimate. Recently, theoretical considerations relating the intertropical convergence zone (ITCZ) position to the demands of global energy balance have emerged in the literature, allowing for a connection to be made between the paleoclimate simulations and recent developments in the understanding of ITCZ dynamics. These energetic considerations aid in explaining the well-known historical, paleoclimatic, and modeling evidence that the ITCZ migrates away from the hemisphere that is energetically deficient in response to asymmetric forcing. Here we use two separate general circulation model (GCM) suites of experiments for the last millennium to relate the ITCZ position to asymmetries in prescribed volcanic sulfate aerosols in the stratosphere and related asymmetric radiative forcing. We discuss the ITCZ shift in the context of atmospheric energetics and discuss the ramifications of transient ITCZ migrations for other sensitive indicators of changes in the tropical hydrologic cycle, including global streamflow. For the first time, we also offer insight into the large-scale fingerprint of water isotopologues in precipitation (δ 18 O p) in response to asymmetries in radiative forcing. The ITCZ shifts away from the hemisphere with greater volcanic forcing. Since the isotopic composition of precipitation in the ITCZ is relatively depleted compared to areas outside this zone, this meridional precipitation migration results in a large-scale enrichment (depletion) in the isotopic composition of tropical precipitation in regions the ITCZ moves away from (toward). Our results highlight the need for careful consideration of the spatial structure of volcanic forcing for interpreting volcanic signals in proxy records and therefore in evaluating the skill of Common Era climate model output.


Introduction
The intertropical convergence zone (ITCZ) is the narrow belt of deep convective clouds and strong precipitation that develops in the rising branch of the Hadley circulation. Migrations in the position of the ITCZ have important consequences for local rainfall availability, drought and river discharge, and the distribution of water isotopologues (e.g., δ 18 O and δD, hereafter simply referred to as water isotopes, with notation developed in Sect. 3.3) that are used to derive inferences of past climate change in the tropics.
Meridional displacements of the ITCZ are constrained by requirements of reaching a consistent energy balance on both sides of the ascending branch of the Hadley circulation (e.g., Kang et al., 2008Kang et al., , 2009Schneider et al., 2014). Although the ITCZ is a convergence zone in near-surface meridional mass flux, it is a divergence zone energetically. The stratification of the tropical atmosphere is such that moist static energy (MSE) is greater aloft than near the surface, compelling Hadley cells to transport energy in the direction of their upper tropospheric flow (Neelin and Held, 1987). If the system is perturbed with preferred heating or cooling in one hemisphere, the anomalous circulation that develops resists the resulting asymmetry by transporting energy from the heated to the cooled hemisphere. Conversely, meridional moisture transport in the Hadley circulation is primarily confined to the low-level equatorward flow, so the response of the tropical circulation to asymmetric heating demands an ITCZ migration away from the hemisphere that is energetically deficient. Since the mean circulation dominates the atmospheric energy transport (AET) in the vicinity of the equator, the recognition that the ITCZ is approximately co-located with the latitude where meridional column-integrated energy fluxes vanish has provided a basis for relating the mean ITCZ position to AET. We note that this perspective focused on atmospheric energetics is distinct from one that emphasizes sea surface temperature gradients across the tropics (Maroon et al., 2016).
This energetic framework has emerged as a central paradigm of climate change problems, providing high explanatory and predictive power for ITCZ migrations across timescales and forcing mechanisms McGee et al., 2014;Schneider et al., 2014). It is also a compelling basis for understanding why the climatological annual-mean ITCZ resides in the Northern Hemisphere (NH); it has been shown that this is associated with ocean heat transport, which in the prevailing climate is directed northward across the equator Marshall et al., 2014). The energetic paradigm also predicts an ITCZ response for asymmetric perturbations that arise from remote extratropical forcing. This phenomenon is exhibited in many numerical experiments, is borne out paleoclimatically, and has gradually matured in its theoretical articulation (Chiang and Bitz, 2005;Broccoli et al., 2006;Kang et al., 2008Kang et al., , 2009Brocolli, 2008, 2009;Chiang and Friedman, 2012;Frierson and Hwang, 2012;Bischoff and Schneider, 2014;Adam et al., 2016).
Thus far, however, little or only very recent attention has been given to the relation between transient ITCZ migrations and explosive volcanism (although see Iles et al., 2013;Liu et al., 2016, Sect. 2). This connection has received recent consideration using carbon isotopes in paleo-records (Ridley et al., 2015) or in the context of volcanic and anthropogenic aerosol forcing in the 20th century (Friedman et al., 2013;Allen et al., 2015;Haywood et al., 2015). The purpose of this paper is to use the energetic paradigm as our vehicle for interpreting the climate response in paleoclimate simulations featuring explosive volcanism of varying spatial structure.
Much of the existing literature highlighting the importance of spatial structure in volcanic forcing focuses on the problem of tropical vs. high-latitude eruptions and dynamical ramifications of changing pole-to-equator temperature gradients (e.g., Robock, 2000;Stenchikov et al., 2002;Shindell et al., 2004;Oman et al., 2005Oman et al., , 2006; Kravitz and Robock, 2011), which is a distinct problem from one focused on interhemispheric asymmetries in the volcanic forcing. Further-more, episodes with preferentially higher aerosol loading in the Southern Hemisphere (SH) have received comparatively little attention, probably due to the greater propensity for both natural or anthropogenic aerosol forcing to be skewed toward the NH.
Here we show that it matters greatly over which hemisphere the aerosol loading is concentrated and that this asymmetry in aerosol forcing has a first-order impact on changes in the tropical hydrologic cycle, atmospheric energetics, and the distribution of the isotopic composition of precipitation.

Methods
To illuminate how the spatial structure of volcanic forcing expresses itself in the climate system, we call upon two stateof-the-art models that were run over the preindustrial part of the last millennium, nominally 850-1850 CE (hereafter, LM), the most recent key interval identified by the Paleoclimate Model Intercomparison Project Phase 3 (PMIP3). An analysis of this time period is motivated by the fact that volcanic forcing is the most important radiative perturbation during the LM (LeGrande and Anchukaitis, 2015;Atwood et al., 2016). Furthermore, the available input data that define volcanic forcing in CMIP5/PMIP3 feature a greater sample of events, larger radiative excursions, and richer diversity in their spatial structure than is available over the historical period. This allows for a robust composite analysis to be performed over this interval.
The two general circulation models (GCMs) that we use as our laboratory are NASA GISS ModelE2-R (hereafter, GISS-E2) and the National Center for Atmospheric Research (NCAR) Community Earth System Model (version 1.1) Last Millennium Ensemble (hereafter, just CESM to describe this set of simulations). The GISS-E2 version used here is the same as the noninteractive atmospheric composition physics version used in the CMIP5 initiative (called "NINT" in Miller et al., 2014). CESM is a community resource that became available in 2015  and consists of several component models each representing different aspects of the Earth system; the atmospheric component is the Community Atmosphere Model version 5 (CAM5, see Hurrell et al., 2013), which in CESM features a 1.9 • latitude × 2.5 • longitude horizontal resolution with 30 vertical levels up to ∼ 2 hPa. The GISS-E2 model is run at a comparable horizontal resolution (2 • × 2.5 • ) and with 40 vertical levels up to 0.1 hPa.
Both GISS-E2 and CESM feature multiple ensemble members that include volcanic forcing. There are only a small number of volcanic eruptions in our different forcing classifications (see below) in each 1000-year realization of the LM, motivating an ensemble approach to sample multiple realizations of each eruption. There are currently 18 members in CESM, including 13 with all transient forcings during the LM and five volcano-only simulations. This number is much higher than the number of ensembles used for participating LM simulations in CMIP5/PMIP3. The volcanic reconstruction is based on Gao et al. (2008, hereafter, G08) and the ensemble spread is generated from round-off differences in the initial atmospheric state (∼ 10 −14 • C changes in the temperature field). Sampling many realizations of internal variability is critical in the context of volcanic eruptions given the different trajectories that can arise in the atmosphere-ocean system in response to a similar forcing (Deser et al., 2012).
For GISS-E2, there exist six available members that include a transient volcanic forcing history. Here, however, we use only the three simulations that utilize the G08 reconstruction. This was done in order to composite over the same dates as the CESM events and because the other volcanic forcing dataset that NASA explored in their suite of simulations (Crowley and Unterman, 2013) only provides data over four latitude bands, complicating inferences concerning hemispheric asymmetry. Taken together, there are 21 000 years of simulation time in which to explore the postvolcanic response while probing both initial-condition sensitivity and the structural uncertainty between two different models. The three GISS-E2 members also differ in the combination of transient solar/land-use histories employed, but since our analysis focuses only on the immediate postvolcanic imprint, the impact of these smaller amplitude and slowly varying forcings is very small. We tested this using the composite methodology developed below on no-volcano simulations with other single forcing runs (in CESM) or with combined forcings (in GISS-E2) and found the results to be indistinguishable from that of a control run (not shown).
In both GISS-E2 and CESM, the model response is a slave to the spatial distribution of the imposed radiative forcing, which was based on the aerosol transport model of G08 rather than the coupled model stratospheric wind field, thus losing potential insight into the seasonal dependence of the response that may arise in the real world. For our purpose, however, this is a more appropriate experimental setup, since the spatial structure of the forcing is implicitly known (Fig. 1).
The original G08 dataset provides sulfate aerosol loading from 9 to 30 km (at 0.5 km resolution) for each 10 • latitude belt. This reconstruction is based on sulfate peaks in ice cores and a model of transport that determines the latitudinal, height, and time distribution of the stratospheric aerosol. In CESM, aerosols are treated as a fixed size distribution in three levels of the stratosphere, which provide a radiative effect, including shortwave scattering and longwave absorption. The GISS-E2 model is forced with prescribed aerosol optical depth (AOD) from 15 to 35 km, based on a linear scaling with the G08-derived column volcanic aerosol mass (Stothers, 1984;Schmidt et al., 2011), with a size distribution as a function of AOD as in Sato et al. (1993) -thus altering the relative longwave and shortwave forcing (Lacis et al., 1992;Lacis, 2015). Note that Samalas is omitted, as discussed in text. The time series is at seasonal (5-month) resolution, and thus multiple points may be associated with a single eruption. The hemispheric contrast (NH minus SH) clear-sky net solar radiation (FSNTC -in W m −2 ) in CESM LME is shown in orange (offset to have zero mean).
We note that the GISS-E2 runs forced with the G08 reconstruction in CMIP5/PMIP3 were mis-scaled to give approximately twice the appropriate AOD forcing, although the spatial structure of forcing in the model is still coherent with G08. For this reason, we emphasize the CESM results in this study. However, we still choose to examine the results from the GISS-E2 model for two reasons. First, we view this error as an opportunity to explore the climate response to a wider range of hemispheric forcing gradients, even though it comes at the expense of not being able to relate the results to actual events during the LM. Secondly, the GISS-E2 LM runs were equipped with interactive water isotopes (Sect. 3.3). A self-consistent simulation of the isotope field in a GCM is important, since it removes a degree of uncertainty in the error-prone conversion of isotopic signals into more fundamental climate variables. To our knowledge, an explicit simulation of the isotopic distribution following asymmetries in volcanic forcing has not previously been reported.
In our analysis, we classify volcanic events as "symmetric" (SYMM), and "asymmetric" (ASYMM X ), where the subscript X refers to a preferred forcing in the NH or SH. Composites are formed from all events within each of the three classifications in order to isolate the volcanic sig- nal. All events must have a global aerosol loading > 8 Tg (1 Tg = 10 12 g) averaged over at least one 5-month period to qualify as an eruption and enter the composite. For comparison, the 1991 Mt. Pinatubo eruption remains elevated at ∼ 20-30 Tg sulfate aerosol in the G08 dataset for about a year and drops off to < 1 Tg after 4-5 years.
Events fall into the SYMM category if they have less than a 25 % difference in aerosol loading between hemispheres, while the ASYMM NH events have at least a 25 % higher loading in the NH relative to the SH. The opposite applies to events falling into the ASYMM SH category. The dates for which these thresholds are satisfied are taken from the original G08 dataset (Table 1), and thus the CESM and GISS-E2 composites are formed using the same events despite the GISS-E2 mis-scaling and other differences in model implementation.
Results are reported for the boreal warm season (averaged over the MJJAS months) and cold season (NDJFM), except for annual-mean results in Figs. 8-9 and for results showing the progression of signals at monthly resolution (Figs. S6, S9-S12 in the Supplement). For each eruption, we identify the post-volcanic response by averaging the number of consecutive seasons during which the above criteria are met, typically 1-3 years. All seasons for an eruption lasting longer than 1 year are first averaged together to avoid overweighting its influence in the composite. Anomalies are given with respect to the corresponding time of year during the 5 years prior to the eruption. For overlapping eruptions, the 5 years prior to the first eruption are used instead. This relatively short reference period allows creating composites that are unaffected by changes in the mean background state due to low-frequency climate change during the LM. Composites for the SYMM, ASYMM NH , and ASYMM SH cases are then obtained for each season and model by averaging over all anomaly fields within the appropriate classification, including all ensemble members. A two-sided Student's t test was applied to all composites in order to identify regions where the anomalous signal is significantly different (p < 0.05) from the mean background conditions.
In no case does the classification of a given eruption change over the duration of the event, with the exception of the largest eruption (Samalas, 1258 CE), which straddles the 25 % asymmetry criterion (SYMM and ASYMM NH ) throughout the years following the event. This eruption would project itself most strongly onto the symmetric composite but may reasonably be classified as ASYMM NH due to the greater absolute aerosol loadings in the NH. Due to this ambiguity, we omit the Samalas event from our main results. We note that there are far more asymmetric eruptions during the LM based on our criteria than SYMM cases, most of which easily meet the two thresholds outlined above. Because of this, the classification assigned to each event is quite robust to slightly different criteria in defining the ratio (or differences) in hemispheric aerosol loading. Since the asymmetric composites are formed from a relatively large number of events, our results are insensitive to the addition or removal of individual eruptions that may be more ambiguous in their degree of asymmetry. However, the SYMM composites are formed from only a few events and are therefore more sensitive to each of the individual eruptions that are included.

MJJAS
we are screening events by spatial structure and since different magnitude eruptions enter into the different composites, a quantitative comparison of the different event classifications (or the two models) is not our primary objective and would require a more controlled experiment. Instead, we are reporting on the different composite responses as they exist in current LM simulations and highlight the emergent structure that arises from different choices of how eruptions are sorted, much of which are shown to be scalable to different eruption sizes and robust to choices of model implementation. Figure 2 illustrates the composite temperature anomaly for each classification and season in the CESM model. In both the ASYMM NH and ASYMM SH cases, the hemisphere that is subjected to the strongest forcing is preferentially cooled. In the ASYMM NH results, the cooling peaks over the Eurasian and North American continents. As expected, there tends to be a much larger response over land, as well as evidence of NH winter warming in the mid-to-high latitudes, a phenomenon previously highlighted in the literature and often associated with increased (decreased) pole-toequator stratospheric (mid-tropospheric) temperature gradi-  (Robock andMao, 1992, 1995;Stenchikov et al., 2002;Shindell et al., 2004;Ortega et al., 2015). This effect is weak in the ASYMM NH composite, likely because the maximal radiative forcing is located in the NH, offsetting any dynamical response, but is present in the SYMM and ASYMM SH composites in both models (see Fig. S2 for the GISS-E2 composite).

Temperature, precipitation, and El Niño-Southern Oscillation (ENSO) response
In the SH, cooling is muted by larger heat capacity associated with a smaller land fraction, with weak responses over the Southern Ocean while still exhibiting statistically significant cooling in South America, South Africa, and Australia in all cases. In fact, the cooling in the ASYMM SH composites is largely confined to the tropics, in contrast to the polar amplified pattern that is common to most climate change experiments. The cooling in all categories is communicated vertically ( Fig. S1 in the Supplement) and across the free tropical troposphere, suggesting AET toward the forced hemisphere (Sect. 3.4) for asymmetric forcing.
The cooling in the GISS-E2 model (Fig. S2 in the Supplement) displays a very similar spatial structure to CESM in all categories but with much greater amplitude due to the larger forcing. We note that the composite-mean forcing is similar between the four asymmetric panels but larger in the symmetric cases. In Fig. 3, we show the hemispheric and global average temperature response for both models after normalizing each event by a common global aerosol mass excursion, thereby accounting for differences in the average forcing among the different eruptions. This is done to highlight spread associated with internal variability and model differ-ences, and it assumes that the response pattern scales linearly with global forcing, which is unlikely to be true across all events and for the two models. Nonetheless, the gross features of the hemispheric contrast and reduction in globalmean temperature are shared between both models.
The CESM precipitation response is shown in Fig. 4  (Fig. S3 in the Supplement for GISS-E2). For both the ASYMM NH and ASYMM SH cases, the ITCZ shows a robust displacement away from the forced hemisphere. The precipitation reduction in the SYMM composites is much less zonally coherent, instead featuring tropical-mean reductions in precipitation and a slight increase toward the subtropics (see also Iles et al., 2013;Iles and Hegerl, 2014). Despite global cooling and reduced global evaporation (not shown), the ITCZ shift in ASYMM NH and ASYMM SH tends to result in precipitation increases in the hemisphere that is least forced (Fig. 5), since the hemispheric-mean precipitation signal is largely influenced by the ITCZ migration itself.
The ensemble spread in precipitation for a selected eruption (1762 CE, NDJFM) is shown in Fig. S4 in the Supplement, corresponding to the Icelandic Laki aerosol loading (a large ASYMM NH event). We note that the Laki eruption in Iceland actually occurred in 1783 CE but is earlier in our composite due to an alignment error in the first version of the G08 dataset. Results are shown for the 1763 CE boreal winter only (the full composite also includes 1762, see Table 1; Fig. S4 also reports the winter 1763 Niño 3.4 anomaly in surface temperature for each ensemble member, and therefore we restrict the anomalous precipitation field to the same season). The ITCZ shift away from the NH is fairly robust across the ensemble members, particularly in the Atlantic basin, although internal variability still leads to large differences in the spatial pattern of precipitation, notably in the central and eastern Pacific.
The monthly time evolution of the composite temperature and precipitation responses for the ASYMM NH and ASYMM SH cases can be viewed in an animation (Figs. S9-S12 in the Supplement). The global and hemispheric difference in aerosol loadings is also shown for each timestep (at monthly resolution) in the animations. When averaged over the individual eruptions within each classification, the global aerosol mass loading remains elevated above 8 Tg for nearly 2 years, coincident with the peak temperature and precipitation response that begins to dampen out gradually and relaxes back to pre-eruption noise levels after ∼ 4-5 years. The seasonal migration of anomalous precipitation in the ITCZ domain occurs in nearly the same way as the meridional movement of climatological rainfall, highlighting important connections between the timing of the eruption relative to the seasonal cycle of rainfall at a given location.
In both CESM and GISS-E2, the ITCZ shift is approximately scalable to eruption size. For both models, we define a precipitation asymmetry index, PAi , in each season as the area-weighted NH tropical precipitation minus SH tropical precipitation (extending to 20 latitude) normalized by the model tropical-mean precipitation, i.e., Supplement Fig. S5 illustrates the relationship between PAi and the AOD gradient between hemispheres (AOD is inferred for the CESM model by dividing the aerosol loading by 75 Tg in each hemisphere, an approximate conversion factor to compare the results with GISS-E2). The misscaling in GISS-E2 results in a wider range of AOD gradients than occurs in CESM. Both models feature more tropical precipitation in the NH (SH) during boreal summer (winter) in their climatology, with more asymmetry in CESM during boreal summer. Interestingly, the most asymmetric events in GISS-E2 (those that result in equatorward precipitation movements) can be sufficient to produce more precipitation in the tropical winter hemisphere, thus competing with the seasonal insolation cycle in determining the seasonal precipitation distribution. The meridional ITCZ shift leads to a number of important tropical climate responses. For example, an intriguing feature of the temperature pattern in Fig. 2 is the El Niño response that is unique to the ASYMM NH composites. This is unlikely to be a residual feature of unforced variability, since there are 288 events in the ASYMM NH composites (16 eruptions in Table 1, multiplied by 18 ensemble members), significantly more than in the other categories. The GISS-E2 temperature composite (Fig. S2 in the Supplement) also features a relatively weak cooling for ASYMM NH , despite the very large radiative forcing. The relationship between ENSO and volcanic eruptions has, historically, been quite complicated due to the problem of separating natural variability from the forced response and due to a limited sample of historical eruptions where ENSO events were already underway prior to the eruption. Older studies have suggested that El Niño events may be more likely 1 to 2 years following a large eruption (e.g., Adams et al., 2003;Mann et al., 2005;Emile-Geay et al., 2008). Our findings are also consistent with recent results (Pausata et al., 2015) that found an El Niño tendency to arise from a Laki-like forcing (in that study, a sequence of aerosol pulses in the high latitudes that was confined to the NH extratropics), and the El Niño response in CESM to different expressions of volcanic forcing was recently explored in Stevensonn et al. (2016). Pausata et al. (2015) attributed the El Niño development directly to a southward ITCZ displacement. Since low-level converging winds are weak in the vicinity of the ITCZ, a southward ITCZ displacement leads to weaker easterly winds (a westerly anomaly) across the central equatorial Pacific. This was shown for a different model (NorESM1-M) and experimental setup but also emerges in the ASYMM NH composite results for CESM. Indeed, a composite anomaly of ∼ 0.5 • C emerges over the Niño 3.4 domain, lasting up to 2 years (Fig. S6 in the Supplement) with peak anomalies in the first two boreal winters after an eruption. Consistent with the sea surface temperature (SST) anomalies, a relaxation of the zonal winds and redistribution of water mass across the Pacific Ocean can be observed in the ASYMM NH composite response (Fig. S7 in the Supplement).
Since the ITCZ shift is a consequence of differential aerosol loading, we argue that the El Niño tendency in CESM is a forced response in ASYMM NH but otherwise depends on the state of internal variability concurrent with a given eruption, as no such ENSO response is associated with the composite SYMM or ASYMM SH , although we note that El Niño does tend to develop in response to the Samalas eruption that was removed from our composite and would strongly influence the interpretation of the SYMM results due to the few events sampled (not shown, though see Stevenson et al., 2016). However, we also caution that this version of CESM exhibits ENSO amplitudes much larger than observations and also features strong El Niño events with amplitudes that are ∼ 2 times larger than strong La Niña events even in noneruption years. Therefore, we choose not to further explore the dependence of our results on ENSO phasing.
Because the ITCZ responds differently to the three eruption classifications, there are implications for best practices in assessing the skill of climate model output against proxy evidence. For example, Anchukaitis et al. (2010) noted discrepancies between well-validated tree-ring proxy reconstructions of eruption-induced drought in the Asian monsoon sector and the precipitation response following volcanic eruptions derived from the NCAR Climate System Model (CSM) 1.4 millennial simulation. However, we note that monsoonal rainfall responds differently to ASYMM NH , ASYMM SH , or SYMM events in both GISS-E2 and CESM. Figure S8 in the Supplement shows a histogram of boreal summer (MJJAS) Asian Pacific rainfall anomalies for all events in both models. ASYMM NH and SYMM eruptions generally lead to reductions in rainfall over the broad re-gion averaged from 65 to 150 • E and 10-40 • N (see also the spatial patterns in Fig. 4 for CESM and Fig. S3 in the Supplement for GISS E2-R). Because of the southward ITCZ shift in ASYMM NH , the most pronounced precipitation reductions occur for events within this category. In contrast, for ASYMM SH events, the northward ITCZ shift and associated monsoon developments are such that precipitation changes are relatively muted, and often the anomalies are positive.

River outflow
An ITCZ shift away from the forced hemisphere will manifest itself in several other components of the tropical hydroclimate system that are important to consider from the standpoint of both impacts as well as the development of testable predictions. One such important component of the hydrologic cycle is global streamflow, a variable that is related to excessive or deficient precipitation over a catchment. Rivers are important for ecosystem integrity, agriculture, industry, power generation, and human consumption. Streamflow anomalies associated with volcanic forcing in observations and models have previously been documented for the historical period (Trenberth and Dai, 2007;Iles and Hegerl, 2015). Here, we discuss this variable in the context of our symmetric and asymmetric composites.
The hydrology module of the land component of CESM simulates surface and subsurface fluxes of water, which serve as input into the CESM River Transport Model (RTM). The RTM was developed to route river runoff downstream to the ocean or marginal seas and enable closure of the hydrologic cycle (Oleson et al., 2010). The RTM is run on a finer grid (0.5 • × 0.5 • ) than the atmospheric component of CESM. Figure 6 shows the river discharge anomalies in our different forcing categories. The southward ITCZ shift in ASYMM NH results in enhanced discharge in central and southern South America, especially in the southern Amazon and Paraná river networks. These territories of South America, along with southern Africa and Australia, are the primary regions where land precipitation increases in the tropics for ASYMM NH , and the river flow in these areas tends to increase. Our results are also consistent with Oman et al. (2006), who argue for a reduced Nile River level (northeastern Africa) following several large high northern latitude eruptions, including Laki and the Katmai (1912 CE) eruption. Their results were viewed through the lens of weakened African and Indian monsoons associated with reduced landocean temperature differences; our composite results suggest that regional precipitation reductions may also be part of a zonally coherent precipitation shift.
In ASYMM SH , the ITCZ moves northward, resulting in reduced river flux in the Amazon sector and increases (reduction) in the Niger of central and western Africa during boreal summer (boreal winter). Interestingly, the Nile flow is also reduced in this case, although to a lesser extent, despite very modest precipitation increases during MJJAS for a Southern Hemisphere biased aerosol forcing. There are also modest discharge increases in southern Asia. However, there is simply very little land in regions where northward ITCZ shifts result in enhanced precipitation, suggesting less opportunity for increases in discharge to a SH biased eruption. For the SYMM eruptions, river discharge is reduced nearly everywhere in the tropics, consistent with the precipitation reductions that occur (Fig. 3). The response is weaker or even reversed in the subtropics, such as in southern South America, where precipitation tends to increase (Iles and Hegerl, 2015).

Water isotopic variability
Another important variable that integrates several aspects of the tropical climate system is the isotopic composition of precipitation. Here, we focus on the relative abundance of 1 H 18 2 O versus the more abundant 1 H 16 2 O, commonly expressed as δ 18 O, such that where O 18 mp and O 16 mp are the moles of oxygen isotope in a sample, in our case precipitation (denoted by the subscript mp). Delta values are given with respect to the isotopic ratio in a standard sample, the Vienna Standard Mean Ocean Water (VSMOW = 2.005 × 10 −3 ). δ 18 O p is a variable that is directly obtained from many paleoclimate proxy records. Therefore, rather than relying on a conversion of the local isotope signal to some climate variable, the explicit simulation of isotopic variability is preferred for generating potentially falsifiable predictions concerning the imprint associated with asymmetric volcanic eruptions. Indeed, δ 18 O p variability is 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.
Water isotope tracers have been incorporated into the GISS-E2 model's atmosphere, land surface, sea ice, and ocean and are advected and tracked through every stage of the hydrologic cycle. A fractionation factor is applied at each phase change and all freshwater fluxes are tagged isotopically. Stable isotope results from the lineage of GISS-E2 models have a long history of being tested against observations and proxy records (e.g., Vuille et al., 2003;Schmidt et al., 2007;Schmidt, 2008, 2009). Figure 7 shows the δ 18 O p response in the GISS-E2 model. Seasonal calculations are weighted by the precipita- , except for ASYMM SH eruptions. Grey envelope corresponds to the total AET anomaly vs. latitude in a control simulation using 50 realizations of a 17-event composite (17 "events" with no external forcing, corresponding to the size of the ensemble). Vertical bars correspond to the range of (aqua) latent and (orange) dry components of cross-equatorial energy transport (AET eq ) in the control composite.
tion amount for each month, although changes in the seasonality of precipitation are not important in driving our results (not shown). The literature on mechanistic explanations for isotope variability has a rich history of being described by several "effects" such as a precipitation amount effect in deep convective regions or a temperature effect at high latitudes (Dansgaard, 1964;Araguás-Araguás et al., 2000), so named as to reflect the most important climatic driver of isotopic variability at a site or climate regime. Notably, δ 18 O p tends to be negatively correlated with precipitation amount in the deep tropics and positively correlated with temperature at high latitudes (see, e.g., Hoffman and Heimann (1997) for a review of mechanisms). However, isotope-climate relations are generally complex. In our experiments, the δ 18 O p spatial pattern in the tropics (Fig. 7) exhibits a similar pattern to precipitation changes induced by the ITCZ shift (Fig. S5 in the Supplement for GISS-E2), particularly over the ocean. The meridional movement of the ITCZ leads to an isotopic signal that is more positive (enriched in heavy isotopes) in the preferentially forced hemisphere. The hemisphere toward which the ITCZ is displaced on the other hand experiences increased tropical rainfall and a relative depletion of the heavy isotope (more negative δ 18 O p ). Thus, the paleoclimatic fin-gerprint of asymmetric volcanic eruptions is characterized by a tropical dipole pattern, with more positive (negative) δ 18 O p associated with reduced (increased) rainfall. Over land, South America stands out as exhibiting a palette of isotopic patterns depending on forcing category and season. The South American monsoon system peaks in austral summer, and the largest precipitation reductions occur in ASYMM SH when the ITCZ moves northward. There is a dipole pattern, characterized by isotopic enrichment (depletion) in 18 O in the northern (southern) tropics of South America in ASYMM NH during NDJFM, while the opposite pattern emerges in ASYMM SH , both associated with Atlantic and east Pacific ITCZ displacements. During the austral winter, climatological South American precipitation peaks in the northern part of the continent, and precipitation in this region is reduced in both the SYMM and ASYMM SH composites, leading to a large increase in δ 18 O p . This is consistent with recent results in Colose et al. (2016), who used the isotope-enabled GISS-E2 model to form a composite of all large (AOD > 0.1) LM tropical volcanic events based on the Crowley and Unterman (2013) dataset. The eruptions analyzed in that study were smaller in amplitude due to differences in the scaling during implementation, as well as the fact that G08 tends to have larger volcanic events in the original dataset to begin with. In regions where tropical South American precipitation does not exhibit very large changes, such as in the NDJFM SYMM composites, temperature may explain much of the isotopic response, again consistent with findings in Colose et al. (2016).

Atmospheric energetics
The overarching purpose of this work was to consider the influence of asymmetric volcanic forcing on the energetic paradigm outlined in Sect. 1. This framework of analyzing ITCZ shifts in the context of asymmetric forcing predicts a net AET anomaly toward the hemisphere that is preferentially forced by explosive volcanism, with anticorrelated dry and latent energy fluxes both contributing to drive the ITCZ away from the forced hemisphere. To examine this relationship in CESM, we first write a zonal-mean energy budget for the atmosphere (Trenberth, 1997;Donohoe and Battisti, 2013): where ASR TOA is the absorbed solar radiation, OLR TOA is outgoing longwave radiation at the top of the atmosphere (TOA), SW ↑ sfc is reflected surface shortwave radiation, SW ↓ sfc is shortwave received by the surface (sfc), LW ↑ sfc is longwave radiation emitted (or reflected) by the surface, LW ↓ sfc is longwave radiation received by the surface, LH is the latent heat flux, SH is the sensible heat flux, Sn is snowfall rate, q is specific humidity, k is kinetic energy, φ is latitude, a is the radius of the Earth, T is temperature, c p is specific heat capacity, L v and L f are the latent heats of vaporization and fusion, p is pressure (p = p s at the surface), and g is the acceleration due to gravity. All terms are defined positive in the atmosphere, and the subscripts denote TOA or surface flux (sfc) diagnostics. Equation (3) effectively calculates MSE transport (Sect. 1) as a residual of energy fluxes in the model.
The last term (∂ ∂t ) on the right side of Eq. (3) is the time-tendency term, representing storage of energy in the atmosphere. (Hereafter, STOR L and STOR D for latent and dry energy, respectively. The time derivative is calculated using finite differencing of the monthly-mean fields. The term in the parentheses is the moist enthalpy, or MSE minus geopotential energy. The kinetic energy is calculated in this study but is several orders of magnitude smaller than other terms, and hereafter is folded into the definition of STOR D .) The tendency term must vanish on timescales of several years or longer but is important in our context. We explicitly write out the snowfall term since CESM (and any CMIP5 model) does not include surface energy changes associated with snowmelt over the ice-free ocean as part of the latent heat diagnostic and must be calculated to close the model energy budget.
Integrating yields an expression for the atmospheric heat transport across a latitude circle: where we have combined the TOA terms into R TOA and the snowfall and surface diagnostics have collapsed into a single variable F sfc . Similarly, the latent heat flux H L across a latitude circle is where P is precipitation in kg m −2 s −1 . We note that transport calculations are presented for CESM and were done for only 17 ensemble members, since there are missing output files for the requisite diagnostics in one run. Figure 8a shows the annual-mean climatological northward heat transport in CESM, as performed by the atmosphere, in addition to the dry and moisture-related components of AET. The total CESM climatological poleward transport is in good agreement with observational estimates (e.g., Trenberth and Caron, 2001;Wunsch, 2005;Fasullo and Trenberth, 2008), peaking at ∼ 5.0 PW and ∼ 5.2 PW in the SH and NH subtropics, respectively (1 PW = 10 15 W). In CESM, the SH receives slightly more net TOA solar radiation than the NH (by ∼ 1.3 W m −2 in the annual-mean), and the NH loses slightly more net TOA longwave radiation to space (by ∼ 0.89 W m −2 ). However, the CESM annual ocean heat transport is northward across the equator (not shown), keeping the NH warmer than the SH by ∼ 0.98 • C. As a consequence, AET is directed southward across the equator (red line). Moisture makes it more difficult for the tropical circulation to transport energy poleward, and the transport of moisture in the low-level equatorward flow is directed northward across the equator and associated with an annual-mean ITCZ approximately co-located with the atmospheric energy flux equator (EFE), the latitude where AET vanishes. This arrangement of the tropical climate is consistent with satellite and reanalysis results for the present climate Kang et al., 2014).
In response to asymmetric volcanic forcing, anomalous AET is directed toward the preferentially forced hemisphere (Fig. 8b, c), along the imposed temperature gradient. Results are shown for the annual-mean AET anomaly in ASYMM NH and ASYMM SH for 1 year beginning with the January after each eruption, although averaging the first 2-3 years yields similar results with slightly smaller amplitudes. The equatorial AET (AET eq ) anomaly averaged over all events and ensemble members for ASYMM NH (ASYMM SH ) is approximately 0.08 (−0.06) PW, defined positive northward, with much larger near-compensating dry and latent components. The anomalous moisture convergence drives the ITCZ shift away from the forced hemisphere. Anomalies in AET eq when considering each unique volcanic event (after averaging over the 17 ensemble members) are strongly anticorrelated with changes in the energy flux equator (r = −0.97, not shown), the latitude where AET vanishes.
The change in cross-equatorial energy transport for the SYMM ensemble/eruption mean (not shown) does not exhibit the coherence of the asymmetric cases for either AET or the individual dry and moist components and in all cases does not emerge from background internal variability.
Quantifying the ITCZ shift is nontrivial, since the precipitation field is less sharply defined than the EFE, and climate models (including the two discussed here) exhibit a bimodal tropical precipitation distribution (often called a "double ITCZ"), often with one mode of higher amplitude in the NH (centered at 8-9 • N in CESM). However, despite pervasive biases that still exist in the climatology of tropical precipitation in CMIP5 (e.g., Oueslati and Bellon, 2015), the anomalous precipitation response is still characterized by a well-defined ITCZ shift (or a shift in the bimodal precipitation distribution; e.g., Fig. 9 in Stevenson et al., 2016), and the gross features presented here are in agreement with theoretical considerations. In our analysis, a movement in the latitude of maximum precipitation is not found to be a persuasive indicator of our ITCZ shift. In fact, the meridional shift is better described as a movement in the center of mass of the precipitation distribution, including changes in the relative amplitude of the two modes (e.g., a heightening of the SH mode for a southward ITCZ shift). Different metrics to describe the shift in the center of mass have been presented in the literature (e.g., Frierson and Hwang, 2012;Adam et al., 2016).
Here, we first adopt the precipitation median φ med (e.g., Frierson and Hwang, 2012) defined as the latitude where area-weighted precipitation from 20 • S to φ med equals the precipitation amount from φ med to 20 • N, i.e., where the following is satisfied: When considering the spread across eruption size (regressing the different events in all three categories together after averaging over ensemble members), we find a ∼ −8.9 • shift in ITCZ latitude per 1 PW of anomalous AET eq (Fig. 9). The sign of this relationship is a robust property of the present climate system, although it is higher than other estimates  that analyzed the ITCZ scaling with AET eq to a number of other time periods and forcing mechanisms (not volcanic), including the seasonal cycle, CO 2 doubling, Last Glacial Maximum, and mid-Holocene. It was argued in  that the ITCZ is "stiff" in the sense that a large AET eq is required to move the ITCZ. However, the sensitivity of this relationship may vary considerably depending on the ITCZ metric considered ( Fig. 9 presents a scaling with different indices), based on the following equation (Adam et al., 2016): Here, N controls the weighting given to the modes in the precipitation distribution. Typically φ ITCZ moves toward the precipitation maximum as N increases, but importantly, the sensitivity of a φ ITCZ migration to a given anomaly in AET eq also changes. Figure 9 shows the regression of anomalous φ med and φ ITCZ (N = 5) against anomalous AET eq (r = −0.94). φ ITCZ (N = 3) yields a high correlation (r = −0.95) and best follows a 1 : 1 line with the EFE (Fig. 9, bottom left). The slope of the relationship between ITCZ location and AET eq may vary by a factor of 4-5 depending on the relationship used. For example, there is approximately a −11.7 • shift in ITCZ latitude per 1 PW of anomalous AET eq using φ ITCZ (N = 3). Thus, we interpret our results as suggesting that, energetically, it is not necessarily difficult to move the ITCZ and urge caution in characterizing past ITCZ shifts as being difficult to reconcile with paleo-forcing estimates . Indeed, as many studies have used a "precipitation centroid" or a similar variant to quantify tropical precipitation migrations, we recommend exploring the sensitivity of ITCZ shifts to different ways of characterizing the movement in precipitation mass unless the community can agree upon a well-defined N that suitably characterizes the precipitation distribution in both climate models and observations.

Conclusions
In this work, we have examined two models -NASA GISS ModelE2-R and the recently completed CESM Last Millennium Ensemble -and stratified volcanic events by their degree of asymmetry between hemispheres. We find a robust ITCZ shift away from the preferentially forced hemisphere as a consequence of adjustments in the Hadley circulation that transports anomalous energy into the cooled hemisphere. An important component of our work was using the GISS-E2 model to explicitly simulate the oxygen isotopic imprint following major volcanic eruptions with asymmetric aerosol forcing. The ITCZ shift following asymmetric forcing leads to a more positive isotopic signal in the tropical regions the ITCZ migrates away from and a relative depletion in heavy isotopes in regions the ITCZ migrates to. These results provide a framework for the search of asymmetric volcanic signals in high-resolution isotopic or other temperature-and precipitation-sensitive proxy data from the tropics.
There is still considerable uncertainty in the timing and magnitude of LM eruptions. Improvements in particle size representation have been identified as a critical target for improved modeling and comparisons to proxy data (e.g., Mann et al., 2015). Here, we argue that the interhemispheric asymmetry of the aerosol forcing also emerges as being of firstorder importance for the expected volcanic response. Future developments in model-proxy comparisons should probe the uncertainty space not just in the global-mean radiative forcing and coincident internal variability at the time of the eruption but also in the spatial structure of the aerosol cloud. For example, simulations that represent volcanic forcing simply as an equivalent reduction in total solar irradiance at the TOA are unrealistic and cannot be expected to be faithful to tropical climate proxy records.
We hope this contribution will help motivate the connection between the spatial structure of volcanic episodes and the expression in tropical hydroclimate as an urgent paleoclimate target in future studies and model intercomparisons. Such investigation also calls for high-resolution and accurately dated tropical proxy networks that reach across hemispheres. Developments in seasonally and annually resolved volcanic reconstructions from both hemispheres (Sigl et al., 2015) are of considerable importance in such assessments. Future modeling efforts that are forced with the explicit injection of volcanic species, while also probing multiple realizations of internal variability that will dictate the spatiotemporal evolution of the volcanic aerosol, are also urgently required as a tool for understanding both past and future volcanic impacts.

Data availability
Volcanic input for 850-1850 CE can be obtained from http: //climate.envsci.rutgers.edu/IVI2/ (original dataset no. 2). ModelE2-R climate model data included in this study can be downloaded from the Earth System Grid Federation (account required); we used three of the past1000 runs described in http://data.giss.nasa.gov/modelE/ar5/ ("GRA" for volcanic forcing). The water isotope variable for the model runs is not available publicly. Written correspondence can be sent to Allegra N. LeGrande (allegra.n.legrande@nasa.gov) to obtain this diagnostic. Diagnostics for the NCAR last millennium simulations are available at https://www.earthsystemgrid. org/dataset/ucar.cgd.ccsm4.CESM_CAM5_LME.html.
The Supplement related to this article is available online at doi:10.5194/esd-7-681-2016-supplement.